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

Diff for /imach/src/imach.c between versions 1.18 and 1.19

version 1.18, 2002/02/20 17:17:09 version 1.19, 2002/02/20 17:19:10
Line 129  int m,nb; Line 129  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, ***mobaverage;
   double dateintmean=0;
   
 double *weight;  double *weight;
 int **s; /* Status */  int **s; /* Status */
Line 1150  void lubksb(double **a, int n, int *indx Line 1151  void lubksb(double **a, int n, int *indx
 }  }
   
 /************ Frequencies ********************/  /************ Frequencies ********************/
 void  freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1,double **mint,double **anint, int boolprev, double dateprev1,double dateprev2)  void  freqsummary(char fileres[], 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)
 {  /* 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;
   double ***freq; /* Frequencies */    double ***freq; /* Frequencies */
   double *pp;    double *pp;
   double pos, k2;    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,130 ,1,8, 1,8);    probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
   strcpy(fileresp,"p");    strcpy(fileresp,"p");
   strcat(fileresp,fileres);    strcat(fileresp,fileres);
   if((ficresp=fopen(fileresp,"w"))==NULL) {    if((ficresp=fopen(fileresp,"w"))==NULL) {
Line 1183  void  freqsummary(char fileres[], int ag Line 1184  void  freqsummary(char fileres[], int ag
          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;
          
           dateintsum=0;
           k2cpt=0;
        for (i=1; i<=imx; i++) {         for (i=1; i<=imx; i++) {
          bool=1;           bool=1;
          if  (cptcovn>0) {           if  (cptcovn>0) {
Line 1192  void  freqsummary(char fileres[], int ag Line 1195  void  freqsummary(char fileres[], int ag
                bool=0;                 bool=0;
          }           }
          if (bool==1) {           if (bool==1) {
            if (boolprev==1){             for(m=firstpass; m<=lastpass; m++){
              for(m=fprev1; m<=lprev1; m++){               k2=anint[m][i]+(mint[m][i]/12.);
                if ((k2>=dateprev1) && (k2<=dateprev2)) {
                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]] += weight[i];                 freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
                freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];                 freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
                  if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {
                    dateintsum=dateintsum+k2;
                    k2cpt++;
                  }
   
              }               }
            }             }
            else {           }
             for(m=firstpass; m<=lastpass; m++){  
              k2=anint[m][i]+(mint[m][i]/12.);  
              if ((k2>=dateprev1) && (k2<=dateprev2)) {  
              if(agev[m][i]==0) agev[m][i]=agemax+1;  
              if(agev[m][i]==1) agev[m][i]=agemax+2;  
              freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];  
              freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];  
              }  
             }  
            }  
           }  
        }         }
         if  (cptcovn>0) {          if  (cptcovn>0) {
          fprintf(ficresp, "\n#********** Variable ");           fprintf(ficresp, "\n#********** Variable ");
Line 1271  void  freqsummary(char fileres[], int ag Line 1269  void  freqsummary(char fileres[], int ag
     }      }
     }      }
  }   }
     dateintmean=dateintsum/k2cpt;
     
   fclose(ficresp);    fclose(ficresp);
   free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);    free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
   
 }  /* End of Freq */    /* End of Freq */
   }
   
 /************ Prevalence ********************/  /************ Prevalence ********************/
 void prevalence(int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax, int fprev1,int lprev1, double **mint,double **anint,int boolprev, double dateprev1, double dateprev2)  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)
 {  /* 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 1288  void prevalence(int agemin, int agemax, Line 1288  void prevalence(int agemin, int agemax,
   double pos, k2;    double pos, k2;
   
   pp=vector(1,nlstate);    pp=vector(1,nlstate);
   probs= ma3x(1,130 ,1,8, 1,8);    probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
     
   freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);    freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
   j1=0;    j1=0;
Line 1303  void prevalence(int agemin, int agemax, Line 1303  void prevalence(int agemin, int agemax,
       for (i=-1; i<=nlstate+ndeath; i++)          for (i=-1; i<=nlstate+ndeath; i++)  
         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;
Line 1311  void prevalence(int agemin, int agemax, Line 1311  void prevalence(int agemin, int agemax,
           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) {
           if (boolprev==1){            for(m=firstpass; m<=lastpass; m++){
             for(m=fprev1; m<=lprev1; m++){              k2=anint[m][i]+(mint[m][i]/12.);
               if ((k2>=dateprev1) && (k2<=dateprev2)) {
               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]] += weight[i];                freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-1/12.)] += weight[i];
               freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];                freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];  
             }  
           }  
           else {  
             for(m=firstpass; m<=lastpass; m++){  
               k2=anint[m][i]+(mint[m][i]/12.);  
               if ((k2>=dateprev1) && (k2<=dateprev2)) {  
                 if(agev[m][i]==0) agev[m][i]=agemax+1;  
                 if(agev[m][i]==1) agev[m][i]=agemax+2;  
                 freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];  
                 freq[s[m][i]][s[m+1][i]][(int) agemax+3] += 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 1368  void prevalence(int agemin, int agemax, Line 1359  void prevalence(int agemin, int agemax,
   free_vector(pp,1,nlstate);    free_vector(pp,1,nlstate);
     
 }  /* End of Freq */  }  /* End of Freq */
   
 /************* Waves Concatenation ***************/  /************* Waves Concatenation ***************/
   
 void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)  void  concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int  firstpass, int lastpass, int imx, int nlstate, int stepm)
