/* $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 ***
/* #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 <math.h>
#include <stdio.h>
}
/**** praxis ****/
+# include <float.h>
void transpose_in_place ( int n, double **a )
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[])
{
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);
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
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");
/* 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);
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) {
/* 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 */
/* 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
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++) {
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);
/* 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
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) {
}
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; */
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;
/*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;
}
<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 */
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,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
}
/* State specific survival functions (period) */
for(cpt=1; cpt<=nlstate;cpt++){
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);
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);
}
/* Period (forward stable) prevalence in each health state */
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,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
}
/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
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), \
- 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. <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," <img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
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)){
#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);
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 ");
}
/* 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;
/* */
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 */
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);
git://henry.ined.fr / .git / commitdiff