# Chunks of weak compositions

Here is slightly modified version of algorithm used in ``Combinatorica`NextComposition`` converted to a `LibraryFunction`.

 Needs["CCompilerDriver`"]
 "
 #include \"WolframLibrary.h\"
 
 DLLEXPORT mint WolframLibrary_getVersion() {
 return WolframLibraryVersion;
 }
 DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) {
 return LIBRARY_NO_ERROR;
 }
 DLLEXPORT void WolframLibrary_uninitialize(WolframLibraryData libData) {}
 
 DLLEXPORT int nextWeakCompositionsChunk(
 WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res
 ) {
 // number which compositions are calculated
 mint n = MArgument_getInteger(Args[0]);
 
 // length of composition
 mint k = MArgument_getInteger(Args[1]);
 
 // \" current \" composition tensor
 MTensor compositionT = MArgument_getMTensor(Args[2]);
 // actual data of \" current \" composition
 mint* composition = libData->MTensor_getIntegerData(compositionT);
 
 // one-dimentional, one-element tensor holding index of last non-zero 
 // element in \"current\" composition
 MTensor lastNonZeroIndexT = MArgument_getMTensor(Args[3]);
 // pointer to index of last non-zero element in \"current\" composition
 mint* lastNonZeroIndex = libData->MTensor_getIntegerData(lastNonZeroIndexT);
 
 // number of compositions in returned chunk
 mint chunkSize = MArgument_getInteger(Args[4]);
 
 // dimensions of returned chunk tensor
 mint chunkDims[2] = {chunkSize, k};
 // 2 dimentional tensor with chunk of compositions to be returned
 MTensor chunkT;
 libData->MTensor_new(MType_Integer, 2, chunkDims, &chunkT);
 // actual data of the chunk tensor
 mint* chunk = libData->MTensor_getIntegerData(chunkT);
 
 // index enumerating compositions in chunk
 mint i;
 // index enumerating elements of composition
 mint j;
 
 for (i = 0; i < chunkSize; i++) {
 // Moving to next composition does basically this:
 // {..., x, y} -> {..., x+1, y-1} if last element is non-zero,
 // {..., x, y, 0, ..., 0} -> {..., x+1, 0, 0, ..., y-1}
 // if last element is zero and y is last non-zero element.
 
 ++(composition[*lastNonZeroIndex - 1]);
 
 if (*lastNonZeroIndex == k -1) {
 --(composition[*lastNonZeroIndex]);
 } else {
 composition[k - 1] = composition[*lastNonZeroIndex] - 1;
 composition[*lastNonZeroIndex] = 0;
 }
 
 // If last element is zero it means we moved last non-zero element to
 // the left. If last element is non zero, then, obviously, it is the
 // last non-zero element.
 if (composition[k - 1] == 0) {
 --(*lastNonZeroIndex);
 } else {
 *lastNonZeroIndex = k - 1;
 }
 
 // Copy \"current\" composition to i-th row of chunk.
 for (j = 0; j < k; j++) {
 chunk[i*k + j] = composition[j];
 }
 }
 
 // Return control over composition and lastNonZeroIndex tensors back to
 // Wolfram Language.
 libData->MTensor_disown(lastNonZeroIndexT);
 libData->MTensor_disown(compositionT);
 
 // Set chunk tensor as returned value of LibraryFunction.
 MArgument_setMTensor(Res, chunkT);
 
 return LIBRARY_NO_ERROR;
 }
 ";
 CreateLibrary[%, "weakCompositions"]
 
 ClearAll[nextWeakCompositionsChunk]
 nextWeakCompositionsChunk = 
 LibraryFunctionLoad[
 %%, 
 "nextWeakCompositionsChunk",
 {Integer, Integer, {Integer, 1, "Shared"}, {Integer, 1, "Shared"}, Integer},
 {Integer, 2}
 ]

It returns chunk of weak `k`-compositions of `n`, in lexicographic order starting from composition after given one.
It accepts 5 arguments: numbers `n` and `k`, list of integers representing "previous" composition, list containing one integer being index (in C sense i.e. starting from 0) of last non-zero element in "previous" composition, last argument is number of compositions in returned chunk.

