From c138dd12b3565ff750795acbe407902955f093f9 Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Wed, 14 Jun 2023 14:55:52 +0000 Subject: [PATCH] * imach.c (Module): Testing if conjugate gradient could be quicker when lot of variables POWELLORIGINCONJUGATE --- src/imach.c | 31 +++++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/src/imach.c b/src/imach.c index 76a9ffb..ccfbedb 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,6 +1,11 @@ /* $Id$ $State$ $Log$ + Revision 1.356 2023/05/23 12:08:43 brouard + Summary: 0.99r46 + + * imach.c (Module): Fixed PROB_r + Revision 1.355 2023/05/22 17:03:18 brouard Summary: 0.99r46 @@ -1287,6 +1292,7 @@ Important routines /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */ /* #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 */ #include #include @@ -2750,7 +2756,9 @@ 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); /* 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 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 */ @@ -2780,9 +2788,8 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret fprintf(ficlog,"\n"); #endif } /* 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 */ - /* New value of last point Pn is not computed, P(n-1) */ for(j=1;j<=n;j++) { if(flatdir[j] >0){ printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); @@ -2920,10 +2927,17 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret } } #endif +#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 */ + xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ + } +#endif #ifdef LINMINORIGINAL #else for (j=1, flatd=0;j<=n;j++) { @@ -2952,6 +2966,15 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret #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); + /* The minimization in direction $\xi_1$ gives $P_1$. From $P_1$ minimization in direction $\xi_2$ gives */ + /* $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 */ + /* direction $\xi_1^{(1)}=\xi_2$. Also second direction is replaced by new direction $\xi^{(1)}_2=P_2-P_0$. */ + + /* At the second iteration, starting from $P_0^{(1)}$, minimization amongst $\xi^{(1)}_1$ gives point $P^{(1)}_1$. */ + /* 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); @@ -4962,7 +4985,7 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub 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++) xi[i][j]=(i==j ? 1.0 : 0.0); printf("Powell\n"); fprintf(ficlog,"Powell\n"); -- 2.43.0