Diff for /imach/src/imach.c between versions 1.212 and 1.215

version 1.212, 2015/11/21 12:47:24 version 1.215, 2015/12/16 08:52:24
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.215  2015/12/16 08:52:24  brouard
     Summary: 0.98r4 working
   
     Revision 1.214  2015/12/16 06:57:54  brouard
     Summary: temporary not working
   
     Revision 1.213  2015/12/11 18:22:17  brouard
     Summary: 0.98r4
   
   Revision 1.212  2015/11/21 12:47:24  brouard    Revision 1.212  2015/11/21 12:47:24  brouard
   Summary: minor typo    Summary: minor typo
   
Line 824  double **matprod2(); /* test */ Line 833  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, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop;  FILE *ficpar, *ficparo,*ficres, *ficresp, *ficresphtm, *ficresphtmfr, *ficrespl, *ficrespij, *ficrest,*ficresf,*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 1386  char *subdirf3(char fileres[], char *pre Line 1395  char *subdirf3(char fileres[], char *pre
   strcat(tmpout,fileres);    strcat(tmpout,fileres);
   return tmpout;    return tmpout;
 }  }
    
   /*************** function subdirfext ***********/
   char *subdirfext(char fileres[], char *preop, char *postop)
   {
     
     strcpy(tmpout,preop);
     strcat(tmpout,fileres);
     strcat(tmpout,postop);
     return tmpout;
   }
   
   /*************** function subdirfext3 ***********/
   char *subdirfext3(char fileres[], char *preop, char *postop)
   {
     
     /* Caution optionfilefiname is hidden */
     strcpy(tmpout,optionfilefiname);
     strcat(tmpout,"/");
     strcat(tmpout,preop);
     strcat(tmpout,fileres);
     strcat(tmpout,postop);
     return tmpout;
   }
    
 char *asc_diff_time(long time_sec, char ascdiff[])  char *asc_diff_time(long time_sec, char ascdiff[])
 {  {
   long sec_left, days, hours, minutes;    long sec_left, days, hours, minutes;
Line 2223  double ***hpxij(double ***po, int nhstep Line 2255  double ***hpxij(double ***po, int nhstep
   double **out, cov[NCOVMAX+1];    double **out, cov[NCOVMAX+1];
   double **newm;    double **newm;
   double agexact;    double agexact;
     double agebegin, ageend;
   
   /* Hstepm could be zero and should return the unit matrix */    /* Hstepm could be zero and should return the unit matrix */
   for (i=1;i<=nlstate+ndeath;i++)    for (i=1;i<=nlstate+ndeath;i++)
Line 2236  double ***hpxij(double ***po, int nhstep Line 2269  double ***hpxij(double ***po, int nhstep
       newm=savm;        newm=savm;
       /* Covariates have to be included here again */        /* Covariates have to be included here again */
       cov[1]=1.;        cov[1]=1.;
       agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM;        agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /* age just before transition */
       cov[2]=agexact;        cov[2]=agexact;
       if(nagesqr==1)        if(nagesqr==1)
         cov[3]= agexact*agexact;          cov[3]= agexact*agexact;
Line 2607  double funcone( double *x) Line 2640  double funcone( double *x)
   int s1, s2;    int s1, s2;
   double bbh, survp;    double bbh, survp;
   double agexact;    double agexact;
     double agebegin, ageend;
   /*extern weight */    /*extern weight */
   /* We are differentiating ll according to initial status */    /* We are differentiating ll according to initial status */
   /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/    /*  for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
Line 2625  double funcone( double *x) Line 2659  double funcone( double *x)
           oldm[ii][j]=(ii==j ? 1.0 : 0.0);            oldm[ii][j]=(ii==j ? 1.0 : 0.0);
           savm[ii][j]=(ii==j ? 1.0 : 0.0);            savm[ii][j]=(ii==j ? 1.0 : 0.0);
         }          }
       for(d=0; d<dh[mi][i]; d++){        
         agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
         ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
         for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */
           /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
             and mw[mi+1][i]. dh depends on stepm.*/
         newm=savm;          newm=savm;
         agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;          agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
         cov[2]=agexact;          cov[2]=agexact;
Line 2673  double funcone( double *x) Line 2712  double funcone( double *x)
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */        /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       if(globpr){        if(globpr){
         fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\          fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \   %11.6f %11.6f %11.6f ", \
                 num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,                  num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
                 2*weight[i]*lli,out[s1][s2],savm[s1][s2]);                  2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
         for(k=1,llt=0.,l=0.; k<=nlstate; k++){          for(k=1,llt=0.,l=0.; k<=nlstate; k++){
           llt +=ll[k]*gipmx/gsw;            llt +=ll[k]*gipmx/gsw;
Line 2713  void likelione(FILE *ficres,double p[], Line 2752  void likelione(FILE *ficres,double p[],
       printf("Problem with resultfile: %s\n", fileresilk);        printf("Problem with resultfile: %s\n", fileresilk);
       fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);        fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
     }      }
     fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");      fprintf(ficresilk, "#individual(line's_record) count ageb ageend s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
     fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");      fprintf(ficresilk, "#num_i ageb agend i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");
     /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */      /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
     for(k=1; k<=nlstate; k++)       for(k=1; k<=nlstate; k++) 
       fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);        fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
Line 3188  void pstamp(FILE *fichier) Line 3227  void pstamp(FILE *fichier)
 }  }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[])  void  freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \
                     int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[],\
                     int firstpass,  int lastpass, int stepm, int weightopt, char model[])
 {  /* Some frequencies */  {  /* Some frequencies */
       
   int i, m, jk, j1, bool, z1,j;    int i, m, jk, j1, bool, z1,j;
     int mi; /* Effective wave */
   int first;    int first;
   double ***freq; /* Frequencies */    double ***freq; /* Frequencies */
   double *pp, **prop;    double *pp, **prop;
   double pos,posprop, k2, dateintsum=0,k2cpt=0;    double pos,posprop, k2, dateintsum=0,k2cpt=0;
   char fileresp[FILENAMELENGTH];    char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH];
       double agebegin, ageend;
       
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
   prop=matrix(1,nlstate,iagemin,iagemax+3);    prop=matrix(1,nlstate,iagemin,iagemax+3);
   strcpy(fileresp,"P_");    strcpy(fileresp,"P_");
   strcat(fileresp,fileresu);    strcat(fileresp,fileresu);
     /*strcat(fileresphtm,fileresu);*/
   if((ficresp=fopen(fileresp,"w"))==NULL) {    if((ficresp=fopen(fileresp,"w"))==NULL) {
     printf("Problem with prevalence resultfile: %s\n", fileresp);      printf("Problem with prevalence resultfile: %s\n", fileresp);
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);      fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
     exit(0);      exit(0);
   }    }
   
     strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
     if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
       printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
       fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
       fflush(ficlog);
       exit(70); 
     }
     else{
       fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
   <hr size=\"2\" color=\"#EC5E5E\"> \n\
   Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
             fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
     }
       fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
       
     strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
     if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
       printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
       fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
       fflush(ficlog);
       exit(70); 
     }
     else{
       fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
   <hr size=\"2\" color=\"#EC5E5E\"> \n\
   Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
             fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
     }
     fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
   
   freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3);    freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3);
   j1=0;    j1=0;
       
