--- imach/src/imachprax.c 2023/06/22 11:22:40 1.2 +++ imach/src/imachprax.c 2023/06/22 12:50:51 1.4 @@ -1,6 +1,12 @@ -/* $Id: imachprax.c,v 1.2 2023/06/22 11:22:40 brouard Exp $ +/* $Id: imachprax.c,v 1.4 2023/06/22 12:50:51 brouard Exp $ $State: Exp $ $Log: imachprax.c,v $ + 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 *** + Revision 1.2 2023/06/22 11:22:40 brouard Summary: with svd but not working yet @@ -1281,6 +1287,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 @@ -1373,12 +1380,12 @@ double gnuplotversion=GNUPLOTVERSION; #define ODIRSEPARATOR '\\' #endif -/* $Id: imachprax.c,v 1.2 2023/06/22 11:22:40 brouard Exp $ */ +/* $Id: imachprax.c,v 1.4 2023/06/22 12:50:51 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022"; -char fullversion[]="$Revision: 1.2 $ $Date: 2023/06/22 11:22:40 $"; +char fullversion[]="$Revision: 1.4 $ $Date: 2023/06/22 12:50:51 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -2613,6 +2620,41 @@ void linmin(double p[], double xi[], int } /**** praxis ****/ + +void transpose_in_place ( int n, double **a ) + +/******************************************************************************/ +/* + Purpose: + + TRANSPOSE_IN_PLACE transposes a square matrix in place. + Licensing: + This code is distributed under the GNU LGPL license. + + Input, int N, the number of rows and columns of the matrix A. + + Input/output, double A[N*N], the matrix to be transposed. +*/ +{ + int i; + int j; + double t; + + /* for ( j = 0; j < n; j++ ){ */ + /* for ( i = 0; i < j; i++ ) { */ + for ( j = 1; j <= n; j++ ){ + for ( i = 1; i < j; i++ ) { + /* t = a[i+j*n]; */ + /* a[i+j*n] = a[j+i*n]; */ + /* a[j+i*n] = t; */ + t = a[i][j]; + a[i][j] = a[j][i]; + a[j][i] = t; + } + } + return; +} + double pythag( double x, double y ) /******************************************************************************/ @@ -2653,9 +2695,16 @@ double pythag( double x, double y ) } /* void svdcmp(double **a, int m, int n, double w[], double **v) */ -void svdcmp(double **a, int m, int n, double w[]) +void svdminfit(double **a, int m, int n, double w[]) { - /* Golub 1970 Algol60 */ + /* From numerical recipes */ + /* Given a matrix a[1..m][1..n], this routine computes its singular value */ + /* decomposition, A = U ·W ·V T . The matrix U replaces a on output. */ + /* The diagonal matrix of singular values W is out- put as a vector w[1..n]. */ + /* The matrix V (not the transpose V T ) is output as v[1..n][1..n]. */ + + /* But in fact from Golub 1970 Algol60 */ + /* Computation of the singular values and complete orthogonal decom- */ /* position of a real rectangular matrix A, */ /* A = U diag (q) V^T, U^T U = V^T V =I , */ @@ -2803,7 +2852,7 @@ void svdcmp(double **a, int m, int n, do } break; } - if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations"); + if (its == 30) nrerror("no convergence in 30 svdcmp iterations"); x=w[l]; /* shift from bottom 2 x 2 minor; */ nm=k-1; y=w[nm]; @@ -2999,7 +3048,9 @@ void powell(double p[], double **xi, int 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 */ @@ -3029,9 +3080,8 @@ void powell(double p[], double **xi, int 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]); @@ -3169,10 +3219,37 @@ void powell(double p[], double **xi, int } } #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 */ + } +/* + 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. + + 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 LINMINORIGINAL #else for (j=1, flatd=0;j<=n;j++) { @@ -3201,6 +3278,15 @@ void powell(double p[], double **xi, int #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); @@ -4811,7 +4897,7 @@ double funcone( double *x) * 3 ncovta=15 +age*V3*V2+age*V2+agev3+ageV4 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 * 3 TvarAVVA[1]@15= itva 3 2 2 3 4 6 7 6 3 7 3 6 4 7 4 * 3 ncovta 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 - * + *?TvarAVVAind[1]@15= V3 is in k=2 1 1 2 3 4 5 4,2 5,2, 4,3 5 3}TvarVVAind[] * TvarAVVAind[1]@15= V3 is in k=6 6 12 13 14 15 16 18 18 19,19, 20,20 21,21}TvarVVAind[] * 3 ncovvta=10 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva] @@ -4825,7 +4911,8 @@ double funcone( double *x) * 2, 3, 4, 6, 7, * 6, 8, 9, 10, 11} * TvarFind[itv] 0 0 0 - * FixedV[itv] 1 1 1 0 1 0 1 0 1 0 1 0 1 0 + * FixedV[itv] 1 1 1 0 1 0 1 0 1 0 0 + *? FixedV[itv] 1 1 1 0 1 0 1 0 1 0 1 0 1 0 * Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4 * Tvar[TvarFind[itv]] [0]=? ?ncovv 1 à ncovvt] * Not a fixed cotvar[mw][itv][i] 6 7 6 2 7, 2, 6, 3, 7, 3, 6, 4, 7, 4} @@ -5210,7 +5297,7 @@ void mlikeli(FILE *ficres,double p[], in 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"); @@ -5279,16 +5366,16 @@ void mlikeli(FILE *ficres,double p[], in #else /* FLATSUP */ /* powell(p,xi,npar,ftol,&iter,&fret,func);*/ /* praxis ( t0, h0, n, prin, x, beale_f ); */ - int prin=4; - double h0=0.25; -#include "praxis.h" +/* int prin=4; */ +/* double h0=0.25; */ +/* #include "praxis.h" */ /* Be careful that praxis start at x[0] and powell start at p[1] */ /* praxis ( ftol, h0, npar, prin, p, func ); */ -p1= (p+1); /* p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */ -printf("Praxis \n"); -fprintf(ficlog, "Praxis \n");fflush(ficlog); -praxis ( ftol, h0, npar, prin, p1, func ); -printf("End Praxis\n"); +/* p1= (p+1); /\* p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] *\/ */ +/* printf("Praxis \n"); */ +/* fprintf(ficlog, "Praxis \n");fflush(ficlog); */ +/* praxis ( ftol, h0, npar, prin, p1, func ); */ +/* printf("End Praxis\n"); */ #endif /* FLATSUP */ #ifdef LINMINORIGINAL