Besides returning chunk of compositions `nextWeakCompositionsChunk` has a side effect.
"Previous" composition and list containing index of last non-zero element are passed as "Shared" tensors.
They are modified inside `nextWeakCompositionsChunk` function and they stay modified.
For example get 3 weak 4-compositions of 5 starting from one after `{1, 3, 0, 1}`:

 ClearAll[currentComposition, lastNonZeroIndex]
 currentComposition = {1, 3, 0, 1} // Developer`ToPackedArray;
 lastNonZeroIndex = {3} // Developer`ToPackedArray;
 nextWeakCompositionsChunk[5, 4, currentComposition, lastNonZeroIndex, 3]
 (* {{1, 3, 1, 0}, {1, 4, 0, 0}, {2, 0, 0, 3}} *)
 currentComposition
 (* {2, 0, 0, 3} *)
 lastNonZeroIndex
 (* {3} *)

To get first composition we need to pass a fake composition that implemented algorithm will change to the first one.
This "fake one before first composition" is always of form `{0, ..., 0, -1, n+1}`.
For example to get all 10 weak 3-compositions of 3 we use: 

 ClearAll[currentComposition, lastNonZeroIndex]
 currentComposition = {0, -1, 4} // Developer`ToPackedArray;
 lastNonZeroIndex = {2} // Developer`ToPackedArray;
 nextWeakCompositionsChunk[3, 3, currentComposition, lastNonZeroIndex, 10]
 (* {{0, 0, 3}, {0, 1, 2}, {0, 2, 1}, {0, 3, 0}, {1, 0, 2}, {1, 1, 1}, {1, 2, 0}, {2, 0, 1}, {2, 1, 0}, {3, 0, 0}} *)
 currentComposition
 (* {3, 0, 0} *)
 lastNonZeroIndex
 (* {0} *)

This is a "low level" function (in package it would be private) that does no tests of consistency of given arguments.
It is intended to be used only through functions defined in next sections.
Passing wrong arguments can result in kernel crash.


# All weak compositions

Using our low level function we can define a function that returns all weak compositions.

 ClearAll[weakCompositions]
 weakCompositions[n_Integer?NonNegative, 0] := {}
 weakCompositions[n_Integer?NonNegative, 1] := {{n}}
 weakCompositions[n_Integer?NonNegative, k_Integer /; k >= 2] :=
 nextWeakCompositionsChunk[
 n,
 k,
 Join[ConstantArray[0, k - 2], {-1, n + 1}] // Developer`ToPackedArray,
 {k - 1} // Developer`ToPackedArray,
 Binomial[n + k - 1, n]
 ]

For example get all weak 3-compositions of 2:

 weakCompositions[2, 3]
 (* {{0, 0, 2}, {0, 1, 1}, {0, 2, 0}, {1, 0, 1}, {1, 1, 0}, {2, 0, 0}} *)

Test that we get same output as sorted `Permutations` of `IntegerPartitions`:

 And @@ Flatten @ Table[
 SameQ[
 weakCompositions[n, k],
 weakCompositionsPermPart[n, k] // Sort
 ],
 {n, 12}, {k, 12}
 ]
 (* True *)

`weakCompositions` is faster and more memory efficient than combination of built-ins.

Below are some benchmarks to prove it.

Since `Join`ing of permutations of partitions can double memory usage, and in some cases is not needed, here we use function that doesn't join them.

 ClearAll[weakCompositionsPermPartNonJoined]
 weakCompositionsPermPartNonJoined[n_, k_] :=
 Permutations /@ IntegerPartitions[n, {k}, Range[0, n]]

Some helper functions:

 ClearAll[$functionsList, fixedNFunctionsList, fixedKFunctionsList]
