[BACK]Return to imach.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.216 and 1.217

version 1.216, 2015/12/18 17:32:11 version 1.217, 2015/12/23 17:18:31
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.217  2015/12/23 17:18:31  brouard
     Summary: Experimental backcast
   
   Revision 1.216  2015/12/18 17:32:11  brouard    Revision 1.216  2015/12/18 17:32:11  brouard
   Summary: 0.98r4 Warning and status=-2    Summary: 0.98r4 Warning and status=-2
   
Line 841  double **matprod2(); /* test */ Line 844  double **matprod2(); /* test */
 double **oldm, **newm, **savm; /* Working pointers to matrices */  double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 /*FILE *fic ; */ /* Used in readdata only */  /*FILE *fic ; */ /* Used in readdata only */
 FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;  FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficresplb,*ficrespij, *ficrespijb, *ficrest,*ficresf, *ficresfb,*ficrespop;
 FILE *ficlog, *ficrespow;  FILE *ficlog, *ficrespow;
 int globpr=0; /* Global variable for printing or not */  int globpr=0; /* Global variable for printing or not */
 double fretone; /* Only one call to likelihood */  double fretone; /* Only one call to likelihood */
Line 864  char fileresv[FILENAMELENGTH]; Line 867  char fileresv[FILENAMELENGTH];
 FILE  *ficresvpl;  FILE  *ficresvpl;
 char fileresvpl[FILENAMELENGTH];  char fileresvpl[FILENAMELENGTH];
 char title[MAXLINE];  char title[MAXLINE];
 char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];  char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH],  fileresplb[FILENAMELENGTH];
 char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];  char plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
 char tmpout[FILENAMELENGTH],  tmpout2[FILENAMELENGTH];   char tmpout[FILENAMELENGTH],  tmpout2[FILENAMELENGTH]; 
 char command[FILENAMELENGTH];  char command[FILENAMELENGTH];
 int  outcmd=0;  int  outcmd=0;
   
 char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filerespijb[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
 char fileresu[FILENAMELENGTH]; /* fileres without r in front */  char fileresu[FILENAMELENGTH]; /* fileres without r in front */
 char filelog[FILENAMELENGTH]; /* Log file */  char filelog[FILENAMELENGTH]; /* Log file */
 char filerest[FILENAMELENGTH];  char filerest[FILENAMELENGTH];
Line 2141  Earliest age to start was %d-%d=%d, ncvl Line 2144  Earliest age to start was %d-%d=%d, ncvl
   return prlim; /* should not reach here */    return prlim; /* should not reach here */
 }  }
   
   
    /**** Back Prevalence limit (stable or period prevalence)  ****************/
   
   double **bprevalim(double **bprlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int *ncvyear, int ij)
   {
     /* Computes the prevalence limit in each live state at age x by left multiplying the unit
        matrix by transitions matrix until convergence is reached with precision ftolpl */
     /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1  = Wx-n Px-n ... Px-2 Px-1 I */
     /* Wx is row vector: population in state 1, population in state 2, population dead */
     /* or prevalence in state 1, prevalence in state 2, 0 */
     /* newm is the matrix after multiplications, its rows are identical at a factor */
     /* Initial matrix pimij */
     /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */
     /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */
     /*  0,                   0                  , 1} */
     /*
      * and after some iteration: */
     /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */
     /*  0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */
     /*  0,                   0                  , 1} */
     /* And prevalence by suppressing the deaths are close to identical rows in prlim: */
     /* {0.51571254859325999, 0.4842874514067399, */
     /*  0.51326036147820708, 0.48673963852179264} */
     /* If we start from prlim again, prlim tends to a constant matrix */
   
     int i, ii,j,k;
     double *min, *max, *meandiff, maxmax,sumnew=0.;
     /* double **matprod2(); */ /* test */
     double **out, cov[NCOVMAX+1], **bmij();
     double **newm;
     double agefin, delaymax=200. ; /* 100 Max number of years to converge */
     int ncvloop=0;
     
     min=vector(1,nlstate);
     max=vector(1,nlstate);
     meandiff=vector(1,nlstate);
   
     for (ii=1;ii<=nlstate+ndeath;ii++)
       for (j=1;j<=nlstate+ndeath;j++){
         oldm[ii][j]=(ii==j ? 1.0 : 0.0);
       }
     
     cov[1]=1.;
     
     /* Even if hstepm = 1, at least one multiplication by the unit matrix */
     /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
     for(agefin=age+stepm/YEARM; agefin<=age+delaymax; agefin=agefin+stepm/YEARM){
       ncvloop++;
       newm=savm;
       /* Covariates have to be included here again */
       cov[2]=agefin;
       if(nagesqr==1)
         cov[3]= agefin*agefin;;
       for (k=1; k<=cptcovn;k++) {
         /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
         cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
         /* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */
       }
       /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
       /* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]*cov[2]; */
       for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
       for (k=1; k<=cptcovprod;k++) /* Useless */
         /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
         cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
       
       /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
       /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
       /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/
       /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
       /* out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /\* Bug Valgrind *\/ */
       out=matprod2(newm, oldm ,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, bmij(pmmij,cov,ncovmodel,x,nlstate)); /* Bug Valgrind */
       
       savm=oldm;
       oldm=newm;
   
       for(j=1; j<=nlstate; j++){
         max[j]=0.;
         min[j]=1.;
       }
         /* sumnew=0; */
         /* for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k]; */
       for(j=1; j<=nlstate; j++){ 
         for(i=1;i<=nlstate;i++){
           /* bprlim[i][j]= newm[i][j]/(1-sumnew); */
           bprlim[i][j]= newm[i][j];
           max[i]=FMAX(max[i],bprlim[i][j]);
           min[i]=FMIN(min[i],bprlim[i][j]);
         }
       }
   
       maxmax=0.;
       for(i=1; i<=nlstate; i++){
         meandiff[i]=(max[i]-min[i])/(max[i]+min[i])*2.; /* mean difference for each column */
         maxmax=FMAX(maxmax,meandiff[i]);
         /* printf("Back age= %d meandiff[%d]=%f, agefin=%d max[%d]=%f min[%d]=%f maxmax=%f\n", (int)age, i, meandiff[i],(int)agefin, i, max[i], i, min[i],maxmax); */
       } /* j loop */
       *ncvyear= -( (int)age- (int)agefin);
       /* printf("Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear); */
       if(maxmax < ftolpl){
         printf("OK Back maxmax=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, ncvloop, (int)age, (int)agefin, *ncvyear);
         free_vector(min,1,nlstate);
         free_vector(max,1,nlstate);
         free_vector(meandiff,1,nlstate);
         return bprlim;
       }
     } /* age loop */
       /* After some age loop it doesn't converge */
     printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
   Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
     /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
     free_vector(min,1,nlstate);
     free_vector(max,1,nlstate);
     free_vector(meandiff,1,nlstate);
     
     return bprlim; /* should not reach here */
   }
   
 /*************** transition probabilities ***************/   /*************** transition probabilities ***************/ 
   
 double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )  double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
