Variance of one-step transitions added and forecasting of prevalences

diff --git a/src/imach.c b/src/imach.c
index c32eb50e85f9ebd287c3dc1ede2a68deebaa3eed..60ac82e56643bc5a9c50de72f85ed548e03d64c7 100644 (file)
--- a/src/imach.c
+++ b/src/imach.c
@@ -83,8 +83,8 @@ int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */
 double jmean; /* Mean space between 2 waves */\r
 double **oldm, **newm, **savm; /* Working pointers to matrices */\r
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */\r
-FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest;\r
-FILE *ficgp, *fichtm;\r
+FILE *fic,*ficpar, *ficparo,*ficres,  *ficrespl, *ficrespij, *ficrest,*ficresf;\r
+FILE *ficgp, *fichtm,*ficresprob;\r
 FILE *ficreseij;\r
   char filerese[FILENAMELENGTH];\r
  FILE  *ficresvij;\r
@@ -127,7 +127,7 @@ int stepm;
 int m,nb;\r
 int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;\r
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;\r
-double **pmmij;\r
+double **pmmij, ***probs, ***mobaverage;\r
 \r
 double *weight;\r
 int **s; /* Status */\r
@@ -760,7 +760,7 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate )
 \r
 double **matprod2(double **out, double **in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double **b)\r
 {\r
-  /* Computes the matric product of in(1,nrh-nrl+1)(1,nch-ncl+1) times\r
+  /* Computes the matrix product of in(1,nrh-nrl+1)(1,nch-ncl+1) times\r
      b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */\r
   /* in, b, out are matrice of pointers which should have been initialized \r
      before: only the contents of out is modified. The function returns\r
@@ -1160,7 +1160,7 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
   char fileresp[FILENAMELENGTH];\r
 \r
   pp=vector(1,nlstate);\r
-\r
+ probs= ma3x(1,130 ,1,8, 1,8);\r
   strcpy(fileresp,"p");\r
   strcat(fileresp,fileres);\r
   if((ficresp=fopen(fileresp,"w"))==NULL) {\r
@@ -1237,8 +1237,11 @@ void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev
       else\r
        printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);\r
       if( i <= (int) agemax){\r
-       if(pos>=1.e-5)\r
+       if(pos>=1.e-5){\r
          fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);\r
+         probs[i][jk][j1]= pp[jk]/pos;\r
+         /*printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);*/\r
+       }\r
       else\r
          fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);\r
       }\r
@@ -1635,7 +1638,110 @@ void varprevlim(char fileres[], double **varpl, double **matcov, double x[], dou
 \r
 }\r
 \r
+/************ Variance of one-step probabilities  ******************/\r
+void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)\r
+{\r
+  int i, j;\r
+  int k=0, cptcode;\r
+  double **dnewm,**doldm;\r
+  double *xp;\r
+  double *gp, *gm;\r
+  double **gradg, **trgradg;\r
+  double age,agelim, cov[NCOVMAX];\r
+  int theta;\r
+  char fileresprob[FILENAMELENGTH];\r
+\r
+  strcpy(fileresprob,"prob"); \r
+  strcat(fileresprob,fileres);\r
+  if((ficresprob=fopen(fileresprob,"w"))==NULL) {\r
+    printf("Problem with resultfile: %s\n", fileresprob);\r
+  }\r
+  printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);\r
+  \r
+\r
+  xp=vector(1,npar);\r
+  dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);\r
+  doldm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,(nlstate+ndeath)*(nlstate+ndeath));\r
+  \r
+  cov[1]=1;\r
+  for (age=bage; age<=fage; age ++){ \r
+    cov[2]=age;\r
+    gradg=matrix(1,npar,1,9);\r
+    trgradg=matrix(1,9,1,npar);\r
+    gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));\r
+    gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));\r
+    \r
+    for(theta=1; theta <=npar; theta++){\r
+      for(i=1; i<=npar; i++)\r
+       xp[i] = x[i] + (i==theta ?delti[theta]:0);\r
+     \r
+      pmij(pmmij,cov,ncovmodel,xp,nlstate);\r
+   \r
+      k=0;\r
+      for(i=1; i<= (nlstate+ndeath); i++){\r
+       for(j=1; j<=(nlstate+ndeath);j++){\r
+          k=k+1;\r
+         gp[k]=pmmij[i][j];\r
+       }\r
+      }\r
+\r
+      for(i=1; i<=npar; i++)\r
+       xp[i] = x[i] - (i==theta ?delti[theta]:0);\r
+    \r
+\r
+      pmij(pmmij,cov,ncovmodel,xp,nlstate);\r
+      k=0;\r
+      for(i=1; i<=(nlstate+ndeath); i++){\r
+       for(j=1; j<=(nlstate+ndeath);j++){\r
+         k=k+1;\r
+         gm[k]=pmmij[i][j];\r
+       }\r
+      }\r
+     \r
+       for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++) \r
+          gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];  \r
+    }\r
+\r
+     for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)\r
+      for(theta=1; theta <=npar; theta++)\r
+      trgradg[j][theta]=gradg[theta][j];\r
+ \r
+     matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);\r
+     matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);\r
+\r
+     pmij(pmmij,cov,ncovmodel,x,nlstate);\r
+\r
+     k=0;\r
+     for(i=1; i<=(nlstate+ndeath); i++){\r
+       for(j=1; j<=(nlstate+ndeath);j++){\r
+        k=k+1;\r
+        gm[k]=pmmij[i][j];\r
+       }\r
+     }\r
+     \r
+     /*printf("\n%d ",(int)age);\r
+     for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){\r
+       \r
+\r
+       printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i]));\r
+     }*/\r
+\r
+  fprintf(ficresprob,"\n%d ",(int)age);\r
+\r
+  for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){\r
+    if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);\r
+if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);\r
+  }\r
 \r
+    free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));\r
+    free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));\r
+    free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);\r
+    free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);\r
+}\r
+ free_vector(xp,1,npar);\r
+fclose(ficresprob);\r
+ exit(0);\r
+}\r
 \r
 /***********************************************/\r
 /**************** Main Program *****************/\r
@@ -1658,7 +1764,7 @@ int main()
   char line[MAXLINE], linepar[MAXLINE];\r
   char title[MAXLINE];\r
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];\r
-  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];\r
+  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];\r
   char filerest[FILENAMELENGTH];\r
   char fileregp[FILENAMELENGTH];\r
   char path[80],pathc[80],pathcd[80],pathtot[80],model[20];\r
@@ -1682,9 +1788,12 @@ int main()
   double ***eij, ***vareij;\r
   double **varpl; /* Variances of prevalence limits by age */\r
   double *epj, vepp;\r
+  double kk1;\r
+\r
   char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";\r
   char *alph[]={"a","a","b","c","d","e"}, str[4];\r
 \r
+\r
   char z[1]="c", occ;\r
 #include <sys/time.h>\r
 #include <time.h>\r
@@ -2300,7 +2409,7 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),c
       fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);   \r
       fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);\r
     } \r
-  }\r
+  }  \r
 \r
   /* proba elementaires */\r
    for(i=1,jk=1; i <=nlstate; i++){\r
@@ -2501,6 +2610,7 @@ fclose(fichtm);
       }\r
     }\r
   fclose(ficrespl);\r
+\r
   /*------------- h Pij x at various ages ------------*/\r
   \r
   strcpy(filerespij,"pij");  strcat(filerespij,fileres);\r
@@ -2510,7 +2620,7 @@ fclose(fichtm);
   printf("Computing pij: result on file '%s' \n", filerespij);\r
   \r
   stepsize=(int) (stepm+YEARM-1)/YEARM;\r
