--- imach/src/imachprax.c 2023/10/09 09:10:01 1.5 +++ imach/src/imachprax.c 2024/04/24 21:10:29 1.6 @@ -1,6 +1,9 @@ -/* $Id: imachprax.c,v 1.5 2023/10/09 09:10:01 brouard Exp $ +/* $Id: imachprax.c,v 1.6 2024/04/24 21:10:29 brouard Exp $ $State: Exp $ $Log: imachprax.c,v $ + 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.5 2023/10/09 09:10:01 brouard Summary: trying to reconsider @@ -1384,12 +1387,12 @@ double gnuplotversion=GNUPLOTVERSION; #define ODIRSEPARATOR '\\' #endif -/* $Id: imachprax.c,v 1.5 2023/10/09 09:10:01 brouard Exp $ */ +/* $Id: imachprax.c,v 1.6 2024/04/24 21:10:29 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; 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: 1.5 $ $Date: 2023/10/09 09:10:01 $"; +char fullversion[]="$Revision: 1.6 $ $Date: 2024/04/24 21:10:29 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1436,7 +1439,8 @@ int *wav; /* Number of waves for this in 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 */ @@ -2623,365 +2627,1512 @@ void linmin(double p[], double xi[], int free_vector(pcom,1,n); } -/**** praxis ****/ -# include +/**** praxis gegen ****/ -void transpose_in_place ( int n, double **a ) +/* 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. */ +/* */ +/*********************************************************************/ -/******************************************************************************/ -/* - Purpose: +#include +#include +#include +#include /* for DBL_EPSILON */ +/* #include "machine.h" */ - TRANSPOSE_IN_PLACE transposes a square matrix in place. - Licensing: - This code is distributed under the GNU LGPL license. - Input, int N, the number of rows and columns of the matrix A. +/* 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; - Input/output, double A[N*N], the matrix to be transposed. -*/ +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; - double t; - - /* for ( j = 0; j < n; j++ ){ */ - /* for ( i = 0; i < j; i++ ) { */ - for ( j = 1; j <= n; j++ ){ - for ( i = 1; i < j; i++ ) { - /* t = a[i+j*n]; */ - /* a[i+j*n] = a[j+i*n]; */ - /* a[j+i*n] = t; */ - t = a[i][j]; - a[i][j] = a[j][i]; - a[j][i] = t; + 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; } - } - return; -} -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 ( fabs ( x ) < fabs ( y ) ) { - a = fabs ( y ); - b = fabs ( x ); - } else { - a = fabs ( x ); - b = fabs ( y ); - } + if ( 1 < ilo ){ + i2lo = ilo; + }else{ + i2lo = 1; + } + if ( m < ihi ){ + i2hi = m; + }else{ + i2hi = ihi; + } + + 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 0) { /* linear search */ + /* for (i=0; i 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 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 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; i0.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 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= 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 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(); - for (j1 = 1; j1 < n; j1++) { +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 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 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 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 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 (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 2) + matprint(" NEW DIRECTIONS vectors:",v,n,n); + /* for (j=0; j 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 z[i]) + s = z[i]; + } + /* for (i=0; 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 1.0) { +#ifdef DEBUGPRAX + printf(" Unscale the axes.\n"); +#endif + /* for (i=0; 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 */ -/* void svdcmp(double **a, int m, int n, double w[], double **v) */ -void svdminfit(double **a, int m, int n, double w[]) -{ - /* From numerical recipes */ - /* Given a matrix a[1..m][1..n], this routine computes its singular value */ - /* decomposition, A = U ·W ·V T . The matrix U replaces a on output. */ - /* The diagonal matrix of singular values W is out- put as a vector w[1..n]. */ - /* The matrix V (not the transpose V T ) is output as v[1..n][1..n]. */ - - /* But in fact from Golub 1970 Algol60 */ - -/* Computation of the singular values and complete orthogonal decom- */ -/* 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 */ - - - 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; - - 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); -} +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 ************************/ /* @@ -3106,8 +4257,9 @@ void powell(double p[], double **xi, int } } 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); @@ -3115,7 +4267,7 @@ void powell(double p[], double **xi, int 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 @@ -3149,6 +4301,7 @@ void powell(double p[], double **xi, int } /* 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]); @@ -3171,6 +4324,28 @@ void powell(double p[], double **xi, int /* 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); @@ -3179,15 +4354,19 @@ void powell(double p[], double **xi, int 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 @@ -3217,194 +4396,59 @@ void powell(double p[], double **xi, int 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++) { @@ -3429,28 +4473,28 @@ void powell(double p[], double **xi, int 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 @@ -3655,7 +4699,9 @@ void powell(double p[], double **xi, int 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); @@ -3690,7 +4736,7 @@ void powell(double p[], double **xi, int /* 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 */ @@ -3957,9 +5003,9 @@ double **pmij(double **ps, double *cov, /* 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) */ @@ -4172,11 +5218,11 @@ double ***hpxij(double ***po, int nhstep */ - 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++) @@ -4353,11 +5399,11 @@ double ***hbxij(double ***po, int nhstep 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 */ @@ -4876,7 +5922,7 @@ double func( double *x) 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; @@ -5426,13 +6472,13 @@ void likelione(FILE *ficres,double p[], 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; @@ -5512,18 +6558,23 @@ void mlikeli(FILE *ficres,double p[], in } 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 @@ -6654,7 +7705,7 @@ void prevalence(double ***probs, double 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; @@ -8541,7 +9592,7 @@ void printinghtml(char fileresu[], char 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,"
  • Result files (first order: no variance)\n \
  • Result files (second order (variance)\n \ @@ -10405,10 +11456,10 @@ void prevforecast(char fileres[], double */ /* 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]; @@ -10556,10 +11607,10 @@ void prevforecast(char fileres[], double 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]; @@ -11242,7 +12293,7 @@ void printinggnuplotmort(char fileresu[] char dirfileres[132],optfileres[132]; - int ng; + /*int ng;*/ /*#ifdef windows */ @@ -11266,7 +12317,7 @@ int readdata(char datafile[], int firsto /*-------- 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 */ @@ -11883,7 +12934,7 @@ int decodemodel( char model[], int lasto */ /* 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]; @@ -13311,7 +14362,7 @@ int hPijx(double *p, int bage, int fage) 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; @@ -13517,7 +14568,7 @@ int main(int argc, char *argv[]) double fret; double dum=0.; /* Dummy variable */ - double ***p3mat; + /* double*** p3mat;*/ /* double ***mobaverage; */ double wald; @@ -13530,7 +14581,7 @@ int main(int argc, char *argv[]) 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; @@ -13565,7 +14616,7 @@ int main(int argc, char *argv[]) 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; @@ -13623,7 +14674,7 @@ int main(int argc, char *argv[]) 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)){ @@ -14525,7 +15576,7 @@ This file: %s
    Tit /* 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"); @@ -15835,7 +16886,7 @@ Please run with mle=-1 to get a correct 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); */