Line 2223  double **pmij(double **ps, double *cov, Line 2343  double **pmij(double **ps, double *cov,
     return ps;      return ps;
 }  }
   
   /*************** transition probabilities ***************/ 
   
   double **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
   {
     /* According to parameters values stored in x and the covariate's values stored in cov,
        computes the probability to be observed in state j being in state i by appying the
        model to the ncovmodel covariates (including constant and age).
        lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc]
        and, according on how parameters are entered, the position of the coefficient xij(nc) of the
        ncth covariate in the global vector x is given by the formula:
        j<i nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel
        j>=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel
        Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation,
        sums on j different of i to get 1-pii/pii, deduces pii, and then all pij.
        Outputs ps[i][j] the probability to be observed in j being in j according to
        the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij]
     */
     double s1, lnpijopii;
     /*double t34;*/
     int i,j, nc, ii, jj;
   
       for(i=1; i<= nlstate; i++){
         for(j=1; j<i;j++){
           for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
             /*lnpijopii += param[i][j][nc]*cov[nc];*/
             lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
   /*       printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
           }
           ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
   /*      printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
         }
         for(j=i+1; j<=nlstate+ndeath;j++){
           for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
             /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
             lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
   /*        printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
           }
           ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
         }
       }
       
       for(i=1; i<= nlstate; i++){
         s1=0;
         for(j=1; j<i; j++){
           s1+=exp(ps[i][j]); /* In fact sums pij/pii */
           /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
         }
         for(j=i+1; j<=nlstate+ndeath; j++){
           s1+=exp(ps[i][j]); /* In fact sums pij/pii */
           /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
         }
         /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
         ps[i][i]=1./(s1+1.);
         /* Computing other pijs */
         for(j=1; j<i; j++)
           ps[i][j]= exp(ps[i][j])*ps[i][i];
         for(j=i+1; j<=nlstate+ndeath; j++)
           ps[i][j]= exp(ps[i][j])*ps[i][i];
         /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
       } /* end i */
       
       for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
         for(jj=1; jj<= nlstate+ndeath; jj++){
           ps[ii][jj]=0;
           ps[ii][ii]=1;
         }
       }
       /* Added for backcast */
       for(jj=1; jj<= nlstate; jj++){
         s1=0.;
         for(ii=1; ii<= nlstate; ii++){
           s1+=ps[ii][jj];
         }
         for(ii=1; ii<= nlstate; ii++){
           ps[ii][jj]=ps[ii][jj]/s1;
         }
       }
        
       /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
       /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
       /*  printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
       /*   } */
       /*   printf("\n "); */
       /* } */
       /* printf("\n ");printf("%lf ",cov[2]);*/
       /*
         for(i=1; i<= npar; i++) printf("%f ",x[i]);
         goto end;*/
       return ps;
   }
   
   
 /**************** Product of 2 matrices ******************/  /**************** Product of 2 matrices ******************/
   
 double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b)  double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b)