-  if (stepm<=24) stepsize=2;\r
+  /*if (stepm<=24) stepsize=2;*/\r
 \r
   agelim=AGESUP;\r
   hstepm=stepsize*YEARM; /* Every year of age */\r
@@ -2549,8 +2659,110 @@ fclose(fichtm);
     }\r
   }\r
 \r
+  /* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/\r
+\r
   fclose(ficrespij);\r
 \r
+  exit(0);\r
+  /*---------- Forecasting ------------------*/\r
+\r
+  strcpy(fileresf,"f"); \r
+  strcat(fileresf,fileres);\r
+  if((ficresf=fopen(fileresf,"w"))==NULL) {\r
+    printf("Problem with forecast resultfile: %s\n", fileresf);goto end;\r
+  }\r
+  printf("Computing forecasting: result on file '%s' \n", fileresf);\r
+\r
+  /* Mobile average */\r
+\r
+  /* for (agedeb=bage; agedeb<=fage; agedeb++)\r
+    for (i=1; i<=nlstate;i++)\r
+      for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++)\r
+      printf("%f %d i=%d j1=%d\n", probs[(int)agedeb][i][cptcod],(int) agedeb,i,cptcod);*/\r
+\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
+         }\r
+         mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;\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
+  agelim=AGESUP;\r
+  hstepm=stepsize*YEARM; /* Every year of age */\r
+  hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ \r
+  hstepm=12;\r
+   k=0;\r
+  for(cptcov=1;cptcov<=i1;cptcov++){\r
+    for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){\r
+      k=k+1;\r
+      fprintf(ficresf,"\n#****** ");\r
+      for(j=1;j<=cptcoveff;j++) {\r
+       fprintf(ficresf,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);\r
+      }\r
+      \r
+      fprintf(ficresf,"******\n");\r
+\r
+      fprintf(ficresf,"# StartingAge FinalAge Horizon(in years)");\r
+      for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);\r
+\r
+      for (agedeb=fage; agedeb>=bage; agedeb--){ \r
+       fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);\r
+       for(j=1; j<=nlstate;j++) \r
+         fprintf(ficresf,"%.3f ",mobaverage[(int)agedeb][j][cptcod]);\r
+      }\r
+      for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");\r
+\r
+      for (cpt=1; cpt<=8;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
+       nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */\r
+       /*printf("stepm=%d hstepm=%d nhstepm=%d \n",stepm,hstepm,nhstepm);*/\r
+\r
+       p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+       oldm=oldms;savm=savms;\r
+       hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  \r
+               \r
+       for (h=0; h<=nhstepm; h++){\r
+       \r
+         if (h*hstepm/YEARM*stepm==cpt)\r
+ fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);\r
+         \r
+         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
+         }\r
+       free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);\r
+       }\r
+      }\r
+    }\r
+  fclose(ficresf);\r
+\r
   /*---------- Health expectancies and variances ------------*/\r
 \r
   strcpy(filerest,"t");\r
@@ -2629,6 +2841,9 @@ fclose(fichtm);
     }\r
   }\r
        \r
+       \r
+\r
+\r
  fclose(ficreseij);\r
  fclose(ficresvij);\r
   fclose(ficrest);\r
@@ -2674,7 +2889,7 @@ strcpy(fileresvpl,"vpl");
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);\r
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);\r
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);\r
-  \r
+ \r
   free_matrix(matcov,1,npar,1,npar);\r
   free_vector(delti,1,npar);\r
   \r
@@ -2692,9 +2907,7 @@ strcpy(fileresvpl,"vpl");
 #ifdef windows\r
  chdir(pathcd);\r
 #endif \r
- /*system("wgnuplot graph.plt");*/\r
- /*system("../gp37mgw/wgnuplot graph.plt");*/\r
- /*system("cd ../gp37mgw");*/\r
+ \r
  system("..\\gp37mgw\\wgnuplot graph.plt");\r
 \r
 #ifdef windows\r