# ---------------------------------------------------------------------
# CONTINGENCY code for Maple (Version of 02-26-07)
# Based on the paper: Counting Magic squares in quasi-polynomial time
#                     by Alexander Barvinok, Alexander Samorodnitsky and
#                        Alexander Yong
# Report errors to Alexander Yong (ayong@math.umn.edu). 
# Maple compatibility: Tested for Maple 7,10.
# ---------------------------------------------------------------------

# -----------------------------------------------------------------------------------------
# QUICK INSTRUCTIONS:
#
# Purpose: This software provides an implementation of the ``simplified''
#          estimator of the number of _magic squares_ with parameters
#          n and t, i.e., positive integer n x n matrices whose row and
#          column sums equal t. (Those familiar with magic squares from
#          recreational math will be familiar with the  extra assumption 
#          that the diagonals sum to the ``magic total'' t. Our terminology 
#          uses the convention of usage in modern combinatorics, following,
#          e.g., MacMahon and Stanley.)
#
#          In fact, this software applies the same techniques to more
#          contingency tables with more general margins.
#
# Usage:   Start maple and load contingency.txt. Type ``contingency();''
#          to run the user interface. The row and column margins are
#          inputted as a vector, e.g. [3,5,7,2,10]. In the case of
#          magic squares, use, e.g., ``same_array(15, 10)'' to create
#          a length 10 vector of 15's: [15,15,15,15,...,15] 
# 
# Example 1: Suppose you want to compute the number of 6x6 matrices
#            with line sum t=36. We do this by expressing the integral
#            as S_1\times\prod_{i=1}^{m-1}S_{i+1}/S_i (see section 3.6
#            of the paper), where t_1=8. Each S_i,i>1 is computed by 
#            hit and run Markov chain Monte Carlo. This is done at
#            each step with three samples derived after waiting for
#            the random walk to ``mix'' 7000 steps (thus a total of 
#            21000 steps at each stage). The integral S_1 is computed by a 
#            sample mean using 100000 samples. Since the double stochasticizaion
#            used in the code is approximate, we choose an error of 0.001.
#
# Then Running contingency(); we answer the
#            questions, as follows:
#
#            ``Number of rows: 6;''
#            ``Number of columns: 6;''
#            ``Row margins, e.g., [3,3,...,3]: same_array(36,6)''
#            ``Column margins, e.g., [3,3,...,3]: same_array(36,6)''
#            ``Number of telescoping steps: 20''
#            ``Vector of hit and run steps: same_array(21000,20);''
#            ``Vector of hit and run steps to wait between samples: same_array(7000,20);''
#            ``Row margins in the final step: same_array(8,6);''
#            ``Column margins in the final step: same_array(8,6);''         
#            ``The error used in the sigma calculation. Example 0.001: 0.001;"
#       
#            The code outputs sum intermediate data (the margins used in the
#            telescoping steps and log of the integral $S_i=\int_{\Delta}\Psi_{i}d\mu$).
#            This is useful for parallel processing.
#            The estimated number of tables is given as [NUMBER_OF_TABLES, 0.2205723997 10^28]
#            in the case of one run we did.
#
# Example 2: Alternatively, one can run Example 1 noninteractively via:
#            markov_telescope_July06(6,6, same_array(36,6),
#                                    same_array(36,6), 20, same_array(21000,20),
#                                    same_array(7000,20), 100000, same_array(8,6), same_array(8,6),
#                                    0.00001,0);
#            (The last entry can be set to 1 if the user wants certain intermediate computations
#            outputted.)
#
# Example 3: In the above example, we estimated the number of magic squares. However, the
#            code also works to give an answer for contingency tables (although we have
#            no proven results in general, this was used to come up with conjectures in
#            the paper). This explains, e.g., why we use same_array(36,6) to describe
#            the line sum t=36. So for example, we could compute the number of 4x4
#            contingency tables with row margins R=[220,215,93,64] and column margins
#            C=[108,286,71,127] by
#            markov_telescope_July06(4,4, [220,215,93,64], [108,286,71,127], 50,
#                                    [op(same_array(20000,50))], [op(same_array(20000,50))],
#                                     100000, same_array(10,4), same_array(10,4), 0.001,0);
#
#            So here, we only used one sample point in each of the intermediate telescoping steps
#            but sampled 100000 points in the final sample mean. This led to a fairly accurate
#            estimate of the correct answer. Note that we telescope to the 4x4 t=10 magic square
#            case. Our experience is that having the final step be a magic squares case with
#            t moderately larger than n leads to better estimates.
# -----------------------------------------------------------------------------------------------
   
# Computes contingency tables via the hit and run markov chain
# Monte Carlo method
# INPUT:   mxn are the dimensions of the matrix, row_marginsin ans column_marginsin
#          are self-explanatory, num_steps is the number of the telescoping ratios
#          to estimate. num_samples is the total number of sample points per step
#          num_iterations_before sample is the length of any random walk before
#          actually taking a point to average with. Note both of these two
#          previous inputs are VECTORS (hence the length and waiting time of walks
#          can vary from step to step. final_steps is the number of
#          naive MC samples to compute the final ratio. rowend_marginsin and
#          colend_marginsin are the margins of the tables that we compute with
#          naive MC at the end. error1 is used in Db_stochize... and plotdata is
#          0 or 1 depending on if the user wants to see shape of f on a chosen
#          line during the hit and run process.
# RETURNS: [log(integral), integral, #tables] estimations

