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

Diff for /imach/src/imachprax.c between versions 1.2 and 1.5

version 1.2, 2023/06/22 11:22:40 version 1.5, 2023/10/09 09:10:01
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.5  2023/10/09 09:10:01  brouard
     Summary: trying to reconsider
   
     Revision 1.4  2023/06/22 12:50:51  brouard
     Summary: stil on going
   
     Revision 1.3  2023/06/22 11:28:07  brouard
     *** empty log message ***
   
   Revision 1.2  2023/06/22 11:22:40  brouard    Revision 1.2  2023/06/22 11:22:40  brouard
   Summary: with svd but not working yet    Summary: with svd but not working yet
   
Line 1281  Important routines Line 1290  Important routines
 /* #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 *\/ */
 /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */  /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
 /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */  /* #define FLATSUP  *//* Suppresses directions where likelihood is flat */
   /* #define POWELLORIGINCONJUGATE  /\* Don't use conjugate but biggest decrease if valuable *\/ */
   /* #define NOTMINFIT */
   
 #include <math.h>  #include <math.h>
 #include <stdio.h>  #include <stdio.h>
Line 2613  void linmin(double p[], double xi[], int Line 2624  void linmin(double p[], double xi[], int
 }   } 
   
 /**** praxis ****/  /**** praxis ****/
   # include <float.h>
   
   void transpose_in_place ( int n, double **a )
   
   /******************************************************************************/
   /*
     Purpose:
   
      TRANSPOSE_IN_PLACE transposes a square matrix in place.
     Licensing:
       This code is distributed under the GNU LGPL license.
   
       Input, int N, the number of rows and columns of the matrix A.
   
       Input/output, double A[N*N], the matrix to be transposed.
   */
   {
     int i;
     int j;
     double t;
   
     /* for ( j = 0; j < n; j++ ){ */
     /*   for ( i = 0; i < j; i++ ) { */
     for ( j = 1; j <= n; j++ ){
       for ( i = 1; i < j; i++ ) {
         /* t        = a[i+j*n]; */
         /* a[i+j*n] = a[j+i*n]; */
         /* a[j+i*n] = t; */
         t        = a[i][j];
         a[i][j] = a[j][i];
         a[j][i] = t;
       }
     }
     return;
   }
   
 double pythag( double x, double y )  double pythag( double x, double y )
   
 /******************************************************************************/  /******************************************************************************/
Line 2652  double pythag( double x, double y ) Line 2699  double pythag( double x, double y )
   return value;    return value;
 }  }
   
   void svsort ( int n, double d[], double **v ) 
   
   /******************************************************************************/
   /*
     Purpose:
   
       SVSORT descending sorts D and adjusts the corresponding columns of V.
   
     Discussion:
       A simple bubble sort is used on D.
       In our application, D contains singular values, and the columns of V are
       the corresponding right singular vectors.
     Author:
       Original FORTRAN77 version by Richard Brent.
       Richard Brent,
       Algorithms for Minimization with Derivatives,
       Prentice Hall, 1973,
       Reprinted by Dover, 2002.
   
     Parameters:
       Input, int N, the length of D, and the order of V.
       Input/output, double D[N], the vector to be sorted.  
       On output, the entries of D are in descending order.
   
       Input/output, double V[N,N], an N by N array to be adjusted 
       as D is sorted.  In particular, if the value that was in D(I) on input is
       moved to D(J) on output, then the input column V(*,I) is moved to
       the output column V(*,J).
   */
   {
     int i, j1, j2, j3;
     double t;
   
     for (j1 = 1; j1 < n; j1++) {
       /*
        * Find J3, the index of the largest entry in D[J1:N-1].
        * MAXLOC apparently requires its output to be an array.
       */
       j3 = j1;
       for (j2 = j1+1; j2 < n; j2++) {
         if (d[j3] < d[j2]) {
           j3 = j2;
         }
       }
       /*
        * If J1 != J3, swap D[J1] and D[J3], and columns J1 and J3 of V.
       */
       if (j1 != j3) {
         t     = d[j1];
         d[j1] = d[j3];
         d[j3] = t;
         for (i = 1; i <= n; i++) {
           t = v[i][j1];
           v[i][j1] = v[i][j3];
           v[i][j3] = t;
         } /* end i */
       } /* end j1 != j3 */
     } /* end j1 */
     return;
   }
   
 /* void svdcmp(double **a, int m, int n, double w[], double **v)  */  /* void svdcmp(double **a, int m, int n, double w[], double **v)  */
 void svdcmp(double **a, int m, int n, double w[])   void svdminfit(double **a, int m, int n, double w[]) 
 {  {
   /* Golub 1970 Algol60 */    /* From numerical recipes */
     /* Given a matrix a[1..m][1..n], this routine computes its singular value */
     /*   decomposition, A = U ·W ·V T . The matrix U replaces a on output. */
     /*   The diagonal matrix of singular values W is out- put as a vector w[1..n]. */
     /*   The matrix V (not the transpose V T ) is output as v[1..n][1..n]. */
     
     /* But in fact from Golub 1970 Algol60 */
   
 /*   Computation of the singular values and complete orthogonal decom- */  /*   Computation of the singular values and complete orthogonal decom- */
 /* position of a real rectangular matrix A, */  /* position of a real rectangular matrix A, */
 /* A = U diag (q) V^T, U^T U = V^T V =I , */  /* A = U diag (q) V^T, U^T U = V^T V =I , */
Line 2803  void svdcmp(double **a, int m, int n, do Line 2918  void svdcmp(double **a, int m, int n, do
         }           } 
         break;           break; 
       }         } 
       if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations");         if (its == 30) nrerror("no convergence in 30 svdcmp iterations"); 
       x=w[l];  /* shift from bottom 2 x 2 minor; */        x=w[l];  /* shift from bottom 2 x 2 minor; */
       nm=k-1;         nm=k-1; 
       y=w[nm];         y=w[nm]; 
