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

Diff for /imach/src/imach.c between versions 1.53 and 1.58

version 1.53, 2002/07/23 23:59:37 version 1.58, 2002/07/26 12:29:55
Line 39 Line 39
   hPijx.    hPijx.
   
   Also this programme outputs the covariance matrix of the parameters but also    Also this programme outputs the covariance matrix of the parameters but also
   of the life expectancies. It also computes the prevalence limits.     of the life expectancies. It also computes the stable prevalence. 
       
   Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).    Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
            Institut national d'études démographiques, Paris.             Institut national d'études démographiques, Paris.
Line 60 Line 60
 /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/  /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/
 #define FILENAMELENGTH 80  #define FILENAMELENGTH 80
 /*#define DEBUG*/  /*#define DEBUG*/
 #define unix  #define windows
 #define GLOCK_ERROR_NOPATH              -1      /* empty path */  #define GLOCK_ERROR_NOPATH              -1      /* empty path */
 #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */  #define GLOCK_ERROR_GETCWD              -2      /* cannot get cwd */
   
Line 83 Line 83
 #define ODIRSEPARATOR '\\'  #define ODIRSEPARATOR '\\'
 #endif  #endif
   
 char version[80]="Imach version 0.8j, July 2002, INED-EUROREVES ";  char version[80]="Imach version 0.8k, July 2002, INED-EUROREVES ";
 int erreur; /* Error number */  int erreur; /* Error number */
 int nvar;  int nvar;
 int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;  int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
Line 163  int estepm; Line 163  int estepm;
 int m,nb;  int m,nb;
 int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;  int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
 double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;  double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
 double **pmmij, ***probs, ***mobaverage;  double **pmmij, ***probs;
 double dateintmean=0;  double dateintmean=0;
   
 double *weight;  double *weight;
Line 701  void powell(double p[], double **xi, int Line 701  void powell(double p[], double **xi, int
         printf("\n");          printf("\n");
         fprintf(ficlog,"\n");          fprintf(ficlog,"\n");
 #endif  #endif
       }         }
     }       } 
   }     } 
 }   } 
   
 /**** Prevalence limit ****************/  /**** Prevalence limit (stable prevalence)  ****************/
   
 double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)  double **prevalim(double **prlim, int nlstate, double x[], double age, double **oldm, double **savm, double ftolpl, int ij)
 {  {
Line 1284  void  freqsummary(char fileres[], int ag Line 1284  void  freqsummary(char fileres[], int ag
             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.);
             if ((k2>=dateprev1) && (k2<=dateprev2)) {              if ((k2>=dateprev1) && (k2<=dateprev2)) {
Line 1575  void  concatwav(int wav[], int **dh, int Line 1575  void  concatwav(int wav[], int **dh, int
 /*********** Tricode ****************************/  /*********** Tricode ****************************/
 void tricode(int *Tvar, int **nbcode, int imx)  void tricode(int *Tvar, int **nbcode, int imx)
 {  {
   int Ndum[20],ij=1, k, j, i;    
     int Ndum[20],ij=1, k, j, i, maxncov=19;
   int cptcode=0;    int cptcode=0;
   cptcoveff=0;     cptcoveff=0; 
     
   for (k=0; k<19; k++) Ndum[k]=0;    for (k=0; k<maxncov; k++) Ndum[k]=0;
   for (k=1; k<=7; k++) ncodemax[k]=0;    for (k=1; k<=7; k++) ncodemax[k]=0;
   
   for (j=1; j<=(cptcovn+2*cptcovprod); j++) {    for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
     for (i=1; i<=imx; i++) {      for (i=1; i<=imx; i++) { /*reads the data file to get the maximum 
       ij=(int)(covar[Tvar[j]][i]);                                 modality*/ 
       Ndum[ij]++;         ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual*/
         Ndum[ij]++; /*store the modality */
       /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/        /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
       if (ij > cptcode) cptcode=ij;         if (ij > cptcode) cptcode=ij; /* getting the maximum of covariable 
                                          Tvar[j]. If V=sex and male is 0 and 
                                          female is 1, then  cptcode=1.*/
     }      }
   
     for (i=0; i<=cptcode; i++) {      for (i=0; i<=cptcode; i++) {
       if(Ndum[i]!=0) ncodemax[j]++;        if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */
     }      }
     ij=1;   
   
   
       ij=1; 
     for (i=1; i<=ncodemax[j]; i++) {      for (i=1; i<=ncodemax[j]; i++) {
       for (k=0; k<=19; k++) {        for (k=0; k<= maxncov; k++) {
         if (Ndum[k] != 0) {          if (Ndum[k] != 0) {
           nbcode[Tvar[j]][ij]=k;             nbcode[Tvar[j]][ij]=k; 
             /* store the modality in an array. k is a modality. If we have model=V1+V1*sex then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
                       
           ij++;            ij++;
         }          }
Line 1608  void tricode(int *Tvar, int **nbcode, in Line 1612  void tricode(int *Tvar, int **nbcode, in
     }       } 
   }      }  
   
  for (k=0; k<19; k++) Ndum[k]=0;   for (k=0; k< maxncov; k++) Ndum[k]=0;
   
  for (i=1; i<=ncovmodel-2; i++) {   for (i=1; i<=ncovmodel-2; i++) { 
      /* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
    ij=Tvar[i];     ij=Tvar[i];
    Ndum[ij]++;      Ndum[ij]++;
  }   }
   
  ij=1;   ij=1;
  for (i=1; i<=10; i++) {   for (i=1; i<= maxncov; i++) {
    if((Ndum[i]!=0) && (i<=ncovcol)){     if((Ndum[i]!=0) && (i<=ncovcol)){
      Tvaraff[ij]=i;        Tvaraff[ij]=i; /*For printing */
      ij++;       ij++;
    }     }
  }   }
     
  cptcoveff=ij-1;   cptcoveff=ij-1; /*Number of simple covariates*/
 }  }
   
 /*********** Health Expectancies ****************/  /*********** Health Expectancies ****************/
Line 1819  void varevsij(char optionfilefiname[], d Line 1824  void varevsij(char optionfilefiname[], d
   double ***mobaverage;    double ***mobaverage;
   int theta;    int theta;
   char digit[4];    char digit[4];
   char digitp[16];    char digitp[25];
   
   char fileresprobmorprev[FILENAMELENGTH];    char fileresprobmorprev[FILENAMELENGTH];
   
   if(popbased==1)    if(popbased==1){
     strcpy(digitp,"-populbased-");      if(mobilav!=0)
   else        strcpy(digitp,"-populbased-mobilav-");
       else strcpy(digitp,"-populbased-nomobil-");
     }
     else 
     strcpy(digitp,"-stablbased-");      strcpy(digitp,"-stablbased-");
   if(mobilav==1)  
     strcat(digitp,"mobilav-");    if (mobilav!=0) {
   else  
     strcat(digitp,"nomobil-");  
   if (mobilav==1) {  
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     movingaverage(probs, bage, fage, mobaverage);      if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
   }    }
   
   strcpy(fileresprobmorprev,"prmorprev");     strcpy(fileresprobmorprev,"prmorprev"); 
Line 1928  void varevsij(char optionfilefiname[], d Line 1936  void varevsij(char optionfilefiname[], d
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
   
       if (popbased==1) {        if (popbased==1) {
         if(mobilav !=1){          if(mobilav ==0){
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             prlim[i][i]=probs[(int)age][i][ij];              prlim[i][i]=probs[(int)age][i][ij];
         }else{ /* mobilav=1 */           }else{ /* mobilav */ 
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             prlim[i][i]=mobaverage[(int)age][i][ij];              prlim[i][i]=mobaverage[(int)age][i][ij];
         }          }
Line 1956  void varevsij(char optionfilefiname[], d Line 1964  void varevsij(char optionfilefiname[], d
       prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);        prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
     
       if (popbased==1) {        if (popbased==1) {
         if(mobilav !=1){          if(mobilav ==0){
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             prlim[i][i]=probs[(int)age][i][ij];              prlim[i][i]=probs[(int)age][i][ij];
         }else{ /* mobilav=1 */           }else{ /* mobilav */ 
           for(i=1; i<=nlstate;i++)            for(i=1; i<=nlstate;i++)
             prlim[i][i]=mobaverage[(int)age][i][ij];              prlim[i][i]=mobaverage[(int)age][i][ij];
         }          }
Line 2025  void varevsij(char optionfilefiname[], d Line 2033  void varevsij(char optionfilefiname[], d
     prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);      prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);
     
     if (popbased==1) {      if (popbased==1) {
       if(mobilav !=1){        if(mobilav ==0){
         for(i=1; i<=nlstate;i++)          for(i=1; i<=nlstate;i++)
           prlim[i][i]=probs[(int)age][i][ij];            prlim[i][i]=probs[(int)age][i][ij];
       }else{ /* mobilav=1 */         }else{ /* mobilav */ 
         for(i=1; i<=nlstate;i++)          for(i=1; i<=nlstate;i++)
           prlim[i][i]=mobaverage[(int)age][i][ij];            prlim[i][i]=mobaverage[(int)age][i][ij];
       }        }
Line 2084  void varevsij(char optionfilefiname[], d Line 2092  void varevsij(char optionfilefiname[], d
   free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);    free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
   free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);    free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);
   free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);    free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
   free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   fclose(ficresprobmorprev);    fclose(ficresprobmorprev);
   fclose(ficgp);    fclose(ficgp);
   fclose(fichtm);    fclose(fichtm);
   
 }  }
   
 /************ Variance of prevlim ******************/  /************ Variance of prevlim ******************/
Line 2106  void varprevlim(char fileres[], double * Line 2113  void varprevlim(char fileres[], double *
   double age,agelim;    double age,agelim;
   int theta;    int theta;
         
   fprintf(ficresvpl,"# Standard deviation of prevalence's limit\n");    fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");
   fprintf(ficresvpl,"# Age");    fprintf(ficresvpl,"# Age");
   for(i=1; i<=nlstate;i++)    for(i=1; i<=nlstate;i++)
       fprintf(ficresvpl," %1d-%1d",i,i);        fprintf(ficresvpl," %1d-%1d",i,i);
Line 2599  void printinggnuplot(char fileres[], dou Line 2606  void printinggnuplot(char fileres[], dou
     fprintf(ficlog,"Problem with file %s",optionfilegnuplot);      fprintf(ficlog,"Problem with file %s",optionfilegnuplot);
   }    }
   
 #ifdef windows    /*#ifdef windows */
     fprintf(ficgp,"cd \"%s\" \n",pathc);      fprintf(ficgp,"cd \"%s\" \n",pathc);
 #endif      /*#endif */
 m=pow(2,cptcoveff);  m=pow(2,cptcoveff);
       
  /* 1eme*/   /* 1eme*/
Line 2614  m=pow(2,cptcoveff); Line 2621  m=pow(2,cptcoveff);
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");         else fprintf(ficgp," \%%*lf (\%%*lf)");
      }       }
      fprintf(ficgp,"\" t\"Stationary prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);       fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);
      for (i=1; i<= nlstate ; i ++) {       for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");         else fprintf(ficgp," \%%*lf (\%%*lf)");
Line 2768  m=pow(2,cptcoveff); Line 2775  m=pow(2,cptcoveff);
   
   
 /*************** Moving average **************/  /*************** Moving average **************/
 void movingaverage(double ***probs, double bage,double fage, double ***mobaverage){  int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){
   
   int i, cpt, cptcod;    int i, cpt, cptcod;
     int modcovmax =1;
     int mobilavrange, mob;
   double age;    double age;
   for (age=bage; age<=fage; age++)  
     for (i=1; i<=nlstate;i++)    modcovmax=2*cptcoveff;/* Max number of modalities. We suppose 
       for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)                             a covariate has 2 modalities */
         mobaverage[(int)age][i][cptcod]=0.;    if (cptcovn<1) modcovmax=1; /* At least 1 pass */
     
   for (age=bage+4; age<=fage; age++){    if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){
     for (i=1; i<=nlstate;i++){      if(mobilav==1) mobilavrange=5; /* default */
       for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){      else mobilavrange=mobilav;
         for (cpt=0;cpt<=4;cpt++){      for (age=bage; age<=fage; age++)
           mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]+probs[(int)age-cpt][i][cptcod];        for (i=1; i<=nlstate;i++)
           for (cptcod=1;cptcod<=modcovmax;cptcod++)
             mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod];
       /* We keep the original values on the extreme ages bage, fage and for 
          fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2
          we use a 5 terms etc. until the borders are no more concerned. 
       */ 
       for (mob=3;mob <=mobilavrange;mob=mob+2){
         for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){
           for (i=1; i<=nlstate;i++){
             for (cptcod=1;cptcod<=modcovmax;cptcod++){
               mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod];
                 for (cpt=1;cpt<=(mob-1)/2;cpt++){
                   mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod];
                   mobaverage[(int)age][i][cptcod] +=probs[(int)age+cpt][i][cptcod];
                 }
               mobaverage[(int)age][i][cptcod]=mobaverage[(int)age][i][cptcod]/mob;
             }
         }          }
         mobaverage[(int)age-2][i][cptcod]=mobaverage[(int)age-2][i][cptcod]/5;        }/* end age */
       }      }/* end mob */
     }    }else return -1;
   }    return 0;
     }/* End movingaverage */
 }  
   
   
 /************** Forecasting ******************/  /************** Forecasting ******************/
Line 2799  prevforecast(char fileres[], double anpr Line 2824  prevforecast(char fileres[], double anpr
   double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;    double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
   double *popeffectif,*popcount;    double *popeffectif,*popcount;
   double ***p3mat;    double ***p3mat;
     double ***mobaverage;
   char fileresf[FILENAMELENGTH];    char fileresf[FILENAMELENGTH];
   
  agelim=AGESUP;   agelim=AGESUP;
 calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;   calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
   
   prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);    prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
     
Line 2818  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2844  calagedate=(anproj1+mproj1/12.+jproj1/36
   
   if (cptcoveff==0) ncodemax[cptcoveff]=1;    if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
   if (mobilav==1) {    if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     movingaverage(probs, ageminpar,fage, mobaverage);      if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
   }    }
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;    stepsize=(int) (stepm+YEARM-1)/YEARM;
Line 2891  calagedate=(anproj1+mproj1/12.+jproj1/36 Line 2920  calagedate=(anproj1+mproj1/12.+jproj1/36
     }      }
   }    }
                 
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   fclose(ficresf);    fclose(ficresf);
 }  }
Line 2903  populforecast(char fileres[], double anp Line 2932  populforecast(char fileres[], double anp
   double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;    double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
   double *popeffectif,*popcount;    double *popeffectif,*popcount;
   double ***p3mat,***tabpop,***tabpopprev;    double ***p3mat,***tabpop,***tabpopprev;
     double ***mobaverage;
   char filerespop[FILENAMELENGTH];    char filerespop[FILENAMELENGTH];
   
   tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
Line 2924  populforecast(char fileres[], double anp Line 2954  populforecast(char fileres[], double anp
   
   if (cptcoveff==0) ncodemax[cptcoveff]=1;    if (cptcoveff==0) ncodemax[cptcoveff]=1;
   
   if (mobilav==1) {    if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     movingaverage(probs, ageminpar, fage, mobaverage);      if (movingaverage(probs, ageminpar, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
   }    }
   
   stepsize=(int) (stepm+YEARM-1)/YEARM;    stepsize=(int) (stepm+YEARM-1)/YEARM;
Line 3039  populforecast(char fileres[], double anp Line 3072  populforecast(char fileres[], double anp
    }      } 
   }    }
     
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   if (popforecast==1) {    if (popforecast==1) {
     free_ivector(popage,0,AGESUP);      free_ivector(popage,0,AGESUP);
Line 3187  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3220  while((c=getc(ficpar))=='#' && c!= EOF){
       
         
   covar=matrix(0,NCOVMAX,1,n);     covar=matrix(0,NCOVMAX,1,n); 
   cptcovn=0;     cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/
   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;    if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
   
   ncovmodel=2+cptcovn;    ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */
   nvar=ncovmodel-1; /* Suppressing age as a basic covariate */    nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
       
   /* Read guess parameters */    /* Read guess parameters */
Line 3228  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3261  while((c=getc(ficpar))=='#' && c!= EOF){
       fprintf(ficparo,"\n");        fprintf(ficparo,"\n");
     }      }
       
     npar= (nlstate+ndeath-1)*nlstate*ncovmodel;      npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
   
   p=param[1][1];    p=param[1][1];
       
Line 3387  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3420  while((c=getc(ficpar))=='#' && c!= EOF){
   Tvard=imatrix(1,15,1,2);    Tvard=imatrix(1,15,1,2);
   Tage=ivector(1,15);          Tage=ivector(1,15);      
         
   if (strlen(model) >1){    if (strlen(model) >1){ /* If there is at least 1 covariate */
     j=0, j1=0, k1=1, k2=1;      j=0, j1=0, k1=1, k2=1;
     j=nbocc(model,'+');      j=nbocc(model,'+'); /* j=Number of '+' */
     j1=nbocc(model,'*');      j1=nbocc(model,'*'); /* j1=Number of '*' */
     cptcovn=j+1;      cptcovn=j+1; 
     cptcovprod=j1;      cptcovprod=j1; /*Number of products */
           
     strcpy(modelsav,model);       strcpy(modelsav,model); 
     if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){      if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
Line 3401  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3434  while((c=getc(ficpar))=='#' && c!= EOF){
       goto end;        goto end;
     }      }
           
       /* This loop fill the array Tvar from the string 'model'.*/
   
     for(i=(j+1); i>=1;i--){      for(i=(j+1); i>=1;i--){
       cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */         cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ 
       if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */        if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */
Line 3450  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3485  while((c=getc(ficpar))=='#' && c!= EOF){
     } /* end of loop + */      } /* end of loop + */
   } /* end model */    } /* end model */
       
     /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products.
       If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/
   
   /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);    /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
   printf("cptcovprod=%d ", cptcovprod);    printf("cptcovprod=%d ", cptcovprod);
   fprintf(ficlog,"cptcovprod=%d ", cptcovprod);    fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
   scanf("%d ",i);*/  
     fclose(fic);    scanf("%d ",i);
     fclose(fic);*/
   
     /*  if(mle==1){*/      /*  if(mle==1){*/
     if (weightopt != 1) { /* Maximisation without weights*/      if (weightopt != 1) { /* Maximisation without weights*/
Line 3525  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3564  while((c=getc(ficpar))=='#' && c!= EOF){
       }        }
     }      }
   
 printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);      printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
  fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);       fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); 
   
     free_vector(severity,1,maxwav);      free_vector(severity,1,maxwav);
     free_imatrix(outcome,1,maxwav+1,1,n);      free_imatrix(outcome,1,maxwav+1,1,n);
Line 3545  printf("Total number of individuals= %d, Line 3584  printf("Total number of individuals= %d,
     /* Concatenates waves */      /* Concatenates waves */
       concatwav(wav, dh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);        concatwav(wav, dh, mw, s, agedc, agev,  firstpass, lastpass, imx, nlstate, stepm);
   
         /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
   
       Tcode=ivector(1,100);        Tcode=ivector(1,100);
       nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);         nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); 
       ncodemax[1]=1;        ncodemax[1]=1;
       if (cptcovn > 0) tricode(Tvar,nbcode,imx);        if (cptcovn > 0) tricode(Tvar,nbcode,imx);
               
    codtab=imatrix(1,100,1,10);        codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of 
                                        the estimations*/
    h=0;     h=0;
    m=pow(2,cptcoveff);     m=pow(2,cptcoveff);
     
Line 3704  printf("Total number of individuals= %d, Line 3745  printf("Total number of individuals= %d,
    ungetc(c,ficpar);     ungetc(c,ficpar);
       
    fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);     fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
    fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,&mobilav);     fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
    fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,&mobilav);     fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
         
    while((c=getc(ficpar))=='#' && c!= EOF){     while((c=getc(ficpar))=='#' && c!= EOF){
      ungetc(c,ficpar);       ungetc(c,ficpar);
Line 3748  while((c=getc(ficpar))=='#' && c!= EOF){ Line 3789  while((c=getc(ficpar))=='#' && c!= EOF){
   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(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);    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 -------------*/
   strcpy(optionfilegnuplot,optionfilefiname);   strcpy(optionfilegnuplot,optionfilefiname);
   strcat(optionfilegnuplot,".gp");   strcat(optionfilegnuplot,".gp");
   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {   if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
     printf("Problem with file %s",optionfilegnuplot);     printf("Problem with file %s",optionfilegnuplot);
   }   }
   fclose(ficgp);   else{
      fprintf(ficgp,"\n# %s\n", version); 
      fprintf(ficgp,"# %s\n", optionfilegnuplot); 
      fprintf(ficgp,"set missing 'NaNq'\n");
   }
    fclose(ficgp);
  printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);   printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
 /*--------- index.htm --------*/  /*--------- index.htm --------*/
   
Line 3792  Interval (in months) between two waves: Line 3839  Interval (in months) between two waves:
  fclose(ficres);   fclose(ficres);
   
   
   /*--------------- Prevalence limit --------------*/    /*--------------- Prevalence limit  (stable prevalence) --------------*/
       
   strcpy(filerespl,"pl");    strcpy(filerespl,"pl");
   strcat(filerespl,fileres);    strcat(filerespl,fileres);
   if((ficrespl=fopen(filerespl,"w"))==NULL) {    if((ficrespl=fopen(filerespl,"w"))==NULL) {
     printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;      printf("Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
     fprintf(ficlog,"Problem with Prev limit resultfile: %s\n", filerespl);goto end;      fprintf(ficlog,"Problem with stable prevalence resultfile: %s\n", filerespl);goto end;
   }    }
   printf("Computing prevalence limit: result on file '%s' \n", filerespl);    printf("Computing stable prevalence: result on file '%s' \n", filerespl);
   fprintf(ficlog,"Computing prevalence limit: result on file '%s' \n", filerespl);    fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);
   fprintf(ficrespl,"#Prevalence limit\n");    fprintf(ficrespl,"#Stable prevalence \n");
   fprintf(ficrespl,"#Age ");    fprintf(ficrespl,"#Age ");
   for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);    for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
   fprintf(ficrespl,"\n");    fprintf(ficrespl,"\n");
Line 3948  Interval (in months) between two waves: Line 3995  Interval (in months) between two waves:
   }    }
   printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);    printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);    fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
   
   calagedate=-1;    calagedate=-1;
   
   prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);    prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
   if (mobilav==1) {  
     if (mobilav!=0) {
     mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);      mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
     movingaverage(probs, ageminpar, fage, mobaverage);      if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){
         fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
         printf(" Error in movingaverage mobilav=%d\n",mobilav);
       }
   }    }
   
   k=0;    k=0;
Line 3994  Interval (in months) between two waves: Line 4047  Interval (in months) between two waves:
       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);
         if (popbased==1) {          if (popbased==1) {
           if(mobilav !=1){            if(mobilav ==0){
             for(i=1; i<=nlstate;i++)              for(i=1; i<=nlstate;i++)
               prlim[i][i]=probs[(int)age][i][k];                prlim[i][i]=probs[(int)age][i][k];
           }else{ /* mobilav=1 */             }else{ /* mobilav */ 
             for(i=1; i<=nlstate;i++)              for(i=1; i<=nlstate;i++)
               prlim[i][i]=mobaverage[(int)age][i][k];                prlim[i][i]=mobaverage[(int)age][i][k];
           }            }
Line 4032  free_matrix(mint,1,maxwav,1,n); Line 4085  free_matrix(mint,1,maxwav,1,n);
   fclose(ficpar);    fclose(ficpar);
   free_vector(epj,1,nlstate+1);    free_vector(epj,1,nlstate+1);
       
   /*------- Variance limit prevalence------*/       /*------- Variance of stable prevalence------*/   
   
   strcpy(fileresvpl,"vpl");    strcpy(fileresvpl,"vpl");
   strcat(fileresvpl,fileres);    strcat(fileresvpl,fileres);
   if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {    if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
     printf("Problem with variance prev lim resultfile: %s\n", fileresvpl);      printf("Problem with variance of stable prevalence  resultfile: %s\n", fileresvpl);
     exit(0);      exit(0);
   }    }
   printf("Computing Variance-covariance of Prevalence limit: file '%s' \n", fileresvpl);    printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl);
   
   k=0;    k=0;
   for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
Line 4075  free_matrix(mint,1,maxwav,1,n); Line 4128  free_matrix(mint,1,maxwav,1,n);
   free_vector(delti,1,npar);    free_vector(delti,1,npar);
   free_matrix(agev,1,maxwav,1,imx);    free_matrix(agev,1,maxwav,1,imx);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);    free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);    if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
   
   fprintf(fichtm,"\n</body>");    fprintf(fichtm,"\n</body>");
   fclose(fichtm);    fclose(fichtm);
Line 4109  free_matrix(mint,1,maxwav,1,n); Line 4162  free_matrix(mint,1,maxwav,1,n);
  strcpy(plotcmd,GNUPLOTPROGRAM);   strcpy(plotcmd,GNUPLOTPROGRAM);
  strcat(plotcmd," ");   strcat(plotcmd," ");
  strcat(plotcmd,optionfilegnuplot);   strcat(plotcmd,optionfilegnuplot);
    printf("Starting: %s\n",plotcmd);fflush(stdout);
  system(plotcmd);   system(plotcmd);
   
 #ifdef windows   /*#ifdef windows*/
   while (z[0] != 'q') {    while (z[0] != 'q') {
     /* chdir(path); */      /* chdir(path); */
     printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");      printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: ");
Line 4121  free_matrix(mint,1,maxwav,1,n); Line 4175  free_matrix(mint,1,maxwav,1,n);
     else if (z[0] == 'g') system(plotcmd);      else if (z[0] == 'g') system(plotcmd);
     else if (z[0] == 'q') exit(0);      else if (z[0] == 'q') exit(0);
   }    }
 #endif     /*#endif */
 }  }
   
   
Removed from v.1.53  
changed lines
  Added in v.1.58

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