/*  gencaph1.c - Generates capture history records.          */

/* Here's how it works:
 
The expected numbers of animals exhibiting each capture history are computed
using a recursive algorithm coded in C. The algorithm follows a population of
animals, exposing them to capture and survival until death or the end of the
study, while saving a vector of codes (0=not captured, 1=captured) indicating
capture history. The process is repeated for each cohort of injected animals
in capture occasions after the first.
 
There are two modes of operation of the program - Deterministic and
Stochastic. Under Deterministic mode, the number of animals is retained
throughout the process as a fractional number and is multiplied by the
parameters to get the final number of animals with each capture-history. This
results in two things: capture-history frequencies are not integers, and
parameter estimates will match input parameters exactly (almost). In
Stochastic mode, the number of animals is achieved via a simulation function
(binomial) on the parameters. In this mode, the number of animals with each
capture-history will be integers, and the data (and resulting estimates) will
be different for each program run.
 
Terminology and Notation for Input Data
N(i) = Number of new animals introduced into the population just before time i 
phi(i) = survival rate from time i to i+1, 
theta(i) = tag-retention rate from age or time i to i+1 
p(i) = capture probability for unmarked animals in time i 
c(i) = capture probability for marked animals in time i 
f(i) = number of new individuals in time i per old individual in time i-1
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#ifdef DOS
 #include <process.h>
 #include <conio.h>
#endif
FILE *fopen(), *fp;
#define eps .00005
#define MXPER 27
long nper, determ=1, dim, mxgrps=1, grp=0;  int thetaflg=0;
double **nx,*pa,*p2,*sa,*tcap,*talive,*th;

//            uniform random number generator
double unif(void) { return(rand()/(double)RAND_MAX); }

//            normal random number generator
double unorm(double xn, double se) { 
    int i; double sumr=-6.; for (i=0; i<12; i++) sumr+=rand()/(double)RAND_MAX; 
    sumr=sumr*se+xn; if (!determ) sumr=floor(sumr+0.5);
    return(sumr); 
}

//            binomial random number generator
double binoml(double n,double p) /* binomial function */
{  double np,tr,b,x,bin,pp; long i,nn; double y;
    bin=0; if (determ)  return(n*p);
    if (p>0)                                      /* compute binomial directly     */
        if (n<=100) { nn=floor(n); for (i=1;i<=nn;bin+=(unif()<=p),i++); }   
        else
        { np=n*p;                                     /* normal approx */
            if (p>=.05 && p<=.95) { for (x=-6,i=1;i<=12;x+=unif(),i++); y=np*(1-p);bin=floor(sqrt(y)*x+np); }
            else {
                pp=p; if (p>.9) pp=1-p;                           /* poisson approx  */
                for ( b=-n*pp, tr=0; tr>=b; bin++) tr+=log(unif());
                if (p>.9) bin=n-bin;
              //  printf("binoml(%f, %f)=%f\n",n,p,bin);
            }
        }
    return(bin);
}

/*         caphis function - recursive routine
                  which follows a cohort of animals
                  through time, keeping track of capture 
                  status at each occasion */
void caphis(long ch,double n,long per, long fper)
{ double nalive, ncap, notcap, ndead, nlosttag, nkepttag, captprob; 
    long i, firstper;
    captprob=pa[per-1]; if (fper>=0) captprob=p2[per-1];
    ncap=binoml(n,captprob); 
    
    tcap[per-1]+=ncap;  notcap=n-ncap;
    if (per<nper) {
        nalive=binoml(ncap,sa[per-1]); talive[per]+=nalive; 
        firstper=fper; if (fper<0) firstper=per; 
        if (thetaflg==0)
            nkepttag=binoml(nalive,th[per-firstper]); 
        else
            nkepttag=binoml(nalive,th[per-1]); 
        nlosttag=nalive-nkepttag;
        ndead=ncap-nkepttag; 
		i=(ch|1<<(nper-per)); 
		nx[grp][i]+=(ndead); 
        if (nkepttag>0) caphis(i,nkepttag,per+1,firstper);
        if (nlosttag>0) caphis(0,nlosttag,per+1,-1);
        nalive=binoml(notcap,sa[per-1]); talive[per]+=nalive; 
        nkepttag=nalive; 
        if (fper>=0) {
            if (thetaflg==0) nkepttag=binoml(nalive,th[per-fper]);
            else nkepttag=binoml(nalive,th[per-1]);
        }
        nlosttag=nalive-nkepttag;
        ndead=notcap-nkepttag; nx[grp][ch]+=ndead;
        if (nkepttag>0) caphis(ch,nkepttag,per+1,fper);
        if (nlosttag>0) caphis(0,nlosttag,per+1,fper);
    }
    else { nx[grp][(ch|1)]+=ncap; nx[grp][ch]+=notcap; }
}

