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

Diff for /imach/src/imach.c between versions 1.201 and 1.202

version 1.201, 2015/09/15 17:34:58 version 1.202, 2015/09/22 19:45:16
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.202  2015/09/22 19:45:16  brouard
     Summary: Adding some overall graph on contribution to likelihood. Might change
   
   Revision 1.201  2015/09/15 17:34:58  brouard    Revision 1.201  2015/09/15 17:34:58  brouard
   Summary: 0.98r0    Summary: 0.98r0
   
Line 644 Line 647
   
 /* #define DEBUG */  /* #define DEBUG */
 /* #define DEBUGBRENT */  /* #define DEBUGBRENT */
   #define DEBUGLINMIN
 #define POWELL /* Instead of NLOPT */  #define POWELL /* Instead of NLOPT */
 #define POWELLF1F3 /* Skip test */  #define POWELLF1F3 /* Skip test */
 /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */  /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
Line 802  char command[FILENAMELENGTH]; Line 806  char command[FILENAMELENGTH];
 int  outcmd=0;  int  outcmd=0;
   
 char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];  char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
 char fileresu[FILENAMELENGTH]; /* Without r in front */  char fileresu[FILENAMELENGTH]; /* fileres without r in front */
 char filelog[FILENAMELENGTH]; /* Log file */  char filelog[FILENAMELENGTH]; /* Log file */
 char filerest[FILENAMELENGTH];  char filerest[FILENAMELENGTH];
 char fileregp[FILENAMELENGTH];  char fileregp[FILENAMELENGTH];
Line 1606  void linmin(double p[], double xi[], int Line 1610  void linmin(double p[], double xi[], int
   nrfunc=func;     nrfunc=func; 
   for (j=1;j<=n;j++) {     for (j=1;j<=n;j++) { 
     pcom[j]=p[j];       pcom[j]=p[j]; 
     xicom[j]=xi[j];       xicom[j]=xi[j]; /* Former scale xi[j] of currrent direction i */
   }     } 
   
   /* axs=0.0; */    /* axs=0.0; */
Line 1630  void linmin(double p[], double xi[], int Line 1634  void linmin(double p[], double xi[], int
   
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
   printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);    printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);
     fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);
 #endif  #endif
   *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/    *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/
   /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */    /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */
Line 1642  void linmin(double p[], double xi[], int Line 1647  void linmin(double p[], double xi[], int
 #endif  #endif
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
   printf("linmin end ");    printf("linmin end ");
     fprintf(ficlog,"linmin end ");
 #endif  #endif
   for (j=1;j<=n;j++) {     for (j=1;j<=n;j++) { 
     /* printf(" before xi[%d]=%12.8f", j,xi[j]); */      /* printf(" before xi[%d]=%12.8f", j,xi[j]); */
Line 1653  void linmin(double p[], double xi[], int Line 1659  void linmin(double p[], double xi[], int
   /* printf("\n"); */    /* printf("\n"); */
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
   printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));    printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
     fprintf(ficlog,"Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p));
   for (j=1;j<=n;j++) {     for (j=1;j<=n;j++) { 
     printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]);      printf(" xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]);
     if(j % ncovmodel == 0)      fprintf(ficlog," xi[%d]= %14.10f p[%d]= %12.7f",j,xi[j],j,p[j]);
       if(j % ncovmodel == 0){
       printf("\n");        printf("\n");
         fprintf(ficlog,"\n");
       }
   }    }
 #endif  #endif
   free_vector(xicom,1,n);     free_vector(xicom,1,n); 
Line 1691  void powell(double p[], double **xi, int Line 1701  void powell(double p[], double **xi, int
   xits=vector(1,n);     xits=vector(1,n); 
   *fret=(*func)(p);     *fret=(*func)(p); 
   for (j=1;j<=n;j++) pt[j]=p[j];     for (j=1;j<=n;j++) pt[j]=p[j]; 
     rcurr_time = time(NULL);      rcurr_time = time(NULL);  
   for (*iter=1;;++(*iter)) {     for (*iter=1;;++(*iter)) { 
     fp=(*fret); /* From former iteration or initial value */      fp=(*fret); /* From former iteration or initial value */
     ibig=0;       ibig=0; 
Line 1836  void powell(double p[], double **xi, int Line 1846  void powell(double p[], double **xi, int
       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */
       t= t- del*SQR(fp-fptt);        t= t- del*SQR(fp-fptt);
 #endif  #endif
       directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */        directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */
 #ifdef DEBUG  #ifdef DEBUG
       printf("t1= %.12lf, t2= %.12lf, t=%.12lf  directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);        printf("t1= %.12lf, t2= %.12lf, t=%.12lf  directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
       fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);        fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
Line 1851  void powell(double p[], double **xi, int Line 1861  void powell(double p[], double **xi, int
       if (t < 0.0) { /* Then we use it for new direction */        if (t < 0.0) { /* Then we use it for new direction */
 #else  #else
       if (directest*t < 0.0) { /* Contradiction between both tests */        if (directest*t < 0.0) { /* Contradiction between both tests */
         printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);          printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
         printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);          printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
         fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);          fprintf(ficlog,"directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
         fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);          fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
       }         } 
       if (directest < 0.0) { /* Then we use it for new direction */        if (directest < 0.0) { /* Then we use it for new direction */
Line 1861  void powell(double p[], double **xi, int Line 1871  void powell(double p[], double **xi, int
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
         printf("Before linmin in direction P%d-P0\n",n);          printf("Before linmin in direction P%d-P0\n",n);
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
           printf("Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);            printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
           if(j % ncovmodel == 0)            fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
             if(j % ncovmodel == 0){
             printf("\n");              printf("\n");
               fprintf(ficlog,"\n");
             }
         }          }
 #endif  #endif
         linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/          linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
           printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);            printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
           if(j % ncovmodel == 0)            fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
             if(j % ncovmodel == 0){
             printf("\n");              printf("\n");
               fprintf(ficlog,"\n");
             }
         }          }
 #endif  #endif
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
Line 1911  double **prevalim(double **prlim, int nl Line 1927  double **prevalim(double **prlim, int nl
   /* double **matprod2(); */ /* test */    /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **pmij();    double **out, cov[NCOVMAX+1], **pmij();
   double **newm;    double **newm;
   double agefin, delaymax=50 ; /* Max number of years to converge */    double agefin, delaymax=100 ; /* Max number of years to converge */
     long int ncvyear=0, ncvloop=0;
       
   for (ii=1;ii<=nlstate+ndeath;ii++)    for (ii=1;ii<=nlstate+ndeath;ii++)
     for (j=1;j<=nlstate+ndeath;j++){      for (j=1;j<=nlstate+ndeath;j++){
Line 1921  double **prevalim(double **prlim, int nl Line 1938  double **prevalim(double **prlim, int nl
   cov[1]=1.;    cov[1]=1.;
       
   /* Even if hstepm = 1, at least one multiplication by the unit matrix */    /* Even if hstepm = 1, at least one multiplication by the unit matrix */
     /* Start at agefin= age, computes the matrix of passage and loops decreasing agefin until convergence is reached */
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){    for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
       ncvloop++;
     newm=savm;      newm=savm;
     /* Covariates have to be included here again */      /* Covariates have to be included here again */
     cov[2]=agefin;      cov[2]=agefin;
Line 1956  double **prevalim(double **prlim, int nl Line 1975  double **prevalim(double **prlim, int nl
         sumnew=0;          sumnew=0;
         for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];          for(k=1; k<=ndeath; k++) sumnew+=newm[i][nlstate+k];
         prlim[i][j]= newm[i][j]/(1-sumnew);          prlim[i][j]= newm[i][j]/(1-sumnew);
         /*printf(" prevalim i=%d, j=%d, prmlim[%d][%d]=%f, agefin=%d \n", i, j, i, j, prlim[i][j],(int)agefin);*/  
         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); */
       }        }
       maxmin=max-min;        maxmin=max-min;
       maxmax=FMAX(maxmax,maxmin);        maxmax=FMAX(maxmax,maxmin);
     } /* j loop */      } /* j loop */
     if(maxmax < ftolpl){      if(maxmax < ftolpl){
         /* printf("maxmax=%lf maxmin=%lf ncvloop=%ld, ncvyear=%d \n", maxmax, maxmin, ncvloop, (int)age-(int)agefin); */
       return prlim;        return prlim;
     }      }
   } /* age loop */    } /* age loop */
     printf("Warning: the stable prevalence did not converge with the required precision ftolpl=6*10^5*ftol=%g. \n\
   Earliest age to start was %d-%d=%d, ncvloop=%ld, ncvyear=%d\n\
   Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin);
   return prlim; /* should not reach here */    return prlim; /* should not reach here */
 }  }
   
Line 2542  double funcone( double *x) Line 2565  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 %6d %2d %2d %1d %1d %3d %11.6f %8.4f\          fprintf(ficresilk,"%9ld %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f\
  %11.6f %11.6f %11.6f ", \   %11.6f %11.6f %11.6f ", \
                 num[i],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],
                 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 2577  void likelione(FILE *ficres,double p[], Line 2600  void likelione(FILE *ficres,double p[],
   
   if(*globpri !=0){ /* Just counts and sums, no printings */    if(*globpri !=0){ /* Just counts and sums, no printings */
     strcpy(fileresilk,"ILK_");       strcpy(fileresilk,"ILK_"); 
     strcat(fileresilk,fileres);      strcat(fileresilk,fileresu);
     if((ficresilk=fopen(fileresilk,"w"))==NULL) {      if((ficresilk=fopen(fileresilk,"w"))==NULL) {
       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) 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, "#num_i i s1 s2 mi mw dh likeli weight 2wlli out sav ");      fprintf(ficresilk, "#num_i age i s1 s2 mi mw dh likeli 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 2593  void likelione(FILE *ficres,double p[], Line 2616  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 (if mle=1): <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));      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));
     fflush(fichtm);       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> \
   <img src=\"%s.png\">",subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"),subdirf2(optionfilefiname,"ILK_"));
       fflush(fichtm);
   }     } 
   return;    return;
 }  }
Line 3841  void varevsij(char optionfilefiname[], d Line 3866  void varevsij(char optionfilefiname[], d
   /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/    /*printf("DIGIT=%s, ij=%d ijr=%-d|\n",digit, ij,ij);*/
   strcat(fileresprobmorprev,digit); /* Tvar to be done */    strcat(fileresprobmorprev,digit); /* Tvar to be done */
   strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */    strcat(fileresprobmorprev,digitp); /* Popbased or not, mobilav or not */
   strcat(fileresprobmorprev,fileres);    strcat(fileresprobmorprev,fileresu);
   if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {    if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprobmorprev);      printf("Problem with resultfile: %s\n", fileresprobmorprev);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
Line 4213  void varprob(char optionfilefiname[], do Line 4238  void varprob(char optionfilefiname[], do
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
   }    }
   strcpy(fileresprobcov,"PROBCOV_");     strcpy(fileresprobcov,"PROBCOV_"); 
   strcat(fileresprobcov,fileres);    strcat(fileresprobcov,fileresu);
   if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {    if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprobcov);      printf("Problem with resultfile: %s\n", fileresprobcov);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
   }    }
   strcpy(fileresprobcor,"PROBCOR_");     strcpy(fileresprobcor,"PROBCOR_"); 
   strcat(fileresprobcor,fileres);    strcat(fileresprobcor,fileresu);
   if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {    if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
     printf("Problem with resultfile: %s\n", fileresprobcor);      printf("Problem with resultfile: %s\n", fileresprobcor);
     fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
Line 4571  fprintf(fichtm," \n<ul><li><b>Graphs</b> Line 4596  fprintf(fichtm," \n<ul><li><b>Graphs</b>
      fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \       fprintf(fichtm,"<br>- Logit model, for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \
 <img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);  <img src=\"%s_%d-1.svg\">",subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);
      /* Pij */       /* Pij */
      fprintf(fichtm,"<br>\n- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \       fprintf(fichtm,"<br>\n- Pij or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \
 <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);       <img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);     
      /* Quasi-incidences */       /* Quasi-incidences */
      fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\       fprintf(fichtm,"<br>\n- Iij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\
Line 4698  void printinggnuplot(char fileresu[], ch Line 4723  void printinggnuplot(char fileresu[], ch
     /*#endif */      /*#endif */
   m=pow(2,cptcoveff);    m=pow(2,cptcoveff);
   
     /* Contribution to likelihood */
     /* Plot the probability implied in the likelihood */
       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,"\nset ter svg size 640, 480"); */ /* Too big for svg */
       fprintf(ficgp,"\nset ter png size 640, 480");
   /* good 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 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.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,"\nreplot  \"%s\" u 2:($3 <= 3 ? -$11 : 1/0):3 t \"First 3 individuals\" with line lc variable", subdirf(fileresilk));
       fprintf(ficgp,"\nset out\n");
       /* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
   
   strcpy(dirfileres,optionfilefiname);    strcpy(dirfileres,optionfilefiname);
   strcpy(optfileres,"vpl");    strcpy(optfileres,"vpl");
  /* 1eme*/   /* 1eme*/
Line 4816  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 4857  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   } /* end covariate */      } /* end covariate */  
   
   /* Survival functions (period) from state i in state j by final state j */    /* Survival functions (period) from state i in state j by final state j */
   for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */    for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
       k=3;        k=3;
       fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);        fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
Line 4849  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 4890  plot [%.f:%.f]  ", ageminpar, agemaxpar)
     } /* end cpt state*/       } /* end cpt state*/ 
   } /* end covariate */      } /* end covariate */  
   
   /* CV preval stable (period) */    /* CV preval stable (period) for each covariate */
   for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */    for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
       k=3;        k=3;
       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);        fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt);
Line 6276  void syscompilerinfo(int logged) Line 6317  void syscompilerinfo(int logged)
   
  }   }
   
 int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){   int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar, double ftolpl){
   /*--------------- Prevalence limit  (period or stable prevalence) --------------*/    /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
   int i, j, k, i1 ;    int i, j, k, i1 ;
   double ftolpl = 1.e-10;    /* double ftolpl = 1.e-10; */
   double age, agebase, agelim;    double age, agebase, agelim;
   
     strcpy(filerespl,"PL_");    strcpy(filerespl,"PL_");
     strcat(filerespl,fileresu);    strcat(filerespl,fileresu);
     if((ficrespl=fopen(filerespl,"w"))==NULL) {    if((ficrespl=fopen(filerespl,"w"))==NULL) {
       printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;      printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
       fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;      fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
     }    }
     printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);    printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
     fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);    fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
     pstamp(ficrespl);    pstamp(ficrespl);
     fprintf(ficrespl,"# Period (stable) prevalence \n");    fprintf(ficrespl,"# Period (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");
       
     /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */      /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
   
Line 6684  int main(int argc, char *argv[]) Line 6725  int main(int argc, char *argv[])
     }      }
     printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);      printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt);
   }    }
     ftolpl=6*ftol*1.e5; /* 6.e-3 make convergences in less than 80 loops for the prevalence limit */
   /* Third parameter line */    /* Third parameter line */
   while(fgets(line, MAXLINE, ficpar)) {    while(fgets(line, MAXLINE, ficpar)) {
     /* If line starts with a # it is a comment */      /* If line starts with a # it is a comment */
Line 7089  Please run with mle=-1 to get a correct Line 7130  Please run with mle=-1 to get a correct
            *    15 i=8 1     2     2     2             *    15 i=8 1     2     2     2
            *    16     2     2     2     2             *    16     2     2     2     2
            */             */
   for(h=1; h <=100 ;h++){     /* /\* for(h=1; h <=100 ;h++){  *\/ */
     /* printf("h=%2d ", h); */    /*   /\* printf("h=%2d ", h); *\/ */
      /* for(k=1; k <=10; k++){ */    /*    /\* for(k=1; k <=10; k++){ *\/ */
        /* printf("k=%d %d ",k,codtabm(h,k)); */    /*      /\* printf("k=%d %d ",k,codtabm(h,k)); *\/ */
      /*   codtab[h][k]=codtabm(h,k); */    /*    /\*   codtab[h][k]=codtabm(h,k); *\/ */
      /* } */    /*    /\* } *\/ */
      /* printf("\n"); */    /*    /\* printf("\n"); *\/ */
   }    /* } */
   /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */    /* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */
   /*   for(i=1; i <=pow(2,cptcoveff-k);i++){ /\* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 *\/  */    /*   for(i=1; i <=pow(2,cptcoveff-k);i++){ /\* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 *\/  */
   /*     for(j=1; j <= ncodemax[k]; j++){ /\* For each modality of this covariate ncodemax=2*\/ */    /*     for(j=1; j <= ncodemax[k]; j++){ /\* For each modality of this covariate ncodemax=2*\/ */
Line 7717  Please run with mle=-1 to get a correct Line 7758  Please run with mle=-1 to get a correct
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/      /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
     /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */      /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */
     prlim=matrix(1,nlstate,1,nlstate);      prlim=matrix(1,nlstate,1,nlstate);
     prevalence_limit(p, prlim,  ageminpar, agemaxpar);      prevalence_limit(p, prlim,  ageminpar, agemaxpar, ftolpl);
     fclose(ficrespl);      fclose(ficrespl);
   
 #ifdef FREEEXIT2  #ifdef FREEEXIT2
Removed from v.1.201  
changed lines
  Added in v.1.202

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