Line 2297  double ***hpxij(double ***po, int nhstep Line 2509  double ***hpxij(double ***po, int nhstep
       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/        /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
       out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,         out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, 
                    pmij(pmmij,cov,ncovmodel,x,nlstate));                     pmij(pmmij,cov,ncovmodel,x,nlstate));
         /* if((int)age == 70){ */
         /*        printf(" Forward hpxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
         /*        for(i=1; i<=nlstate+ndeath; i++) { */
         /*          printf("%d pmmij ",i); */
         /*          for(j=1;j<=nlstate+ndeath;j++) { */
         /*            printf("%f ",pmmij[i][j]); */
         /*          } */
         /*          printf(" oldm "); */
         /*          for(j=1;j<=nlstate+ndeath;j++) { */
         /*            printf("%f ",oldm[i][j]); */
         /*          } */
         /*          printf("\n"); */
         /*        } */
         /* } */
       savm=oldm;        savm=oldm;
       oldm=newm;        oldm=newm;
     }      }
Line 2311  double ***hpxij(double ***po, int nhstep Line 2537  double ***hpxij(double ***po, int nhstep
   return po;    return po;
 }  }
   
   /************* Higher Back Matrix Product ***************/
   
   double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij )
   {
     /* Computes the transition matrix starting at age 'age' over 
        'nhstepm*hstepm*stepm' months (i.e. until
        age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying 
        nhstepm*hstepm matrices. 
        Output is stored in matrix po[i][j][h] for h every 'hstepm' step 
        (typically every 2 years instead of every month which is too big 
        for the memory).
        Model is determined by parameters x and covariates have to be 
        included manually here. 
   
        */
   
     int i, j, d, h, k;
     double **out, cov[NCOVMAX+1];
     double **newm;
     double agexact;
     double agebegin, ageend;
   
     /* Hstepm could be zero and should return the unit matrix */
     for (i=1;i<=nlstate+ndeath;i++)
       for (j=1;j<=nlstate+ndeath;j++){
         oldm[i][j]=(i==j ? 1.0 : 0.0);
         po[i][j][0]=(i==j ? 1.0 : 0.0);
       }
     /* Even if hstepm = 1, at least one multiplication by the unit matrix */
     for(h=1; h <=nhstepm; h++){
       for(d=1; d <=hstepm; d++){
         newm=savm;
         /* Covariates have to be included here again */
         cov[1]=1.;
         agexact=age-((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
         /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */
         cov[2]=agexact;
         if(nagesqr==1)
           cov[3]= agexact*agexact;
         for (k=1; k<=cptcovn;k++) 
           cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];
           /* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */
         for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
           /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
           cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
           /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */
         for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
           cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
           /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
   
   
         /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
         /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
         out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, 
                     oldm);
         /* if((int)age == 70){ */
         /*        printf(" Backward hbxij age=%d agexact=%f d=%d nhstepm=%d hstepm=%d\n", (int) age, agexact, d, nhstepm, hstepm); */
         /*        for(i=1; i<=nlstate+ndeath; i++) { */
         /*          printf("%d pmmij ",i); */
         /*          for(j=1;j<=nlstate+ndeath;j++) { */
         /*            printf("%f ",pmmij[i][j]); */
         /*          } */
         /*          printf(" oldm "); */
         /*          for(j=1;j<=nlstate+ndeath;j++) { */
         /*            printf("%f ",oldm[i][j]); */
         /*          } */
         /*          printf("\n"); */
         /*        } */
         /* } */
         savm=oldm;
         oldm=newm;
       }
       for(i=1; i<=nlstate+ndeath; i++)
         for(j=1;j<=nlstate+ndeath;j++) {
           po[i][j][h]=newm[i][j];
           /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
         }
       /*printf("h=%d ",h);*/
     } /* end h */
   /*     printf("\n H=%d \n",h); */
     return po;
   }
   
   
 #ifdef NLOPT  #ifdef NLOPT
   double  myfunc(unsigned n, const double *p1, double *grad, void *pd){    double  myfunc(unsigned n, const double *p1, double *grad, void *pd){
   double fret;    double fret;
Line 2694  double funcone( double *x) Line 3004  double funcone( double *x)
               
       s1=s[mw[mi][i]][i];        s1=s[mw[mi][i]][i];
       s2=s[mw[mi+1][i]][i];        s2=s[mw[mi+1][i]][i];
       if(s2==-1){        /* if(s2==-1){ */
         printf(" s1=%d, s2=%d i=%d \n", s1, s2, i);        /*        printf(" s1=%d, s2=%d i=%d \n", s1, s2, i); */
         /* exit(1); */        /*        /\* exit(1); *\/ */
       }        /* } */
       bbh=(double)bh[mi][i]/(double)stepm;         bbh=(double)bh[mi][i]/(double)stepm; 
       /* bias is positive if real duration        /* bias is positive if real duration
        * is higher than the multiple of stepm and negative otherwise.         * is higher than the multiple of stepm and negative otherwise.
Line 3614  void  concatwav(int wav[], int **dh, int Line 3924  void  concatwav(int wav[], int **dh, int
   int i, mi, m;    int i, mi, m;
   /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;    /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
      double sum=0., jmean=0.;*/       double sum=0., jmean=0.;*/
   int first, firstwo;    int first, firstwo, firsthree;
   int j, k=0,jk, ju, jl;    int j, k=0,jk, ju, jl;
   double sum=0.;    double sum=0.;
   first=0;    first=0;
   firstwo=0;    firstwo=0;
     firsthree=0;
   jmin=100000;    jmin=100000;
   jmax=-1;    jmax=-1;
   jmean=0.;    jmean=0.;
Line 3631  void  concatwav(int wav[], int **dh, int Line 3942  void  concatwav(int wav[], int **dh, int
       }        }
       if(m >=lastpass){        if(m >=lastpass){
         if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){          if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
           printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);            if(firsthree == 0){
           fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);              printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
               fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
               firsthree=1;
             }else{
               fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
             }
           mw[++mi][i]=m;            mw[++mi][i]=m;
         }          }
         if(s[m][i]==-2){ /* Vital status is really unknown */          if(s[m][i]==-2){ /* Vital status is really unknown */
Line 4998  To be simple, these graphs help to under Line 5314  To be simple, these graphs help to under
 void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \  void printinghtml(char fileresu[], char title[], char datafile[], int firstpass, \
                   int lastpass, int stepm, int weightopt, char model[],\                    int lastpass, int stepm, int weightopt, char model[],\
                   int imx,int jmin, int jmax, double jmeanint,char rfileres[],\                    int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
                   int popforecast, int prevfcast, int estepm ,          \                    int popforecast, int prevfcast, int backcast, int estepm , \
                   double jprev1, double mprev1,double anprev1, double dateprev1, \                    double jprev1, double mprev1,double anprev1, double dateprev1, \
                   double jprev2, double mprev2,double anprev2, double dateprev2){                    double jprev2, double mprev2,double anprev2, double dateprev2){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
Line 5016  void printinghtml(char fileresu[], char Line 5332  void printinghtml(char fileresu[], char
  - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",   - Estimated transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
            stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_"));             stepm,subdirf2(fileresu,"PIJ_"),subdirf2(fileresu,"PIJ_"));
    fprintf(fichtm,"\     fprintf(fichtm,"\
    - Estimated back transition probabilities over %d (stepm) months: <a href=\"%s\">%s</a><br>\n ",
              stepm,subdirf2(fileresu,"PIJB_"),subdirf2(fileresu,"PIJB_"));
      fprintf(fichtm,"\
  - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",   - Period (stable) prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
            subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));             subdirf2(fileresu,"PL_"),subdirf2(fileresu,"PL_"));
    fprintf(fichtm,"\     fprintf(fichtm,"\
    - Period (stable) back prevalence in each health state: <a href=\"%s\">%s</a> <br>\n",
              subdirf2(fileresu,"PLB_"),subdirf2(fileresu,"PLB_"));
      fprintf(fichtm,"\
  - (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \   - (a) Life expectancies by health status at initial age, e<sub>i.</sub> (b) health expectancies by health status at initial age, e<sub>ij</sub> . If one or more covariates are included, specific tables for each value of the covariate are output in sequences within the same file (estepm=%2d months): \
    <a href=\"%s\">%s</a> <br>\n",     <a href=\"%s\">%s</a> <br>\n",
            estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));             estepm,subdirf2(fileresu,"E_"),subdirf2(fileresu,"E_"));
Line 5070  divided by h: <sub>h</sub>P<sub>ij</sub> Line 5392  divided by h: <sub>h</sub>P<sub>ij</sub>
      }       }
      /* Period (stable) prevalence in each health state */       /* Period (stable) prevalence in each health state */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \         fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);  <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
      }       }
       if(backcast==1){
        /* Period (stable) back prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){
          fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \
   <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1);
        }
       }
     if(prevfcast==1){      if(prevfcast==1){
       /* Projection of prevalence up to period (stable) prevalence in each health state */        /* Projection of prevalence up to period (stable) prevalence in each health state */
       for(cpt=1; cpt<=nlstate;cpt++){        for(cpt=1; cpt<=nlstate;cpt++){
Line 5360  unset log y\n\ Line 5689  unset log y\n\
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  plot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3;        k=3;
       for (i=1; i<= nlstate ; i ++){        for (i=1; i<= nlstate ; i ++){
         if(i==1)          if(i==1){
           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));            fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
         else          }else{
           fprintf(ficgp,", '' ");            fprintf(ficgp,", '' ");
           }
         l=(nlstate+ndeath)*(i-1)+1;          l=(nlstate+ndeath)*(i-1)+1;
         fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);          fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
         for (j=2; j<= nlstate+ndeath ; j ++)          for (j=2; j<= nlstate+ndeath ; j ++)
Line 5453  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 5783  plot [%.f:%.f]  ", ageminpar, agemaxpar)
     } /* end cpt state*/       } /* end cpt state*/ 
   } /* end covariate */      } /* end covariate */  
   
       /* CV back preval stable (period) for each covariate */
     for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
       for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
         fprintf(ficgp,"\n#\n#\n#CV Back preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
         for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
           lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[lv]][lv];
           fprintf(ficgp," V%d=%d ",k,vlv);
         }
         fprintf(ficgp,"\n#\n");
   
         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
   set ter svg size 640, 480\n\
   unset log y\n\
   plot [%.f:%.f]  ", ageminpar, agemaxpar);
         k=3; /* Offset */
         for (i=1; i<= nlstate ; i ++){
           if(i==1)
             fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
           else
             fprintf(ficgp,", '' ");
           l=(nlstate+ndeath)*(i-1)+1;
           fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /* a vérifier */
           for (j=2; j<= nlstate ; j ++)
             fprintf(ficgp,"+$%d",k+l+j-1);
           fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
         } /* nlstate */
         fprintf(ficgp,"\nset out\n");
       } /* end cpt state*/ 
     } /* end covariate */  
   
   if(prevfcast==1){    if(prevfcast==1){
   /* Projection from cross-sectional to stable (period) for each covariate */    /* Projection from cross-sectional to stable (period) for each covariate */
   
Line 5620  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 5985  plot [%.f:%.f]  ", ageminpar, agemaxpar)
                else                 else
                  fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);                   fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
              }               }
              if(ng != 1){             }else{
                fprintf(ficgp,")/(1");               i=i-ncovmodel;
                if(ng !=1 ) /* For logit formula of log p11 is more difficult to get */
                  fprintf(ficgp," (1.");
              }
              
              if(ng != 1){
                fprintf(ficgp,")/(1");
                             
                for(k1=1; k1 <=nlstate; k1++){                for(k1=1; k1 <=nlstate; k1++){ 
                  if(nagesqr==0)                 if(nagesqr==0)
                    fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);                   fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
                  else /* nagesqr =1 */                 else /* nagesqr =1 */
                    fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);                   fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr);
                                   
                  ij=1;                 ij=1;
                  for(j=3; j <=ncovmodel-nagesqr; j++){                 for(j=3; j <=ncovmodel-nagesqr; j++){
                    if(ij <=cptcovage) { /* Bug valgrind */                   if(ij <=cptcovage) { /* Bug valgrind */
                      if((j-2)==Tage[ij]) { /* Bug valgrind */                     if((j-2)==Tage[ij]) { /* Bug valgrind */
                        fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);                       fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
                        /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */                       /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); */
                        ij++;                       ij++;
                      }  
                    }                     }
                    else  
                      fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);  
                  }                   }
                  fprintf(ficgp,")");                   else
                      fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]);
                }                 }
                fprintf(ficgp,")");                 fprintf(ficgp,")");
                if(ng ==2)  
                  fprintf(ficgp," t \"p%d%d\" ", k2,k);  
                else /* ng= 3 */  
                  fprintf(ficgp," t \"i%d%d\" ", k2,k);  
              }else{ /* end ng <> 1 */  
                fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);  
              }               }
              if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");               fprintf(ficgp,")");
              i=i+ncovmodel;               if(ng ==2)
                  fprintf(ficgp," t \"p%d%d\" ", k2,k);
                else /* ng= 3 */
                  fprintf(ficgp," t \"i%d%d\" ", k2,k);
              }else{ /* end ng <> 1 */
                if( k !=k2) /* logit p11 is hard to draw */
                  fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);
            }             }
              if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
                fprintf(ficgp,",");
              if (ng == 1 && k!=k2 && (k+k2)!= (nlstate*2+ndeath))
                fprintf(ficgp,",");
              i=i+ncovmodel;
          } /* end k */           } /* end k */
        } /* end k2 */         } /* end k2 */
        fprintf(ficgp,"\n set out\n");         fprintf(ficgp,"\n set out\n");
