Correction of initialisations in case of several death states.

diff --git a/src/imach.c b/src/imach.c
index b3ac7208e029d4ed9a002ab000fa9087c68b78a5..d578d95eec74b06463d0d4142f9301fde3d89088 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -919,7 +919,7 @@ double func( double *x)
   double sw; /* Sum of weights */
   double lli; /* Individual log likelihood */
   int s1, s2;
-  double bbh;
+  double bbh, survp;
   long ipmx;
   /*extern weight */
   /* We are differentiating ll according to initial status */
@@ -972,6 +972,14 @@ double func( double *x)
         * is higher than the multiple of stepm and negative otherwise.
         */
        /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/
+        /* if s2=-2 lli=out[1][1]+out[1][2];*/
+       if (s2==-2) {
+         for (j=1,survp=0. ; j<=nlstate; j++) 
+           survp += out[s1][j];
+         lli= survp;
+       }
+
+       else
        lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));  /* linear interpolation */
        /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
        /*if(lli ==000.0)*/
@@ -1416,7 +1424,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
     fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
     exit(0);
   }
-  freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
+  freq= ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3);
   j1=0;
   
   j=cptcoveff;
@@ -1531,8 +1539,8 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
          }
        }
        
-       for(jk=-1; jk <=nlstate+ndeath; jk++)
-         for(m=-1; m <=nlstate+ndeath; m++)
+       for(jk=-2; jk <=nlstate+ndeath; jk++)
+         for(m=-2; m <=nlstate+ndeath; m++)
            if(freq[jk][m][i] !=0 ) {
            if(first==1)
              printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
@@ -1549,7 +1557,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
   dateintmean=dateintsum/k2cpt; 
  
   fclose(ficresp);
-  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
+  free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);
   
   /* End of Freq */
@@ -1566,7 +1574,7 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
 
   pp=vector(1,nlstate);
   
-  freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
+  freq=ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,agemin,agemax+3);
   j1=0;
   
   j=cptcoveff;
@@ -1576,8 +1584,8 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
     for(i1=1; i1<=ncodemax[k1];i1++){
       j1++;
       
-      for (i=-1; i<=nlstate+ndeath; i++)  
-       for (jk=-1; jk<=nlstate+ndeath; jk++)  
+      for (i=-2; i<=nlstate+ndeath; i++)  
+       for (jk=-2; jk<=nlstate+ndeath; jk++)  
          for(m=agemin; m <= agemax+3; m++)
            freq[i][jk][m]=0;
      
@@ -1634,7 +1642,7 @@ void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, i
   } /* end k1 */
 
   
-  free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
+  free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);
   
 }  /* End of Freq */
@@ -1664,7 +1672,7 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
     mi=0;
     m=firstpass;
     while(s[m][i] <= nlstate){
-      if(s[m][i]>=1)
+      if(s[m][i]>=1 || s[m][i]==-2)
        mw[++mi][i]=m;
       if(m >=lastpass)
        break;
@@ -1704,10 +1712,12 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
          if (j <= jmin) jmin=j;
          sum=sum+j;
          /*if (j<0) printf("j=%d num=%d \n",j,i); */
+         /*      printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
          }
        }
        else{
          j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
+         /*      printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
          k=k+1;
          if (j >= jmax) jmax=j;
          else if (j <= jmin)jmin=j;
@@ -1742,7 +1752,6 @@ void  concatwav(int wav[], int **dh, int **bh,  int **mw, int **s, double *agedc
            bh[mi][i]=ju; /* At least one step */
            printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);
          }
-         if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);
        }
       } /* end if mle */
     } /* end wave */
@@ -2136,7 +2145,7 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
          as a weighted average of prlim.
       */
       for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){
-       for(i=1; i<= nlstate; i++)
+       for(i=1,gpp[j]=0.; i<= nlstate; i++)
          gpp[j] += prlim[i][i]*p3mat[i][j][1];
       }    
       /* end probability of death */
