5
$\begingroup$

Firstly I wanted to find all the non-negative integer solutions for the following equation given the value of $k_{max}$

$$ k_1 +k_2 +k_3+k_4 = k_{max} $$

For this I use the following code:-

findAll[kmax_]:=Partition[
 Flatten[Permutations /@ IntegerPartitions[kmax, {4}, Range[0, 
 kmax]]], 4]

Secondly I want to do in some sense the inverse of the above problem which is to find the position of the $k_i$ state given a set of $\{k_1,k_2,k_3,k_4\}$ and $k_{max}$.

For this part I use the code below.

findIndex[kmax_, state_] := Position[Partition[
Flatten[Permutations /@ IntegerPartitions[kmax, {4}, Range[0, kmax]]], 
4], val_ /; val == state]

findIndex[4, {0, 1, 1, 2}]
(*{{34}}*)

I really don't like this method for the second part.

Can this be achieved without going over all the permutations?

$\endgroup$

4 Answers 4

7
$\begingroup$

I think you'll find this considerably more efficient than the accepted answer:

f= With[{s = Accumulate@Reverse@# + 1, r = Range[Length@# - 1]},
         Tr[(Pochhammer[Rest@s, r] - Pochhammer[Most@s, r])/r!] + 1] &;

Example:

test={77620, 58215, 38810, 19405};

f@test // RepeatedTiming

Block[{$RecursionLimit = 100000}, index[Tr@test, Most@test] // RepeatedTiming]

{0.000017, 954814548275041}

{0.7, 954814548275041}

About 4 orders of magnitude faster.

Additionally, no need to mess with $RecursionLimit for larger cases (which even with it, can crash kernel using index, at least this happens on my machine).

$\endgroup$
5
  • $\begingroup$ ciao, yours is based on the same Sort[ ] of Hubble07's data, correct? $\endgroup$
    – MikeY
    Dec 19, 2018 at 14:15
  • $\begingroup$ @MikeY - "same" as in what? It does not use the order he generates - because that's a kind of weird order - it uses the sort of that result. Looking at your code, it's the same order used there. Is that what you are asking? If not, let me know... $\endgroup$
    – ciao
    Dec 19, 2018 at 19:13
  • $\begingroup$ That's it...you got it. $\endgroup$
    – MikeY
    Dec 19, 2018 at 19:17
  • $\begingroup$ @ciao Can you please have a look at link . It is a more general version of this solution which also includes negative integers. $\endgroup$
    – Hubble07
    Dec 28, 2018 at 18:02
  • $\begingroup$ @Hubble07 - sure - but it may be a day or two before I can devote any thinking time to it - a quick glance gives me the gut feeling that may be a very hard problem. $\endgroup$
    – ciao
    Dec 29, 2018 at 7:04
4
$\begingroup$

You can replace your homegrown findAll[ ] function with FrobeniusSolve[ ] and get a single call that creates the list in a nice ordered ranking (that you can also get by Sort[ ]-ing the results of findAll[ ])

 n=3; res=FrobeniusSolve[{1,1,1,1},n];

 (* {{0, 0, 0, 3}, {0, 0, 1, 2}, {0, 0, 2, 1}, {0, 0, 3, 0}, {0, 1, 0, 2},
     {0, 1, 1, 1}, {0, 1, 2, 0}, {0, 2, 0, 1}, {0, 2, 1, 0}, {0, 3, 0, 0},
     {1, 0, 0, 2}, {1, 0, 1, 1}, {1, 0, 2, 0}, {1, 1, 0, 1}, {1, 1, 1, 0}, 
     {1, 2, 0, 0}, {2, 0, 0, 1}, {2, 0, 1, 0}, {2, 1, 0, 0}, {3, 0, 0, 0}} *)

add an index to the list and display...

 resa = Transpose[Join[{Range[Length[res]]}, (res // Transpose)]];
 resa // TableForm

$\left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & 3 \\ 2 & 0 & 0 & 1 & 2 \\ 3 & 0 & 0 & 2 & 1 \\ 4 & 0 & 0 & 3 & 0 \\ 5 & 0 & 1 & 0 & 2 \\ 6 & 0 & 1 & 1 & 1 \\ 7 & 0 & 1 & 2 & 0 \\ 8 & 0 & 2 & 0 & 1 \\ 9 & 0 & 2 & 1 & 0 \\ 10 & 0 & 3 & 0 & 0 \\ 11 & 1 & 0 & 0 & 2 \\ 12 & 1 & 0 & 1 & 1 \\ 13 & 1 & 0 & 2 & 0 \\ 14 & 1 & 1 & 0 & 1 \\ 15 & 1 & 1 & 1 & 0 \\ 16 & 1 & 2 & 0 & 0 \\ 17 & 2 & 0 & 0 & 1 \\ 18 & 2 & 0 & 1 & 0 \\ 19 & 2 & 1 & 0 & 0 \\ 20 & 3 & 0 & 0 & 0 \\ \end{array} \right)$

You can see the well-structured order in the list. The 4th column is redundant, since its value is determined by the other three, so only work with $k1$, $k2$, and $k3$. Define the following functions, which do not require you to construct the list of results in order to get the index position:

 inc[n_, k_, d_] := Binomial[n - k + 1 + d - 1, d - 1] + inc[n, k - 1, d];
 inc[n_, k_, 1] := k;
 inc[n_, 0, d_] := 0;
 index[n_, {k1_, k2_, k3_}] := inc[n, k1, 3] + inc[n - k1, k2, 2] + inc[n, k3, 1] + 1

Now try it on a real problem...

 n = 80; res = FrobeniusSolve[{1, 1, 1, 1}, n];

Pick off a point

 resa[[91000]]
 (* {64, 6, 0, 10} *)

Work it backwards

 index[80,{64,6,0}]
 (* 91000 *)

A little bigger...

 n = 400; res = FrobeniusSolve[{1, 1, 1, 1}, n]; 
 res[[800000]]
 (* {10, 31, 258, 101} *)

 index[400,{10,31,258}]
 (* 800000 *)

The function is essentially instantaneous

 index[8*10^20 , {64, 6, 0}]
 (* 20479999999999999998468800000000000000039320 *)
$\endgroup$
1
  • $\begingroup$ Note that resource usage of index function is linear in k1 and k2. index[1022, {1022, 0, 0}] hits default recursion limit. Block[{$RecursionLimit = Infinity}, index[15000, {15000, 0, 0}]] crashes my kernel. $\endgroup$
    – jkuczm
    Dec 18, 2018 at 21:09
4
$\begingroup$

If you don't really care about specific order given by Permutations of IntegerPartitions and can accept lexical order, then there's relatively easy algorithm calculating composition index.

Let's start with function generating all weak compositions in lexicographic order:

weakCompositionsLex // ClearAll
weakCompositionsLex[n_Integer, k_Integer?NonNegative] :=
  Sort[Join @@ (Permutations /@ IntegerPartitions[n, {k}, Range[0, n]])]

Index of weak $k$-composition $c = (c_1, c_2, \ldots, c_k)$, in lexicographic order, can be calculated using following formula:

(formula = Binomial[Sum[Indexed[c, j], {j, 1, k}] + k - 1, k - 1]
  - Sum[Binomial[Sum[Indexed[c, j], {j, i + 1, k}] - 1 + k - i, k-i], {i, 1, k - 1}]
) // TraditionalForm

$$ \binom{\sum_{j=1}^k c_j + k-1}{k-1}-\sum_{i=1}^{k-1} \binom{\sum_{j=i+1}^k c_j-1 + k-i}{k-i} $$

Let's check that it gives correct results for all compositions with n and k up to 9:

tmp = Evaluate@formula & /. {Indexed[c, i_] :> #[[i]]};
Table[tmp /@ weakCompositionsLex[n, k] === Range@Binomial[n + k - 1, k - 1], {n, 0, 9}, {k, 0, 9}];
And @@ Flatten@%
(* True *)

Above formula can be implemented as following top level function:

weakCompositionLexIndexTop = If[# === {},
  1
(* else *),
  With[{acc = Accumulate@Reverse@#, range = Range[Length@# - 1]},
    Binomial[Last@acc + Length@# - 1, Length@# - 1] -
    Total@Binomial[Most@acc + range - 1, range]
  ]
]&;

Basic tests:

Table[weakCompositionLexIndexTop /@ weakCompositionsLex[n, k] === Range@Binomial[n + k - 1, k - 1], {n, 0, 10}, {k, 0, 10}];
And @@ Flatten@%
(* True *)

We can also implement it as compiled function:

weakCompositionLexIndexC = Hold@Compile[{{comp, _Integer, 1}},
  Module[{k, m, index, sum, tmpSum, jMin, jMax, bin, binPrev, l, i, iZeroSum, d, el},
    k = Length@comp;
    If[k <= 1, Return@1];
    index = binPrev = comp[[k]];
    sum = index - 1;
    i = k - 1;
    m = 1;
    If[sum === -1,
      While[comp[[i]] === 0 && i >= 2, --i];
      m += k - 1 - i;
    ];
    tmpSum = sum + comp[[i]];
    If[tmpSum === 0 && i >= 2,
      iZeroSum = i;
      While[tmpSum === 0 && i >= 2, tmpSum += comp[[--i]]];
      d = iZeroSum - i;
      m += d;
      index += d;
      sum = 0;
      binPrev = 1;
    ];
    While[i >= 2,
      ++m;
      el = comp[[i]];
      tmpSum = sum;
      sum += el;
      jMin = 2;
      If[m < sum,
        jMax = m;
        bin = l = sum + 1;
      (* else *),
        jMax = sum;
        bin = l = m + 1;
      ];
      If[el < jMax - jMin && binPrev > 0,
        jMin = tmpSum + 1;
        jMax = sum;
        l = tmpSum + m;
        bin = Quotient[binPrev l, m];
      ];
      Do[bin = Quotient[bin (++l), j], {j, jMin, jMax}];
      binPrev = bin;
      index += bin;
      --i;
    ];
    sum += First@comp + 1;
    If[sum === 0, Return@1];
    jMax = If[m < sum, m, sum];
    bin = l = sum + m;
    Do[bin = Quotient[bin (--l), j], {j, 2, jMax}];
    bin - index
  ],
  RuntimeOptions -> {"Speed", "CatchMachineIntegerOverflow" -> True, "WarningMessages" -> False},
  CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True
] /.
  Part -> Compile`GetElement //. HoldPattern[Compile`GetElement@x__ = y_] :> (Part@x = y) //.
  HoldPattern[Plus][pre___, x_, HoldPattern[Times][-1, y_], post___] :> Plus[pre, Subtract[x, y], post] //
  ReleaseHold;

Test all compositions with n and k up to 12:

Table[weakCompositionLexIndexC@weakCompositionsLex[n, k] === Range@Binomial[n + k - 1, k - 1], {n, 0, 12}, {k, 12}];
And @@ Flatten@%
(* True *)

Final function that uses compiled function when possible, and falls back to top level version on integer overflow and for symbolic arguments.

With[{cf = ReplacePart[weakCompositionLexIndexC, 7 -> weakCompositionLexIndexTop]},
  weakCompositionLexIndex // ClearAll;
  weakCompositionLexIndex@l_List := cf@l;
]

Tests:

weakCompositionLexIndex@{}
(* 1 *)

weakCompositionLexIndex@{a, b, c, d}
(* 1 + c - 1/2 (1 + c + d) (2 + c + d) + 1/2 (1 + b + c + d) (2 + b + c + d) - 1/6 (1 + b + c + d) (2 + b + c + d) (3 + b + c + d) + 1/6 (1 + a + b + c + d) (2 + a + b + c + d) (3 + a + b + c + d) *)

% /. {a -> 2, b -> 0, c -> 1, d -> 5}
(* 83 *)

weakCompositionLexIndex@{2, 0, 1, 5}
(* 83 *)

weakCompositionLexIndex@Join[{1}, ConstantArray[0, 10^6]]
(* 1000001 *)

weakCompositionLexIndex@{210, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
(* 2620460268444674457614627940 *)

weakCompositionLexIndex@{64, 6, 0, 8*10^20 - 64 - 6}
(* 20479999999999999998468800000000000000039320 *)

weakCompositionLexIndex@c
(* weakCompositionLexIndex[c] *)

Table[weakCompositionLexIndex@weakCompositionsLex[n, k] === Range@Binomial[n + k - 1, k - 1], {n, 0, 12}, {k, 12}];
And @@ Flatten@%
(* True *)

Let's check time and memory usage on data analogous to one used in answer by Henrik Schumacher.

compositions = weakCompositionsLex[400, 4];
rCompositions = RandomChoice[compositions, 100000];
weakCompositionLexIndex@rCompositions; // MaxMemoryUsed // RepeatedTiming
(* {0.021, 801872} *)

We can index all ten millions compositions and check that we get correct consecutive numbers:

res = weakCompositionLexIndex@compositions; // MaxMemoryUsed // RepeatedTiming
res === Range@Binomial[400 + 4 - 1, 4 - 1]
(* {2.36, 86621048} *)
(* True *)
$\endgroup$
2
  • 1
    $\begingroup$ This returns incorrect results, e.g. weakCompositionLexIndex@{210, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} gives 211148762779630179, it s/b 2620460268444674457614627940... $\endgroup$
    – ciao
    Dec 12, 2018 at 23:46
  • $\begingroup$ @ciao Thanks, added overflow caching. $\endgroup$
    – jkuczm
    Dec 13, 2018 at 14:07
3
$\begingroup$
kmax = 400;
partitions = Join @@ (Permutations /@ IntegerPartitions[kmax, {4}, Range[0, kmax]]);
rpartitions = RandomChoice[partitions, 100000];

If you have to do many lookups, you can create a lookup table stored in an Association. The advantage is that the average complexity of a lookup within an association should grow logarithmically with the length of the Association while the average complexity of Position grows linearly with the length of list.

createLookupTable[partitions_] := AssociationThread[partitions -> Range[Length[partitions]]];
findIndex[tab_Association, state_?MatrixQ] := Lookup[tab, state];
findIndex[tab_Association, state_?VectorQ] := Lookup[tab, {state}][[1]];

tab = createLookupTable[partitions]; // AbsoluteTiming // First
idx = findIndex[tab, rpartitions]; // AbsoluteTiming // First

12.6645

0.147966

2681434480

Alternatively, you may use a SparseArray as lookup table; it is quite different performance characteristics (it is built faster and needs considerably less memory, but lookup operations are much slower):

createLookupTable2[partitions_] := SparseArray[partitions + 1 -> Range[Length[partitions]], {1, 1, 1, 1} (1 + Max[partitions])];
findIndex2[tab_, state_] := Extract[tab, state + 1];

tab2 = createLookupTable2[partitions]; // AbsoluteTiming // First
idx2 = findIndex2[tab2, rpartitions]; // AbsoluteTiming // First

3.67315

3.80017

346480928

A third variant can be provided by Nearest. Since you data can be interpreted as points $\mathbb{R}^4$, it allows to use $k$D-trees to make quick lookups (logarithmic complexity in terms of length of the list to look up in). Moreover, the preparation time is the lowest so far and the storage requirements are as low as the SparseArray method.

createLookupTable3[partitions_] := Nearest[partitions -> Automatic];
findIndex3[tab_, state_] := Join @@ tab3[state, {1, 0}];

tab3 = createLookupTable3[partitions]; // AbsoluteTiming // First
idx3 = findIndex3[tab3, rpartitions]; // AbsoluteTiming // First
ByteCount[tab3]

2.02313

0.060386

346477184

Of course, all three methods return the same result:

idx == idx2 == idx3

True

$\endgroup$
2
  • $\begingroup$ Do you think its possible to do this without using the command Permutations. I mean the second part $\endgroup$
    – Hubble07
    Dec 10, 2018 at 17:34
  • 1
    $\begingroup$ Yes, there is certainly a way. But notice that Permutations is called precisely once for constructing the lookup table. The lookups do not need a call to Permutations. $\endgroup$ Dec 10, 2018 at 17:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.