/* $Id$
$State$
$Log$
+ Revision 1.5 2023/10/09 09:10:01 brouard
+ Summary: trying to reconsider
+
Revision 1.4 2023/06/22 12:50:51 brouard
Summary: stil on going
int maxwav=0; /* Maxim number of waves */
int jmin=0, jmax=0; /* min, max spacing between 2 waves */
int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */
-int gipmx=0, gsw=0; /* Global variables on the number of contributions
+int gipmx = 0;
+double gsw = 0; /* Global variables on the number of contributions
to the likelihood and the sum of weights (done by funcone)*/
int mle=1, weightopt=0;
int **mw; /* mw[mi][i] is number of the mi wave for this individual */
free_vector(pcom,1,n);
}
-/**** praxis ****/
-# include <float.h>
-
-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.
+/**** praxis gegen ****/
+
+/* This has been tested by Visual C from Microsoft and works */
+/* meaning tha valgrind could be wrong */
+/*********************************************************************/
+/* f u n c t i o n p r a x i s */
+/* */
+/* praxis is a general purpose routine for the minimization of a */
+/* function in several variables. the algorithm used is a modifi- */
+/* cation of conjugate gradient search method by powell. the changes */
+/* are due to r.p. brent, who gives an algol-w program, which served */
+/* as a basis for this function. */
+/* */
+/* references: */
+/* - powell, m.j.d., 1964. an efficient method for finding */
+/* the minimum of a function in several variables without */
+/* calculating derivatives, computer journal, 7, 155-162 */
+/* - brent, r.p., 1973. algorithms for minimization without */
+/* derivatives, prentice hall, englewood cliffs. */
+/* */
+/* problems, suggestions or improvements are always wellcome */
+/* karl gegenfurtner 07/08/87 */
+/* c - version */
+/*********************************************************************/
+/* */
+/* usage: min = praxis(tol, macheps, h, n, prin, x, func) */
+/* macheps has been suppressed because it is replaced by DBL_EPSILON */
+/* and if it was an argument of praxis (as it is in original brent) */
+/* it should be declared external */
+/* usage: min = praxis(tol, h, n, prin, x, func) */
+/* was min = praxis(fun, x, n); */
+/* */
+/* fun the function to be minimized. fun is called from */
+/* praxis with x and n as arguments */
+/* x a double array containing the initial guesses for */
+/* the minimum, which will contain the solution on */
+/* return */
+/* n an integer specifying the number of unknown */
+/* parameters */
+/* min praxis returns the least calculated value of fun */
+/* */
+/* some additional global variables control some more aspects of */
+/* the inner workings of praxis. setting them is optional, they */
+/* are all set to some reasonable default values given below. */
+/* */
+/* prin controls the printed output from the routine. */
+/* 0 -> no output */
+/* 1 -> print only starting and final values */
+/* 2 -> detailed map of the minimization process */
+/* 3 -> print also eigenvalues and vectors of the */
+/* search directions */
+/* the default value is 1 */
+/* tol is the tolerance allowed for the precision of the */
+/* solution. praxis returns if the criterion */
+/* 2 * ||x[k]-x[k-1]|| <= sqrt(macheps) * ||x[k]|| + tol */
+/* is fulfilled more than ktm times. */
+/* the default value depends on the machine precision */
+/* ktm see just above. default is 1, and a value of 4 leads */
+/* to a very(!) cautious stopping criterion. */
+/* h0 or step is a steplength parameter and should be set equal */
+/* to the expected distance from the solution. */
+/* exceptionally small or large values of step lead to */
+/* slower convergence on the first few iterations */
+/* the default value for step is 1.0 */
+/* scbd is a scaling parameter. 1.0 is the default and */
+/* indicates no scaling. if the scales for the different */
+/* parameters are very different, scbd should be set to */
+/* a value of about 10.0. */
+/* illc should be set to true (1) if the problem is known to */
+/* be ill-conditioned. the default is false (0). this */
+/* variable is automatically set, when praxis finds */
+/* the problem to be ill-conditioned during iterations. */
+/* maxfun is the maximum number of calls to fun allowed. praxis */
+/* will return after maxfun calls to fun even when the */
+/* minimum is not yet found. the default value of 0 */
+/* indicates no limit on the number of calls. */
+/* this return condition is only checked every n */
+/* iterations. */
+/* */
+/*********************************************************************/
- Input, int N, the number of rows and columns of the matrix A.
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h> /* for DBL_EPSILON */
+/* #include "machine.h" */
+
+
+/* extern void minfit(int n, double eps, double tol, double **ab, double q[]); */
+/* extern void minfit(int n, double eps, double tol, double ab[N][N], double q[]); */
+/* control parameters */
+/* control parameters */
+#define SQREPSILON 1.0e-19
+/* #define EPSILON 1.0e-8 */ /* in main */
+
+double tol = SQREPSILON,
+ scbd = 1.0,
+ step = 1.0;
+int ktm = 1,
+ /* prin = 2, */
+ maxfun = 0,
+ illc = 0;
+
+/* some global variables */
+static int i, j, k, k2, nl, nf, kl, kt;
+/* static double s; */
+double sl, dn, dmin,
+ fx, f1, lds, ldt, sf, df,
+ qf1, qd0, qd1, qa, qb, qc,
+ m2, m4, small_windows, vsmall, large,
+ vlarge, ldfac, t2;
+/* static double d[N], y[N], z[N], */
+/* q0[N], q1[N], v[N][N]; */
+
+static double *d, *y, *z;
+static double *q0, *q1, **v;
+double *tflin; /* used in flin: return (*fun)(tflin, n); */
+double *e; /* used in minfit, don't konw how to free memory and thus made global */
+/* static double s, sl, dn, dmin, */
+/* fx, f1, lds, ldt, sf, df, */
+/* qf1, qd0, qd1, qa, qb, qc, */
+/* m2, m4, small, vsmall, large, */
+/* vlarge, ldfac, t2; */
+/* static double d[N], y[N], z[N], */
+/* q0[N], q1[N], v[N][N]; */
+
+/* these will be set by praxis to point to it's arguments */
+static int prin; /* added */
+static int n;
+static double *x;
+static double (*fun)();
+/* static double (*fun)(double *x, int n); */
+
+/* these will be set by praxis to the global control parameters */
+/* static double h, macheps, t; */
+extern double macheps;
+static double h;
+static double t;
+
+static double
+drandom() /* return random no between 0 and 1 */
+{
+ return (double)(rand()%(8192*2))/(double)(8192*2);
+}
- Input/output, double A[N*N], the matrix to be transposed.
-*/
+static void sort() /* d and v in descending order */
{
- int i;
- int j;
- double t;
+ int k, i, j;
+ double s;
+
+ for (i=1; i<=n-1; i++) {
+ k = i; s = d[i];
+ for (j=i+1; j<=n; j++) {
+ if (d[j] > s) {
+ k = j;
+ s = d[j];
+ }
+ }
+ if (k > i) {
+ d[k] = d[i];
+ d[i] = s;
+ for (j=1; j<=n; j++) {
+ s = v[j][i];
+ v[j][i] = v[j][k];
+ v[j][k] = s;
+ }
+ }
+ }
+}
- /* 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;
+double randbrent ( int *naught )
+{
+ double ran1, ran3[127], half;
+ int ran2, q, r, i, j;
+ int init=0; /* false */
+ double rr;
+ /* REAL*8 RAN1,RAN3(127),HALF */
+
+ /* INTEGER RAN2,Q,R */
+ /* LOGICAL INIT */
+ /* DATA INIT/.FALSE./ */
+ /* IF (INIT) GO TO 3 */
+ if(!init){
+/* R = MOD(NAUGHT,8190) + 1 *//* 1804289383 rand () */
+ r = *naught % 8190 + 1;/* printf(" naught r %d %d",*naught,r); */
+ ran2=127;
+ for(i=ran2; i>0; i--){
+/* RAN2 = 128 */
+/* DO 2 I=1,127 */
+ ran2 = ran2-1;
+/* RAN2 = RAN2 - 1 */
+ ran1 = -pow(2.0,55);
+/* RAN1 = -2.D0**55 */
+/* DO 1 J=1,7 */
+ for(j=1; j<=7;j++){
+/* R = MOD(1756*R,8191) */
+ r = (1756*r) % 8191;/* printf(" i=%d (1756*r)%8191=%d",j,r); */
+ q=r/32;
+/* Q = R/32 */
+/* 1 RAN1 = (RAN1 + Q)*(1.0D0/256) */
+ ran1 =(ran1+q)*(1.0/256);
+ }
+/* 2 RAN3(RAN2) = RAN1 */
+ ran3[ran2] = ran1; /* printf(" ran2=%d ran1=%.7g \n",ran2,ran1); */
}
+/* INIT = .TRUE. */
+ init=1;
+/* 3 IF (RAN2.EQ.1) RAN2 = 128 */
}
- return;
+ if(ran2 == 0) ran2 = 126;
+ else ran2 = ran2 -1;
+ /* RAN2 = RAN2 - 1 */
+ /* RAN1 = RAN1 + RAN3(RAN2) */
+ ran1 = ran1 + ran3[ran2];/* printf("BIS ran2=%d ran1=%.7g \n",ran2,ran1); */
+ half= 0.5;
+ /* HALF = .5D0 */
+ /* IF (RAN1.GE.0.D0) HALF = -HALF */
+ if(ran1 >= 0.) half =-half;
+ ran1 = ran1 +half;
+ ran3[ran2] = ran1;
+ rr= ran1+0.5;
+ /* RAN1 = RAN1 + HALF */
+ /* RAN3(RAN2) = RAN1 */
+ /* RANDOM = RAN1 + .5D0 */
+/* r = ( ( double ) ( *seed ) ) * 4.656612875E-10; */
+ return rr;
}
+static void matprint(char *s, double **v, int m, int n)
+/* char *s; */
+/* double v[N][N]; */
+{
+#define INCX 8
+ int i;
+
+ int i2hi;
+ int ihi;
+ int ilo;
+ int i2lo;
+ int jlo=1;
+ int j;
+ int j2hi;
+ int jhi;
+ int j2lo;
+ ilo=1;
+ ihi=n;
+ jlo=1;
+ jhi=n;
+
+ printf ("\n" );
+ printf ("%s\n", s );
+ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
+ {
+ j2hi = j2lo + INCX - 1;
+ if ( n < j2hi )
+ {
+ j2hi = n;
+ }
+ if ( jhi < j2hi )
+ {
+ j2hi = jhi;
+ }
-double pythag( double x, double y )
+ /* fprintf ( ficlog, "\n" ); */
+ printf ("\n" );
+/*
+ For each column J in the current range...
-/******************************************************************************/
+ Write the header.
+*/
+ /* fprintf ( ficlog, " Col: "); */
+ printf ("Col:");
+ for ( j = j2lo; j <= j2hi; j++ )
+ {
+ /* fprintf ( ficlog, " %7d ", j - 1 ); */
+ /* printf (" %9d ", j - 1 ); */
+ printf (" %9d ", j );
+ }
+ /* fprintf ( ficlog, "\n" ); */
+ /* fprintf ( ficlog, " Row\n" ); */
+ /* fprintf ( ficlog, "\n" ); */
+ printf ("\n" );
+ printf (" Row\n" );
+ printf ("\n" );
/*
- Purpose:
- R8_HYPOT returns the value of sqrt ( X^2 + Y^2 ).
- Licensing:
- This code is distributed under the GNU LGPL license.
- Modified:
- 26 March 2012
- Author:
- John Burkardt
- Parameters:
- Input, double X, Y, the arguments.
- Output, double R8_HYPOT, the value of sqrt ( X^2 + Y^2 ).
+ Determine the range of the rows in this strip.
*/
-{
- double a;
- double b;
- double value;
+ if ( 1 < ilo ){
+ i2lo = ilo;
+ }else{
+ i2lo = 1;
+ }
+ if ( m < ihi ){
+ i2hi = m;
+ }else{
+ i2hi = ihi;
+ }
- if ( fabs ( x ) < fabs ( y ) ) {
- a = fabs ( y );
- b = fabs ( x );
- } else {
- a = fabs ( x );
- b = fabs ( y );
- }
+ for ( i = i2lo; i <= i2hi; i++ ){
/*
- A contains the larger value.
+ Print out (up to) 5 entries in row I, that lie in the current strip.
*/
- if ( a == 0.0 ) {
- value = 0.0;
- } else {
- value = a * sqrt ( 1.0 + ( b / a ) * ( b / a ) );
+ /* fprintf ( ficlog, "%5d:", i - 1 ); */
+ /* printf ("%5d:", i - 1 ); */
+ printf ("%5d:", i );
+ for ( j = j2lo; j <= j2hi; j++ )
+ {
+ /* fprintf ( ficlog, " %14g", a[i-1+(j-1)*m] ); */
+ /* printf ("%14.7g ", a[i-1+(j-1)*m] ); */
+ /* printf("%14.7f ", v[i-1][j-1]); */
+ printf("%14.7f ", v[i][j]);
+ /* fprintf ( stdout, " %14g", a[i-1+(j-1)*m] ); */
+ }
+ /* fprintf ( ficlog, "\n" ); */
+ printf ("\n" );
+ }
}
- return value;
+
+ /* printf("%s\n", s); */
+ /* for (k=0; k<n; k++) { */
+ /* for (i=0; i<n; i++) { */
+ /* /\* printf("%20.10e ", v[k][i]); *\/ */
+ /* } */
+ /* printf("\n"); */
+ /* } */
+#undef INCX
}
-void svsort ( int n, double d[], double **v )
+void vecprint(char *s, double *x, int n)
+/* char *s; */
+/* double x[N]; */
+{
+ int i=0;
+
+ printf(" %s", s);
+ /* for (i=0; i<n; i++) */
+ for (i=1; i<=n; i++)
+ printf (" %14.7g", x[i] );
+ /* printf(" %8d: %14g\n", i, x[i]); */
+ printf ("\n" );
+}
-/******************************************************************************/
-/*
- Purpose:
+static void print() /* print a line of traces */
+{
+
+
+ printf("\n");
+ /* printf("... chi square reduced to ... %20.10e\n", fx); */
+ /* printf("... after %u function calls ...\n", nf); */
+ /* printf("... including %u linear searches ...\n", nl); */
+ printf("%10d %10d%14.7g",nl, nf, fx);
+ vecprint("... current values of x ...", x, n);
+}
+/* static void print2(int n, double *x, int prin, double fx, int nf, int nl) */ /* print a line of traces */
+static void print2() /* print a line of traces */
+{
+ int i; double fmin=0.;
+
+ /* printf("\n"); */
+ /* printf("... chi square reduced to ... %20.10e\n", fx); */
+ /* printf("... after %u function calls ...\n", nf); */
+ /* printf("... including %u linear searches ...\n", nl); */
+ /* printf("%10d %10d%14.7g",nl, nf, fx); */
+ printf ( "\n" );
+ printf ( " Linear searches %d", nl );
+ /* printf ( " Linear searches %d\n", nl ); */
+ /* printf ( " Function evaluations %d\n", nf ); */
+ /* printf ( " Function value FX = %g\n", fx ); */
+ printf ( " Function evaluations %d", nf );
+ printf ( " Function value FX = %.12lf\n", fx );
+#ifdef DEBUGPRAX
+ printf("n=%d prin=%d\n",n,prin);
+#endif
+ if(fx <= fmin) printf(" UNDEFINED "); else printf("%14.7g",log(fx-fmin));
+ if ( n <= 4 || 2 < prin )
+ {
+ /* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */
+ for(i=1;i<=n;i++)printf("%14.7g",x[i]);
+ /* r8vec_print ( n, x, " X:" ); */
+ }
+ printf("\n");
+ }
- SVSORT descending sorts D and adjusts the corresponding columns of V.
- Discussion:
- A simple bubble sort is used on D.
- In our application, D contains singular values, and the columns of V are
- the corresponding right singular vectors.
- Author:
- Original FORTRAN77 version by Richard Brent.
- Richard Brent,
- Algorithms for Minimization with Derivatives,
- Prentice Hall, 1973,
- Reprinted by Dover, 2002.
-
- Parameters:
- Input, int N, the length of D, and the order of V.
- Input/output, double D[N], the vector to be sorted.
- On output, the entries of D are in descending order.
-
- Input/output, double V[N,N], an N by N array to be adjusted
- as D is sorted. In particular, if the value that was in D(I) on input is
- moved to D(J) on output, then the input column V(*,I) is moved to
- the output column V(*,J).
-*/
+/* #ifdef MSDOS */
+/* static double tflin[N]; */
+/* #endif */
+
+static double flin(double l, int j)
+/* double l; */
{
- int i, j1, j2, j3;
- double t;
+ int i;
+ /* #ifndef MSDOS */
+ /* double tflin[N]; */
+ /* #endif */
+ /* double *tflin; */ /* Be careful to put tflin on a vector n */
+
+ /* j is used from 0 to n-1 and can be -1 for parabolic search */
+
+ /* if (j != -1) { /\* linear search *\/ */
+ if (j > 0) { /* linear search */
+ /* for (i=0; i<n; i++){ */
+ for (i=1; i<=n; i++){
+ tflin[i] = x[i] + l *v[i][j];
+#ifdef DEBUGPRAX
+ /* printf(" flin i=%14d t=%14.7f x=%14.7f l=%14.7f v[%d,%d]=%14.7f nf=%14d\n",i+1, tflin[i],x[i],l,i,j,v[i][j],nf); */
+ printf(" flin i=%14d t=%14.7f x=%14.7f l=%14.7f v[%d,%d]=%14.7f nf=%14d\n",i, tflin[i],x[i],l,i,j,v[i][j],nf);
+#endif
+ }
+ }
+ else { /* search along parabolic space curve */
+ qa = l*(l-qd1)/(qd0*(qd0+qd1));
+ qb = (l+qd0)*(qd1-l)/(qd0*qd1);
+ qc = l*(l+qd0)/(qd1*(qd0+qd1));
+#ifdef DEBUGPRAX
+ printf(" search along a parabolic space curve. j=%14d nf=%14d l=%14.7f qd0=%14.7f qd1=%14.7f\n",j,nf,l,qd0,qd1);
+#endif
+ /* for (i=0; i<n; i++){ */
+ for (i=1; i<=n; i++){
+ tflin[i] = qa*q0[i]+qb*x[i]+qc*q1[i];
+#ifdef DEBUGPRAX
+ /* printf(" parabole i=%14d t(i)=%14.7f q0=%14.7f x=%14.7f q1=%14.7f\n",i+1,tflin[i],q0[i],x[i],q1[i]); */
+ printf(" parabole i=%14d t(i)=%14.7e q0=%14.7e x=%14.7e q1=%14.7e\n",i,tflin[i],q0[i],x[i],q1[i]);
+#endif
+ }
+ }
+ nf++;
- for (j1 = 1; j1 < n; j1++) {
- /*
- * Find J3, the index of the largest entry in D[J1:N-1].
- * MAXLOC apparently requires its output to be an array.
- */
- j3 = j1;
- for (j2 = j1+1; j2 < n; j2++) {
- if (d[j3] < d[j2]) {
- j3 = j2;
+#ifdef NR_SHIFT
+ return (*fun)((tflin-1), n);
+#else
+ /* return (*fun)(tflin, n);*/
+ return (*fun)(tflin);
+#endif
+}
+
+void minny(int j, int nits, double *d2, double *x1, double f1, int fk)
+/* double *d2, *x1, f1; */
+{
+/* here j is from 0 to n-1 and can be -1 for parabolic search */
+ /* MINIMIZES F FROM X IN THE DIRECTION V(*,J) */
+ /* UNLESS J<1, WHEN A QUADRATIC SEARCH IS DONE */
+ /* IN THE PLANE DEFINED BY Q0, Q1 AND X. */
+ /* D2 AN APPROXIMATION TO HALF F'' (OR ZERO), */
+ /* X1 AN ESTIMATE OF DISTANCE TO MINIMUM, */
+ /* RETURNED AS THE DISTANCE FOUND. */
+ /* IF FK = TRUE THEN F1 IS FLIN(X1), OTHERWISE */
+ /* X1 AND F1 ARE IGNORED ON ENTRY UNLESS FINAL */
+ /* FX > F1. NITS CONTROLS THE NUMBER OF TIMES */
+ /* AN ATTEMPT IS MADE TO HALVE THE INTERVAL. */
+ /* SIDE EFFECTS: USES AND ALTERS X, FX, NF, NL. */
+ /* IF J < 1 USES VARIABLES Q... . */
+ /* USES H, N, T, M2, M4, LDT, DMIN, MACHEPS; */
+ int k, i, dz;
+ double x2, xm, f0, f2, fm, d1, t2, sf1, sx1;
+ double s;
+ double macheps;
+ macheps=pow(16.0,-13.0);
+ sf1 = f1; sx1 = *x1;
+ k = 0; xm = 0.0; fm = f0 = fx; dz = *d2 < macheps;
+ /* h=1.0;*/ /* To be revised */
+#ifdef DEBUGPRAX
+ /* printf("min macheps=%14g h=%14g step=%14g t=%14g fx=%14g\n",macheps,h, step,t, fx); */
+ /* Where is fx coming from */
+ printf(" min macheps=%14g h=%14g t=%14g fx=%.9lf dirj=%d\n",macheps, h, t, fx, j);
+ matprint(" min vectors:",v,n,n);
+#endif
+ /* find step size */
+ s = 0.;
+ /* for (i=0; i<n; i++) s += x[i]*x[i]; */
+ for (i=1; i<=n; i++) s += x[i]*x[i];
+ s = sqrt(s);
+ if (dz)
+ t2 = m4*sqrt(fabs(fx)/dmin + s*ldt) + m2*ldt;
+ else
+ t2 = m4*sqrt(fabs(fx)/(*d2) + s*ldt) + m2*ldt;
+ s = s*m4 + t;
+ if (dz && t2 > s) t2 = s;
+ if (t2 < small_windows) t2 = small_windows;
+ if (t2 > 0.01*h) t2 = 0.01 * h;
+ if (fk && f1 <= fm) {
+ xm = *x1;
+ fm = f1;
+ }
+#ifdef DEBUGPRAX
+ printf(" additional flin X1=%14.7f t2=%14.7f *f1=%14.7f fm=%14.7f fk=%d\n",*x1,t2,f1,fm,fk);
+#endif
+ if (!fk || fabs(*x1) < t2) {
+ *x1 = (*x1 >= 0 ? t2 : -t2);
+ /* *x1 = (*x1 > 0 ? t2 : -t2); */ /* kind of error */
+#ifdef DEBUGPRAX
+ printf(" additional flin X1=%16.10e dirj=%d fk=%d\n",*x1, j, fk);
+#endif
+ f1 = flin(*x1, j);
+#ifdef DEBUGPRAX
+ printf(" after flin f1=%18.12e dirj=%d fk=%d\n",f1, j,fk);
+#endif
+ }
+ if (f1 <= fm) {
+ xm = *x1;
+ fm = f1;
+ }
+L0: /*L0 loop or next */
+/*
+ Evaluate FLIN at another point and estimate the second derivative.
+*/
+ if (dz) {
+ x2 = (f0 < f1 ? -(*x1) : 2*(*x1));
+#ifdef DEBUGPRAX
+ printf(" additional second flin x2=%14.8e x1=%14.8e f0=%14.8e f1=%18.12e dirj=%d\n",x2,*x1,f0,f1,j);
+#endif
+ f2 = flin(x2, j);
+#ifdef DEBUGPRAX
+ printf(" additional second flin x2=%16.10e x1=%16.10e f1=%18.12e f0=%18.10e f2=%18.10e fm=%18.10e\n",x2, *x1, f1,f0,f2,fm);
+#endif
+ if (f2 <= fm) {
+ xm = x2;
+ fm = f2;
}
- }
- /*
- * If J1 != J3, swap D[J1] and D[J3], and columns J1 and J3 of V.
+ /* d2 is the curvature or double difference f1 doesn't seem to be accurately computed */
+ *d2 = (x2*(f1-f0) - (*x1)*(f2-f0))/((*x1)*x2*((*x1)-x2));
+#ifdef DEBUGPRAX
+ double d11,d12;
+ d11=(f1-f0)/(*x1);d12=(f2-f0)/x2;
+ printf(" d11=%18.12e d12=%18.12e d11-d12=%18.12e x1-x2=%18.12e (d11-d12)/(x2-(*x1))=%18.12e\n", d11 ,d12, d11-d12, x2-(*x1), (d11-d12)/(x2-(*x1)));
+ printf(" original computing f1=%18.12e *d2=%16.10e f0=%18.12e f1-f0=%16.10e f2-f0=%16.10e\n",f1,*d2,f0,f1-f0, f2-f0);
+ double ff1=7.783920622852e+04;
+ double f1mf0=9.0344736236e-05;
+ *d2 = (f1mf0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2);
+ /* *d2 = (ff1-f0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); */
+ printf(" simpliff computing *d2=%16.10e f1mf0=%18.12e,f1=f0+f1mf0=%18.12e\n",*d2,f1mf0,f0+f1mf0);
+ *d2 = ((f1-f0)/ (*x1) - (f2-f0)/x2)/((*x1)-x2);
+ printf(" overlifi computing *d2=%16.10e\n",*d2);
+#endif
+ *d2 = ((f1-f0)/ (*x1) - (f2-f0)/x2)/((*x1)-x2);
+ }
+#ifdef DEBUGPRAX
+ printf(" additional second flin xm=%14.8e fm=%14.8e *d2=%14.8e\n",xm, fm,*d2);
+#endif
+ /*
+ Estimate the first derivative at 0.
+ */
+ d1 = (f1-f0)/(*x1) - *x1**d2; dz = 1;
+ /*
+ Predict the minimum.
*/
- if (j1 != j3) {
- t = d[j1];
- d[j1] = d[j3];
- d[j3] = t;
- for (i = 1; i <= n; i++) {
- t = v[i][j1];
- v[i][j1] = v[i][j3];
- v[i][j3] = t;
- } /* end i */
- } /* end j1 != j3 */
- } /* end j1 */
- return;
+ if (*d2 <= small_windows) {
+ x2 = (d1 < 0 ? h : -h);
+ }
+ else {
+ x2 = - 0.5*d1/(*d2);
+ }
+#ifdef DEBUGPRAX
+ printf(" AT d1=%14.8e d2=%14.8e small=%14.8e dz=%d x1=%14.8e x2=%14.8e\n",d1,*d2,small_windows,dz,*x1,x2);
+#endif
+ if (fabs(x2) > h)
+ x2 = (x2 > 0 ? h : -h);
+L1: /* L1 or try loop */
+#ifdef DEBUGPRAX
+ printf(" AT predicted minimum flin x2=%14.8e x1=%14.8e K=%14d NITS=%14d dirj=%d\n",x2,*x1,k,nits,j);
+#endif
+ f2 = flin(x2, j); /* x[i]+x2*v[i][j] */
+#ifdef DEBUGPRAX
+ printf(" after flin f0=%14.8e f1=%14.8e f2=%14.8e fm=%14.8e\n",f0,f1,f2, fm);
+#endif
+ if ((k < nits) && (f2 > f0)) {
+#ifdef DEBUGPRAX
+ printf(" NO SUCCESS SO TRY AGAIN;\n");
+#endif
+ k++;
+ if ((f0 < f1) && (*x1*x2 > 0.0))
+ goto L0; /* or next */
+ x2 *= 0.5;
+ goto L1;
+ }
+ nl++;
+#ifdef DEBUGPRAX
+ printf(" bebeBE end of min x1=%14.8e x2=%14.8e f1=%14.8e f2=%14.8e f0=%14.8e fm=%14.8e d2=%14.8e\n",*x1, x2, f1, f2, f0, fm, *d2);
+#endif
+ if (f2 > fm) x2 = xm; else fm = f2;
+ if (fabs(x2*(x2-*x1)) > small_windows) {
+ *d2 = (x2*(f1-f0) - *x1*(fm-f0))/(*x1*x2*(*x1-x2));
+ }
+ else {
+ if (k > 0) *d2 = 0;
+ }
+#ifdef DEBUGPRAX
+ printf(" bebe end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
+#endif
+ if (*d2 <= small_windows) *d2 = small_windows;
+ *x1 = x2; fx = fm;
+ if (sf1 < fx) {
+ fx = sf1;
+ *x1 = sx1;
+ }
+ /*
+ Update X for linear search.
+ */
+#ifdef DEBUGPRAX
+ printf(" end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
+#endif
+
+ /* if (j != -1) */
+ /* for (i=0; i<n; i++) */
+ /* x[i] += (*x1)*v[i][j]; */
+ if (j > 0)
+ for (i=1; i<=n; i++)
+ x[i] += (*x1)*v[i][j];
}
-/* void svdcmp(double **a, int m, int n, double w[], double **v) */
-void svdminfit(double **a, int m, int n, double w[])
+void quad() /* look for a minimum along the curve q0, q1, q2 */
{
- /* 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 */
+ int i;
+ double l, s;
+
+ s = fx; fx = qf1; qf1 = s; qd1 = 0.0;
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ s = x[i]; l = q1[i]; x[i] = l; q1[i] = s;
+ qd1 = qd1 + (s-l)*(s-l);
+ }
+ s = 0.0; qd1 = sqrt(qd1); l = qd1;
+#ifdef DEBUGPRAX
+ printf(" QUAD after sqrt qd1=%14.8e \n",qd1);
+#endif
+
+ if (qd0>0.0 && qd1>0.0 &&nl>=3*n*n) {
+#ifdef DEBUGPRAX
+ printf(" QUAD before min value=%14.8e \n",qf1);
+#endif
+ /* min(-1, 2, &s, &l, qf1, 1); */
+ minny(0, 2, &s, &l, qf1, 1);
+ qa = l*(l-qd1)/(qd0*(qd0+qd1));
+ qb = (l+qd0)*(qd1-l)/(qd0*qd1);
+ qc = l*(l+qd0)/(qd1*(qd0+qd1));
+ }
+ else {
+ fx = qf1; qa = qb = 0.0; qc = 1.0;
+ }
+#ifdef DEBUGPRAX
+ printf("after eventual min qd0=%14.8e qd1=%14.8e nl=%d\n",qd0, qd1,nl);
+#endif
+ qd0 = qd1;
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ s = q0[i]; q0[i] = x[i];
+ x[i] = qa*s + qb*x[i] + qc*q1[i];
+ }
+#ifdef DEBUGQUAD
+ vecprint ( " X after QUAD:" , x, n );
+#endif
+}
-/* 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 , */
-/* where the arrays a[1:m, 1:n], u[1:m, 1 :n], v[1:n, 1:n], q[1:n] re- */
-/* present A, U, V, q respectively. The actual parameters corresponding */
-/* to a, u, v may all be identical unless withu=withv = true . In this */
-/* case, the actual parameters corresponding to u and v must differ. */
-/* m >= n is assumed; */
-
- /* Simplified (as in praxis) in order to output only V (replacing A), w is the diagonal matrix q */
+/* void minfit(int n, double eps, double tol, double ab[N][N], double q[]) */
+void minfit(int n, double eps, double tol, double **ab, double q[])
+/* int n; */
+/* double eps, tol, ab[N][N], q[N]; */
+{
+ int l, kt, l2, i, j, k;
+ double c, f, g, h, s, x, y, z;
+ /* double eps; */
+/* #ifndef MSDOS */
+/* double e[N]; /\* plenty of stack on a vax *\/ */
+/* #endif */
+ /* double *e; */
+ /* e=vector(0,n-1); /\* should be freed somewhere but gotos *\/ */
+
+ /* householder's reduction to bidiagonal form */
+
+ if(n==1){
+ /* q[1-1]=ab[1-1][1-1]; */
+ /* ab[1-1][1-1]=1.0; */
+ q[1]=ab[1][1];
+ ab[1][1]=1.0;
+ return; /* added from hardt */
+ }
+ /* eps=macheps; */ /* added */
+ x = g = 0.0;
+#ifdef DEBUGPRAX
+ matprint (" HOUSE holder:", ab, n, n);
+#endif
+
+ /* for (i=0; i<n; i++) { /\* FOR I := 1 UNTIL N DO *\/ */
+ for (i=1; i<=n; i++) { /* FOR I := 1 UNTIL N DO */
+ e[i] = g; s = 0.0; l = i+1;
+ /* for (j=i; j<n; j++) /\* FOR J := I UNTIL N DO S := S*AB(J,I)**2; *\/ /\* not correct *\/ */
+ for (j=i; j<=n; j++) /* FOR J := I UNTIL N DO S := S*AB(J,I)**2; */ /* not correct */
+ s += ab[j][i] * ab[j][i];
+#ifdef DEBUGPRAXFIN
+ printf("i=%d s=%d %.7g tol=%.7g",i,s,tol);
+#endif
+ if (s < tol) {
+ g = 0.0;
+ }
+ else {
+ /* f = ab[i][i]; */
+ f = ab[i][i];
+ if (f < 0.0)
+ g = sqrt(s);
+ else
+ g = -sqrt(s);
+ /* h = f*g - s; ab[i][i] = f - g; */
+ h = f*g - s; ab[i][i] = f - g;
+ /* for (j=l; j<n; j++) { */ /* FOR J := L UNTIL N DO */ /* wrong */
+ for (j=l; j<=n; j++) {
+ f = 0.0;
+ /* for (k=i; k<n; k++) /\* FOR K := I UNTIL N DO *\/ /\* wrong *\/ */
+ for (k=i; k<=n; k++) /* FOR K := I UNTIL N DO */
+ /* f += ab[k][i] * ab[k][j]; */
+ f += ab[k][i] * ab[k][j];
+ f /= h;
+ for (k=i; k<=n; k++) /* FOR K := I UNTIL N DO */
+ /* for (k=i; k<n; k++)/\* FOR K := I UNTIL N DO *\/ /\* wrong *\/ */
+ ab[k][j] += f * ab[k][i];
+ /* ab[k][j] += f * ab[k][i]; */
+#ifdef DEBUGPRAX
+ printf("Holder J=%d F=%.7g",j,f);
+#endif
+ }
+ } /* end s */
+ /* q[i] = g; s = 0.0; */
+ q[i] = g; s = 0.0;
+#ifdef DEBUGPRAX
+ printf(" I Q=%d %.7g",i,q[i]);
+#endif
+
+ /* if (i < n) */
+ /* if (i <= n) /\* I is always lower or equal to n wasn't in golub reinsch*\/ */
+ /* for (j=l; j<n; j++) */
+ for (j=l; j<=n; j++)
+ s += ab[i][j] * ab[i][j];
+ /* s += ab[i][j] * ab[i][j]; */
+ if (s < tol) {
+ g = 0.0;
+ }
+ else {
+ if(i<n)
+ /* f = ab[i][i+1]; */ /* Brent golub overflow */
+ f = ab[i][i+1];
+ if (f < 0.0)
+ g = sqrt(s);
+ else
+ g = - sqrt(s);
+ h = f*g - s;
+ /* h = f*g - s; ab[i][i+1] = f - g; */ /* Overflow for i=n Error in Golub too but not Burkardt*/
+ /* for (j=l; j<n; j++) */
+ /* e[j] = ab[i][j]/h; */
+ if(i<n){
+ ab[i][i+1] = f - g;
+ for (j=l; j<=n; j++)
+ e[j] = ab[i][j]/h;
+ /* for (j=l; j<n; j++) { */
+ for (j=l; j<=n; j++) {
+ s = 0.0;
+ /* for (k=l; k<n; k++) s += ab[j][k]*ab[i][k]; */
+ for (k=l; k<=n; k++) s += ab[j][k]*ab[i][k];
+ /* for (k=l; k<n; k++) ab[j][k] += s * e[k]; */
+ for (k=l; k<=n; k++) ab[j][k] += s * e[k];
+ } /* END J */
+ } /* END i <n */
+ } /* end s */
+ /* y = fabs(q[i]) + fabs(e[i]); */
+ y = fabs(q[i]) + fabs(e[i]);
+ if (y > x) x = y;
+#ifdef DEBUGPRAX
+ printf(" I Y=%d %.7g",i,y);
+#endif
+#ifdef DEBUGPRAX
+ printf(" i=%d e(i) %.7g",i,e[i]);
+#endif
+ } /* end i */
+ /*
+ Accumulation of right hand transformations */
+ /* for (i=n-1; i >= 0; i--) { */ /* FOR I := N STEP -1 UNTIL 1 DO */
+ /* We should avoid the overflow in Golub */
+ /* ab[n-1][n-1] = 1.0; */
+ /* g = e[n-1]; */
+ ab[n][n] = 1.0;
+ g = e[n];
+ l = n;
+
+ /* for (i=n; i >= 1; i--) { */
+ for (i=n-1; i >= 1; i--) { /* n-1 loops, different from brent and golub*/
+ if (g != 0.0) {
+ /* h = ab[i-1][i]*g; */
+ h = ab[i][i+1]*g;
+ for (j=l; j<=n; j++) ab[j][i] = ab[i][j] / h;
+ for (j=l; j<=n; j++) {
+ /* h = ab[i][i+1]*g; */
+ /* for (j=l; j<n; j++) ab[j][i] = ab[i][j] / h; */
+ /* for (j=l; j<n; j++) { */
+ s = 0.0;
+ /* for (k=l; k<n; k++) s += ab[i][k] * ab[k][j]; */
+ /* for (k=l; k<n; k++) ab[k][j] += s * ab[k][i]; */
+ for (k=l; k<=n; k++) s += ab[i][k] * ab[k][j];
+ for (k=l; k<=n; k++) ab[k][j] += s * ab[k][i];
+ }/* END J */
+ }/* END G */
+ /* for (j=l; j<n; j++) */
+ /* ab[i][j] = ab[j][i] = 0.0; */
+ /* ab[i][i] = 1.0; g = e[i]; l = i; */
+ for (j=l; j<=n; j++)
+ ab[i][j] = ab[j][i] = 0.0;
+ ab[i][i] = 1.0; g = e[i]; l = i;
+ }/* END I */
+#ifdef DEBUGPRAX
+ matprint (" HOUSE accumulation:",ab,n, n );
+#endif
+
+ /* diagonalization to bidiagonal form */
+ eps *= x;
+ /* for (k=n-1; k>= 0; k--) { */
+ for (k=n; k>= 1; k--) {
+ kt = 0;
+TestFsplitting:
+#ifdef DEBUGPRAX
+ printf(" TestFsplitting: k=%d kt=%d\n",k,kt);
+ /* for(i=1;i<=n;i++)printf(" e(%d)=%.14f",i,e[i]);printf("\n"); */
+#endif
+ kt = kt+1;
+/* TestFsplitting: */
+ /* if (++kt > 30) { */
+ if (kt > 30) {
+ e[k] = 0.0;
+ fprintf(stderr, "\n+++ MINFIT - Fatal error\n");
+ fprintf ( stderr, " The QR algorithm failed to converge.\n" );
+ }
+ /* for (l2=k; l2>=0; l2--) { */
+ for (l2=k; l2>=1; l2--) {
+ l = l2;
+#ifdef DEBUGPRAX
+ printf(" l e(l)< eps %d %.7g %.7g ",l,e[l], eps);
+#endif
+ /* if (fabs(e[l]) <= eps) */
+ if (fabs(e[l]) <= eps)
+ goto TestFconvergence;
+ /* if (fabs(q[l-1]) <= eps)*/ /* missing if ( 1 < l ){ *//* printf(" q(l-1)< eps %d %.7g %.7g ",l-1,q[l-2], eps); */
+ if (fabs(q[l-1]) <= eps)
+ break; /* goto Cancellation; */
+ }
+ Cancellation:
+#ifdef DEBUGPRAX
+ printf(" Cancellation:\n");
+#endif
+ c = 0.0; s = 1.0;
+ for (i=l; i<=k; i++) {
+ f = s * e[i]; e[i] *= c;
+ /* f = s * e[i]; e[i] *= c; */
+ if (fabs(f) <= eps)
+ goto TestFconvergence;
+ /* g = q[i]; */
+ g = q[i];
+ if (fabs(f) < fabs(g)) {
+ double fg = f/g;
+ h = fabs(g)*sqrt(1.0+fg*fg);
+ }
+ else {
+ double gf = g/f;
+ h = (f!=0.0 ? fabs(f)*sqrt(1.0+gf*gf) : 0.0);
+ }
+ /* COMMENT: THE ABOVE REPLACES Q(I):=H:=LONGSQRT(G*G+F*F) */
+ /* WHICH MAY GIVE INCORRECT RESULTS IF THE */
+ /* SQUARES UNDERFLOW OR IF F = G = 0; */
+
+ /* q[i] = h; */
+ q[i] = h;
+ if (h == 0.0) { h = 1.0; g = 1.0; }
+ c = g/h; s = -f/h;
+ }
+TestFconvergence:
+ #ifdef DEBUGPRAX
+ printf(" TestFconvergence: l=%d k=%d\n",l,k);
+#endif
+ /* z = q[k]; */
+ z = q[k];
+ if (l == k)
+ goto Convergence;
+ /* shift from bottom 2x2 minor */
+ /* x = q[l]; y = q[k-l]; g = e[k-1]; h = e[k]; */ /* Error */
+ x = q[l]; y = q[k-1]; g = e[k-1]; h = e[k];
+ f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);
+ g = sqrt(f*f+1.0);
+ if (f <= 0.0)
+ f = ((x-z)*(x+z) + h*(y/(f-g)-h))/x;
+ else
+ f = ((x-z)*(x+z) + h*(y/(f+g)-h))/x;
+ /* next qr transformation */
+ s = c = 1.0;
+ for (i=l+1; i<=k; i++) {
+#ifdef DEBUGPRAXQR
+ printf(" Before Mid TestFconvergence: l+1=%d i=%d k=%d h=%.6e e(i)=%14.8f e(i-1)=%14.8f\n",l+1,i,k, h, e[i],e[i-1]);
+#endif
+ /* g = e[i]; y = q[i]; h = s*g; g *= c; */
+ g = e[i]; y = q[i]; h = s*g; g *= c;
+ if (fabs(f) < fabs(h)) {
+ double fh = f/h;
+ z = fabs(h) * sqrt(1.0 + fh*fh);
+ }
+ else {
+ double hf = h/f;
+ z = (f!=0.0 ? fabs(f)*sqrt(1.0+hf*hf) : 0.0);
+ }
+ /* e[i-1] = z; */
+ e[i-1] = z;
+#ifdef DEBUGPRAXQR
+ printf(" Mid TestFconvergence: l+1=%d i=%d k=%d h=%.6e e(i)=%14.8f e(i-1)=%14.8f\n",l+1,i,k, h, e[i],e[i-1]);
+#endif
+ if (z == 0.0)
+ f = z = 1.0;
+ c = f/z; s = h/z;
+ f = x*c + g*s; g = - x*s + g*c; h = y*s;
+ y *= c;
+ /* for (j=0; j<n; j++) { */
+ /* x = ab[j][i-1]; z = ab[j][i]; */
+ /* ab[j][i-1] = x*c + z*s; */
+ /* ab[j][i] = - x*s + z*c; */
+ /* } */
+ for (j=1; j<=n; j++) {
+ x = ab[j][i-1]; z = ab[j][i];
+ ab[j][i-1] = x*c + z*s;
+ ab[j][i] = - x*s + z*c;
+ }
+ if (fabs(f) < fabs(h)) {
+ double fh = f/h;
+ z = fabs(h) * sqrt(1.0 + fh*fh);
+ }
+ else {
+ double hf = h/f;
+ z = (f!=0.0 ? fabs(f)*sqrt(1.0+hf*hf) : 0.0);
+ }
+#ifdef DEBUGPRAXQR
+ printf(" qr transformation z f h=%.7g %.7g %.7g i=%d k=%d\n",z,f,h, i, k);
+#endif
+ q[i-1] = z;
+ if (z == 0.0)
+ z = f = 1.0;
+ c = f/z; s = h/z;
+ f = c*g + s*y; /* f can be very small */
+ x = - s*g + c*y;
+ }
+ /* e[l] = 0.0; e[k] = f; q[k] = x; */
+ e[l] = 0.0; e[k] = f; q[k] = x;
+#ifdef DEBUGPRAXQR
+ printf(" aftermid loop l=%d k=%d e(l)=%7g e(k)=%.7g q(k)=%.7g x=%.7g\n",l,k,e[l],e[k],q[k],x);
+#endif
+ goto TestFsplitting;
+ Convergence:
+#ifdef DEBUGPRAX
+ printf(" Convergence:\n");
+#endif
+ if (z < 0.0) {
+ /* q[k] = - z; */
+ /* for (j=0; j<n; j++) ab[j][k] = - ab[j][k]; */
+ q[k] = - z;
+ for (j=1; j<=n; j++) ab[j][k] = - ab[j][k];
+ }/* END Z */
+ }/* END K */
+} /* END MINFIT */
+
+
+double praxis(double tol, double macheps, double h0, int _n, int _prin, double *_x, double (*_fun)(double *_x))
+/* double praxis(double tol, double macheps, double h0, int _n, int _prin, double *_x, double (*_fun)(double *_x, int _n)) */
+/* double praxis(double (*_fun)(), double _x[], int _n) */
+/* double (*_fun)(); */
+/* double _x[N]; */
+/* double (*_fun)(); */
+/* double _x[N]; */
+{
+ /* init global extern variables and parameters */
+ /* double *d, *y, *z, */
+ /* *q0, *q1, **v; */
+ /* double *tflin; /\* used in flin: return (*fun)(tflin, n); *\/ */
+ /* double *e; /\* used in minfit, don't konw how to free memory and thus made global *\/ */
- double pythag(double a, double b);
- int flag,i,its,j,jj,k,l,nm;
- double anorm,c,f,g,h,s,scale,x,y,z,*rv1;
+ int seed; /* added */
+ int biter=0;
+ double r;
+ double randbrent( int (*));
+ double s, sf;
- rv1=vector(1,n);
- /* Householder's reduction to bidiagonal form; */
- g=scale=anorm=0.0;
- for (i=1;i<=n;i++) {
- l=i+1;
- rv1[i]=scale*g;
- g=s=scale=0.0;
- if (i <= m) {
- for (k=i;k<=m;k++) scale += fabs(a[k][i]);
- if (scale) {
- for (k=i;k<=m;k++) {
- a[k][i] /= scale;
- s += a[k][i]*a[k][i];
- }
- f=a[i][i];
- g = -SIGN(sqrt(s),f);
- h=f*g-s;
- a[i][i]=f-g;
- for (j=l;j<=n;j++) {
- for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
- f=s/h;
- for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
- }
- for (k=i;k<=m;k++) a[k][i] *= scale;
- }
- }
- w[i]=scale *g; /* p 411*/
- g=s=scale=0.0;
- if (i <= m && i != n) {
- for (k=l;k<=n;k++) scale += fabs(a[i][k]);
- if (scale) {
- for (k=l;k<=n;k++) {
- a[i][k] /= scale;
- s += a[i][k]*a[i][k];
- }
- f=a[i][l];
- g = -SIGN(sqrt(s),f);
- h=f*g-s;
- a[i][l]=f-g;
- for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
- for (j=l;j<=m;j++) {
- for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
- for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
- }
- for (k=l;k<=n;k++) a[i][k] *= scale;
- }
- }
- /* y : = abs(q[i]) +abs(e[i]); if y > x then x : = y */
- /* anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i]))); */
- anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
- }
- /* Comment accumulation of right-hand transformations */
- for (i=n;i>=1;i--) {
- if (i < n) {
- if (g) {
- for (j=l;j<=n;j++) a[j][i]=(a[i][j]/a[i][l])/g; /* Double division to avoid possible underflow. */
- /* for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g; */
- for (j=l;j<=n;j++) {
- /* for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j]; */
- for (s=0.0,k=l;k<=n;k++) s += a[i][k]*a[k][j];
- /* for (k=l;k<=n;k++) v[k][j] += s*v[k][i]; */
- for (k=l;k<=n;k++) a[k][j] += s*a[k][i];
- }
- }
- /* for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0; */
- for (j=l;j<=n;j++) a[i][j]=a[j][i]=0.0;
- }
- /* v[i][i]=1.0; */
- a[i][i]=1.0;
- g=rv1[i];
- l=i;
- }
- /* Comment accumulation of left-hand transformations; */
- /* for (i=IMIN(m,n);i>=1;i--) { */
- for (i=FMIN(m,n);i>=1;i--) {
- l=i+1;
- g=w[i];
- for (j=l;j<=n;j++) a[i][j]=0.0;
- if (g) {
- g=1.0/g;
- for (j=l;j<=n;j++) {
- for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
- f=(s/a[i][i])*g;
- for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
- }
- for (j=i;j<=m;j++) a[j][i] *= g;
- } else for (j=i;j<=m;j++) a[j][i]=0.0;
- ++a[i][i];
- }
-/* comment diagonalization of the bidiagonal form; */
- for (k=n;k>=1;k--) {
- for (its=1;its<=30;its++) { /* Loop over singular values, and over allowed iterations. */
- flag=1;
- for (l=k;l>=1;l--) { /* Test for splitting. */
- nm=l-1; /* Note that rv1[1] is always zero. */
- if ((double)(fabs(rv1[l])+anorm) == anorm) {
- flag=0;
- break;
- }
- if ((double)(fabs(w[nm])+anorm) == anorm) break;
- }
- if (flag) {
- c=0.0; /* Cancellation of rv1[l], if l > 1. */
- s=1.0;
- for (i=l;i<=k;i++) {
- f=s*rv1[i];
- rv1[i]=c*rv1[i];
- if ((double)(fabs(f)+anorm) == anorm) break;
- g=w[i];
- h=pythag(f,g);
- w[i]=h;
- h=1.0/h;
- c=g*h;
- s = -f*h;
- for (j=1;j<=m;j++) {
- y=a[j][nm];
- z=a[j][i];
- a[j][nm]=y*c+z*s;
- a[j][i]=z*c-y*s;
- }
- }
- }
- z=w[k];
- if (l == k) { /* Convergence */
- if (z < 0.0) { /* Singular value is made nonnegative. */
- w[k] = -z;
- /* for (j=1;j<=n;j++) v[j][k] = -v[j][k]; */
- for (j=1;j<=n;j++) a[j][k] = -a[j][k];
- }
- break;
- }
- 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];
- g=rv1[nm];
- h=rv1[k];
- f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
- g=pythag(f,1.0);
- /* g=dpythag(f,1.0); */
- f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
- c=s=1.0; /* next QR transformation */
- for (j=l;j<=nm;j++) {
- i=j+1;
- g=rv1[i];
- y=w[i];
- h=s*g;
- g=c*g;
- /* z=dpythag(f,h); */
- z=pythag(f,h);
- rv1[j]=z;
- c=f/z;
- s=h/z;
- f=x*c+g*s;
- g = g*c-x*s;
- h=y*s;
- y *= c;
-/* if withv then for j:= 1 step 1 until n do */
- for (jj=1;jj<=n;jj++) {
- /* x=v[jj][j]; */
- /* z=v[jj][i]; */
- x=a[jj][j];
- z=a[jj][i];
- /* v[jj][j]=x*c+z*s; */
- /* v[jj][i]=z*c-x*s; */
- a[jj][j]=x*c+z*s;
- a[jj][i]=z*c-x*s;
- }
- /* z=dpythag(f,h); */
- z=pythag(f,h);
- w[j]=z;
- if (z) {
- z=1.0/z;
- c=f*z;
- s=h*z;
- }
- f=c*g+s*y;
- x=c*y-s*g;
- /* if withu then for j:=1 step 1 until m do */
- /* for (jj=1;jj<=m;jj++) { */
- /* y=a[jj][j]; */
- /* z=a[jj][i]; */
- /* a[jj][j]=y*c+z*s; */
- /* a[jj][i]=z*c-y*s; */
- /* } */
- }
- rv1[l]=0.0;
- rv1[k]=f;
- w[k]=x;
- }
- }
- free_vector(rv1,1,n);
-}
+ h = h0; /* step; */
+ t = tol;
+ scbd = 1.0;
+ illc = 0;
+ ktm = 1;
+
+ macheps = DBL_EPSILON;
+ /* prin=4; */
+#ifdef DEBUGPRAX
+ printf("Praxis macheps=%14g h=%14g step=%14g tol=%14g\n",macheps,h, h0,tol);
+#endif
+ n = _n;
+ x = _x;
+ prin = _prin;
+ fun = _fun;
+ d=vector(1, n);
+ y=vector(1, n);
+ z=vector(1, n);
+ q0=vector(1, n);
+ q1=vector(1, n);
+ e=vector(1, n);
+ tflin=vector(1, n);
+ v=matrix(1, n, 1, n);
+ for(i=1;i<=n;i++){d[i]=y[i]=z[i]=q0[0]=e[i]=tflin[i]=0.;}
+ small_windows = (macheps) * (macheps); vsmall = small_windows*small_windows;
+ large = 1.0/small_windows; vlarge = 1.0/vsmall;
+ m2 = sqrt(macheps); m4 = sqrt(m2);
+ seed = 123456789; /* added */
+ ldfac = (illc ? 0.1 : 0.01);
+ for(i=1;i<=n;i++) z[i]=0.; /* Was missing in Gegenfurtner as well as Brent's algol or fortran */
+ nl = kt = 0; nf = 1;
+#ifdef NR_SHIFT
+ fx = (*fun)((x-1), n);
+#else
+ fx = (*fun)(x);
+#endif
+ qf1 = fx;
+ t2 = small_windows + fabs(t); t = t2; dmin = small_windows;
+#ifdef DEBUGPRAX
+ printf("praxis2 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);
+#endif
+ if (h < 100.0*t) h = 100.0*t;
+#ifdef DEBUGPRAX
+ printf("praxis3 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);
+#endif
+ ldt = h;
+ /* for (i=0; i<n; i++) for (j=0; j<n; j++) */
+ for (i=1; i<=n; i++) for (j=1; j<=n; j++)
+ v[i][j] = (i == j ? 1.0 : 0.0);
+ d[1] = 0.0; qd0 = 0.0;
+ /* for (i=0; i<n; i++) q1[i] = x[i]; */
+ for (i=1; i<=n; i++) q1[i] = x[i];
+ if (prin > 1) {
+ printf("\n------------- enter function praxis -----------\n");
+ printf("... current parameter settings ...\n");
+ printf("... scaling ... %20.10e\n", scbd);
+ printf("... tol ... %20.10e\n", t);
+ printf("... maxstep ... %20.10e\n", h);
+ printf("... illc ... %20u\n", illc);
+ printf("... ktm ... %20u\n", ktm);
+ printf("... maxfun ... %20u\n", maxfun);
+ }
+ if (prin) print2();
+
+mloop:
+ biter++; /* Added to count the loops */
+ /* sf = d[0]; */
+ /* s = d[0] = 0.0; */
+ printf("\n Big iteration %d \n",biter);
+ fprintf(ficlog,"\n Big iteration %d \n",biter);
+ sf = d[1];
+ s = d[1] = 0.0;
+
+ /* minimize along first direction V(*,1) */
+#ifdef DEBUGPRAX
+ printf(" Minimize along the first direction V(*,1). illc=%d\n",illc);
+ /* fprintf(ficlog," Minimize along the first direction V(*,1).\n"); */
+#endif
+#ifdef DEBUGPRAX2
+ printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);
+#endif
+ /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */
+ minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global */
+#ifdef DEBUGPRAX
+ printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx);
+#endif
+ if (s <= 0.0)
+ /* for (i=0; i < n; i++) */
+ for (i=1; i <= n; i++)
+ v[i][1] = -v[i][1];
+ /* if ((sf <= (0.9 * d[0])) || ((0.9 * sf) >= d[0])) */
+ if ((sf <= (0.9 * d[1])) || ((0.9 * sf) >= d[1]))
+ /* for (i=1; i<n; i++) */
+ for (i=2; i<=n; i++)
+ d[i] = 0.0;
+ /* for (k=1; k<n; k++) { */
+ for (k=2; k<=n; k++) {
+ /*
+ The inner loop starts here.
+ */
+#ifdef DEBUGPRAX
+ printf(" The inner loop here from k=%d to n=%d.\n",k,n);
+ /* fprintf(ficlog," The inner loop here from k=%d to n=%d.\n",k,n); */
+#endif
+ /* for (i=0; i<n; i++) */
+ for (i=1; i<=n; i++)
+ y[i] = x[i];
+ sf = fx;
+#ifdef DEBUGPRAX
+ printf(" illc=%d and kt=%d and ktm=%d\n", illc, kt, ktm);
+#endif
+ illc = illc || (kt > 0);
+next:
+ kl = k;
+ df = 0.0;
+ if (illc) { /* random step to get off resolution valley */
+#ifdef DEBUGPRAX
+ printf(" A random step follows, to avoid resolution valleys.\n");
+ matprint(" before rand, vectors:",v,n,n);
+#endif
+ for (i=1; i<=n; i++) {
+#ifdef NOBRENTRAND
+ r = drandom();
+#else
+ seed=i;
+ /* seed=i+1; */
+#ifdef DEBUGRAND
+ printf(" Random seed=%d, brent i=%d",seed,i); /* YYYY i=5 j=1 vji= -0.0001170073 */
+#endif
+ r = randbrent ( &seed );
+#endif
+#ifdef DEBUGRAND
+ printf(" Random r=%.7g \n",r);
+#endif
+ z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (r - 0.5);
+ /* z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (drandom() - 0.5); */
+
+ s = z[i];
+ for (j=1; j <= n; j++)
+ x[j] += s * v[j][i];
+ }
+#ifdef DEBUGRAND
+ matprint(" after rand, vectors:",v,n,n);
+#endif
+#ifdef NR_SHIFT
+ fx = (*fun)((x-1), n);
+#else
+ fx = (*fun)(x, n);
+#endif
+ /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/
+ nf++;
+ }
+ /* minimize along non-conjugate directions */
+#ifdef DEBUGPRAX
+ printf(" Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n);
+ /* fprintf(ficlog," Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n); */
+#endif
+ /* for (k2=k; k2<n; k2++) { /\* Be careful here k2 <=n ? *\/ */
+ for (k2=k; k2<=n; k2++) { /* Be careful here k2 <=n ? */
+ sl = fx;
+ s = 0.0;
+#ifdef DEBUGPRAX
+ printf(" Minimize along the 'NON-CONJUGATE' true direction k2=%14d fx=%14.7f\n",k2, fx);
+ matprint(" before min vectors:",v,n,n);
+#endif
+ /* min(k2, 2, &d[k2], &s, fx, 0); */
+ /* jsearch=k2-1; */
+ /* min(jsearch, 2, &d[jsearch], &s, fx, 0); */
+ minny(k2, 2, &d[k2], &s, fx, 0);
+#ifdef DEBUGPRAX
+ printf(" . D(%d)=%14.7f d[k2]=%14.7f z[k2]=%14.7f illc=%14d fx=%14.7f\n",k2,d[k2],d[k2],z[k2],illc,fx);
+#endif
+ if (illc) {
+ /* double szk = s + z[k2]; */
+ /* s = d[k2] * szk*szk; */
+ double szk = s + z[k2];
+ s = d[k2] * szk*szk;
+ }
+ else
+ s = sl - fx;
+ /* if (df < s) { */
+ if (df <= s) {
+ df = s;
+ kl = k2;
+#ifdef DEBUGPRAX
+ printf(" df=%.7g and choose kl=%d \n",df,kl); /* UUUU */
+#endif
+ }
+ } /* end loop k2 */
+ /*
+ If there was not much improvement on the first try, set
+ ILLC = true and start the inner loop again.
+ */
+#ifdef DEBUGPRAX
+ printf(" If there was not much improvement on the first try, set ILLC = true and start the inner loop again. illc=%d\n",illc);
+ /* fprintf(ficlog," If there was not much improvement on the first try, set ILLC = true and start the inner loop again.\n"); */
+#endif
+ if (!illc && (df < fabs(100.0 * (macheps) * fx))) {
+#ifdef DEBUGPRAX
+ printf("\n NO SUCCESS because DF is small, starts inner loop with same K(=%d), fabs( 100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e > df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc);
+#endif
+ illc = 1;
+ goto next;
+ }
+#ifdef DEBUGPRAX
+ printf("\n SUCCESS, BREAKS inner loop K(=%d) because DF is big, fabs( 100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e <= df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc);
+#endif
+
+ /* if ((k == 1) && (prin > 1)){ /\* be careful k=2 *\/ */
+ if ((k == 2) && (prin > 1)){ /* be careful k=2 */
+#ifdef DEBUGPRAX
+ printf(" NEW D The second difference array d:\n" );
+ /* fprintf(ficlog, " NEW D The second difference array d:\n" ); */
+#endif
+ vecprint(" NEW D The second difference array d:",d,n);
+ }
+ /* minimize along conjugate directions */
+ /*
+ Minimize along the "conjugate" directions V(*,1),...,V(*,K-1).
+ */
+#ifdef DEBUGPRAX
+ printf("Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1);
+ /* fprintf(ficlog,"Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1); */
+#endif
+ /* for (k2=0; k2<=k-1; k2++) { */
+ for (k2=1; k2<=k-1; k2++) {
+ s = 0.0;
+ /* min(k2-1, 2, &d[k2-1], &s, fx, 0); */
+ minny(k2, 2, &d[k2], &s, fx, 0);
+ }
+ f1 = fx;
+ fx = sf;
+ lds = 0.0;
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ sl = x[i];
+ x[i] = y[i];
+ y[i] = sl - y[i];
+ sl = y[i];
+ lds = lds + sl*sl;
+ }
+ lds = sqrt(lds);
+#ifdef DEBUGPRAX
+ printf("Minimization done 'conjugate', shifted all points, computed lds=%.8f\n",lds);
+#endif
+ /*
+ Discard direction V(*,kl).
+
+ If no random step was taken, V(*,KL) is the "non-conjugate"
+ direction along which the greatest improvement was made.
+ */
+ if (lds > small_windows) {
+#ifdef DEBUGPRAX
+ printf("lds big enough to throw direction V(*,kl=%d). If no random step was taken, V(*,KL) is the 'non-conjugate' direction along which the greatest improvement was made.\n",kl);
+ matprint(" before shift new conjugate vectors:",v,n,n);
+#endif
+ for (i=kl-1; i>=k; i--) {
+ /* for (j=0; j < n; j++) */
+ for (j=1; j <= n; j++)
+ /* v[j][i+1] = v[j][i]; */ /* This is v[j][i+1]=v[j][i] i=kl-1 to k */
+ v[j][i+1] = v[j][i]; /* This is v[j][i+1]=v[j][i] i=kl-1 to k */
+ /* v[j][i+1] = v[j][i]; */
+ /* d[i+1] = d[i];*/ /* last is d[k+1]= d[k] */
+ d[i+1] = d[i]; /* last is d[k]= d[k-1] */
+ }
+#ifdef DEBUGPRAX
+ matprint(" after shift new conjugate vectors:",v,n,n);
+#endif /* d[k] = 0.0; */
+ d[k] = 0.0;
+ for (i=1; i <= n; i++)
+ v[i][k] = y[i] / lds;
+ /* v[i][k] = y[i] / lds; */
+#ifdef DEBUGPRAX
+ printf("Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector: (new x) - (old x). d2=%14.7g lds=%.10f\n",k,d[k],lds);
+ /* fprintf(ficlog,"Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector: (new x) - (old x).\n",k); */
+ matprint(" before min new conjugate vectors:",v,n,n);
+#endif
+ /* min(k-1, 4, &d[k-1], &lds, f1, 1); */
+ minny(k, 4, &d[k], &lds, f1, 1);
+#ifdef DEBUGPRAX
+ printf(" after min d(k)=%d %.7g lds=%14f\n",k,d[k],lds);
+ matprint(" after min vectors:",v,n,n);
+#endif
+ if (lds <= 0.0) {
+ lds = -lds;
+#ifdef DEBUGPRAX
+ printf(" lds changed sign lds=%.14f k=%d\n",lds,k);
+#endif
+ /* for (i=0; i<n; i++) */
+ /* v[i][k] = -v[i][k]; */
+ for (i=1; i<=n; i++)
+ v[i][k] = -v[i][k];
+ }
+ }
+ ldt = ldfac * ldt;
+ if (ldt < lds)
+ ldt = lds;
+ if (prin > 0){
+#ifdef DEBUGPRAX
+ printf(" k=%d",k);
+ /* fprintf(ficlog," k=%d",k); */
+#endif
+ print2();/* n, x, prin, fx, nf, nl ); */
+ }
+ t2 = 0.0;
+ /* for (i=0; i<n; i++) */
+ for (i=1; i<=n; i++)
+ t2 += x[i]*x[i];
+ t2 = m2 * sqrt(t2) + t;
+ /*
+ See whether the length of the step taken since starting the
+ inner loop exceeds half the tolerance.
+ */
+#ifdef DEBUGPRAX
+ printf("See if step length exceeds half the tolerance.\n"); /* ZZZZZ */
+ /* fprintf(ficlog,"See if step length exceeds half the tolerance.\n"); */
+#endif
+ if (ldt > (0.5 * t2))
+ kt = 0;
+ else
+ kt++;
+#ifdef DEBUGPRAX
+ printf("if kt=%d >? ktm=%d gotoL2 loop\n",kt,ktm);
+#endif
+ if (kt > ktm){
+ if ( 0 < prin ){
+ /* printf("\nr8vec_print\n X:\n"); */
+ /* fprintf(ficlog,"\nr8vec_print\n X:\n"); */
+ vecprint ("END X:", x, n );
+ }
+ goto fret;
+ }
+#ifdef DEBUGPRAX
+ matprint(" end of L2 loop vectors:",v,n,n);
+#endif
+
+ }
+ /* printf("The inner loop ends here.\n"); */
+ /* fprintf(ficlog,"The inner loop ends here.\n"); */
+ /*
+ The inner loop ends here.
+
+ Try quadratic extrapolation in case we are in a curved valley.
+ */
+#ifdef DEBUGPRAX
+ printf("Try QUAD ratic extrapolation in case we are in a curved valley.\n");
+#endif
+ /* try quadratic extrapolation in case */
+ /* we are stuck in a curved valley */
+ quad();
+ dn = 0.0;
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ d[i] = 1.0 / sqrt(d[i]);
+ if (dn < d[i])
+ dn = d[i];
+ }
+ if (prin > 2)
+ matprint(" NEW DIRECTIONS vectors:",v,n,n);
+ /* for (j=0; j<n; j++) { */
+ for (j=1; j<=n; j++) {
+ s = d[j] / dn;
+ /* for (i=0; i < n; i++) */
+ for (i=1; i <= n; i++)
+ v[i][j] *= s;
+ }
+
+ if (scbd > 1.0) { /* scale axis to reduce condition number */
+#ifdef DEBUGPRAX
+ printf("Scale the axes to try to reduce the condition number.\n");
+#endif
+ /* fprintf(ficlog,"Scale the axes to try to reduce the condition number.\n"); */
+ s = vlarge;
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ sl = 0.0;
+ /* for (j=0; j < n; j++) */
+ for (j=1; j <= n; j++)
+ sl += v[i][j]*v[i][j];
+ z[i] = sqrt(sl);
+ if (z[i] < m4)
+ z[i] = m4;
+ if (s > z[i])
+ s = z[i];
+ }
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ sl = s / z[i];
+ z[i] = 1.0 / sl;
+ if (z[i] > scbd) {
+ sl = 1.0 / scbd;
+ z[i] = scbd;
+ }
+ }
+ }
+ for (i=1; i<=n; i++)
+ /* for (j=0; j<=i-1; j++) { */
+ /* for (j=1; j<=i; j++) { */
+ for (j=1; j<=i-1; j++) {
+ s = v[i][j];
+ v[i][j] = v[j][i];
+ v[j][i] = s;
+ }
+#ifdef DEBUGPRAX
+ printf(" Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n");
+#endif
+ /*
+ 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.
+ */
+ #ifdef DEBUGPRAX
+ printf(" MINFIT 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");
+#endif
+
+ minfit(n, macheps, vsmall, v, d);
+ /* for(i=0; i<n;i++)printf(" %14.7g",d[i]); */
+ /* v is overwritten with R. */
+ /*
+ Unscale the axes.
+ */
+ if (scbd > 1.0) {
+#ifdef DEBUGPRAX
+ printf(" Unscale the axes.\n");
+#endif
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ s = z[i];
+ /* for (j=0; j<n; j++) */
+ for (j=1; j<=n; j++)
+ v[i][j] *= s;
+ }
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ s = 0.0;
+ /* for (j=0; j<n; j++) */
+ for (j=1; j<=n; j++)
+ s += v[j][i]*v[j][i];
+ s = sqrt(s);
+ d[i] *= s;
+ s = 1.0 / s;
+ /* for (j=0; j<n; j++) */
+ for (j=1; j<=n; j++)
+ v[j][i] *= s;
+ }
+ }
+ /* for (i=0; i<n; i++) { */
+ double dni; /* added for compatibility with buckhardt but not brent */
+ for (i=1; i<=n; i++) {
+ dni=dn*d[i]; /* added for compatibility with buckhardt but not brent */
+ if ((dn * d[i]) > large)
+ d[i] = vsmall;
+ else if ((dn * d[i]) < small_windows)
+ d[i] = vlarge;
+ else
+ d[i] = 1.0 / dni / dni; /* added for compatibility with buckhardt but not brent */
+ /* d[i] = pow(dn * d[i],-2.0); */
+ }
+#ifdef DEBUGPRAX
+ vecprint ("\n Before sort Eigenvalues of a:",d,n );
+#endif
+
+ sort(); /* the new eigenvalues and eigenvectors */
+#ifdef DEBUGPRAX
+ vecprint( " After sort the eigenvalues ....\n", d, n);
+ matprint( " After sort the eigenvectors....\n", v, n,n);
+#endif
+#ifdef DEBUGPRAX
+ printf(" Determine the smallest eigenvalue.\n");
+#endif
+ /* dmin = d[n-1]; */
+ dmin = d[n];
+ if (dmin < small_windows)
+ dmin = small_windows;
+ /*
+ The ratio of the smallest to largest eigenvalue determines whether
+ the system is ill conditioned.
+ */
+
+ /* illc = (m2 * d[0]) > dmin; */
+ illc = (m2 * d[1]) > dmin;
+#ifdef DEBUGPRAX
+ printf(" The ratio of the smallest to largest eigenvalue determines whether\n the system is ill conditioned=%d . dmin=%.10lf < m2=%.10lf * d[1]=%.10lf \n",illc, dmin,m2, d[1]);
+#endif
+
+ if ((prin > 2) && (scbd > 1.0))
+ vecprint("\n The scale factors:",z,n);
+ if (prin > 2)
+ vecprint(" Principal values (EIGEN VALUES OF A) of the quadratic form:",d,n);
+ if (prin > 2)
+ matprint(" The principal axes (EIGEN VECTORS OF A:",v,n, n);
+
+ if ((maxfun > 0) && (nf > maxfun)) {
+ if (prin)
+ printf("\n... maximum number of function calls reached ...\n");
+ goto fret;
+ }
+#ifdef DEBUGPRAX
+ printf("Goto main loop\n");
+#endif
+ goto mloop; /* back to main loop */
+
+fret:
+ if (prin > 0) {
+ vecprint("\n X:", x, n);
+ /* printf("\n... ChiSq reduced to %20.10e ...\n", fx); */
+ /* printf("... after %20u function calls.\n", nf); */
+ }
+ free_vector(d, 1, n);
+ free_vector(y, 1, n);
+ free_vector(z, 1, n);
+ free_vector(q0, 1, n);
+ free_vector(q1, 1, n);
+ free_matrix(v, 1, n, 1, n);
+ /* double *d, *y, *z, */
+ /* *q0, *q1, **v; */
+ free_vector(tflin, 1, n);
+ /* double *tflin; /\* used in flin: return (*fun)(tflin, n); *\/ */
+ free_vector(e, 1, n);
+ /* double *e; /\* used in minfit, don't konw how to free memory and thus made global *\/ */
+
+ return(fx);
+}
-/* end praxis */
+/* end praxis gegen */
/*************** powell ************************/
/*
}
}
for (i=1;i<=n;i++) { /* For each direction i, maximisation after loading directions */
- for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */
- fptt=(*fret);
+ for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales. xi is not changed but one dim xit */
+
+ fptt=(*fret); /* Computes likelihood for parameters xit */
#ifdef DEBUG
printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret);
printf("%d",i);fflush(stdout); /* print direction (parameter) i */
fprintf(ficlog,"%d",i);fflush(ficlog);
#ifdef LINMINORIGINAL
- linmin(p,xit,n,fret,func); /* New point i minimizing in direction xit i has coordinates p[j].*/
+ linmin(p,xit,n,fret,func); /* New point i minimizing in direction xit, 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
} /* end loop on each direction i */
/* 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]);
/* the scales of the directions and the directions, because the are reset to canonical directions */
/* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */
/* under the tolerance value. If the tolerance is very small 1.e-9, it could last long. */
+#ifdef DEBUG
+ int k[2],l;
+ k[0]=1;
+ k[1]=-1;
+ printf("Max: %.12e",(*func)(p));
+ fprintf(ficlog,"Max: %.12e",(*func)(p));
+ for (j=1;j<=n;j++) {
+ printf(" %.12e",p[j]);
+ fprintf(ficlog," %.12e",p[j]);
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
+ for(l=0;l<=1;l++) {
+ for (j=1;j<=n;j++) {
+ ptt[j]=p[j]+(p[j]-pt[j])*k[l];
+ printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
+ fprintf(ficlog,"l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
+ }
+ printf("func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
+ fprintf(ficlog,"func(ptt)=%.12e, deriv=%.12e\n",(*func)(ptt),(ptt[j]-p[j])/((*func)(ptt)-(*func)(p)));
+ }
+#endif
free_vector(xit,1,n);
free_vector(xits,1,n);
return;
} /* enough precision */
if (*iter == ITMAX*n) nrerror("powell exceeding maximum iterations.");
- for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0)=2Pn-P0 */
+ for (j=1;j<=n;j++) { /* Computes the extrapolated point and value f3, P_0 + 2 (P_n-P_0)=2Pn-P0 and xit is direction Pn-P0 */
ptt[j]=2.0*p[j]-pt[j];
- xit[j]=p[j]-pt[j]; /* Coordinate j of last direction xi_n=P_n-_0 */
+ xit[j]=p[j]-pt[j]; /* Coordinate j of last direction xi_n=P_n-P_0 */
+#ifdef DEBUG
printf("\n %d xit=%12.7g p=%12.7g pt=%12.7g ",j,xit[j],p[j],pt[j]);
- pt[j]=p[j];
+#endif
+ pt[j]=p[j]; /* New P0 is Pn */
}
+#ifdef DEBUG
printf("\n");
+#endif
fptt=(*func)(ptt); /* f_3 */
-#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */
+#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in directions until some iterations are done */
if (*iter <=4) {
#else
#endif
t= t- del*SQR(fp-fptt);
#endif
directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */
- printf(" t=%g, directest=%g\n",t, directest);
-#ifdef POWELLNOTTRUECONJUGATE /* Searching for IBIG and testing for replacement */
+#ifdef DEBUG
+ printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
+ fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
+ printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
+ (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
+ fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
+ (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
+ printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
+ fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
+#endif
#ifdef POWELLORIGINAL
if (t < 0.0) { /* Then we use it for new direction */
-#else /* Not POWELLOriginal but Brouard's */
+#else
if (directest*t < 0.0) { /* Contradiction between both tests */
printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
}
- if (directest < 0.0) { /* Then we use (P0, Pn) for new direction Xi_n or Xi_iBig */
-#endif /* end POWELLOriginal */
-#endif /* POWELLNOTTRUECONJUGATE else means systematic replacement by new direction P_0P_n */
+ if (directest < 0.0) { /* Then we use it for new direction */
+#endif
+#ifdef DEBUGLINMIN
+ printf("Before linmin in direction P%d-P0\n",n);
+ for (j=1;j<=n;j++) {
+ printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
+ printf("\n");
+ fprintf(ficlog,"\n");
+ }
+ }
+#endif
#ifdef LINMINORIGINAL
- /* xit[j]=p[j]-pt[j] */
- printf("\n Computes min on P_0, P_n direction iter=%d Bigter=%d\n",*iter,Bigter);
- linmin(p,xit,n,fret,func); /* computes minimum on P_0,P_n direction: changes p and rescales xit.*/
-#else /* NOT LINMINORIGINAL but with searching for flat directions */
- printf("\n Flat Computes min on P_0, P_n direction iter=%d Bigter=%d\n",*iter,Bigter);
+ linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
+#else
linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
flatdir[i]=flat; /* Function is vanishing in that direction i */
#endif
-#ifdef POWELLNOTTRUECONJUGATE
-#else
-#ifdef POWELLORIGINCONJUGATE
+#ifdef DEBUGLINMIN
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 (i=1;i<=n-1;i++) {
- for (j=1;j<=n;j++) {
- xi[j][i]=xi[j][i+1]; /* Standard method of conjugate directions, not Powell who changes the nth direction by p0 pn . */
+ printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
+ if(j % ncovmodel == 0){
+ printf("\n");
+ fprintf(ficlog,"\n");
}
}
+#endif
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 */
}
-#endif /* POWELLORIGINCONJUGATE*/
-#endif /*POWELLNOTTRUECONJUGATE*/
- printf(" Standard method of conjugate directions\n");
- printf("\n#A Before prax Bigter=%d model= 1 + age ", Bigter);
- for(j=1;j<=n;j++){
- printf("%d \n",j);
- for(i=1;i<=n;i++){
- printf(" %f",xi[j][i]);
- }
- }
- printf("\n");
-
-#ifdef NOTMINFIT
-#else
- if(*iter >n){
- /* if(Bigter >n){ */
- printf("\n#Before prax Bigter=%d model= 1 + age ", Bigter);
- printf("\n");
- for(j=1;j<=n;j++){
- printf("%d \n",j);
- for(i=1;i<=n;i++){
- printf(" %f",xi[j][i]);
- }
- }
- printf("\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(" Bigter=%d Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n",Bigter);
- 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 *d; /* eigenvalues of principal directions */
- d=vector(1,n);
-
-
- svdminfit (xi, n, n, d ); /* In Brent's notation find d such that V=Q Diagonal(d) R, and Lambda=d^(-1/2) */
-
- printf("\n#After prax model= 1 + age ");
- fprintf(ficlog,"\n#model= 1 + age ");
-
- if(nagesqr==1){
- printf(" + age*age ");
- fprintf(ficlog," + age*age ");
- }
- for(j=1;j <=ncovmodel-2;j++){
- if(Typevar[j]==0) {
- printf(" + V%d ",Tvar[j]);
- fprintf(ficlog," + V%d ",Tvar[j]);
- }else if(Typevar[j]==1) {
- printf(" + V%d*age ",Tvar[j]);
- fprintf(ficlog," + V%d*age ",Tvar[j]);
- }else if(Typevar[j]==2) {
- printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
- fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
- }else if(Typevar[j]==3) {
- printf(" + V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
- fprintf(ficlog," + V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
- }
- }
- printf("\n");
- fprintf(ficlog,"\n");
- for(i=1,jk=1; i <=nlstate; i++){
- for(k=1; k <=(nlstate+ndeath); k++){
- if (k != i) {
- printf("%d%d ",i,k);
- fprintf(ficlog,"%d%d ",i,k);
- for(j=1; j <=ncovmodel; j++){
- printf("%12.7f ",p[jk]);
- fprintf(ficlog,"%12.7f ",p[jk]);
- jk++;
- }
- printf("\n");
- fprintf(ficlog,"\n");
- }
- }
- }
- /* minfit ( n, vsmall, v, d ); */
- /* v is overwritten with R. */
- /*
- Heuristic numbers:
- If the axes may be badly scaled (which is to be avoided if
- possible), then set SCBD = 10. Otherwise set SCBD = 1.
-
- If the problem is known to be ill-conditioned, initialize ILLC = true.
- KTM is the number of iterations without improvement before the
- algorithm terminates. KTM = 4 is very cautious; usually KTM = 1
- is satisfactory.
- */
- double machep, small;
- double dmin;
- int illc=0; /* Local, int ILLC, is TRUE if the system is ill-conditioned. */
- machep = DBL_EPSILON;
- small = machep * machep;
- /* m2 = dsqrt(machep); */
-
- /*
- * Sort the eigenvalues and eigenvectors.
- */
- printf(" Sort the eigenvalues and eigenvectors....\n");
- svsort ( n, d, xi );
- printf("Sorted Eigenvalues:\n");
- for(i=1; i<=n;i++){
- printf(" d[%d]=%g",i,d[i]);
- }
- printf("\n");
- /*
- * Determine the smallest eigenvalue.
- */
- printf(" Determine the smallest eigenvalue.\n");
- dmin = fmax ( d[n], small );
- /*
- * The ratio of the smallest to largest eigenvalue determines whether
- * the system is ill conditioned.
- */
- if ( dmin < sqrt(machep) * d[1] ) { /* m2*d[0] */
- illc = 1;
- } else {
- illc = 0;
- }
- printf(" The ratio of the smallest to largest eigenvalue determines whether\n \
- the system is ill conditioned=%d . dmin=%.12lf < m2=%.12lf * d[1]=%.12lf \n",illc, dmin,sqrt(machep), d[1]);
- /* if ( 1.0 < scbd ) { */
- /* r8vec_print ( n, z, " The scale factors:" ); */
- /* } */
- /* r8vec_print ( n, d, " Principal values of the quadratic form:" ); */
- /* } */
- /* if ( 3 < prin ) { */
- /* r8mat_print ( n, n, v, " The principal axes:" ); */
- /* } */
- free_vector(d,1,n);
- /*
- * The main loop ends here.
- */
-
- /* if ( 0 < prin ) */
- /* { */
- /* r8vec_print ( n, x, " X:" ); */
- /* } */
- }
-#endif /* NOTMINFIT */
#ifdef LINMINORIGINAL
#else
for (j=1, flatd=0;j<=n;j++) {
free_vector(pt,1,n);
return;
#endif
- } /* endif(flatd >0) */
+ }
#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 POWELLNOTTRUECONJUGATE
- } /* end of t or directest negative */
+
+#ifdef DEBUG
+ printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
+ fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
+ for(j=1;j<=n;j++){
+ printf(" %lf",xit[j]);
+ fprintf(ficlog," %lf",xit[j]);
+ }
+ printf("\n");
+ fprintf(ficlog,"\n");
#endif
+ } /* end of t or directest negative */
+ printf(" Directest is positive, P_n-P_0 does not increase the conjugacy. n=%d\n",n);
+ fprintf(ficlog," Directest is positive, P_n-P_0 does not increase the conjugacy. n=%d\n",n);
#ifdef POWELLNOF3INFF1TEST
#else
} /* end if (fptt < fp) */
#endif
-#ifdef POWELLORIGINCONJUGATE
- } /* if (t < 0.0) */
-#endif
#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */
} /*NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done */
#else
first++;
}
- /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
+ /* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl,
+ * (int)age, (int)delaymax, (int)agefin, ncvloop,
+ * (int)age-(int)agefin); */
free_vector(min,1,nlstate);
free_vector(max,1,nlstate);
free_vector(meandiff,1,nlstate);
/* 0.51326036147820708, 0.48673963852179264} */
/* If we start from prlim again, prlim tends to a constant matrix */
- int i, ii,j,k, k1;
+ int i, ii,j, k1;
int first=0;
double *min, *max, *meandiff, maxmax,sumnew=0.;
/* double **matprod2(); */ /* test */
/* Computes the backward probability at age agefin, cov[2], and covariate combination 'ij'. In fact cov is already filled and x too.
* Call to pmij(cov and x), call to cross prevalence, sums and inverses, left multiply, and returns in **ps as well as **bmij.
*/
- int i, ii, j,k;
+ int ii, j;
- double **out, **pmij();
+ double **pmij();
double sumnew=0.;
double agefin;
double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */
*/
- int i, j, d, h, k, k1;
+ int i, j, d, h, k1;
double **out, cov[NCOVMAX+1];
double **newm;
double agexact;
- double agebegin, ageend;
+ /*double agebegin, ageend;*/
/* Hstepm could be zero and should return the unit matrix */
for (i=1;i<=nlstate+ndeath;i++)
The addresss of po (p3mat allocated to the dimension of nhstepm) should be stored for output
*/
- int i, j, d, h, k, k1;
+ int i, j, d, h, k1;
double **out, cov[NCOVMAX+1], **bmij();
double **newm, ***newmm;
double agexact;
- double agebegin, ageend;
+ /*double agebegin, ageend;*/
double **oldm, **savm;
newmm=po; /* To be saved */
double funcone( double *x)
{
/* Same as func but slower because of a lot of printf and if */
- int i, ii, j, k, mi, d, kk, kv=0, kf=0;
+ int i, ii, j, k, mi, d, kv=0, kf=0;
int ioffset=0;
int ipos=0,iposold=0,ncovv=0;
void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double []))
{
- int i,j,k, jk, jkk=0, iter=0;
+ int i,j, jkk=0, iter=0;
double **xi;
- double fret;
- double fretone; /* Only one call to likelihood */
+ /*double fret;*/
+ /*double fretone;*/ /* Only one call to likelihood */
/* char filerespow[FILENAMELENGTH];*/
- double * p1; /* Shifted parameters from 0 instead of 1 */
+ /*double * p1;*/ /* Shifted parameters from 0 instead of 1 */
#ifdef NLOPT
int creturn;
nlopt_opt opt;
}
powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
#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 ); */
-/* int prin=4; */
-/* double h0=0.25; */
+ int prin=1;
+ double h0=0.25;
+ double macheps;
+ double fmin;
+ macheps=pow(16.0,-13.0);
/* #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); */
+/* p1= (p+1); */ /* p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */
+printf("Praxis Gegenfurtner \n");
+fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog);
/* praxis ( ftol, h0, npar, prin, p1, func ); */
-/* printf("End Praxis\n"); */
+ /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */
+ fmin = praxis(ftol,macheps, h0, npar, prin, p, func);
+printf("End Praxis\n");
#endif /* FLATSUP */
#ifdef LINMINORIGINAL
int i, m, jk, j1, bool, z1,j, iv;
int mi; /* Effective wave */
int iage;
- double agebegin, ageend;
+ double agebegin; /*, ageend;*/
double **prop;
double posprop;
int popforecast, int mobilav, int prevfcast, int mobilavproj, int prevbcast, int estepm , \
double jprev1, double mprev1,double anprev1, double dateprev1, double dateprojd, double dateback1, \
double jprev2, double mprev2,double anprev2, double dateprev2, double dateprojf, double dateback2){
- int jj1, k1, i1, cpt, k4, nres;
+ int jj1, k1, cpt, nres;
/* In fact some results are already printed in fichtm which is open */
fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
<li><a href='#secondorder'>Result files (second order (variance)</a>\n \
*/
/* double anprojd, mprojd, jprojd; */
/* double anprojf, mprojf, jprojf; */
- int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+ int yearp, stepsize, hstepm, nhstepm, j, k, i, h, nres=0;
double agec; /* generic age */
- double agelim, ppij, yp,yp1,yp2;
- double *popeffectif,*popcount;
+ double agelim, ppij;
+ /*double *popcount;*/
double ***p3mat;
/* double ***mobaverage; */
char fileresf[FILENAMELENGTH];
anback2 year of end of backprojection (same day and month as back1).
prevacurrent and prev are prevalences.
*/
- int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
+ int yearp, stepsize, hstepm, nhstepm, j, k, i, h, nres=0;
double agec; /* generic age */
- double agelim, ppij, ppi, yp,yp1,yp2; /* ,jintmean,mintmean,aintmean;*/
- double *popeffectif,*popcount;
+ double agelim, ppij, ppi; /* ,jintmean,mintmean,aintmean;*/
+ /*double *popcount;*/
double ***p3mat;
/* double ***mobaverage; */
char fileresfb[FILENAMELENGTH];
char dirfileres[132],optfileres[132];
- int ng;
+ /*int ng;*/
/*#ifdef windows */
/*-------- data file ----------*/
FILE *fic;
char dummy[]=" ";
- int i=0, j=0, n=0, iv=0, v;
+ int i = 0, j = 0, n = 0, iv = 0;/* , v;*/
int lstra;
int linei, month, year,iout;
int noffset=0; /* This is the offset if BOM data file */
*/
/* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1, Tage[1]=2 */
{
- int i, j, k, ks, v;
+ int i, j, k, ks;/* , v;*/
int n,m;
int j1, k1, k11, k12, k2, k3, k4;
char modelsav[300];
int agelim;
int hstepm;
int nhstepm;
- int h, i, i1, j, k, k4, nres=0;
+ int h, i, i1, j, k, nres=0;
double agedeb;
double ***p3mat;
double fret;
double dum=0.; /* Dummy variable */
- double ***p3mat;
+ /* double*** p3mat;*/
/* double ***mobaverage; */
double wald;
char pathr[MAXLINE], pathimach[MAXLINE];
char *tok, *val; /* pathtot */
/* int firstobs=1, lastobs=10; /\* nobs = lastobs-firstobs declared globally ;*\/ */
- int c, h , cpt, c2;
+ int c, h; /* c2; */
int jl=0;
int i1, j1, jk, stepsize=0;
int count=0;
double ***delti3; /* Scale */
double *delti; /* Scale */
double ***eij, ***vareij;
- double **varpl; /* Variances of prevalence limits by age */
+ //double **varpl; /* Variances of prevalence limits by age */
double *epj, vepp;
getcwd(pathcd, size);
#endif
syscompilerinfo(0);
- printf("\nIMaCh prax version minfit %s, %s\n%s",version, copyright, fullversion);
+ printf("\nIMaCh prax version %s, %s\n%s",version, copyright, fullversion);
if(argc <=1){
printf("\nEnter the parameter file name: ");
if(!fgets(pathr,FILENAMELENGTH,stdin)){
/* Calculates basic frequencies. Computes observed prevalence at single age
and for any valid combination of covariates
and prints on file fileres'p'. */
- freqsummary(fileres, p, pstart, agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
+ freqsummary(fileres, p, pstart, (double)agemin, agemax, s, agev, nlstate, imx, Tvaraff, invalidvarcomb, nbcode, ncodemax,mint,anint,strstart, \
firstpass, lastpass, stepm, weightopt, model);
fprintf(fichtm,"\n");
free_ivector(TmodelInvind,1,NCOVMAX);
free_ivector(TmodelInvQind,1,NCOVMAX);
- free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/
+ /* free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /\* Could be elsewhere ?*\/ */
free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
/* free_imatrix(codtab,1,100,1,10); */
git://henry.ined.fr / .git / commitdiff