Line 2918  void powell(double p[], double **xi, int Line 3033  void powell(double p[], double **xi, int
     curr_time = *localtime(&rcurr_time);      curr_time = *localtime(&rcurr_time);
     /* printf("\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); */      /* printf("\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); */
     /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); */      /* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); */
     Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /* Big iteration, i.e on ncovmodel cycle */      /* Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /\* Big iteration, i.e on ncovmodel cycle *\/ */
       Bigter=(*iter - (*iter-1) % n)/n +1; /* Big iteration, i.e on ncovmodel cycle */
     printf("\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);      printf("\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
     fprintf(ficlog,"\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);      fprintf(ficlog,"\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
     fprintf(ficrespow,"%d %d %.12f %d",*iter,Bigter, *fret,curr_time.tm_sec-start_time.tm_sec);      fprintf(ficrespow,"%d %d %.12f %d",*iter,Bigter, *fret,curr_time.tm_sec-start_time.tm_sec);
Line 2989  void powell(double p[], double **xi, int Line 3105  void powell(double p[], double **xi, int
         fprintf(ficlog,"   - if your program needs %d BIG iterations  (%d iterations) to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);          fprintf(ficlog,"   - if your program needs %d BIG iterations  (%d iterations) to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
       }        }
     }      }
     for (i=1;i<=n;i++) { /* For each direction i */      for (i=1;i<=n;i++) { /* For each direction i, maximisation after loading directions */
       for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */        for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */
       fptt=(*fret);         fptt=(*fret); 
 #ifdef DEBUG  #ifdef DEBUG
Line 2999  void powell(double p[], double **xi, int Line 3115  void powell(double p[], double **xi, int
       printf("%d",i);fflush(stdout); /* print direction (parameter) i */        printf("%d",i);fflush(stdout); /* print direction (parameter) i */
       fprintf(ficlog,"%d",i);fflush(ficlog);        fprintf(ficlog,"%d",i);fflush(ficlog);
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
       linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/        linmin(p,xit,n,fret,func); /* New point i minimizing in direction xit i has coordinates p[j].*/
         /* xit[j] gives the n coordinates of direction i as input.*/
         /* *fret gives the maximum value on direction xit */
 #else  #else
       linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/        linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
                         flatdir[i]=flat; /* Function is vanishing in that direction i */        flatdir[i]=flat; /* Function is vanishing in that direction i */
 #endif  #endif
                         /* Outputs are fret(new point p) p is updated and xit rescaled */        /* Outputs are fret(new point p) p is updated and xit rescaled */
       if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions */        if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions */
                                 /* because that direction will be replaced unless the gain del is small */          /* because that direction will be replaced unless the gain del is small */
                                 /* in comparison with the 'probable' gain, mu^2, with the last average direction. */          /* in comparison with the 'probable' gain, mu^2, with the last average direction. */
                                 /* Unless the n directions are conjugate some gain in the determinant may be obtained */          /* Unless the n directions are conjugate some gain in the determinant may be obtained */
                                 /* with the new direction. */          /* with the new direction. */
                                 del=fabs(fptt-(*fret));           del=fabs(fptt-(*fret)); 
                                 ibig=i;           ibig=i; 
       }         } 
 #ifdef DEBUG  #ifdef DEBUG
       printf("%d %.12e",i,(*fret));        printf("%d %.12e",i,(*fret));
       fprintf(ficlog,"%d %.12e",i,(*fret));        fprintf(ficlog,"%d %.12e",i,(*fret));
       for (j=1;j<=n;j++) {        for (j=1;j<=n;j++) {
                                 xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);          xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
                                 printf(" x(%d)=%.12e",j,xit[j]);          printf(" x(%d)=%.12e",j,xit[j]);
                                 fprintf(ficlog," x(%d)=%.12e",j,xit[j]);          fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
       }        }
       for(j=1;j<=n;j++) {        for(j=1;j<=n;j++) {
                                 printf(" p(%d)=%.12e",j,p[j]);          printf(" p(%d)=%.12e",j,p[j]);
                                 fprintf(ficlog," p(%d)=%.12e",j,p[j]);          fprintf(ficlog," p(%d)=%.12e",j,p[j]);
       }        }
       printf("\n");        printf("\n");
       fprintf(ficlog,"\n");        fprintf(ficlog,"\n");
 #endif  #endif
     } /* end loop on each direction i */      } /* end loop on each direction i */
     /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */       /* Convergence test will use last linmin estimation (fret) and compare to former iteration (fp) */ 
     /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */      /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */
     /* New value of last point Pn is not computed, P(n-1) */  
     for(j=1;j<=n;j++) {      for(j=1;j<=n;j++) {
       if(flatdir[j] >0){        if(flatdir[j] >0){
         printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);          printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
Line 3054  void powell(double p[], double **xi, int Line 3171  void powell(double p[], double **xi, int
       /* the scales of the directions and the directions, because the are reset to canonical directions */        /* the scales of the directions and the directions, because the are reset to canonical directions */
       /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */        /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */
       /* under the tolerance value. If the tolerance is very small 1.e-9, it could last long.  */        /* under the tolerance value. If the tolerance is very small 1.e-9, it could last long.  */
 #ifdef DEBUG  
       int k[2],l;  
       k[0]=1;  
       k[1]=-1;  
       printf("Max: %.12e",(*func)(p));  
       fprintf(ficlog,"Max: %.12e",(*func)(p));  
       for (j=1;j<=n;j++) {  
         printf(" %.12e",p[j]);  
         fprintf(ficlog," %.12e",p[j]);  
       }  
       printf("\n");  
       fprintf(ficlog,"\n");  
       for(l=0;l<=1;l++) {  
         for (j=1;j<=n;j++) {  
           ptt[j]=p[j]+(p[j]-pt[j])*k[l];  
           printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);  
           fprintf(ficlog,"l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);  
         }  
         printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));  
         fprintf(ficlog,"func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));  
       }  
 #endif  
   
       free_vector(xit,1,n);         free_vector(xit,1,n); 
       free_vector(xits,1,n);         free_vector(xits,1,n); 
Line 3084  void powell(double p[], double **xi, int Line 3179  void powell(double p[], double **xi, int
       return;         return; 
     } /* enough precision */       } /* enough precision */ 
     if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations.");       if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations."); 
     for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */      for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0)=2Pn-P0 */
       ptt[j]=2.0*p[j]-pt[j];         ptt[j]=2.0*p[j]-pt[j]; 
       xit[j]=p[j]-pt[j];         xit[j]=p[j]-pt[j]; /* Coordinate j of last direction xi_n=P_n-_0 */
         printf("\n %d xit=%12.7g p=%12.7g pt=%12.7g ",j,xit[j],p[j],pt[j]);
       pt[j]=p[j];         pt[j]=p[j]; 
     }       }
       printf("\n");
     fptt=(*func)(ptt); /* f_3 */      fptt=(*func)(ptt); /* f_3 */
 #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in drections until some iterations are done */  #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in drections until some iterations are done */
                 if (*iter <=4) {                  if (*iter <=4) {
Line 3109  void powell(double p[], double **xi, int Line 3206  void powell(double p[], double **xi, int
       /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */        /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
       /*  Even if f3 <f1, directest can be negative and t >0 */        /*  Even if f3 <f1, directest can be negative and t >0 */
       /* mu² and del² are equal when f3=f1 */        /* mu² and del² are equal when f3=f1 */
                         /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */        /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */
                         /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */        /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */
                         /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0  */        /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0  */
                         /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0  */        /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0  */
 #ifdef NRCORIGINAL  #ifdef NRCORIGINAL
       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/
 #else  #else
Line 3120  void powell(double p[], double **xi, int Line 3217  void powell(double p[], double **xi, int
       t= t- del*SQR(fp-fptt);        t= t- del*SQR(fp-fptt);
 #endif  #endif
       directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta 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        printf(" t=%g, directest=%g\n",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);  #ifdef POWELLNOTTRUECONJUGATE   /* Searching for IBIG and testing for replacement */  
       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);  
       printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),  
              (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));  
       fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),  
              (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));  
       printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);  
       fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);  
 #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 /* end POWELLOriginal */
 #ifdef DEBUGLINMIN  #endif  /* POWELLNOTTRUECONJUGATE  else means systematic replacement by new direction P_0P_n */
         printf("Before linmin in direction P%d-P0\n",n);  
         for (j=1;j<=n;j++) {  
           printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);  
           fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);  
           if(j % ncovmodel == 0){  
             printf("\n");  
             fprintf(ficlog,"\n");  
           }  
         }  
 #endif  
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
         linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/          /* xit[j]=p[j]-pt[j] */
 #else          printf("\n Computes min on P_0, P_n direction iter=%d Bigter=%d\n",*iter,Bigter);
           linmin(p,xit,n,fret,func); /* computes minimum on P_0,P_n direction: changes p and rescales xit.*/
   #else  /* NOT LINMINORIGINAL but with searching for flat directions */
           printf("\n Flat Computes min on P_0, P_n direction iter=%d Bigter=%d\n",*iter,Bigter);
         linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/          linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
         flatdir[i]=flat; /* Function is vanishing in that direction i */          flatdir[i]=flat; /* Function is vanishing in that direction i */
 #endif  #endif
                   
 #ifdef DEBUGLINMIN  #ifdef POWELLNOTTRUECONJUGATE
   #else
   #ifdef POWELLORIGINCONJUGATE
         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]);            xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
           fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);            xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
           if(j % ncovmodel == 0){          }
             printf("\n");  #else
             fprintf(ficlog,"\n");          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 . */
           }            }
         }          }
 #endif  
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
           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 */
         }          }
   #endif  /* POWELLORIGINCONJUGATE*/
   #endif  /*POWELLNOTTRUECONJUGATE*/
           printf(" Standard method of conjugate directions\n");
           printf("\n#A Before prax Bigter=%d model=  1      +     age ", Bigter);
           for(j=1;j<=n;j++){
             printf("%d \n",j);
             for(i=1;i<=n;i++){
               printf(" %f",xi[j][i]);
             }
           }
           printf("\n");
   
   #ifdef NOTMINFIT
   #else 
           if(*iter >n){
           /* if(Bigter >n){ */
             printf("\n#Before prax Bigter=%d model=  1      +     age ", Bigter);
             printf("\n");
             for(j=1;j<=n;j++){
               printf("%d \n",j);
               for(i=1;i<=n;i++){
                 printf(" %f",xi[j][i]);
               }
             }
             printf("\n");
             /*
              * Calculate a new set of orthogonal directions before repeating
              * the main loop.
              * Transpose V for SVD (minfit) (because minfit returns the right V in ULV=A):
              */
             printf(" Bigter=%d Calculate a new set of orthogonal directions before repeating  the main loop.\n  Transpose V for MINFIT:...\n",Bigter);
             transpose_in_place ( n, xi );
             /*
               SVD/MINFIT finds the singular value decomposition of V.
               
               This gives the principal values and principal directions of the
               approximating quadratic form without squaring the condition number.
             */
             printf(" SVDMINFIT finds the singular value decomposition of V. \n This gives the principal values and principal directions of the\n  approximating quadratic form without squaring the condition number...\n");
             double *d; /* eigenvalues of principal directions */
             d=vector(1,n);
             
             
             svdminfit (xi, n, n, d ); /* In Brent's notation find d such that V=Q Diagonal(d) R, and Lambda=d^(-1/2) */
           
             printf("\n#After prax model=  1      +     age ");
             fprintf(ficlog,"\n#model=  1      +     age ");
             
             if(nagesqr==1){
               printf("  + age*age  ");
               fprintf(ficlog,"  + age*age  ");
             }
             for(j=1;j <=ncovmodel-2;j++){
               if(Typevar[j]==0) {
                 printf("  +      V%d  ",Tvar[j]);
                 fprintf(ficlog,"  +      V%d  ",Tvar[j]);
               }else if(Typevar[j]==1) {
                 printf("  +    V%d*age ",Tvar[j]);
                 fprintf(ficlog,"  +    V%d*age ",Tvar[j]);
               }else if(Typevar[j]==2) {
                 printf("  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
                 fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
               }else if(Typevar[j]==3) {
                 printf("  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
                 fprintf(ficlog,"  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
               }
             }
             printf("\n");
             fprintf(ficlog,"\n");
             for(i=1,jk=1; i <=nlstate; i++){
               for(k=1; k <=(nlstate+ndeath); k++){
                 if (k != i) {
                   printf("%d%d ",i,k);
                   fprintf(ficlog,"%d%d ",i,k);
                   for(j=1; j <=ncovmodel; j++){
                     printf("%12.7f ",p[jk]);
                     fprintf(ficlog,"%12.7f ",p[jk]);
                     jk++; 
                   }
                   printf("\n");
                   fprintf(ficlog,"\n");
                 }
               }
             }
             /* minfit ( n, vsmall, v, d ); */
             /* v is overwritten with R. */
             /*
               Heuristic numbers:
               If the axes may be badly scaled (which is to be avoided if
               possible), then set SCBD = 10.  Otherwise set SCBD = 1.
               
               If the problem is known to be ill-conditioned, initialize ILLC = true.
               KTM is the number of iterations without improvement before the
               algorithm terminates.  KTM = 4 is very cautious; usually KTM = 1
               is satisfactory.
             */
             double machep, small;
             double dmin;
             int illc=0; /*    Local, int ILLC, is TRUE if the system is ill-conditioned. */
             machep = DBL_EPSILON;
             small = machep * machep;
             /* m2 = dsqrt(machep); */
             
             /*
              * Sort the eigenvalues and eigenvectors.
              */
             printf(" Sort the eigenvalues and eigenvectors....\n");
             svsort ( n, d, xi );
             printf("Sorted Eigenvalues:\n");
             for(i=1; i<=n;i++){
               printf(" d[%d]=%g",i,d[i]);
             }
             printf("\n");
             /*
              * Determine the smallest eigenvalue.
              */
             printf("  Determine the smallest eigenvalue.\n");
             dmin = fmax ( d[n], small );
             /*
              * The ratio of the smallest to largest eigenvalue determines whether
              * the system is ill conditioned.
              */
             if ( dmin < sqrt(machep) * d[1] ) {  /* m2*d[0] */
               illc = 1;
             } else  {
               illc = 0;
             }
             printf("  The ratio of the smallest to largest eigenvalue determines whether\n \
     the system is ill conditioned=%d . dmin=%.12lf < m2=%.12lf * d[1]=%.12lf \n",illc, dmin,sqrt(machep), d[1]);
             /*    if ( 1.0 < scbd ) { */
             /*      r8vec_print ( n, z, "  The scale factors:" ); */
             /*    }  */
             /*    r8vec_print ( n, d, "  Principal values of the quadratic form:" ); */
             /* } */
             /* if ( 3 < prin ) { */
             /*    r8mat_print ( n, n, v, "  The principal axes:" ); */
             /* } */
             free_vector(d,1,n);   
             /*
              * The main loop ends here.
              */
             
             /* if ( 0 < prin ) */
             /* { */
             /*   r8vec_print ( n, x, "  X:" ); */
             /* } */
           }       
   #endif  /* NOTMINFIT */
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
 #else  #else
         for (j=1, flatd=0;j<=n;j++) {          for (j=1, flatd=0;j<=n;j++) {
Line 3197  void powell(double p[], double **xi, int Line 3429  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
         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);
             /* The minimization in direction $\xi_1$ gives $P_1$. From $P_1$ minimization in direction $\xi_2$ gives */
 #ifdef DEBUG    /* $P_2$. Minimization of line $P_2-P_1$ gives new starting point $P^{(1)}_0$ and direction $\xi_1$ is dropped and replaced by second */
         printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);    /* direction $\xi_1^{(1)}=\xi_2$. Also second direction is replaced by new direction $\xi^{(1)}_2=P_2-P_0$. */
         fprintf(ficlog,"Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);  
         for(j=1;j<=n;j++){    /* At the second iteration, starting from $P_0^{(1)}$, minimization amongst $\xi^{(1)}_1$ gives point $P^{(1)}_1$. */
           printf(" %lf",xit[j]);    /* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$.  As $P^{(2)}_1$ and */
           fprintf(ficlog," %lf",xit[j]);    /* $P^{(1)}_0$ are minimizing in the same direction $P^{(1)}_2 - P^{(1)}_1= P_2-P_0$, directions $P^{(1)}_2-P^{(1)}_0$ */
         }    /* and $P_2-P_0$ (parallel to $\xi$ and $\xi^c$) are conjugate.  } */
         printf("\n");  #ifdef POWELLNOTTRUECONJUGATE  
         fprintf(ficlog,"\n");  
 #endif  
       } /* end of t or directest negative */        } /* end of t or directest negative */
   #endif
 #ifdef POWELLNOF3INFF1TEST  #ifdef POWELLNOF3INFF1TEST
 #else  #else
       } /* end if (fptt < fp)  */        } /* end if (fptt < fp)  */
 #endif  #endif
   #ifdef POWELLORIGINCONJUGATE
       } /* if (t < 0.0) */
   #endif
 #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in drections until some iterations are done */  #ifdef NODIRECTIONCHANGEDUNTILNITER  /* No change in drections until some iterations are done */
     } /*NODIRECTIONCHANGEDUNTILNITER  No change in drections until some iterations are done */      } /*NODIRECTIONCHANGEDUNTILNITER  No change in drections until some iterations are done */
 #else  #else
