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

Diff for /imach/src/imach.c between versions 1.203 and 1.207

version 1.203, 2015/09/30 17:45:14 version 1.207, 2015/10/27 17:36:57
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.207  2015/10/27 17:36:57  brouard
     *** empty log message ***
   
     Revision 1.206  2015/10/24 07:14:11  brouard
     *** empty log message ***
   
     Revision 1.205  2015/10/23 15:50:53  brouard
     Summary: 0.98r3 some clarification for graphs on likelihood contributions
   
     Revision 1.204  2015/10/01 16:20:26  brouard
     Summary: Some new graphs of contribution to likelihood
   
   Revision 1.203  2015/09/30 17:45:14  brouard    Revision 1.203  2015/09/30 17:45:14  brouard
   Summary: looking at better estimation of the hessian    Summary: looking at better estimation of the hessian
   
Line 749  typedef struct { Line 761  typedef struct {
 /* $State$ */  /* $State$ */
 #include "version.h"  #include "version.h"
 char version[]=__IMACH_VERSION__;  char version[]=__IMACH_VERSION__;
 char copyright[]="September 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";  char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 940  static int split( char *path, char *dirc Line 952  static int split( char *path, char *dirc
     }      }
     /* got dirc from getcwd*/      /* got dirc from getcwd*/
     printf(" DIRC = %s \n",dirc);      printf(" DIRC = %s \n",dirc);
   } else {                              /* strip direcotry from path */    } else {                              /* strip directory from path */
     ss++;                               /* after this, the filename */      ss++;                               /* after this, the filename */
     l2 = strlen( ss );                  /* length of filename */      l2 = strlen( ss );                  /* length of filename */
     if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );      if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
Line 1954  double **prevalim(double **prlim, int nl Line 1966  double **prevalim(double **prlim, int nl
 {  {
   /* Computes the prevalence limit in each live state at age x by left multiplying the unit    /* Computes the prevalence limit in each live state at age x by left multiplying the unit
      matrix by transitions matrix until convergence is reached with precision ftolpl */       matrix by transitions matrix until convergence is reached with precision ftolpl */
       /* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1  = Wx-n Px-n ... Px-2 Px-1 I */
     /* Wx is row vector: population in state 1, population in state 2, population dead */
     /* or prevalence in state 1, prevalence in state 2, 0 */
     /* newm is the matrix after multiplications, its rows are identical at a factor */
     /* Initial matrix pimij */
     /* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */
     /* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */
     /*  0,                   0                  , 1} */
     /*
      * and after some iteration: */
     /* {0.45504275246439968, 0.42731458730878791, 0.11764266022681241, */
     /*  0.45201005341706885, 0.42865420071559901, 0.11933574586733192, */
     /*  0,                   0                  , 1} */
     /* And prevalence by suppressing the deaths are close to identical rows in prlim: */
     /* {0.51571254859325999, 0.4842874514067399, */
     /*  0.51326036147820708, 0.48673963852179264} */
     /* If we start from prlim again, prlim tends to a constant matrix */
   
   int i, ii,j,k;    int i, ii,j,k;
   double min, max, maxmin, maxmax,sumnew=0.;    double min, max, maxmin, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    /* double **matprod2(); */ /* test */
Line 2010  double **prevalim(double **prlim, int nl Line 2039  double **prevalim(double **prlim, int nl
         prlim[i][j]= newm[i][j]/(1-sumnew);          prlim[i][j]= newm[i][j]/(1-sumnew);
         max=FMAX(max,prlim[i][j]);          max=FMAX(max,prlim[i][j]);
         min=FMIN(min,prlim[i][j]);          min=FMIN(min,prlim[i][j]);
         /* printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min); */          printf(" age= %d prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d max=%f min=%f\n", (int)age, i, j, i, j, prlim[i][j],(int)agefin, max, min);
       }        }
       maxmin=(max-min)/(max+min)*2;        maxmin=(max-min)/(max+min)*2;
       maxmax=FMAX(maxmax,maxmin);        maxmax=FMAX(maxmax,maxmin);
     } /* j loop */      } /* j loop */
     *ncvyear= (int)age- (int)agefin;      *ncvyear= (int)age- (int)agefin;
     /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */      printf("maxmax=%lf maxmin=%lf ncvloop=%d, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear);
     if(maxmax < ftolpl){      if(maxmax < ftolpl){
       /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */        /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, age=%d, agefin=%d ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age, (int)agefin, *ncvyear); */
       return prlim;        return prlim;
Line 2600  double funcone( double *x) Line 2629  double funcone( double *x)
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;        ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */        /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       if(globpr){        if(globpr){
         fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f\          fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \   %11.6f %11.6f %11.6f ", \
                 num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],                  num[i], agexact, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
                 2*weight[i]*lli,out[s1][s2],savm[s1][s2]);                  2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
         for(k=1,llt=0.,l=0.; k<=nlstate; k++){          for(k=1,llt=0.,l=0.; k<=nlstate; k++){
           llt +=ll[k]*gipmx/gsw;            llt +=ll[k]*gipmx/gsw;
Line 2640  void likelione(FILE *ficres,double p[], Line 2669  void likelione(FILE *ficres,double p[],
       printf("Problem with resultfile: %s\n", fileresilk);        printf("Problem with resultfile: %s\n", fileresilk);
       fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);        fprintf(ficlog,"Problem with resultfile: %s\n", fileresilk);
     }      }
     fprintf(ficresilk, "#individual(line's_record) s1 s2 wave# effective_wave# number_of_matrices_product pij weight -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");      fprintf(ficresilk, "#individual(line's_record) count age s1 s2 wave# effective_wave# number_of_matrices_product pij weight weight/gpw -2ln(pij)*weight 0pij_x 0pij_(x-stepm) cumulating_loglikeli_by_health_state(reweighted=-2ll*weightXnumber_of_contribs/sum_of_weights) and_total\n");
     fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight 2wlli out sav ");      fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli weight %%weight 2wlli out sav ");
     /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */      /*  i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],2*weight[i]*lli,out[s1][s2],savm[s1][s2]); */
     for(k=1; k<=nlstate; k++)       for(k=1; k<=nlstate; k++) 
       fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);        fprintf(ficresilk," -2*gipw/gsw*weight*ll[%d]++",k);
Line 2651  void likelione(FILE *ficres,double p[], Line 2680  void likelione(FILE *ficres,double p[],
   *fretone=(*funcone)(p);    *fretone=(*funcone)(p);
   if(*globpri !=0){    if(*globpri !=0){
     fclose(ficresilk);      fclose(ficresilk);
     fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with initial parameters and mle >= 1. You should at least run with mle >= 1 and starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));      if (mle ==0)
     fprintf(fichtm,"<br>- The first 3 individuals are drawn with lines. The function drawn is -2Log(L) in log scale: <a href=\"%s.png\">%s.png</a><br> \        fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with initial parameters and mle = %d.",mle);
 <img src=\"%s.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));      else if(mle >=1)
         fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);
       fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
       
       fprintf(fichtm,"<br>- The function drawn is -2Log(L) in Log scale: by state of origin <a href=\"%s-ori.png\">%s-ori.png</a><br> \
   <img src=\"%s-ori.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
       fprintf(fichtm,"<br>- and by state of destination <a href=\"%s-dest.png\">%s-dest.png</a><br> \
   <img src=\"%s-dest.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
     fflush(fichtm);      fflush(fichtm);
   }         
       for (k=1; k<= nlstate ; k++) {
         fprintf(fichtm,"<br>- Probability p%dj by origin %d and destination j <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
   <img src=\"%s-p%dj.png\">",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k);
       }
     }
   return;    return;
 }  }
   
Line 4264  void cvevsij(double ***eij, double x[], Line 4305  void cvevsij(double ***eij, double x[],
 /************ Variance of prevlim ******************/  /************ Variance of prevlim ******************/
  void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[])   void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyear, int ij, char strstart[])
 {  {
   /* Variance of prevalence limit */    /* Variance of prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/
   
   double **dnewm,**doldm;    double **dnewm,**doldm;
Line 4323  void cvevsij(double ***eij, double x[], Line 4364  void cvevsij(double ***eij, double x[],
   
     for(i=1;i<=nlstate;i++)      for(i=1;i<=nlstate;i++)
       varpl[i][(int)age] =0.;        varpl[i][(int)age] =0.;
       if((int)age==67 ||(int)age== 66 ){
       matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
       matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
       }else{
     matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);      matprod2(dnewm,trgradg,1,nlstate,1,npar,1,npar,matcov);
     matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);      matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg);
       }
     for(i=1;i<=nlstate;i++)      for(i=1;i<=nlstate;i++)
       varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */        varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */
   
Line 4755  divided by h: hPij/h : <a href=\"%s_%d-3 Line 4801  divided by h: hPij/h : <a href=\"%s_%d-3
 <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);  <img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);
      }       }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s%d%d.svg\">%s%d%d.svg</a> <br> \         fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \
 <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);  <img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);
      }       }
    /* } /\* end i1 *\/ */     /* } /\* end i1 *\/ */
Line 4825  See page 'Matrix of variance-covariance Line 4871  See page 'Matrix of variance-covariance
      }       }
      for(cpt=1; cpt<=nlstate;cpt++) {       for(cpt=1; cpt<=nlstate;cpt++) {
        fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \         fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \
 prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg <br>\  prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d%d.svg\"> %s_%d-%d.svg <br>\
 <img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);    <img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);  
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \       fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \  health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
 true period expectancies (those weighted with period prevalences are also\  true period expectancies (those weighted with period prevalences are also\
  drawn in addition to the population based expectancies computed using\   drawn in addition to the population based expectancies computed using\
  observed and cahotic prevalences: %s_%d.svg<br>\   observed and cahotic prevalences:  <a href=\"%s_%d.svg\">%s_%d.svg<br>\
 <img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);  <img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);
    /* } /\* end i1 *\/ */     /* } /\* end i1 *\/ */
  }/* End k1 */   }/* End k1 */
  fprintf(fichtm,"</ul>");   fprintf(fichtm,"</ul>");
Line 4862  void printinggnuplot(char fileresu[], ch Line 4908  void printinggnuplot(char fileresu[], ch
     fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");      fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
     fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");      fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Likelihood (-2Log(L))\";");
     /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */      /* fprintf(ficgp,"\nset ter svg size 640, 480"); */ /* Too big for svg */
     fprintf(ficgp,"\nset ter png size 640, 480");      fprintf(ficgp,"\nset ter pngcairo size 640, 480");
 /* good for mle=4 plot by number of matrix products.  /* nice for mle=4 plot by number of matrix products.
    replot  "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */     replot  "rrtest1/toto.txt" u 2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with point lc 1 */
 /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)"  */  /* replot exp(p1+p2*x)/(1+exp(p1+p2*x)+exp(p3+p4*x)+exp(p5+p6*x)) t "p12(x)"  */
     /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */      /* fprintf(ficgp,"\nset out \"%s.svg\";",subdirf2(optionfilefiname,"ILK_")); */
     fprintf(ficgp,"\nset out \"%s.png\";",subdirf2(optionfilefiname,"ILK_"));      fprintf(ficgp,"\nset out \"%s-dest.png\";",subdirf2(optionfilefiname,"ILK_"));
     fprintf(ficgp,"\nplot  \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk));      fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$12):5 t \"All sample, transitions colored by destination\" with dots lc variable; set out;\n",subdirf(fileresilk));
     fprintf(ficgp,"\nreplot  \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk));      fprintf(ficgp,"\nset out \"%s-ori.png\";",subdirf2(optionfilefiname,"ILK_"));
       fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$12):4 t \"All sample, transitions colored by origin\" with dots lc variable; set out;\n\n",subdirf(fileresilk));
       for (i=1; i<= nlstate ; i ++) {
         fprintf(ficgp,"\nset out \"%s-p%dj.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i);
         fprintf(ficgp,"unset log;\n# plot weighted, mean weight should have point size of 0.5\n plot  \"%s\"",subdirf(fileresilk));
         fprintf(ficgp,"  u  2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable \\\n",i,1,i,1);
         for (j=2; j<= nlstate+ndeath ; j ++) {
           fprintf(ficgp,",\\\n \"\" u  2:($4 == %d && $5==%d ? $9 : 1/0):($11/4.):5 t \"p%d%d\" with points pointtype 7 ps variable lc variable ",i,j,i,j);
         }
         fprintf(ficgp,";\nset out; unset ylabel;\n"); 
       }
       /* unset log; plot  "rrtest1_sorted_4/ILK_rrtest1_sorted_4.txt" u  2:($4 == 1 && $5==2 ? $9 : 1/0):5 t "p12" with points lc variable */              
       /* fprintf(ficgp,"\nset log y;plot  \"%s\" u 2:(-$11):3 t \"All sample, all transitions\" with dots lc variable",subdirf(fileresilk)); */
       /* fprintf(ficgp,"\nreplot  \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk)); */
     fprintf(ficgp,"\nset out;unset log\n");      fprintf(ficgp,"\nset out;unset log\n");
     /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */      /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
   
Line 6718  int main(int argc, char *argv[]) Line 6777  int main(int argc, char *argv[])
   printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion);    printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion);
   if(argc <=1){    if(argc <=1){
     printf("\nEnter the parameter file name: ");      printf("\nEnter the parameter file name: ");
     fgets(pathr,FILENAMELENGTH,stdin);      if(!fgets(pathr,FILENAMELENGTH,stdin)){
         printf("ERROR Empty parameter file name\n");
         goto end;
       }
     i=strlen(pathr);      i=strlen(pathr);
     if(pathr[i-1]=='\n')      if(pathr[i-1]=='\n')
       pathr[i-1]='\0';        pathr[i-1]='\0';
     i=strlen(pathr);      i=strlen(pathr);
     if(pathr[i-1]==' ') /* This may happen when dragging on oS/X! */      if(i >= 1 && pathr[i-1]==' ') {/* This may happen when dragging on oS/X! */
       pathr[i-1]='\0';        pathr[i-1]='\0';
    for (tok = pathr; tok != NULL; ){      }
       i=strlen(pathr);
       if( i==0 ){
         printf("ERROR Empty parameter file name\n");
         goto end;
       }
       for (tok = pathr; tok != NULL; ){
       printf("Pathr |%s|\n",pathr);        printf("Pathr |%s|\n",pathr);
       while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');        while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
       printf("val= |%s| pathr=%s\n",val,pathr);        printf("val= |%s| pathr=%s\n",val,pathr);
Line 7319  Please run with mle=-1 to get a correct Line 7387  Please run with mle=-1 to get a correct
     printf("Problem with file %s",optionfilegnuplot);      printf("Problem with file %s",optionfilegnuplot);
   }    }
   else{    else{
     fprintf(ficgp,"\n# %s\n", version);       fprintf(ficgp,"\n# IMaCh-%s\n", version); 
     fprintf(ficgp,"# %s\n", optionfilegnuplot);       fprintf(ficgp,"# %s\n", optionfilegnuplot); 
     //fprintf(ficgp,"set missing 'NaNq'\n");      //fprintf(ficgp,"set missing 'NaNq'\n");
     fprintf(ficgp,"set datafile missing 'NaNq'\n");      fprintf(ficgp,"set datafile missing 'NaNq'\n");
Line 7346  Please run with mle=-1 to get a correct Line 7414  Please run with mle=-1 to get a correct
   else{    else{
   fprintf(fichtmcov,"<html><head>\n<title>IMaCh Cov %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \    fprintf(fichtmcov,"<html><head>\n<title>IMaCh Cov %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n",\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
           optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);            optionfilehtmcov,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
   }    }
   
   fprintf(fichtm,"<html><head>\n<title>IMaCh %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \    fprintf(fichtm,"<html><head>\n<head>\n<meta charset=\"utf-8\"/><meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n<title>IMaCh %s</title></head>\n <body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n<font size=\"3\">Sponsored by Copyright (C)  2002-2015 <a href=http://www.ined.fr>INED</a>-EUROREVES-Institut de longévité-Japan Society for the Promotion of Sciences 日本学術振興会 (<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - <a href=https://software.intel.com/en-us>Intel Software 2015</a></font><br>  \
   <hr size=\"2\" color=\"#EC5E5E\"> \n\
   <font size=\"2\">IMaCh-%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n\
 \n\  \n\
 <hr  size=\"2\" color=\"#EC5E5E\">\  <hr  size=\"2\" color=\"#EC5E5E\">\
  <ul><li><h4>Parameter files</h4>\n\   <ul><li><h4>Parameter files</h4>\n\
Line 7619  Please run with mle=-1 to get a correct Line 7689  Please run with mle=-1 to get a correct
     free_matrix(ximort,1,NDIM,1,NDIM);      free_matrix(ximort,1,NDIM,1,NDIM);
 #endif  #endif
   } /* Endof if mle==-3 mortality only */    } /* Endof if mle==-3 mortality only */
   /* Standard maximisation */    /* Standard  */
   else{ /* For mle !=- 3 */    else{ /* For mle !=- 3, could be 0 or 1 or 4 etc. */
     globpr=0;/* debug */      globpr=0;/* Computes sum of likelihood for globpr=1 and funcone */
     /* Computes likelihood for initial parameters */      /* Computes likelihood for initial parameters, uses funcone to compute gpimx and gsw */
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */      likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
     printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);      printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
     for (k=1; k<=npar;k++)      for (k=1; k<=npar;k++)
       printf(" %d %8.5f",k,p[k]);        printf(" %d %8.5f",k,p[k]);
     printf("\n");      printf("\n");
     globpr=1; /* again, to print the contributions */      if(mle>=1){ /* Could be 1 or 2, Real Maximization */
         /* mlikeli uses func not funcone */
         mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
       }
       if(mle==0) {/* No optimization, will print the likelihoods for the datafile */
         globpr=0;/* Computes sum of likelihood for globpr=1 and funcone */
         /* Computes likelihood for initial parameters, uses funcone to compute gpimx and gsw */
         likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
       }
       globpr=1; /* again, to print the individual contributions using computed gpimx and gsw */
     likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */      likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
     printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);      printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
     for (k=1; k<=npar;k++)      for (k=1; k<=npar;k++)
       printf(" %d %8.5f",k,p[k]);        printf(" %d %8.5f",k,p[k]);
     printf("\n");      printf("\n");
     if(mle>=1){ /* Could be 1 or 2, Real Maximisation */  
       mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);  
     }  
           
     /*--------- results files --------------*/      /*--------- results files --------------*/
     fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);      fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model);
Line 8089  Please run with mle=-1 to get a correct Line 8165  Please run with mle=-1 to get a correct
             }              }
                   
             fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);              fprintf(ficrest," %4.0f %d %d",age, vpopbased, mobilav);
               /* fprintf(ficrest," %4.0f %d %d %d %d",age, vpopbased, mobilav,Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); */ /* to be done */
             /* printf(" age %4.0f ",age); */              /* printf(" age %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++) {
Removed from v.1.203  
changed lines
  Added in v.1.207

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