Forecasting

diff --git a/src/imach.c b/src/imach.c
index 60ac82e56643bc5a9c50de72f85ed548e03d64c7..f3f7921a5eb16e4c688a0f21052f5675a90704d5 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -73,6 +73,7 @@ int npar=NPARMAX;
 int nlstate=2; /* Number of live states */\r
 int ndeath=1; /* Number of dead states */\r
 int ncovmodel, ncov;     /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */\r
+int popbased=0, fprev,lprev;\r
 \r
 int *wav; /* Number of waves for this individuual 0 is possible */\r
 int maxwav; /* Maxim number of waves */\r
@@ -1191,7 +1192,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
               bool=0;\r
         }\r
          if (bool==1) {\r
-          for(m=firstpass; m<=lastpass-1; m++){\r
+          for(m=fprev; m<=lprev; m++){\r
             if(agev[m][i]==0) agev[m][i]=agemax+1;\r
             if(agev[m][i]==1) agev[m][i]=agemax+2;\r
             freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];\r
@@ -1215,7 +1216,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
       printf("Age %d", i);\r
     for(jk=1; jk <=nlstate ; jk++){\r
       for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)\r
-       pp[jk] += freq[jk][m][i];\r
+       pp[jk] += freq[jk][m][i]; \r
     }\r
     for(jk=1; jk <=nlstate ; jk++){\r
       for(m=-1, pos=0; m <=0 ; m++)\r
@@ -1225,10 +1226,12 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
       else\r
         printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);\r
     }\r
-    for(jk=1; jk <=nlstate ; jk++){\r
-      for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)\r
+\r
+     for(jk=1; jk <=nlstate ; jk++){\r
+      for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)\r
        pp[jk] += freq[jk][m][i];\r
-    }\r
+     }\r
+\r
     for(jk=1,pos=0; jk <=nlstate ; jk++)\r
       pos += pp[jk];\r
     for(jk=1; jk <=nlstate ; jk++){\r
@@ -1492,6 +1495,12 @@ void varevsij(char fileres[], double ***vareij, double **matcov, double x[], dou
       }\r
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  \r
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);\r
+\r
+      if (popbased==1) {\r
+       for(i=1; i<=nlstate;i++)\r
+         prlim[i][i]=probs[(int)age][i][ij];\r
+      }\r
+      \r
       for(j=1; j<= nlstate; j++){\r
        for(h=0; h<=nhstepm; h++){\r
          for(i=1, gp[h][j]=0.;i<=nlstate;i++)\r
@@ -1503,12 +1512,19 @@ void varevsij(char fileres[], double ***vareij, double **matcov, double x[], dou
        xp[i] = x[i] - (i==theta ?delti[theta]:0);\r
       hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);  \r
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);\r
+\r
+      if (popbased==1) {\r
+       for(i=1; i<=nlstate;i++)\r
+         prlim[i][i]=probs[(int)age][i][ij];\r
+      }\r
+\r
       for(j=1; j<= nlstate; j++){\r
        for(h=0; h<=nhstepm; h++){\r
          for(i=1, gm[h][j]=0.;i<=nlstate;i++)\r
            gm[h][j] += prlim[i][i]*p3mat[i][j][h];\r
        }\r
       }\r
+\r
       for(j=1; j<= nlstate; j++)\r
        for(h=0; h<=nhstepm; h++){\r
          gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];\r
@@ -1773,9 +1789,10 @@ int main()
   int c,  h , cpt,l;\r
   int ju,jl, mi;\r
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;\r
-  int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;\r
- \r
+  int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; \r
+  int mobilav=0, fprevfore=1, lprevfore=1;\r
   int hstepm, nhstepm;\r
+\r
   double bage, fage, age, agelim, agebase;\r
   double ftolpl=FTOL;\r
   double **prlim;\r
@@ -1852,7 +1869,25 @@ split(pathtot, path,optionfile);
   fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncov, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);\r
   printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate,ndeath, maxwav, mle, weightopt,model);\r
   fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncov,nlstate,ndeath,maxwav, mle, weightopt,model);\r
-\r
+while((c=getc(ficpar))=='#' && c!= EOF){\r
+    ungetc(c,ficpar);\r
+    fgets(line, MAXLINE, ficpar);\r
+    puts(line);\r
+    fputs(line,ficparo);\r
+  }\r
+  ungetc(c,ficpar);\r
+  \r
+  fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);\r
+ while((c=getc(ficpar))=='#' && c!= EOF){\r
+    ungetc(c,ficpar);\r
+    fgets(line, MAXLINE, ficpar);\r
+    puts(line);\r
+    fputs(line,ficparo);\r
+  }\r
+  ungetc(c,ficpar);\r
+  \r
+  fscanf(ficpar,"fprevalence=%d lprevalence=%d mob_average=%d\n",&fprevfore,&lprevfore,&mobilav);\r
+ \r
   covar=matrix(0,NCOVMAX,1,n); \r
   cptcovn=0; \r
   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;\r