Line 4811  double funcone( double *x) Line 5045  double funcone( double *x)
         * 3 ncovta=15    +age*V3*V2+age*V2+agev3+ageV4 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4          * 3 ncovta=15    +age*V3*V2+age*V2+agev3+ageV4 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
         * 3 TvarAVVA[1]@15= itva 3 2    2      3    4        6       7        6 3         7 3         6 4         7 4           * 3 TvarAVVA[1]@15= itva 3 2    2      3    4        6       7        6 3         7 3         6 4         7 4 
         * 3 ncovta             1 2      3      4    5        6       7        8 9       10 11       12 13        14 15          * 3 ncovta             1 2      3      4    5        6       7        8 9       10 11       12 13        14 15
         *           *?TvarAVVAind[1]@15= V3 is in k=2 1 1  2    3        4       5        4,2         5,2,      4,3           5 3}TvarVVAind[]
         * TvarAVVAind[1]@15= V3 is in k=6 6 12  13   14      15      16       18 18       19,19,     20,20        21,21}TvarVVAind[]          * TvarAVVAind[1]@15= V3 is in k=6 6 12  13   14      15      16       18 18       19,19,     20,20        21,21}TvarVVAind[]
         * 3 ncovvta=10     +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4          * 3 ncovvta=10     +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4
         * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva]          * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva]