$functionsList = {weakCompositionsPermPartNonJoined, weakCompositions};
 fixedNFunctionsList[n_] := Function[func, func[n, #] &] /@ $functionsList
fixedKFunctionsList[k_] := Function[func, func[#, k] &] /@ $functionsList

We're benchmarking both time and memory usage using my [`timeAndMemoryBenchmarkPlots` function](http://mathematica.stackexchange.com/a/90034/14303).

Fixed `n`, variable `k`:

 timeAndMemoryBenchmarkPlots[fixedNFunctionsList[2], Identity, {1, 2, 3, 4, 5, 8, 11, 15, 21, 29, 41, 58, 81, 114, 160, 225, 315, 442, 620, 870}, TimeConstraint -> Infinity]
 timeAndMemoryBenchmarkPlots[fixedNFunctionsList[8], Identity, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 17, 19, 21, 23, 25}, TimeConstraint -> Infinity]
 timeAndMemoryBenchmarkPlots[fixedNFunctionsList[32], Identity, Range[8], TimeConstraint -> Infinity]
 
[![benchmarks weakCompositions fixed n](https://i.sstatic.net/S0qQX.png)](https://i.sstatic.net/S0qQX.png)


Fixed `k`, variable `n`:

 timeAndMemoryBenchmarkPlots[fixedKFunctionsList[2], Identity, {1, 2, 4, 7, 13, 25, 47, 90, 171, 324, 617, 1172, 2229, 4237, 8054, 15312, 29109, 55338, 105203, 200000}, TimeConstraint -> Infinity]
 timeAndMemoryBenchmarkPlots[fixedKFunctionsList[8], Identity, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 17, 19, 22, 25, 28, 32, 36}, TimeConstraint -> Infinity]
 timeAndMemoryBenchmarkPlots[fixedKFunctionsList[32], Identity, Range[6], TimeConstraint -> Infinity]
 
[![benchmarks weakCompositions fixed k](https://i.sstatic.net/9ludl.png)](https://i.sstatic.net/9ludl.png)
 

# Lazy weak compositions

And now the lazy approach:

 Needs["Streaming`"]

 ClearAll[lazyWeakCompositions]

 lazyWeakCompositions[
 n_Integer?NonNegative, 0, chunkSize : _Integer?Positive : 10^5,
 opts : OptionsPattern[]
 ] :=
 LazyListCreate[{}, chunkSize, opts]

 lazyWeakCompositions[
 n_Integer?NonNegative, 1, chunkSize : _Integer?Positive : 10^5,
 opts : OptionsPattern[]
 ] :=
 LazyListCreate[{{n}}, chunkSize, opts]

 lazyWeakCompositions[
 n_Integer?NonNegative, k_Integer /; k >= 2,
 chunkSize : _Integer?Positive : 10^5, opts : OptionsPattern[]
 ] :=
 With[{length = Binomial[n + k - 1, n]},
 Module[
 {
 active = False,
 compositionsLeft = length,
 prevComposition = 
 Join[ConstantArray[0, k - 2], {-1, n + 1}] // 
 Developer`ToPackedArray,
 lastNonZeroIndex = {k - 1} // Developer`ToPackedArray
 }
 ,
 LazyListCreate[
 IteratorCreate[
 ListIterator,
 (active = True) &,
 Module[{realChunkSize = Min[chunkSize, compositionsLeft]},
 If[realChunkSize === 0,
 {}
 ,(* else *)
 compositionsLeft -= realChunkSize;
 nextWeakCompositionsChunk[
 n, k, prevComposition, lastNonZeroIndex,
 realChunkSize
 ]
 ]
 ] &,
 TrueQ[active] &,
 Remove[
 active, compositionsLeft,
 prevComposition, lastNonZeroIndex
 ] &
 ],
 chunkSize,
 opts,
 "Length" -> length,
 "FiniteList" -> True
 ]
 ]
 ]

Again all weak 3-compositions of 2, but this time taken in chunks of 2:

 lazyWeakCompositions[2, 3, 2]
 (* « LazyList[{0,0,2},{0,1,1},...] » *)
 Scan[Print, %]
 (* {0,0,2}
 {0,1,1}
 {0,2,0}
 {1,0,1}
 {1,1,0}
 {2,0,0} *)

Test that `Normal`ized lazy compositions are the same as sorted `Permutations` of `IntegerPartitions`:

 And @@ Flatten @ Table[
 SameQ[
 lazyWeakCompositions[n, k] // Normal,
 weakCompositionsPermPart[n, k] // Sort
 ],
 {n, 12}, {k, 12}
 ]
 (* True *)

We will benchmark scanning of all weak compositions in three ways:

1. All compositions generated at once:

 ClearAll[scanWeakCompositionsPermPartNonJoined]
 scanWeakCompositionsPermPartNonJoined[n_, k_] := (
 Scan[Identity, weakCompositionsPermPartNonJoined[n, k], {2}];
 DeleteFile@FileNames["*.mx", $StreamingCacheBase]
 )

2. Compositions generated by permuting lazy partitions.
 Partitions are taken one by one, so in single chunk we have all permutations of given partition.

 ClearAll[lazyWeakCompositionsPermPart, scanLazyWeakCompositionsPermPart]
 lazyWeakCompositionsPermPart[n_, k_] := 
 Permutations /@ lazyIntegerPartitions[n, {k}, Range[0, n], 1]

 scanLazyWeakCompositionsPermPart[n_, k_] := (
 Scan[Scan[Identity], lazyWeakCompositionsPermPart[n, k]];
 DeleteFile@FileNames["*.mx", $StreamingCacheBase]
 )

3. Compositions generated using our `lazyWeakCompositions` function in chunks of default length: 10^5.

 ClearAll[scanLazyWeakCompositions]
 scanLazyWeakCompositions[n_, k_] := (
 Scan[Identity, lazyWeakCompositions[n, k]];
 DeleteFile@FileNames["*.mx", $StreamingCacheBase]
 )

Since lazy lists cache their values and benchmark measures running time multiple times, we clear cache after each usage.

 timeAndMemoryBenchmarkPlots[fixedKFunctionsList[16], Identity, Range[16], TimeConstraint -> Infinity, MemoryConstraint -> 2^31]
 timeAndMemoryBenchmarkPlots[fixedKFunctionsList[32], Identity, Range[8], TimeConstraint -> Infinity, MemoryConstraint -> 2^31]
 
[![enter image description here](https://i.sstatic.net/SBVYU.png)](https://i.sstatic.net/SBVYU.png)

We see that speed for large values is comparable.
Memory usage when generating all compositions at once grows fast.
When using lazy integer partitions memory usage grows slower, but still grows with number of permutations.
When using `lazyWeakCompositions` memory usage grows until it reaches size of one chunk, then it's (almost) constant.


# C-free version

For those who don't have C compiler here's version using `CompiledFunction` instead of `LibraryFunction`.
It can be compiled to virtual machine.
`weakCompositions` and `lazyWeakCompositions` functions work exactly as above, but are slower and less memory efficient.

 ClearAll[nextWeakCompositionsChunk]
 nextWeakCompositionsChunk = 
 Compile[{{n, _Integer}, {k, _Integer}, {prevComposition, _Integer, 1}, {chunkSize, _Integer}},
 Module[{lastNonZeroPosition = k, composition = prevComposition},
 While[composition[[lastNonZeroPosition]] == 0, --lastNonZeroPosition];
 Table[
 ++composition[[lastNonZeroPosition - 1]];
 If[lastNonZeroPosition == k,
 --composition[[-1]]
 ,(* else *)
 composition[[-1]] = composition[[lastNonZeroPosition]] - 1;
 composition[[lastNonZeroPosition]] = 0
 ];
 If[composition[[-1]] == 0,
 --lastNonZeroPosition
 ,(* else *)
 lastNonZeroPosition = k
 ];
 composition
 ,
 chunkSize
 ]
 ],
 RuntimeOptions -> "Speed",
 CompilationTarget -> "WVM"
 ];

 ClearAll[weakCompositions]
 weakCompositions[n_Integer?NonNegative, 0] := {}
 weakCompositions[n_Integer?NonNegative, 1] := {{n}}
 weakCompositions[n_Integer?NonNegative, k_Integer /; k >= 2] :=
 nextWeakCompositionsChunk[n, k, Join[ConstantArray[0, k - 2], {-1, n + 1}], Binomial[n + k - 1, n]]

 
 Needs["Streaming`"]

 ClearAll[lazyWeakCompositions]
 lazyWeakCompositions[n_Integer?NonNegative, 0, chunkSize : _Integer?Positive : 10^5, opts : OptionsPattern[]] := LazyListCreate[{}, chunkSize, opts]
 lazyWeakCompositions[n_Integer?NonNegative, 1, chunkSize : _Integer?Positive : 10^5, opts : OptionsPattern[]] := LazyListCreate[{{n}}, chunkSize, opts]
 lazyWeakCompositions[n_Integer?NonNegative, k_Integer /; k >= 2, chunkSize : _Integer?Positive : 10^5, opts : OptionsPattern[]] :=
 With[{length = Binomial[n + k - 1, n]},
 Module[{active = False, compositionsLeft = length, prevComposition = Join[ConstantArray[0, k - 2], {-1, n + 1}]},
 LazyListCreate[
 IteratorCreate[
 ListIterator,
 (active = True) &,
 Module[{realChunkSize = Min[chunkSize, compositionsLeft], taken},
 If[realChunkSize === 0,
 {}
 ,(* else *)
 compositionsLeft -= realChunkSize;
 taken = nextWeakCompositionsChunk[n, k, prevComposition, realChunkSize];
 prevComposition = Last[taken];
 taken
 ]
 ] &,
 TrueQ[active] &,
 Remove[active, compositionsLeft, prevComposition] &
 ],
 chunkSize,
 opts,
 "Length" -> length,
 "FiniteList" -> True
 ]
 ]
 ]