version 1.2, 2023/06/22 11:22:40
|
version 1.4, 2023/06/22 12:50:51
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$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 *** |
|
|
Revision 1.2 2023/06/22 11:22:40 brouard |
Revision 1.2 2023/06/22 11:22:40 brouard |
Summary: with svd but not working yet |
Summary: with svd but not working yet |
|
|
Line 1281 Important routines
|
Line 1287 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 2613 void linmin(double p[], double xi[], int
|
Line 2620 void linmin(double p[], double xi[], int
|
} |
} |
|
|
/**** praxis ****/ |
/**** 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 ) |
double pythag( double x, double y ) |
|
|
/******************************************************************************/ |
/******************************************************************************/ |
Line 2653 double pythag( double x, double y )
|
Line 2695 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[], 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- */ |
/* Computation of the singular values and complete orthogonal decom- */ |
/* position of a real rectangular matrix A, */ |
/* position of a real rectangular matrix A, */ |
/* A = U diag (q) V^T, U^T U = V^T V =I , */ |
/* A = U diag (q) V^T, U^T U = V^T V =I , */ |
Line 2803 void svdcmp(double **a, int m, int n, do
|
Line 2852 void svdcmp(double **a, int m, int n, do
|
} |
} |
break; |
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; */ |
x=w[l]; /* shift from bottom 2 x 2 minor; */ |
nm=k-1; |
nm=k-1; |
y=w[nm]; |
y=w[nm]; |
Line 2999 void powell(double p[], double **xi, int
|
Line 3048 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 3029 void powell(double p[], double **xi, int
|
Line 3080 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 3169 void powell(double p[], double **xi, int
|
Line 3219 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 */ |
|
} |
|
/* |
|
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 |
#ifdef LINMINORIGINAL |
#else |
#else |
for (j=1, flatd=0;j<=n;j++) { |
for (j=1, flatd=0;j<=n;j++) { |
Line 3201 void powell(double p[], double **xi, int
|
Line 3278 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 4811 double funcone( double *x)
|
Line 4897 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 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 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 |
* 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[] |
* 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 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] |
* 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva] |
Line 4825 double funcone( double *x)
|
Line 4911 double funcone( double *x)
|
* 2, 3, 4, 6, 7, |
* 2, 3, 4, 6, 7, |
* 6, 8, 9, 10, 11} |
* 6, 8, 9, 10, 11} |
* TvarFind[itv] 0 0 0 |
* 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[ncovf]]=[1]=2 [2]=3 [4]=4 |
* Tvar[TvarFind[itv]] [0]=? ?ncovv 1 à ncovvt] |
* 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} |
* Not a fixed cotvar[mw][itv][i] 6 7 6 2 7, 2, 6, 3, 7, 3, 6, 4, 7, 4} |
Line 5210 void mlikeli(FILE *ficres,double p[], in
|
Line 5297 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"); |
Line 5279 void mlikeli(FILE *ficres,double p[], in
|
Line 5366 void mlikeli(FILE *ficres,double p[], in
|
#else /* FLATSUP */ |
#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 ); */ |
/* praxis ( t0, h0, n, prin, x, beale_f ); */ |
int prin=4; |
/* int prin=4; */ |
double h0=0.25; |
/* double h0=0.25; */ |
#include "praxis.h" |
/* #include "praxis.h" */ |
/* Be careful that praxis start at x[0] and powell start at p[1] */ |
/* Be careful that praxis start at x[0] and powell start at p[1] */ |
/* praxis ( ftol, h0, npar, prin, p, func ); */ |
/* praxis ( ftol, h0, npar, prin, p, func ); */ |
p1= (p+1); /* p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */ |
/* p1= (p+1); /\* p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] *\/ */ |
printf("Praxis \n"); |
/* printf("Praxis \n"); */ |
fprintf(ficlog, "Praxis \n");fflush(ficlog); |
/* fprintf(ficlog, "Praxis \n");fflush(ficlog); */ |
praxis ( ftol, h0, npar, prin, p1, func ); |
/* praxis ( ftol, h0, npar, prin, p1, func ); */ |
printf("End Praxis\n"); |
/* printf("End Praxis\n"); */ |
#endif /* FLATSUP */ |
#endif /* FLATSUP */ |
|
|
#ifdef LINMINORIGINAL |
#ifdef LINMINORIGINAL |