####################################################################### TestTupel := proc() local n, LowerLimits, UpperLimits, Result; n:=3: LowerLimits[1]:=1; LowerLimits[2]:=1; LowerLimits[3]:=1; UpperLimits[1]:=4; UpperLimits[2]:=4; UpperLimits[3]:=4; Result := Tupel(n, LowerLimits, UpperLimits); print ("Result: ",Result); end proc; Tupel := proc(n::integer,LowerLimits,UpperLimits) # The Maple program Tupel returns a list of all integer n-tupels [t[1], t[2],...,t[i],...,t[n]] # with LowerLimits[i] <= t[i] <= UpperLimits[i]. # Tupel is a Maple implementation of Algorithm M (mixed radix generation) from Donald E. Knuth. # Source: # Donald E. Knuth # The Art of Computer Programming # Volume 4, Fascile 3 # Generating All Tuples and Permutations. # With the extension that the lower limits may be different from zero. # Begonnen am: 27. Januar 2011 # Letzte Änderung am 29. Januar 2011 # thomas wieder@t-online.de # Written for Maple 13. # Usage: #> n:=3; LowerLimits[1]:=1; LowerLimits[2]:=1; LowerLimits[3]:=1; UpperLimits[1]:=4; UpperLimits[2]:=4; UpperLimits[3]:=4; #> Tupel(n, LowerLimits, UpperLimits); #"Result: ", [1, 1, 1], [1, 1, 2], [1, 1, 3], [1, 2, 1], [1, 2, 2], [1, 2, 3], # [1, 3, 1], [1, 3, 2], [1, 3, 3], [2, 1, 1], [2, 1, 2], [2, 1, 3], # [2, 2, 1], [2, 2, 2], [2, 2, 3], [2, 3, 1], [2, 3, 2], [2, 3, 3], # [3, 1, 1], [3, 1, 2], [3, 1, 3], [3, 2, 1], [3, 2, 2], [3, 2, 3], # [3, 3, 1], [3, 3, 2], [3, 3, 3] local iverbose, alu, i, j, ListOfTuples; iverbose := 0; if iverbose = 1 then print("Tupel: Input: n, LowerLimits, UpperLimits:",n, LowerLimits, UpperLimits); end if; ListOfTuples := NULL; # M1. for j from 1 to n do alu[j] := LowerLimits[j]; end do; UpperLimits[0] := 2; do #M2. ListOfTuples := ListOfTuples, [seq(alu[i],i=1..n)]; if iverbose = 1 then print("Tupel: ListOfTuples = ",ListOfTuples); end if; #M3. j := n; #M4. do if alu[j] = UpperLimits[j] - 1 then alu[j] := LowerLimits[j]; j := j - 1; else break end if; end do; #M5. if j = 0 or LowerLimits[j] >= UpperLimits[j] then break end if; alu[j] := alu[j]+1; end do; return ListOfTuples; end proc; ####################################################################### # NEXT PROGRAM ####################################################################### frtotx2 := proc (s) # Calculate the Cartesian product among the lists of a list of lists. # Maple 13 source code. # Author: xavier, see http://www.mapleprimes.com/users/xavier # and http://www.mapleprimes.com/questions/41503-Combinations-Of-Numbers-In-Sets. # With modifications by Thomas Wieder. # Begonnen am: 16.7.2010 # Letzte Änderung: 16.7.2010 # thomas.wieder@t-online.de #> frtotx2([[1, 2, 3], [a, b, c], [A, B]]); # [1, a, A], [1, a, B], [1, b, A], [1, b, B], [1, c, A], [1, c, B], [2, a, A], # [2, a, B], [2, b, A], [2, b, B], [2, c, A], [2, c, B], [3, a, A], [3, a, B], # [3, b, A], [3, b, B], [3, c, A], [3, c, B] local iverbose,liste,i,j,ki,Result; iverbose:=0; liste := [seq([s[1][i]], i = 1 .. nops(s[1]))]; if iverbose=1 then print("liste:",liste) fi; for ki from 2 to nops(s) do Result:=NULL; for i from 1 to nops(liste) do for j from 1 to nops(s[ki]) do Result := Result,[op(liste[i]), s[ki][j]]; if iverbose=1 then print("Result:",Result) fi; # End of do loop j. end do; # End of do loop i. end do; liste:=[Result]; # End of do loop ki. end do; return Result end proc; ####################################################################### # NEXT PROGRAM ####################################################################### KP:=proc(L::list) # Begonnen am: 28.8.2007 # Letzte Änderung am: 28.8.2007 # thomas.wieder@t-online.de # Berechne das kartesische Produkt einer Liste von Listen. # Usage: # > L1 := [a, b, c]; L2 := [1, 2, 3]; L3 := [A, B, C]; L := [L1, L2, L3]; # [[a, b, c], [1, 2, 3], [A, B, C]] # > KP(L); # [[a, 1, A], [a, 1, B], [a, 1, C], [a, 2, A], [a, 2, B], [a, 2, C], [a, 3, A], # [a, 3, B], [a, 3, C], [b, 1, A], [b, 1, B], [b, 1, C], [b, 2, A], [b, 2, B], # [b, 2, C], [b, 3, A], [b, 3, B], [b, 3, C], [c, 1, A], [c, 1, B], [c, 1, C], # [c, 2, A], [c, 2, B], [c, 2, C], [c, 3, A], [c, 3, B], [c, 3, C]] local l, KP, ZahlElemente; global Resultat; if nops(L)=1 then Resultat:=L; end if; ZahlElemente:=nops(L); LP(L[1],L[2]); KP:=Resultat; for l from 3 to ZahlElemente do LP(KP,L[l]); KP:=Resultat; end do; Resultat:=KP; end proc; LP:=proc(L1,L2) # Begonnen am: 28.8.2007 # Letzte Änderung am: 28.8.2007 # Berechne das einfache Produkt zweier Listen. local iverbose, l1, l2, L; global Resultat; iverbose:=0; L:=NULL; if not type(L1,integer) and not type(L1,list) then print("Argument not an integer or a list:",L1); stop; end if; if not type(L2,integer) and not type(L2,list) then print("Argument not an integer or a list:",L2); stop; end if; for l1 in L1 do for l2 in L2 do L:=L,[op(l1),l2]; end do; end do; L:=[L]; Resultat:=L; end proc; ####################################################################### # NEXT PROGRAM ####################################################################### simplex_lattice_point_next := proc( n::integer, t::integer, more:: boolean, x_in::list) # John Burkardt's Fortran program simplex_lattice_point_next # translated to Maple programming language by Thomas Wieder. # start: 8.2.2011 # last change: 9.2.2011 # thomas.wieder@t-online.de # Source: # John Burkardt # http://people.sc.fsu.edu/~jburkardt/f77_src/asa299/asa299.f # Last revised on 10 January 2007. # # Usage with Maple worksheet: #> read "C:/Users/Wieder/Documents/EDV/Maple/Latgen.mpl"; #> Result := 'Result'; Output := NULL; #> n := 3; t := 4; Result[1] := false; Result[2] := [seq(0, i = 1 .. n)]; #> do Result := simplex_lattice_point_next(n, t, Result[1], Result[2]); Output := Output, Result[2]; if Result[1] = false then break end if end do; Output; #[4, 0, 0], [3, 1, 0], [3, 0, 1], [2, 2, 0], [2, 1, 1], [2, 0, 2], [1, 3, 0], [1, 2, 1], [1, 1, 2], [1, 0, 3], [0, 4, 0], [0, 3, 1], [0, 2, 2], [0, 1, 3], [0, 0, 4] # #********************************************************************* # # SIMPLEX_LATTICE_POINT_NEXT generates lattice points in a simplex. # # Discussion: # # The simplex is defined by N-dimensional points X such that: # # 0 <= X(1:N) # # and # # sum ( X(1:N) ) <= T # # where T is an integer. # # Lattice points are points X which satisfy the simplex conditions and # for which all the components are integers. # # This routine generates all the lattice points in a given simplex, one at # a time, in a reverse lexicographic order. # # To use the routine, initialize by setting N and T to appropriate values, # and MORE to FALSE. The initial value of X is not important. # # Call the routine. On return, X will contain the first lattice point in # the simplex. If MORE is TRUE, then the routine may be called again to # get the next point. In fact, as long as the output value of MORE is # TRUE, there is at least one more lattice point that can be found by # making another call. When MORE is returned as FALSE, then there are no # more lattice points; the value of X returned at that time is the # "last" such point. # # During the computation of a sequence of lattice points, the user should # not change the values of N, T, MORE or X. # # The output for N = 3, T = 4 would be: # # 1 4 0 0 # 2 3 1 0 # 3 3 0 1 # 4 2 2 0 # 5 2 1 1 # 6 2 0 2 # 7 1 3 0 # 8 1 2 1 # 9 1 1 2 # 10 1 0 3 # 11 0 4 0 # 12 0 3 1 # 13 0 2 2 # 14 0 1 3 # 15 0 0 4 # # The number of lattice points is (T+N-1)! / ( T! * (N-1)! ) # # Thus, for N = 3, T = 4, we have # # number = 6! / ( 4! * 2! ) = 720 / ( 24 * 2 ) = 15 # # Modified: # # 10 January 2007 # # Author: # # John Burkardt # # Reference: # # Scott Chasalow, Richard Brand, # Algorithm AS 299: # Generation of Simplex Lattice Points, # Applied Statistics, # Volume 44, Number 4, 1995, pages 534-545. # # Albert Nijenhuis, Herbert Wilf, # Combinatorial Algorithms for Computers and Calculators, # Second Edition, # Academic Press, 1978, # ISBN: 0-12-519260-6, # LC: QA164.N54. # # Parameters: # # Input, integer N, the spatial dimension. # N must be positive. # # Input, integer T, the characteristic of the simplex. # T must be nonnegative. # # Input/output, logical MORE, initialized to FALSE by the user to # begin a sequence of calculations, returned by the routine as TRUE, # if there are more values of X that can be calculated, or FALSE # if the accompanying value of X is the last one for this sequence. # # Input/output, integer X(N), not initialized by the user, but not # changed by the user on subsequent calls. The routine returns # a new point on each call, and on subsequent calls uses the input # value (old point) to compute the output value (next point). # # implicit none # # integer n # # integer i # integer j # integer k # logical more # integer t # integer x(n) # local iverbose, i::integer, j::integer, k::integer, x::list, more_out::boolean; iverbose := 0; for i from 1 to n do x[i]:=x_in[i]; end do; more_out := true; if ( more = false) then if iverbose=1 then print("Latgen: more=false!"); print(seq(x[i], i = 1 .. n)); end if; if ( n < 1 ) then print (" SLPN: "); print (" SLPN: SIMPLEX_LATTICE_POINT_NEXT - Fatal error!" ); print (" SLPN: N < 1." ); more_out := false; return [more_out, [seq(x[i],i=1..n)]] end if; if ( t < 0 ) then print (" SLPN: " ); print (" SLPN: SIMPLEX_LATTICE_POINT_NEXT - Fatal error!" ); print (" SLPN: T < 0." ); more_out := false; return [more_out, [seq(x[i],i=1..n)]] end if; more_out := true; j := 1; x[1] := t; for i from 2 to n do x[i] := 0; end do; # # The first point can actually also be the last! # if ( n = 1 ) then more_out := false; end if; else if iverbose=1 then print("Latgen: more=true!"); print(seq(x[i], i = 1 .. n)); end if; # # Search X(N-1 down to 1) for the first nonzero element. # If none, then terminate. (This should not happen.) # Otherwise, set J to this index. # Decrement X(J) by 1. # Set X(J+1:N) to (T-X(1:J),0,0,...0). # j := n; for i from n-1 to 1 by -1 do if ( 0 < x[i] ) then j := i; x[j] := x[j] - 1; x[j+1] := t; for k from 1 to j do x[j+1] := x[j+1] - x[k]; end do; for k from j+2 to n do x[k] := 0; end do; if ( x[n] = t ) then more_out := false; end if; return [more_out, [seq(x[i],i=1..n)]] end if; end do; if ( j = n ) then print (" SLPN: " ); print (" SLPN: SIMPLEX_LATTICE_POINT_NEXT - Fatal error!" ); print (" SLPN: The input X vector is nonpositive in" ); print (" SLPN: all entries except possibly the last" ); print (" SLPN: one." ), print (" SLPN: " ); print (" SLPN: Perhaps the user has miscalled the" ); print (" SLPN: routine or altered data between calls." ); print (" SLPN: " ); print (" SLPN: ABNORMAL TERMINATION." ); more := false; return more_out end if; end if; if iverbose=1 then print("Latgen: End!"); print(seq(x[i], i = 1 .. n)); end if; return [more_out, [seq(x[i],i=1..n)]]; end proc; ####################################################################### # NEXT PROGRAM ####################################################################### enum := proc( r::integer, c::integer, n::(list(integer)), m::(list(integer)) ) # Ian Saunders' Fortran program enum.f as modified by John Burkardt # translated to Maple's programming language by Thomas Wieder. # start: 10.2.2011 # last change: 31.3.2011 # thomas.wieder@t-online.de # Source: # John Burkardt # http://people.sc.fsu.edu/~jburkardt/f77_src/asa205/asa205.f # http://people.sc.fsu.edu/~jburkardt/f77_src/asa205/asa205.html # Modified: 29 November 2006. # See the source code below for more information on the authors. # # Usage with Maple worksheet: # #> read "C:/Users/Wieder/Documents/EDV/Maple/enum.mpl"; #> ttt := enum(2, 3, [1, 3], [1, 1, 2]); # " " # " Table" # 1, 0, 0 # 0, 1, 2 # " " # " Table" # 0, 0, 1 # 1, 1, 1 # " " # "EVAL summary" # " Number of cases (ignoring multiplicity):", 2 # " Number of cases (allowing multiplicity):", 3 # " Probability sum = ", 1.000000000 #ttt := [Matrix(2, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 2}), # Matrix(2, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (2, 1) = 1, (2, 2) = 1, (2, 3) = 1})] #********************************************************************* # # ENUM generates contingency tables with given shape and row and column sums. # # Discussion: # # The routine enumerates all M by N contingency tables with given row # and column totals, and calculates the hypergeometric probability of # each table. # # For tables having two or more row sums repeated, equivalent # tables differing only by a row permutation are not separately # enumerated. A representative of each equivalence class is enumerated # and the multiplicity of each class calculated. # <-- No, I have commented out this feature, # --> All solutions will be displyed. # --> See comment lines marked by the date 31.3.2011, Thomas Wieder: # # For each table enumerated, subroutine EVAL is called to carry out # calculations on the table. # # Note that the entries in the column sum and row sum vectors will # be (implicitly) sorted into ascending order, and results will be # returned as though these orderings were being used! # # Modified: # # 29 November 2006 # # Author: # # Ian Saunders # Modifications by John Burkardt # # Reference: # # Ian Saunders, # Algorithm AS 205: # Enumeration of R x C Tables with Repeated Row Totals, # Journal of the Royal Statistical Society: Series C (Applied Statistics), # Volume 33, Number 3, London, 1984, pages 340-352. # # Parameters: # # Input, integer R, the number of rows. # # Input, integer C, the number of columns. # # Input, integer N(R), the row sums. # # Input, integer M(C), the column sums. # # Input, external EVAL, the name of the user supplied routine # that will be called each time a new table is determined. # The routine has the form: # # subroutine eval ( iflag, table, r, c, n, m, prob, mult ) # # See the dummy version of EVAL for an explanation of the meaning # of the arguments. # # Output, integer IFAULT, error flag. # 0, no error occurred. # # Local parameters: # # table(I,j) - (I,J)-th entry of current table # # NTOTAL - total of table entries # # bound(I,J) - current upper bound on table(i,j) to satisfy # row and column totals # # rept(I) - LOGICAL = true if row totals N(I), N(I-1) are equal # = false otherwise # # reps(I) - number of previous rows equal to row I # # MAX - maximum dimension of TABLE # # NMAX - maximum number of observations + 1 # # mult2(I) - maximum number of equivalent tables given first i rows # # Z(J) - lower bound on sum of entries used by algorithm C # # PROB2(I,J) - partial sum of terms in log(P) # # implicit none # # integer c # integer mmax # parameter ( mmax = 10 ) # integer nmax # parameter ( nmax = 201 ) # integer r # # integer bound(mmax,mmax) # integer c2 # integer cm # external eval # double precision factlm(nmax) # double precision flogm(nmax) # integer i # logical ieqim # integer ifault # integer iflag # integer ii # integer iim # integer iip # integer ij # integer j # integer jj # integer jjm # integer jkeep # integer jm # integer jnext # integer k # integer keep # integer left # integer m(c) # integer m2(mmax) # integer maxrc # integer mult # integer mult2(mmax) # integer multc # integer multr # integer mult0 # integer n(r) # integer n2(mmax) # integer ntot # integer ntotal # double precision prob # double precision prob0 # double precision prob2(mmax,mmax) # integer r2 # integer reps(mmax) # integer repsc # integer repsr # logical rept(mmax) # logical reptc(mmax) # integer rm # integer rowbnd # integer rowsum # integer table(mmax,mmax) # integer table2(mmax,mmax) # integer z(mmax) local iverbose, mmax, nmax, bound::Array, c2, cm, factlm::list, flogm::list, i, ieqim, ifault, iflag, ii, iim, iip, ij, j, jj, jjm, jkeep, jm, jnext, k, keep, left, m2::list, maxrc, mult, mult2::list, multc, multr, mult0, n2::list, ntot, ntotal, prob, prob0, prob2::Array, r2, reps::list, repsc, repsr, rept::list, reptc::(list,boolean), rm, rowbnd, rowsum, table1::Array, table2::Array, z::list, goto_label::string, Result::list, ResultList::list; global label_123, label_200, label_210, label_220, label_230, label_300, label_330, lable360, label370, label_380, label_393, label_400, label_500; iverbose := 0; # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # The parameters nmax and mmax limit the size of the problem! # You can try to increase them if you run into trouble with enum. # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! mmax := 200; nmax := 200; Result := [0,0,0.0]; ResultList := NULL; # # Check the input values. # ifault := 0; # print(r,c,n,m); # if (r = c) and (r = 1) then # ResultList := n[1]; # return [ResultList]; # end if; if ( mmax < r ) then ifault := 1; error "enum: mmax < r, ifault =", ifault end if; if ( mmax < c ) then ifault := 2; error "enum: mmax < c ,ifault =", ifault end if; if ( r <= 0 ) then ifault := 3; error "enum: r <= 0, ifault =", ifault end if; if ( c <= 0 ) then ifault := 4; error "enum: c <= 0, ifault =", ifault end if; if ( r <> nops(n) ) then ifault := 5; error "enum: r <> rowdim(n), ifault =", ifault end if; if ( c <> nops(m) ) then ifault := 6; error "enum: c <> coldim(m), ifault =", ifault end if; ntotal := 0; for i from 1 to r do if ( n[i] < 0 ) then ifault := 7; error "enum: n[i] <= 0!,i,n[i]:",i,n[i], "ifault =", ifault end if; ntotal := ntotal + n[i]; end do; ntot := 0; for j from 1 to c do if ( m[j] < 0 ) then ifault := 8; error "enum: m[j] <= 0! j, m[j]:",j,m[j], "ifault =", ifault end if; ntot := ntot + m[j]; end do; if ( ntot <> ntotal ) then ifault := 9; error "enum: ntot <> ntotal, ifault =", ifault end if; if ( nmax <= ntotal ) then ifault := 10; error "enum: nmax <= ntotal, ifault =", ifault end if; # Addition by Thomas Wieder to cover the case that a vector is given instead of a matrix. if r = 1 then ResultList := m; return [ResultList]; elif c = 1 then ResultList := n; return [ResultList]; end if; # # We need to copy the input parameters. # One reason is that we may end up "transposing" the matrix, # so we need to use arrays that are large enough to work with # either the row or column dimension. # c2 := c; r2 := r; for i from 1 to r2 do n2[i] := n[i]; end do; for j from 1 to c2 do m2[j] := m[j]; end do; # # The first call to EVAL is simply to allow EVAL to initialize. # iflag := 1; prob := 0.0; mult := 0; for j from 1 to c do for i from 1 to r do table1[i,j] := 0; end do; end do; Result := eval_enum ( iflag, table1, r, c, n, m, prob, mult, Result ); # # Initialize FLOGM(K) = LOG(K-1), FACTLM(K) = LOG(K-1 FACTORIAL) # flogm[1] := 0.0; factlm[1] := 0.0; for k from 1 to ntotal do flogm[k+1] := evalf(ln( k )); factlm[k+1] := factlm[k] + flogm[k+1]; end do; # # Constants. # rm := r2 - 1; cm := c2 - 1; # # Sort rows and columns into ascending order. # 13.3.2011: # Sorting of rows and columns commented out by Thomas Wieder. # The desired consequence is that rows and columns are not permuted # in comparison to the initial order of the marginals. # Try for example (with and without sorting) # >ttt := NULL; ttt := enum(3, 3, [1, 1, 4], [1, 2, 3]); # rtable(1 .. 3, 1 .. 3, ttt[1]); # # for i from 1 to r2 - 1 do # for ii from i+1 to r2 do # if ( n2[ii] < n2[i] ) then # keep := n2[i]; # n2[i] := n2[ii]; # n2[ii] := keep; # end if; # end do; # end do; # # for j from 1 to c2 - 1 do # for jj from j+1 to c2 do # if ( m2[jj] < m2[j] ) then # keep := m2[j]; # m2[j] := m2[jj]; # m2[jj] := keep; # end if; # end do; # end do; # Calculate multiplicities of rows and columns. # # REPTC(J) = TRUE if columns J and J-1 have the same total. # REPT(I) = TRUE if rows I and I-1 have the same total. # multc := 1; repsc := 1; reptc[1] := false; for j from 2 to c2 do # 31.3.2011: Commented out by Thomas Wieder # --> All permutated solutions will be displayed. # reptc[j] := is( m2[j] = m2[j-1] ); reptc[j] := false; if ( reptc[j] ) then repsc := repsc + 1; multc := multc * repsc; else repsc := 1; end if; end do; multr := 1; repsr := 1; rept[1] := false; for i from 2 to r2 do # 31.3.2011: Commented out by Thomas Wieder # --> All permutated solutions will be displayed. # rept[i] := is( n2[i] = n2[i-1] ); rept[i] := false; if ( rept[i] ) then repsr := repsr + 1; multr := multr * repsr; else repsr := 1; end if; end do; # # If column multiplicity exceeds row multiplicity, transpose the table. # if ( multr < multc ) then maxrc := max ( r2, c2 ); for ij from 1 to maxrc do keep := n2[ij]; n2[ij] := m2[ij]; m2[ij] := keep; end do; keep := r2; r2 := c2; c2 := keep; rm := r2 - 1; cm := c2 - 1; for i from 1 to r2 do rept[i] := reptc[i]; end do; multr := multc; end if; # # Set up the initial table. # # Maximum multiplicity. # mult2[1] := multr; reps[1] := 1; # # Constant term in probability. # prob0 := -factlm[ntotal+1]; for i from 1 to r2 do ii := n2[i]; prob0 := prob0 + factlm[ii+1]; end do; for j from 1 to c2 do jj := m2[j]; prob0 := prob0 + factlm[jj+1]; end do; # # Calculate bounds on row 1. # for j from 1 to c2 do bound[1,j] := m2[j]; end do; # # For each I, find the greatest I-TH row satisfying the bounds. # for i from 1 to r2 do if ( i <> 1 ) then prob0 := prob2[i-1,c2]; end if; left := n2[i]; # # Elements of row I. # ieqim := rept[i]; for j from 1 to cm do ij := min ( left, bound[i,j] ); table2[i,j] := ij; if ( j = 1 ) then prob2[i,j] := prob0 - factlm[ij+1]; else prob2[i,j] := prob2[i,j-1] - factlm[ij+1]; end if; left := left - table2[i,j]; if ( i < r2 ) then bound[i+1,j] := bound[i,j] - table2[i,j]; end if; if ( left = 0 ) then for jj from j+1 to c2 do table2[i,jj] := 0; prob2[i,jj] := prob2[i,jj-1]; bound[i+1,jj] := bound[i,jj]; end do; #Remarks on the following command goto(label_123): #1.) Maple does not allow a goto which breaks a loop. #2.) The jump to label label_123 causes a disturbance in the output. #3.) Fortunately, a simple break instaed of the goto seems to work well. #goto(label_123) break; end if; if ( ieqim ) then ieqim := is(table2[i,j] = table2[i-1,j]); end if; end do; table2[i,c2] := left; prob2[i,c2] := prob2[i,cm] - factlm[left+1]; if ( i < r2 ) then bound[i+1,c2] := bound[i,c2] - left; end if; label_123: if ( i <> 1 ) then mult2[i] := mult2[i-1]; reps[i] := 1; if ( ieqim ) then reps[i] := reps[i-1] + 1; mult2[i] := mult2[i] / reps[i]; end if; end if; end do; # # Call EVAL for the first table. # iflag := 2; prob := prob2[r2,c2]; mult := mult2[r2]; if ( r = r2 and c = c2 ) then for j from 1 to c do for i from 1 to r do table1[i,j] := table2[i,j]; end do; end do; else for j from 1 to c do for i from 1 to r do table1[i,j] := table2[j,i]; end do; end do; end if; # The command # ResultList := ResultList,table1; # does not provide the tables table1 on return from enum. # The tables table1 get lost on return, just the name "table1" is returned. ResultList := eval(ResultList),rtable(1..r,1..c,map(eval,table1)); #print("enum test 1:",table1); Result := eval_enum ( iflag, table1, r, c, n, m, prob, mult, Result ); # # Commence enumeration of remaining tables. # label_200: i := r2; label_210: i := i - 1; # # If I = 0, no more tables are possible. # if ( i = 0 ) then iflag := 3; prob := 0.0; mult := 0; for j from 1 to c do for i from 1 to r do table1[i,j] := 0; end do; end do; Result := eval_enum ( iflag, table1, r, c, n, m, prob, mult, Result ); if ifault <> 0 then return "Error return from enum: ifault =", ifault else # Return to the calling program: #print("enum test 2:",ResultList); return [ResultList]; fi; end if; j := cm; left := table2[i,c2]; rowbnd := bound[i,c2]; # # Try to decrease element (I,J). # label_220: #print(r,c,n,m,table2[i,j],left,rowbnd); if ( 0 < table2[i,j] and left < rowbnd ) then if iverbose=1 then print("goto(label_230)") fi; goto(label_230) end if; # # Element (I,J) cannot be decreased. Try (I,J-1). # if ( j = 1 ) then if iverbose=1 then print("goto(label_210)") fi; goto(label_210) end if; left := left + table2[i,j]; rowbnd := rowbnd + bound[i,j]; j := j - 1; if iverbose=1 then print("goto(label_220)") fi; goto(label_220); # # Decrease element (I,J). # label_230: ij := table2[i,j]; prob2[i,j] := prob2[i,j] + flogm[ij+1]; table2[i,j] := table2[i,j] - 1; bound[i+1,j] := bound[i+1,j] + 1; # # If row I was the same as row I-1, it is no longer. # if ( reps[i] <> 1 ) then reps[i] := 1; mult2[i] := mult2[i-1]; end if; # # Complete row I with the largest possible values. # ii := i; iip := ii + 1; iim := ii - 1; jnext := j + 1; left := left + 1; if iverbose=1 then print("goto(label_380)") fi; goto(label_380); # # Fill up remaining rows. # label_300: ii := ii + 1; # # The last row is treated separately. # if ( ii = r2 ) then if iverbose=1 then print("goto(label_400)") fi; goto(label_400) end if; iip := ii + 1; iim := ii - 1; # # Row total N(II) is not a repeat. Make row II as large as possible. # if is( not rept[ii] ) then left := n2[ii]; jnext := 1; goto(label_380) end if; # # Repeated row totals. # # (I) If row II-1 satisfies the bounds on row II, repeat it. # for j from 1 to c2 do if ( bound[ii,j] < table2[iim,j] ) then goto(label_330) end if; ij := table2[iim,j]; table2[ii,j] := ij; bound[iip,j] := bound[ii,j] - table2[ii,j]; if ( j = 1 ) then prob2[ii,j] := prob2[iim,c2] - factlm[ij+1]; else prob2[ii,j] := prob2[ii,j-1] - factlm[ij+1]; end if; end do; # # Row II is a repeat of row II-1. # reps[ii] := reps[iim] + 1; mult2[ii] := mult2[iim] / reps[ii]; if iverbose=1 then print("goto(label_300)") fi; goto(label_300); # # Element J of row II-1 was too big. # # Construct the sequence Z(J) of lower bounds. # label_330: # # If J = 1, the bounds are satisfied automatically. # if ( j = 1 ) then ij := bound[ii,1]; table2[ii,1] := ij; prob2[ii,1] := prob2[iim,c2] - factlm[ij+1]; jnext := 2; left := n2[ii] - table2[ii,1]; bound[iip,1] := 0; if iverbose=1 then print("goto(label_380)") fi; goto(label_380) end if; z[j] := n2[ii]; jm := j - 1; for jj from j+1 to c2 do z[j] := z[j] - bound[ii,jj]; end do; for jjm from 1 to jm do jj := j - jjm; z[jj] := z[jj+1] - bound[ii,jj+1]; end do; # # (II) If the cumulative totals of row II-1 all exceed the bounds Z(J), # make element (II,J) equal to its bound. # rowsum := 0; jkeep := 0; for jj from 1 to jm do rowsum := rowsum + table2[iim,jj]; if ( rowsum < z[jj] ) then if iverbose=1 then print("goto(label_360)") fi; goto(label_360) end if; if ( z[jj] < rowsum and 0 < table2[iim,jj] ) then jkeep := jj; end if; end do; table2[ii,j] := bound[ii,j]; bound[iip,j] := 0; ij := table2[ii,j]; prob2[ii,j] := prob2[ii,jm] - factlm[ij+1]; reps[ii] := 1; mult2[ii] := mult2[iim]; # # Complete row II with the largest possible elements. # jnext := j + 1; left := n2[ii]; for jj from 1 to j do left := left - table2[ii,jj]; end do; goto(label_380); # # (III) The cumulative sums violate the bounds. # If no element of row II-1 can be changed to satisfy the bounds, # no suitable row II is possible. # In that case, go back and try decreasing row II-1. # label_360: if ( jkeep = 0 ) then i := ii; if iverbose=1 then print("goto(label_210)") fi; goto(label_210) end if; # # Element (II,JKEEP) can be decreased. # label_370: bound[iip,jkeep] := bound[iip,jkeep] + 1; ij := table2[ii,jkeep]; prob2[ii,jkeep] := prob2[ii,jkeep] + flogm[ij+1]; table2[ii,jkeep] := table2[ii,jkeep] - 1; # # Complete the row. # jnext := jkeep + 1; left := n2[ii]; for jj from 1 to jkeep do left := left - table2[ii,jj]; end do; # # Row II is complete up to element JNEXT-1. # Make the remaining elements as large as possible. # (This section of code is used for every row, repeated or not.) # label_380: for j from jnext to cm do table2[ii,j] := min ( left, bound[ii,j] ); left := left - table2[ii,j]; bound[iip,j] := bound[ii,j] - table2[ii,j]; ij := table2[ii,j]; if ( j = 1 ) then prob2[ii,j] := prob2[iim,c2] - factlm[ij+1]; else prob2[ii,j] := prob2[ii,j-1] - factlm[ij+1]; end if; if ( left = 0 ) then for jj from j+1 to c2 do table2[ii,jj] := 0; prob2[ii,jj] := prob2[ii,jj-1]; bound[iip,jj] := bound[ii,jj]; end do; if iverbose=1 then print("goto(label_393)") fi; goto(label_393) end if; end do; table2[ii,c2] := left; prob2[ii,c2] := prob2[ii,cm] - factlm[left+1]; bound[iip,c2] := bound[ii,c2] - left; label_393: reps[ii] := 1; if ( 1 < ii ) then mult2[ii] := mult2[iim]; end if; if iverbose=1 then print("goto(label_300)") fi; goto(label_300); # # The final row. # label_400: if ( not rept[r2] ) then # # Not a repeat. Set row R equal to its bounds. # ij := bound[r2,1]; table2[r2,1] := ij; prob2[r2,1] := prob2[rm,c2] - factlm[ij+1]; for j from 2 to c2 do ij := bound[r2,j]; table2[r2,j] := ij; prob2[r2,j] := prob2[r2,j-1] - factlm[ij+1]; end do; mult2[r2] := mult2[rm]; if iverbose=1 then print("goto(label_500)") fi; goto(label_500) end if; # # Row total R is a repeat. Ensure that it is less than row R-1. # for j from 1 to c2 do # # If row R would be bigger than row R-1 go back and try # decreasing row R-2. # if ( table2[rm,j] < bound[r2,j] ) then i := rm; if iverbose=1 then print("goto(label_210)") fi; goto(label_210) end if; ij := bound[r2,j]; table2[r2,j] := ij; if ( j = 1 ) then prob2[r2,j] := prob2[rm,c2] - factlm[ij+1]; else prob2[r2,j] := prob2[r2,j-1] - factlm[ij+1]; end if; # # Row R is already less than row R-1, so no more checks are needed. # if ( table2[r2,j] <> table2[rm,j] ) then for jj from j+1 to c2 do ij := bound[r2,jj]; table2[r2,jj] := ij; prob2[r2,jj] := prob2[r2,jj-1] - factlm[ij+1]; end do; mult2[r2] := mult2[rm]; if iverbose=1 then print("goto(label_500)") fi; goto(label_500) end if; end do; # # Row R is a repeat of row R-1. # reps[r2] := reps[rm] + 1; mult2[r2] := mult2[rm] / reps[r2]; # # The table is complete. # label_500: iflag := 2; prob := prob2[r2,c2]; mult := mult2[r2]; if (( r = r2) and (c = c2 )) then for j from 1 to c do for i from 1 to r do table1[i,j] := table2[i,j]; end do; end do; else for j from 1 to c do for i from 1 to r do table1[i,j] := table2[j,i]; end do; end do; end if; # The command # ResultList := ResultList,table1; # does not provide the tables table1 on return from enum. # The tables table1 get lost on return, just the name "table1" is returned. ResultList := eval(ResultList),rtable(1..r,1..c,map(eval,table1)); #print("enum test 3:",table1); Result := eval_enum ( iflag, table1, r, c, n, m, prob, mult, Result ); if iverbose=1 then print("goto(label_200)") fi; goto(label_200) end proc; eval_enum := proc( iflag, table1::table, m, n, rowsum::list, colsum::list, prob, mult, count_in::list ) #********************************************************************* # # EVAL is called by ENUM every time a new contingency table is determined. # # Discussion: # # This is a dummy version of the routine. # # The user might wish to print out each contingency table, or collect # some statistics. # # Modified: # # 27 November 2006 # # Author: # # John Burkardt # # Reference: # # Ian Saunders, # Algorithm AS 205, # Enumeration of R x C Tables with Repeated Row Totals, # Applied Statistics, # Volume 33, Number 3, pages 340-352, 1984. # # Parameters: # # Input, integer IFLAG, input flag. # 1, this is the first call. No table is input. # 2, this is a call with a new table. # 3, this is the last call. No table is input. # # Input, integer TABLE(M,N), the current contingency table. # # Input, integer M, the number of rows. # # Input, integer N, the number of columns. # # Input, integer ROWSUM(M), the row sums. # # Input, integer COLSUM(N), the column sums. # # Input, double precision PROB, the logarithm of the hypergeometric probability # of this table. # # Input, integer MULT, the multiplicity of this table, that is, # the number of different tables that still have the same set of # entries, but differ by a permutation of some rows and columns. # # implicit none # # integer mmax # parameter ( mmax = 10 ) # # integer m # integer n # # integer colsum(n) # integer count1 # integer count2 # integer iflag # integer mult # double precision prob # double precision psum # integer rowsum(m) # integer table(mmax,mmax) # save count1 # save count2 # save psum # # data count1 / 0 / # data count2 / 0 / # data psum / 0.0D+00 / local mmax, count1, count2, psum, do_print; # 10.3.2011 do_print := 0; mmax := 10; count1 := op(1,count_in); count2 := op(2,count_in); psum := op(3,count_in); # # First call, no table, initialize. # if ( iflag = 1 ) then count1 := 0; count2 := 0; psum := 0.0; # # Call with a new table. # elif ( iflag = 2 ) then count1 := count1 + 1; count2 := count2 + mult; psum := psum + mult * evalf(exp( prob )); if do_print = 1 then i4mat_print ( m, n, table1, " Table" ); fi; # # Last call, no table. # elif ( iflag = 3 ) then if do_print = 1 then print(" "); print ( "EVAL summary" ); print ( " Number of cases (ignoring multiplicity):", count1); print ( " Number of cases (allowing multiplicity):", count2); print ( " Probability sum = ", psum); fi; end if; return [count1, count2, psum] end proc; i4mat_print := proc( m, n, a::table, title::string ) #********************************************************************* # # I4MAT_PRINT prints an I4MAT. # # Discussion: # # An I4MAT is a rectangular array of integer values. # # Modified: # # 30 June 2003 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows in A. # # Input, integer N, the number of columns in A. # # Input, integer A(M,N), the matrix to be printed. # # Input, character*(*) TITLE, a title to be printed first. # TITLE may be blank. # # implicit none # # integer mmax # parameter ( mmax = 10 ) # # integer m # integer n # # integer a(mmax,mmax) # integer ihi # integer ilo # integer jhi # integer jlo # character*(*) title local mmax, ihi, ilo, jhi, jlo; mmax := 10; ilo := 1; ihi := m; jlo := 1; jhi := n; i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ); return end proc; i4mat_print_some := proc( m, n, a::table, ilo, jlo, ihi, jhi, title::string ) #********************************************************************* # # I4MAT_PRINT_SOME prints some of an I4MAT. # # Discussion: # # An I4MAT is a rectangular array of integer values. # # Modified: # # 04 November 2003 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, N, the number of rows and columns. # # Input, integer A(M,N), an M by N matrix to be printed. # # Input, integer ILO, JLO, the first row and column to print. # # Input, integer IHI, JHI, the last row and column to print. # # Input, character*(*) TITLE, an optional title. # # implicit none # # integer mmax # parameter ( mmax = 10 ) # # integer incx # parameter ( incx = 10 ) # integer m # integer n # # integer a(mmax,mmax) # character*(8) ctemp(incx) # integer i # integer i2hi # integer i2lo # integer ihi # integer ilo # integer inc # integer j # integer j2 # integer j2hi # integer j2lo # integer jhi # integer jlo # integer s_len_trim # character*(*) title local mmax, ctemp::(list,string), i, i2hi, i2lo, inc, incx, j, j2, j2hi, j2lo; mmax := 10; incx := 10; if ( 0 < length ( title ) ) then print(" "); print ( title ); end if; for j2lo from max ( jlo, 1 ) to min ( jhi, n ) by incx do j2hi := j2lo + incx - 1; j2hi := min ( j2hi, n ); j2hi := min ( j2hi, jhi ); inc := j2hi + 1 - j2lo; for j from j2lo to j2hi do j2 := j + 1 - j2lo; ctemp[j2] := j; end do; #print(" ", seq(ctemp[j], j=1..inc) ); i2lo := max ( ilo, 1 ); i2hi := min ( ihi, m ); for i from i2lo to i2hi do #print ( i," ", seq(a[i,j2lo - 1 + j2],j2=1..inc) ); print ( seq(a[i,j2lo - 1 + j2],j2=1..inc) ); end do; end do; return end proc; ####################################################################### # NEXT PROGRAM ####################################################################### TestLatgen := proc() local i::integer, n::integer, L::list, U::list, Result::list; n := 3; L := NULL; U := NULL; Result := NULL; for i from 1 to n do L := L,0; U := U,n; end do; L := [L]; U := [U]; print("L, U, n:",L, U, n); Result := Latgen(L, U, n); return(Result); end proc; Latgen := proc(L, U, n) # Begonnen am: 8.2.2011 # Letzte Änderung am: 20.2.2011 # thomas.wieder@t-online.de # # Latgen calculates all n-tuples with lower bounds L and upper bounds U. # # Source: # Derek O´Connor # A Mimimum-change Algorithm for Generating Lattice Points # University College, Dublin, Ireland # June 2006. # # Written for Maple 13. # This subroutine generates a sequence of n-tuples # (x(1), x(2),..., x(n)) with l(i) <= x(i) <= u(i) # and successive n-tuples differ by one component. # Usage with a Maple worksheet: #> n := 3; L := NULL; U := NULL; Result := NULL; #> for i to n do L := L, 0; U := U, n end do; L := [L]; U := [U]; # [0, 0, 0] # [3, 3, 3] #> Result := Latgen(L, U, n); #[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [3, 1, 0], [2, 1, 0], [1, 1, 0], # [0, 1, 0], [0, 2, 0], [1, 2, 0], [2, 2, 0], [3, 2, 0], [3, 3, 0], [2, 3, 0], # [1, 3, 0], [0, 3, 0], [0, 3, 1], [1, 3, 1], [2, 3, 1], [3, 3, 1], [3, 2, 1], # [2, 2, 1], [1, 2, 1], [0, 2, 1], [0, 1, 1], [1, 1, 1], [2, 1, 1], [3, 1, 1], # [3, 0, 1], [2, 0, 1], [1, 0, 1], [0, 0, 1], [0, 0, 2], [1, 0, 2], [2, 0, 2], # [3, 0, 2], [3, 1, 2], [2, 1, 2], [1, 1, 2], [0, 1, 2], [0, 2, 2], [1, 2, 2], # [2, 2, 2], [3, 2, 2], [3, 3, 2], [2, 3, 2], [1, 3, 2], [0, 3, 2], [0, 3, 3], # [1, 3, 3], [2, 3, 3], [3, 3, 3], [3, 2, 3], [2, 2, 3], [1, 2, 3], [0, 2, 3], # [0, 1, 3], [1, 1, 3], [2, 1, 3], [3, 1, 3], [3, 0, 3], [2, 0, 3], [1, 0, 3], # [0, 0, 3]] #> nops(Result); # 64 local iverbose, i, x::list, sign::list, Lattice::list; global label_20, label_30, label_40; # Initialization iverbose := 0; Lattice := NULL; Lattice := Lattice,[seq(L[j],j=1..n)]; for i from 1 to n do x[i] := L[i]; sign[i] := +1; end do; # start the algorithm goto(label_40); label_20: sign[i] := -sign[i]; i := i+1; if (i = n+1) then return [Lattice]; fi; label_30: if (x[i] = L[i]) and (sign[i] = -1) then goto(label_20) fi; if (x[i] = U[i]) and (sign[i] = +1) then goto(label_20) fi; x[i] := x[i] + sign[i]; Lattice := Lattice,[seq(x[j],j=1..n)]; label_40: if iverbose = 1 then print ("Latgen:",Lattice); fi; i := 1; goto(label_30): end proc;