Line 5788  void prevforecast(char fileres[], double Line 6162  void prevforecast(char fileres[], double
           fprintf(ficresf," p%d%d",i,j);            fprintf(ficresf," p%d%d",i,j);
         fprintf(ficresf," p.%d",j);          fprintf(ficresf," p.%d",j);
       }        }
       for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {         for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
         fprintf(ficresf,"\n");          fprintf(ficresf,"\n");
         fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);             fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);   
   
         for (agec=fage; agec>=(ageminpar-1); agec--){           for (agec=fage; agec>=(ageminpar-1); agec--){ 
           nhstepm=(int) rint((agelim-agec)*YEARM/stepm);             nhstepm=(int) rint((agelim-agec)*YEARM/stepm); 
           nhstepm = nhstepm/hstepm;             nhstepm = nhstepm/hstepm; 
           p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);            p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           oldm=oldms;savm=savms;            oldm=oldms;savm=savms;
           hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);              hpxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);
                   
           for (h=0; h<=nhstepm; h++){            for (h=0; h<=nhstepm; h++){
             if (h*hstepm/YEARM*stepm ==yearp) {              if (h*hstepm/YEARM*stepm ==yearp) {
Line 5837  void prevforecast(char fileres[], double Line 6210  void prevforecast(char fileres[], double
   
 }  }
   
   /************** Back Forecasting ******************/
   void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
     /* back1, year, month, day of starting backection 
        agemin, agemax range of age
        dateprev1 dateprev2 range of dates during which prevalence is computed
        anback2 year of en of backection (same day and month as back1).
     */
     int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1;
     double agec; /* generic age */
     double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
     double *popeffectif,*popcount;
     double ***p3mat;
     double ***mobaverage;
     char fileresfb[FILENAMELENGTH];
   
     agelim=AGESUP;
     /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
        in each health status at the date of interview (if between dateprev1 and dateprev2).
        We still use firstpass and lastpass as another selection.
     */
     /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
     /*          firstpass, lastpass,  stepm,  weightopt, model); */
     prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
    
     strcpy(fileresfb,"FB_"); 
     strcat(fileresfb,fileresu);
     if((ficresfb=fopen(fileresfb,"w"))==NULL) {
       printf("Problem with back forecast resultfile: %s\n", fileresfb);
       fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb);
     }
     printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb);
     fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb);
   
     if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
     if (mobilav!=0) {
       mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
       if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
     }
   
     stepsize=(int) (stepm+YEARM-1)/YEARM;
     if (stepm<=12) stepsize=1;
     if(estepm < stepm){
       printf ("Problem %d lower than %d\n",estepm, stepm);
     }
     else  hstepm=estepm;   
   
     hstepm=hstepm/stepm; 
     yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and
                                  fractional in yp1 */
     anprojmean=yp;
     yp2=modf((yp1*12),&yp);
     mprojmean=yp;
     yp1=modf((yp2*30.5),&yp);
     jprojmean=yp;
     if(jprojmean==0) jprojmean=1;
     if(mprojmean==0) jprojmean=1;
   
     i1=cptcoveff;
     if (cptcovn < 1){i1=1;}
     
     fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2); 
     
     fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
   
   /*            if (h==(int)(YEARM*yearp)){ */
     for(cptcov=1, k=0;cptcov<=i1;cptcov++){
       for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
         k=k+1;
         fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#");
         for(j=1;j<=cptcoveff;j++) {
           fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
         }
         fprintf(ficresfb," yearbproj age");
         for(j=1; j<=nlstate+ndeath;j++){ 
           for(i=1; i<=nlstate;i++)              
             fprintf(ficresfb," p%d%d",i,j);
           fprintf(ficresfb," p.%d",j);
         }
         for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) { 
        /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {  */
           fprintf(ficresfb,"\n");
           fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);   
           for (agec=fage; agec>=(ageminpar-1); agec--){ 
             nhstepm=(int) rint((agelim-agec)*YEARM/stepm); 
             nhstepm = nhstepm/hstepm; 
             p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
             oldm=oldms;savm=savms;
             hbxij(p3mat,nhstepm,agec,hstepm,p,nlstate,stepm,oldm,savm, k);  
           
             for (h=0; h<=nhstepm; h++){
               if (h*hstepm/YEARM*stepm ==yearp) {
                 fprintf(ficresfb,"\n");
                 for(j=1;j<=cptcoveff;j++) 
                   fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
                 fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm);
               } 
               for(j=1; j<=nlstate+ndeath;j++) {
                 ppij=0.;
                 for(i=1; i<=nlstate;i++) {
                   if (mobilav==1) 
                     ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
                   else {
                     ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
                   }
                   if (h*hstepm/YEARM*stepm== yearp) {
                     fprintf(ficresfb," %.3f", p3mat[i][j][h]);
                   }
                 } /* end i */
                 if (h*hstepm/YEARM*stepm==yearp) {
                   fprintf(ficresfb," %.3f", ppij);
                 }
               }/* end j */
             } /* end h */
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           } /* end agec */
         } /* end yearp */
       } /* end cptcod */
     } /* end  cptcov */
          
     if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
     fclose(ficresfb);
     printf("End of Computing Back forecasting \n");
     fprintf(ficlog,"End of Computing Back forecasting\n");
   
   }
   
 /************** Forecasting *****not tested NB*************/  /************** Forecasting *****not tested NB*************/
 void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){  void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
       
Line 7020  void syscompilerinfo(int logged) Line 7524  void syscompilerinfo(int logged)
         return 0;          return 0;
 }  }
   
    int back_prevalence_limit(double *p, double **bprlim, double ageminpar, double agemaxpar, double ftolpl, int *ncvyearp){
     /*--------------- Back Prevalence limit  (period or stable prevalence) --------------*/
     int i, j, k, i1 ;
     /* double ftolpl = 1.e-10; */
     double age, agebase, agelim;
     double tot;
   
     strcpy(fileresplb,"PLB_");
     strcat(fileresplb,fileresu);
     if((ficresplb=fopen(fileresplb,"w"))==NULL) {
       printf("Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1;
       fprintf(ficlog,"Problem with period (stable) back prevalence resultfile: %s\n", fileresplb);return 1;
     }
     printf("Computing period (stable) back prevalence: result on file '%s' \n", fileresplb);
     fprintf(ficlog,"Computing period (stable) back prevalence: result on file '%s' \n", fileresplb);
     pstamp(ficresplb);
     fprintf(ficresplb,"# Period (stable) back prevalence. Precision given by ftolpl=%g \n", ftolpl);
     fprintf(ficresplb,"#Age ");
     for(i=1; i<=nlstate;i++) fprintf(ficresplb,"%d-%d ",i,i);
     fprintf(ficresplb,"\n");
     
       /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
   
       agebase=ageminpar;
       agelim=agemaxpar;
   
       i1=pow(2,cptcoveff);
       if (cptcovn < 1){i1=1;}
   
       for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
         //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
           k=k+1;
           /* to clean */
           //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov));
           fprintf(ficresplb,"#******");
           printf("#******");
           fprintf(ficlog,"#******");
           for(j=1;j<=cptcoveff;j++) {
             fprintf(ficresplb," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
             printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
             fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
           }
           fprintf(ficresplb,"******\n");
           printf("******\n");
           fprintf(ficlog,"******\n");
   
           fprintf(ficresplb,"#Age ");
           for(j=1;j<=cptcoveff;j++) {
             fprintf(ficresplb,"V%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
           }
           for(i=1; i<=nlstate;i++) fprintf(ficresplb,"  %d-%d   ",i,i);
           fprintf(ficresplb,"Total Years_to_converge\n");
           
           for (age=agebase; age<=agelim; age++){
           /* for (age=agebase; age<=agebase; age++){ */
             bprevalim(bprlim, nlstate, p, age, oldm, savm, ftolpl, ncvyearp, k);
             fprintf(ficresplb,"%.0f ",age );
             for(j=1;j<=cptcoveff;j++)
               fprintf(ficresplb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
             tot=0.;
             for(i=1; i<=nlstate;i++){
               tot +=  bprlim[i][i];
               fprintf(ficresplb," %.5f", bprlim[i][i]);
             }
             fprintf(ficresplb," %.3f %d\n", tot, *ncvyearp);
           } /* Age */
           /* was end of cptcod */
       } /* cptcov */
           return 0;
   }
   
 int hPijx(double *p, int bage, int fage){  int hPijx(double *p, int bage, int fage){
     /*------------- h Pij x at various ages ------------*/      /*------------- h Pij x at various ages ------------*/
   
Line 7090  int hPijx(double *p, int bage, int fage) Line 7666  int hPijx(double *p, int bage, int fage)
         return 0;          return 0;
 }  }
   
    int hBijx(double *p, int bage, int fage){
       /*------------- h Bij x at various ages ------------*/
   
     int stepsize;
     int agelim;
     int hstepm;
     int nhstepm;
     int h, i, i1, j, k;
   
     double agedeb;
     double ***p3mat;
   
       strcpy(filerespijb,"PIJB_");  strcat(filerespijb,fileresu);
       if((ficrespijb=fopen(filerespijb,"w"))==NULL) {
         printf("Problem with Pij back resultfile: %s\n", filerespijb); return 1;
         fprintf(ficlog,"Problem with Pij back resultfile: %s\n", filerespijb); return 1;
       }
       printf("Computing pij back: result on file '%s' \n", filerespijb);
       fprintf(ficlog,"Computing pij back: result on file '%s' \n", filerespijb);
     
       stepsize=(int) (stepm+YEARM-1)/YEARM;
       /*if (stepm<=24) stepsize=2;*/
   
       agelim=AGESUP;
       hstepm=stepsize*YEARM; /* Every year of age */
       hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
   
       /* hstepm=1;   aff par mois*/
       pstamp(ficrespijb);
       fprintf(ficrespijb,"#****** h Pij x Back Probability to be in state i at age x-h being in j at x ");
       i1= pow(2,cptcoveff);
      /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
      /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
      /*   k=k+1;  */
       for (k=1; k <= (int) pow(2,cptcoveff); k++){
         fprintf(ficrespijb,"\n#****** ");
         for(j=1;j<=cptcoveff;j++) 
           fprintf(ficrespijb,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
         fprintf(ficrespijb,"******\n");
         
         /* for (agedeb=fage; agedeb>=bage; agedeb--){ /\* If stepm=6 months *\/ */
         for (agedeb=bage; agedeb<=fage; agedeb++){ /* If stepm=6 months */
           nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
           nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
           
           /*        nhstepm=nhstepm*YEARM; aff par mois*/
           
           p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           oldm=oldms;savm=savms;
           hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
           fprintf(ficrespijb,"# Cov Agex agex-h hpijx with i,j=");
           for(i=1; i<=nlstate;i++)
             for(j=1; j<=nlstate+ndeath;j++)
               fprintf(ficrespijb," %1d-%1d",i,j);
           fprintf(ficrespijb,"\n");
           for (h=0; h<=nhstepm; h++){
             /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
             fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb - h*hstepm/YEARM*stepm );
             /* fprintf(ficrespijb,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); */
             for(i=1; i<=nlstate;i++)
               for(j=1; j<=nlstate+ndeath;j++)
                 fprintf(ficrespijb," %.5f", p3mat[i][j][h]);
             fprintf(ficrespijb,"\n");
           }
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           fprintf(ficrespijb,"\n");
         }
         /*}*/
       }
           return 0;
   }
   
   
 /***********************************************/  /***********************************************/
 /**************** Main Program *****************/  /**************** Main Program *****************/
Line 7141  int main(int argc, char *argv[]) Line 7789  int main(int argc, char *argv[])
   
   int *tab;     int *tab; 
   int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */    int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
     int backcast=0;
   int mobilav=0,popforecast=0;    int mobilav=0,popforecast=0;
   int hstepm=0, nhstepm=0;    int hstepm=0, nhstepm=0;
   int agemortsup;    int agemortsup;
Line 7151  int main(int argc, char *argv[]) Line 7800  int main(int argc, char *argv[])
   double bage=0, fage=110., age, agelim=0., agebase=0.;    double bage=0, fage=110., age, agelim=0., agebase=0.;
   double ftolpl=FTOL;    double ftolpl=FTOL;
   double **prlim;    double **prlim;
     double **bprlim;
   double ***param; /* Matrix of parameters */    double ***param; /* Matrix of parameters */
   double  *p;    double  *p;
   double **matcov; /* Matrix of covariance */    double **matcov; /* Matrix of covariance */
Line 7162  int main(int argc, char *argv[]) Line 7812  int main(int argc, char *argv[])
   double *epj, vepp;    double *epj, vepp;
   
   double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;    double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
     double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000;
   
   double **ximort;    double **ximort;
   char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";    char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
   int *dcwave;    int *dcwave;
Line 8464  Please run with mle=-1 to get a correct Line 9116  Please run with mle=-1 to get a correct
     fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);      fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
     /* day and month of proj2 are not used but only year anproj2.*/      /* day and month of proj2 are not used but only year anproj2.*/
           
       while((c=getc(ficpar))=='#' && c!= EOF){
         ungetc(c,ficpar);
         fgets(line, MAXLINE, ficpar);
         fputs(line,stdout);
         fputs(line,ficparo);
       }
       ungetc(c,ficpar);
       
       fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
       fscanf(ficparo,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
       fscanf(ficlog,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
       fscanf(ficres,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
       /* day and month of proj2 are not used but only year anproj2.*/
           
           
      /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */       /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
Line 8481  Please run with mle=-1 to get a correct Line 9146  Please run with mle=-1 to get a correct
       printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p);        printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, pathc,p);
           
     printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\      printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt,\
                  model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,estepm, \                   model,imx,jmin,jmax,jmean,rfileres,popforecast,prevfcast,backcast, estepm, \
                  jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);                   jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
               
    /*------------ free_vector  -------------*/     /*------------ free_vector  -------------*/
Line 8508  Please run with mle=-1 to get a correct Line 9173  Please run with mle=-1 to get a correct
     prevalence_limit(p, prlim,  ageminpar, agemaxpar, ftolpl, &ncvyear);      prevalence_limit(p, prlim,  ageminpar, agemaxpar, ftolpl, &ncvyear);
     fclose(ficrespl);      fclose(ficrespl);
   
       /*--------------- Back Prevalence limit  (period or stable prevalence) --------------*/
       /*#include "prevlim.h"*/  /* Use ficresplb, ficlog */
       bprlim=matrix(1,nlstate,1,nlstate);
       back_prevalence_limit(p, bprlim,  ageminpar, agemaxpar, ftolpl, &ncvyear);
       fclose(ficresplb);
   
       
 #ifdef FREEEXIT2  #ifdef FREEEXIT2
 #include "freeexit2.h"  #include "freeexit2.h"
 #endif  #endif
Line 8517  Please run with mle=-1 to get a correct Line 9189  Please run with mle=-1 to get a correct
     hPijx(p, bage, fage);      hPijx(p, bage, fage);
     fclose(ficrespij);      fclose(ficrespij);
   
       hBijx(p, bage, fage);
       fclose(ficrespijb);
   
   /*-------------- Variance of one-step probabilities---*/    /*-------------- Variance of one-step probabilities---*/
     k=1;      k=1;
     varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);      varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
Line 8533  Please run with mle=-1 to get a correct Line 9208  Please run with mle=-1 to get a correct
     if(prevfcast==1){      if(prevfcast==1){
       /*    if(stepm ==1){*/        /*    if(stepm ==1){*/
       prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);        prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
       /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/  
       /*      }  */  
       /*      else{ */  
       /*        erreur=108; */  
       /*        printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */  
       /*        fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */  
       /*      } */  
     }      }
       if(backcast==1){
         prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
       }
       /* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
       /*      }  */
       /*      else{ */
       /*        erreur=108; */
       /*        printf("Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
       /*        fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
       /*      } */
       
     
     /* ------ Other prevalence ratios------------ */      /* ------ Other prevalence ratios------------ */
   
Removed from v.1.216  
changed lines
  Added in v.1.217

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>