/* $Id$
$State$
$Log$
- Revision 1.357 2023/06/14 14:55:52 brouard
- * imach.c (Module): Testing if conjugate gradient could be quicker when lot of variables POWELLORIGINCONJUGATE
+ Revision 1.6 2024/04/24 21:10:29 brouard
+ Summary: First IMaCh version using Brent Praxis software based on Buckhardt and Gegenfürtner C codes
- Revision 1.356 2023/05/23 12:08:43 brouard
- Summary: 0.99r46
+ Revision 1.5 2023/10/09 09:10:01 brouard
+ Summary: trying to reconsider
- * imach.c (Module): Fixed PROB_r
+ Revision 1.4 2023/06/22 12:50:51 brouard
+ Summary: stil on going
- Revision 1.355 2023/05/22 17:03:18 brouard
- Summary: 0.99r46
-
- * imach.c (Module): In the ILK....txt file, the number of columns
- before the covariates values is dependent of the number of states (16+nlstate): 0.99r46
+ Revision 1.3 2023/06/22 11:28:07 brouard
+ *** empty log message ***
- Revision 1.354 2023/05/21 05:05:17 brouard
- Summary: Temporary change for imachprax
+ Revision 1.2 2023/06/22 11:22:40 brouard
+ Summary: with svd but not working yet
Revision 1.353 2023/05/08 18:48:22 brouard
*** empty log message ***
/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
/* #define FLATSUP *//* Suppresses directions where likelihood is flat */
-#define POWELLORIGINCONJUGATE /* Don't use conjugate but biggest decrease if valuable */
+/* #define POWELLORIGINCONJUGATE /\* Don't use conjugate but biggest decrease if valuable *\/ */
+/* #define NOTMINFIT */
#include <math.h>
#include <stdio.h>
/* $State$ */
#include "version.h"
char version[]=__IMACH_VERSION__;
-char copyright[]="Testing conjugate April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
+char copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
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 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. */
+/* */
+/*********************************************************************/
+
+#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);
+}
+
+static void sort() /* d and v in descending order */
+{
+ 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;
+ }
+ }
+ }
+}
+
+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 */
+ }
+ 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;
+ }
+
+ /* 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" );
+/*
+ Determine the range of the rows in this strip.
+*/
+ if ( 1 < ilo ){
+ i2lo = ilo;
+ }else{
+ i2lo = 1;
+ }
+ if ( m < ihi ){
+ i2hi = m;
+ }else{
+ i2hi = ihi;
+ }
+
+ for ( i = i2lo; i <= i2hi; i++ ){
+/*
+ Print out (up to) 5 entries in row I, that lie in the current strip.
+*/
+ /* 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" );
+ }
+ }
+
+ /* 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 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" );
+}
+
+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");
+ }
+
+
+/* #ifdef MSDOS */
+/* static double tflin[N]; */
+/* #endif */
+
+static double flin(double l, int j)
+/* double l; */
+{
+ 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++;
+
+#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;
+ }
+ /* 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 (*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 quad() /* look for a minimum along the curve q0, q1, q2 */
+{
+ 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
+}
+
+/* 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 *\/ */
+
+
+ int seed; /* added */
+ int biter=0;
+ double r;
+ double randbrent( int (*));
+ double s, sf;
+
+ 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 gegen */
/*************** powell ************************/
/*
curr_time = *localtime(&rcurr_time);
/* printf("\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); */
/* fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",*iter,*fret, fp-*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); */
- Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /* Big iteration, i.e on ncovmodel cycle */
+ /* Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /\* Big iteration, i.e on ncovmodel cycle *\/ */
+ Bigter=(*iter - (*iter-1) % n)/n +1; /* Big iteration, i.e on ncovmodel cycle */
printf("\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
fprintf(ficlog,"\nPowell iter=%d Big Iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,Bigter,*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
fprintf(ficrespow,"%d %d %.12f %d",*iter,Bigter, *fret,curr_time.tm_sec-start_time.tm_sec);
fprintf(ficlog," - if your program needs %d BIG iterations (%d iterations) to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
}
}
- for (i=1;i<=n;i++) { /* For each direction i */
- for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */
- fptt=(*fret);
+ 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. 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 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
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 */
#endif
- /* Outputs are fret(new point p) p is updated and xit rescaled */
+ /* Outputs are fret(new point p) p is updated and xit rescaled */
if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions */
- /* because that direction will be replaced unless the gain del is small */
- /* in comparison with the 'probable' gain, mu^2, with the last average direction. */
- /* Unless the n directions are conjugate some gain in the determinant may be obtained */
- /* with the new direction. */
- del=fabs(fptt-(*fret));
- ibig=i;
+ /* because that direction will be replaced unless the gain del is small */
+ /* in comparison with the 'probable' gain, mu^2, with the last average direction. */
+ /* Unless the n directions are conjugate some gain in the determinant may be obtained */
+ /* with the new direction. */
+ del=fabs(fptt-(*fret));
+ ibig=i;
}
#ifdef DEBUG
printf("%d %.12e",i,(*fret));
fprintf(ficlog,"%d %.12e",i,(*fret));
for (j=1;j<=n;j++) {
- xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
- printf(" x(%d)=%.12e",j,xit[j]);
- fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
+ xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
+ printf(" x(%d)=%.12e",j,xit[j]);
+ fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
}
for(j=1;j<=n;j++) {
- printf(" p(%d)=%.12e",j,p[j]);
- fprintf(ficlog," p(%d)=%.12e",j,p[j]);
+ printf(" p(%d)=%.12e",j,p[j]);
+ fprintf(ficlog," p(%d)=%.12e",j,p[j]);
}
printf("\n");
fprintf(ficlog,"\n");
} /* 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]);
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) */
+ 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];
- pt[j]=p[j];
- }
+ 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]);
+#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=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
/* Even if f3 <f1, directest can be negative and t >0 */
/* mu² and del² are equal when f3=f1 */
- /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */
- /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */
- /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0 */
- /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */
+ /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */
+ /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */
+ /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0 */
+ /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */
#ifdef NRCORIGINAL
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/
#else
if (t < 0.0) { /* Then we use it for new direction */
#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("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);
}
}
#endif
-#ifdef POWELLORIGINCONJUGATE
for (j=1;j<=n;j++) {
xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */
}
-#else
- for (j=1;j<=n-1;j++) {
- xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */
- xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */
- }
-#endif
#ifdef LINMINORIGINAL
#else
for (j=1, flatd=0;j<=n;j++) {
#endif
printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
- /* The minimization in direction $\xi_1$ gives $P_1$. From $P_1$ minimization in direction $\xi_2$ gives */
- /* $P_2$. Minimization of line $P_2-P_1$ gives new starting point $P^{(1)}_0$ and direction $\xi_1$ is dropped and replaced by second */
- /* direction $\xi_1^{(1)}=\xi_2$. Also second direction is replaced by new direction $\xi^{(1)}_2=P_2-P_0$. */
-
- /* At the second iteration, starting from $P_0^{(1)}$, minimization amongst $\xi^{(1)}_1$ gives point $P^{(1)}_1$. */
- /* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$. As $P^{(2)}_1$ and */
- /* $P^{(1)}_0$ are minimizing in the same direction $P^{(1)}_2 - P^{(1)}_1= P_2-P_0$, directions $P^{(1)}_2-P^{(1)}_0$ */
- /* and $P_2-P_0$ (parallel to $\xi$ and $\xi^c$) are conjugate. } */
-
#ifdef DEBUG
printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
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) */
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;
for (i=1;i<=npar;i++) /* Starting with canonical directions j=1,n xi[i=1,n][j] */
for (j=1;j<=npar;j++)
xi[i][j]=(i==j ? 1.0 : 0.0);
- printf("Powell\n"); fprintf(ficlog,"Powell\n");
+ printf("Powell-prax\n"); fprintf(ficlog,"Powell-prax\n");
strcpy(filerespow,"POW_");
strcat(filerespow,fileres);
if((ficrespow=fopen(filerespow,"w"))==NULL) {
}
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=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 Gegenfurtner \n");
+fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog);
+/* praxis ( ftol, h0, npar, prin, p1, func ); */
+ /* 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;
if(j==0) j=1; /* Survives at least one month after exam */
else if(j<0){
nberr++;
- printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
j=1; /* Temporary Dangerous patch */
printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
- fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
}
k=k+1;
/*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
if(j<0){
nberr++;
- printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
- fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ printf("Error! Negative delay (%d) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
+ fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld (around line %d) who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
}
sum=sum+j;
}
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 \
<img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Survival functions (period) in state j */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
+ fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages.
Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
}
/* State specific survival functions (period) */
for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\
- And probability to be observed in various states (up to %d) being in state %d at different ages. \
+ And probability to be observed in various states (up to %d) being in state %d at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. \
<a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
}
/* Period (forward stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
- fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
+ fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be
alive in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
}
/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
- from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \
- account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
+ from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d). Randomness in cross-sectional prevalence is not taken into \
+ account but can visually be appreciated. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_"));
fprintf(fichtm," <img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
*/
/* 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 version %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");
#ifdef GSL
printf("GSL optimization\n"); fprintf(ficlog,"Powell\n");
#else
- printf("Powell\n"); fprintf(ficlog,"Powell\n");
+ printf("Powell-mort\n"); fprintf(ficlog,"Powell-mort\n");
#endif
strcpy(filerespow,"POW-MORT_");
strcat(filerespow,fileresu);
for(i=1; i <=NDIM; i++)
for(j=i+1;j<=NDIM;j++)
- matcov[i][j]=matcov[j][i];
+ matcov[i][j]=matcov[j][i];
printf("\nCovariance matrix\n ");
fprintf(ficlog,"\nCovariance matrix\n ");
}
/* Results */
- /* Value of covariate in each resultine will be compututed (if product) and sorted according to model rank */
+ /* Value of covariate in each resultine will be computed (if product) and sorted according to model rank */
/* It is precov[] because we need the varying age in order to compute the real cov[] of the model equation */
precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1);
endishere=0;
/* */
if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
continue;
- printf("\n# model %s \n#****** Result for:", model); /* HERE model is empty */
- fprintf(ficrest,"\n# model %s \n#****** Result for:", model);
- fprintf(ficlog,"\n# model %s \n#****** Result for:", model);
+ printf("\n# model=1+age+%s \n#****** Result for:", model); /* HERE model is empty */
+ fprintf(ficrest,"\n# model=1+age+%s \n#****** Result for:", model);
+ fprintf(ficlog,"\n# model=1+age+%s \n#****** Result for:", model);
/* It might not be a good idea to mix dummies and quantitative */
/* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model *\/ */
for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
} /* mle==-3 arrives here for freeing */
/* endfree:*/
+ if(mle!=-3) free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);
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