[BACK]Return to imach.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.27 and 1.33

version 1.27, 2002/02/28 17:49:07 version 1.33, 2002/03/12 22:13:38
Line 1185  void  freqsummary(char fileres[], int ag Line 1185  void  freqsummary(char fileres[], int ag
   double pos, k2, dateintsum=0,k2cpt=0;    double pos, k2, dateintsum=0,k2cpt=0;
   FILE *ficresp;    FILE *ficresp;
   char fileresp[FILENAMELENGTH];    char fileresp[FILENAMELENGTH];
    
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
   probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);    probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   strcpy(fileresp,"p");    strcpy(fileresp,"p");
Line 1307  void  freqsummary(char fileres[], int ag Line 1307  void  freqsummary(char fileres[], int ag
 }  }
   
 /************ Prevalence ********************/  /************ Prevalence ********************/
 void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedate)  void prevalence(int agemin, float agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2, double calagedate)
 {  /* Some frequencies */  {  /* Some frequencies */
     
   int i, m, jk, k1, i1, j1, bool, z1,z2,j;    int i, m, jk, k1, i1, j1, bool, z1,z2,j;
Line 1332  void prevalence(int agemin, int agemax, Line 1332  void prevalence(int agemin, int agemax,
         for (jk=-1; jk<=nlstate+ndeath; jk++)            for (jk=-1; jk<=nlstate+ndeath; jk++)  
           for(m=agemin; m <= agemax+3; m++)            for(m=agemin; m <= agemax+3; m++)
             freq[i][jk][m]=0;              freq[i][jk][m]=0;
             
       for (i=1; i<=imx; i++) {        for (i=1; i<=imx; i++) {
         bool=1;          bool=1;
         if  (cptcovn>0) {          if  (cptcovn>0) {
           for (z1=1; z1<=cptcoveff; z1++)            for (z1=1; z1<=cptcoveff; z1++)
             if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])              if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
               bool=0;                bool=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.);              k2=anint[m][i]+(mint[m][i]/12.);
Line 1347  void prevalence(int agemin, int agemax, Line 1347  void prevalence(int agemin, int agemax,
               if(agev[m][i]==0) agev[m][i]=agemax+1;                if(agev[m][i]==0) agev[m][i]=agemax+1;
               if(agev[m][i]==1) agev[m][i]=agemax+2;                if(agev[m][i]==1) agev[m][i]=agemax+2;
               freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];                freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
               freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];                  /* freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];  */
             }              }
           }            }
         }          }
       }        }
        
         for(i=(int)agemin; i <= (int)agemax+3; i++){          for(i=(int)agemin; i <= (int)agemax+3; i++){
           for(jk=1; jk <=nlstate ; jk++){            for(jk=1; jk <=nlstate ; jk++){
             for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)              for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
Line 1444  void  concatwav(int wav[], int **dh, int Line 1443  void  concatwav(int wav[], int **dh, int
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           if (j <= jmin) jmin=j;            if (j <= jmin) jmin=j;
           sum=sum+j;            sum=sum+j;
           /* if (j<10) printf("j=%d num=%d ",j,i); */            /*if (j<0) printf("j=%d num=%d \n",j,i); */
           }            }
         }          }
         else{          else{
Line 1452  void  concatwav(int wav[], int **dh, int Line 1451  void  concatwav(int wav[], int **dh, int
           k=k+1;            k=k+1;
           if (j >= jmax) jmax=j;            if (j >= jmax) jmax=j;
           else if (j <= jmin)jmin=j;            else if (j <= jmin)jmin=j;
           /*   if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */            /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
           sum=sum+j;            sum=sum+j;
         }          }
         jk= j/stepm;          jk= j/stepm;
Line 1528  void tricode(int *Tvar, int **nbcode, in Line 1527  void tricode(int *Tvar, int **nbcode, in
 void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij)  void evsij(char fileres[], double ***eij, double x[], int nlstate, int stepm, int bage, int fage, double **oldm, double **savm, int ij)
 {  {
   /* Health expectancies */    /* Health expectancies */
   int i, j, nhstepm, hstepm, h;    int i, j, nhstepm, hstepm, h, nstepm, k;
   double age, agelim,hf;    double age, agelim,hf;
   double ***p3mat;    double ***p3mat;
     
Line 1539  void evsij(char fileres[], double ***eij Line 1538  void evsij(char fileres[], double ***eij
       fprintf(ficreseij," %1d-%1d",i,j);        fprintf(ficreseij," %1d-%1d",i,j);
   fprintf(ficreseij,"\n");    fprintf(ficreseij,"\n");
   
   hstepm=1*YEARM; /*  Every j years of age (in month) */    k=1;             /* For example stepm=6 months */
   hstepm=hstepm/stepm; /* Typically in stepm units, if j= 2 years, = 2/6 months = 4 */    hstepm=k*YEARM; /* (a) Every k years of age (in months), for example every k=2 years 24 m */
     hstepm=stepm;   /* or (b) We decided to compute the life expectancy with the smallest unit */
     /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
        nhstepm is the number of hstepm from age to agelim
        nstepm is the number of stepm from age to agelin.
        Look at hpijx to understand the reason of that which relies in memory size
        and note for a fixed period like k years */
     /* We decided (b) to get a life expectancy respecting the most precise curvature of the
        survival function given by stepm (the optimization length). Unfortunately it
        means that if the survival funtion is printed only each two years of age and if
        you sum them up and add 1 year (area under the trapezoids) you won't get the same
        results. So we changed our mind and took the option of the best precision.
     */
     hstepm=hstepm/stepm; /* Typically in stepm units, if k= 2 years, = 2/6 months = 4 */
   
   agelim=AGESUP;    agelim=AGESUP;
   for (age=bage; age<=fage; age ++){ /* If stepm=6 months */    for (age=bage; age<=fage; age ++){ /* If stepm=6 months */
     /* nhstepm age range expressed in number of stepm */      /* nhstepm age range expressed in number of stepm */
     nhstepm=(int) rint((agelim-age)*YEARM/stepm);      nstepm=(int) rint((agelim-age)*YEARM/stepm);
     /* Typically if 20 years = 20*12/6=40 stepm */      /* Typically if 20 years nstepm = 20*12/6=40 stepm */
     if (stepm >= YEARM) hstepm=1;      /* if (stepm >= YEARM) hstepm=1;*/
     nhstepm = nhstepm/hstepm;/* Expressed in hstepm, typically 40/4=10 */      nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
     /* Computed by stepm unit matrices, product of hstepm matrices, stored      /* Computed by stepm unit matrices, product of hstepm matrices, stored
        in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */         in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
     hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);        hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm, savm, ij);  
       hf=hstepm*stepm/YEARM;  /* Duration of hstepm expressed in year unit. */
   
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++)        for(j=1; j<=nlstate;j++)
         for (h=0, eij[i][j][(int)age]=0; h<=nhstepm; h++){          for (h=0, eij[i][j][(int)age]=0; h<=nhstepm-1; h++){
           eij[i][j][(int)age] +=p3mat[i][j][h];            eij[i][j][(int)age] += (p3mat[i][j][h]+p3mat[i][j][h+1])/2.0*hf;
             /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
         }          }
          fprintf(ficreseij,"%3.0f",age );
     hf=1;  
     if (stepm >= YEARM) hf=stepm/YEARM;  
     fprintf(ficreseij,"%.0f",age );  
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){        for(j=1; j<=nlstate;j++){
         fprintf(ficreseij," %.4f", hf*eij[i][j][(int)age]);          fprintf(ficreseij," %9.4f", eij[i][j][(int)age]);
       }        }
     fprintf(ficreseij,"\n");      fprintf(ficreseij,"\n");
     free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
Line 1883  fclose(ficresprob); Line 1892  fclose(ficresprob);
 }  }
   
 /******************* Printing html file ***********/  /******************* Printing html file ***********/
 void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[] ){  void printinghtml(char fileres[], char title[], char datafile[], int firstpass, int lastpass, int stepm, int weightopt, char model[],int imx,int jmin, int jmax, double jmeanint,char optionfile[],char optionfilehtm[],char rfileres[] ){
   int jj1, k1, i1, cpt;    int jj1, k1, i1, cpt;
   FILE *fichtm;    FILE *fichtm;
   /*char optionfilehtm[FILENAMELENGTH];*/    /*char optionfilehtm[FILENAMELENGTH];*/
Line 1894  void printinghtml(char fileres[], char t Line 1903  void printinghtml(char fileres[], char t
     printf("Problem with %s \n",optionfilehtm), exit(0);      printf("Problem with %s \n",optionfilehtm), exit(0);
   }    }
   
  fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\">   fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.71e </font> <hr size=\"2\" color=\"#EC5E5E\">
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
   
 Total number of observations=%d <br>  Total number of observations=%d <br>
Line 1910  Interval (in months) between two waves: Line 1919  Interval (in months) between two waves:
         - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>          - Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
         - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>          - Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
         - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>          - Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
         - Prevalences and population forecasting: <a href=\"f%s\">f%s</a> <br>          - Prevalences forecasting: <a href=\"f%s\">f%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,fileres,fileres);          - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>
           <br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
     
 fprintf(fichtm," <li>Graphs</li><p>");  fprintf(fichtm," <li>Graphs</li><p>");
   
Line 1953  fclose(fichtm); Line 1963  fclose(fichtm);
 }  }
   
 /******************* Gnuplot file **************/  /******************* Gnuplot file **************/
 void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemax, double fage , char pathc[], double p[]){  void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double agemin, double agemaxpar, double fage , char pathc[], double p[]){
   
   int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;    int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
   
Line 2048  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2058  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
     for (k1=1; k1<= m ; k1 ++) {      for (k1=1; k1<= m ; k1 ++) {
     for (cpt=1; cpt<nlstate ; cpt ++) {      for (cpt=1; cpt<nlstate ; cpt ++) {
       k=3;        k=3;
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",agemin,agemax,fileres,k1,k+cpt+1,k+1);        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter gif small size 400,300\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",agemin,agemaxpar,fileres,k1,k+cpt+1,k+1);
   
       for (i=1; i< nlstate ; i ++)        for (i=1; i< nlstate ; i ++)
         fprintf(ficgp,"+$%d",k+i+1);          fprintf(ficgp,"+$%d",k+i+1);
Line 2080  fprintf(ficgp,"\nset out \"v%s%d%d.gif\" Line 2090  fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
     }      }
   
     for(jk=1; jk <=m; jk++) {      for(jk=1; jk <=m; jk++) {
   fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemax);    fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot  [%.f:%.f] ",agemin,agemaxpar);
    i=1;     i=1;
    for(k2=1; k2<=nlstate; k2++) {     for(k2=1; k2<=nlstate; k2++) {
      k3=i;       k3=i;
Line 2159  prevforecast(char fileres[], double anpr Line 2169  prevforecast(char fileres[], double anpr
   
  agelim=AGESUP;   agelim=AGESUP;
 calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;  calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
    
   prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);  
   
     prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
    
     
   strcpy(fileresf,"f");    strcpy(fileresf,"f");
   strcat(fileresf,fileres);    strcat(fileresf,fileres);
Line 2210  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2220  calagedate=(anproj1+mproj1/12.+jproj1/36
       for (cpt=0; cpt<=(anproj2-anproj1);cpt++) {        for (cpt=0; cpt<=(anproj2-anproj1);cpt++) {
         fprintf(ficresf,"\n");          fprintf(ficresf,"\n");
         fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);            fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);  
          
         for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){          for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
           nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);            nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
           nhstepm = nhstepm/hstepm;            nhstepm = nhstepm/hstepm;
Line 2250  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2260  calagedate=(anproj1+mproj1/12.+jproj1/36
   fclose(ficresf);    fclose(ficresf);
 }  }
 /************** Forecasting ******************/  /************** Forecasting ******************/
 populforecast(char fileres[], double anproj1,double mproj1,double jproj1,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){  populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double agemin, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
     
   int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;    int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
   int *popage;    int *popage;
Line 2259  populforecast(char fileres[], double anp Line 2269  populforecast(char fileres[], double anp
   double ***p3mat,***tabpop,***tabpopprev;    double ***p3mat,***tabpop,***tabpopprev;
   char filerespop[FILENAMELENGTH];    char filerespop[FILENAMELENGTH];
   
 tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
 tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
  agelim=AGESUP;    agelim=AGESUP;
 calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;    calagedate=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;
     
   prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);    prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
    
     
   strcpy(filerespop,"pop");    strcpy(filerespop,"pop");
   strcat(filerespop,fileres);    strcat(filerespop,fileres);
   if((ficrespop=fopen(filerespop,"w"))==NULL) {    if((ficrespop=fopen(filerespop,"w"))==NULL) {
Line 2312  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2322  calagedate=(anproj1+mproj1/12.+jproj1/36
         fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);          fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       }        }
       fprintf(ficrespop,"******\n");        fprintf(ficrespop,"******\n");
       fprintf(ficrespop,"# StartingAge FinalAge");        fprintf(ficrespop,"# Age");
       for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespop," P.%d",j);        for(j=1; j<=nlstate+ndeath;j++) fprintf(ficrespop," P.%d",j);
       if (popforecast==1)  fprintf(ficrespop," [Population]");        if (popforecast==1)  fprintf(ficrespop," [Population]");
             
       for (cpt=0; cpt<=0;cpt++) {        for (cpt=0; cpt<=0;cpt++) {
         fprintf(ficrespop,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);            fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);  
                 
         for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){          for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
           nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);            nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
Line 2329  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2339  calagedate=(anproj1+mproj1/12.+jproj1/36
                 
           for (h=0; h<=nhstepm; h++){            for (h=0; h<=nhstepm; h++){
             if (h==(int) (calagedate+YEARM*cpt)) {              if (h==(int) (calagedate+YEARM*cpt)) {
               fprintf(ficrespop,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);                fprintf(ficrespop,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm);
             }              }
             for(j=1; j<=nlstate+ndeath;j++) {              for(j=1; j<=nlstate+ndeath;j++) {
               kk1=0.;kk2=0;                kk1=0.;kk2=0;
Line 2349  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2359  calagedate=(anproj1+mproj1/12.+jproj1/36
             for(i=1; i<=nlstate;i++){              for(i=1; i<=nlstate;i++){
               kk1=0.;                kk1=0.;
                 for(j=1; j<=nlstate;j++){                  for(j=1; j<=nlstate;j++){
                   kk1= kk1+Tabpop[(int)(agedeb)][j][cptcod];                    kk1= kk1+tabpop[(int)(agedeb)][j][cptcod];
                 }                  }
                   tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedate+12*cpt)*hstepm/YEARM*stepm-1)];                    tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1*popeffectif[(int)(agedeb+(calagedate+12*cpt)*hstepm/YEARM*stepm-1)];
             }              }
   
 if (h==(int)(calagedate+12*cpt)) for(j=1; j<=nlstate;j++) fprintf(ficrespop," %.3f",tabpopprev[(int)(agedeb+1)][j][cptcod]);              if (h==(int)(calagedate+12*cpt)) for(j=1; j<=nlstate;j++)
                 fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
           }            }
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);            free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         }          }
Line 2362  if (h==(int)(calagedate+12*cpt)) for(j=1 Line 2373  if (h==(int)(calagedate+12*cpt)) for(j=1
     
   /******/    /******/
   
       for (cpt=1; cpt<=4;cpt++) {        for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {
         fprintf(ficrespop,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);            fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);  
         for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){          for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(agemin-((int)calagedate %12)/12.); agedeb--){
           nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);            nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
           nhstepm = nhstepm/hstepm;            nhstepm = nhstepm/hstepm;
Line 2373  if (h==(int)(calagedate+12*cpt)) for(j=1 Line 2384  if (h==(int)(calagedate+12*cpt)) for(j=1
           hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);              hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
           for (h=0; h<=nhstepm; h++){            for (h=0; h<=nhstepm; h++){
             if (h==(int) (calagedate+YEARM*cpt)) {              if (h==(int) (calagedate+YEARM*cpt)) {
               fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);                fprintf(ficresf,"\n %3.f ",agedeb+h*hstepm/YEARM*stepm);
             }              }
             for(j=1; j<=nlstate+ndeath;j++) {              for(j=1; j<=nlstate+ndeath;j++) {
               kk1=0.;kk2=0;                kk1=0.;kk2=0;
               for(i=1; i<=nlstate;i++) {                              for(i=1; i<=nlstate;i++) {              
                 kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod];                      kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod];    
               }                }
               if (h==(int)(calagedate+12*cpt)) fprintf(ficresf," %.3f", kk1);                  if (h==(int)(calagedate+12*cpt)) fprintf(ficresf," %15.2f", kk1);
             }              }
           }            }
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);            free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
Line 2410  int main(int argc, char *argv[]) Line 2421  int main(int argc, char *argv[])
   
   int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;    int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
   double agedeb, agefin,hf;    double agedeb, agefin,hf;
   double agemin=1.e20, agemax=-1.e20;    double agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;
   
   double fret;    double fret;
   double **xi,tmp,delta;    double **xi,tmp,delta;
Line 2423  int main(int argc, char *argv[]) Line 2434  int main(int argc, char *argv[])
   char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];    char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH];
   char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];    char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
     
   char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH];    char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
   
   char filerest[FILENAMELENGTH];    char filerest[FILENAMELENGTH];
   char fileregp[FILENAMELENGTH];    char fileregp[FILENAMELENGTH];
Line 2437  int main(int argc, char *argv[]) Line 2448  int main(int argc, char *argv[])
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;    int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
   int mobilav=0,popforecast=0;    int mobilav=0,popforecast=0;
   int hstepm, nhstepm;    int hstepm, nhstepm;
   double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2;    double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1;
   
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
Line 2455  int main(int argc, char *argv[]) Line 2466  int main(int argc, char *argv[])
   double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;    double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;
     
   
   char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";    char version[80]="Imach version 0.71e, March 2002, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
   
   
Line 2623  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2634  while((c=getc(ficpar))=='#' && c!= EOF){
   printf("\n");    printf("\n");
   
   
     /*-------- data file ----------*/      /*-------- Rewriting paramater file ----------*/
     if((ficres =fopen(fileres,"w"))==NULL) {       strcpy(rfileres,"r");    /* "Rparameterfile */
       printf("Problem with resultfile: %s\n", fileres);goto end;       strcat(rfileres,optionfilefiname);    /* Parameter file first name*/
        strcat(rfileres,".");    /* */
        strcat(rfileres,optionfilext);    /* Other files have txt extension */
       if((ficres =fopen(rfileres,"w"))==NULL) {
         printf("Problem writing new parameter file: %s\n", fileres);goto end;
     }      }
     fprintf(ficres,"#%s\n",version);      fprintf(ficres,"#%s\n",version);
         
       /*-------- data file ----------*/
     if((fic=fopen(datafile,"r"))==NULL)    {      if((fic=fopen(datafile,"r"))==NULL)    {
       printf("Problem with datafile: %s\n", datafile);goto end;        printf("Problem with datafile: %s\n", datafile);goto end;
     }      }
Line 2901  printf("Total number of individuals= %d, Line 2917  printf("Total number of individuals= %d,
     }      }
         
     /*--------- results files --------------*/      /*--------- results files --------------*/
     fprintf(ficres,"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);      fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, weightopt,model);
     
   
    jk=1;     jk=1;
Line 2973  printf("Total number of individuals= %d, Line 2989  printf("Total number of individuals= %d,
     }      }
     ungetc(c,ficpar);      ungetc(c,ficpar);
     
     fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);      fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemaxpar, &bage, &fage);
         
     if (fage <= 2) {      if (fage <= 2) {
       bage = agemin;        bage = agemin;
       fage = agemax;        fage = agemaxpar;
     }      }
         
     fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");      fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
     fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);      fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage);
     fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);      fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemaxpar,bage,fage);
     
     while((c=getc(ficpar))=='#' && c!= EOF){      while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
Line 3009  printf("Total number of individuals= %d, Line 3025  printf("Total number of individuals= %d,
    dateprev2=anprev2+mprev2/12.+jprev2/365.;     dateprev2=anprev2+mprev2/12.+jprev2/365.;
   
   fscanf(ficpar,"pop_based=%d\n",&popbased);    fscanf(ficpar,"pop_based=%d\n",&popbased);
    fprintf(ficparo,"pop_based=%d\n",popbased);      fprintf(ficparo,"pop_based=%d\n",popbased);  
    fprintf(ficres,"pop_based=%d\n",popbased);      fprintf(ficres,"pop_based=%d\n",popbased);  
    
   while((c=getc(ficpar))=='#' && c!= EOF){    while((c=getc(ficpar))=='#' && c!= EOF){
     ungetc(c,ficpar);      ungetc(c,ficpar);
     fgets(line, MAXLINE, ficpar);      fgets(line, MAXLINE, ficpar);
Line 3019  printf("Total number of individuals= %d, Line 3035  printf("Total number of individuals= %d,
     fputs(line,ficparo);      fputs(line,ficparo);
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
   fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mov_average=%d\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav);  
 fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);    fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mov_average=%d\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilav);
 fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);  fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
   fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mov_average=%d\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilav);
   
   
   while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
   
     fscanf(ficpar,"popforecast=%d popfile=%s popfiledate=%lf/%lf/%lf last-popfiledate=%lf/%lf/%lf\n",&popforecast,popfile,&jpyram,&mpyram,&anpyram,&jpyram1,&mpyram1,&anpyram1);
     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,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
   
 /*------------ gnuplot -------------*/  /*------------ gnuplot -------------*/
  printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemax,fage, pathc,p);   printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, agemin,agemaxpar,fage, pathc,p);
     
 /*------------ free_vector  -------------*/  /*------------ free_vector  -------------*/
  chdir(path);   chdir(path);