@@ -2167,8 +2176,8 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
          as a weighted average of prlim.
       */
       for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){
-       for(i=1; i<= nlstate; i++)
-         gmp[j] += prlim[i][i]*p3mat[i][j][1];
+       for(i=1,gmp[j]=0.; i<= nlstate; i++)
+         gmp[j] += prlim[i][i]*p3mat[i][j][1];
       }    
       /* end probability of death */
 
@@ -2176,6 +2185,7 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
        for(h=0; h<=nhstepm; h++){
          gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
        }
+
       for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
        gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
       }
@@ -2190,8 +2200,9 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
          trgradg[h][j][theta]=gradg[h][theta][j];
 
     for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
-      for(theta=1; theta <=npar; theta++)
+      for(theta=1; theta <=npar; theta++) {
        trgradgp[j][theta]=gradgp[theta][j];
+      }
 
     hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
     for(i=1;i<=nlstate;i++)
@@ -2211,9 +2222,12 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
     /* pptj */
     matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
     matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
-    for(j=nlstate+1;j<=nlstate+ndeath;j++)
-      for(i=nlstate+1;i<=nlstate+ndeath;i++)
+ 
+   for(j=nlstate+1;j<=nlstate+ndeath;j++)
+     for(i=nlstate+1;i<=nlstate+ndeath;i++){
        varppt[j][i]=doldmp[j][i];
+     }
+
     /* end ppptj */
     /*  x centered again */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);  
@@ -2233,8 +2247,8 @@ void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double
        computed over hstepm (estepm) matrices product = hstepm*stepm months) 
        as a weighted average of prlim.
     */
-    for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){
-      for(i=1; i<= nlstate; i++)
+    for(j=nlstate+1;j<=nlstate+ndeath;j++){
+      for(i=1,gmp[j]=0.;i<= nlstate; i++) 
        gmp[j] += prlim[i][i]*p3mat[i][j][1]; 
     }    
     /* end probability of death */
@@ -3704,11 +3718,10 @@ int main(int argc, char *argv[])
   for (i=1; i<=imx; i++)  {
     agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
     for(m=1; (m<= maxwav); m++){
-      if(s[m][i] >0){
+      if(s[m][i] >0 || s[m][i]==-2){
        if (s[m][i] >= nlstate+1) {
          if(agedc[i]>0)
-           if(moisdc[i]!=99 && andc[i]!=9999)
-             agev[m][i]=agedc[i];
+           if(moisdc[i]!=99 && andc[i]!=9999) agev[m][i]=agedc[i];
          /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
            else {
              if (andc[i]!=9999){
@@ -3770,7 +3783,8 @@ int main(int argc, char *argv[])
   dh=imatrix(1,lastpass-firstpass+1,1,imx);
   bh=imatrix(1,lastpass-firstpass+1,1,imx);
   mw=imatrix(1,lastpass-firstpass+1,1,imx);
-   
+  
+
   /* Concatenates waves */
   concatwav(wav, dh, bh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
 
@@ -3979,6 +3993,7 @@ int main(int argc, char *argv[])
   fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
   fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
 
+
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
 
   /*------------ gnuplot -------------*/
@@ -4141,7 +4156,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
 
 
   /*---------- Forecasting ------------------*/
-  if((stepm == 1) && (strcmp(model,".")==0)){
+   if((stepm == 1) && (strcmp(model,".")==0)){
     prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
     if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);
   } 
@@ -4149,7 +4164,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
     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);
-  }
+    }
   
 
   /*---------- Health expectancies and variances ------------*/
@@ -4173,6 +4188,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
   printf("Computing Health Expectancies: result on file '%s' \n", filerese);
   fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);
 
+
   strcpy(fileresv,"v");
   strcat(fileresv,fileres);
   if((ficresvij=fopen(fileresv,"w"))==NULL) {