Line 3215  void  freqsummary(char fileres[], int ia Line 3290  void  freqsummary(char fileres[], int ia
   
   first=1;    first=1;
   
   /* for(k1=1; k1<=j ; k1++){ */  /* Loop on covariates */    for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */
   /*  for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */  
   /*    j1++; */  
   for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){  
       /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
         scanf("%d", i);*/          scanf("%d", i);*/
       for (i=-5; i<=nlstate+ndeath; i++)          for (i=-5; i<=nlstate+ndeath; i++)  
Line 3232  void  freqsummary(char fileres[], int ia Line 3304  void  freqsummary(char fileres[], int ia
               
       dateintsum=0;        dateintsum=0;
       k2cpt=0;        k2cpt=0;
       for (i=1; i<=imx; i++) {        for (i=1; i<=imx; i++) { /* For each individual i */
         bool=1;          bool=1;
         if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */          if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
           for (z1=1; z1<=cptcoveff; z1++)                   for (z1=1; z1<=cptcoveff; z1++)       
Line 3245  void  freqsummary(char fileres[], int ia Line 3317  void  freqsummary(char fileres[], int ia
               /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/                /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
             }               } 
         } /* cptcovn > 0 */          } /* cptcovn > 0 */
    
         if (bool==1){          if (bool==1){
           for(m=firstpass; m<=lastpass; m++){            /* for(m=firstpass; m<=lastpass; m++){ */
             k2=anint[m][i]+(mint[m][i]/12.);            for(mi=1; mi<wav[i];mi++){
             /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/              m=mw[mi][i];
               if(agev[m][i]==0) agev[m][i]=iagemax+1;              /* dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective (mi) waves m=mw[mi][i]
               if(agev[m][i]==1) agev[m][i]=iagemax+2;                 and mw[mi+1][i]. dh depends on stepm. */
               if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i];              agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
               ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /* Age at end of wave and transition */
               if(m >=firstpass && m <=lastpass){
                 k2=anint[m][i]+(mint[m][i]/12.);
                 /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
                 if(agev[m][i]==0) agev[m][i]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */
                 if(agev[m][i]==1) agev[m][i]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */
                 if (s[m][i]>0 && s[m][i]<=nlstate)  /* If status at wave m is known and a live state */
                   prop[s[m][i]][(int)agev[m][i]] += weight[i];  /* At age of beginning of transition, where status is known */
               if (m<lastpass) {                if (m<lastpass) {
                 freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];                  /* if(s[m][i]==4 && s[m+1][i]==4) */
                 freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i];                  /*   printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i]); */
                   if(s[m][i]==-1)
                     printf(" num=%ld m=%d, i=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[i], m, i,s[m][i],s[m+1][i], (int)agev[m][i],agebegin, ageend, (int)((agebegin+ageend)/2.));
                   freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i]; /* At age of beginning of transition, where status is known */
                   /* freq[s[m][i]][s[m+1][i]][(int)((agebegin+ageend)/2.)] += weight[i]; */
                   freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
               }                }
                             }  
               if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) {              if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3)) && (anint[m][i]!=9999) && (mint[m][i]!=99)) {
                 dateintsum=dateintsum+k2;                dateintsum=dateintsum+k2;
                 k2cpt++;                k2cpt++;
                 /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */                /* printf("i=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",i, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
               }              }
               /*}*/              /*}*/
           } /* end m */            } /* end m */
         } /* end bool */          } /* end bool */
       } /* end i = 1 to imx */        } /* end i = 1 to imx */
Line 3272  void  freqsummary(char fileres[], int ia Line 3357  void  freqsummary(char fileres[], int ia
       pstamp(ficresp);        pstamp(ficresp);
       if  (cptcovn>0) {        if  (cptcovn>0) {
         fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, "\n#********** Variable "); 
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);          fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
         fprintf(ficresp, "**********\n#");          fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
           for (z1=1; z1<=cptcoveff; z1++){
             fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
             fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
             fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
           }
             fprintf(ficresp, "**********\n#");
           fprintf(ficresphtm, "**********</h3>\n");
           fprintf(ficresphtmfr, "**********</h3>\n");
         fprintf(ficlog, "\n#********** Variable ");           fprintf(ficlog, "\n#********** Variable "); 
         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
         fprintf(ficlog, "**********\n#");          fprintf(ficlog, "**********\n");
       }        }
       for(i=1; i<=nlstate;i++)         fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
         for(i=1; i<=nlstate;i++) {
         fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);          fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
           fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
         }
       fprintf(ficresp, "\n");        fprintf(ficresp, "\n");
         fprintf(ficresphtm, "\n");
               
         /* Header of frequency table by age */
         fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
         fprintf(ficresphtmfr,"<th>Age</th> ");
         for(jk=-1; jk <=nlstate+ndeath; jk++){
           for(m=-1; m <=nlstate+ndeath; m++){
             if(jk!=0 && m!=0)
               fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
           }
         }
         fprintf(ficresphtmfr, "\n");
         
         /* For each age */
       for(i=iagemin; i <= iagemax+3; i++){        for(i=iagemin; i <= iagemax+3; i++){
         if(i==iagemax+3){          fprintf(ficresphtm,"<tr>");
           if(i==iagemax+1){
             fprintf(ficlog,"1");
             fprintf(ficresphtmfr,"<tr><th>0</th> ");
           }else if(i==iagemax+2){
             fprintf(ficlog,"0");
             fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
           }else if(i==iagemax+3){
           fprintf(ficlog,"Total");            fprintf(ficlog,"Total");
             fprintf(ficresphtmfr,"<tr><th>Total</th> ");
         }else{          }else{
           if(first==1){            if(first==1){
             first=0;              first=0;
             printf("See log file for details...\n");              printf("See log file for details...\n");
           }            }
             fprintf(ficresphtmfr,"<tr><th>%d</th> ",i);
           fprintf(ficlog,"Age %d", i);            fprintf(ficlog,"Age %d", i);
         }          }
         for(jk=1; jk <=nlstate ; jk++){          for(jk=1; jk <=nlstate ; jk++){
Line 3332  void  freqsummary(char fileres[], int ia Line 3450  void  freqsummary(char fileres[], int ia
           if( i <= iagemax){            if( i <= iagemax){
             if(pos>=1.e-5){              if(pos>=1.e-5){
               fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop);                fprintf(ficresp," %d %.5f %.0f %.0f",i,prop[jk][i]/posprop, prop[jk][i],posprop);
                 fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",i,prop[jk][i]/posprop, prop[jk][i],posprop);
               /*probs[i][jk][j1]= pp[jk]/pos;*/                /*probs[i][jk][j1]= pp[jk]/pos;*/
               /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/                /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/
             }              }
             else              else{
               fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop);                fprintf(ficresp," %d NaNq %.0f %.0f",i,prop[jk][i],posprop);
                 fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",i, prop[jk][i],posprop);
               }
           }            }
         }          }
                   
         for(jk=-1; jk <=nlstate+ndeath; jk++)          for(jk=-1; jk <=nlstate+ndeath; jk++){
           for(m=-1; m <=nlstate+ndeath; m++)            for(m=-1; m <=nlstate+ndeath; m++){
             if(freq[jk][m][i] !=0 ) {              if(freq[jk][m][i] !=0 ) { /* minimizing output */
             if(first==1)                if(first==1){
               printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);                  printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
                 }
               fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);                fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);
             }              }
         if(i <= iagemax)              if(jk!=0 && m!=0)
                 fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][i]);
             }
           }
           fprintf(ficresphtmfr,"</tr>\n ");
           if(i <= iagemax){
           fprintf(ficresp,"\n");            fprintf(ficresp,"\n");
             fprintf(ficresphtm,"</tr>\n");
           }
         if(first==1)          if(first==1)
           printf("Others in log...\n");            printf("Others in log...\n");
         fprintf(ficlog,"\n");          fprintf(ficlog,"\n");
       } /* end loop i */        } /* end loop i */
         fprintf(ficresphtm,"</table>\n");
         fprintf(ficresphtmfr,"</table>\n");
       /*}*/        /*}*/
   } /* end j1 */    } /* end j1 */
   dateintmean=dateintsum/k2cpt;     dateintmean=dateintsum/k2cpt; 
     
   fclose(ficresp);    fclose(ficresp);
     fclose(ficresphtm);
     fclose(ficresphtmfr);
   free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);    free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
   free_matrix(prop,1,nlstate,iagemin, iagemax+3);    free_matrix(prop,1,nlstate,iagemin, iagemax+3);
Line 3373  void prevalence(double ***probs, double Line 3506  void prevalence(double ***probs, double
   */    */
     
   int i, m, jk, j1, bool, z1,j;    int i, m, jk, j1, bool, z1,j;
     int mi; /* Effective wave */
     int iage;
     double agebegin, ageend;
   
   double **prop;    double **prop;
   double posprop;     double posprop; 
Line 3392  void prevalence(double ***probs, double Line 3528  void prevalence(double ***probs, double
       
   first=1;    first=1;
   for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){    for(j1=1; j1<= (int) pow(2,cptcoveff);j1++){
     /*for(i1=1; i1<=ncodemax[k1];i1++){      for (i=1; i<=nlstate; i++)  
       j1++;*/        for(iage=iagemin; iage <= iagemax+3; iage++)
                 prop[i][iage]=0.0;
       for (i=1; i<=nlstate; i++)        
         for(m=iagemin; m <= iagemax+3; m++)      for (i=1; i<=imx; i++) { /* Each individual */
           prop[i][m]=0.0;        bool=1;
              if  (cptcovn>0) {  /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
       for (i=1; i<=imx; i++) { /* Each individual */          for (z1=1; z1<=cptcoveff; z1++) 
         bool=1;            if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) 
         if  (cptcovn>0) {              bool=0;
           for (z1=1; z1<=cptcoveff; z1++)         } 
             if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)])         if (bool==1) { 
               bool=0;          /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */
         }           for(mi=1; mi<wav[i];mi++){
         if (bool==1) {             m=mw[mi][i];
           for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/            agebegin=agev[m][i]; /* Age at beginning of wave before transition*/
             /* ageend=agev[m][i]+(dh[m][i])*stepm/YEARM; /\* Age at end of wave and transition *\/ */
             if(m >=firstpass && m <=lastpass){
             y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */              y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
             if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */              if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
               if(agev[m][i]==0) agev[m][i]=iagemax+1;                if(agev[m][i]==0) agev[m][i]=iagemax+1;
               if(agev[m][i]==1) agev[m][i]=iagemax+2;                if(agev[m][i]==1) agev[m][i]=iagemax+2;
               if((int)agev[m][i] <iagemin || (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m);                 if((int)agev[m][i] <iagemin || (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m); 
               if (s[m][i]>0 && s[m][i]<=nlstate) {                 if (s[m][i]>0 && s[m][i]<=nlstate) { 
                 /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/                  /*if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);*/
                 prop[s[m][i]][(int)agev[m][i]] += weight[i];                  prop[s[m][i]][(int)agev[m][i]] += weight[i];/* At age of beginning of transition, where status is known */
                 prop[s[m][i]][iagemax+3] += weight[i];                   prop[s[m][i]][iagemax+3] += weight[i]; 
               }                 } /* end valid statuses */ 
             }              } /* end selection of dates */
           } /* end selection of waves */            } /* end selection of waves */
         }          } /* end effective waves */
       }        } /* end bool */
       for(i=iagemin; i <= iagemax+3; i++){        }
         for(jk=1,posprop=0; jk <=nlstate ; jk++) {       for(i=iagemin; i <= iagemax+3; i++){  
           posprop += prop[jk][i];         for(jk=1,posprop=0; jk <=nlstate ; jk++) { 
         }           posprop += prop[jk][i]; 
                 } 
         for(jk=1; jk <=nlstate ; jk++){             
           if( i <=  iagemax){         for(jk=1; jk <=nlstate ; jk++){       
             if(posprop>=1.e-5){           if( i <=  iagemax){ 
               probs[i][jk][j1]= prop[jk][i]/posprop;            if(posprop>=1.e-5){ 
             } else{              probs[i][jk][j1]= prop[jk][i]/posprop;
               if(first==1){            } else{
                 first=0;              if(first==1){
                 printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);                first=0;
               }                printf("Warning Observed prevalence probs[%d][%d][%d]=%lf because of lack of cases\nSee others on log file...\n",jk,i,j1,probs[i][jk][j1]);
             }              }
           }             }
         }/* end jk */           } 
       }/* end i */         }/* end jk */ 
       }/* end i */ 
     /*} *//* end i1 */      /*} *//* end i1 */
   } /* end j1 */    } /* end j1 */
       
Line 3462  void  concatwav(int wav[], int **dh, int Line 3601  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;    int first, firstwo;
   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;
   jmin=100000;    jmin=100000;
   jmax=-1;    jmax=-1;
   jmean=0.;    jmean=0.;
   for(i=1; i<=imx; i++){    for(i=1; i<=imx; i++){  /* For simple cases and if state is death */
     mi=0;      mi=0;
     m=firstpass;      m=firstpass;
     while(s[m][i] <= nlstate){      while(s[m][i] <= nlstate){  /* a live state */
       if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5)        if(s[m][i]>=1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5)
         mw[++mi][i]=m;          mw[++mi][i]=m;
       if(m >=lastpass)        if(m >=lastpass)
Line 3480  void  concatwav(int wav[], int **dh, int Line 3620  void  concatwav(int wav[], int **dh, int
       else        else
         m++;          m++;
     }/* end while */      }/* end while */
     if (s[m][i] > nlstate){      if (s[m][i] > nlstate){  /* In a death state */
       mi++;     /* Death is another wave */        mi++;     /* Death is another wave */
       /* if(mi==0)  never been interviewed correctly before death */        /* if(mi==0)  never been interviewed correctly before death */
          /* Only death is a correct wave */           /* Only death is a correct wave */
       mw[mi][i]=m;        mw[mi][i]=m;
       }else if (andc[i] != 9999) { /* A death occured after lastpass */
         m++;
         mi++;
         s[m][i]=nlstate+1;  /* We are setting the status to the last of non live state */
         mw[mi][i]=m;
         nbwarn++;
         if(firstwo==0){
           printf("Warning! Death for individual %ld line=%d  occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\nOthers in log file only\n",num[i],i,lastpass,nlstate+1, m);
           fprintf(ficlog,"Warning! Death for individual %ld line=%d  occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m);
           firstwo=1;
         }
         if(firstwo==1){
           fprintf(ficlog,"Warning! Death for individual %ld line=%d  occurred after last wave %d. Since 0.98r4 we considered a status %d at wave %d\n",num[i],i,lastpass,nlstate+1, m);
         }
     }      }
   
     wav[i]=mi;      wav[i]=mi;
     if(mi==0){      if(mi==0){
       nbwarn++;        nbwarn++;
Line 3499  void  concatwav(int wav[], int **dh, int Line 3652  void  concatwav(int wav[], int **dh, int
       }        }
     } /* end mi==0 */      } /* end mi==0 */
   } /* End individuals */    } /* End individuals */
     /* wav and mw are no more changed */
   
     
   for(i=1; i<=imx; i++){    for(i=1; i<=imx; i++){
     for(mi=1; mi<wav[i];mi++){      for(mi=1; mi<wav[i];mi++){
       if (stepm <=0)        if (stepm <=0)
Line 4815  void printinghtml(char fileresu[], char Line 4970  void printinghtml(char fileresu[], char
                   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 estepm ,          \
                   double jprev1, double mprev1,double anprev1, \                    double jprev1, double mprev1,double anprev1, double dateprev1, \
                   double jprev2, double mprev2,double anprev2){                    double jprev2, double mprev2,double anprev2, double dateprev2){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
   
    fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \     fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
    <li><a href='#secondorder'>Result files (second order (variance)</a>\n \     <li><a href='#secondorder'>Result files (second order (variance)</a>\n \
 </ul>");  </ul>");
    fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \     fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n");
  - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",     fprintf(fichtm,"<li>- Observed frequency between two states (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file)<br/>\n",
            jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));             jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTMFR_",".htm"),subdirfext3(optionfilefiname,"PHTMFR_",".htm"));
      fprintf(fichtm,"<li> - Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> (html file) ",
              jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirfext3(optionfilefiname,"PHTM_",".htm"),subdirfext3(optionfilefiname,"PHTM_",".htm"));
      fprintf(fichtm,",  <a href=\"%s\">%s</a> (text file) <br>\n",subdirf2(fileresu,"P_"),subdirf2(fileresu,"P_"));
    fprintf(fichtm,"\     fprintf(fichtm,"\
  - 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_"));
Line 4889  divided by h: <sub>h</sub>P<sub>ij</sub> Line 5047  divided by h: <sub>h</sub>P<sub>ij</sub>
     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++){
         fprintf(fichtm,"<br>\n- Projection of prevalece up 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- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up 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,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);  <img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);
       }        }
     }      }
   
Line 5009  true period expectancies (those weighted Line 5167  true period expectancies (those weighted
 /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)"  */  /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)"  */
     /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */      /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
     fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));      fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
     fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$12):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));      fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$13):6 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
     fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));      fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
     fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$12):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));      fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$13):5 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
     for (i=1; i<= nlstate ; i ++) {      for (i=1; i<= nlstate ; i ++) {
       fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);        fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
       fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot  \"%s\"",subdirf(fileresilk));        fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot  \"%s\"",subdirf(fileresilk));
       fprintf(ficgp,"  u  2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);        fprintf(ficgp,"  u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
       for (j=2; j<= nlstate+ndeath ; j ++) {        for (j=2; j<= nlstate+ndeath ; j ++) {
         fprintf(ficgp,",\\\n \"\" u  2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);          fprintf(ficgp,",\\\n \"\" u  2:($5 == %d && $6==%d ? $10 : 1/0):($12/4.):6 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
       }        }
       fprintf(ficgp,";\nset out; unset ylabel;\n");         fprintf(ficgp,";\nset out; unset ylabel;\n"); 
     }      }
Line 5539  void prevforecast(char fileres[], double Line 5697  void prevforecast(char fileres[], double
      in each health status at the date of interview (if between dateprev1 and dateprev2).       in each health status at the date of interview (if between dateprev1 and dateprev2).
      We still use firstpass and lastpass as another selection.       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);    prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
     
   strcpy(fileresf,"F_");     strcpy(fileresf,"F_"); 
Line 5547  void prevforecast(char fileres[], double Line 5707  void prevforecast(char fileres[], double
     printf("Problem with forecast resultfile: %s\n", fileresf);      printf("Problem with forecast resultfile: %s\n", fileresf);
     fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);      fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
   }    }
   printf("Computing forecasting: result on file '%s' \n", fileresf);    printf("Computing forecasting: result on file '%s', please wait... \n", fileresf);
   fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf);    fprintf(ficlog,"Computing forecasting: result on file '%s', please wait... \n", fileresf);
   
   if (cptcoveff==0) ncodemax[cptcoveff]=1;    if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
Line 5643  void prevforecast(char fileres[], double Line 5803  void prevforecast(char fileres[], double
   if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   fclose(ficresf);    fclose(ficresf);
     printf("End of Computing forecasting \n");
     fprintf(ficlog,"End of Computing forecasting\n");
   
 }  }
   
 /************** Forecasting *****not tested NB*************/  /************** Forecasting *****not tested NB*************/
Line 6482  int calandcheckages(int imx, int maxwav, Line 6645  int calandcheckages(int imx, int maxwav,
   for (i=1; i<=imx; i++)  {    for (i=1; i<=imx; i++)  {
     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);      agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
     for(m=firstpass; (m<= lastpass); m++){      for(m=firstpass; (m<= lastpass); m++){
       if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){        if(s[m][i] >0  || s[m][i]==-1 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ /* What if s[m][i]=-1 */
         if (s[m][i] >= nlstate+1) {          if (s[m][i] >= nlstate+1) {
           if(agedc[i]>0){            if(agedc[i]>0){
             if((int)moisdc[i]!=99 && (int)andc[i]!=9999){              if((int)moisdc[i]!=99 && (int)andc[i]!=9999){
               agev[m][i]=agedc[i];                agev[m][i]=agedc[i];
           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/                /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
             }else {              }else {
               if ((int)andc[i]!=9999){                if ((int)andc[i]!=9999){
                 nbwarn++;                  nbwarn++;
Line 6497  int calandcheckages(int imx, int maxwav, Line 6660  int calandcheckages(int imx, int maxwav,
               }                }
             }              }
           } /* agedc > 0 */            } /* agedc > 0 */
         }          } /* end if */
         else if(s[m][i] !=9){ /* Standard case, age in fractional          else if(s[m][i] !=9){ /* Standard case, age in fractional
                                  years but with the precision of a month */                                   years but with the precision of a month */
           agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);            agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
Line 6513  int calandcheckages(int imx, int maxwav, Line 6676  int calandcheckages(int imx, int maxwav,
           }            }
           /*agev[m][i]=anint[m][i]-annais[i];*/            /*agev[m][i]=anint[m][i]-annais[i];*/
           /*     agev[m][i] = age[i]+2*m;*/            /*     agev[m][i] = age[i]+2*m;*/
         }          } /* en if 9*/
         else { /* =9 */          else { /* =9 */
             /* printf("Debug num[%d]=%ld s[%d][%d]=%d\n",i,num[i], m,i, s[m][i]); */
           agev[m][i]=1;            agev[m][i]=1;
           s[m][i]=-1;            s[m][i]=-1;
         }          }
       }        }
       else /*= 0 Unknown */        else if(s[m][i]==0) /*= 0 Unknown */
         agev[m][i]=1;          agev[m][i]=1;
     }        else{
               printf("Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]); 
           fprintf(ficlog, "Warning, num[%d]=%ld, s[%d][%d]=%d\n", i, num[i], m, i,s[m][i]); 
           agev[m][i]=0;
         }
       } /* End for lastpass */
   }    }
       
   for (i=1; i<=imx; i++)  {    for (i=1; i<=imx; i++)  {
     for(m=firstpass; (m<=lastpass); m++){      for(m=firstpass; (m<=lastpass); m++){
       if (s[m][i] > (nlstate+ndeath)) {        if (s[m][i] > (nlstate+ndeath)) {
Line 7522  Please run with mle=-1 to get a correct Line 7691  Please run with mle=-1 to get a correct
   free_vector(annais,1,n);    free_vector(annais,1,n);
   /* free_matrix(mint,1,maxwav,1,n);    /* free_matrix(mint,1,maxwav,1,n);
      free_matrix(anint,1,maxwav,1,n);*/       free_matrix(anint,1,maxwav,1,n);*/
   free_vector(moisdc,1,n);    /* free_vector(moisdc,1,n); */
   free_vector(andc,1,n);    /* free_vector(andc,1,n); */
   /* */    /* */
       
   wav=ivector(1,imx);    wav=ivector(1,imx);
   dh=imatrix(1,lastpass-firstpass+1,1,imx);    /* dh=imatrix(1,lastpass-firstpass+1,1,imx); */
   bh=imatrix(1,lastpass-firstpass+1,1,imx);    /* bh=imatrix(1,lastpass-firstpass+1,1,imx); */
   mw=imatrix(1,lastpass-firstpass+1,1,imx);    /* mw=imatrix(1,lastpass-firstpass+1,1,imx); */
     dh=imatrix(1,lastpass-firstpass+2,1,imx); /* We are adding a wave if status is unknown at last wave but death occurs after last wave.*/
     bh=imatrix(1,lastpass-firstpass+2,1,imx);
     mw=imatrix(1,lastpass-firstpass+2,1,imx);
         
   /* Concatenates waves */    /* Concatenates waves */
     /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i.
        Death is a valid wave (if date is known).
        mw[mi][i] is the number of (mi=1 to wav[i]) effective wave out of mi of individual i
        dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
        and mw[mi+1][i]. dh depends on stepm.
     */
   
   concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);    concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
   /* */    /* */
     
     free_vector(moisdc,1,n);
     free_vector(andc,1,n);
   
   /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */    /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
   
   nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);     nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); 
Line 7708  Title=%s <br>Datafile=%s Firstpass=%d La Line 7890  Title=%s <br>Datafile=%s Firstpass=%d La
           optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);            optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }    }
   
   fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015</a></font><br>  \    fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-2013-2016-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br>  \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 <font size=\"2\">IMaCh-%s <br> %s</font> \  <font size=\"2\">IMaCh-%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
Line 7738  Title=%s <br>Datafile=%s Firstpass=%d La Line 7920  Title=%s <br>Datafile=%s Firstpass=%d La
       
   /* Calculates basic frequencies. Computes observed prevalence at single age    /* Calculates basic frequencies. Computes observed prevalence at single age
      and prints on file fileres'p'. */       and prints on file fileres'p'. */
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart);    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\
                 firstpass, lastpass,  stepm,  weightopt, model);
   
   fprintf(fichtm,"\n");    fprintf(fichtm,"\n");
   fprintf(fichtm,"<br>Total number of observations=%d <br>\n\    fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
Line 8269  Please run with mle=-1 to get a correct Line 8452  Please run with mle=-1 to get a correct
           
     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,estepm, \
                  jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);                   jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);
               
    /*------------ free_vector  -------------*/     /*------------ free_vector  -------------*/
    /*  chdir(path); */     /*  chdir(path); */
     
     free_ivector(wav,1,imx);      /* free_ivector(wav,1,imx); */  /* Moved after last prevalence call */
     free_imatrix(dh,1,lastpass-firstpass+1,1,imx);      /* free_imatrix(dh,1,lastpass-firstpass+2,1,imx); */
     free_imatrix(bh,1,lastpass-firstpass+1,1,imx);      /* free_imatrix(bh,1,lastpass-firstpass+2,1,imx); */
     free_imatrix(mw,1,lastpass-firstpass+1,1,imx);         /* free_imatrix(mw,1,lastpass-firstpass+2,1,imx);    */
     free_lvector(num,1,n);      free_lvector(num,1,n);
     free_vector(agedc,1,n);      free_vector(agedc,1,n);
     /*free_matrix(covar,0,NCOVMAX,1,n);*/      /*free_matrix(covar,0,NCOVMAX,1,n);*/
Line 8337  Please run with mle=-1 to get a correct Line 8520  Please run with mle=-1 to get a correct
     /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\      /*  printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d,  mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
         ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);          ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
     */      */
       free_ivector(wav,1,imx);
       free_imatrix(dh,1,lastpass-firstpass+2,1,imx);
       free_imatrix(bh,1,lastpass-firstpass+2,1,imx);
       free_imatrix(mw,1,lastpass-firstpass+2,1,imx);   
   
   
     if (mobilav!=0) {      if (mobilav!=0) {
       mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);        mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);

Removed from v.1.212  
changed lines
  Added in v.1.215


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