`markov_telescope_July06`:=proc(m,n, row_marginsin, column_marginsin, num_steps, num_samples, 
                                num_iterations_before_sample, final_steps, rowend_marginsin, 
                                colend_marginsin, error1, plotdata)

	local aa,ii,jj, row_descending_marginslist, column_descending_marginslist, 
              row_dropvector, column_dropvector, NN, ans, row_margins, column_margins, 
              numtables, log_factorial_product, rowend_margins, colend_margins;

	row_margins:=[]:
	column_margins:=[]:
	rowend_margins:=[]:
	colend_margins:=[]:
	
	for ii from 1 to nops(row_marginsin) do
	    row_margins:=[op(row_margins), op(ii, row_marginsin)]:
	od:
	for jj from 1 to nops(column_marginsin) do
	    column_margins:=[op(column_margins), op(jj, column_marginsin)]:
	od:
	for ii from 1 to nops(rowend_marginsin) do
	    rowend_margins:=[op(rowend_margins), op(ii, rowend_marginsin)]:
	od:
	for jj from 1 to nops(colend_marginsin) do
	    colend_margins:=[op(colend_margins), op(jj, colend_marginsin)]:
	od:

	row_descending_marginslist:=[row_margins]:
	column_descending_marginslist:=[column_margins]:
	row_dropvector:=[]:
	column_dropvector:=[]:
	NN:=sumarrayelts(row_margins, nops(row_margins)):

	# two symmetric cases, depending on which side of the rectangle is shorter
	if(m<=n) then
		 for ii from 1 to m do
		     row_dropvector:=[op(row_dropvector),
                                      (op(ii,row_margins)-op(ii,rowend_margins))/num_steps]:
		 od:	 	 		 
		 for jj from 1 to n do
		     if(jj<=m) then
		          column_dropvector:=[op(column_dropvector), 
                                              (op(jj, column_margins)-op(jj,colend_margins))/num_steps]:
		     else
		          column_dropvector:=[op(column_dropvector), (op(jj,column_margins))/num_steps]:
		     fi:
		 od:

		 #print(ROW_DROP_VECTOR, row_dropvector, COLUMN_DROP_VECTOR, column_dropvector);

		 for aa from 1 to num_steps do
		     for ii from 1 to m do
			 row_margins[ii]:=row_margins[ii]-op(ii,row_dropvector):
		     od:    
		     row_descending_marginslist:=[op(row_descending_marginslist), row_margins]:
		   

		     for jj from 1 to n do
			 column_margins[jj]:=column_margins[jj]-op(jj,column_dropvector):
		     od:
		     column_descending_marginslist:=[op(column_descending_marginslist), column_margins]:
		 od:
	else
		 for jj from 1 to n do
		     column_dropvector:=[op(column_dropvector),
                                         (op(jj,column_margins)-op(jj,colend_margins))/num_steps]:
		 od:	 	 		 
		 for ii from 1 to m do
		     if(ii<=n) then
		          row_dropvector:=[op(row_dropvector), 
                                           (op(ii, row_margins)-op(ii,rowend_margins))/num_steps]:
		     else
		          row_dropvector:=[op(row_dropvector), (op(ii, row_margins))/num_steps]:
		     fi:
		 od:

		 #print(ROW_DROP_VECTOR, row_dropvector, COLUMN_DROP_VECTOR, column_dropvector);

		 for aa from 1 to num_steps do
		     for jj from 1 to n do
			 column_margins[jj]:=column_margins[jj]-op(jj,column_dropvector):
		     od:    
		     column_descending_marginslist:=[op(column_descending_marginslist), 
                                                                          column_margins]:

		     for ii from 1 to m do
			 row_margins[ii]:=row_margins[ii]-op(ii,row_dropvector):
		     od:
		     row_descending_marginslist:=[op(row_descending_marginslist), row_margins]:
		 od:
	fi:
	print([row_descending_marginslist, column_descending_marginslist]);

	ans:=0;
	for ii from 1 to nops(row_descending_marginslist)-1 do
	    ans:=ans+(hitandrun_July06(m,n, op(ii, row_descending_marginslist), 
                                       op(ii, column_descending_marginslist), 
                                       op(ii+1, row_descending_marginslist), 
                                       op(ii+1, column_descending_marginslist),  
                                       op(ii,num_samples), op(ii,num_iterations_before_sample), 
                                       error1, plotdata));
	    print(ii, ans);
	od:

    	ans:=ans+(op(1, genzoom_sigma(m,n, op(nops(row_descending_marginslist), 
                                      row_descending_marginslist),
                                      op(nops(column_descending_marginslist), 
                                         column_descending_marginslist), 
                                      final_steps, error1, 1, 0))):

	# factors for the number of tables
	NN:=sumarrayelts(row_marginsin, m):
	log_factorial_product:=0:
	for ii from 1 to m do
	    log_factorial_product:=log_factorial_product+log(op(ii, row_marginsin)!*1.0):
	od:
	for ii from 1 to n do
	    log_factorial_product:=log_factorial_product+log(op(ii, column_marginsin)!*1.0):
	od:

	numtables:=evalf(log((NN*1.0)!))-NN*log(NN*1.0)+evalf(log((NN*1.0+m*n-1.0)!))
                                   -evalf(log((m*n-1.0)!*1.0))-log_factorial_product:
	numtables:=numtables+ans;
	RETURN([NUMBER_OF_TABLES, exp(numtables)]);
end:

# hit and run
`hitandrun_July06`:=proc(m,n, row_margins, column_margins, row_marginsprimed, column_marginsprimed, 
                         num_samples, num_iterations_before_sample, error1, plotdata)

	local aa, bb, ii, jj, tt, pos, gamma_matrix, dividing_total, random_list, gamma_list,
              ratioslist, mu_matrix, min_mu, t_plus, t_minus, tau_minus, tau_plus, stopflag, 
              max_value, tset, fset, tau_left, tau_right, f_tau_minus, f_tau_left, f_tau_right,
              f_tau_plus, pairsets, new_max_value, p_total, pset, ptotal, randnum, currentp,
              choosej, taufinal, ans, L_value, actualratioslist, pick, newstartpoint,
              pick_forplotdata, ansset, stored_minus, stored_plus, f_tau_minus_temp,
              f_tau_plus_temp, f_tau_left_temp, f_tau_right_temp, plotset, vvv;

	ratioslist:=[];   	
	newstartpoint:=1;  
	pick_forplotdata:=1:
	pick:=0:

	for tt from 1 to num_samples do
	   if(newstartpoint=1) then
		# perform one sample of GAMMA: sample initial GAMMA at random from
		# uniform distribution on the simplex, and improve by iterating hit and run
		# sample mn independent random variables
		gamma_matrix:=array(1..m, 1..n):
		random_list:=[stats[random, exponential[1,0]](m*n)]:	
		pos:=1:
		for aa from 1 to m do
		    for bb from 1 to n do
			gamma_matrix[aa,bb]:=op(pos,random_list):
			pos:=pos+1:
		    od:
		od:

		# now warm up,..., essentially, the zoom-in transform with \beta=1
		dividing_total:=0;
		for aa from 1 to m do
		    for bb from 1 to n do
			gamma_matrix[aa,bb]:=gamma_matrix[aa,bb]*op(aa,row_margins)*op(bb, column_margins):
			dividing_total:=dividing_total+gamma_matrix[aa,bb]:
		    od:
		od:
		for aa from 1 to m do
		    for bb from 1 to n do
			gamma_matrix[aa,bb]:=gamma_matrix[aa,bb]/dividing_total:
		    od:
		od:
		newstartpoint:=0:	
	   fi:

	   pick:=pick+1: # counter that counts up to num_iterations_before_sample

	   # To make one iteration of hit and run, we choose a random line \ell through
	   # the current point \gamma\in\Delta, and improve \gamma\to better \gamma

	   # mu_matrix defines the direction to hit and run from gamma_matrix
	   mu_matrix:=array(1..m, 1..n):
	   dividing_total:=0:
	   random_list:=[stats[random, normald[0,1]](m*n)]:	
	   pos:=1:
	   for aa from 1 to m do
	      for bb from 1 to n do
		mu_matrix[aa,bb]:=op(pos,random_list):
		dividing_total:=dividing_total+mu_matrix[aa,bb]:
		pos:=pos+1:
	      od:
	   od:
	
	   # balance
	   for aa from 1 to m do
	      for bb from 1 to n do
		mu_matrix[aa,bb]:=mu_matrix[aa,bb]-(1/(m*n*1.0))*dividing_total:
	      od:
	   od:

   	   # now we sample a point from \Delta\cap\ell	
	   # First, figure out the intersection \Delta\cap\ell
	   min_mu:=infinity;
	   for aa from 1 to m do
	      for bb from 1 to n do
		if(mu_matrix[aa,bb]<0) then
			if((gamma_matrix[aa,bb]/(-mu_matrix[aa,bb]))<min_mu) then
				min_mu:=gamma_matrix[aa,bb]/(-mu_matrix[aa,bb])
			fi:
		fi:		       
	      od:
	   od:	
	   t_plus:=min_mu:
	   
	   min_mu:=infinity;
	   for aa from 1 to m do
	      for bb from 1 to n do
		if(mu_matrix[aa,bb]>0) then
			if(gamma_matrix[aa,bb]/mu_matrix[aa,bb]<min_mu) then
				min_mu:=gamma_matrix[aa,bb]/mu_matrix[aa,bb]:
			fi:
		fi:		       
	      od:
	   od:	
	   t_minus:=-min_mu:

	   # next we are going to find the maximum of \sigma on the interval
	   # \Delta\cap\ell
	   # f(t)=ln\sigma(\gamma+t\mu)

	   tau_minus:=t_minus;
	   tau_plus:=t_plus;
	   stopflag:=0; 

	   max_value:=infinity;
	   tset:=[tau_minus,tau_plus];
	   
	   ansset:=[sigma_hitrun_July06(m,n, gamma_matrix, mu_matrix, tau_minus, row_marginsprimed,
                                      column_marginsprimed, same_array(1,m), same_array(1,n), error1), 
                    sigma_hitrun_July06(m,n, gamma_matrix, mu_matrix, tau_plus, row_marginsprimed,
                                      column_marginsprimed, same_array(1,m), same_array(1,n), error1)]:

	   fset:=[op(1,op(1,ansset)), op(1,op(2,ansset))]:

	   # These two variables hold the info for tau_minus and tau_plus, for use
	   # in the next round of computation, thus making a bit more savings
	   stored_minus:=op(1,ansset):
	   stored_plus:=op(2,ansset):

	   while(stopflag=0) do

	      tau_left:=tau_minus+((tau_plus-tau_minus)/3):
	      tau_right:=tau_plus-((tau_plus-tau_minus)/3):
	      f_tau_minus_temp:=stored_minus:
	      f_tau_plus_temp:=stored_plus:
	      
	      f_tau_left_temp:=sigma_hitrun_July06(m,n, gamma_matrix, mu_matrix, tau_left, 
                                      row_marginsprimed, column_marginsprimed, 
                                      op(2,f_tau_minus_temp), op(3,f_tau_minus_temp),
                                      error1):
	      f_tau_right_temp:=sigma_hitrun_July06(m,n, gamma_matrix, mu_matrix, tau_right, 
                                      row_marginsprimed, column_marginsprimed, 
                                      op(2,f_tau_plus_temp), op(3,f_tau_plus_temp), 
                                      error1):	    	
	      f_tau_minus:=op(1,f_tau_minus_temp):
	      f_tau_left:=op(1,f_tau_left_temp):
	      f_tau_right:=op(1,f_tau_right_temp):
	      f_tau_plus:=op(1,f_tau_plus_temp):

	      pairsets:=genzoom_list_sorter(tset, fset, tau_left, f_tau_left, infinity);
	      tset:=op(1,pairsets);
	      fset:=op(2,pairsets);	    

	      pairsets:=genzoom_list_sorter(tset, fset, tau_right, f_tau_right, infinity);
	      tset:=op(1,pairsets);
	      fset:=op(2,pairsets);	    

	      if((f_tau_minus<=f_tau_left) and (f_tau_left<=f_tau_right) and 
                 (f_tau_right<=f_tau_plus)) then
	           tau_minus:=tau_right:
		   stored_minus:=f_tau_right_temp:
	      else if((f_tau_minus<=f_tau_left) and (f_tau_left<=f_tau_right) and
                     (f_tau_right>=f_tau_plus)) then
		        tau_minus:=tau_left:
		        stored_minus:=f_tau_left_temp:
		   else if((f_tau_minus<=f_tau_left) and (f_tau_left>=f_tau_right)) then
			     tau_plus:=tau_right:
			     stored_plus:=f_tau_right_temp:
		        else if(f_tau_minus>=f_tau_left) then
                                tau_plus:=tau_left:
			        stored_plus:=f_tau_left_temp:
			else
			        print(THISSHOULDNTHAPPEN):	
		               fi:
                        fi:
		   fi:
              fi:
		  
	      new_max_value:=max(f_tau_minus, f_tau_left, f_tau_right, f_tau_plus):
	      if(abs(exp(max_value) - exp(new_max_value))<0.000000001) then
		  stopflag:=1;
	      fi:
	      max_value:=new_max_value:
	  od:

	  # give data to plot tset and fset
	  if(plotdata=1) then
		plotset:=[]:
	        for vvv from 1 to nops(tset) do
	           plotset:=[op(plotset), [op(vvv,tset), op(vvv,fset)]]:
	        od:
	        print(PLOTSET, plotset);
	  fi:
	
	  # Finally, let us sample a new improved \gamma...
	  # so at this point tset should give t_minus=t_0<t_1<...<t_k=t_plus
          # and fset should give the corresponding values
	  p_total:=0;
	  pset:=[]:
	  for ii from 1 to nops(tset)-1 do
	      if(op(ii,fset)=op(ii+1,fset)) then
		pset:=[op(pset), (op(ii+1,tset)-op(ii,tset))*exp(op(ii+1,fset))]:
	      else
	        pset:=[op(pset), ((op(ii+1,tset)-op(ii,tset))/
                                 (op(ii+1,fset)-op(ii,fset)))*(exp(op(ii+1,fset))-exp(op(ii,fset)))]:
	      fi:
	  od:

	  ptotal:=0;
	  # normalize pset so the sum is 1
	  for ii from 1 to nops(pset) do
	      ptotal:=ptotal+op(ii,pset):
	  od:
	  for ii from 1 to nops(pset) do
	      pset[ii]:=pset[ii]/ptotal:
	  od:

	  # now pick an integer j\in\{1,...,k\} with probability op(jj, pset)
	  randnum:=stats[random,uniform[0,1]](1):
	  currentp:=0;
	  ii:=1;
	  stopflag:=0;

	  while((stopflag=0) and (ii<=nops(pset))) do
	    currentp:=currentp+op(ii,pset):
	    if(randnum<=currentp) then
		choosej:=ii;
		stopflag:=1;		
	    fi:
	    ii:=ii+1;
	  od:	

	  randnum:=stats[random,uniform[0,1]](1):
	  if(op(choosej+1,fset)>op(choosej,fset)) then
              taufinal:=op(choosej,tset)+((op(choosej+1,tset)-op(choosej,tset))/
                                          (op(choosej+1,fset)-op(choosej,fset)))*
                                           log(randnum*(exp(op(choosej+1,fset)-op(choosej,fset))-1)+1):
	  fi:
	  if(op(choosej+1,fset)<op(choosej,fset)) then
              taufinal:=op(choosej+1,tset)+((op(choosej,tset)-op(choosej+1,tset))/
                                            (op(choosej,fset)-op(choosej+1,fset)))*
                                             log(randnum*(exp(op(choosej,fset)-op(choosej+1,fset))-1)+1):
	  fi:	
	  if(op(choosej+1,fset)=op(choosej,fset)) then
              taufinal:=op(choosej,tset)+randnum*(op(choosej+1,tset)-op(choosej,tset)):
	  fi:		
	
	  if(plotdata=1) then
		print(POINTSELECTEDWAS, taufinal, SIGMAVALUE, op(1,sigma_hitrun_July06(m,n, gamma_matrix, 
                                                                                  mu_matrix, taufinal, 
                                                                                  row_marginsprimed, 
                                                                                  column_marginsprimed,
										  same_array(1,m), 
										  same_array(1,n), 
                                                                                  error1)));
		if(pick_forplotdata=num_iterations_before_sample) then
			print(ABOVEPOINTUSED);
			pick_forplotdata:=1;
		else
			pick_forplotdata:=pick_forplotdata+1;
		fi:
	  fi:

	  #finally, update gamma_matrix
	  for ii from 1 to m do
	      for jj from 1 to n do
		gamma_matrix[ii,jj]:=gamma_matrix[ii,jj]+taufinal*mu_matrix[ii,jj]:
	      od:
	  od:
	
	  if(pick=num_iterations_before_sample) then
	     ratioslist:=[op(ratioslist),op(2,Db_stochize_adjust_numerics(gamma_matrix, row_margins, 
                          column_margins, error1))-op(2,Db_stochize_adjust_numerics(gamma_matrix, 
                                                              row_marginsprimed, column_marginsprimed, 
                                                              error1))];   
	     if(plotdata=1) then
		print(ratioslist):
	     fi:
	     pick:=0:
	     newstartpoint:=1:
	  fi:  
        od:

	print(ratioslist);
	L_value:=-infinity;
	for ii from 1 to nops(ratioslist) do
	    L_value:=max(L_value, op(ii, ratioslist)):
	od:
	
	ans:=0;
	for ii from 1 to nops(ratioslist) do
	    ans:=ans+exp(op(ii, ratioslist)-L_value):
	od:
	ans:=evalf(log((ans/nops(ratioslist))*1.0))+L_value:
	RETURN(ans);
