/* $Id$
$State$
$Log$
+ Revision 1.182 2015/02/12 08:19:57 brouard
+ Summary: Trying to keep directest which seems simpler and more general
+ Author: Nicolas Brouard
+
Revision 1.181 2015/02/11 23:22:24 brouard
Summary: Comments on Powell added
*/
#define POWELL /* Instead of NLOPT */
-#define POWELLDIRECT /* Directest to decide new direction instead of Powell test */
+/* #define POWELLORIGINAL */ /* Don't use Directest to decide new direction but original Powell test */
#include <math.h>
#include <stdio.h>
/* $Id$ */
/* $State$ */
-char version[]="Imach version 0.98p, Février 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
+char version[]="Imach version 0.98q0, March 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
#define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
#define rint(a) floor(a+0.5)
/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */
-/* #define mytinydouble 1.0e-16 */
+#define mytinydouble 1.0e-16
/* #define DEQUAL(a,b) (fabs((a)-(b))<mytinydouble) */
/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/mynrutils_8h-source.html */
/* static double dsqrarg; */
if (fu <= fx) {
if (u >= x) a=x; else b=x;
SHFT(v,w,x,u)
- SHFT(fv,fw,fx,fu)
- } else {
- if (u < x) a=u; else b=u;
- if (fu <= fw || w == x) {
- v=w;
- w=u;
- fv=fw;
- fw=fu;
- } else if (fu <= fv || v == x || v == w) {
- v=u;
- fv=fu;
- }
- }
+ SHFT(fv,fw,fx,fu)
+ } else {
+ if (u < x) a=u; else b=u;
+ if (fu <= fw || w == x) {
+ v=w;
+ w=u;
+ fv=fw;
+ fw=fu;
+ } else if (fu <= fv || v == x || v == w) {
+ v=u;
+ fv=fu;
+ }
+ }
}
nrerror("Too many iterations in brent");
*xmin=x;
void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
double (*func)(double))
-{
+{ /* Given a function func , and given distinct initial points ax and bx , this routine searches in
+the downhill direction (defined by the function as evaluated at the initial points) and returns
+new points ax , bx , cx that bracket a minimum of the function. Also returned are the function
+values at the three points, fa, fb , and fc such that fa > fb and fb < fc.
+ */
double ulim,u,r,q, dum;
double fu;
*fb=(*func)(*bx);
if (*fb > *fa) {
SHFT(dum,*ax,*bx,dum)
- SHFT(dum,*fb,*fa,dum)
- }
+ SHFT(dum,*fb,*fa,dum)
+ }
*cx=(*bx)+GOLD*(*bx-*ax);
*fc=(*func)(*cx);
- while (*fb > *fc) { /* Declining fa, fb, fc */
+#ifdef DEBUG
+ printf("mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc);
+ fprintf(ficlog,"mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc);
+#endif
+ while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc */
r=(*bx-*ax)*(*fb-*fc);
q=(*bx-*cx)*(*fb-*fa);
u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
- (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscisse of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */
- ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */
- if ((*bx-u)*(u-*cx) > 0.0) { /* if u between b and c */
+ (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscissa of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */
+ ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscissa where function should be evaluated */
+ if ((*bx-u)*(u-*cx) > 0.0) { /* if u_p is between b and c */
fu=(*func)(u);
#ifdef DEBUG
/* f(x)=A(x-u)**2+f(u) */
fparabu= *fa - A*(*ax-u)*(*ax-u);
printf("mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);
fprintf(ficlog, "mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);
+ /* And thus,it can be that fu > *fc even if fparabu < *fc */
+ /* mnbrak (*ax=7.666299858533, *fa=299039.693133272231), (*bx=8.595447774979, *fb=298976.598289369489),
+ (*cx=10.098840694817, *fc=298946.631474258087), (*u=9.852501168332, fu=298948.773013752128, fparabu=298945.434711494134) */
+ /* In that case, there is no bracket in the output! Routine is wrong with many consequences.*/
#endif
+#ifdef MNBRAKORI
+#else
+ if (fu > *fc) {
+#ifdef DEBUG
+ printf("mnbrak4 fu > fc \n");
+ fprintf(ficlog, "mnbrak4 fu > fc\n");
+#endif
+ /* SHFT(u,*cx,*cx,u) /\* ie a=c, c=u and u=c; in that case, next SHFT(a,b,c,u) will give a=b=b, b=c=u, c=u=c and *\/ */
+ /* SHFT(*fa,*fc,fu,*fc) /\* (b, u, c) is a bracket while test fb > fc will be fu > fc will exit *\/ */
+ dum=u; /* Shifting c and u */
+ u = *cx;
+ *cx = dum;
+ dum = fu;
+ fu = *fc;
+ *fc =dum;
+ } else { /* end */
+#ifdef DEBUG
+ printf("mnbrak3 fu < fc \n");
+ fprintf(ficlog, "mnbrak3 fu < fc\n");
+#endif
+ dum=u; /* Shifting c and u */
+ u = *cx;
+ *cx = dum;
+ dum = fu;
+ fu = *fc;
+ *fc =dum;
+ }
+#endif
} else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */
+#ifdef DEBUG
+ printf("mnbrak2 u after c but before ulim\n");
+ fprintf(ficlog, "mnbrak2 u after c but before ulim\n");
+#endif
fu=(*func)(u);
if (fu < *fc) {
+#ifdef DEBUG
+ printf("mnbrak2 u after c but before ulim AND fu < fc\n");
+ fprintf(ficlog, "mnbrak2 u after c but before ulim AND fu <fc \n");
+#endif
SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
- SHFT(*fb,*fc,fu,(*func)(u))
- }
+ SHFT(*fb,*fc,fu,(*func)(u))
+ }
} else if ((u-ulim)*(ulim-*cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */
+#ifdef DEBUG
+ printf("mnbrak2 u outside ulim (verifying that ulim is beyond c)\n");
+ fprintf(ficlog, "mnbrak2 u outside ulim (verifying that ulim is beyond c)\n");
+#endif
u=ulim;
fu=(*func)(u);
- } else {
+ } else { /* u could be left to b (if r > q parabola has a maximum) */
+#ifdef DEBUG
+ printf("mnbrak2 u could be left to b (if r > q parabola has a maximum)\n");
+ fprintf(ficlog, "mnbrak2 u could be left to b (if r > q parabola has a maximum)\n");
+#endif
u=(*cx)+GOLD*(*cx-*bx);
fu=(*func)(u);
- }
+ } /* end tests */
SHFT(*ax,*bx,*cx,u)
- SHFT(*fa,*fb,*fc,fu)
- }
+ SHFT(*fa,*fb,*fc,fu)
+#ifdef DEBUG
+ printf("mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu);
+ fprintf(ficlog, "mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu);
+#endif
+ } /* end while; ie return (a, b, c, fa, fb, fc) such that a < b < c with f(a) > f(b) and fb < f(c) */
}
/*************** linmin ************************/
#endif
printf("%d",i);fflush(stdout);
fprintf(ficlog,"%d",i);fflush(ficlog);
- linmin(p,xit,n,fret,func);
+ linmin(p,xit,n,fret,func); /* xit[n] has been loaded for direction i */
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.
/* will reach at f3 = fm + h^2/2 f"m ; f" = (f1 -2f2 +f3 ) / h**2 */
/* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del */
/* 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); /* Intel compiler doesn't work on one line */
+#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=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */
t= t- del*SQR(fp-fptt);
+#endif
directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del 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);
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 POWELLDIRECT
+#ifdef POWELLORIGINAL
+ if (t < 0.0) { /* Then we use it for new direction */
+#else
if (directest*t < 0.0) { /* Contradiction between both tests */
printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt);
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,"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 */
-#else
- if (t < 0.0) { /* Then we use it for new direction */
#endif
linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/
for (j=1;j<=n;j++) {
which slows down the processing. The difference can be up to 10%
lower mortality.
*/
- lli=log(out[s1][s2] - savm[s1][s2]);
-
+ /* If, at the beginning of the maximization mostly, the
+ cumulative probability or probability to be dead is
+ constant (ie = 1) over time d, the difference is equal to
+ 0. out[s1][3] = savm[s1][3]: probability, being at state
+ s1 at precedent wave, to be dead a month before current
+ wave is equal to probability, being at state s1 at
+ precedent wave, to be dead at mont of the current
+ wave. Then the observed probability (that this person died)
+ is null according to current estimated parameter. In fact,
+ it should be very low but not zero otherwise the log go to
+ infinity.
+ */
+/* #ifdef INFINITYORIGINAL */
+/* lli=log(out[s1][s2] - savm[s1][s2]); */
+/* #else */
+/* if ((out[s1][s2] - savm[s1][s2]) < mytinydouble) */
+/* lli=log(mytinydouble); */
+/* else */
+/* lli=log(out[s1][s2] - savm[s1][s2]); */
+/* #endif */
+ lli=log(out[s1][s2] - savm[s1][s2]);
} else if (s2==-2) {
for (j=1,survp=0. ; j<=nlstate; j++)
ipmx +=1;
sw += weight[i];
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
+ /* if (lli < log(mytinydouble)){ */
+ /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */
+ /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */
+ /* } */
} /* end of wave */
} /* end of individual */
} else if(mle==2){
pstamp(ficrespij);
fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
i1= pow(2,cptcoveff);
- for(cptcov=1,k=0;cptcov<=i1;cptcov++){
- /*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
- k=k+1;
- /* for (k=1; k <= (int) pow(2,cptcoveff); k++){*/
- fprintf(ficrespij,"\n#****** ");
- for(j=1;j<=cptcoveff;j++)
- fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
- fprintf(ficrespij,"******\n");
+ /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
+ /* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
+ /* k=k+1; */
+ for (k=1; k <= (int) pow(2,cptcoveff); k++){
+ fprintf(ficrespij,"\n#****** ");
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
+ fprintf(ficrespij,"******\n");
+
+ for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
+ nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
+ nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
- for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
- nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
- nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
-
- /* nhstepm=nhstepm*YEARM; aff par mois*/
-
- p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
- oldm=oldms;savm=savms;
- hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
- fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
+ /* nhstepm=nhstepm*YEARM; aff par mois*/
+
+ p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ oldm=oldms;savm=savms;
+ hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
+ fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
+ for(i=1; i<=nlstate;i++)
+ for(j=1; j<=nlstate+ndeath;j++)
+ fprintf(ficrespij," %1d-%1d",i,j);
+ fprintf(ficrespij,"\n");
+ for (h=0; h<=nhstepm; h++){
+ /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
+ fprintf(ficrespij,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm );
for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate+ndeath;j++)
- fprintf(ficrespij," %1d-%1d",i,j);
- fprintf(ficrespij,"\n");
- for (h=0; h<=nhstepm; h++){
- /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
- fprintf(ficrespij,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm );
- for(i=1; i<=nlstate;i++)
- for(j=1; j<=nlstate+ndeath;j++)
- fprintf(ficrespij," %.5f", p3mat[i][j][h]);
- fprintf(ficrespij,"\n");
- }
- free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ fprintf(ficrespij," %.5f", p3mat[i][j][h]);
fprintf(ficrespij,"\n");
}
+ free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
+ fprintf(ficrespij,"\n");
+ }
/*}*/
}
}
git://henry.ined.fr / .git / commitdiff