@@ -2014,8 +2049,8 @@ split(pathtot, path,optionfile);
     if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3;\r
     if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;\r
     if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;\r
-  }\r
-  for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/\r
+    }\r
+    for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/\r
 \r
   /* Calculation of the number of parameter from char model*/\r
   Tvar=ivector(1,15); \r
@@ -2525,7 +2560,9 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
         - Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>\r
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>\r
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>\r
-        - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br><br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);\r
+        - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\r
+        - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\r
+<br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);\r
 \r
  fprintf(fichtm," <li>Graphs</li><p>");\r
 \r
@@ -2663,7 +2700,6 @@ fclose(fichtm);
 \r
   fclose(ficrespij);\r
 \r
-  exit(0);\r
   /*---------- Forecasting ------------------*/\r
 \r
   strcpy(fileresf,"f"); \r
@@ -2682,26 +2718,25 @@ fclose(fichtm);
 \r
   if (cptcoveff==0) ncodemax[cptcoveff]=1;\r
 \r
-  mobaverage= ma3x(1,130 ,1,8, 1,8);\r
-  for (agedeb=bage+3; agedeb<=fage-2; agedeb++)\r
-    for (i=1; i<=nlstate;i++)\r
-      for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)\r
-       mobaverage[(int)agedeb][i][cptcod]=0.;\r
-  \r
-  for (agedeb=bage+4; agedeb<=fage; agedeb++){\r
-    for (i=1; i<=nlstate;i++){\r
-      for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
-       for (cpt=0;cpt<=4;cpt++){\r
-         mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];\r
+  if (mobilav==1) {\r
+    mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
+    for (agedeb=bage+3; agedeb<=fage-2; agedeb++)\r
+      for (i=1; i<=nlstate;i++)\r
+       for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)\r
+         mobaverage[(int)agedeb][i][cptcod]=0.;\r
+    \r
+    for (agedeb=bage+4; agedeb<=fage; agedeb++){\r
+      for (i=1; i<=nlstate;i++){\r
+       for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
+         for (cpt=0;cpt<=4;cpt++){\r
+           mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];\r
          }\r
          mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;\r
+       }\r
       }\r
-    }\r
+    }   \r
   }\r
 \r
-/* if (cptcod==2) printf("m=%f p=%f %d age=%d ",mobaverage[(int)agedeb-2][i][cptcod],probs[(int)agedeb-cpt][i][cptcod],cpt,(int)agedeb-2);*/\r
-\r
-\r
   stepsize=(int) (stepm+YEARM-1)/YEARM;\r
   if (stepm<=24) stepsize=2;\r
 \r
@@ -2725,12 +2760,18 @@ fclose(fichtm);
 \r
       for (agedeb=fage; agedeb>=bage; agedeb--){ \r
        fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);\r
+       if (mobilav==1) {\r
        for(j=1; j<=nlstate;j++) \r
-         fprintf(ficresf,"%.3f ",mobaverage[(int)agedeb][j][cptcod]);\r
-      }\r
+         fprintf(ficresf,"%.5f ",mobaverage[(int)agedeb][j][cptcod]);\r
+       }\r
+       else {\r
+         for(j=1; j<=nlstate;j++) \r
+         fprintf(ficresf,"%.5f ",probs[(int)agedeb][j][cptcod]);\r
+       }\r
+       \r
       for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");\r
-\r
-      for (cpt=1; cpt<=8;cpt++)   \r
+      }\r
+      for (cpt=1; cpt<=NCOVMAX;cpt++)   \r
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */\r
        \r
        nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ \r
@@ -2749,20 +2790,19 @@ fclose(fichtm);
          for(j=1; j<=nlstate+ndeath;j++) {\r
            kk1=0.;\r
            for(i=1; i<=nlstate;i++) {        \r
-             /*   kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];*/\r
-               kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];\r
-           }\r
-           \r
-           if (h*hstepm/YEARM*stepm==cpt) \r
-             fprintf(ficresf," %.5f ", kk1);\r
-         }\r
+             if (mobilav==1)\r
+             kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];\r
+             else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];\r
+           }    \r
+           if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1);\r
          }\r
-       free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
        }\r
+       free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
       }\r
     }\r
+  }\r
+  if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);\r
   fclose(ficresf);\r
-\r
   /*---------- Health expectancies and variances ------------*/\r
 \r
   strcpy(filerest,"t");\r
@@ -2822,6 +2862,11 @@ fclose(fichtm);
       epj=vector(1,nlstate+1);\r
       for(age=bage; age <=fage ;age++){\r
        prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);\r
+       if (popbased==1) {\r
+         for(i=1; i<=nlstate;i++)\r
+           prlim[i][i]=probs[(int)age][i][k];\r
+       }\r
+       \r
        fprintf(ficrest," %.0f",age);\r
        for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){\r
          for(i=1, epj[j]=0.;i <=nlstate;i++) {\r