end:

# computes sigma
`sigma_hitrun_July06`:=proc(m, n, gamma_matrix, mu_matrix, tau, row_margins, column_margins, 
mu_array, lambda_array, error1)

	local ii, jj, the_matrix, ans, correctlnsigmaA, correctmu, correctlambda;

	the_matrix:=array(1..m,1..n):
	for ii from 1 to m do
	    for jj from 1 to n do
		the_matrix[ii,jj]:=gamma_matrix[ii,jj]+tau*mu_matrix[ii,jj]:
	    od:
	od:

	# Now adjust the_matrix
	for ii from 1 to m do
	    for jj from 1 to n do
		the_matrix[ii,jj]:=the_matrix[ii,jj]/(op(ii,mu_array)*op(jj, lambda_array));
	    od:
	od:	

	ans:=Db_stochize_adjust_numerics_July06(the_matrix, row_margins, column_margins, error1):
	correctlnsigmaA:=0;
	for ii from 1 to m do
	    correctlnsigmaA:=correctlnsigmaA+op(ii,row_margins)*evalf(log(op(ii,mu_array)*1.0)):
	od:
	for jj from 1 to n do
	    correctlnsigmaA:=correctlnsigmaA
                                     +op(jj,column_margins)*evalf(log(op(jj,lambda_array)*1.0)):
	od:
	correctlnsigmaA:=correctlnsigmaA+op(2,ans):	

	correctmu:=[]:
	correctlambda:=[]:

	for ii from 1 to m do
	    correctmu:=[op(correctmu), op(ii,op(3,ans))*op(ii,mu_array)]:
	od:
	for jj from 1 to n do
	    correctlambda:=[op(correctlambda), op(jj,op(4,ans))*op(jj,lambda_array)]:
	od:

	RETURN([correctlnsigmaA, correctmu, correctlambda]);
