[BACK]Return to imachprax.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imachprax.c between versions 1.2 and 1.4

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
Removed from v.1.2  
changed lines
  Added in v.1.4

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