/*  main routine - collects input, calls caphis routine, prints output.
 
 *                 Optionally generates input statements for
 *                 programs JOLLY, MARK or CAPTURE.
 * 
 *                 Options are specified on the command line.
 *                   NOPROMPT - If you double-click in Windows, the program
 *                              will prompt you to hit a key to continue
 *                              (allowing you to view output messages).
 *                              This option prevents this prompt so the
 *                              program can be run in batch mode.
 *                   JOLLY - Generate an input file suitable for JOLLY
 *                   CAPTURE - Generate an input file suitable for CAPTURE
 *                   MARK - Generate an input file suitable for MARK
 *                   EXP - Generate expected values (not simulated)
 *                   DETERM - Generate expected values (not simulated)
 *                   STOCHASTIC - Generate simulated values
 *                   THETA=TIME - Make THETA (tag-retention rates) time specific
 *                                instead of the default (age-specific)
*/   
int main(int argc, char **argv) {
    FILE *F; double *na,*fa,*sef,sm; char infile[512]="gencaph1.txt",outfile[512]="gencaph1.out";
    long i,j,k,tst,jlyflg=0,capflg=0,mrkflg=0,tagflg=0,promptflg=1,recruit=0,klmt,seed=12345,USEF1=0; double sumn;
    printf("usage: gencaph1 options\n");
    printf("  input file:gencaph1.txt\n  format:\n");
    printf("   n1,b2,b3,...     (init pop size, recruitment)\n");
    printf("   s1,s2,s3,...     (survival rates)\n");
    printf("   th1,th2,th3,...  (tag-retension rates)\n");
    printf("   p1,p2,p3,...\n   (capt probs (unmarked)");
    printf("   c1,c2,c3,...  (capt probs (marked)\n");
    printf("   f1,f2,f3,... (if RECRUIT= specified)\n");
    printf("   sef1,sef2,sef3,... (if RECRUIT= specified)\n");
    printf("   n1,n2,n3,... (group 2) or -1,0,0,... (if only 1 group)\n");
	printf(" options:\n");
	printf("  {CAPTURE|JOLLY|MARK} : create output file in program format\n");
	printf("  THETA={AGE|TIME}     : Make THETA (tag-retention rates) time specific\n");
    printf("                                instead of the default (age-specific)\n");
    printf("  EXP                  : Generate expected values instead of simulated\n");
	printf("  STOCHASTIC           : Generate simulated values\n");
	printf("  RECRUIT={BINOMIAL|NORMAL}: specify function for recruitment\n");
	printf("  USEF1                : Use all f's instead of computing f(1) using n2 and n1\n");
	printf("  IN=filename          : specify alternate input file instead of gencaph1.txt\n");
	printf("  OUT=filename         : specify alternate output file instead of gencaph1.out\n");
    for (i=1,determ=1; i<argc; i++) {
        if (strstr(argv[i],"NOPROMPT")!=0) promptflg=0;
        if (strstr(argv[i],"JOLLY")!=0) jlyflg=1;
        if (strstr(argv[i],"CAPTURE")!=0) capflg=1;
        if (strstr(argv[i],"MARK")!=0) mrkflg=1;
        if (strstr(argv[i],"EXP")!=0) determ=1;
        if (strstr(argv[i],"DETERM")!=0) determ=1;
        if (strstr(argv[i],"STOCHASTIC")!=0) determ=0;
        if (strstr(argv[i],"RECRUIT=BINOMIAL")!=0) recruit=1;
        if (strstr(argv[i],"RECRUIT=NORMAL")!=0) recruit=2;
        if (strstr(argv[i],"USEF1")!=0) USEF1=1;
        if (strstr(argv[i],"GROUPS=")!=0) { mxgrps=atol(argv[i]+7); printf("max groups=%ld\n",mxgrps);}
        if (strstr(argv[i],"THETA=TIME")!=0) thetaflg=1;
        if (strstr(argv[i],"IN=")!=0) strcpy(infile,argv[i]+3);
        if (strstr(argv[i],"OUT=")!=0) strcpy(outfile,argv[i]+4);
        if (strstr(argv[i],"SEED=")!=0) { seed=atol(argv[i]+5); printf("SEED=%ld\n",seed);}    
        printf("@@%s\n",argv[i]);
    }
    printf("opening %s...\n",infile);
    if (strcmp(infile,"stdin")==0) F=stdin;
    else {
      if ((F=fopen(infile,"r"))==NULL) { 
        printf(" ** error opening %s\n\n",infile); 
#ifdef DOS
        printf("hit <Enter>...");
        while( !_kbhit() );
#endif
        exit(1);
      }
    }
    if (seed==12345) srand( (unsigned)time( NULL ) ); else srand(seed);
    if (strcmp(outfile,"stdout")==0) fp=stdout;
    else fp=fopen(outfile,"w");  
    fscanf(F,"%d",&nper); 
    dim=1<<nper; 
    if (nper>MXPER) {
        printf("\nERROR: This program is limited to %d capture occasions!\n",MXPER);
        exit(1);
    }
    if (determ) printf("\nExpected value data generated\n"); else printf("\nStochastic data generated\n");
    if (jlyflg==1) fprintf(fp,"TITLE=SIMULATED DATA\nPERIODS=%d\nFORMAT=(%dI1,I6)\n",nper,nper); 
    if (capflg==2) 
    { fprintf(fp,"title='simulated data'\ntask read captures occasions=%d x matrix\n",nper);
        fprintf(fp,"format='(a8,1x,%df1.0)'\nread input data\n",nper);
    }
    if (mrkflg==1) {
        fprintf(fp,"proc title exp value data w/ het.;\n");
        fprintf(fp,"proc chmatrix occasions=%d groups=1 etype=Live hist=%d nohist;\n",nper,dim);
        fprintf(fp,"   glabel(1)=Group 1;\n");
        fprintf(fp,"   time interval "); for (i=0; i<nper-1; i++) fprintf(fp," 1"); fprintf(fp,";\n");
    }
    
    if ((na=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for n[%d]\n",nper); exit(1);}
    if ((tcap=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for tcap[%d]\n",nper); exit(1);}
    if ((talive=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for talive[%d]\n",nper); exit(1);}
    if ((pa=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for p[%d]\n",nper); exit(1);}
    if ((p2=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for p[%d]\n",nper); exit(1);}
    if ((sa=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for s[%d]\n",nper); exit(1);}
    if ((fa=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for f[%d]\n",nper); exit(1);}
    if ((sef=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for sef[%d]\n",nper); exit(1);}
    if ((th=malloc(nper*sizeof(double)))==NULL) { printf("error allocating memory for th[%d]\n",nper); exit(1);}
    for (i=0; i<nper-1;i++) th[i]=1;
    
    if ((nx=malloc(mxgrps*sizeof(double*)))==NULL) { printf("error allocating memory for nx\n"); exit(1); }
    for (i=0; i<mxgrps; i++)
      if ((nx[i]=malloc(dim*sizeof(double)))==NULL) { printf("error allocating memory for nx[%d]\n",dim); exit(1);}
    //printf("\nHit Enter...");  while( !_kbhit() );
    for (tst=1; tst>0;) {
        for (sumn=0,i=0; i<nper; sumn+=na[i++]) {
            na[i]=tcap[i]=talive[i]=0;
            if ((tst=fscanf(F,"%lf",&na[i]))<1) break;
        }
        if (tst<1 || sumn<=0) break;
        printf("n="); for (i=0; i<nper; i++) printf(" %lf",na[i]); printf("\n");
        for (i=0; i<nper-1; i++) fscanf(F,"%lf",&sa[i]);
        printf("srv="); for (i=0; i<nper-1; i++) printf(" %lf",sa[i]); printf("\n"); 
        for (i=0; i<nper-1; i++) fscanf(F,"%lf",&th[i]);
        printf("theta="); for (i=0; i<nper-1; i++) printf(" %lf",th[i]); printf("\n"); 
        for (i=0; i<nper; i++) fscanf(F,"%lf",&pa[i]);
        printf("p="); for (i=0; i<nper; i++) printf(" %lf",pa[i]); printf("\n");
        for (i=0; i<nper; i++) fscanf(F,"%lf",&p2[i]);
        printf("c="); for (i=0; i<nper; i++) printf(" %lf",p2[i]); printf("\n");
        if (recruit>0) {
            for (i=0; i<nper-1; i++) fscanf(F,"%lf",&fa[i]); 
            if (USEF1==0) { fa[0]=na[1]/na[0]-sa[0]; na[0]=na[0]/pa[0];} 
			if (fa[0]<0.) fa[0]=0.;
            printf("f="); for (i=0; i<nper-1; i++) printf(" %lf",fa[i]); printf("\n");
            if (recruit==2) {
                for (i=0; i<nper-1; i++) fscanf(F,"%lf",&sef[i]); 
                printf("se(f)="); sef[0]=sef[1];
                for (i=0; i<nper-1; i++) printf(" %lf",sef[i]); printf("\n");
            }
        }
        for (i=0;i<dim;i++) nx[grp][i]=0; 
        
        for (k=0; k<nper; k++) {
            printf("N[%d]=%f",k,na[k]);talive[k]+=na[k];
            if (na[k]>0) caphis(0,na[k],k+1,-1); printf(" tcap,talive[%d]=%f %f\n",k,tcap[k],talive[k]);
            if (recruit>0) na[k+1]=binoml(talive[k],fa[k]);
            if (recruit==2) na[k+1]=unorm(talive[k]*fa[k],talive[k]*sef[k]);
            if (na[k+1]<0) na[k+1]=0;
        }
        if (mxgrps==1) {
            printf("n="); for (k=0; k<nper; k++) printf(" %f",na[k]); printf("\n");
            for (i=1;i<dim ;i++) if (nx[grp][i]>eps) {
                klmt=floor(.5+nx[grp][i]); if (capflg==0) klmt=1;
                for (k=0; k<klmt; k++) {
                    if (capflg==1) fprintf(fp,"%4.4d%4.4d ",i,k+1);
                    for (j=(nper-1);j>=0;j--) fprintf(fp,"%d",(i&(1<<j))>0); 
                    fprintf(fp," %12.4lf;\n",nx[grp][i]);
                }
            }
        }
        else
            grp++;
    }
    if (mxgrps>1) {
        for (i=1;i<dim ;i++) {
            for (j=0,sm=0; j<mxgrps; j++) sm+=nx[j][i];
            if (sm>eps) {
                for (j=(nper-1);j>=0;j--) fprintf(fp,"%d",(i&(1<<j))>0); 
                for (j=0; j<mxgrps; j++) fprintf(fp," %12.4lf",nx[j][i]);
                fprintf(fp,";\n");
            }
        }
    }
    if (capflg==2)  {
        fprintf(fp,"task closure test\n");
        fprintf(fp,"task model selection\n");
        fprintf(fp,"task population estimate NULL JACKKNIFE MT\n");
    }
    if (mrkflg==1)  {
        fprintf(fp,"proc estimate link=Sin NOLOOP varest=2ndPart ;\n");
        fprintf(fp,"model={Phi(t) p(t) };\n");
        fprintf(fp,"   fixed=1;\n");
        fprintf(fp,"   parm(%d)=1;\n",nper-1);
        fprintf(fp,"   group=1 Phi rows=%d cols=%d Triang;\n",nper-1,nper-1);
        for (j=1; j<nper; j++)  {  for (i=j; i<nper; i++) fprintf(fp,"%4d",i); fprintf(fp,";\n");  }
        fprintf(fp,"   group=1 p rows=%d cols=%d Triang;\n",nper-1,nper-1);
        for (j=1; j<nper; j++)  {  for (i=j; i<nper; i++) fprintf(fp,"%4d",i+nper-1); fprintf(fp,";\n");  }
        fprintf(fp,"   design matrix constraints=%d covariates=%d identity;\n",2*nper-2,2*nper-2);
        for (j=1; j<nper; j++) fprintf(fp,"   blabel(%d)=Phi;\n",j);
        for (j=1; j<nper; j++) fprintf(fp,"   blabel(%d)=p;\n",j+nper-1);
        for (j=1; j<nper; j++) fprintf(fp,"   rlabel(%d)=Phi;\n",j);
        for (j=1; j<nper; j++) fprintf(fp,"   rlabel(%d)=p;\n",j+nper-1);
        fprintf(fp,"proc stop;\n");
    }
    i=fclose(fp); fclose(F);  printf("\noutput saved to 'gencaph1.out'\n");
    if (promptflg>0) { 
#ifdef DOS
		printf("\nHit Enter...");  while( !_kbhit() );
#endif
	}
}