Re: random matrix from row and column sums
- To: mathgroup at smc.vnet.net
- Subject: [mg53698] Re: random matrix from row and column sums
- From: "Ray Koopman" <koopman at sfu.ca>
- Date: Mon, 24 Jan 2005 03:37:19 -0500 (EST)
- References: <csinpt$njk$1@smc.vnet.net><csntqf$4jh$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
The algorithm in my previous post is suitable only for tables with
small counts. It runs out of space very quickly as the counts grow.
Here is something that uses less space and is also faster,
faster even than using Reduce. It is built around P[n,m],
which returns a list of all the decompositions of an integer n
into the sum of m nonincreasing nonnegative integers.
In[1]:=
P[0,m_] := {Table[0,{m}]};
P[n_,m_]/;n<m := With[{z=Table[0,{m-n}]},Join[#,z]&/@P[n,n]];
P[n_,1] := {{n}};
P[n_,2] := Table[{n-i,i},{i,0,n/2}];
P[n_,m_] := ToExpression[ "Block[" <>
ToString[Table[SequenceForm["n",k],{k,m-2}]] <>
",Flatten[Table[" <>
ToString[Table[If[k==m,SequenceForm["n",m-2,"-i",m-1],
SequenceForm["i",k]],{k,m,1,-1}]] <> "," <>
StringTake[ToString[Table[Which[
k==1,StringForm["{i1,0,``/``}",n,m],
k==2,StringForm["{i2,i1,(n1=``-i1)/``}",n,m-1],
True,StringForm["{i`3`,i`2`,(n`2`=n`1`-i`2`)/`4`}",
k-2,k-1,k,m+1-k]], {k,m-1}]], {2,-2}] <> "]," <>
ToString[m-2] <> "]]" ]
In[6]:=
allfxy[fx_,fy_] := Block[{ky = Length@fy, f, fxy, gy},
fxy = Reap[ Scan[ Scan[
If[And@@MapThread[LessEqual,{#,fy}],Sow[{#}]]&, Permutations@#]&,
P[fx[[1]],ky]] ][[2,1]];
Do[f = Reap[ Scan[ Scan[
If[And@@MapThread[LessEqual,{#,fy}],Sow[#]]&, Permutations@#]&,
P[fx[[ix]],ky]] ][[2,1]];
fxy = Reap[Scan[Function[gxy, gy = Plus@@gxy;
Scan[If[And@@MapThread[LessEqual,{#+gy,fy}],Sow[Append[gxy,#]]]&,
f]], fxy]][[2,1]], {ix,2,Length@fx-1}];
Append[#,fy-Plus@@#]&/@fxy]
In[7]:=
NonnegativeIntegerMatrices[rows_, cols_] := Module[
{mat, r, c, a, vars, cond, red, m = Length[rows], n = Length[cols]},
mat = Table[a[i, j], {i, m}, {j, n}];
r = Thread[Plus @@ Transpose[mat] == rows];
c = Thread[Plus @@ mat == cols];
vars = Flatten[mat];
cond = Thread[Flatten[mat] >= 0];
red = {ToRules[Reduce[Join[r, c, cond], vars, Integers]]};
If[red != {}, mat /. red, {}]]
In[8]:=
Timing@Length[m = allfxy[{7,3,2},{2,9,1}]]
Timing@Length[mat = NonnegativeIntegerMatrices[{7,3,2},{2,9,1}]]
Sort@m === Sort@mat
Out[8]= {0.01 Second,17}
Out[9]= {0.03 Second,17}
Out[10]= True
In[11]:=
Timing@Length[m = allfxy[Range@4,Range@4]]
Timing@Length[mat = NonnegativeIntegerMatrices[Range@4,Range@4]]
Sort@m === Sort@mat
Out[11]= {0.02 Second,261}
Out[12]= {0.44 Second,261}
Out[13]= True
In[14]:=
Timing@Length[m = allfxy[Range@5,Range@5]]
Timing@Length[mat = NonnegativeIntegerMatrices[Range@5,Range@5]]
Sort@m === Sort@mat
Out[14]= {2.38 Second,22645}
Out[15]= {57.49 Second,22645}
Out[16]= True