end:

# Double stochasticize a matrix
`Db_stochize_adjust_numerics_July06`:=proc(AA, row_margins, col_margins, error1)
	local flag, ii, jj, current_total, numrows, numcols, log_sigma_A, max_diff,
	      maxerror, diag, A, mu_array, lambda_array;

	A:=matrix(AA):
	numrows:=linalg[rowdim](A):
	numcols:=linalg[coldim](A):

	mu_array:=[]:     # row scalars
	lambda_array:=[]: # column scalars

	for ii from 1 to numrows do
	    mu_array:=[op(mu_array),1]:
	od:
	
	for jj from 1 to numcols do
	    lambda_array:=[op(lambda_array),1]:
	od:
 
	flag:=0:
	while(flag=0) do 		     
		      for ii from 1 to numrows do
			  current_total:=0:
			  for jj from 1 to numcols do
			      current_total:=1.0*current_total+A[ii,jj]:
			  od:
			  for jj from 1 to numcols do
			      A[ii,jj]:=(A[ii,jj]/current_total)*op(ii,row_margins):
			  od:
			  mu_array[ii]:=mu_array[ii]*(current_total/op(ii,row_margins)):
		      od:

		     
		      for jj from 1 to numcols do
			  current_total:=0:
			  for ii from 1 to numrows do
			      current_total:=1.0*current_total+A[ii,jj]:
			  od:
			  for ii from 1 to numrows do
			      A[ii,jj]:=(A[ii,jj]/current_total)*op(jj, col_margins):
			  od:
			  lambda_array[jj]:=lambda_array[jj]*(current_total/op(jj, col_margins)):
		      od:

		      
		      # determine if we are within the error bounds
		      maxerror:=-1;
		      for ii from 1 to numrows do
			  current_total:=0:
			  for jj from 1 to numcols do
			      current_total:=1.0*current_total+A[ii,jj]:
			  od:
			  maxerror:=max(maxerror, abs(current_total-op(ii,row_margins))):
		      od:

		      for jj from 1 to numcols do
			  current_total:=0:
			  for ii from 1 to numrows do
			      current_total:=1.0*current_total+A[ii,jj]:
			  od:
			  maxerror:=max(maxerror, abs(current_total-op(jj, col_margins))):
		      od:
		      if(maxerror<error1) then
			  flag:=1:
		      fi:
	od:

	log_sigma_A:=0:
	for ii from 1 to numrows do
	    log_sigma_A:=log_sigma_A+(row_margins[ii])*log((mu_array[ii]*row_margins[ii])):
	od:
	
	for jj from 1 to numcols do
	    log_sigma_A:=log_sigma_A+(col_margins[jj])*log(lambda_array[jj]*col_margins[jj]):
	od:	
	RETURN([A, log_sigma_A, mu_array, lambda_array]):
end:

# Estimate the number of contingency tables using the sigma estimator
`genzoom_sigma`:=proc(m,n, row_margins, column_margins, num_samples, error1, beta, verbosity)
	local ii, jj, aa, bb, rr, current_A, entry_total, monte_carlo_sum, NN, log_factorial_product,
              lower_ans, upper_ans, t_val, ss, tt, new_avg, L_value, sigma_list, temp_val,
	      log_sum, beta_sum, non_scaled_new_avg, Rmax_matrix, max_log_sum, gamma_log_sum,
              newupper_ans;


	# some precalculations used in the averaging process
	NN:=sumarrayelts(row_margins, m):
	log_factorial_product:=0:
	for ii from 1 to m do
	    log_factorial_product:=log_factorial_product+log(op(ii, row_margins)!*1.0):
	od:
	for ii from 1 to n do
	    log_factorial_product:=log_factorial_product+log(op(ii, column_margins)!*1.0):
	od:
	t_val:=t_value(m,n, row_margins, column_margins);

	Rmax_matrix:=array(1..m,1..n):
	for ii from 1 to m do
	    for jj from 1 to n do
		Rmax_matrix[ii,jj]:=row_margins[ii]*column_margins[jj]:
	    od:
	od:

	sigma_list:=[]:	
	L_value:=-infinity:
	# get the create_A_matrix and rescale
	for ii from 1 to num_samples do
	    current_A:=create_A_matrix_simplex_beta_tricky_060220(m,n, beta, Rmax_matrix):
	    # compute log(sigma(A))
	    temp_val:=op(2,Db_stochize_adjust_numerics(current_A, row_margins, column_margins, error1)):
	    # scale the temp_value
	    max_log_sum:=0;
	    gamma_log_sum:=0;
	    beta_sum:=0;
	    for aa from 1 to m do
		for bb from 1 to n do
		    max_log_sum:=max_log_sum+evalf(log(Rmax_matrix[aa,bb]*1.0));
		    gamma_log_sum:=gamma_log_sum+evalf(log(current_A[aa,bb]*1.0));
		    beta_sum:=beta_sum+((Rmax_matrix[aa,bb]*1.0)^(-1/beta))*((current_A[aa,bb]*1.0)^(1/beta));
		od:
	    od:
            temp_val:=temp_val+(m*n-1)*evalf(log(beta*1.0))+(1/beta)*max_log_sum+(1-(1/beta))*gamma_log_sum+(m*n)*evalf(log(beta_sum*1.0));
	    sigma_list:=[op(sigma_list), temp_val]:
            L_value:=max(L_value,temp_val):
	    
	    # report on intermediate data if the user desires it
	    if(verbosity=1) then
		new_avg:=0;
		for rr from 1 to nops(sigma_list) do
		    new_avg:=new_avg+exp(sigma_list[rr]-L_value);
		od:
		new_avg:=log(new_avg/nops(sigma_list))+L_value:
		new_avg:=new_avg-log_factorial_product:
		lower_ans:=evalf(log((NN*1.0)!))-NN*log(NN*1.0)+evalf(log((NN*1.0+m*n-1.0)!))
-evalf(log((m*n-1.0)!*1.0))+new_avg:
		upper_ans:=new_avg+(NN/t_val)*evalf(log(t_val!*1.0))-NN*1.0*evalf(log(t_val))
                   +evalf(log((NN+m*n-1.0)!))-evalf(log((m*n-1.0)!*1.0)):
		newupper_ans:=new_avg+log_nicer_upper_bound(row_margins, column_margins)
                   +evalf(log((NN+m*n-1.0)!))-evalf(log((m*n-1.0)!*1.0)):   
		printf("num_sample, log(lower bound), log(upper bound), log(new upper bound): %d, %e, %e, %e\n", ii, lower_ans, upper_ans, newupper_ans); 
	    fi:
	od:

	# Now deal with issue 4 above.
	new_avg:=0;
	for rr from 1 to nops(sigma_list) do
	    new_avg:=new_avg+exp(sigma_list[rr]-L_value);
	od:
	new_avg:=log(new_avg/nops(sigma_list))+L_value:
	non_scaled_new_avg:=evalf(new_avg*1.0);	
	# Now the rest is done with this "more refined tropical" average
	new_avg:=new_avg-log_factorial_product:
	lower_ans:=evalf(log((NN*1.0)!))-NN*log(NN*1.0)+evalf(log((NN*1.0+m*n-1.0)!))
-evalf(log((m*n-1.0)!*1.0))+new_avg:		
	upper_ans:=new_avg+(NN/t_val)*evalf(log(t_val!*1.0))-NN*1.0*evalf(log(t_val))
                   +evalf(log((NN+m*n-1.0)!))-evalf(log((m*n-1.0)!*1.0)):

	newupper_ans:=new_avg+log_nicer_upper_bound(row_margins, column_margins)
                      +evalf(log((NN+m*n-1.0)!))-evalf(log((m*n-1.0)!*1.0)):
	RETURN([non_scaled_new_avg, new_avg, lower_ans, upper_ans, newupper_ans, beta]):	
end:


# Log of the factor for a "nicer upper bound" for the number of contingency tables,
# returns: min{\prod_i r_{i}^{-r_{i}} r_{i}!, \prod_j c_{j}^{-c_{j}} c_{j}!}  
`log_nicer_upper_bound`:=proc(row_margins, column_margins)

	local ii, jj, rowans, colans;
	rowans:=0;
	colans:=0;
	for ii from 1 to nops(row_margins) do
	    rowans:=rowans-row_margins[ii]*evalf(log(row_margins[ii]*1.0))+evalf(log(1.0*row_margins[ii]!)):
	od:
	for jj from 1 to nops(column_margins) do
	    colans:=colans-column_margins[jj]*evalf(log(column_margins[jj]*1.0))+evalf(log(1.0*column_margins[jj]!)):
	od:
	RETURN(min(rowans,colans)):
end:

# Double stochasticize a matrix (similar to the July06 version, with slight differences)
`Db_stochize_adjust_numerics`:=proc(AA, row_margins, col_margins, error1)
	local flag, ii, jj, current_total, numrows, numcols, log_sigma_A, max_diff,
	      maxerror, diag, A, mu_array, lambda_array, count;

	A:=matrix(AA):
	numrows:=linalg[rowdim](A):
	numcols:=linalg[coldim](A):

	mu_array:=[]:     # row scalars
	lambda_array:=[]: # column scalars

	for ii from 1 to numrows do
	    mu_array:=[op(mu_array),1]:
	od:
	
	for jj from 1 to numcols do
	    lambda_array:=[op(lambda_array),1]:
	od:
 
	flag:=0:
	count:=0:
	while(flag=0) do 		     
		      for ii from 1 to numrows do
			  current_total:=0:
			  for jj from 1 to numcols do
			      current_total:=1.0*current_total+A[ii,jj]:
			  od:
			  for jj from 1 to numcols do
			      A[ii,jj]:=(A[ii,jj]/current_total)*op(ii,row_margins):
			  od:
			  mu_array[ii]:=mu_array[ii]*(current_total/op(ii,row_margins)):
		      od:

		     
		      for jj from 1 to numcols do
			  current_total:=0:
			  for ii from 1 to numrows do
			      current_total:=1.0*current_total+A[ii,jj]:
			  od:
			  for ii from 1 to numrows do
			      A[ii,jj]:=(A[ii,jj]/current_total)*op(jj, col_margins):
			  od:
			  lambda_array[jj]:=lambda_array[jj]*(current_total/op(jj, col_margins)):
		      od:

		      
		      # determine if we are within the error bounds
		      maxerror:=-1;
		      for ii from 1 to numrows do
			  current_total:=0:
			  for jj from 1 to numcols do
			      current_total:=1.0*current_total+A[ii,jj]:
			  od:
			  maxerror:=max(maxerror, abs(current_total-op(ii,row_margins))):
		      od:

		      for jj from 1 to numcols do
			  current_total:=0:
			  for ii from 1 to numrows do
			      current_total:=1.0*current_total+A[ii,jj]:
			  od:
			  maxerror:=max(maxerror, abs(current_total-op(jj, col_margins))):
		      od:
		      if(maxerror<error1) then
			  flag:=1:
		      fi:
		      count:=count+1:
	od:

	log_sigma_A:=0:
	for ii from 1 to numrows do
	    log_sigma_A:=log_sigma_A+(row_margins[ii])*log((mu_array[ii]*row_margins[ii])):
	od:
	
	for jj from 1 to numcols do
	    log_sigma_A:=log_sigma_A+(col_margins[jj])*log(lambda_array[jj]*col_margins[jj]):
	od:	

	# print(DBSTOCHEDNUMTIMES, count):
	RETURN([A, log_sigma_A]):
end:

# Generates a random matrix and scales it into the simplex 
`create_A_matrix_simplex_beta_tricky_060220`:=proc(m, n, beta, Rmax_matrix) 

	local pos, aa, bb, ss, tt, total_entries, random_list, pre_A_matrix;

	pre_A_matrix:=array(1..m, 1..n):
	random_list:=[stats[random, exponential[1,0]](m*n)]:
	pos:=1:
	for aa from 1 to m do
	    for bb from 1 to n do
		pre_A_matrix[aa,bb]:=op(pos,random_list):
		pos:=pos+1:
	    od:
	od:

	# sum up entries in pre_A_matrix and divide 
	total_entries:=0:
	for ss from 1 to m do
		for tt from 1 to n do
		    total_entries:=total_entries+Rmax_matrix[ss,tt]*pre_A_matrix[ss,tt]^(beta):
		od:
	od: 

	for ss from 1 to m do
		for tt from 1 to n do
		    pre_A_matrix[ss,tt]:=(Rmax_matrix[ss,tt]*pre_A_matrix[ss,tt]^(beta))/total_entries:
		od:
	od:
	RETURN(pre_A_matrix);
end:


# Miscellaneous utilities below
# -----------------------------

# a sorting algorithm we need
`genzoom_list_sorter`:=proc(sigma_listin, per_listin, poss_sigma_entry, poss_per_list_entry, k)

	local ii, sigma_list, per_list;
	
	sigma_list:=[];
	per_list:=[];
	for ii from 1 to nops(sigma_listin) do
	    sigma_list:=[op(sigma_list), sigma_listin[ii]]:
	    per_list:=[op(per_list), per_listin[ii]]:
	od:

	# first deal with case the lists are empty
	if(nops(sigma_list)=0) then
		sigma_list:=[op(sigma_list),poss_sigma_entry]:
		per_list:=[op(per_list),poss_per_list_entry]:
		RETURN([sigma_list,per_list]):
	fi:

	# now the case that the poss_sigma_entry is too small
	if(poss_sigma_entry<min(op(sigma_list))) then
		if(nops(sigma_list)=k) then
		    RETURN([sigma_list, per_list]):
		else
		    sigma_list:=[poss_sigma_entry, op(sigma_list)]:
		    per_list:=[poss_per_list_entry, op(per_list)]:
		    RETURN([sigma_list,per_list]):
		fi:
	fi:

	# now the case that the lists are of 0< size <k, and we will insert 
	# poss_sigma_entry somewhere 
	if(nops(sigma_list)<k) then
		if(max(op(sigma_list))<poss_sigma_entry) then
			sigma_list:=[op(sigma_list), poss_sigma_entry]:
			per_list:=[op(per_list), poss_per_list_entry]:
			RETURN([sigma_list,per_list]):  			
		fi:
		for ii from 1 to nops(sigma_list) do
			if(op(ii,sigma_list)>=poss_sigma_entry) then
					sigma_list:=[op(sigma_list[1..(ii-1)]), poss_sigma_entry, op(sigma_list[ii..nops(sigma_list)])]:
					per_list:=[op(per_list[1..(ii-1)]), poss_per_list_entry, op(per_list[ii..nops(per_list)])]:
					RETURN([sigma_list,per_list]):
			fi:      
		od:
	else  # now we have exactly k entries
		if(max(op(sigma_list))<poss_sigma_entry) then
			sigma_list:=[op(sigma_list[2..nops(sigma_list)]), poss_sigma_entry]:
			per_list:=[op(per_list[2..nops(sigma_list)]), poss_per_list_entry]:
			RETURN([sigma_list,per_list]):  			
		fi:
		for ii from 1 to nops(sigma_list) do
			if(op(ii,sigma_list)>=poss_sigma_entry) then
					sigma_list:=[op(sigma_list[2..(ii-1)]), poss_sigma_entry, op(sigma_list[ii..nops(sigma_list)])]:
					per_list:=[op(per_list[2..(ii-1)]), poss_per_list_entry, op(per_list[ii..nops(per_list)])]:
					RETURN([sigma_list,per_list]):
			fi:      
		od:	        
	fi:
end:

# Code to compute the Diaconis-Efron heuristic approximation formula
# "which seems to work with m,n small and N large".
`Diaconis_Efron`:=proc(row_vector,column_vector)
	local mm, nn, NN, ww, rr_bar_vector, cc_bar_vector, rr_bar_sq,
             cc_bar_sq, kk, ii, ans, rr_bar_prod, cc_bar_prod;

	mm:=nops(row_vector):
	nn:=nops(column_vector):
	NN:=sumarrayelts(row_vector, mm):
	ww:=1/(1+(mm*nn)/(2*NN)):
	rr_bar_vector:=[]:
	for ii from 1 to mm do
	    rr_bar_vector:=[op(rr_bar_vector), ((1-ww)/mm)+(ww*op(ii, row_vector)/NN)]:
	od:
	cc_bar_vector:=[]:
	for ii from 1 to nn do
	    cc_bar_vector:=[op(cc_bar_vector), ((1-ww)/nn)+(ww*op(ii, column_vector)/NN)]:
	od:	
	rr_bar_sq:=0:
	for ii from 1 to mm do
	    rr_bar_sq:=rr_bar_sq+(op(ii,rr_bar_vector)^2):
	od:
	kk:=((nn+1)/(nn*rr_bar_sq))-(1/nn):

	ans:=(((2*NN+mm*nn)/2)^((mm-1)*(nn-1)))*GAMMA(nn*kk)/((GAMMA(nn)^mm)*(GAMMA(kk)^nn)):
	rr_bar_prod:=1:
	for ii from 1 to mm do
	    rr_bar_prod:=rr_bar_prod*op(ii, rr_bar_vector):
	od:
	cc_bar_prod:=1:
	for ii from 1 to nn do
	    cc_bar_prod:=cc_bar_prod*op(ii, cc_bar_vector):
	od:	
	ans:=evalf(ans*(rr_bar_prod^(nn-1))*(cc_bar_prod^(kk-1))):
	RETURN(ans):	
end:

# Compute t=max{min_{i=1..m}r_i, min_{j=1..n}c_j} 
`t_value`:=proc(m,n, row_margins, column_margins)
	local ii, r_min, c_min, ans;

	r_min:=infinity;
	for ii from 1 to m do
	    if(op(ii, row_margins)<r_min) then
		     r_min:=op(ii, row_margins):
	    fi:
	od:
	c_min:=infinity;
	for ii from 1 to n do
	    if(op(ii, column_margins)<c_min) then
		     c_min:=op(ii, column_margins):
	    fi:
	od:	
	ans:=max(r_min, c_min);
	RETURN(ans):
end:

# Create an array of all length n, all having entry "same_array"
`same_array`:=proc(same_entry, n)
	local the_array, ii;

	the_array:=[]:
	for ii from 1 to n do
	    the_array:=[op(the_array),same_entry]:	    
	od:
	RETURN(the_array):
end:

# Sum the first k elements of an array. 
`sumarrayelts`:=proc(array_in, kk)
			       
	local ii, sumtotal;
	
	sumtotal:=0:
	for ii from 1 to kk do
	    sumtotal:=sumtotal+op(ii,array_in):
	od:
	RETURN(sumtotal):
end:

# A user interface.
`contingency`:=proc()
	local ans, m,n, num_steps, num_samples, num_iterations_before_sample,
              rowend_marginsin, colend_marginsin, final_steps, error1, YN,
              weightmatrix, row_margins, column_margins:

printf("*************************************************\n");
printf("Contingency table enumeration estimator\n");
printf("*************************************************\n");
printf("Remember to end each input with a semicolon.\n");
m:=readstat("Number of rows m: ");
n:=readstat("Number of columns n: ");
row_margins:=readstat("Row margins, e.g., [3,3,...,3]: ");
column_margins:=readstat("Column margins, e.g., [3,3,...,3]: ");
num_steps:=readstat("Number of telescoping steps: ");
num_samples:=readstat("Vector of hit and run steps: ");
num_iterations_before_sample:=readstat("Vector of hit and run steps to wait between samples: ");
rowend_marginsin:=readstat("Row margins of the final step: ");
colend_marginsin:=readstat("Column margins of the final step: ");
final_steps:=readstat("Number of samples for the sample mean for the final step: ");
error1:=readstat("The error used in the sigma calculation. Example 0.001: ");
ans:=markov_telescope_July06(m,n, row_margins, column_margins, 
                                   num_steps, num_samples, 
                                   num_iterations_before_sample, final_steps, 
                                   rowend_marginsin, colend_marginsin, 
                                   error1, 0);
print(ans);
end: