function [res]=listplans(N,m,gamma,delta,beta,plan,plantype)
%generate results for all possible SP and AP plans (no leveraging)
%output gives the best type A and B augmented plans, as well as the
%standard plans with n=1 and n=2
%assume m fixed operators, if beta=1 this simplifies to the no interactio model case
%N is the total number of measurements
%gamma, delta and beta are as defined in the paper all must be between 0 and 1
%gamma is sqrt((so2+sd2+sm2)/(sp2+so2+sd2+sm2))
%plan and plantype are optional parameters if you want to compare a specific plan
%plan should be a vector giving [k,n]
%plantype should be S A or B, for augmented plans ka or kb is then choosen to give a total of N measurements
%res ranks the SP and AP A or B plans by asymptotic standard deviation for gamma
%must choose N as a multiple of m the number of operators!

%example calls 
%[res]=listplans(60,2,.8,.5,.5);
%[res]=listplans(60,3,.9,.5,.5,[10,2],'S');
%run without ending semicolon to see results for all possible plans
%note that the plans types have a numerical code 'S'=0 'A'=100 and 'B'=200

%translate gamma, delta and beta to sp2, so2, sd2 and sm2
sm2=delta*gamma^2;
so2=beta*(1-delta)*gamma^2;
sd2=(1-beta)*(1-delta)*gamma^2;   %interaction
sp2=1-so2-sd2-sm2;   %assume sum is 1

if N/m~=floor(N/m),  disp('choose N as a multiple of number of operators'), end;

%check that plan and plantype make sense, if it is given
if nargin>=6,
    switch plantype
        case 'S'
          plantype=0;
          if m*plan(1)*plan(2)~=N   %standard plan
            disp(['ERROR - the total number of measurements for your standard plan must be ',num2str(N)])
            disp(['your plan has ',num2str(m*Plan(1)*Plan(2)),' measurements'])
          end;
        case 'A'
          plantype=100;
          ka=(N-m*plan(1)*plan(2));
          if ka/m~=round(ka/m)
            disp(['ERROR - the number of measurements in the augmented part must be divisable by the number of operators'])
          end;
        case 'B'
          plantype=200;
          kb=(N-m*plan(1)*plan(2))/m;
          if kb~=round(kb)
            disp(['ERROR - the number of measurements in the augmented part must be divisable by the number of operators'])
          end;
        otherwise
          disp('ERROR - plantype must be S A or B')
    end;
end;
    
%find feasible values of the mus that match so2
%assume close to balance for the relative biases among the operators, this
%this makes no difference as the Fisher information does not depend on so2 or muj
%it does change the change of variables matrix, but D*D' is the same 
muj=zeros(m,1);
if m>1,  %i.e. there is more than a single operator
    if (m/2)==round(m/2),  %m even
        mu=sqrt(so2);
    else
        mu=sqrt(m*so2/(m-1));
    end;
    muj(1:floor(m/2))=mu;
    muj(ceil(m/2)+1:end)=-mu;
end;

if m==1 && so2>0, disp('WARNING, when there is only one operator the so2 value is ignored, it is set to zero'), end;

%find change of variables matrix
D=cmatrix(muj,sp2,sd2,sm2);

res=[];  %save all the results [k,n, ..., se(gamma)]
%k = # parts,   n = # measurements on each part
for k=2:N/m,  %loop through all possible number of parts, avoid having only 1 part in SP
   nmax=floor(N/(m*k));  %largest possible value for n
   if (N/(m*k))==nmax,   %SP possible
       if nmax>1 || (nmax==1 && sd2==0  && m>=2),
          Jmat=J_kmn(k,m,nmax,sp2,sd2,sm2);
          varmat=inv(Jmat);  %diagonals give variance for sp2, sd2 and sm2
          SEsp=sqrt(varmat(m+1,m+1))/(2*sqrt(sp2));      %determine the SE for sp  (not sp2)
          temp=D*inv(Jmat)*D';
          if sd2~=0,  
              res=[res;k,nmax,0,0,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(4,4)),sqrt(temp(1,1))];
          else
              if m>=2,  
                  res=[res;k,nmax,0,0,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(1,1))]; 
              else
                  res=[res;k,nmax,0,0,SEsp,sqrt(temp(2,2)),sqrt(temp(1,1))];  
              end;
          end;
          nmax=nmax-1;  %don't use an AP with nmax since that is an SP
       end;
   end;
   for n=2:nmax,   %loop through looking for all possible type A and B APs
      leftmeas=N-m*n*k;  %total # measurement left for Stage 2
      if floor(leftmeas/m)==leftmeas/m,  %lmeas is divisible by m and we can find APs
         u=leftmeas/m;   %number of measurements for each operator in Stage 2
         %type A AP
         Jmat=J_kmn(k,m,n,sp2,sd2,sm2)+J_typeA(m,u,sp2,sd2,sm2);  %add the two information matrices
         varmat=inv(Jmat);  %diagonals give variance for sp2, sd2 and sm2
         SEsp=sqrt(varmat(m+1,m+1))/(2*sqrt(sp2));      %determine the SE for sp  (not sp2)
         temp=D*inv(Jmat)*D';
         %ka = u*m for a AP type A
         if sd2~=0,   
             res=[res;k,n,u*m,100,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(4,4)),sqrt(temp(1,1))];
         else
             if m>=2,  
                 res=[res;k,n,u*m,100,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(1,1))];
             else
                 res=[res;k,n,u*m,100,SEsp,sqrt(temp(2,2)),sqrt(temp(1,1))]; 
             end;
         end;
         %type B AP
         Jmat=J_kmn(k,m,n,sp2,sd2,sm2)+J_kmn(u,m,1,sp2,sd2,sm2);  %add the two information matrices
         varmat=inv(Jmat);  %diagonals give variance for sp2, sd2 and sm2
         SEsp=sqrt(varmat(m+1,m+1))/(2*sqrt(sp2));      %determine the SE for sp  (not sp2)
         temp=D*inv(Jmat)*D';
         %kb=u for a AP type B
         if sd2~=0,   
             res=[res;k,n,u,200,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(4,4)),sqrt(temp(1,1))];
         else
             if m>=2,  
                 res=[res;k,n,u,200,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(1,1))];  
             else
                 res=[res;k,n,u,200,SEsp,sqrt(temp(2,2)),sqrt(temp(1,1))];  
             end;
         end;
      end;
   end;  
end;

%also add the case of an SP with n=1 when sd>=0, i.e. interaction model), this plan allows estimation of gamma,
%but does not allow estimation of sigmm2 and sigmad2 seperately
%to get results for this case we use the no interaction code, i.e. set sd2=0, but increase sm2 to inlcude the interaction variance
if sd2>0, 
    D=cmatrix(muj,sp2,0,sd2+sm2);  %change of variables matrix
    Jmat=J_kmn(N/m,m,1,sp2,0,sd2+sm2);
    varmat=inv(Jmat);  %diagonals give variance for sp2, sd2 and sm2
    SEsp=sqrt(varmat(m+1,m+1))/(2*sqrt(sp2));      %determine the SE for sp  (not sp2)
    temp=D*inv(Jmat)*D';
    %temp=inv(Jfinal);   %diagonals give SE for sp2, sd2 and sm2
    if sd2~=0,   
        res=[res;N/m,1,0,0,SEsp,sqrt(temp(2,2)),0,sqrt(temp(3,3)),sqrt(temp(1,1))];
    else
        res=[res;N/m,1,0,0,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(1,1))];
    end;
end;

%rank the results by standard error of gamma
[temp,ind]=sort(res(:,end));
res=res(ind,:);

gam=sqrt((so2+sd2+sm2)/(sp2+so2+sd2+sm2));
%display results
disp(['number of operators = ',num2str(m),'  total number of measurements = ',num2str(N)])
disp(['  gamma = ',num2str(gam),'  delta = ',num2str(delta),'  beta = ',num2str(beta)])
if beta==1,
    disp('use the no interaction model (no sd2) since beta = 1')
    disp(['sigmap2 = ',num2str(sp2),'  sigmao2 = ',num2str(so2),'  sigmam2 = ',num2str(sm2)])
else
    disp('use the interaction model')
    disp(['sigmap2 = ',num2str(sp2),'  sigmao2 = ',num2str(so2),'  sigmad2 = ',num2str(sd2),'  sigmam2 = ',num2str(sm2)])
end;
if sd2~=0,   disp('      k       n     ka/b    plan      SEsp      SEso      SEsd      SEsm    SEgamma   SEratio')
else
    if m>=2,  disp('      k       n     ka/b    plan      SEsp      SEso      SEsm    SEgamma   SEratio') 
    else  disp('      k       n     ka/b     plan      SEsp      SEsm    SEgamma   SEratio'),   end;
end
nplans=min(20,length(res(:,1)));
topSP=find(res(:,4)==0,1);    %find top SP, will have n=1
topSPs=find(res(:,4)==0,2);    %find all SP
topSP2=topSPs(2);  %top SP other than the one with n=1
topAPA=find(res(:,4)==100,1);    %find top AP A
topAPB=find(res(:,4)==200,1);    %find top AP B
%display the top plan in order based on SE(gamma)
topplans=res([topAPB,topAPA,topSP2,topSP],:);   %all the top plans

if nargin>=6,   %also find specific plan
    temp=find(res(:,4)==plantype);  
    temp2=find(res(:,1)==plan(1));  
    temp3=find(res(:,2)==plan(2));  
    temp4=intersect(temp,intersect(temp2,temp3));   %should be one value
    topplans=[topplans;res(temp4,:)];   %add specific plan to list
end;
topplans=[topplans,res(topSP,end)./topplans(:,end)];   %also calculate efficiency ratio (compare to to Standard plan)

[ntop,temp]=size(topplans);
for ii=1:ntop,   %look at each plan seperately
    switch topplans(ii,4)
        case 0, pp='S';
        case 100, pp='A';
        case 200, pp='B';
    end;
    disp(['      ',num2str(topplans(ii,1:3),'%8d'),'      ',pp,'      ',num2str(topplans(ii,5:end),'%10.4f')])
end;


function [Jmat]=J_kmn(k,m,n,sp2,sd2,sm2)
%find the information matrix for the SP given k, m, n
if sd2 ~=0,  Jmat=zeros(m+3,m+3);   else   Jmat=zeros(m+2,m+2);   end;

%Note that the information matrix Jmat is always block diagonal with terms
%related to the mujs and another block for the rest

b1 = 1/sm2;
b2 = -sd2 / (sm2 * (sm2 + n * sd2));
b3 = -sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2);

%setup top bit related to the mus
Jmat(1:m,1:m)=-k*n*n*b3;
for i=1:m,   %replace the diagonal
    Jmat(i,i)=-k*n*(b1+n*(b2+b3));
end;

%expected vlaues of the various sums of squares
A=k*m*n*(sp2+sd2+sm2);
B=k*m*n*(n*sp2+n*sd2+sm2);
C=k*m*n*(m*n*sp2+n*sd2+sm2);

%now look at bottom right 3x3 part related to sp2, sd2 and sm2 in that order!
if sd2~=0, 
    Jmat(m+1,m+1) = -0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m ^ 2 * n ^ 2 + k * m ^ 2 * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+1,m+2) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C * n / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n ^ 2 / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n ^ 2 + k * m * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+2,m+1) = Jmat(m+1,m+2);
    Jmat(m+1,m+3) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n + k * m * n / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+3,m+1) = Jmat(m+1,m+3);
    Jmat(m+2,m+2) = -(1 / sm2 / (sm2 + n * sd2) ^ 2 * B * n) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B * n ^ 2) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C * n ^ 2) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n ^ 2) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * n ^ 2) - (k * (2 * n ^ 2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 * n ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * (n * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * n) / 0.2e1 + (k * (n * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) * n / (sm2 + n * sd2)) / 0.2e1;
    Jmat(m+2,m+3) = -0.1e1 / sm2 ^ 2 / (sm2 + n * sd2) * B / 0.2e1 - 0.1e1 / sm2 / (sm2 + n * sd2) ^ 2 * B / 0.2e1 + sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B * n / 0.2e1 + sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B * n + sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C * n + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * n - k * (0.2e1 * n * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) ^ 2 / (sm2 + n * sd2) ^ 2 * n * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) * (m - 0.1e1) / (sm2 + n * sd2) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) * (m * n - m) / sm2 / 0.2e1;
    Jmat(m+3,m+2) = Jmat(m+2,m+3);
    Jmat(m+3,m+3) = -(1 / sm2 ^ 3 * A) + (sd2 / sm2 ^ 3 / (sm2 + n * sd2) * B) + (sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C) - (k * (2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) ^ 2 / sm2 ^ 2 - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 ^ 2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) / (sm2 + n * sd2)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m * n - m) / sm2) / 0.2e1;
