## 26 May 2013
# V9

with(LinearAlgebra):

interface(rtablesize=infinity);

intM:=proc(K::integer, C::integer)

description "Creat interaction terms (matrix form).";

local 
	M:=Matrix(K,C,symbol=mix),
	tolt,cols,rows;

# Set constraints

	M(..,1):=0;
	M(1..2,..):=0;
	
	tolt:=convert(M,`+`);
	cols:=MTM:-sum(M,1);
	rows:=MTM:-sum(M,2);
	
	M(2,..):=-cols;
	M(..,1):=-rows;
	M(2,1):=tolt;
	
	M;
	
end proc;


indprob:=proc(x,
	K::integer,
	C::integer,
	tvarphi::truefalse:=false,
	tvarp::truefalse:=false,
	phiMix::truefalse:=false
	)

description "Compute the Probability for given capture history $x$.";

local
	# phi,p, # Vectors of phi and p
	int0,  # Matrix of interaction terms
	ind,   # Find indices of where the '1's are
	fcap,  # First capture
	li,    # Last capture
	ans;   # Desired probability

global phi,p; # Vectors of phi and p;

# take care of phi (nu)
# SINGLE CLASS
	
	if tvarphi then
		phi:=Array(1..K,[seq(1/(1+exp(-nu[j])),j=1..K)]);
		# phi:=Array(1..K,[seq(nu[j],j=1..K)]);
	else
		phi:=Array(1..K,[seq(1/(1+exp(-nu)),j=1..K)]);
		# phi:=Array(1..K,[seq(nu,j=1..K)]);
	end if;
	phi(K):=0;
	
# take care of p

	if tvarp then
		if phiMix then
			int0:=intM(K,C);
			p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j]+eta[c]+int0[j,c]))));
			# p:=Array(1..K,1..C,(j,c)->mu+tau[j]+eta[c]+int0[j,c]);
		else
			if C=1 then
				p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j]))));
				# p:=Array(1..K,1..C,(j,c)->mu+tau[j]);
			else
				p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+tau[j]+eta[c]))));
				# p:=Array(1..K,1..C,(j,c)->mu+tau[j]+eta[c]);
			end if;
		end if;
	else
			p:=Array(1..K,1..C,(j,c)->1/(1+exp(-(mu+eta[c]))));
			# p:=Array(1..K,1..C,(j,c)->mu+eta[c]);
	end if;
	
	p(1,..):=0;
	p:=subs({tau[2]=0,eta[1]=0},p);

	ind:=ArrayTools[SearchArray](Array(x));
	fcap,li:=ind[1],ind[-1];
	
	ans:=
		add(
			w[c]*
			add(
				(1-phi[d])*
				mul(phi[j],j=fcap..(d-1))
				*
				mul(p[j,c]^(x[j])*(1-p[j,c])^(1-x[j]),j=(fcap+1)..d), 
				d=li..nops(x)
			),
		c=1..C);
	
	ans;
	
end proc;


myP:=proc(
	K::integer,
	C::integer,
	tvarphi::truefalse:=false,
	tvarp::truefalse:=false,
	phiMix::truefalse:=false,
	dosum::truefalse:=false,
	n::integer:=30)

description "Test for Parameter Redundancy (for mixtures of $p$).";	

local
	CH,      # All Capture Histories from $K$ sample
	kappa,   # Exhaustive Summary
	allpars, # Complete Set of Model Parameters
	DD1,     # Derivative Matrix (Jacobian) of Kappa
	
	## Numerical ##
	sys,     # System of eqns
	myrand,  # Random number generator

	phitrue, # Independent Set of phi
	phinum,  # Random points
	phiparM, # Vector of phi
	
	ptrue,   # Independent Set of p
	pnum,    # Random points
	pparM,   # Vector of p
	
	ws,      # Random points
	wsnum,   # Vector of w
	
	num,     # Complete Set of Random Points
	DD1num,  # Numerical Derivative Matrix

	rS,      # Symbolic Rank
	rN       # Numerical Rank
	;
	
# Compute Capture Histories

	CH:=combinat[permute]([seq(s$K,s=[0,1])],K);
	CH:=CH[2..-1];

# Get Exhaustive Summary

	kappa:=map(indprob,CH,K,C,tvarphi,tvarp,phiMix);
	kappa:=subs({w[C]=1-add(w[i],i=1..(C-1))},kappa);
	# kappa:=subs({seq(w[i]=1/(1+exp(-wnew[i])),i=1..(C-1))},kappa);
	kappa:=kappa[2..-1];
	
# Get Jacobian
	allpars:=indets(kappa,name);
	# pars:=allpars minus {seq(w[i],i=1..(C-1))};
	DD1:=VectorCalculus[Jacobian](Vector(kappa),convert(allpars,list));

	# rS:=Rank(DD1);

# Initialise parameters

	myrand:=evalf[2](rand(10..90)/100);
	
	# Take care of phi
	phitrue:=convert(phi,set) minus {0};
	sys:={op(phitrue)=~seq(myrand(),i=1..nops(phitrue))};
	phinum:=solve(sys,indets(phitrue,name));
	phiparM:=subs(phinum,phi);
	phiparM:=evalf(phiparM);
	phiparM:=phiparM[1..-2];

	# Take care of p
	ptrue:=convert(p,set):
	if tvarp then
		if phiMix then
			ptrue:=ptrue minus {0};
		else
			ptrue:=ptrue[2..(K+C-1)];
		end if;
	else
		ptrue:=ptrue minus {0};
	end if;
	
	sys:=ptrue=~seq(myrand(),i=1..nops(indets(ptrue,name)));
	pnum:=solve({sys},indets(ptrue,name));
	pparM:=subs(pnum,p);
	pparM:=evalf(pparM);
	pparM:=pparM(2..,..);
	
	# Take care of weights
	wsnum:=[seq(w[i]=myrand(),i=1..(C-1))];
	ws:=convert(rhs~(wsnum),Vector);
	wsnum:=convert(wsnum,set);
	
	# Get num into DD1
	num:=phinum union pnum union wsnum;
	DD1num:=eval(DD1,num);
	rN:=Rank(DD1num);

# Print Summary #
if dosum then

	printf("Vector for phi \n");
	printf("%.4f\n\n",phiparM);
	
	printf("Parameter Matrix for p \n");
	printf("%.4f\n\n",pparM);
	
	
	printf("With vector of weights \n");
	printf("%.4f\n\n",ws);
	
	#printf("Number of Parameters = %d, Numerical Rank = %d, defficiency = %d.\n", nops(allpars), rN, nops(allpars)-rN);
	printf("Number of Parameters = %d.\n",nops(allpars));
	printf("Numerical Rank = %d.\n",rN);
	printf("Defficiency = %d.\n",nops(allpars)-rN);
	printf("---------------------------------------------------------------\n");
end if;
# Print Summary # -End- #

	return rN;
	
end proc;