Finally, I parse the candidate partition sets, splitting them at the points where the running "stuttering" totals equal the target coefficients, tally these, extract the results using an undocumented form of extractundocumented form of extract I discovered for speed, and applying the adaptive bin binomial coefficient extractors. Lastly, I check if the result is empty (coefficient set was invalid) and return 0 (invalid) or the total of the coefficients (valid).
replaced http://mathematica.stackexchange.com/ with https://mathematica.stackexchange.com/
getCoeff[ipow_, epow_, expset_] :=
Module[{f, fastRF, myIP, px = {{}}, target = DeleteCases[expset, 0],
set = Join @@ MapThread[ConstantArray, {ipow, epow}], nres, x, bin,
tset, breaks, accs, splits, seqs, res},
(* fugly patch for now *)
If[(Tr[ipow*epow] < Tr@target) || (Min@target < Min@ipow), Return[0]];
fastRF[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
o = Ordering[c];
x = 1 - 2 UnitStep[-1 - Length[b] + o];
x = FoldList[Max[#, 0] + #2 &, x];
x[[o]] = x;
Pick[c, x, -1]];
myIP[set_, used_, left_, n_] :=
With[{tmp = fastRF[set, used]},
(* Another fugly patch for non-conformance *)
If[Length@tmp < left, Return[{{}}]];
IntegerPartitions[n, Length@tmp - left, DeleteDuplicates@tmp]];
tset = Tally[set];
Module[{cnt = #[[2]], tmp, reset = #[[2]]},
bin[#[[1]], x_] := (tmp = Binomial[cnt, x]; cnt -= x;
If[cnt == 0, cnt = reset]; tmp)] & /@ tset;
(f[#[[1]], z_] = z <= #[[2]]) & /@ Tally[set];
nres = Module[{z = 0},
Nest[With[{zz = #}, z++;
Join @@ (Module[{mip = myIP[set, #, Length@target - z, target[[z]]], tmp = #},
mip = Join[tmp, #] & /@ mip;
Pick[mip, Apply[And, Apply[f, Tally /@ mip, {2}], {1}]]] & /@
zz)] &, px, Length@target]];
breaks = Accumulate@target;
accs = Accumulate /@ nres;
Reap[Do[
x = Join @@ Position[accs[[z]], Alternatives @@ breaks];
x = Differences[Prepend[x, 0]];
splits =
Transpose@{Span @@@ Transpose[{Accumulate[Prepend[Most@x, 0]] + 1, Accumulate[x]}]};
segs = Tally /@ Rest@Extract[nres[[z]], Prepend[splits, {{}}]];
res = Times @@ Join @@ Apply[bin, segs, {2}];
Sow[res];, {z, Length@nres}]] //
If[#[[2]] == {}, 0, Tr@#[[2, 1]]] &
]
Update 25/02/2015 per comments:
Looking at your test harness, there's an obvious issue - you can't just Quiet something and expect that errors sans messages means the same performance. Even though using Quiet does usually reduce timings when errors cause messages, there's still overhead in Mathematica handling the issue, and of course the likelihood that there's a cascade of problems presented to the code... also on large result lists, using rules can be quite inefficient for replacements. Look into Dispatch and using Replace directly with a limited level specification if you have tasks that need to do such things in the future.
Since the original question, comments, Math.SE question, etc. all implied the feeding of well-conditioned exponent lists, this was built with that in mind. I've added another ugly patch for now, but be advised, it just short-ciruits the problem, but it also means there's work, sometimes a lot, getting done where it need not be.
That said, I took your code, and just changed one line to be:
getCoeff[Sequence @@ Transpose@Tally@list1, #] & /@ list2
in your partialONk function. No other changes made, so there's extra work converting the argument format, and the extra work being done in my function because bad sets are passed to it. Nonetheless, the following is the result (netbook timings, I ran out of patience at 12):
Both grow quickly (expected from the problem), but the original coeff explodes.
There appears to be much overhead from the rest of your code, based on just timing the coefficient functions. I'd venture re-thinking it, and in particular generating only valid exponent sets, will go a long way in improving timing. As it stands, much of the boost in finding the coefficients is getting swamped in other overhead. By example, the two functions, over a range of 1 to 5 components, each with random powers between 2 and 5 for internal/external power, and a random exponent set of appropriate size. I gave up at six components - getCoeff finished in under a second, I aborted coeff after twenty minutes. At five components, getCoeff was ~2800 times faster... N.B. This had to be done with a Log scale so getCoeff would even show up on the plot:
getCoeff[ipow_, epow_, expset_] :=
Module[{f, fastRF, myIP, px = {{}}, target = DeleteCases[expset, 0],
set = Join @@ MapThread[ConstantArray, {ipow, epow}], nres, x, bin,
tset, breaks, accs, splits, seqs, res},
(* fugly patch for now *)
If[(Tr[ipow*epow] < Tr@target) || (Min@target < Min@ipow), Return[0]];
fastRF[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
o = Ordering[c];
x = 1 - 2 UnitStep[-1 - Length[b] + o];
x = FoldList[Max[#, 0] + #2 &, x];
x[[o]] = x;
Pick[c, x, -1]];
myIP[set_, used_, left_, n_] :=
With[{tmp = fastRF[set, used]},
IntegerPartitions[n, Length@tmp - left, DeleteDuplicates@tmp]];
tset = Tally[set];
Module[{cnt = #[[2]], tmp, reset = #[[2]]},
bin[#[[1]], x_] := (tmp = Binomial[cnt, x]; cnt -= x;
If[cnt == 0, cnt = reset]; tmp)] & /@ tset;
(f[#[[1]], z_] = z <= #[[2]]) & /@ Tally[set];
nres = Module[{z = 0},
Nest[With[{zz = #}, z++;
Join @@ (Module[{mip = myIP[set, #, Length@target - z, target[[z]]], tmp = #},
mip = Join[tmp, #] & /@ mip;
Pick[mip, Apply[And, Apply[f, Tally /@ mip, {2}], {1}]]] & /@
zz)] &, px, Length@target]];
breaks = Accumulate@target;
accs = Accumulate /@ nres;
Reap[Do[
x = Join @@ Position[accs[[z]], Alternatives @@ breaks];
x = Differences[Prepend[x, 0]];
splits =
Transpose@{Span @@@ Transpose[{Accumulate[Prepend[Most@x, 0]] + 1, Accumulate[x]}]};
segs = Tally /@ Rest@Extract[nres[[z]], Prepend[splits, {{}}]];
res = Times @@ Join @@ Apply[bin, segs, {2}];
Sow[res];, {z, Length@nres}]] //
If[#[[2]] == {}, 0, Tr@#[[2, 1]]] &
]
getCoeff[ipow_, epow_, expset_] :=
Module[{f, fastRF, myIP, px = {{}}, target = DeleteCases[expset, 0],
set = Join @@ MapThread[ConstantArray, {ipow, epow}], nres, x, bin,
tset, breaks, accs, splits, seqs, res},
(* fugly patch for now *)
If[(Tr[ipow*epow] < Tr@target) || (Min@target < Min@ipow), Return[0]];
fastRF[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
o = Ordering[c];
x = 1 - 2 UnitStep[-1 - Length[b] + o];
x = FoldList[Max[#, 0] + #2 &, x];
x[[o]] = x;
Pick[c, x, -1]];
myIP[set_, used_, left_, n_] :=
With[{tmp = fastRF[set, used]},
(* Another fugly patch for non-conformance *)
If[Length@tmp < left, Return[{{}}]];
IntegerPartitions[n, Length@tmp - left, DeleteDuplicates@tmp]];
tset = Tally[set];
Module[{cnt = #[[2]], tmp, reset = #[[2]]},
bin[#[[1]], x_] := (tmp = Binomial[cnt, x]; cnt -= x;
If[cnt == 0, cnt = reset]; tmp)] & /@ tset;
(f[#[[1]], z_] = z <= #[[2]]) & /@ Tally[set];
nres = Module[{z = 0},
Nest[With[{zz = #}, z++;
Join @@ (Module[{mip = myIP[set, #, Length@target - z, target[[z]]], tmp = #},
mip = Join[tmp, #] & /@ mip;
Pick[mip, Apply[And, Apply[f, Tally /@ mip, {2}], {1}]]] & /@
zz)] &, px, Length@target]];
breaks = Accumulate@target;
accs = Accumulate /@ nres;
Reap[Do[
x = Join @@ Position[accs[[z]], Alternatives @@ breaks];
x = Differences[Prepend[x, 0]];
splits =
Transpose@{Span @@@ Transpose[{Accumulate[Prepend[Most@x, 0]] + 1, Accumulate[x]}]};
segs = Tally /@ Rest@Extract[nres[[z]], Prepend[splits, {{}}]];
res = Times @@ Join @@ Apply[bin, segs, {2}];
Sow[res];, {z, Length@nres}]] //
If[#[[2]] == {}, 0, Tr@#[[2, 1]]] &
]
Update 25/02/2015 per comments:
Looking at your test harness, there's an obvious issue - you can't just Quiet something and expect that errors sans messages means the same performance. Even though using Quiet does usually reduce timings when errors cause messages, there's still overhead in Mathematica handling the issue, and of course the likelihood that there's a cascade of problems presented to the code... also on large result lists, using rules can be quite inefficient for replacements. Look into Dispatch and using Replace directly with a limited level specification if you have tasks that need to do such things in the future.
Since the original question, comments, Math.SE question, etc. all implied the feeding of well-conditioned exponent lists, this was built with that in mind. I've added another ugly patch for now, but be advised, it just short-ciruits the problem, but it also means there's work, sometimes a lot, getting done where it need not be.
That said, I took your code, and just changed one line to be:
getCoeff[Sequence @@ Transpose@Tally@list1, #] & /@ list2
in your partialONk function. No other changes made, so there's extra work converting the argument format, and the extra work being done in my function because bad sets are passed to it. Nonetheless, the following is the result (netbook timings, I ran out of patience at 12):
Both grow quickly (expected from the problem), but the original coeff explodes.
There appears to be much overhead from the rest of your code, based on just timing the coefficient functions. I'd venture re-thinking it, and in particular generating only valid exponent sets, will go a long way in improving timing. As it stands, much of the boost in finding the coefficients is getting swamped in other overhead. By example, the two functions, over a range of 1 to 5 components, each with random powers between 2 and 5 for internal/external power, and a random exponent set of appropriate size. I gave up at six components - getCoeff finished in under a second, I aborted coeff after twenty minutes. At five components, getCoeff was ~2800 times faster... N.B. This had to be done with a Log scale so getCoeff would even show up on the plot: