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

Diff for /imach/src/imach.c between versions 1.359 and 1.364

version 1.359, 2024/04/24 21:21:17 version 1.364, 2024/06/28 12:27:05
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.364  2024/06/28 12:27:05  brouard
     * imach.c (Module): fixing some bugs in gnuplot and quantitative variables, but not completely solved
   
     Revision 1.363  2024/06/28 09:31:55  brouard
     Summary: Adding log lines too
   
     Revision 1.362  2024/06/28 08:00:31  brouard
     Summary: 0.99s6
   
     * imach.c (Module): s6 errors with age*age (harmless).
   
     Revision 1.361  2024/05/12 20:29:32  brouard
     Summary: Version 0.99s5
   
     * src/imach.c Version 0.99s5 In fact, the covariance of total life
     expectancy e.. with a partial life expectancy e.j is high,
     therefore the complete matrix of variance covariance has to be
     included in the formula of the standard error of the proportion of
     total life expectancy spent in a specific state:
     var(X/Y)=mu_x^2/mu_y^2*(sigma_x^2/mu_x^2 -2
     sigma_xy/mu_x/mu_y+sigma^2/mu_y^2).  Also an error with mle=-3
     made the program core dump. It is fixed in this version.
   
     Revision 1.360  2024/04/30 10:59:22  brouard
     Summary: Version 0.99s4 and estimation of std of e.j/e..
   
   Revision 1.359  2024/04/24 21:21:17  brouard    Revision 1.359  2024/04/24 21:21:17  brouard
   Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes    Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes
   
Line 1225  Important routines Line 1251  Important routines
 - Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities)  - Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities)
   and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually.    and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually.
 - printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables  - printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables
   o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, éliminating 1 1 if    o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, eliminating 1 1 if
     race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless.      race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless.
   
   
Line 1394  double gnuplotversion=GNUPLOTVERSION; Line 1420  double gnuplotversion=GNUPLOTVERSION;
 /* $State$ */  /* $State$ */
 #include "version.h"  #include "version.h"
 char version[]=__IMACH_VERSION__;  char version[]=__IMACH_VERSION__;
 char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";  char copyright[]="April 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024";
 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 1516  char *endptr; Line 1542  char *endptr;
 long lval;  long lval;
 double dval;  double dval;
   
   /* This for praxis gegen */
     /* int prin=1; */
     double h0=0.25;
     double macheps;
     double ffmin;
   
 #define NR_END 1  #define NR_END 1
 #define FREE_ARG char*  #define FREE_ARG char*
 #define FTOL 1.0e-10  #define FTOL 1.0e-10
Line 2992  static void print2() /* print a line of Line 3024  static void print2() /* print a line of
    /* printf("... after %u function calls ...\n", nf); */     /* printf("... after %u function calls ...\n", nf); */
    /* printf("... including %u linear searches ...\n", nl); */     /* printf("... including %u linear searches ...\n", nl); */
    /* printf("%10d    %10d%14.7g",nl, nf, fx); */     /* printf("%10d    %10d%14.7g",nl, nf, fx); */
   printf ( "\n" );    /* printf ( "\n" ); */
   printf ( "  Linear searches      %d", nl );    printf ( "  Linear searches      %d", nl );
     fprintf (ficlog, "  Linear searches      %d", nl );
   /* printf ( "  Linear searches      %d\n", nl ); */    /* printf ( "  Linear searches      %d\n", nl ); */
   /* printf ( "  Function evaluations %d\n", nf ); */    /* printf ( "  Function evaluations %d\n", nf ); */
   /* printf ( "  Function value FX = %g\n", fx ); */    /* printf ( "  Function value FX = %g\n", fx ); */
   printf ( "  Function evaluations %d", nf );    printf ( "  Function evaluations %d", nf );
   printf ( "  Function value FX = %.12lf\n", fx );    printf ( "  Function value FX = %.12lf\n", fx );
     fprintf (ficlog, "  Function evaluations %d", nf );
     fprintf (ficlog, "  Function value FX = %.12lf\n", fx );
 #ifdef DEBUGPRAX  #ifdef DEBUGPRAX
    printf("n=%d prin=%d\n",n,prin);     printf("n=%d prin=%d\n",n,prin);
 #endif  #endif
    if(fx <= fmin) printf(" UNDEFINED "); else  printf("%14.7g",log(fx-fmin));     /* if(fx <= fmin) printf(" UNDEFINED "); else  printf("%14.7g",log(fx-fmin)); */
    if ( n <= 4 || 2 < prin )     if ( n <= 4 || 2 < prin )
    {     {
      /* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */       /* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */
      for(i=1;i<=n;i++)printf("%14.7g",x[i]);       for(i=1;i<=n;i++){
          printf(" %14.7g",x[i]);
          fprintf(ficlog," %14.7g",x[i]);
        }
      /* r8vec_print ( n, x, "  X:" ); */       /* r8vec_print ( n, x, "  X:" ); */
    }     }
    printf("\n");     printf("\n");
      fprintf(ficlog,"\n");
  }   }
   
   
Line 3216  L1:  /* L1 or try loop */ Line 3255  L1:  /* L1 or try loop */
       if (k > 0) *d2 = 0;        if (k > 0) *d2 = 0;
    }     }
 #ifdef DEBUGPRAX  #ifdef DEBUGPRAX
    printf(" bebe end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);     printf(" bebe end of min x1 might be very wrong x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
 #endif  #endif
    if (*d2 <= small_windows) *d2 = small_windows;     if (*d2 <= small_windows) *d2 = small_windows;
    *x1 = x2; fx = fm;     *x1 = x2; fx = fm;
Line 3692  mloop: Line 3731  mloop:
    printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);      printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t); 
 #endif  #endif
    /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */     /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */
    minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global */     minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global it seems that fx doesn't correspond to f(s=*x1) */
 #ifdef DEBUGPRAX  #ifdef DEBUGPRAX
    printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx);      printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx); 
 #endif  #endif
Line 4203  void powell(double p[], double **xi, int Line 4242  void powell(double p[], double **xi, int
         printf("  + age*age  ");          printf("  + age*age  ");
         fprintf(ficlog,"  + age*age  ");          fprintf(ficlog,"  + age*age  ");
     }      }
     for(j=1;j <=ncovmodel-2;j++){      for(j=1;j <=ncovmodel-2-nagesqr;j++){
       if(Typevar[j]==0) {        if(Typevar[j]==0) {
         printf("  +      V%d  ",Tvar[j]);          printf("  +      V%d  ",Tvar[j]);
         fprintf(ficlog,"  +      V%d  ",Tvar[j]);          fprintf(ficlog,"  +      V%d  ",Tvar[j]);
Line 4411  void powell(double p[], double **xi, int Line 4450  void powell(double p[], double **xi, int
 #endif  #endif
 #ifdef POWELLORIGINAL  #ifdef POWELLORIGINAL
       if (t < 0.0) { /* Then we use it for new direction */        if (t < 0.0) { /* Then we use it for new direction */
 #else  #else  /* Not POWELLOriginal but Brouard's */
       if (directest*t < 0.0) { /* Contradiction between both tests */        if (directest*t < 0.0) { /* Contradiction between both tests */
         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("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 (if directest<0 or t<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,"directest= %.12lf (if directest<0 or t<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 (P0, Pn) for new direction Xi_n or Xi_iBig */
 #endif  #endif
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
         printf("Before linmin in direction P%d-P0\n",n);          printf("Before linmin in direction P%d-P0\n",n);
Line 4452  void powell(double p[], double **xi, int Line 4491  void powell(double p[], double **xi, int
           xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */            xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
           xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */            xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
         }          }
   
   /* #else */
   /*      for (i=1;i<=n-1;i++) {  */
   /*        for (j=1;j<=n;j++) {  */
   /*          xi[j][i]=xi[j][i+1]; /\* Standard method of conjugate directions, not Powell who changes the nth direction by p0 pn . *\/ */
   /*        } */
   /*      } */
   /*      for (j=1;j<=n;j++) {  */
   /*        xi[j][n]=xit[j];      /\* and this nth direction by the by the average p_0 p_n *\/ */
   /*      } */
   /*      /\* for (j=1;j<=n-1;j++) {  *\/ */
   /*      /\*   xi[j][1]=xi[j][j+1]; /\\* Standard method of conjugate directions *\\/ *\/ */
   /*      /\*   xi[j][n]=xit[j];      /\\* and this nth direction by the by the average p_0 p_n *\\/ *\/ */
   /*      /\* } *\/ */
   /* #endif */
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
 #else  #else
         for (j=1, flatd=0;j<=n;j++) {          for (j=1, flatd=0;j<=n;j++) {
Line 4476  void powell(double p[], double **xi, int Line 4530  void powell(double p[], double **xi, int
           free_vector(pt,1,n);             free_vector(pt,1,n); 
           return;            return;
 #endif  #endif
         }          }  /* endif(flatd >0) */
 #endif  #endif /* LINMINORIGINAL */
         printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
         fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
                   
Line 6563  void mlikeli(FILE *ficres,double p[], in Line 6617  void mlikeli(FILE *ficres,double p[], in
 #else  /* FLATSUP */  #else  /* FLATSUP */
 /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/  /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/
 /*   praxis ( t0, h0, n, prin, x, beale_f ); */  /*   praxis ( t0, h0, n, prin, x, beale_f ); */
   int prin=1;   int prin=4;
   double h0=0.25;    /* double h0=0.25; */
   double macheps;    /* double macheps; */
   double fmin;    /* double fmin; */
   macheps=pow(16.0,-13.0);    macheps=pow(16.0,-13.0);
 /* #include "praxis.h" */  /* #include "praxis.h" */
   /* Be careful that praxis start at x[0] and powell start at p[1] */    /* Be careful that praxis start at x[0] and powell start at p[1] */
Line 6576  printf("Praxis Gegenfurtner \n"); Line 6630  printf("Praxis Gegenfurtner \n");
 fprintf(ficlog, "Praxis  Gegenfurtner\n");fflush(ficlog);  fprintf(ficlog, "Praxis  Gegenfurtner\n");fflush(ficlog);
 /* praxis ( ftol, h0, npar, prin, p1, func ); */  /* praxis ( ftol, h0, npar, prin, p1, func ); */
   /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */    /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */
   fmin = praxis(ftol,macheps, h0, npar, prin, p, func);    ffmin = praxis(ftol,macheps, h0, npar, prin, p, func);
 printf("End Praxis\n");  printf("End Praxis\n");
 #endif  /* FLATSUP */  #endif  /* FLATSUP */
   
Line 8553  void  concatwav(int wav[], int **dh, int Line 8607  void  concatwav(int wav[], int **dh, int
 /************ Variance ******************/  /************ Variance ******************/
  void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres)   void varevsij(char optionfilefiname[], double ***vareij, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres)
  {   {
    /** Variance of health expectancies      /** Computes the matrix of variance covariance of health expectancies e.j= sum_i w_i e_ij where w_i depends of popbased,
       * either cross-sectional or implied.
       * return vareij[i][j][(int)age]=cov(e.i,e.j)=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20
     *  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 **newm;      * double **newm;
     * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)       * int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav) 
Line 8570  void  concatwav(int wav[], int **dh, int Line 8626  void  concatwav(int wav[], int **dh, int
    double ***gradg, ***trgradg; /**< for var eij */     double ***gradg, ***trgradg; /**< for var eij */
    double **gradgp, **trgradgp; /**< for var p point j */     double **gradgp, **trgradgp; /**< for var p point j */
    double *gpp, *gmp; /**< for var p point j */     double *gpp, *gmp; /**< for var p point j */
    double **varppt; /**< for var p point j nlstate to nlstate+ndeath */     double **varppt; /**< for var p.3 p.death nlstate+1 to nlstate+ndeath */
    double ***p3mat;     double ***p3mat;
    double age,agelim, hf;     double age,agelim, hf;
    /* double ***mobaverage; */     /* double ***mobaverage; */
Line 8638  void  concatwav(int wav[], int **dh, int Line 8694  void  concatwav(int wav[], int **dh, int
    fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");     fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
    fprintf(fichtm,"\n<br>%s  <br>\n",digitp);     fprintf(fichtm,"\n<br>%s  <br>\n",digitp);
   
    varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);     varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath); /* In fact, currently a double */
    pstamp(ficresvij);     pstamp(ficresvij);
    fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");     fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n#  (weighted average of eij where weights are ");
    if(popbased==1)     if(popbased==1)
Line 8707  void  concatwav(int wav[], int **dh, int Line 8763  void  concatwav(int wav[], int **dh, int
              prlim[i][i]=mobaverage[(int)age][i][ij];               prlim[i][i]=mobaverage[(int)age][i][ij];
          }           }
        }         }
        /**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h.         /**< Computes the shifted plus (gp) transition matrix \f$ {}{h}_p^{ij}x\f$ at horizon h.
         */                                */                      
        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Returns p3mat[i][j][h] for h=0 to nhstepm */         hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Returns p3mat[i][j][h] for h=0 to nhstepm */
        /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability         /**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}x\f$, which are the probability
Line 8716  void  concatwav(int wav[], int **dh, int Line 8772  void  concatwav(int wav[], int **dh, int
        for(j=1; j<= nlstate; j++){         for(j=1; j<= nlstate; j++){
          for(h=0; h<=nhstepm; h++){           for(h=0; h<=nhstepm; h++){
            for(i=1, gp[h][j]=0.;i<=nlstate;i++)             for(i=1, gp[h][j]=0.;i<=nlstate;i++)
              gp[h][j] += prlim[i][i]*p3mat[i][j][h];               gp[h][j] += prlim[i][i]*p3mat[i][j][h]; /* gp[h][j]= w_i h_pij */
          }           }
        }         }
        /* Next for computing shifted+ probability of death (h=1 means         /* Next for computing shifted+ probability of death (h=1 means
           computed over hstepm matrices product = hstepm*stepm months)             computed over hstepm matrices product = hstepm*stepm months) 
           as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 .            as a weighted average of prlim(i) * p(i,j) p.3=w1*p13 + w2*p23 .
        */         */
        for(j=nlstate+1;j<=nlstate+ndeath;j++){         for(j=nlstate+1;j<=nlstate+ndeath;j++){ /* Currently only once for theta plus  p.3(age) Sum_i wi pi3*/
          for(i=1,gpp[j]=0.; i<= nlstate; i++)           for(i=1,gpp[j]=0.; i<= nlstate; i++)
            gpp[j] += prlim[i][i]*p3mat[i][j][1];             gpp[j] += prlim[i][i]*p3mat[i][j][1];
        }         }
Line 8745  void  concatwav(int wav[], int **dh, int Line 8801  void  concatwav(int wav[], int **dh, int
          }           }
        }         }
                                                   
        hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);           hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres);  /* Still minus */
                                                   
        for(j=1; j<= nlstate; j++){  /* Sum of wi * eij = e.j */         for(j=1; j<= nlstate; j++){  /* gm[h][j]= Sum_i of wi * pij =  h_p.j */
          for(h=0; h<=nhstepm; h++){           for(h=0; h<=nhstepm; h++){
            for(i=1, gm[h][j]=0.;i<=nlstate;i++)             for(i=1, gm[h][j]=0.;i<=nlstate;i++)
              gm[h][j] += prlim[i][i]*p3mat[i][j][h];               gm[h][j] += prlim[i][i]*p3mat[i][j][h];
Line 8755  void  concatwav(int wav[], int **dh, int Line 8811  void  concatwav(int wav[], int **dh, int
        }         }
        /* This for computing probability of death (h=1 means         /* This for computing probability of death (h=1 means
           computed over hstepm matrices product = hstepm*stepm months)             computed over hstepm matrices product = hstepm*stepm months) 
           as a weighted average of prlim.            as a weighted average of prlim. j is death. gmp[3]=sum_i w_i*p_i3=p.3 minus theta
        */         */
        for(j=nlstate+1;j<=nlstate+ndeath;j++){         for(j=nlstate+1;j<=nlstate+ndeath;j++){  /* Currently only once theta_minus  p.3=Sum_i wi pi3*/
          for(i=1,gmp[j]=0.; i<= nlstate; i++)           for(i=1,gmp[j]=0.; i<= nlstate; i++)
            gmp[j] += prlim[i][i]*p3mat[i][j][1];             gmp[j] += prlim[i][i]*p3mat[i][j][1];
        }             }    
        /* end shifting computations */         /* end shifting computations */
   
        /**< Computing gradient matrix at horizon h          /**< Computing gradient of p.j matrix at horizon h and still for one parameter of vector theta
           * equation 31 and 32
         */          */
        for(j=1; j<= nlstate; j++) /* vareij */         for(j=1; j<= nlstate; j++) /* computes grad p.j(x, over each  h) where p.j is Sum_i w_i*pij(x over h)
                                     * equation 24 */
          for(h=0; h<=nhstepm; h++){           for(h=0; h<=nhstepm; h++){
            gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];             gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
          }           }
        /**< Gradient of overall mortality p.3 (or p.j)          /**< Gradient of overall mortality p.3 (or p.death) 
         */          */
        for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */         for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* computes grad of p.3 from wi+pi3 grad p.3 (theta) */
          gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];           gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
        }         }
                                                   
      } /* End theta */       } /* End theta */
             
      /* We got the gradient matrix for each theta and state j */                       /* We got the gradient matrix for each theta and each state j of gradg(h]theta][j)=grad(_hp.j(theta) */            
      trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */       trgradg =ma3x(0,nhstepm,1,nlstate,1,npar);
                                   
      for(h=0; h<=nhstepm; h++) /* veij */       for(h=0; h<=nhstepm; h++) /* veij */ /* computes the transposed of grad  (_hp.j(theta)*/
        for(j=1; j<=nlstate;j++)         for(j=1; j<=nlstate;j++)
          for(theta=1; theta <=npar; theta++)           for(theta=1; theta <=npar; theta++)
            trgradg[h][j][theta]=gradg[h][theta][j];             trgradg[h][j][theta]=gradg[h][theta][j];
                                   
      for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */       for(j=nlstate+1; j<=nlstate+ndeath;j++) /* computes transposed of grad p.3 (theta)*/
        for(theta=1; theta <=npar; theta++)         for(theta=1; theta <=npar; theta++)
          trgradgp[j][theta]=gradgp[theta][j];           trgradgp[j][theta]=gradgp[theta][j];
      /**< as well as its transposed matrix        /**< as well as its transposed matrix 
Line 8797  void  concatwav(int wav[], int **dh, int Line 8855  void  concatwav(int wav[], int **dh, int
          vareij[i][j][(int)age] =0.;           vareij[i][j][(int)age] =0.;
   
      /* Computing trgradg by matcov by gradg at age and summing over h       /* Computing trgradg by matcov by gradg at age and summing over h
       * and k (nhstepm) formula 15 of article        * and k (nhstepm) formula 32 of article
       * Lievre-Brouard-Heathcote        * Lievre-Brouard-Heathcote so that for each j, computes the cov(e.j,e.k) (formula 31).
         * for given h and k computes trgradg[h](i,j) matcov (theta) gradg(k)(i,j) into vareij[i][j] which is
         cov(e.i,e.j) and sums on h and k
         * including the covariances.
       */        */
             
      for(h=0;h<=nhstepm;h++){       for(h=0;h<=nhstepm;h++){
Line 8807  void  concatwav(int wav[], int **dh, int Line 8868  void  concatwav(int wav[], int **dh, int
          matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);           matprod2(doldm,dnewm,1,nlstate,1,npar,1,nlstate,gradg[k]);
          for(i=1;i<=nlstate;i++)           for(i=1;i<=nlstate;i++)
            for(j=1;j<=nlstate;j++)             for(j=1;j<=nlstate;j++)
              vareij[i][j][(int)age] += doldm[i][j]*hf*hf;               vareij[i][j][(int)age] += doldm[i][j]*hf*hf; /* This is vareij=sum_h sum_k trgrad(h_pij) V(theta) grad(k_pij)
                                                                including the covariances of e.j */
        }         }
      }       }
                                   
      /* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of       /* Mortality: pptj is p.3 or p.death = trgradgp by cov by gradgp, variance of
       * p.j overall mortality formula 49 but computed directly because        * p.3=1-p..=1-sum i p.i  overall mortality computed directly because
       * we compute the grad (wix pijx) instead of grad (pijx),even if        * we compute the grad (wix pijx) instead of grad (pijx),even if
       * wix is independent of theta.        * wix is independent of theta. 
       */        */
      matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);       matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
      matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);       matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
      for(j=nlstate+1;j<=nlstate+ndeath;j++)       for(j=nlstate+1;j<=nlstate+ndeath;j++)
        for(i=nlstate+1;i<=nlstate+ndeath;i++)         for(i=nlstate+1;i<=nlstate+ndeath;i++)
          varppt[j][i]=doldmp[j][i];           varppt[j][i]=doldmp[j][i];  /* This is the variance of p.3 */
      /* end ppptj */       /* end ppptj */
      /*  x centered again */       /*  x centered again */
                                   
Line 8843  void  concatwav(int wav[], int **dh, int Line 8905  void  concatwav(int wav[], int **dh, int
      hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres);         hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij, nres);  
      for(j=nlstate+1;j<=nlstate+ndeath;j++){       for(j=nlstate+1;j<=nlstate+ndeath;j++){
        for(i=1,gmp[j]=0.;i<= nlstate; i++)          for(i=1,gmp[j]=0.;i<= nlstate; i++) 
          gmp[j] += prlim[i][i]*p3mat[i][j][1];            gmp[j] += prlim[i][i]*p3mat[i][j][1]; /* gmp[j] is p.3 */
      }           }    
      /* end probability of death */       /* end probability of death */
                                   
      fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);       fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
      for(j=nlstate+1; j<=(nlstate+ndeath);j++){       for(j=nlstate+1; j<=(nlstate+ndeath);j++){
        fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));         fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));/* p.3 (STD p.3) */
        for(i=1; i<=nlstate;i++){         for(i=1; i<=nlstate;i++){
          fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);           fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]); /* wi, pi3 */
        }         }
      }        } 
      fprintf(ficresprobmorprev,"\n");       fprintf(ficresprobmorprev,"\n");