else  %remove bits related to sd2
    Jmat(m+1,m+1) = -0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m ^ 2 * n ^ 2 + k * m ^ 2 * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+1,m+2) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n + k * m * n / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+2,m+1) = Jmat(m+1,m+2);
    Jmat(m+2,m+2) = -(1 / sm2 ^ 3 * A) + (sd2 / sm2 ^ 3 / (sm2 + n * sd2) * B) + (sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C) - (k * (2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) ^ 2 / sm2 ^ 2 - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 ^ 2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) / (sm2 + n * sd2)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m * n - m) / sm2) / 0.2e1;
end;
%change the sign since we want the negative of the expected value instead of the expected value
Jmat=-Jmat;


function [Jmat]=J_typeA(m,u,sp2,sd2,sm2)
%find the special extra information matrix for the type A AP given m and u
%with AP plan A m operators each measure u *different* parts once each
%total number of parts in extra AP bit is m*u

if sd2 ~=0,  Jmat=zeros(m+3,m+3);
else   Jmat=zeros(m+2,m+2);   end;

%setup top bit related to the mus
Jmat(1:m,1:m)=-u/(sp2+sd2+sm2)*eye(m);  

%n=1; k=u;   %special case
%now look at bottom right 3x3 part related to sp2, sd2, sm2
A=m*u*(sp2+sd2+sm2);
ka=m*u;
if sd2~=0, 
    Jmat(m+1:m+3,m+1:m+3) = -A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;   %3x3 matrix all the same
    %off diagonal elements also all the same
    Jmat(m+1,m+2)=-A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;
    Jmat(m+1,m+3)=Jmat(m+1,m+2);
    Jmat(m+2,m+3)=Jmat(m+1,m+2);
    Jmat(m+2,m+1)=Jmat(m+1,m+2);
    Jmat(m+3,m+2)=Jmat(m+1,m+2);
    Jmat(m+3,m+1)=Jmat(m+1,m+2);
else  %remove bits related to sd2
    Jmat(m+1:m+2,m+1:m+2) = -A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;   %3x3 matrix all the same
    %off diagonal elements also all the same
    Jmat(m+1,m+2)=-A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;
    Jmat(m+2,m+1)=Jmat(m+1,m+2);
end;
%change the sign since we want the negative of the expected value instead of the expected value
Jmat=-Jmat;


function [mat]=cmatrix(mu,sp2,sd2,sm2)
%returns the change of variables matrix based on the partial derivatives
%go from muis sp2, sd2 and sm2 (note variance scale) to gamma, so, sd, sm in that order
%skip row and column for sd if we assume the no interaction model
%also remove the so column if there is only one operat, i.e. m=1

m=length(mu);  %number of operators
if sd2~=0,   %interaction model, this doesn't make sense if there is only one operator
    mat=zeros(4,m+3);   %four rows
    if m==1, disp('WARNING, using the interaction model when there is only one operator does not make sense, select sd2=0'), end;
    mat(3,m+2)=1/(2*sqrt(sd2));  %only element in third row for sd
    mat(4,m+3)=1/(2*sqrt(sm2));  %only element in fourth row for sm
else    %no interaction model
    if m>=2
       mat=zeros(3,m+2);   %three rows
       mat(3,m+2)=1/(2*sqrt(sm2));  %only element in third row for sm
    else  %only one operator case
       mat=zeros(2,m+2);   %two rows, only for gamma and sm
       mat(2,m+2)=1/(2*sqrt(sm2));  %only element in second row for sm
    end;
end;

muiminusmubar=mu-mean(mu);
so=sqrt(sum(muiminusmubar.^2)/m);
if m>=2,  mat(2,1:m)=muiminusmubar/(m*so);   end; %for so

%can't include sp in change of variable matrix, otherwise there is a redundancy
%mat(2,m+1)=1/(2*sqrt(sp2));   %only element in second row for sp

%first row for gamma
gam=sqrt((so^2+sd2+sm2)/(sp2+so^2+sd2+sm2));
mat(1,1:m)=muiminusmubar*(1-gam^2)/gam/(sp2+so^2+sd2+sm2)/m;
so2=so^2;  sp=sqrt(sp2); sd=sqrt(sd2);  sm=sqrt(sm2);

%change of variables on the variance scale
mat(1,m+1) = -((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2 / 0.2e1;
if sd2~=0, 
   mat(1,m+2) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
   mat(1,m+3) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
else
   mat(1,m+2) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
end;