Diff for /imach/src/imach.c between versions 1.266 and 1.267

version 1.266, 2017/05/13 07:26:12 version 1.267, 2017/05/13 10:25:05
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.267  2017/05/13 10:25:05  brouard
     Summary: temporary save for backprojection
   
   Revision 1.266  2017/05/13 07:26:12  brouard    Revision 1.266  2017/05/13 07:26:12  brouard
   Summary: Version 0.99r13 (improvements and bugs fixed)    Summary: Version 0.99r13 (improvements and bugs fixed)
   
Line 818  Back prevalence and projections: Line 821  Back prevalence and projections:
    p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);     p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
                         oldm=oldms;savm=savms;                          oldm=oldms;savm=savms;
   
    - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);     - hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k, nres);
      Computes the transition matrix starting at age 'age' over       Computes the transition matrix starting at age 'age' over
      'nhstepm*hstepm*stepm' months (i.e. until       'nhstepm*hstepm*stepm' months (i.e. until
      age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying       age (in years)  age+nhstepm*hstepm*stepm/12) by multiplying
Line 3174  double ***hpxij(double ***po, int nhstep Line 3177  double ***hpxij(double ***po, int nhstep
     }      }
     for(i=1; i<=nlstate+ndeath; i++)      for(i=1; i<=nlstate+ndeath; i++)
       for(j=1;j<=nlstate+ndeath;j++) {        for(j=1;j<=nlstate+ndeath;j++) {
                                 po[i][j][h]=newm[i][j];          po[i][j][h]=newm[i][j];
                                 /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/          /*if(h==nhstepm) printf("po[%d][%d][%d]=%f ",i,j,h,po[i][j][h]);*/
       }        }
     /*printf("h=%d ",h);*/      /*printf("h=%d ",h);*/
   } /* end h */    } /* end h */
         /*     printf("\n H=%d \n",h); */    /*     printf("\n H=%d \n",h); */
   return po;    return po;
 }  }
   
 /************* Higher Back Matrix Product ***************/  /************* Higher Back Matrix Product ***************/
 /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */  /* double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, double **oldm, double **savm, double **dnewm, double **doldm, double **dsavm, int ij ) */
 double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij )  double ***hbxij(double ***po, int nhstepm, double age, int hstepm, double *x, double ***prevacurrent, int nlstate, int stepm, int ij, int nres )
 {  {
   /* For a combination of dummy covariate ij, computes the transition matrix starting at age 'age' over    /* For a combination of dummy covariate ij, computes the transition matrix starting at age 'age' over
      'nhstepm*hstepm*stepm' months (i.e. until       'nhstepm*hstepm*stepm' months (i.e. until
Line 3230  double ***hbxij(double ***po, int nhstep Line 3233  double ***hbxij(double ***po, int nhstep
       /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */        /* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; *\/ */
         cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];          cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
         /* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */          /* printf("hbxij Dummy agexact=%.0f combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov[%d]=%lf codtabm(%d,Tvar[%d])=%d \n",agexact,ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],2+nagesqr+TvarsDind[k],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
   
       }        }
       for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */        for (k=1; k<=nsq;k++) { /* For single varying covariates only */
         /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */          /* Here comes the value of quantitative after renumbering k with single quantitative covariates */
         cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];          cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k]; 
       /* cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k]])]*cov[2]; */          /* printf("hPxij Quantitative k=%d  TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
       for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */        }
         for (k=1; k<=cptcovage;k++){ /* Should start at cptcovn+1 */
           if(Dummy[Tvar[Tage[k]]]){
             cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
           } else{
             cov[2+nagesqr+Tage[k]]=Tqresult[nres][k]; 
           }
           /* printf("hBxij Age combi=%d k=%d  Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
         }
         for (k=1; k<=cptcovprod;k++){ /* Useless because included in cptcovn */
         cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];          cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)]*nbcode[Tvard[k][2]][codtabm(ij,k)];
       /* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */        }                 
                           
       /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/        /*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
       /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/        /*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
   
       /* Careful transposed matrix */        /* Careful transposed matrix */
       /* age is in cov[2], prevacurrent at beginning of transition. */        /* age is in cov[2], prevacurrent at beginning of transition. */
       /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */        /* out=matprod2(newm, bmij(pmmij,cov,ncovmodel,x,nlstate,prevacurrent, dnewm, doldm, dsavm,ij),\ */
Line 7806  set ter svg size 640, 480\nunset log y\n Line 7817  set ter svg size 640, 480\nunset log y\n
     
   
 /************** Forecasting ******************/  /************** Forecasting ******************/
  void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){  void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
   /* proj1, year, month, day of starting projection     /* proj1, year, month, day of starting projection 
      agemin, agemax range of age       agemin, agemax range of age
      dateprev1 dateprev2 range of dates during which prevalence is computed       dateprev1 dateprev2 range of dates during which prevalence is computed
      anproj2 year of en of projection (same day and month as proj1).       anproj2 year of en of projection (same day and month as proj1).
   */    */
    int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;    int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
   double agec; /* generic age */    double agec; /* generic age */
   double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;    double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
   double *popeffectif,*popcount;    double *popeffectif,*popcount;
Line 7935  set ter svg size 640, 480\nunset log y\n Line 7946  set ter svg size 640, 480\nunset log y\n
 }  }
   
 /* /\************** Back Forecasting ******************\/ */  /* /\************** Back Forecasting ******************\/ */
 /* void prevbackforecast(char fileres[], double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){ */  void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
 /*   /\* back1, year, month, day of starting backection  */    /* back1, year, month, day of starting backection
 /*      agemin, agemax range of age */       agemin, agemax range of age
 /*      dateprev1 dateprev2 range of dates during which prevalence is computed */       dateprev1 dateprev2 range of dates during which prevalence is computed
 /*      anback2 year of en of backection (same day and month as back1). */       anback2 year of en of backection (same day and month as back1).
 /*   *\/ */    */
 /*   int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1; */    int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
 /*   double agec; /\* generic age *\/ */    double agec; /* generic age */
 /*   double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; */    double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
 /*   double *popeffectif,*popcount; */    double *popeffectif,*popcount;
 /*   double ***p3mat; */    double ***p3mat;
 /*   /\* double ***mobaverage; *\/ */    /* double ***mobaverage; */
 /*   char fileresfb[FILENAMELENGTH]; */    char fileresfb[FILENAMELENGTH];
            
 /*   agelim=AGESUP; */    agelim=AGESUP;
 /*   /\* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people */    /* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
 /*      in each health status at the date of interview (if between dateprev1 and dateprev2). */       in each health status at the date of interview (if between dateprev1 and dateprev2).
 /*      We still use firstpass and lastpass as another selection. */       We still use firstpass and lastpass as another selection.
 /*   *\/ */    */
 /*   /\* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ *\/ */    /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart,\ */
 /*   /\*              firstpass, lastpass,  stepm,  weightopt, model); *\/ */    /*          firstpass, lastpass,  stepm,  weightopt, model); */
 /*   prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */  
             /*Do we need to compute prevalence again?*/
 /*   strcpy(fileresfb,"FB_");  */  
 /*   strcat(fileresfb,fileresu); */    /* prevalence(probs, ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); */
 /*   if((ficresfb=fopen(fileresfb,"w"))==NULL) { */  
 /*     printf("Problem with back forecast resultfile: %s\n", fileresfb); */  
 /*     fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb); */  
 /*   } */  
 /*   printf("Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */  
 /*   fprintf(ficlog,"Computing back forecasting: result on file '%s', please wait... \n", fileresfb); */  
           
 /*   if (cptcoveff==0) ncodemax[cptcoveff]=1; */  
           
 /*   /\* if (mobilav!=0) { *\/ */  
 /*   /\*   mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */  
 /*   /\*   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; */  
 /*   if (stepm<=12) stepsize=1; */  
 /*   if(estepm < stepm){ */  
 /*     printf ("Problem %d lower than %d\n",estepm, stepm); */  
 /*   } */  
 /*   else  hstepm=estepm;    */  
           
 /*   hstepm=hstepm/stepm;  */  
 /*   yp1=modf(dateintmean,&yp);/\* extracts integral of datemean in yp  and */  
 /*                                fractional in yp1 *\/ */  
 /*   anprojmean=yp; */  
 /*   yp2=modf((yp1*12),&yp); */  
 /*   mprojmean=yp; */  
 /*   yp1=modf((yp2*30.5),&yp); */  
 /*   jprojmean=yp; */  
 /*   if(jprojmean==0) jprojmean=1; */  
 /*   if(mprojmean==0) jprojmean=1; */  
           
 /*   i1=cptcoveff; */  
 /*   if (cptcovn < 1){i1=1;} */  
       
 /*   fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);  */    strcpy(fileresfb,"FB_");
     strcat(fileresfb,fileresu);
     if((ficresfb=fopen(fileresfb,"w"))==NULL) {
       printf("Problem with back forecast resultfile: %s\n", fileresfb);
       fprintf(ficlog,"Problem with back forecast resultfile: %s\n", fileresfb);
     }
     printf("\nComputing back forecasting: result on file '%s', please wait... \n", fileresfb);
     fprintf(ficlog,"\nComputing back forecasting: result on file '%s', please wait... \n", fileresfb);
       
 /*   fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); */    if (cptcoveff==0) ncodemax[cptcoveff]=1;
             
 /*      /\*           if (h==(int)(YEARM*yearp)){ *\/ */     
 /*   for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */    stepsize=(int) (stepm+YEARM-1)/YEARM;
 /*     for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */    if (stepm<=12) stepsize=1;
 /*       k=k+1; */    if(estepm < stepm){
 /*       fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#"); */      printf ("Problem %d lower than %d\n",estepm, stepm);
 /*       for(j=1;j<=cptcoveff;j++) { */    }
 /*                              fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */    else  hstepm=estepm;
 /*       } */    
 /*       fprintf(ficresfb," yearbproj age"); */    hstepm=hstepm/stepm;
 /*       for(j=1; j<=nlstate+ndeath;j++){  */    yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and
 /*                              for(i=1; i<=nlstate;i++)               */                                 fractional in yp1 */
 /*           fprintf(ficresfb," p%d%d",i,j); */    anprojmean=yp;
 /*                              fprintf(ficresfb," p.%d",j); */    yp2=modf((yp1*12),&yp);
 /*       } */    mprojmean=yp;
 /*       for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {  */    yp1=modf((yp2*30.5),&yp);
 /*                              /\* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {  *\/ */    jprojmean=yp;
 /*                              fprintf(ficresfb,"\n"); */    if(jprojmean==0) jprojmean=1;
 /*                              fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);    */    if(mprojmean==0) jprojmean=1;
 /*                              for (agec=fage; agec>=(ageminpar-1); agec--){  */    
 /*                                      nhstepm=(int) rint((agelim-agec)*YEARM/stepm);  */    i1=pow(2,cptcoveff);
 /*                                      nhstepm = nhstepm/hstepm;  */    if (cptcovn < 1){i1=1;}
 /*                                      p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */    
 /*                                      oldm=oldms;savm=savms; */    fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);
 /*                                      hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm,oldm,savm, dnewm, doldm, dsavm, k);       */    
 /*                                      for (h=0; h<=nhstepm; h++){ */    fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
 /*                                              if (h*hstepm/YEARM*stepm ==yearp) { */    
 /*               fprintf(ficresfb,"\n"); */    /*          if (h==(int)(YEARM*yearp)){ */
 /*               for(j=1;j<=cptcoveff;j++)  */    /* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
 /*                 fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */    /*   for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
 /*                                                      fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm); */    /*     k=k+1; */
 /*                                              }  */    for(nres=1; nres <= nresult; nres++) /* For each resultline */
 /*                                              for(j=1; j<=nlstate+ndeath;j++) { */    for(k=1; k<=i1;k++){
 /*                                                      ppij=0.; */      if(i1 != 1 && TKresult[nres]!= k)
 /*                                                      for(i=1; i<=nlstate;i++) { */        continue;
 /*                                                              if (mobilav==1)  */      if(invalidvarcomb[k]){
 /*                                                                      ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod]; */        printf("\nCombination (%d) projection ignored because no cases \n",k); 
 /*                                                              else { */        continue;
 /*                                                                      ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod]; */      }
 /*                                                              } */      fprintf(ficresfb,"\n#****** hbijx=probability over h years, hp.jx is weighted by observed prev \n#");
 /*                                                              if (h*hstepm/YEARM*stepm== yearp) { */      for(j=1;j<=cptcoveff;j++) {
 /*                                                                      fprintf(ficresfb," %.3f", p3mat[i][j][h]); */        fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
 /*                                                              } */      }
 /*                                                      } /\* end i *\/ */      for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
 /*                                                      if (h*hstepm/YEARM*stepm==yearp) { */        fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
 /*                                                              fprintf(ficresfb," %.3f", ppij); */      }
 /*                                                      } */      fprintf(ficresfb," yearbproj age");
 /*                                              }/\* end j *\/ */      for(j=1; j<=nlstate+ndeath;j++){
 /*                                      } /\* end h *\/ */        for(i=1; i<=nlstate;i++)
 /*                                      free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); */          fprintf(ficresfb," p%d%d",i,j);
 /*                              } /\* end agec *\/ */        fprintf(ficresfb," p.%d",j);
 /*       } /\* end yearp *\/ */      }
 /*     } /\* end cptcod *\/ */      for (yearp=0; yearp>=(anback2-anback1);yearp -=stepsize) {
 /*   } /\* end  cptcov *\/ */        /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {  */
                 fprintf(ficresfb,"\n");
 /*   /\* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); *\/ */        fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
                 /* for (agec=fage; agec>=(ageminpar-1); agec--){ */
 /*   fclose(ficresfb); */        /*        nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */
 /*   printf("End of Computing Back forecasting \n"); */        for (agec=fage; agec>=fage-20; agec--){  /* testing up to 10 */
 /*   fprintf(ficlog,"End of Computing Back forecasting\n"); */          nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
           nhstepm = nhstepm/hstepm;
           p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           oldm=oldms;savm=savms;
           hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
   
           for (h=0; h<=nhstepm; h++){
             if (h*hstepm/YEARM*stepm ==yearp) {
               fprintf(ficresfb,"\n");
               for(j=1;j<=cptcoveff;j++)
                 fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
               fprintf(ficresfb,"%.f %.f ",anback1+yearp,agec+h*hstepm/YEARM*stepm);
             }
             for(j=1; j<=nlstate+ndeath;j++) {
               ppij=0.;
               for(i=1; i<=nlstate;i++) {
                 /* if (mobilav==1) */
                   ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];
                 /* else { */
                 /*        ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
                 /* } */
                 if (h*hstepm/YEARM*stepm== yearp) {
                   fprintf(ficresfb," %.3f", p3mat[i][j][h]);
                 }
               } /* end i */
               if (h*hstepm/YEARM*stepm==yearp) {
                 fprintf(ficresfb," %.3f", ppij);
               }
             }/* end j */
           } /* end h */
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         } /* end agec */
       } /* end yearp */
     } /* end k */
     
     /* if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); */
     
     fclose(ficresfb);
     printf("End of Computing Back forecasting \n");
     fprintf(ficlog,"End of Computing Back forecasting\n");
                   
 /* } */  }
   
 /************** Forecasting *****not tested NB*************/  /************** Forecasting *****not tested NB*************/
 /* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */  /* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */
Line 10071  int hPijx(double *p, int bage, int fage) Line 10091  int hPijx(double *p, int bage, int fage)
   
         /* oldm=oldms;savm=savms; */          /* oldm=oldms;savm=savms; */
         /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */          /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);   */
         hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k);          hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
         /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */          /* hbxij(p3mat,nhstepm,agedeb,hstepm,p,prevacurrent,nlstate,stepm,oldm,savm, dnewm, doldm, dsavm, k); */
         fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");          fprintf(ficrespijb,"# Cov Agex agex-h hbijx with i,j=");
         for(i=1; i<=nlstate;i++)          for(i=1; i<=nlstate;i++)
Line 11732  Please run with mle=-1 to get a correct Line 11752  Please run with mle=-1 to get a correct
       fclose(ficrespijb);        fclose(ficrespijb);
       free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */        free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
   
       /* prevbackforecast(fileresu, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,        prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,
          bage, fage, firstpass, lastpass, anback2, p, cptcoveff); */           bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
       free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);        free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
       free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);        free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
       free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);        free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);

Removed from v.1.266  
changed lines
  Added in v.1.267


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