Line 9910  prevalence (with 95%% confidence interva Line 9972  prevalence (with 95%% confidence interva
        fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"V_"), cpt,k1,nres);
      }       }
      fprintf(fichtm,"\n<br>- Total life expectancy by age and \       fprintf(fichtm,"\n<br>- Total life expectancy by age and \
 health expectancies in each live states (1 to %d). If popbased=1 the smooth (due to the model) \  health expectancies in each live state (1 to %d) with confidence intervals \
   on left y-scale as well as proportions of time spent in each live state \
   (with confidence intervals) on right y-scale 0 to 100%%.\
    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:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);   observed and cahotic prevalences:  <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>",nlstate, subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
Line 10264  void printinggnuplot(char fileresu[], ch Line 10329  void printinggnuplot(char fileresu[], ch
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
         fprintf(ficgp,"\nset label \"popbased %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",vpopbased,gplotlabel);          fprintf(ficgp,"\nset label \"popbased %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",vpopbased,gplotlabel);
         if(vpopbased==0){          if(vpopbased==0){
           fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);            fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nunset ytics; unset y2tics; set ytics nomirror; set y2tics 0,10,100;set y2range [0:100];\nplot [%.f:%.f] ",ageminpar,fage);
         }else          }else
           fprintf(ficgp,"\nreplot ");            fprintf(ficgp,"\nreplot ");
         for (i=1; i<= nlstate+1 ; i ++) {          for (i=1; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */
           k=2*i;            k=2*i;
           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased);            fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ?$4 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */
           for (j=1; j<= nlstate+1 ; j ++) {            for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/
             if (j==i) fprintf(ficgp," %%lf (%%lf)");              if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */
             else fprintf(ficgp," %%*lf (%%*lf)");              else fprintf(ficgp," %%*lf (%%*lf)");  /* skipping that field with a star */
           }               }   
           if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);            if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l lt %d, \\\n",i);
           else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1);            else fprintf(ficgp,"\" t\"LE in state (%d)\" w l lt %d, \\\n",i-1,i+1); /* state=i-1=1 to nlstate  */
           fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);            fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && $4!=0 ? $4-$5*2 : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
           for (j=1; j<= nlstate+1 ; j ++) {            for (j=1; j<= nlstate+1 ; j ++) {
             if (j==i) fprintf(ficgp," %%lf (%%lf)");              if (j==i) fprintf(ficgp," %%lf (%%lf)");
Line 10287  void printinggnuplot(char fileresu[], ch Line 10352  void printinggnuplot(char fileresu[], ch
             if (j==i) fprintf(ficgp," %%lf (%%lf)");              if (j==i) fprintf(ficgp," %%lf (%%lf)");
             else fprintf(ficgp," %%*lf (%%*lf)");              else fprintf(ficgp," %%*lf (%%*lf)");
           }               }   
           if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");            if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0,\\\n"); /* ,\\\n added for th percentage graphs */
           else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");            else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
         } /* state */          } /* state */
           /* again for the percentag spent in state i-1=1 to i-1=nlstate */
           for (i=2; i<= nlstate+1 ; i ++) { /* For state i-1=0 is LE, while i-1=1 to nlstate are origin state */
             k=2*i;
             fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d &&  ($4)<=1 && ($4)>=0 ?($4)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1, vpopbased); /* for fixed variables age, popbased, mobilav */
             for (j=1; j<= nlstate ; j ++)
               fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
             for (j=1; j<= nlstate+1 ; j ++) { /* e.. e.1 e.2 again j-1 is the state of end, wlim_i eij*/
               if (j==i) fprintf(ficgp," %%lf (%%lf)"); /* We want to read e.. i=1,j=1, e.1 i=2,j=2, e.2 i=3,j=3 */
               else fprintf(ficgp," %%*lf (%%*lf)");  /* skipping that field with a star */
             }   
             if (i== 1) fprintf(ficgp,"\" t\"%%TLE\" w l lt %d axis x1y2, \\\n",i); /* Not used */
             else fprintf(ficgp,"\" t\"%%LE in state (%d)\" w l lw 2 lt %d axis x1y2, \\\n",i-1,i+1); /* state=i-1=1 to nlstate  */
             fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4-$5*2)<=1 && ($4-$5*2)>=0? ($4-$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
             for (j=1; j<= nlstate ; j ++)
               fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
             for (j=1; j<= nlstate+1 ; j ++) {
               if (j==i) fprintf(ficgp," %%lf (%%lf)");
               else fprintf(ficgp," %%*lf (%%*lf)");
             }   
             fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,");
             fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2==%d && ($4+$5*2)<=1 && ($4+$5*2)>=0 ? ($4+$5*2)*100. : 1/0) \"%%lf %%lf %%lf",subdirf2(fileresu,"T_"),nres-1,nres-1,vpopbased);
             for (j=1; j<= nlstate ; j ++)
               fprintf(ficgp," %%*lf (%%*lf)"); /* Skipping TLE and LE to read %LE only */
             for (j=1; j<= nlstate+1 ; j ++) {
               if (j==i) fprintf(ficgp," %%lf (%%lf)");
               else fprintf(ficgp," %%*lf (%%*lf)");
             }   
             if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2");
             else fprintf(ficgp,"\" t\"\" w l lt 0 axis x1y2,\\\n");
           } /* state for percent */
       } /* vpopbased */        } /* vpopbased */
       fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; unset label;\n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */        fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; unset label;\n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */
     } /* end nres */      } /* end nres */
Line 10696  set ter svg size 640, 480\nunset log y\n Line 10791  set ter svg size 640, 480\nunset log y\n
               /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */                /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */
               /* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */                /* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */
               kl++;                kl++;
                 /* Problem with quantitative variables TinvDoQresult[nres] */
               /* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */                /* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */
               sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,lv, kl+1, vlv );                sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%lg " ,kl,lv, kl+1, vlv );/* Solved but quantitative must be shifted */
               kl++;                kl++;
               if(k <cptcovs && cptcovs>1)                if(k <cptcovs && cptcovs>1)
                 sprintf(gplotcondition+strlen(gplotcondition)," && ");                  sprintf(gplotcondition+strlen(gplotcondition)," && ");
Line 10711  set ter svg size 640, 480\nunset log y\n Line 10807  set ter svg size 640, 480\nunset log y\n
               fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \                fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt );
               fprintf(ficgp,",\\\n '' ");                fprintf(ficgp,",\\\n '' ");
               fprintf(ficgp," u %d:(",iagec);                 fprintf(ficgp," u %d:(",iagec); /* Below iyearc should be increades if quantitative variable in the reult line */
                 /* $7==6 && $8==2.47 ) && (($9-$10) == 1953 ) ? $12/(1.-$24) : 1/0):7 with labels center not */
                 /* but was  && $7==6 && $8==2 ) && (($7-$8) == 1953 ) ? $12/(1.-$24) : 1/0):7 with labels center not */
               fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \                fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \
                       iyearc, iagec, offyear,                           \                        iyearc, iagec, offyear,                           \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc );
Line 12212  double gompertz(double x[]) Line 12310  double gompertz(double x[])
          A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));           A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
        } else if (cens[i] == 0){         } else if (cens[i] == 0){
         A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))          A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
           +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);            +log(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
           /* +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); */  /* To be seen */
       } else        } else
          printf("Gompertz cens[%d] neither 1 nor 0\n",i);           printf("Gompertz cens[%d] neither 1 nor 0\n",i);
       /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */        /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
Line 14569  int main(int argc, char *argv[]) Line 14668  int main(int argc, char *argv[])
   double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;    double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW;
   double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */    double ageminout=-AGEOVERFLOW,agemaxout=AGEOVERFLOW; /* Smaller Age range redefined after movingaverage */
   
     double stdpercent; /* for computing the std error of percent e.i: e.i/e.. */
   double fret;    double fret;
   double dum=0.; /* Dummy variable */    double dum=0.; /* Dummy variable */
   /* double*** p3mat;*/    /* double*** p3mat;*/
Line 15765  Interval (in months) between two waves: Line 15865  Interval (in months) between two waves:
     gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1)  */      gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1)  */
 #endif  #endif
 #ifdef POWELL  #ifdef POWELL
      powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);  #ifdef LINMINORIGINAL
 #endif    #else /* LINMINORIGINAL */
     
     flatdir=ivector(1,npar); 
     for (j=1;j<=npar;j++) flatdir[j]=0; 
   #endif /*LINMINORIGINAL */
       /* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */
     /* double h0=0.25; */
     macheps=pow(16.0,-13.0);
     printf("Praxis Gegenfurtner mle=%d\n",mle);
     fprintf(ficlog, "Praxis  Gegenfurtner mle=%d\n", mle);fflush(ficlog);
      /* ffmin = praxis(ftol,macheps, h0, npar, prin, p, gompertz); */
     /* For the Gompertz we use only two parameters */
     int _npar=2;
      ffmin = praxis(ftol,macheps, h0, _npar, 4, p, gompertz);
     printf("End Praxis\n");
     fclose(ficrespow);      fclose(ficrespow);
       #ifdef LINMINORIGINAL
   #else
         free_ivector(flatdir,1,npar); 
   #endif  /* LINMINORIGINAL*/
   #endif /* POWELL */   
     hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz);       hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz); 
   
     for(i=1; i <=NDIM; i++)      for(i=1; i <=NDIM; i++)
Line 15902  Please run with mle=-1 to get a correct Line 16020  Please run with mle=-1 to get a correct
       fprintf(ficlog,"  + age*age  ");        fprintf(ficlog,"  + age*age  ");
       fprintf(fichtm, "<th>+ age*age</th>");        fprintf(fichtm, "<th>+ age*age</th>");
     }      }
     for(j=1;j <=ncovmodel-2;j++){      for(j=1;j <=ncovmodel-2-nagesqr;j++){
       if(Typevar[j]==0) {        if(Typevar[j]==0) {
         printf("  +      V%d  ",Tvar[j]);          printf("  +      V%d  ",Tvar[j]);
         fprintf(ficres,"  +      V%d  ",Tvar[j]);          fprintf(ficres,"  +      V%d  ",Tvar[j]);
Line 15973  Please run with mle=-1 to get a correct Line 16091  Please run with mle=-1 to get a correct
         fprintf(ficlog,"  + age*age  ");          fprintf(ficlog,"  + age*age  ");
         fprintf(fichtm, "<th>+ age*age</th>");          fprintf(fichtm, "<th>+ age*age</th>");
       }        }
       for(j=1;j <=ncovmodel-2;j++){        for(j=1;j <=ncovmodel-2-nagesqr;j++){
         if(Typevar[j]==0) {          if(Typevar[j]==0) {
           printf("  +      V%d  ",Tvar[j]);            printf("  +      V%d  ",Tvar[j]);
           fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);            fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);
Line 16740  Please run with mle=-1 to get a correct Line 16858  Please run with mle=-1 to get a correct
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
         oldm=oldms;savm=savms; /* ZZ Segmentation fault */          oldm=oldms;savm=savms; /* ZZ Segmentation fault */
         cptcod= 0; /* To be deleted */          cptcod= 0; /* To be deleted */
         printf("varevsij vpopbased=%d \n",vpopbased);          printf("varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
         fprintf(ficlog, "varevsij vpopbased=%d \n",vpopbased);          fprintf(ficlog, "varevsij vpopbased=%d popbased=%d \n",vpopbased,popbased);
           /* Call to varevsij to get cov(e.i, e.j)= vareij[i][j][(int)age]=sum_h sum_k trgrad(h_p.i) V(theta) grad(k_p.k) Equation 20 */
           /* Depending of popbased which changes the prevalences, either cross-sectional or period */
         varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* cptcod not initialized Intel */          varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart, nres); /* cptcod not initialized Intel */
         fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n#  (weighted average of eij where weights are ");          fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each state\n\
   #  (these are weighted average of eij where weights are ");
         if(vpopbased==1)          if(vpopbased==1)
           fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);            fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally)\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
         else          else
           fprintf(ficrest,"the age specific forward period (stable) prevalences in each health state \n");            fprintf(ficrest,"the age specific forward period (stable) prevalences in each state) \n");
           fprintf(ficrest,"# with proportions of time spent in each state with standard error (on the right of the table.\n ");
         fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */          fprintf(ficrest,"# Age popbased mobilav e.. (std) "); /* Adding covariate values? */
         for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);          for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
           for (i=1;i<=nlstate;i++) fprintf(ficrest," %% e.%d/e.. (std) ",i);
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
         /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */          /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
         printf("Computing age specific forward period (stable) prevalences in each health state \n");          printf("Computing age specific forward period (stable) prevalences in each health state \n");
Line 16775  Please run with mle=-1 to get a correct Line 16898  Please run with mle=-1 to get a correct
               /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/                /*ZZZ  printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/
               /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */                /* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */
             }              }
             epj[nlstate+1] +=epj[j];              epj[nlstate+1] +=epj[j]; /* epp=sum_j epj = sum_j sum_i w_i e_ij */
           }            }
           /* printf(" age %4.0f \n",age); */            /* printf(" age %4.0f \n",age); */
                       
           for(i=1, vepp=0.;i <=nlstate;i++)            for(i=1, vepp=0.;i <=nlstate;i++)  /* Variance of total life expectancy e.. */
             for(j=1;j <=nlstate;j++)              for(j=1;j <=nlstate;j++)
               vepp += vareij[i][j][(int)age];                vepp += vareij[i][j][(int)age]; /* sum_i sum_j cov(e.i, e.j) = var(e..) */
           fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));            fprintf(ficrest," %7.3f (%7.3f)", epj[nlstate+1],sqrt(vepp));
             /* vareij[i][j] is the covariance  cov(e.i, e.j) and vareij[j][j] is the variance  of e.j  */
           for(j=1;j <=nlstate;j++){            for(j=1;j <=nlstate;j++){
             fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));              fprintf(ficrest," %7.3f (%7.3f)", epj[j],sqrt(vareij[j][j][(int)age]));
           }            }
             /* And proportion of time spent in state j */
             /* $$ E[r(X,Y)-E(r(X,Y))]^2=[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]' Var(X,Y)[\frac{1}{\mu_y} -\frac{\mu_x}{{\mu_y}^2}]$$ */
             /* \frac{\mu_x^2}{\mu_y^2} ( \frac{\sigma^2_x}{\mu_x^2}-2\frac{\sigma_{xy}}{\mu_x\mu_y} +\frac{\sigma^2_y}{\mu_y^2}) */
             /* \frac{e_{.i}^2}{e_{..}^2} ( \frac{\Var e_{.i}}{e_{.i}^2}-2\frac{\Var e_{.i} + \sum_{j\ne i} \Cov e_{.j},e_{.i}}{e_{.i}e_{..}} +\frac{\Var e_{..}}{e_{..}^2})*/
             /*\mu_x = epj[j], \sigma^2_x = vareij[j][j][(int)age] and \mu_y=epj[nlstate+1], \sigma^2_y=vepp \sigmaxy= */
             /* vareij[j][j][(int)age]/epj[nlstate+1]^2 + vepp/epj[nlstate+1]^4 */
             for(j=1;j <=nlstate;j++){
               /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */
               /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[j]/epj[j]/epj[j]/epj[j] )); */
               
               for(i=1,stdpercent=0.;i<=nlstate;i++){ /* Computing cov(e..,e.j)=cov(sum_i e.i,e.j)=sum_i cov(e.i, e.j) */
                 stdpercent += vareij[i][j][(int)age];
               }
               stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]* (vareij[j][j][(int)age]/epj[j]/epj[j]-2.*stdpercent/epj[j]/epj[nlstate+1]+ vepp/epj[nlstate+1]/epj[nlstate+1]);
               /* stdpercent= epj[j]*epj[j]/epj[nlstate+1]/epj[nlstate+1]*(vareij[j][j][(int)age]/epj[j]/epj[j] + vepp/epj[nlstate+1]/epj[nlstate+1]); */ /* Without covariance */
               /* fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt( vareij[j][j][(int)age]/epj[nlstate+1]/epj[nlstate+1] + epj[j]*epj[j]*vepp/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1]/epj[nlstate+1] )); */
               fprintf(ficrest," %7.3f (%7.3f)", epj[j]/epj[nlstate+1], sqrt(stdpercent));
             }
           fprintf(ficrest,"\n");            fprintf(ficrest,"\n");
         }          }
       } /* End vpopbased */        } /* End vpopbased */
Removed from v.1.359  
changed lines
  Added in v.1.364

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