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