Line 1894  int main() Line 1886  int main()
   int ju,jl, mi;    int ju,jl, mi;
   int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;    int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
   int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;    int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
   int mobilav=0, fprev, lprev ,fprevfore=1, lprevfore=1,nforecast,popforecast=0;    int mobilav=0,popforecast=0;
   int hstepm, nhstepm;    int hstepm, nhstepm;
   int *popage,boolprev=0;/*boolprev=0 if date and zero if wave*/    int *popage;/*boolprev=0 if date and zero if wave*/
     double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2;
   
   double bage, fage, age, agelim, agebase;    double bage, fage, age, agelim, agebase;
   double ftolpl=FTOL;    double ftolpl=FTOL;
Line 1912  int main() Line 1905  int main()
   double *epj, vepp;    double *epj, vepp;
   double kk1, kk2;    double kk1, kk2;
   double *popeffectif,*popcount;    double *popeffectif,*popcount;
   double dateprev1, dateprev2;    double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate;
     double yp,yp1,yp2;
   
   char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";    char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";
   char *alph[]={"a","a","b","c","d","e"}, str[4];    char *alph[]={"a","a","b","c","d","e"}, str[4];
Line 1922  int main() Line 1916  int main()
 #include <sys/time.h>  #include <sys/time.h>
 #include <time.h>  #include <time.h>
   char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];    char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80];
   char strfprev[10], strlprev[10];   
   char strfprevfore[10], strlprevfore[10];  
   /* long total_usecs;    /* long total_usecs;
   struct timeval start_time, end_time;    struct timeval start_time, end_time;
     
Line 1986  while((c=getc(ficpar))=='#' && c!= EOF){ Line 1979  while((c=getc(ficpar))=='#' && c!= EOF){
   }    }
   ungetc(c,ficpar);    ungetc(c,ficpar);
     
   fscanf(ficpar,"fprevalence=%s lprevalence=%s pop_based=%d\n",strfprev,strlprev,&popbased);     
   fprintf(ficparo,"fprevalence=%s lprevalence=%s pop_based=%d\n",strfprev,strlprev,popbased);  
    
   /* printf("%s %s",strfprev,strlprev);  
      exit(0);*/  
  while((c=getc(ficpar))=='#' && c!= EOF){  
     ungetc(c,ficpar);  
     fgets(line, MAXLINE, ficpar);  
     puts(line);  
     fputs(line,ficparo);  
   }  
   ungetc(c,ficpar);  
    
   fscanf(ficpar,"fprevalence=%s lprevalence=%s nforecast=%d mob_average=%d\n",strfprevfore,strlprevfore,&nforecast,&mobilav);  
   fprintf(ficparo,"fprevalence=%s lprevalence=%s nforecast=%d mob_average=%d\n",strfprevfore,strlprevfore,nforecast,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\n",&popforecast,popfile);  
    
   covar=matrix(0,NCOVMAX,1,n);    covar=matrix(0,NCOVMAX,1,n);
   cptcovn=0;    cptcovn=0;
   if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;    if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
Line 2175  while((c=getc(ficpar))=='#' && c!= EOF){ Line 2142  while((c=getc(ficpar))=='#' && c!= EOF){
     if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;      if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
     if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;      if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
     }      }
     for (i=1; i<=imx; i++) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/  
       for (i=1; i<=imx; i++)
       if (covar[1][i]==0) printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]),  (mint[2][i]), (anint[2][i]), (s[2][i]),  (mint[3][i]), (anint[3][i]), (s[3][i]),  (mint[4][i]), (anint[4][i]), (s[4][i]));*/
   
   /* Calculation of the number of parameter from char model*/    /* Calculation of the number of parameter from char model*/
   Tvar=ivector(1,15);    Tvar=ivector(1,15);