Line 4825  double funcone( double *x) Line 5059  double funcone( double *x)
         *                   2, 3, 4, 6, 7,          *                   2, 3, 4, 6, 7,
         *                     6, 8, 9, 10, 11}          *                     6, 8, 9, 10, 11}
         * TvarFind[itv]                        0      0       0          * TvarFind[itv]                        0      0       0
         * FixedV[itv]                          1      1       1  0      1 0       1 0       1 0      1 0     1 0          * FixedV[itv]                          1      1       1  0      1 0       1 0       1 0       0
           *? FixedV[itv]                          1      1       1  0      1 0       1 0       1 0      1 0     1 0
         * Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4          * Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4
         * Tvar[TvarFind[itv]]                [0]=?      ?ncovv 1 à ncovvt]          * Tvar[TvarFind[itv]]                [0]=?      ?ncovv 1 à ncovvt]
         *   Not a fixed cotvar[mw][itv][i]     6       7      6  2      7, 2,     6, 3,     7, 3,     6, 4,     7, 4}          *   Not a fixed cotvar[mw][itv][i]     6       7      6  2      7, 2,     6, 3,     7, 3,     6, 4,     7, 4}
Line 5210  void mlikeli(FILE *ficres,double p[], in Line 5445  void mlikeli(FILE *ficres,double p[], in
   
   
   xi=matrix(1,npar,1,npar);    xi=matrix(1,npar,1,npar);
   for (i=1;i<=npar;i++)    for (i=1;i<=npar;i++)  /* Starting with canonical directions j=1,n xi[i=1,n][j] */
     for (j=1;j<=npar;j++)      for (j=1;j<=npar;j++)
       xi[i][j]=(i==j ? 1.0 : 0.0);        xi[i][j]=(i==j ? 1.0 : 0.0);
   printf("Powell\n");  fprintf(ficlog,"Powell\n");    printf("Powell-prax\n");  fprintf(ficlog,"Powell-prax\n");
   strcpy(filerespow,"POW_");     strcpy(filerespow,"POW_"); 
   strcat(filerespow,fileres);    strcat(filerespow,fileres);
   if((ficrespow=fopen(filerespow,"w"))==NULL) {    if((ficrespow=fopen(filerespow,"w"))==NULL) {
Line 5277  void mlikeli(FILE *ficres,double p[], in Line 5512  void mlikeli(FILE *ficres,double p[], in
   }    }
   powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);    powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
 #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=4;  /*   int prin=4; */
   double h0=0.25;  /*   double h0=0.25; */
 #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] */
    /* praxis ( ftol, h0, npar, prin, p, func ); */     /* praxis ( ftol, h0, npar, prin, p, func ); */
 p1= (p+1); /*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */  /* p1= (p+1); /\*  p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] *\/ */
 printf("Praxis \n");  /* printf("Praxis \n"); */
 fprintf(ficlog, "Praxis \n");fflush(ficlog);  /* fprintf(ficlog, "Praxis \n");fflush(ficlog); */
 praxis ( ftol, h0, npar, prin, p1, func );  /* praxis ( ftol, h0, npar, prin, p1, func ); */
 printf("End Praxis\n");  /* printf("End Praxis\n"); */
 #endif  /* FLATSUP */  #endif  /* FLATSUP */
   
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
Line 6658  void  concatwav(int wav[], int **dh, int Line 6893  void  concatwav(int wav[], int **dh, int
             if(j==0) j=1;  /* Survives at least one month after exam */              if(j==0) j=1;  /* Survives at least one month after exam */
             else if(j<0){              else if(j<0){
               nberr++;                nberr++;
               printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
               j=1; /* Temporary Dangerous patch */                j=1; /* Temporary Dangerous patch */
               printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);                printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
               fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
               fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);                fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
             }              }
             k=k+1;              k=k+1;
Line 6695  void  concatwav(int wav[], int **dh, int Line 6930  void  concatwav(int wav[], int **dh, int
           /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/            /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
           if(j<0){            if(j<0){
             nberr++;              nberr++;
             printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);              printf("Error! Negative delay (%d) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
             fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);              fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
           }            }
           sum=sum+j;            sum=sum+j;
         }          }
Line 8443  divided by h: <sub>h</sub>P<sub>ij</sub> Line 8678  divided by h: <sub>h</sub>P<sub>ij</sub>
 <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);   <img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); 
      /* Survival functions (period) in state j */       /* Survival functions (period) in state j */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);         fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
        fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));         fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
        fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
      }       }
      /* State specific survival functions (period) */       /* State specific survival functions (period) */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\         fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\
  And probability to be observed in various states (up to %d) being in state %d at different ages.       \   And probability to be observed in various states (up to %d) being in state %d at different ages.  Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. \
  <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);   <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
        fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));         fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
        fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);         fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
      }       }
      /* Period (forward stable) prevalence in each health state */       /* Period (forward stable) prevalence in each health state */
      for(cpt=1; cpt<=nlstate;cpt++){       for(cpt=1; cpt<=nlstate;cpt++){
        fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);         fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be alive in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
        fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));         fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
       fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);        fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
      }       }