Line 3039  fprintf(ficres,"popforecast=%d popfile=% Line 3069  fprintf(ficres,"popforecast=%d popfile=%
  /*free_matrix(covar,1,NCOVMAX,1,n);*/   /*free_matrix(covar,1,NCOVMAX,1,n);*/
  fclose(ficparo);   fclose(ficparo);
  fclose(ficres);   fclose(ficres);
    
   /* Reads comments: lines beginning with '#' */  
   while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);  
     fgets(line, MAXLINE, ficpar);  
     puts(line);  
     fputs(line,ficparo);  
   }  
   ungetc(c,ficpar);  
    
   fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf\n",&agemin,&agemax, &bage, &fage);  
   printf("agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax, bage, fage);  
   fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);  
 /*--------- index.htm --------*/  /*--------- index.htm --------*/
   
   printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm);    printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,optionfile,optionfilehtm,rfileres);
   
     
   /*--------------- Prevalence limit --------------*/    /*--------------- Prevalence limit --------------*/
Line 3078  fprintf(ficres,"popforecast=%d popfile=% Line 3096  fprintf(ficres,"popforecast=%d popfile=%
   oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */    oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */
   k=0;    k=0;
   agebase=agemin;    agebase=agemin;
   agelim=agemax;    agelim=agemaxpar;
   ftolpl=1.e-10;    ftolpl=1.e-10;
   i1=cptcoveff;    i1=cptcoveff;
   if (cptcovn < 1){i1=1;}    if (cptcovn < 1){i1=1;}
Line 3157  fprintf(ficres,"popforecast=%d popfile=% Line 3175  fprintf(ficres,"popforecast=%d popfile=%
   
   
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
   if(stepm == 1) {    if((stepm == 1) && (strcmp(model,".")==0)){
     prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);      prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);
 populforecast(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);
     free_matrix(mint,1,maxwav,1,n);      free_matrix(mint,1,maxwav,1,n);
     free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);      free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n);
     free_vector(weight,1,n);}      free_vector(weight,1,n);}
   else{    else{
     erreur=108;      erreur=108;
     printf("Error %d!! You can only forecast the prevalences if the optimization\n  has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm);      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);
   }    }
     
   
Line 3225  populforecast(fileres, anproj1,mproj1,jp Line 3243  populforecast(fileres, anproj1,mproj1,jp
       for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);        for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
       fprintf(ficrest,"\n");        fprintf(ficrest,"\n");
   
       hf=1;  
       if (stepm >= YEARM) hf=stepm/YEARM;  
       epj=vector(1,nlstate+1);        epj=vector(1,nlstate+1);
       for(age=bage; age <=fage ;age++){        for(age=bage; age <=fage ;age++){
         prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);          prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
Line 3235  populforecast(fileres, anproj1,mproj1,jp Line 3251  populforecast(fileres, anproj1,mproj1,jp
             prlim[i][i]=probs[(int)age][i][k];              prlim[i][i]=probs[(int)age][i][k];
         }          }
                 
         fprintf(ficrest," %.0f",age);          fprintf(ficrest," %4.0f",age);
         for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){          for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
           for(i=1, epj[j]=0.;i <=nlstate;i++) {            for(i=1, epj[j]=0.;i <=nlstate;i++) {
             epj[j] += prlim[i][i]*hf*eij[i][j][(int)age];              epj[j] += prlim[i][i]*eij[i][j][(int)age];
           }            }
           epj[nlstate+1] +=epj[j];            epj[nlstate+1] +=epj[j];
         }          }
         for(i=1, vepp=0.;i <=nlstate;i++)          for(i=1, vepp=0.;i <=nlstate;i++)
           for(j=1;j <=nlstate;j++)            for(j=1;j <=nlstate;j++)
             vepp += vareij[i][j][(int)age];              vepp += vareij[i][j][(int)age];
         fprintf(ficrest," %.2f (%.2f)", epj[nlstate+1],hf*sqrt(vepp));          fprintf(ficrest," %7.2f (%7.2f)", epj[nlstate+1],sqrt(vepp));
         for(j=1;j <=nlstate;j++){          for(j=1;j <=nlstate;j++){
           fprintf(ficrest," %.2f (%.2f)", epj[j],hf*sqrt(vareij[j][j][(int)age]));            fprintf(ficrest," %7.2f (%7.2f)", epj[j],sqrt(vareij[j][j][(int)age]));
         }          }
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
       }        }
Removed from v.1.27  
changed lines
  Added in v.1.33

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