imach/src/imach.c - diff

Return to imach.c CVS log

Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.356 and 1.357

version 1.356, 2023/05/23 12:08:43	version 1.357, 2023/06/14 14:55:52
Line 1	Line 1
/* $Id$	/* $Id$
$State$	$State$
$Log$	$Log$
	Revision 1.357 2023/06/14 14:55:52 brouard
	* imach.c (Module): Testing if conjugate gradient could be quicker when lot of variables POWELLORIGINCONJUGATE

Revision 1.356 2023/05/23 12:08:43 brouard	Revision 1.356 2023/05/23 12:08:43 brouard
Summary: 0.99r46	Summary: 0.99r46

Line 1292 Important routines	Line 1295 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 */

#include <math.h>	#include <math.h>
#include <stdio.h>	#include <stdio.h>
Line 2755 void powell(double p[], double **xi, int	Line 2759 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 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 */
Line 2785 void powell(double p[], double **xi, int	Line 2791 void powell(double p[], double **xi, int
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 2925 void powell(double p[], double **xi, int	Line 2930 void powell(double p[], double **xi, int
}	}
}	}
#endif	#endif
	#ifdef POWELLORIGINCONJUGATE
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][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 */
}	}
	#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	#ifdef LINMINORIGINAL
#else	#else
for (j=1, flatd=0;j<=n;j++) {	for (j=1, flatd=0;j<=n;j++) {
Line 2957 void powell(double p[], double **xi, int	Line 2969 void powell(double p[], double **xi, int
#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 */
	/* $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	#ifdef DEBUG
printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);	printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
Line 4967 void mlikeli(FILE *ficres,double p[], in	Line 4988 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\n"); fprintf(ficlog,"Powell\n");

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

Removed from v.1.356
changed lines
	Added in v.1.357