Line 8482  divided by h: <sub>h</sub>P<sub>ij</sub> Line 8717  divided by h: <sub>h</sub>P<sub>ij</sub>
       /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */        /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){         for(cpt=1; cpt<=nlstate;cpt++){
          fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \           fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
  from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \   from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d). Randomness in cross-sectional prevalence is not taken into \
  account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \   account but can visually be appreciated. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
 with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);  with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
          fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_"));           fprintf(fichtm," (data from text file  <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_"));
          fprintf(fichtm," <img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);           fprintf(fichtm," <img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
Line 13388  int main(int argc, char *argv[]) Line 13623  int main(int argc, char *argv[])
   getcwd(pathcd, size);    getcwd(pathcd, size);
 #endif  #endif
   syscompilerinfo(0);    syscompilerinfo(0);
   printf("\nIMaCh prax version %s, %s\n%s",version, copyright, fullversion);    printf("\nIMaCh prax version minfit %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: ");
     if(!fgets(pathr,FILENAMELENGTH,stdin)){      if(!fgets(pathr,FILENAMELENGTH,stdin)){
Line 14381  Interval (in months) between two waves: Line 14616  Interval (in months) between two waves:
 #ifdef GSL  #ifdef GSL
     printf("GSL optimization\n");  fprintf(ficlog,"Powell\n");      printf("GSL optimization\n");  fprintf(ficlog,"Powell\n");
 #else  #else
     printf("Powell\n");  fprintf(ficlog,"Powell\n");      printf("Powell-mort\n");  fprintf(ficlog,"Powell-mort\n");
 #endif  #endif
     strcpy(filerespow,"POW-MORT_");       strcpy(filerespow,"POW-MORT_"); 
     strcat(filerespow,fileresu);      strcat(filerespow,fileresu);
Line 14484  Interval (in months) between two waves: Line 14719  Interval (in months) between two waves:
   
     for(i=1; i <=NDIM; i++)      for(i=1; i <=NDIM; i++)
       for(j=i+1;j<=NDIM;j++)        for(j=i+1;j<=NDIM;j++)
                                 matcov[i][j]=matcov[j][i];          matcov[i][j]=matcov[j][i];
           
     printf("\nCovariance matrix\n ");      printf("\nCovariance matrix\n ");
     fprintf(ficlog,"\nCovariance matrix\n ");      fprintf(ficlog,"\nCovariance matrix\n ");
Line 14938  Please run with mle=-1 to get a correct Line 15173  Please run with mle=-1 to get a correct
     }      }
             
     /* Results */      /* Results */
     /* Value of covariate in each resultine will be compututed (if product) and sorted according to model rank */      /* Value of covariate in each resultine will be computed (if product) and sorted according to model rank */
     /* It is precov[] because we need the varying age in order to compute the real cov[] of the model equation */        /* It is precov[] because we need the varying age in order to compute the real cov[] of the model equation */  
     precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1);      precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1);
     endishere=0;      endishere=0;