Line 2358  printf("Total number of individuals= %d, Line 2327  printf("Total number of individuals= %d,
    /* Calculates basic frequencies. Computes observed prevalence at single age     /* Calculates basic frequencies. Computes observed prevalence at single age
        and prints on file fileres'p'. */         and prints on file fileres'p'. */
   
     if ((nbocc(strfprev,'/')==1) && (nbocc(strlprev,'/')==1)){     
      boolprev=0;  
      cutv(stra,strb,strfprev,'/');  
      dateprev1=(double)(atoi(strb)+atoi(stra)/12.);  
      cutv(stra,strb,strlprev,'/');  
      dateprev2=(double)(atoi(strb)+atoi(stra)/12.);  
    }  
      
    else if ((nbocc(strfprev,'/')==0) &&(nbocc(strlprev,'/')==0)){  
      boolprev=1;  
      fprev=atoi(strfprev); lprev=atoi(strlprev);  
    }  
     else {  
       printf("Error in statement lprevalence or fprevalence\n");  
       goto end;  
     }  
         
   freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprev, lprev,mint,anint,boolprev,dateprev1,dateprev2);      pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
    
   pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */  
     oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
     savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */      savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
Line 2393  printf("Total number of individuals= %d, Line 2345  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=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);
    fprintf(ficres,"fprevalence=%s lprevalence=%s pop_based=%d\n",strfprev,strlprev,popbased);   
    fprintf(ficres,"fprevalence=%s lprevalence=%s nforecast=%d mob_average=%d\n",strfprevfore,strlprevfore,nforecast,mobilav);  
   
    jk=1;     jk=1;
    fprintf(ficres,"# Parameters\n");     fprintf(ficres,"# Parameters\n");
Line 2472  printf("Total number of individuals= %d, Line 2423  printf("Total number of individuals= %d,
       fage = agemax;        fage = agemax;
     }      }
   
     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.\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,agemax,bage,fage);
       fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
    
       while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
     fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mob_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
     fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
    fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
        
     while((c=getc(ficpar))=='#' && c!= EOF){
       ungetc(c,ficpar);
       fgets(line, MAXLINE, ficpar);
       puts(line);
       fputs(line,ficparo);
     }
     ungetc(c,ficpar);
    
   
         dateprev1=anprev1+mprev1/12.+jprev1/365.;
 /*------------ gnuplot -------------*/     dateprev2=anprev2+mprev2/12.+jprev2/365.;
   
     fscanf(ficpar,"pop_based=%d\n",&popbased);
      fprintf(ficparo,"pop_based=%d\n",popbased);  
      fprintf(ficres,"pop_based=%d\n",popbased);  
   
     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 starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2);
   fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
   fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
   
    freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2);
   
    /*------------ gnuplot -------------*/
 chdir(pathcd);  chdir(pathcd);
   if((ficgp=fopen("graph.plt","w"))==NULL) {    if((ficgp=fopen("graph.plt","w"))==NULL) {
     printf("Problem with file graph.gp");goto end;      printf("Problem with file graph.gp");goto end;
Line 2832  fclose(fichtm); Line 2825  fclose(fichtm);
   fclose(ficrespij);    fclose(ficrespij);
   
   /*---------- Forecasting ------------------*/    /*---------- Forecasting ------------------*/
     calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
   
     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 2840  fclose(fichtm); Line 2837  fclose(fichtm);
   }    }
   printf("Computing forecasting: result on file '%s' \n", fileresf);    printf("Computing forecasting: result on file '%s' \n", fileresf);
   
  if ((nbocc(strfprevfore,'/')==1) && (nbocc(strlprevfore,'/')==1)){  
      boolprev=0;  
      cutv(stra,strb,strfprevfore,'/');  
      dateprev1=(double)(atoi(strb)+atoi(stra)/12.);  
      cutv(stra,strb,strlprevfore,'/');  
      dateprev2=(double)(atoi(strb)+atoi(stra)/12.);  
    }  
      
    else if ((nbocc(strfprevfore,'/')==0) &&(nbocc(strlprevfore,'/')==0)){  
      boolprev=1;  
      fprev=atoi(strfprevfore); lprev=atoi(strlprevfore);  
    }  
     else {  
       printf("Error in statement lprevalence or fprevalence\n");  
       goto end;  
     }  
   
   prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax, fprevfore, lprevfore,mint,anint,boolprev,dateprev1,dateprev2);  
    
   free_matrix(mint,1,maxwav,1,n);    free_matrix(mint,1,maxwav,1,n);
   free_matrix(anint,1,maxwav,1,n);    free_matrix(anint,1,maxwav,1,n);
   free_matrix(agev,1,maxwav,1,imx);    free_matrix(agev,1,maxwav,1,imx);
