version 1.3, 2023/06/22 11:28:07
|
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 |
Revision 1.3 2023/06/22 11:28:07 brouard |
*** empty log message *** |
*** empty log message *** |
|
|
Line 2617 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 2657 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 2807 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 3184 void powell(double p[], double **xi, int
|
Line 3229 void powell(double p[], double **xi, int
|
xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */ |
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 */ |
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 |
#endif |
#ifdef LINMINORIGINAL |
#ifdef LINMINORIGINAL |
#else |
#else |
Line 5301 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 |