Line 15341  Please run with mle=-1 to get a correct Line 15576  Please run with mle=-1 to get a correct
       /* */        /* */
       if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */        if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
         continue;          continue;
       printf("\n# model %s \n#****** Result for:", model);  /* HERE model is empty */        printf("\n# model=1+age+%s \n#****** Result for:", model);  /* HERE model is empty */
       fprintf(ficrest,"\n# model %s \n#****** Result for:", model);        fprintf(ficrest,"\n# model=1+age+%s \n#****** Result for:", model);
       fprintf(ficlog,"\n# model %s \n#****** Result for:", model);        fprintf(ficlog,"\n# model=1+age+%s \n#****** Result for:", model);
       /* It might not be a good idea to mix dummies and quantitative */        /* It might not be a good idea to mix dummies and quantitative */
       /* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model *\/ */        /* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model *\/ */
       for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */        for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */
Line 15538  Please run with mle=-1 to get a correct Line 15773  Please run with mle=-1 to get a correct
     free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);      free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
   }  /* mle==-3 arrives here for freeing */    }  /* mle==-3 arrives here for freeing */
   /* endfree:*/    /* endfree:*/
     if(mle!=-3) free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/
   free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
   free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);    free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
Removed from v.1.2  
changed lines
  Added in v.1.5

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