Line 2889  fclose(fichtm); Line 2867  fclose(fichtm);
   if (stepm<=12) stepsize=1;    if (stepm<=12) stepsize=1;
   
   agelim=AGESUP;    agelim=AGESUP;
   hstepm=stepsize*YEARM; /* Every year of age */    /*hstepm=stepsize*YEARM; *//* Every year of age */
     hstepm=1;
   hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */    hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */
      yp1=modf(dateintmean,&yp);
     anprojmean=yp;
     yp2=modf((yp1*12),&yp);
     mprojmean=yp;
     yp1=modf((yp2*30.5),&yp);
     jprojmean=yp;
     fprintf(ficresf,"Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean);
   
   if (popforecast==1) {    if (popforecast==1) {
     if((ficpop=fopen(popfile,"r"))==NULL)    {      if((ficpop=fopen(popfile,"r"))==NULL)    {
       printf("Problem with population file : %s\n",popfile);goto end;        printf("Problem with population file : %s\n",popfile);goto end;
Line 2906  fclose(fichtm); Line 2892  fclose(fichtm);
         i=i+1;          i=i+1;
       }        }
     imx=i;      imx=i;
       
   for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];      for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
   }    }
   
   for(cptcov=1;cptcov<=i1;cptcov++){    for(cptcov=1;cptcov<=i1;cptcov++){
     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){      for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
       k=k+1;        k=k+1;
       fprintf(ficresf,"\n#****** ");        fprintf(ficresf,"\n#******");
       for(j=1;j<=cptcoveff;j++) {        for(j=1;j<=cptcoveff;j++) {
         fprintf(ficresf,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);          fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
       }        }
       fprintf(ficresf,"******\n");        fprintf(ficresf,"******\n");
       fprintf(ficresf,"# StartingAge FinalAge Horizon(in years)");        fprintf(ficresf,"# StartingAge FinalAge");
       for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);        for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
       if (popforecast==1)  fprintf(ficresf," [Population]");        if (popforecast==1)  fprintf(ficresf," [Population]");
    
       for (agedeb=fage; agedeb>=bage; agedeb--){        for (cpt=0; cpt<=1;cpt++) {
         fprintf(ficresf,"\n%.f %.f 0",agedeb, agedeb);  
        if (mobilav==1) {  
         for(j=1; j<=nlstate;j++)  
           fprintf(ficresf," %.3f",mobaverage[(int)agedeb][j][cptcod]);  
         }  
         else {  
           for(j=1; j<=nlstate;j++)  
           fprintf(ficresf," %.3f",probs[(int)agedeb][j][cptcod]);  
         }    
   
        for(j=1; j<=ndeath;j++) fprintf(ficresf," 0.00000");  
        if (popforecast==1) fprintf(ficresf," [%.f] ",popeffectif[(int)agedeb]);  
       }  
        
       for (cpt=1; cpt<=nforecast;cpt++) {  
         fprintf(ficresf,"\n");          fprintf(ficresf,"\n");
       for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */    fprintf(ficresf,"\nForecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);  
         for (agedeb=(fage-(1/12.)); agedeb>=(bage-(1/12.)); agedeb--){ /* If stepm=6 months */
         nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);          nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
         nhstepm = nhstepm/hstepm;          nhstepm = nhstepm/hstepm;
         /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/          /*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/
Line 2949  fclose(fichtm); Line 2921  fclose(fichtm);
         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+12*cpt)) {
          if (h*hstepm/YEARM*stepm==cpt)              fprintf(ficresf,"h=%d ", h);
             fprintf(ficresf,"\n%.f %.f %.f",agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);              fprintf(ficresf,"\n %f %f ",agedeb,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++) {                        if (mobilav==1)
              if (mobilav==1)  
                 kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];                  kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];
              else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];                else {
              if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];                  kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
                   /*  fprintf(ficresf," p3=%.3f p=%.3f ", p3mat[i][j][h],probs[(int)(agedeb)+1][i][cptcod]);*/
                 }
   
                 if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];
               }
            
               if (h==(int)(calagedate+12*cpt)){
                 fprintf(ficresf," %.3f", kk1);
                
                 if (popforecast==1) fprintf(ficresf," [%.f]", kk2);
             }              }
            if (h*hstepm/YEARM*stepm==cpt) {  
              fprintf(ficresf," %.3f", kk1);  
                if (popforecast==1) fprintf(ficresf," [%.f]", kk2);  
            }  
           }            }
         }          }
         free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);          free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
          
       }        }
       }        }
     }      }
Removed from v.1.18  
changed lines
  Added in v.1.19

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