From: N. Brouard Date: Mon, 9 Oct 2023 09:10:01 +0000 (+0000) Subject: Summary: trying to reconsider X-Git-Tag: imach-099s7~34 X-Git-Url: https://henry.ined.fr/git/?a=commitdiff_plain;h=08f1caf3e4643f505a8ff77017d121b5b25bfc87;p=.git Summary: trying to reconsider --- diff --git a/src/imachprax.c b/src/imachprax.c index ae5bade..06bea6f 100644 --- a/src/imachprax.c +++ b/src/imachprax.c @@ -1,6 +1,9 @@ /* $Id$ $State$ $Log$ + 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 *** @@ -1285,6 +1288,7 @@ Important routines /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */ /* #define FLATSUP *//* Suppresses directions where likelihood is flat */ /* #define POWELLORIGINCONJUGATE /\* Don't use conjugate but biggest decrease if valuable *\/ */ +/* #define NOTMINFIT */ #include #include @@ -2617,6 +2621,7 @@ void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [ } /**** praxis ****/ +# include void transpose_in_place ( int n, double **a ) @@ -2691,6 +2696,67 @@ double pythag( double x, double y ) 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 svdminfit(double **a, int m, int n, double w[]) { @@ -2964,7 +3030,8 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret 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); */ /* 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); 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); @@ -3035,7 +3102,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret 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 */ fptt=(*fret); #ifdef DEBUG @@ -3045,33 +3112,33 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret printf("%d",i);fflush(stdout); /* print direction (parameter) i */ fprintf(ficlog,"%d",i);fflush(ficlog); #ifdef LINMINORIGINAL - linmin(p,xit,n,fret,func); /* New point i minimizing in direction i has coordinates p[j].*/ + 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 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 - /* 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 */ - /* 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. */ - /* Unless the n directions are conjugate some gain in the determinant may be obtained */ - /* with the new direction. */ - del=fabs(fptt-(*fret)); - ibig=i; + /* 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. */ + /* Unless the n directions are conjugate some gain in the determinant may be obtained */ + /* with the new direction. */ + del=fabs(fptt-(*fret)); + ibig=i; } #ifdef DEBUG printf("%d %.12e",i,(*fret)); fprintf(ficlog,"%d %.12e",i,(*fret)); for (j=1;j<=n;j++) { - xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5); - printf(" x(%d)=%.12e",j,xit[j]); - fprintf(ficlog," x(%d)=%.12e",j,xit[j]); + xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5); + printf(" x(%d)=%.12e",j,xit[j]); + fprintf(ficlog," x(%d)=%.12e",j,xit[j]); } for(j=1;j<=n;j++) { - printf(" p(%d)=%.12e",j,p[j]); - fprintf(ficlog," p(%d)=%.12e",j,p[j]); + printf(" p(%d)=%.12e",j,p[j]); + fprintf(ficlog," p(%d)=%.12e",j,p[j]); } printf("\n"); fprintf(ficlog,"\n"); @@ -3101,28 +3168,6 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret /* 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 */ /* 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(xits,1,n); @@ -3131,11 +3176,13 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret return; } /* enough precision */ 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]; - 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]; - } + } + printf("\n"); fptt=(*func)(ptt); /* f_3 */ #ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */ if (*iter <=4) { @@ -3156,10 +3203,10 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */ /* Even if f3 0 */ /* mu² and del² are equal when f3=f1 */ - /* 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 : 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 : 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 : lambda² < mu^2 < del then t is negative and directest >0 */ + /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */ #ifdef NRCORIGINAL t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/ #else @@ -3167,86 +3214,194 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret t= t- del*SQR(fp-fptt); #endif directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */ -#ifdef DEBUG - printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); - 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 + printf(" t=%g, directest=%g\n",t, directest); +#ifdef POWELLNOTTRUECONJUGATE /* Searching for IBIG and testing for replacement */ #ifdef POWELLORIGINAL 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 */ - 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); 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); } - if (directest < 0.0) { /* Then we use it for new direction */ -#endif -#ifdef DEBUGLINMIN - 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 + if (directest < 0.0) { /* Then we use (P0, Pn) for new direction Xi_n or Xi_iBig */ +#endif /* end POWELLOriginal */ +#endif /* POWELLNOTTRUECONJUGATE else means systematic replacement by new direction P_0P_n */ #ifdef LINMINORIGINAL - linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ -#else + /* xit[j]=p[j]-pt[j] */ + 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.*/ flatdir[i]=flat; /* Function is vanishing in that direction i */ #endif -#ifdef DEBUGLINMIN - for (j=1;j<=n;j++) { - printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); - fprintf(ficlog,"After 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 POWELLNOTTRUECONJUGATE +#else #ifdef POWELLORIGINCONJUGATE 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 */ } #else - for (j=1;j<=n-1;j++) { - xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */ + 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 */ } -/* - 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(" Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n"); - transpose_in_place ( n, xi ); -/* - SVD/MINFIT finds the singular value decomposition of V. +#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"); - 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 *w; /* eigenvalues of principal directions */ - w=vector(1,n); - svdminfit (xi, n, n, w ); /* In Brent's notation find d such that V=Q Diagonal(d) R, and Lambda=d^(-1/2) */ - /* minfit ( n, vsmall, v, d ); */ - /* v is overwritten with R. */ - free_vector(w,1,n); -#endif +#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 #else for (j=1, flatd=0;j<=n;j++) { @@ -3271,7 +3426,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret free_vector(pt,1,n); return; #endif - } + } /* endif(flatd >0) */ #endif 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); @@ -3283,23 +3438,16 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret /* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$. As $P^{(2)}_1$ and */ /* $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. } */ - - -#ifdef DEBUG - printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); - 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++){ - printf(" %lf",xit[j]); - fprintf(ficlog," %lf",xit[j]); - } - printf("\n"); - fprintf(ficlog,"\n"); -#endif +#ifdef POWELLNOTTRUECONJUGATE } /* end of t or directest negative */ +#endif #ifdef POWELLNOF3INFF1TEST #else } /* end if (fptt < fp) */ #endif +#ifdef POWELLORIGINCONJUGATE + } /* if (t < 0.0) */ +#endif #ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */ } /*NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done */ #else @@ -5297,7 +5445,7 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub for (i=1;i<=npar;i++) /* Starting with canonical directions j=1,n xi[i=1,n][j] */ for (j=1;j<=npar;j++) 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_"); strcat(filerespow,fileres); if((ficrespow=fopen(filerespow,"w"))==NULL) { @@ -5361,7 +5509,7 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub } powell(p,xi,npar,ftol,&iter,&fret,flatdir,func); #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 ); */ /* int prin=4; */ /* double h0=0.25; */ @@ -6742,10 +6890,10 @@ void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc if(j==0) j=1; /* Survives at least one month after exam */ else if(j<0){ 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 */ 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); } k=k+1; @@ -6779,8 +6927,8 @@ void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc /*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){ 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]); - 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]); + 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 (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; } @@ -8527,21 +8675,21 @@ divided by h: hPij/h : %s_%d-3- ",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres); /* Survival functions (period) in state j */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. %s_%d-%d-%d.svg
", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); + fprintf(fichtm,"
\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. %s_%d-%d-%d.svg
", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_")); fprintf(fichtm,"",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres); } /* State specific survival functions (period) */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
\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. \ %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_")); fprintf(fichtm,"",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres); } /* Period (forward stable) prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ - fprintf(fichtm,"
\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. %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); + fprintf(fichtm,"
\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. %s_%d-%d-%d.svg
", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); fprintf(fichtm," (data from text file %s.txt)\n
",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_")); fprintf(fichtm,"" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres); } @@ -8566,8 +8714,8 @@ divided by h: hPij/h : %s_%d-3- /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */ for(cpt=1; cpt<=nlstate;cpt++){ fprintf(fichtm,"
\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 \ - 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) \ + 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) \ with weights corresponding to observed prevalence at different ages. %s_%d-%d-%d.svg", 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 %s.txt)\n
",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_")); fprintf(fichtm," ", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres); @@ -13472,7 +13620,7 @@ int main(int argc, char *argv[]) getcwd(pathcd, size); #endif 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){ printf("\nEnter the parameter file name: "); if(!fgets(pathr,FILENAMELENGTH,stdin)){ @@ -14465,7 +14613,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ #ifdef GSL printf("GSL optimization\n"); fprintf(ficlog,"Powell\n"); #else - printf("Powell\n"); fprintf(ficlog,"Powell\n"); + printf("Powell-mort\n"); fprintf(ficlog,"Powell-mort\n"); #endif strcpy(filerespow,"POW-MORT_"); strcat(filerespow,fileresu); @@ -14568,7 +14716,7 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ for(i=1; i <=NDIM; i++) for(j=i+1;j<=NDIM;j++) - matcov[i][j]=matcov[j][i]; + matcov[i][j]=matcov[j][i]; printf("\nCovariance matrix\n "); fprintf(ficlog,"\nCovariance matrix\n "); @@ -15022,7 +15170,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa } /* 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 */ precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1); endishere=0; @@ -15425,9 +15573,9 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa /* */ if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */ continue; - printf("\n# model %s \n#****** Result for:", model); /* HERE model is empty */ - fprintf(ficrest,"\n# model %s \n#****** Result for:", model); - fprintf(ficlog,"\n# model %s \n#****** Result for:", model); + printf("\n# model=1+age+%s \n#****** Result for:", model); /* HERE model is empty */ + fprintf(ficrest,"\n# model=1+age+%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 */ /* 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 */ @@ -15622,6 +15770,7 @@ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpa free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); } /* mle==-3 arrives here for freeing */ /* 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(newms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);