--- imach/src/imachprax.c 2023/01/31 09:24:19 1.1 +++ imach/src/imachprax.c 2023/06/22 11:22:40 1.2 @@ -1,8 +1,23 @@ -/* $Id: imachprax.c,v 1.1 2023/01/31 09:24:19 brouard Exp $ +/* $Id: imachprax.c,v 1.2 2023/06/22 11:22:40 brouard Exp $ $State: Exp $ $Log: imachprax.c,v $ - Revision 1.1 2023/01/31 09:24:19 brouard - Summary: version s1 with praxis instead of Powell for large models with age and difficulties to converge + 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 *** + + Revision 1.352 2023/04/29 10:46:21 brouard + *** empty log message *** + + Revision 1.351 2023/04/29 10:43:47 brouard + Summary: 099r45 + + Revision 1.350 2023/04/24 11:38:06 brouard + *** empty log message *** + + Revision 1.349 2023/01/31 09:19:37 brouard + Summary: Improvements in models with age*Vn*Vm Revision 1.347 2022/09/18 14:36:44 brouard Summary: version 0.99r42 @@ -1358,12 +1373,12 @@ double gnuplotversion=GNUPLOTVERSION; #define ODIRSEPARATOR '\\' #endif -/* $Id: imachprax.c,v 1.1 2023/01/31 09:24:19 brouard Exp $ */ +/* $Id: imachprax.c,v 1.2 2023/06/22 11:22:40 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="January 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.1 $ $Date: 2023/01/31 09:24:19 $"; +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.2 $ $Date: 2023/06/22 11:22:40 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1474,6 +1489,7 @@ extern time_t time(); struct tm start_time, end_time, curr_time, last_time, forecast_time; time_t rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */ +time_t rlast_btime; /* raw time */ struct tm tm; char strcurr[80], strfor[80]; @@ -1597,6 +1613,11 @@ int **nbcode, *Tvar; /**< model=V2 => Tv /* Tprod[i]=k 1 2 */ /* Position in model of the ith prod without age */ /* cptcovage 1 2 3 */ /* Counting cov*age in the model equation */ /* Tage[cptcovage]=k 5 8 10 */ /* Position in the model of ith cov*age */ +/* model="V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/ +/* p Tvard[1][1]@21 = {6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}*/ +/* p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 } */ +/* p Tvardk[1][1]@24 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0}*/ +/* p Tvardk[1][1]@22 = {0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0} */ /* Tvard[1][1]@4={4,3,1,2} V4*V3 V1*V2 */ /* Position in model of the ith prod without age */ /* Tvardk[4][1]=4;Tvardk[4][2]=3;Tvardk[7][1]=1;Tvardk[7][2]=2 */ /* Variables of a prod at position in the model equation*/ /* TvarF TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */ @@ -1795,6 +1816,20 @@ char *trimbb(char *out, char *in) return s; } +char *trimbtab(char *out, char *in) +{ /* Trim blanks or tabs in line but keeps first blanks if line starts with blanks */ + char *s; + s=out; + while (*in != '\0'){ + while( (*in == ' ' || *in == '\t')){ /* && *(in+1) != '\0'){*/ + in++; + } + *out++ = *in++; + } + *out='\0'; + return s; +} + /* char *substrchaine(char *out, char *in, char *chain) */ /* { */ /* /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */ @@ -2578,2583 +2613,259 @@ void linmin(double p[], double xi[], int } /**** praxis ****/ -# include -/* # include */ -/* # include */ -/* # include */ -/* # include */ -/* # include */ - -# include "praxis.h" - -/******************************************************************************/ - -double flin ( int n, int jsearch, double l, double (*func) ( double [] ), - double x[], int *nf, double v[], double q0[], double q1[], double *qd0, - double *qd1, double *qa, double *qb, double *qc ) -/* double flin ( int n, int jsearch, double l, double f ( double x[], int n ), */ -/* double x[], int *nf, double v[], double q0[], double q1[], double *qd0, */ -/* double *qd1, double *qa, double *qb, double *qc ) */ +double pythag( double x, double y ) /******************************************************************************/ /* Purpose: - - FLIN is the function of one variable to be minimized by MINNY. - - Discussion: - - F(X) is a scalar function of a vector argument X. - - A minimizer of F(X) is sought along a line or parabola. - + R8_HYPOT returns the value of sqrt ( X^2 + Y^2 ). Licensing: - This code is distributed under the GNU LGPL license. - Modified: - - 28 July 2016 - + 26 March 2012 Author: - - Original FORTRAN77 version by Richard Brent. - C version by John Burkardt. - - Reference: - - Richard Brent, - Algorithms for Minimization with Derivatives, - Prentice Hall, 1973, - Reprinted by Dover, 2002. - + John Burkardt Parameters: - - Input, int N, the number of variables. - - Input, int JSEARCH, indicates the kind of search. - If JSEARCH is a legal column index, linear search along V(*,JSEARCH). - If JSEARCH is -1, then the search is parabolic, based on X, Q0 and Q1. - - Input, double L, is the parameter determining the particular - point at which F is to be evaluated. - For a linear search, L is the step size. - For a quadratic search, L is a parameter which specifies - a point in the plane of X, Q0 and Q1. - - Input, double F ( double X[], int N ), the function to be minimized. - - Input, double X[N], the base point of the search. - - Input/output, int *NF, the function evaluation counter. - - Input, double V[N,N], a matrix whose columns constitute - search directions. - - Input, double Q0[N], Q1[N], two auxiliary points used to - determine the plane when a quadratic search is performed. - - Input, double *QD0, *QD1, values needed to compute the - coefficients QA, QB, QC. - - Output, double *QA, *QB, *QC, coefficients used to combine - Q0, X, and A1 if a quadratic search is used. - - Output, double FLIN, the value of the function at the - minimizing point. + Input, double X, Y, the arguments. + Output, double R8_HYPOT, the value of sqrt ( X^2 + Y^2 ). */ { - int i; - double *t; + double a; + double b; double value; - t = ( double * ) malloc ( n * sizeof ( double ) ); -/* - The search is linear. -*/ - if ( 0 <= jsearch ) - { - for ( i = 0; i < n; i++ ) - { - t[i] = x[i] + l * v[i+jsearch*n]; - } + if ( fabs ( x ) < fabs ( y ) ) { + a = fabs ( y ); + b = fabs ( x ); + } else { + a = fabs ( x ); + b = fabs ( y ); } /* - The search is along a parabolic space curve. + A contains the larger value. */ - else - { - *qa = l * ( l - *qd1 ) / ( *qd0 + *qd1 ) / *qd0; - *qb = - ( l + *qd0 ) * ( l - *qd1 ) / *qd1 / *qd0; - *qc = ( l + *qd0 ) * l / *qd1 / ( *qd0 + *qd1 ); - - for ( i = 0; i < n; i++ ) - { - t[i] = *qa * q0[i] + *qb * x[i] + *qc * q1[i]; - } + if ( a == 0.0 ) { + value = 0.0; + } else { + value = a * sqrt ( 1.0 + ( b / a ) * ( b / a ) ); } -/* - The function evaluation counter NF is incremented. -*/ - *nf = *nf + 1; -/* - Evaluate the function. -*/ - value = (*func) ( (t-1) );/* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/ - /* value = f ( t, n ); */ - - free ( t ); - return value; } -/******************************************************************************/ - -void minfit ( int n, double tol, double a[], double q[] ) - -/******************************************************************************/ -/* - Purpose: - - MINFIT computes the singular value decomposition of an N by N array. - - Discussion: - - This is an improved version of the EISPACK routine MINFIT - restricted to the case M = N and P = 0. - - The singular values of the array A are returned in Q. A is - overwritten with the orthogonal matrix V such that U * diag(Q) = A * V, - where U is another orthogonal matrix. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 30 July 2016 - - Author: - - Original FORTRAN77 version by Richard Brent. - C version by John Burkardt. - - Reference: - - Richard Brent, - Algorithms for Minimization with Derivatives, - Prentice Hall, 1973, - Reprinted by Dover, 2002. - - James Wilkinson, Christian Reinsch, - Handbook for Automatic Computation, - Volume II, Linear Algebra, Part 2, - Springer Verlag, 1971. - - Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, Yasuhiko Ikebe, - Virginia Klema, Cleve Moler, - Matrix Eigensystem Routines, EISPACK Guide, - Lecture Notes in Computer Science, Volume 6, - Springer Verlag, 1976, - ISBN13: 978-3540075462, - LC: QA193.M37. - - Parameters: - - Input, int N, the order of the matrix A. - - Input, double TOL, a tolerance which determines when a vector - (a column or part of a column of the matrix) may be considered - "essentially" equal to zero. - Input/output, double A[N,N]. On input, an N by N array whose - singular value decomposition is desired. On output, the - SVD orthogonal matrix factor V. - - Input/output, double Q[N], the singular values. -*/ +/* void svdcmp(double **a, int m, int n, double w[], double **v) */ +void svdcmp(double **a, int m, int n, double w[]) { - double c; - double *e; - double eps; - double f; - double g; - double h; - int i; - int ii; - int j; - int jj; - int k; - int kt; - const int kt_max = 30; - int l; - int l2; - double s; - int skip; - double temp; - double x; - double y; - double z; -/* - Householder's reduction to bidiagonal form. -*/ - if ( n == 1 ) - { - q[0] = a[0+0*n]; - a[0+0*n] = 1.0; - return; - } - - e = ( double * ) malloc ( n * sizeof ( double ) ); - - eps = DBL_EPSILON; - g = 0.0; - x = 0.0; - - for ( i = 1; i <= n; i++ ) - { - e[i-1] = g; - l = i + 1; - - s = 0.0; - for ( ii = i; ii <= n; ii++ ) - { - s = s + a[ii-1+(i-1)*n] * a[ii-1+(i-1)*n]; - } - - g = 0.0; - - if ( tol <= s ) - { - f = a[i-1+(i-1)*n]; - - g = sqrt ( s ); - - if ( 0.0 <= f ) - { - g = - g; - } - - h = f * g - s; - a[i-1+(i-1)*n] = f - g; - - for ( j = l; j <= n; j++ ) - { - f = 0.0; - for ( ii = i; ii <= n; ii++ ) - { - f = f + a[ii-1+(i-1)*n] * a[ii-1+(j-1)*n]; - } - f = f / h; - - for ( ii = i; ii <= n; ii++ ) - { - a[ii-1+(j-1)*n] = a[ii-1+(j-1)*n] + f * a[ii-1+(i-1)*n]; - } - } - } - - q[i-1] = g; - - s = 0.0; - for ( j = l; j <= n; j++ ) - { - s = s + a[i-1+(j-1)*n] * a[i-1+(j-1)*n]; - } - - g = 0.0; - - if ( tol <= s ) - { - if ( i < n ) - { - f = a[i-1+i*n]; - } - - g = sqrt ( s ); - - if ( 0.0 <= f ) - { - g = - g; - } - - h = f * g - s; - - if ( i < n ) - { - a[i-1+i*n] = f - g; - for ( jj = l; jj <= n; jj++ ) - { - e[jj-1] = a[i-1+(jj-1)*n] / h; - } - - for ( j = l; j <= n; j++ ) - { - s = 0.0; - for ( jj = l; jj <= n; jj++ ) - { - s = s + a[j-1+(jj-1)*n] * a[i-1+(jj-1)*n]; - } - for ( jj = l; jj <= n; jj++ ) - { - a[j-1+(jj-1)*n] = a[j-1+(jj-1)*n] + s * e[jj-1]; - } - } - } - } - - y = fabs ( q[i-1] ) + fabs ( e[i-1] ); - - x = fmax ( x, y ); - } -/* - Accumulation of right-hand transformations. -*/ - a[n-1+(n-1)*n] = 1.0; - g = e[n-1]; - l = n; - - for ( i = n - 1; 1 <= i; i-- ) - { - if ( g != 0.0 ) - { - h = a[i-1+i*n] * g; - - for ( ii = l; ii <= n; ii++ ) - { - a[ii-1+(i-1)*n] = a[i-1+(ii-1)*n] / h; - } - - for ( j = l; j <= n; j++ ) - { - s = 0.0; - for ( jj = l; jj <= n; jj++ ) - { - s = s + a[i-1+(jj-1)*n] * a[jj-1+(j-1)*n]; - } - - for ( ii = l; ii <= n; ii++ ) - { - a[ii-1+(j-1)*n] = a[ii-1+(j-1)*n] + s * a[ii-1+(i-1)*n]; - } - } - } - - for ( jj = l; jj <= n; jj++ ) - { - a[i-1+(jj-1)*n] = 0.0; - } - - for ( ii = l; ii <= n; ii++ ) - { - a[ii-1+(i-1)*n] = 0.0; - } - - a[i-1+(i-1)*n] = 1.0; - - g = e[i-1]; - - l = i; - } -/* - Diagonalization of the bidiagonal form. -*/ - eps = eps * x; - - for ( k = n; 1 <= k; k-- ) - { - kt = 0; - - for ( ; ; ) - { - kt = kt + 1; - - if ( kt_max < kt ) - { - e[k-1] = 0.0; - fprintf ( stderr, "\n" ); - fprintf ( stderr, "MINFIT - Fatal error!\n" ); - fprintf ( stderr, " The QR algorithm failed to converge.\n" ); - exit ( 1 ); - } - - skip = 0; - - for ( l2 = k; 1 <= l2; l2-- ) - { - l = l2; - - if ( fabs ( e[l-1] ) <= eps ) - { - skip = 1; - break; - } - - if ( 1 < l ) - { - if ( fabs ( q[l-2] ) <= eps ) - { - break; - } - } - } -/* - Cancellation of E(L) if 1 < L. -*/ - if ( ! skip ) - { - c = 0.0; - s = 1.0; - - for ( i = l; i <= k; i++ ) - { - f = s * e[i-1]; - e[i-1] = c * e[i-1]; - if ( fabs ( f ) <= eps ) - { - break; - } - g = q[i-1]; -/* - q(i) = h = sqrt(g*g + f*f). -*/ - h = r8_hypot ( f, g ); - - q[i-1] = h; - - if ( h == 0.0 ) - { - g = 1.0; - h = 1.0; - } - - c = g / h; - s = - f / h; - } - } -/* - Test for convergence for this index K. -*/ - z = q[k-1]; - - if ( l == k ) - { - if ( z < 0.0 ) - { - q[k-1] = - z; - for ( i = 1; i <= n; i++ ) - { - a[i-1+(k-1)*n] = - a[i-1+(k-1)*n]; - } - } - break; - } -/* - Shift from bottom 2*2 minor. -*/ - x = q[l-1]; - y = q[k-2]; - g = e[k-2]; - h = e[k-1]; - f = ( ( y - z ) * ( y + z ) + ( g - h ) * ( g + h ) ) / ( 2.0 * h * y ); - - g = r8_hypot ( f, 1.0 ); - - if ( f < 0.0 ) - { - temp = f - g; - } - else - { - temp = f + g; - } - - f = ( ( x - z ) * ( x + z ) + h * ( y / temp - h ) ) / x; -/* - Next QR transformation. -*/ - c = 1.0; - s = 1.0; - - for ( i = l + 1; i <= k; i++ ) - { - g = e[i-1]; - y = q[i-1]; - h = s * g; - g = g * c; - - z = r8_hypot ( f, h ); - - e[i-2] = z; - - if ( z == 0.0 ) - { - f = 1.0; - z = 1.0; - } - - c = f / z; - s = h / z; - f = x * c + g * s; - g = - x * s + g * c; - h = y * s; - y = y * c; - - for ( j = 1; j <= n; j++ ) - { - x = a[j-1+(i-2)*n]; - z = a[j-1+(i-1)*n]; - a[j-1+(i-2)*n] = x * c + z * s; - a[j-1+(i-1)*n] = - x * s + z * c; - } - - z = r8_hypot ( f, h ); - - q[i-2] = z; - - if ( z == 0.0 ) - { - f = 1.0; - z = 1.0; - } - - c = f / z; - s = h / z; - f = c * g + s * y; - x = - s * g + c * y; - } - - e[l-1] = 0.0; - e[k-1] = f; - q[k-1] = x; - } - } - - free ( e ); - - return; -} -/******************************************************************************/ - -void minny ( int n, int jsearch, int nits, double *d2, double *x1, double *f1, - int fk, double (*func) ( double []), double x[], double t, double h, - double v[], double q0[], double q1[], int *nl, int *nf, double dmin, - double ldt, double *fx, double *qa, double *qb, double *qc, double *qd0, - double *qd1 ) -/* void minny ( int n, int jsearch, int nits, double *d2, double *x1, double *f1, */ -/* int fk, double f ( double x[], int n ), double x[], double t, double h, */ -/* double v[], double q0[], double q1[], int *nl, int *nf, double dmin, */ -/* double ldt, double *fx, double *qa, double *qb, double *qc, double *qd0, */ -/* double *qd1 ) */ - -/******************************************************************************/ -/* - Purpose: - - MINNY minimizes a scalar function of N variables along a line. - - Discussion: - - MINNY minimizes F along the line from X in the direction V(*,JSEARCH) - or else using a quadratic search in the plane defined by Q0, Q1 and X. - - If FK = true, then F1 is FLIN(X1). Otherwise X1 and F1 are ignored - on entry unless final FX is greater than F1. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 03 August 2016 - - Author: - - Original FORTRAN77 version by Richard Brent. - C version by John Burkardt. - - Reference: - - Richard Brent, - Algorithms for Minimization with Derivatives, - Prentice Hall, 1973, - Reprinted by Dover, 2002. - - Parameters: - - Input, int N, the number of variables. - - Input, int JSEARCH, indicates the kind of search. - If J is a legal columnindex, linear search in the direction of V(*,JSEARCH). - Otherwise, the search is parabolic, based on X, Q0 and Q1. - - Input, int NITS, the maximum number of times the interval - may be halved to retry the calculation. - - Input/output, double *D2, is either zero, or an approximation to - the value of (1/2) times the second derivative of F. - - Input/output, double *X1, on entry, an estimate of the - distance from X to the minimum along V(*,JSEARCH), or a curve. - On output, the distance between X and the minimizer that was found. - - Input/output, double *F1, ? - - Input, int FK; if FK is TRUE, then on input F1 contains - the value FLIN(X1). - - Input, double F ( double X[], int N ), is the name of the function to - be minimized. - - Input/output, double X[N], ? - - Input, double T, ? - - Input, double H, ? - - Input, double V[N,N], a matrix whose columns are direction - vectors along which the function may be minimized. - - ?, double Q0[N], ? - - ?, double Q1[N], ? - - Input/output, int *NL, the number of linear searches. - - Input/output, int *NF, the number of function evaluations. - - Input, double DMIN, an estimate for the smallest eigenvalue. - - Input, double LDT, the length of the step. - - Input/output, double *FX, the value of F(X,N). - - Input/output, double *QA, *QB, *QC; - - Input/output, double *QD0, *QD1, ?. -*/ -{ - double d1; - int dz; - double f0; - double f2; - double fm; - int i; - int k; - double m2; - double m4; - double machep; - int ok; - double s; - double sf1; - double small; - double sx1; - double t2; - double temp; - double x2; - double xm; - - machep = DBL_EPSILON; - small = machep * machep; - m2 = sqrt ( machep ); - m4 = sqrt ( m2 ); - sf1 = *f1; - sx1 = *x1; - k = 0; - xm = 0.0; - fm = *fx; - f0 = *fx; - dz = ( *d2 < machep ); -/* - Find the step size. -*/ - s = r8vec_norm ( n, x ); - - if ( dz ) - { - temp = dmin; - } - else - { - temp = *d2; - } - - t2 = m4 * sqrt ( fabs ( *fx ) / temp + s * ldt ) + m2 * ldt; - s = m4 * s + t; - if ( dz && s < t2 ) - { - t2 = s; - } - - t2 = fmax ( t2, small ); - t2 = fmin ( t2, 0.01 * h ); - - if ( fk && *f1 <= fm ) - { - xm = *x1; - fm = *f1; - } - - if ( ( ! fk ) || fabs ( *x1 ) < t2 ) - { - if ( 0.0 <= *x1 ) - { - temp = 1.0; - } - else - { - temp = - 1.0; - } - - *x1 = temp * t2; - *f1 = flin ( n, jsearch, *x1, func, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); - /* *f1 = flin ( n, jsearch, *x1, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); */ - } - - if ( *f1 <= fm ) - { - xm = *x1; - fm = *f1; - } -/* - Evaluate FLIN at another point and estimate the second derivative. -*/ - for ( ; ; ) - { - if ( dz ) - { - if ( *f1 <= f0 ) - { - x2 = 2.0 * *x1; - } - else - { - x2 = - *x1; - } - - f2 = flin ( n, jsearch, x2, func, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); - /* f2 = flin ( n, jsearch, x2, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); */ - - if ( f2 <= fm ) - { - xm = x2; - fm = f2; - } - - *d2 = ( x2 * ( *f1 - f0 ) - *x1 * ( f2 - f0 ) ) - / ( ( *x1 * x2 ) * ( *x1 - x2 ) ); - } -/* - Estimate the first derivative at 0. -*/ - d1 = ( *f1 - f0 ) / *x1 - *x1 * *d2; - dz = 1; -/* - Predict the minimum. -*/ - if ( *d2 <= small ) - { - if ( 0.0 <= d1 ) - { - x2 = - h; - } - else - { - x2 = h; - } - } - else - { - x2 = ( - 0.5 * d1 ) / *d2; - } - - if ( h < fabs ( x2 ) ) - { - if ( x2 <= 0.0 ) - { - x2 = - h; - } - else - { - x2 = h; - } - } -/* - Evaluate F at the predicted minimum. -*/ - ok = 1; - - for ( ; ; ) - { - f2 = flin ( n, jsearch, x2, func, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); - /* f2 = flin ( n, jsearch, x2, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); */ - - if ( nits <= k || f2 <= f0 ) - { - break; - } - - k = k + 1; - - if ( f0 < *f1 && 0.0 < *x1 * x2 ) - { - ok = 0; - break; - } - x2 = 0.5 * x2; - } - - if ( ok ) - { - break; - } - } -/* - Increment the one-dimensional search counter. -*/ - *nl = *nl + 1; - - if ( fm < f2 ) - { - x2 = xm; - } - else - { - fm = f2; - } -/* - Get a new estimate of the second derivative. -*/ - if ( small < fabs ( x2 * ( x2 - *x1 ) ) ) - { - *d2 = ( x2 * ( *f1 - f0 ) - *x1 * ( fm - f0 ) ) - / ( ( *x1 * x2 ) * ( *x1 - x2 ) ); - } - else - { - if ( 0 < k ) - { - *d2 = 0.0; - } - } - - *d2 = fmax ( *d2, small ); - - *x1 = x2; - *fx = fm; - - if ( sf1 < *fx ) - { - *fx = sf1; - *x1 = sx1; - } -/* - Update X for linear search. -*/ - if ( 0 <= jsearch ) - { - for ( i = 0; i < n; i++ ) - { - x[i] = x[i] + *x1 * v[i+jsearch*n]; - } - } - - return; -} -/******************************************************************************/ - -/* double praxis ( double t0, double h0, int n, int prin, double x[], */ -/* double f ( double x[], int n ) ) */ -double praxis ( double t0, double h0, int n, int prin, double x[], - double (*func) ( double [] )) - -/******************************************************************************/ -/* - Purpose: - - PRAXIS seeks an N-dimensional minimizer X of a scalar function F(X). - - Discussion: - - PRAXIS returns the minimum of the function F(X,N) of N variables - using the principal axis method. The gradient of the function is - not required. - - The approximating quadratic form is - - Q(x") = F(x,n) + (1/2) * (x"-x)" * A * (x"-x) - - where X is the best estimate of the minimum and - - A = inverse(V") * D * inverse(V) - - V(*,*) is the matrix of search directions; - D(*) is the array of second differences. - - If F(X) has continuous second derivatives near X0, then A will tend - to the hessian of F at X0 as X approaches X0. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 03 August 2016 - - Author: - - Original FORTRAN77 version by Richard Brent. - C version by John Burkardt. - - Reference: - - Richard Brent, - Algorithms for Minimization with Derivatives, - Prentice Hall, 1973, - Reprinted by Dover, 2002. - - Parameters: - - Input, double T0, is a tolerance. PRAXIS attempts to return - praxis = f(x) such that if X0 is the true local minimum near X, then - norm ( x - x0 ) < T0 + sqrt ( EPSILON ) * norm ( X ), - where EPSILON is the machine precision. - - Input, double H0, is the maximum step size. H0 should be - set to about the maximum distance from the initial guess to the minimum. - If H0 is set too large or too small, the initial rate of - convergence may be slow. - - Input, int N, the number of variables. - - Input, int PRIN, controls printing intermediate results. - 0, nothing is printed. - 1, F is printed after every n+1 or n+2 linear minimizations. - final X is printed, but intermediate X is printed only - if N is at most 4. - 2, the scale factors and the principal values of the approximating - quadratic form are also printed. - 3, X is also printed after every few linear minimizations. - 4, the principal vectors of the approximating quadratic form are - also printed. - - Input/output, double X[N], is an array containing on entry a - guess of the point of minimum, on return the estimated point of minimum. - - Input, double F ( double X[], int N ), is the name of the function to be - minimized. - - Output, double PRAXIS, the function value at the minimizer. - - Local parameters: - - Local, double DMIN, an estimate for the smallest eigenvalue. - - Local, double FX, the value of F(X,N). - - Local, int ILLC, is TRUE if the system is ill-conditioned. - - Local, double LDT, the length of the step. - - Local, int NF, the number of function evaluations. - - Local, int NL, the number of linear searches. -*/ -{ - int biter=0; /* Added to count the loops */ - double *d; - double d2; - double df; - double dmin; - double dn; - double dni; - double f1; - int fk; - double fx; - double h; - int i; - int illc; - int j; - int jsearch; - int k; - int k2; - int kl; - int kt; - int ktm; - double large; - double ldfac; - double lds; - double ldt; - double m2; - double m4; - double machep; - int nits; - int nl; - int nf; - double *q0; - double *q1; - double qa; - double qb; - double qc; - double qd0; - double qd1; - double qf1; - double r; - double s; - double scbd; - int seed; - double sf; - double sl; - double small; - double t; - double temp; - double t2; - double *v; - double value; - double vlarge; - double vsmall; - double *y; - double *z; -/* - Allocation. -*/ - d = ( double * ) malloc ( n * sizeof ( double ) ); - q0 = ( double * ) malloc ( n * sizeof ( double ) ); - q1 = ( double * ) malloc ( n * sizeof ( double ) ); - v = ( double * ) malloc ( n * n * sizeof ( double ) ); - y = ( double * ) malloc ( n * sizeof ( double ) ); - z = ( double * ) malloc ( n * sizeof ( double ) ); -/* - Initialization. -*/ - machep = DBL_EPSILON; - small = machep * machep; - vsmall = small * small; - large = 1.0 / small; - vlarge = 1.0 / vsmall; - m2 = sqrt ( machep ); - m4 = sqrt ( m2 ); - seed = 123456789; -/* - 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. -*/ - scbd = 1.0; - illc = 0; - ktm = 1; - - if ( illc ) - { - ldfac = 0.1; - } - else - { - ldfac = 0.01; - } - - kt = 0; - nl = 0; - nf = 1; - /* fx = f ( x, n ); */ - fx = (*func) ( (x-1) );/* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/ - qf1 = fx; - t = small + fabs ( t0 ); - t2 = t; - dmin = small; - h = h0; - h = fmax ( h, 100.0 * t ); - ldt = h; -/* - The initial set of search directions V is the identity matrix. -*/ - for ( j = 0; j < n; j++ ) - { - for ( i = 0; i < n; i++ ) - { - v[i+j*n] = 0.0; - } - v[j+j*n] = 1.0; - } - - for ( i = 0; i < n; i++ ) - { - d[i] = 0.0; - } - qa = 0.0; - qb = 0.0; - qc = 0.0; - qd0 = 0.0; - qd1 = 0.0; - r8vec_copy ( n, x, q0 ); - r8vec_copy ( n, x, q1 ); - - if ( 0 < prin ) - { - print2 ( n, x, prin, fx, nf, nl ); - } -/* - The main loop starts here. -*/ - for ( ; ; ) - { - biter++; /* Added to count the loops */ - printf("\n Big iteration %d \n",biter); - sf = d[0]; - d[0] = 0.0; -/* - Minimize along the first direction V(*,1). -*/ - jsearch = 0; - nits = 2; - d2 = d[0]; - s = 0.0; - value = fx; - fk = 0; - - minny ( n, jsearch, nits, &d2, &s, &value, fk, func, x, t, - h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); - /* minny ( n, jsearch, nits, &d2, &s, &value, fk, func, x, t, */ - /* h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); */ - - d[0] = d2; - - if ( s <= 0.0 ) - { - for ( i = 0; i < n; i++ ) - { - v[i+0*n] = - v[i+0*n]; - } - } - - if ( sf <= 0.9 * d[0] || d[0] <= 0.9 * sf ) - { - for ( i = 1; i < n; i++ ) - { - d[i] = 0.0; - } - } -/* - The inner loop starts here. -*/ - for ( k = 2; k <= n; k++ ) - { - r8vec_copy ( n, x, y ); - - sf = fx; - - if ( 0 < kt ) - { - illc = 1; - } - - for ( ; ; ) - { - kl = k; - df = 0.0; -/* - A random step follows, to avoid resolution valleys. -*/ - if ( illc ) - { - for ( j = 0; j < n; j++ ) - { - r = r8_uniform_01 ( &seed ); - s = ( 0.1 * ldt + t2 * pow ( 10.0, kt ) ) * ( r - 0.5 ); - z[j] = s; - for ( i = 0; i < n; i++ ) - { - x[i] = x[i] + s * v[i+j*n]; - } - } - - fx = (*func) ( (x-1) );/* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/ - /* fx = f ( x, n ); */ - nf = nf + 1; - } -/* - Minimize along the "non-conjugate" directions V(*,K),...,V(*,N). -*/ - for ( k2 = k; k2 <= n; k2++ ) - { - sl = fx; - - jsearch = k2 - 1; - nits = 2; - d2 = d[k2-1]; - s = 0.0; - value = fx; - fk = 0; - - minny ( n, jsearch, nits, &d2, &s, &value, fk, func, x, t, - h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); - /* minny ( n, jsearch, nits, &d2, &s, &value, fk, f, x, t, */ - /* h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); */ - - d[k2-1] = d2; - - if ( illc ) - { - s = d[k2-1] * pow ( s + z[k2-1], 2 ); - } - else - { - s = sl - fx; - } - - if ( df <= s ) - { - df = s; - kl = k2; - } - } -/* - If there was not much improvement on the first try, set - ILLC = true and start the inner loop again. -*/ - if ( illc ) - { - break; - } - printf("\n fabs( 100.0 * machep(=%.12lf) * fx(=%.12lf) ) <=? df(=%.12lf)\n", machep, fx, df); - if ( fabs ( 100.0 * machep * fx ) <= df ) - { - break; - } - illc = 1; - } - - if ( k == 2 && 1 < prin ) - { - r8vec_print ( n, d, " The second difference array:" ); - } -/* - Minimize along the "conjugate" directions V(*,1),...,V(*,K-1). -*/ - for ( k2 = 1; k2 < k; k2++ ) - { - jsearch = k2 - 1; - nits = 2; - d2 = d[k2-1]; - s = 0.0; - value = fx; - fk = 0; - - minny ( n, jsearch, nits, &d2, &s, &value, fk, func, x, t, - h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); - /* minny ( n, jsearch, nits, &d2, &s, &value, fk, f, x, t, */ - /* h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); */ - - d[k2-1] = d2; - } - - f1 = fx; - fx = sf; - - for ( i = 0; i < n; i++ ) - { - temp = x[i]; - x[i] = y[i]; - y[i] = temp - y[i]; - } - - lds = r8vec_norm ( n, y ); -/* - 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 ( small < lds ) - { - for ( j = kl - 1; k <= j; j-- ) - { - for ( i = 1; i <= n; i++ ) - { - v[i-1+j*n] = v[i-1+(j-1)*n]; - } - d[j] = d[j-1]; - } - - d[k-1] = 0.0; - - for ( i = 1; i <= n; i++ ) - { - v[i-1+(k-1)*n] = y[i-1] / lds; - } -/* - Minimize along the new "conjugate" direction V(*,k), which is - the normalized vector: (new x) - (old x). -*/ - jsearch = k - 1; - nits = 4; - d2 = d[k-1]; - value = f1; - fk = 1; - - minny ( n, jsearch, nits, &d2, &lds, &value, fk, func, x, t, - h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); - /* minny ( n, jsearch, nits, &d2, &lds, &value, fk, f, x, t, */ - /* h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); */ - - d[k-1] = d2; - - if ( lds <= 0.0 ) - { - lds = - lds; - for ( i = 1; i <= n; i++ ) - { - v[i-1+(k-1)*n] = - v[i-1+(k-1)*n]; - } - } - } - - ldt = ldfac * ldt; - ldt = fmax ( ldt, lds ); - - if ( 0 < prin ) - { - printf(" k=%d",k); - print2 ( n, x, prin, fx, nf, nl ); - } - - t2 = r8vec_norm ( n, x ); - - t2 = m2 * t2 + t; -/* - See whether the length of the step taken since starting the - inner loop exceeds half the tolerance. -*/ - if ( 0.5 * t2 < ldt ) - { - kt = - 1; - } - - kt = kt + 1; - - if ( ktm < kt ) - { - if ( 0 < prin ) - { - r8vec_print ( n, x, " X:" ); - } - - free ( d ); - free ( q0 ); - free ( q1 ); - free ( v ); - free ( y ); - free ( z ); - - return fx; - } - } -/* - The inner loop ends here. - - Try quadratic extrapolation in case we are in a curved valley. -*/ - quad ( n, func, x, t, h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qf1, - &qa, &qb, &qc, &qd0, &qd1 ); - /* quad ( n, f, x, t, h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qf1, */ - /* &qa, &qb, &qc, &qd0, &qd1 ); */ - - for ( j = 0; j < n; j++ ) - { - d[j] = 1.0 / sqrt ( d[j] ); - } - - dn = r8vec_max ( n, d ); - - if ( 3 < prin ) - { - r8mat_print ( n, n, v, " The new direction vectors:" ); - } - - for ( j = 0; j < n; j++ ) - { - for ( i = 0; i < n; i++ ) - { - v[i+j*n] = ( d[j] / dn ) * v[i+j*n]; - } - } -/* - Scale the axes to try to reduce the condition number. -*/ - if ( 1.0 < scbd ) - { - for ( i = 0; i < n; i++ ) - { - s = 0.0; - for ( j = 0; j < n; j++ ) - { - s = s + v[i+j*n] * v[i+j*n]; - } - s = sqrt ( s ); - z[i] = fmax ( m4, s ); - } - - s = r8vec_min ( n, z ); - - for ( i = 0; i < n; i++ ) - { - sl = s / z[i]; - z[i] = 1.0 / sl; - - if ( scbd < z[i] ) - { - sl = 1.0 / scbd; - z[i] = scbd; - } - for ( j = 0; j < n; j++ ) - { - v[i+j*n] = sl * v[i+j*n]; - } - } - } -/* - Calculate a new set of orthogonal directions before repeating - the main loop. - - Transpose V for MINFIT: -*/ - printf(" Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n"); - r8mat_transpose_in_place ( n, v ); -/* - 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(" 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"); - minfit ( n, vsmall, v, d ); -/* - Unscale the axes. -*/ - printf(" Unscale the axes.\n"); - if ( 1.0 < scbd ) - { - for ( i = 0; i < n; i++ ) - { - for ( j = 0; j < n; j++ ) - { - v[i+j*n] = z[i] * v[i+j*n]; - } - } - - for ( j = 0; j < n; j++ ) - { - s = 0.0; - for ( i = 0; i < n; i++ ) - { - s = s + v[i+j*n] * v[i+j*n]; - } - s = sqrt ( s ); - - d[j] = s * d[j]; - for ( i = 0; i < n; i++ ) - { - v[i+j*n] = v[i+j*n] / s; - } - } - } - - for ( i = 0; i < n; i++ ) - { - dni = dn * d[i]; - - if ( large < dni ) - { - d[i] = vsmall; - } - else if ( dni < small ) - { - d[i] = vlarge; - } - else - { - d[i] = 1.0 / dni / dni; - } - } -/* - Sort the eigenvalues and eigenvectors. -*/ - printf(" Sort the eigenvalues and eigenvectors....\n"); - svsort ( n, d, v ); -/* - Determine the smallest eigenvalue. -*/ - printf(" Determine the smallest eigenvalue.\n"); - dmin = fmax ( d[n-1], small ); -/* - The ratio of the smallest to largest eigenvalue determines whether - the system is ill conditioned. -*/ - - if ( dmin < 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[0]=%.12lf \n",illc, dmin,m2, d[0]); - - if ( 1 < prin ) - { - 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:" ); - } -/* - The main loop ends here. -*/ - } - - if ( 0 < prin ) - { - r8vec_print ( n, x, " X:" ); - } -/* - Free memory. -*/ - free ( d ); - free ( q0 ); - free ( q1 ); - free ( v ); - free ( y ); - free ( z ); - - return fx; -} -/******************************************************************************/ - -void print2 ( int n, double x[], int prin, double fx, int nf, int nl ) - -/******************************************************************************/ -/* - Purpose: - - PRINT2 prints certain data about the progress of the iteration. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 28 July 2016 - - Author: - - Original FORTRAN77 version by Richard Brent. - C version by John Burkardt. - - Reference: - - Richard Brent, - Algorithms for Minimization with Derivatives, - Prentice Hall, 1973, - Reprinted by Dover, 2002. - - Parameters: - - Input, int N, the number of variables. - - Input, double X[N], the current estimate of the minimizer. - - Input, int PRIN, the user-specifed print level. - 0, nothing is printed. - 1, F is printed after every n+1 or n+2 linear minimizations. - final X is printed, but intermediate X is printed only - if N is at most 4. - 2, the scale factors and the principal values of the approximating - quadratic form are also printed. - 3, X is also printed after every few linear minimizations. - 4, the principal vectors of the approximating quadratic form are - also printed. - - Input, double FX, the smallest value of F(X) found so far. - - Input, int NF, the number of function evaluations. - - Input, int NL, the number of linear searches. -*/ -{ - 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 ); - - if ( n <= 4 || 2 < prin ) - { - r8vec_print ( n, x, " X:" ); - } - - return; -} -/******************************************************************************/ - -void quad ( int n, double (*func) ( double [] ), double x[], double t, - double h, double v[], double q0[], double q1[], int *nl, int *nf, double dmin, - double ldt, double *fx, double *qf1, double *qa, double *qb, double *qc, - double *qd0, double *qd1 ) -/* void quad ( int n, double f ( double x[], int n ), double x[], double t, */ -/* double h, double v[], double q0[], double q1[], int *nl, int *nf, double dmin, */ -/* double ldt, double *fx, double *qf1, double *qa, double *qb, double *qc, */ -/* double *qd0, double *qd1 ) */ - -/******************************************************************************/ -/* - Purpose: - - QUAD seeks to minimize the scalar function F along a particular curve. - - Discussion: - - The minimizer to be sought is required to lie on a curve defined - by Q0, Q1 and X. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 30 July 2016 - - Author: - - Original FORTRAN77 version by Richard Brent. - C version by John Burkardt. - - Reference: - - Richard Brent, - Algorithms for Minimization with Derivatives, - Prentice Hall, 1973, - Reprinted by Dover, 2002. - - Parameters: - - Input, int N, the number of variables. - - Input, double F ( double X[], int N ), the name of the function to - be minimized. - - Input/output, double X[N], ? - - Input, double T, ? - - Input, double H, ? - - Input, double V[N,N], the matrix of search directions. - - Input/output, double Q0[N], Q1[N], two auxiliary points used to define - a curve through X. - - Input/output, int *NL, the number of linear searches. - - Input/output, int *NF, the number of function evaluations. - - Input, double DMIN, an estimate for the smallest eigenvalue. - - Input, double LDT, the length of the step. - - Input/output, double *FX, the value of F(X,N). - - Input/output, double *QF1, *QA, *QB, *QC, *QD0, *QD1 ? -*/ -{ - int fk; - int i; - int jsearch; - double l; - int nits; - double s; - double temp; - double value; - - temp = *fx; - *fx = *qf1; - *qf1 = temp; - - for ( i = 0; i < n; i++ ) - { - temp = x[i]; - x[i] = q1[i]; - q1[i] = temp; - } - - *qd1 = 0.0; - for ( i = 0; i < n; i++ ) - { - *qd1 = *qd1 + ( x[i] - q1[i] ) * ( x[i] - q1[i] ); - } - *qd1 = sqrt ( *qd1 ); - - if ( *qd0 <= 0.0 || *qd1 <= 0.0 || *nl < 3 * n * n ) - { - *fx = *qf1; - *qa = 0.0; - *qb = 0.0; - *qc = 1.0; - s = 0.0; - } - else - { - jsearch = - 1; - nits = 2; - s = 0.0; - l = *qd1; - value = *qf1; - fk = 1; - - minny ( n, jsearch, nits, &s, &l, &value, fk, func, x, t, - h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 ); - /* minny ( n, jsearch, nits, &s, &l, &value, fk, f, x, t, */ - /* h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 ); */ - - *qa = l * ( l - *qd1 ) / ( *qd0 + *qd1 ) / *qd0; - *qb = - ( l + *qd0 ) * ( l - *qd1 ) / *qd1 / *qd0; - *qc = ( l + *qd0 ) * l / *qd1 / ( *qd0 + *qd1 ); - } - - *qd0 = *qd1; - - for ( i = 0; i < n; i++ ) - { - s = q0[i]; - q0[i] = x[i]; - x[i] = *qa * s + *qb * x[i] + *qc * q1[i]; - } - - return; -} -/******************************************************************************/ - -double r8_hypot ( double x, double y ) - -/******************************************************************************/ -/* - 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 ). -*/ -{ - double a; - double b; - double value; - - if ( fabs ( x ) < fabs ( y ) ) - { - a = fabs ( y ); - b = fabs ( x ); - } - else - { - a = fabs ( x ); - b = fabs ( y ); - } -/* - A contains the larger value. -*/ - if ( a == 0.0 ) - { - value = 0.0; - } - else - { - value = a * sqrt ( 1.0 + ( b / a ) * ( b / a ) ); - } - - return value; -} -/******************************************************************************/ - -double r8_uniform_01 ( int *seed ) - -/******************************************************************************/ -/* - Purpose: - - R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1]. - - Discussion: - - This routine implements the recursion - - seed = 16807 * seed mod ( 2^31 - 1 ) - r8_uniform_01 = seed / ( 2^31 - 1 ) - - The integer arithmetic never requires more than 32 bits, - including a sign bit. - - If the initial seed is 12345, then the first three computations are - - Input Output R8_UNIFORM_01 - SEED SEED - - 12345 207482415 0.096616 - 207482415 1790989824 0.833995 - 1790989824 2035175616 0.947702 - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 11 August 2004 - - Author: - - John Burkardt - - Reference: - - Paul Bratley, Bennett Fox, Linus Schrage, - A Guide to Simulation, - Springer Verlag, pages 201-202, 1983. - - Pierre L'Ecuyer, - Random Number Generation, - in Handbook of Simulation - edited by Jerry Banks, - Wiley Interscience, page 95, 1998. - - Bennett Fox, - Algorithm 647: - Implementation and Relative Efficiency of Quasirandom - Sequence Generators, - ACM Transactions on Mathematical Software, - Volume 12, Number 4, pages 362-376, 1986. - - P A Lewis, A S Goodman, J M Miller, - A Pseudo-Random Number Generator for the System/360, - IBM Systems Journal, - Volume 8, pages 136-143, 1969. - - Parameters: - - Input/output, int *SEED, the "seed" value. Normally, this - value should not be 0. On output, SEED has been updated. - - Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between - 0 and 1. -*/ -{ - const int i4_huge = 2147483647; - int k; - double r; - - if ( *seed == 0 ) - { - fprintf ( stderr, "\n" ); - fprintf ( stderr, "R8_UNIFORM_01 - Fatal error!\n" ); - fprintf ( stderr, " Input value of SEED = 0\n" ); - exit ( 1 ); - } - - k = *seed / 127773; - - *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; - - if ( *seed < 0 ) - { - *seed = *seed + i4_huge; - } - - r = ( ( double ) ( *seed ) ) * 4.656612875E-10; - - return r; -} -/******************************************************************************/ - -void r8mat_print ( int m, int n, double a[], char *title ) - -/******************************************************************************/ -/* - Purpose: - - R8MAT_PRINT prints an R8MAT. - - Discussion: - - An R8MAT is a doubly dimensioned array of R8 values, stored as a vector - in column-major order. - - Entry A(I,J) is stored as A[I+J*M] - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 28 May 2008 - - Author: - - John Burkardt - - Parameters: - - Input, int M, the number of rows in A. - - Input, int N, the number of columns in A. - - Input, double A[M*N], the M by N matrix. - - Input, char *TITLE, a title. -*/ -{ - r8mat_print_some ( m, n, a, 1, 1, m, n, title ); - - return; -} -/******************************************************************************/ - -void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, - int jhi, char *title ) - -/******************************************************************************/ -/* - Purpose: - - R8MAT_PRINT_SOME prints some of an R8MAT. - - Discussion: - - An R8MAT is a doubly dimensioned array of R8 values, stored as a vector - in column-major order. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 26 June 2013 - - Author: - - John Burkardt - - Parameters: - - Input, int M, the number of rows of the matrix. - M must be positive. - - Input, int N, the number of columns of the matrix. - N must be positive. - - Input, double A[M*N], the matrix. - - Input, int ILO, JLO, IHI, JHI, designate the first row and - column, and the last row and column to be printed. - - Input, char *TITLE, a title. -*/ -{ -# define INCX 5 - - int i; - int i2hi; - int i2lo; - int j; - int j2hi; - int j2lo; - - fprintf ( stdout, "\n" ); - fprintf ( stdout, "%s\n", title ); - - if ( m <= 0 || n <= 0 ) - { - fprintf ( stdout, "\n" ); - fprintf ( stdout, " (None)\n" ); - return; - } -/* - Print the columns of the matrix, in strips of 5. -*/ - for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) - { - j2hi = j2lo + INCX - 1; - if ( n < j2hi ) - { - j2hi = n; - } - if ( jhi < j2hi ) - { - j2hi = jhi; - } - - fprintf ( stdout, "\n" ); -/* - For each column J in the current range... - - Write the header. -*/ - fprintf ( stdout, " Col: "); - for ( j = j2lo; j <= j2hi; j++ ) - { - fprintf ( stdout, " %7d ", j - 1 ); - } - fprintf ( stdout, "\n" ); - fprintf ( stdout, " Row\n" ); - fprintf ( stdout, "\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 ( stdout, "%5d:", i - 1 ); - for ( j = j2lo; j <= j2hi; j++ ) - { - fprintf ( stdout, " %14g", a[i-1+(j-1)*m] ); - } - fprintf ( stdout, "\n" ); - } - } - - return; -# undef INCX -} -/******************************************************************************/ - -void r8mat_transpose_in_place ( int n, double a[] ) - -/******************************************************************************/ -/* - Purpose: - - R8MAT_TRANSPOSE_IN_PLACE transposes a square matrix in place. - - Discussion: - - An R8MAT is a doubly dimensioned array of R8 values, stored as a vector - in column-major order. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 26 June 2008 - - Author: - - John Burkardt - - Parameters: - - Input, int N, the number of rows and columns of the matrix A. - - Input/output, double A[N*N], the matrix to be transposed. -*/ -{ - int i; - int j; - double t; - - for ( j = 0; j < n; j++ ) - { - for ( i = 0; i < j; i++ ) - { - t = a[i+j*n]; - a[i+j*n] = a[j+i*n]; - a[j+i*n] = t; - } - } - return; -} -/******************************************************************************/ - -void r8vec_copy ( int n, double a1[], double a2[] ) - -/******************************************************************************/ -/* - Purpose: - - R8VEC_COPY copies an R8VEC. - - Discussion: - - An R8VEC is a vector of R8's. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 03 July 2005 - - Author: - - John Burkardt - - Parameters: - - Input, int N, the number of entries in the vectors. - - Input, double A1[N], the vector to be copied. - - Input, double A2[N], the copy of A1. -*/ -{ - int i; - - for ( i = 0; i < n; i++ ) - { - a2[i] = a1[i]; - } - return; -} -/******************************************************************************/ - -double r8vec_max ( int n, double r8vec[] ) - -/******************************************************************************/ -/* - Purpose: - - R8VEC_MAX returns the value of the maximum element in a R8VEC. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 05 May 2006 - - Author: - - John Burkardt - - Parameters: - - Input, int N, the number of entries in the array. - - Input, double R8VEC[N], a pointer to the first entry of the array. - - Output, double R8VEC_MAX, the value of the maximum element. This - is set to 0.0 if N <= 0. -*/ -{ - int i; - double value; - - if ( n <= 0 ) - { - value = 0.0; - return value; - } - - value = r8vec[0]; - - for ( i = 1; i < n; i++ ) - { - if ( value < r8vec[i] ) - { - value = r8vec[i]; - } - } - return value; -} -/******************************************************************************/ - -double r8vec_min ( int n, double r8vec[] ) - -/******************************************************************************/ -/* - Purpose: - - R8VEC_MIN returns the value of the minimum element in a R8VEC. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 05 May 2006 - - Author: - - John Burkardt - - Parameters: - - Input, int N, the number of entries in the array. - - Input, double R8VEC[N], the array to be checked. - - Output, double R8VEC_MIN, the value of the minimum element. -*/ -{ - int i; - double value; - - value = r8vec[0]; - - for ( i = 1; i < n; i++ ) - { - if ( r8vec[i] < value ) - { - value = r8vec[i]; - } - } - return value; -} -/******************************************************************************/ - -double r8vec_norm ( int n, double a[] ) - -/******************************************************************************/ -/* - Purpose: - - R8VEC_NORM returns the L2 norm of an R8VEC. - - Discussion: - - The vector L2 norm is defined as: - - R8VEC_NORM = sqrt ( sum ( 1 <= I <= N ) A(I)^2 ). - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 01 March 2003 - - Author: - - John Burkardt - - Parameters: - - Input, int N, the number of entries in A. - - Input, double A[N], the vector whose L2 norm is desired. - - Output, double R8VEC_NORM, the L2 norm of A. -*/ -{ - int i; - double v; - - v = 0.0; - - for ( i = 0; i < n; i++ ) - { - v = v + a[i] * a[i]; - } - v = sqrt ( v ); - - return v; -} -/******************************************************************************/ - -void r8vec_print ( int n, double a[], char *title ) - -/******************************************************************************/ -/* - Purpose: - - R8VEC_PRINT prints an R8VEC. - - Discussion: - - An R8VEC is a vector of R8's. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 08 April 2009 - - Author: - - John Burkardt - - Parameters: - - Input, int N, the number of components of the vector. - - Input, double A[N], the vector to be printed. - - Input, char *TITLE, a title. -*/ -{ - int i,j, jk, k; - - double *p; - - p=(a-1); /* So that a[0]=p[1] */ - /* for (i=1;i<=n;i++) { */ - /* fprintf(ficrespow," %.12lf", p[i]); */ - /* } */ - /* fprintf(ficrespow,"\n");fflush(ficrespow); */ - printf("\n#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"); -/* printf("12 47.0114589 0.0154322 33.2424412 0.3279905 2.3731903 */ -/* 13 -21.5392400 0.1118147 1.2680506 1.2973408 -1.0663662 */ - 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"); - } - } - } - /* fprintf ( stdout, "\n" ); */ - /* fprintf ( stdout, " %s\n", title ); */ - fprintf ( stdout, " %s", title ); - /* fprintf ( stdout, "\n" ); */ - for ( i = 0; i < n; i++ ) - { - /* fprintf ( stdout, " %8d: %14g", i+1, a[i] ); */ - fprintf ( stdout, " %.12lf", a[i] ); - } - fprintf ( stdout, "\n" ); - /* for ( i = 0; i < n; i++ ) */ - /* { */ - /* fprintf ( stdout, " %8d: %14g\n", i, a[i] ); */ - /* } */ - - return; -} -/******************************************************************************/ - -void svsort ( int n, double d[], double v[] ) - -/******************************************************************************/ -/* - Purpose: - - SVSORT descending sorts D and adjusts the corresponding columns of V. - - Discussion: - - A simple bubble sort is used on D. - - In our application, D contains singular values, and the columns of V are - the corresponding right singular vectors. - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 28 July 2016 - - Author: - - Original FORTRAN77 version by Richard Brent. - C version by John Burkardt. - - Reference: - - Richard Brent, - Algorithms for Minimization with Derivatives, - Prentice Hall, 1973, - Reprinted by Dover, 2002. - - Parameters: - - Input, int N, the length of D, and the order of V. - - Input/output, double D[N], the vector to be sorted. - On output, the entries of D are in descending order. - - Input/output, double V[N,N], an N by N array to be adjusted - as D is sorted. In particular, if the value that was in D(I) on input is - moved to D(J) on output, then the input column V(*,I) is moved to - the output column V(*,J). -*/ -{ - int i; - int j1; - int j2; - int j3; - double t; - - for ( j1 = 0; j1 < n - 1; j1++ ) - { -/* - Find J3, the index of the largest entry in D[J1:N-1]. - MAXLOC apparently requires its output to be an array. -*/ - j3 = j1; - for ( j2 = j1 + 1; j2 < n; j2++ ) - { - if ( d[j3] < d[j2] ) - { - j3 = j2; - } - } -/* - If J1 != J3, swap D[J1] and D[J3], and columns J1 and J3 of V. -*/ - if ( j1 != j3 ) - { - t = d[j1]; - d[j1] = d[j3]; - d[j3] = t; - for ( i = 0; i < n; i++ ) - { - t = v[i+j1*n]; - v[i+j1*n] = v[i+j3*n]; - v[i+j3*n] = t; - } - } - } - - return; -} -/******************************************************************************/ - -void timestamp ( ) - -/******************************************************************************/ -/* - Purpose: - - TIMESTAMP prints the current YMDHMS date as a time stamp. - - Example: - - 31 May 2001 09:45:54 AM - - Licensing: - - This code is distributed under the GNU LGPL license. - - Modified: - - 24 September 2003 - - Author: - - John Burkardt - - Parameters: - - None -*/ -{ -# define TIME_SIZE 40 - - static char time_buffer[TIME_SIZE]; - const struct tm *tm; - time_t now; - - now = time ( NULL ); - tm = localtime ( &now ); - - strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); - - fprintf ( stdout, "%s\n", time_buffer ); + /* 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 dsvdcmp 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); +} - return; -# undef TIME_SIZE -} /* end praxis */ /*************** powell ************************/ @@ -5187,7 +2898,8 @@ void powell(double p[], double **xi, int double fp,fptt; double *xits; int niterf, itmp; - + int Bigter=0, nBigterf=1; + pt=vector(1,n); ptt=vector(1,n); xit=vector(1,n); @@ -5200,14 +2912,16 @@ void powell(double p[], double **xi, int ibig=0; del=0.0; rlast_time=rcurr_time; + rlast_btime=rcurr_time; /* (void) gettimeofday(&curr_time,&tzp); */ rcurr_time = time(NULL); 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); */ - printf("\nPowell iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,*fret,fp-*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); - fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f gain=%.3lg %ld sec. %ld sec.",*iter,*fret,fp-*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); -/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ + Bigter=(*iter - *iter % ncovmodel)/ncovmodel +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); fp=(*fret); /* From former iteration or initial value */ for (i=1;i<=n;i++) { fprintf(ficrespow," %.12lf", p[i]); @@ -5262,15 +2976,17 @@ void powell(double p[], double **xi, int strcurr[itmp-1]='\0'; printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); - for(niterf=10;niterf<=30;niterf+=10){ + for(nBigterf=1;nBigterf<=31;nBigterf+=10){ + niterf=nBigterf*ncovmodel; + /* rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); */ rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); forecast_time = *localtime(&rforecast_time); strcpy(strfor,asctime(&forecast_time)); itmp = strlen(strfor); if(strfor[itmp-1]=='\n') strfor[itmp-1]='\0'; - printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); - fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); + printf(" - 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); + 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 */ @@ -6568,94 +4284,46 @@ double func( double *x) for(k=1; k<=nlstate; k++) ll[k]=0.; ioffset=0; - for (i=1,ipmx=0, sw=0.; i<=imx; i++){ - /* Computes the values of the ncovmodel covariates of the model - depending if the covariates are fixed or varying (age dependent) and stores them in cov[] - Then computes with function pmij which return a matrix p[i][j] giving the elementary probability - to be observed in j being in i according to the model. - */ - ioffset=2+nagesqr ; - /* Fixed */ - for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummy or quant or prod */ - /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */ - /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ - /* TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */ - /* TvarFind; TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */ - cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/ - /* V1*V2 (7) TvarFind[2]=7, TvarFind[3]=9 */ - } - /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] - is 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6 - has been calculated etc */ - /* For an individual i, wav[i] gives the number of effective waves */ - /* We compute the contribution to Likelihood of each effective transition - mw[mi][i] is real wave of the mi th effectve wave */ - /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i]; - s2=s[mw[mi+1][i]][i]; - And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] because now is moved after nvocol+nqv - But if the variable is not in the model TTvar[iv] is the real variable effective in the model: - meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] - */ - for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */ + if(mle==1){ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + /* Computes the values of the ncovmodel covariates of the model + depending if the covariates are fixed or varying (age dependent) and stores them in cov[] + Then computes with function pmij which return a matrix p[i][j] giving the elementary probability + to be observed in j being in i according to the model. + */ + ioffset=2+nagesqr ; + /* Fixed */ + for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummy or quant or prod */ + /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */ + /* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + /* TvarF[1]=Tvar[6]=2, TvarF[2]=Tvar[7]=7, TvarF[3]=Tvar[9]=1 ID of fixed covariates or product V2, V1*V2, V1 */ + /* TvarFind; TvarFind[1]=6, TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod) */ + cov[ioffset+TvarFind[kf]]=covar[Tvar[TvarFind[kf]]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (TvarFind[1]=6)*/ + /* V1*V2 (7) TvarFind[2]=7, TvarFind[3]=9 */ + } + /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] + is 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6 + has been calculated etc */ + /* For an individual i, wav[i] gives the number of effective waves */ + /* We compute the contribution to Likelihood of each effective transition + mw[mi][i] is real wave of the mi th effectve wave */ + /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i]; + s2=s[mw[mi+1][i]][i]; + And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] because now is moved after nvocol+nqv + But if the variable is not in the model TTvar[iv] is the real variable effective in the model: + meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i] + */ + for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */ /* Wave varying (but not age varying) */ - /* for(k=1; k <= ncovv ; k++){ /\* Varying covariates in the model (single and product but no age )"V5+V4+V3+V4*V3+V5*age+V1*age+V1" +TvarVind 1,2,3,4(V4*V3) Tvar[1]@7{5, 4, 3, 6, 5, 1, 1 ; 6 because the created covar is after V5 and is 6, minus 1+1, 3,2,1,4 positions in cotvar*\/ */ - /* /\* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; but where is the crossproduct? *\/ */ - /* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; */ - /* } */ - for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Varying covariates (single and product but no age )*/ - itv=TvarVV[ncovv]; /* TvarVV={3, 1, 3} gives the name of each varying covariate */ - ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/ - if(FixedV[itv]!=0){ /* Not a fixed covariate */ - cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i]; /* cotvar[wav][ncovcol+nqv+iv][i] */ - }else{ /* fixed covariate */ - cotvarv=covar[itv][i]; /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */ - } - if(ipos!=iposold){ /* Not a product or first of a product */ - cotvarvold=cotvarv; - }else{ /* A second product */ - cotvarv=cotvarv*cotvarvold; - } - iposold=ipos; - cov[ioffset+ipos]=cotvarv; - } - /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)*\/ */ - /* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff *\/ */ - /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ - /* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model *\/ */ - /* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */ - /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ - /* } */ - /* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates *\/ */ - /* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff *\/ */ - /* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); *\/ */ - /* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */ - /* } */ - /* for products of time varying to be done */ - for (ii=1;ii<=nlstate+ndeath;ii++) - for (j=1;j<=nlstate+ndeath;j++){ - oldm[ii][j]=(ii==j ? 1.0 : 0.0); - savm[ii][j]=(ii==j ? 1.0 : 0.0); - } - - agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */ - ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */ - for(d=0; d 1 the results are less biased than in previous versions. - */ - s1=s[mw[mi][i]][i]; - s2=s[mw[mi+1][i]][i]; - bbh=(double)bh[mi][i]/(double)stepm; - /* bias bh is positive if real duration - * is higher than the multiple of stepm and negative otherwise. - */ - /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ - if( s2 > nlstate){ - /* i.e. if s2 is a death state and if the date of death is known - then the contribution to the likelihood is the probability to - die between last step unit time and current step unit time, - which is also equal to probability to die before dh - minus probability to die before dh-stepm . - In version up to 0.92 likelihood was computed - as if date of death was unknown. Death was treated as any other - health state: the date of the interview describes the actual state - and not the date of a change in health state. The former idea was - to consider that at each interview the state was recorded - (healthy, disable or death) and IMaCh was corrected; but when we - introduced the exact date of death then we should have modified - the contribution of an exact death to the likelihood. This new - contribution is smaller and very dependent of the step unit - stepm. It is no more the probability to die between last interview - and month of death but the probability to survive from last - interview up to one month before death multiplied by the - probability to die within a month. Thanks to Chris - Jackson for correcting this bug. Former versions increased - mortality artificially. The bad side is that we add another loop - which slows down the processing. The difference can be up to 10% - lower mortality. - */ - /* If, at the beginning of the maximization mostly, the - cumulative probability or probability to be dead is - constant (ie = 1) over time d, the difference is equal to - 0. out[s1][3] = savm[s1][3]: probability, being at state - s1 at precedent wave, to be dead a month before current - wave is equal to probability, being at state s1 at - precedent wave, to be dead at mont of the current - wave. Then the observed probability (that this person died) - is null according to current estimated parameter. In fact, - it should be very low but not zero otherwise the log go to - infinity. - */ + /* But now since version 0.9 we anticipate for bias at large stepm. + * If stepm is larger than one month (smallest stepm) and if the exact delay + * (in months) between two waves is not a multiple of stepm, we rounded to + * the nearest (and in case of equal distance, to the lowest) interval but now + * we keep into memory the bias bh[mi][i] and also the previous matrix product + * (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the + * probability in order to take into account the bias as a fraction of the way + * from savm to out if bh is negative or even beyond if bh is positive. bh varies + * -stepm/2 to stepm/2 . + * For stepm=1 the results are the same as for previous versions of Imach. + * For stepm > 1 the results are less biased than in previous versions. + */ + s1=s[mw[mi][i]][i]; + s2=s[mw[mi+1][i]][i]; + bbh=(double)bh[mi][i]/(double)stepm; + /* bias bh is positive if real duration + * is higher than the multiple of stepm and negative otherwise. + */ + /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ + if( s2 > nlstate){ + /* i.e. if s2 is a death state and if the date of death is known + then the contribution to the likelihood is the probability to + die between last step unit time and current step unit time, + which is also equal to probability to die before dh + minus probability to die before dh-stepm . + In version up to 0.92 likelihood was computed + as if date of death was unknown. Death was treated as any other + health state: the date of the interview describes the actual state + and not the date of a change in health state. The former idea was + to consider that at each interview the state was recorded + (healthy, disable or death) and IMaCh was corrected; but when we + introduced the exact date of death then we should have modified + the contribution of an exact death to the likelihood. This new + contribution is smaller and very dependent of the step unit + stepm. It is no more the probability to die between last interview + and month of death but the probability to survive from last + interview up to one month before death multiplied by the + probability to die within a month. Thanks to Chris + Jackson for correcting this bug. Former versions increased + mortality artificially. The bad side is that we add another loop + which slows down the processing. The difference can be up to 10% + lower mortality. + */ + /* If, at the beginning of the maximization mostly, the + cumulative probability or probability to be dead is + constant (ie = 1) over time d, the difference is equal to + 0. out[s1][3] = savm[s1][3]: probability, being at state + s1 at precedent wave, to be dead a month before current + wave is equal to probability, being at state s1 at + precedent wave, to be dead at mont of the current + wave. Then the observed probability (that this person died) + is null according to current estimated parameter. In fact, + it should be very low but not zero otherwise the log go to + infinity. + */ /* #ifdef INFINITYORIGINAL */ /* lli=log(out[s1][s2] - savm[s1][s2]); */ /* #else */ @@ -6739,49 +4444,194 @@ double func( double *x) /* else */ /* lli=log(out[s1][s2] - savm[s1][s2]); */ /* #endif */ - lli=log(out[s1][s2] - savm[s1][s2]); + lli=log(out[s1][s2] - savm[s1][s2]); + + } else if ( s2==-1 ) { /* alive */ + for (j=1,survp=0. ; j<=nlstate; j++) + survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; + /*survp += out[s1][j]; */ + lli= log(survp); + } + /* else if (s2==-4) { */ + /* for (j=3,survp=0. ; j<=nlstate; j++) */ + /* survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */ + /* lli= log(survp); */ + /* } */ + /* else if (s2==-5) { */ + /* for (j=1,survp=0. ; j<=2; j++) */ + /* survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */ + /* lli= log(survp); */ + /* } */ + else{ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ + } + /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ + /*if(lli ==000.0)*/ + /* printf("num[i], i=%d, bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + /* if (lli < log(mytinydouble)){ */ + /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */ + /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */ + /* } */ + } /* end of wave */ + } /* end of individual */ + } else if(mle==2){ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + ioffset=2+nagesqr ; + for (k=1; k<=ncovf;k++) + cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i]; + for(mi=1; mi<= wav[i]-1; mi++){ + for(k=1; k <= ncovv ; k++){ + cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ + } + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } + for(d=0; d<=dh[mi][i]; d++){ + newm=savm; + agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; + cov[2]=agexact; + if(nagesqr==1) + cov[3]= agexact*agexact; + for (kk=1; kk<=cptcovage;kk++) { + cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; + } + out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, + 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); + savm=oldm; + oldm=newm; + } /* end mult */ + + s1=s[mw[mi][i]][i]; + s2=s[mw[mi+1][i]][i]; + bbh=(double)bh[mi][i]/(double)stepm; + lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + } /* end of wave */ + } /* end of individual */ + } else if(mle==3){ /* exponential inter-extrapolation */ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; + for(mi=1; mi<= wav[i]-1; mi++){ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } + for(d=0; d1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + } /* end of wave */ + } /* end of individual */ + }else if (mle==4){ /* ml=4 no inter-extrapolation */ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; + for(mi=1; mi<= wav[i]-1; mi++){ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } + for(d=0; d(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ - } else if(mle==2){ - lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ - } else if(mle==3){ /* exponential inter-extrapolation */ - lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ - } else if (mle==4){ /* mle=4 no inter-extrapolation */ - lli=log(out[s1][s2]); /* Original formula */ - } else{ /* mle=0 back to 1 */ - lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ - /*lli=log(out[s1][s2]); */ /* Original formula */ - } - /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ - /*if(lli ==000.0)*/ - /* printf("num[i], i=%d, bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ - ipmx +=1; - sw += weight[i]; - ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; - /* if (lli < log(mytinydouble)){ */ - /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */ - /* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */ - /* } */ - } /* end of wave */ - } /* end of individual */ + out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, + 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); + savm=oldm; + oldm=newm; + } /* end mult */ + + s1=s[mw[mi][i]][i]; + s2=s[mw[mi+1][i]][i]; + if( s2 > nlstate){ + lli=log(out[s1][s2] - savm[s1][s2]); + } else if ( s2==-1 ) { /* alive */ + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= log(survp); + }else{ + lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ + } + ipmx +=1; + sw += weight[i]; + ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; + /* printf("num[i]=%09ld, i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */ + } /* end of wave */ + } /* end of individual */ + }else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ + for (i=1,ipmx=0, sw=0.; i<=imx; i++){ + for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; + for(mi=1; mi<= wav[i]-1; mi++){ + for (ii=1;ii<=nlstate+ndeath;ii++) + for (j=1;j<=nlstate+ndeath;j++){ + oldm[ii][j]=(ii==j ? 1.0 : 0.0); + savm[ii][j]=(ii==j ? 1.0 : 0.0); + } + for(d=0; d- Probability p%dj by origin %d and destination j. Dot's sizes are related to corresponding weight: %s-p%dj.png
\n \ \n",k,k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k); for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */ - /* kvar=Tvar[TvarFind[kf]]; */ /* variable */ - fprintf(fichtm,"
- Probability p%dj by origin %d and destination j with colored covariate V%d. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): %s-p%dj.png
\ -",k,k,Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,subdirf2(optionfilefiname,"ILK_"),k,Tvar[TvarFind[kf]]); + kvar=Tvar[TvarFind[kf]]; /* variable */ + fprintf(fichtm,"
- Probability p%dj by origin %d and destination j with colored covariate V%d. Same dot size of all points but with a different color for transitions with dummy variable V%d=1 at beginning of transition (keeping former color for V%d=0): ",k,k,Tvar[TvarFind[kf]],Tvar[TvarFind[kf]],Tvar[TvarFind[kf]]); + fprintf(fichtm,"%s-p%dj-%d.png
",subdirf2(optionfilefiname,"ILK_"),k,kvar,subdirf2(optionfilefiname,"ILK_"),k,kvar); + fprintf(fichtm,"",subdirf2(optionfilefiname,"ILK_"),k,Tvar[TvarFind[kf]]); } for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Loop on the time varying extended covariates (with extension of Vn*Vm */ ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */ @@ -10107,7 +7888,7 @@ void varprob(char optionfilefiname[], do double ***varpij; strcpy(fileresprob,"PROB_"); - strcat(fileresprob,fileres); + strcat(fileresprob,fileresu); if((ficresprob=fopen(fileresprob,"w"))==NULL) { printf("Problem with resultfile: %s\n", fileresprob); fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob); @@ -10926,7 +8707,9 @@ void printinggnuplot(char fileresu[], ch for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */ kvar=Tvar[TvarFind[kf]]; /* variable name */ /* k=18+Tvar[TvarFind[kf]];/\*offset because there are 18 columns in the ILK_ file but could be placed else where *\/ */ - k=18+kf;/*offset because there are 18 columns in the ILK_ file */ + /* k=18+kf;/\*offset because there are 18 columns in the ILK_ file *\/ */ + /* k=19+kf;/\*offset because there are 19 columns in the ILK_ file *\/ */ + k=16+nlstate+kf;/*offset because there are 19 columns in the ILK_ file, first cov Vn on col 21 with 4 living states */ for (i=1; i<= nlstate ; i ++) { fprintf(ficgp,"\nset out \"%s-p%dj-%d.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i,kvar); fprintf(ficgp,"unset log;\n# For each simple dummy covariate of the model \n plot \"%s\"",subdirf(fileresilk)); @@ -11612,18 +9395,22 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp," u %d:(",ioffset); kl=0; strcpy(gplotcondition,"("); - for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */ + /* for (k=1; k<=cptcoveff; k++){ /\* For each covariate writing the chain of conditions *\/ */ /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */ - lv=codtabm(k1,TnsdVar[Tvaraff[k]]); + for (k=1; k<=cptcovs; k++){ /* For each covariate k get corresponding value lv for combination k1 */ + /* lv=codtabm(k1,TnsdVar[Tvaraff[k]]); */ + lv=Tvresult[nres][k]; + vlv=TinvDoQresult[nres][Tvresult[nres][k]]; /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */ - vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; + /* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */ kl++; - sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); + /* sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */ + sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,lv, kl+1, vlv ); kl++; - if(k 1) + if(k 1) sprintf(gplotcondition+strlen(gplotcondition)," && "); } strcpy(gplotcondition+strlen(gplotcondition),")"); @@ -11707,7 +9494,8 @@ set ter svg size 640, 480\nunset log y\n }else{ fprintf(ficgp,",\\\n '' "); } - if(cptcoveff ==0){ /* No covariate */ + /* if(cptcoveff ==0){ /\* No covariate *\/ */ + if(cptcovs ==0){ /* No covariate */ ioffset=2; /* Age is in 2 */ /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ @@ -11731,7 +9519,7 @@ set ter svg size 640, 480\nunset log y\n /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ iyearc=ioffset-1; iagec=ioffset; - fprintf(ficgp," u %d:(",ioffset); /* PROBLEM HERE VERIFY */ + fprintf(ficgp," u %d:(",ioffset); kl=0; strcpy(gplotcondition,"("); for (k=1; k<=cptcovs; k++){ /* For each covariate k of the resultline, get corresponding value lv for combination k1 */ @@ -11819,7 +9607,8 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n"); fprintf(ficgp,"#model=1+age+%s \n",model); fprintf(ficgp,"# Type of graphic ng=%d\n",ng); - fprintf(ficgp,"# k1=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */ + /* fprintf(ficgp,"# k1=1 to 2^%d=%d\n",cptcoveff,m);/\* to be checked *\/ */ + fprintf(ficgp,"# k1=1 to 2^%d=%d\n",cptcovs,m);/* to be checked */ /* for(k1=1; k1 <=m; k1++) /\* For each combination of covariate *\/ */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ /* k1=nres; */ @@ -11948,20 +9737,20 @@ set ter svg size 640, 480\nunset log y\n if(cptcovdageprod >0){ /* if(j==Tprod[ijp]) { */ /* not necessary */ /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */ - if(ijp <=cptcovprod) { /* Product */ - if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */ - if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */ + if(ijp <=cptcovprod) { /* Product Vn*Vm and age*VN*Vm*/ + if(DummyV[Tvardk[ijp][1]]==0){/* Vn is dummy */ + if(DummyV[Tvardk[ijp][2]]==0){/* Vn and Vm are dummy */ /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */ fprintf(ficgp,"+p%d*%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); }else{ /* Vn is dummy and Vm is quanti */ /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */ - fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]); } - }else{ /* Vn*Vm Vn is quanti */ + }else{ /* age* Vn*Vm Vn is quanti HERE */ if(DummyV[Tvard[ijp][2]]==0){ - fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); + fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[ijp][2]],Tqinvresult[nres][Tvardk[ijp][1]]); }else{ /* Both quanti */ - fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]); } } ijp++; @@ -12060,22 +9849,22 @@ set ter svg size 640, 480\nunset log y\n /* if(j==Tprod[ijp]) { /\* *\/ */ /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */ if(ijp <=cptcovprod) { /* Product */ - if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */ - if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are dummy */ + if(DummyV[Tvardk[ijp][1]]==0){/* Vn is dummy */ + if(DummyV[Tvardk[ijp][2]]==0){/* Vn and Vm are dummy */ /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */ - fprintf(ficgp,"+p%d*%d*%d*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); + fprintf(ficgp,"+p%d*%d*%d*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][1]],Tinvresult[nres][Tvardk[ijp][2]]); /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); */ }else{ /* Vn is dummy and Vm is quanti */ /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */ - fprintf(ficgp,"+p%d*%d*%f*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + fprintf(ficgp,"+p%d*%d*%f*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]); /* fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */ } }else{ /* Vn*Vm Vn is quanti */ - if(DummyV[Tvard[ijp][2]]==0){ - fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); + if(DummyV[Tvardk[ijp][2]]==0){ + fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][2]],Tqinvresult[nres][Tvardk[ijp][1]]); /* fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); */ }else{ /* Both quanti */ - fprintf(ficgp,"+p%d*%f*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); + fprintf(ficgp,"+p%d*%f*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]); /* fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */ } } @@ -12437,30 +10226,36 @@ void prevforecast(char fileres[], double /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */ /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */ /* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */ - i1=pow(2,cptcoveff); - if (cptcovn < 1){i1=1;} + /* i1=pow(2,cptcoveff); */ + /* if (cptcovn < 1){i1=1;} */ fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); fprintf(ficresf,"#****** Routine prevforecast **\n"); /* if (h==(int)(YEARM*yearp)){ */ - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ /* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) */ - if(i1 != 1 && TKresult[nres]!= k) - continue; - if(invalidvarcomb[k]){ - printf("\nCombination (%d) projection ignored because no cases \n",k); - continue; - } + for(nres=1; nres <= nresult; nres++){ /* For each resultline */ + k=TKresult[nres]; + if(TKresult[nres]==0) k=1; /* To be checked for noresult */ + /* for(k=1; k<=i1;k++){ /\* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) *\/ */ + /* if(i1 != 1 && TKresult[nres]!= k) */ + /* continue; */ + /* if(invalidvarcomb[k]){ */ + /* printf("\nCombination (%d) projection ignored because no cases \n",k); */ + /* continue; */ + /* } */ fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); - for(j=1;j<=cptcoveff;j++) { - /* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); */ - fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - } - for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ - fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); + for(j=1;j<=cptcovs;j++){ + /* for(j=1;j<=cptcoveff;j++) { */ + /* /\* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); *\/ */ + /* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + /* } */ + /* for (k4=1; k4<= nsq; k4++){ /\* For each selected (single) quantitative value *\/ */ + /* fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */ + /* } */ + fprintf(ficresf," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); } + fprintf(ficresf," yearproj age"); for(j=1; j<=nlstate+ndeath;j++){ for(i=1; i<=nlstate;i++) @@ -12485,9 +10280,11 @@ void prevforecast(char fileres[], double } } fprintf(ficresf,"\n"); - for(j=1;j<=cptcoveff;j++) + /* for(j=1;j<=cptcoveff;j++) */ + for(j=1;j<=cptcovs;j++) + fprintf(ficresf,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); /\* Tvaraff not correct *\/ */ - fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /* TnsdVar[Tvaraff] correct */ + /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /\* TnsdVar[Tvaraff] correct *\/ */ fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm); for(j=1; j<=nlstate+ndeath;j++) { @@ -12579,29 +10376,35 @@ void prevforecast(char fileres[], double /* if(jintmean==0) jintmean=1; */ /* if(mintmean==0) jintmean=1; */ - i1=pow(2,cptcoveff); - if (cptcovn < 1){i1=1;} + /* i1=pow(2,cptcoveff); */ + /* if (cptcovn < 1){i1=1;} */ fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ - if(i1 != 1 && TKresult[nres]!= k) - continue; - if(invalidvarcomb[k]){ - printf("\nCombination (%d) projection ignored because no cases \n",k); - continue; - } + for(nres=1; nres <= nresult; nres++){ /* For each resultline */ + k=TKresult[nres]; + if(TKresult[nres]==0) k=1; /* To be checked for noresult */ + /* for(k=1; k<=i1;k++){ */ + /* if(i1 != 1 && TKresult[nres]!= k) */ + /* continue; */ + /* if(invalidvarcomb[k]){ */ + /* printf("\nCombination (%d) projection ignored because no cases \n",k); */ + /* continue; */ + /* } */ fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#"); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - } - for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ - fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); + for(j=1;j<=cptcovs;j++){ + /* for(j=1;j<=cptcoveff;j++) { */ + /* fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + /* } */ + fprintf(ficresfb," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); } + /* fprintf(ficrespij,"******\n"); */ + /* for (k4=1; k4<= nsq; k4++){ /\* For each selected (single) quantitative value *\/ */ + /* fprintf(ficresfb," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */ + /* } */ fprintf(ficresfb," yearbproj age"); for(j=1; j<=nlstate+ndeath;j++){ for(i=1; i<=nlstate;i++) @@ -12632,8 +10435,10 @@ void prevforecast(char fileres[], double } } fprintf(ficresfb,"\n"); - for(j=1;j<=cptcoveff;j++) - fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); + /* for(j=1;j<=cptcoveff;j++) */ + for(j=1;j<=cptcovs;j++) + fprintf(ficresfb,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); + /* fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm); for(i=1; i<=nlstate+ndeath;i++) { ppij=0.;ppi=0.; @@ -13597,14 +11402,15 @@ int decoderesult( char resultline[], int if (strlen(resultsav) >1){ j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */ } - if(j == 0){ /* Resultline but no = */ + if(j == 0 && cptcovs== 0){ /* Resultline but no = and no covariate in the model */ TKresult[nres]=0; /* Combination for the nresult and the model */ return (0); } if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ - printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model); - fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model); - /* return 1;*/ + fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(ficlog); + printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(stdout); + if(j==0) + return 1; } for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */ if(nbocc(resultsav,'=') >1){ @@ -13787,20 +11593,29 @@ int decoderesult( char resultline[], int precov[nres][k1]=Tvalsel[k3q]; /* printf("Decoderesult Quantitative nres=%d,precov[nres=%d][k1=%d]=%.f V(k2q=V%d)= Tvalsel[%d]=%d, Tvarsel[%d]=%f\n",nres, nres, k1,precov[nres][k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); */ k4q++;; - }else if( Dummy[k1]==2 ){ /* For dummy with age product */ - /* Tvar[k1]; */ /* Age variable */ + }else if( Dummy[k1]==2 ){ /* For dummy with age product "V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/ + /* Tvar[k1]; */ /* Age variable */ /* 17 age*V6*V2 ?*/ /* Wrong we want the value of variable name Tvar[k1] */ - - k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/ - k2=(int)Tvarsel[k3]; /* nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/ - TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][4]=1 */ - precov[nres][k1]=Tvalsel[k3]; + if(Typevar[k1]==2 || Typevar[k1]==3 ){ /* For product quant or dummy (with or without age) */ + precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]]; + /* printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]); */ + }else{ + k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/ + k2=(int)Tvarsel[k3]; /* nres=1 k1=2=>k3=1 Tvarsel[resultmodel[2]]= Tvarsel[1] = 4=k2 (V4); k1=3=>k3=2 Tvarsel[2]=3 (V3)*/ + TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][4]=1 */ + precov[nres][k1]=Tvalsel[k3]; + } /* printf("Decoderesult Dummy with age k=%d, k1=%d precov[nres=%d][k1=%d]=%.f Tvar[%d]=V%d k2=Tvarsel[%d]=%d Tvalsel[%d]=%d\n",k, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k3,(int)Tvarsel[k3], k3, (int)Tvalsel[k3]); */ }else if( Dummy[k1]==3 ){ /* For quant with age product */ - k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */ - k2q=(int)Tvarsel[k3q]; /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */ - TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */ - precov[nres][k1]=Tvalsel[k3q]; + if(Typevar[k1]==2 || Typevar[k1]==3 ){ /* For product quant or dummy (with or without age) */ + precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]]; + /* printf("Decoderesult Quantitative or Dummy (not with age) nres=%d k1=%d precov[nres=%d][k1=%d]=%.f V%d(=%.f) * V%d(=%.f) \n",nres, k1, nres, k1,precov[nres][k1], Tvardk[k1][1], TinvDoQresult[nres][Tvardk[k1][1]], Tvardk[k1][2], TinvDoQresult[nres][Tvardk[k1][2]]); */ + }else{ + k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */ + k2q=(int)Tvarsel[k3q]; /* Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */ + TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */ + precov[nres][k1]=Tvalsel[k3q]; + } /* printf("Decoderesult Quantitative with age nres=%d, k1=%d, precov[nres=%d][k1=%d]=%f Tvar[%d]=V%d V(k2q=%d)= Tvarsel[%d]=%d, Tvalsel[%d]=%f\n",nres, k1, nres, k1,precov[nres][k1], k1, Tvar[k1], k2q, k3q, Tvarsel[k3q], k3q, Tvalsel[k3q]); */ }else if(Typevar[k1]==2 || Typevar[k1]==3 ){ /* For product quant or dummy (with or without age) */ precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]]; @@ -13883,17 +11698,17 @@ int decodemodel( char model[], int lasto if (strlen(modelsav) >1){ /* V2 +V3 +V4 +V6 +V7 +V6*V2 +V7*V2 +V6*V3 +V7*V3 +V6*V4 +V7*V4 +age*V2 +age*V3 +age*V4 +age*V6 +age*V7 +age*V6*V2 +V7*V2 +age*V6*V3 +age*V7*V3 +age*V6*V4 +age*V7*V4 */ j=nbocc(modelsav,'+'); /**< j=Number of '+' */ j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */ - /* cptcovs=j+1-j1; */ /* is wrong , see after */ + cptcovs=0; /**< Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age => V1 + V3 =4+1-3=2 Wrong */ cptcovt= j+1; /* Number of total covariates in the model, not including * cst, age and age*age * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/ /* including age products which are counted in cptcovage. * but the covariates which are products must be treated * separately: ncovn=4- 2=2 (V1+V3). */ - cptcovprod=0; /**< Number of products and single product with age V1*V2 +v3*age = 2 */ - cptcovdageprod=0; /* Number of double products with age age*Vn*VM or Vn*age*Vm or Vn*Vm*age */ + cptcovprod=0; /**< Number of products V1*V2 +v3*age = 2 */ + cptcovdageprod=0; /* Number of doouble products with age age*Vn*VM or Vn*age*Vm or Vn*Vm*age */ cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */ - cptcovprodage=0; /**< Number of varying covariate products with age: age*V6(v)*V3(f) =1 */ + cptcovprodage=0; /* cptcovprodage=nboccstr(modelsav,"age");*/ /* Design @@ -13948,6 +11763,22 @@ int decodemodel( char model[], int lasto Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0; } cptcovage=0; + + /* First loop in order to calculate */ + /* for age*VN*Vm + * Provides, Typevar[k], Tage[cptcovage], existcomb[n][m], FixedV[ncovcolt+k12] + * Tprod[k1]=k Tposprod[k]=k1; Tvard[k1][1] =m; + */ + /* Needs FixedV[Tvardk[k][1]] */ + /* For others: + * Sets Typevar[k]; + * Tvar[k]=ncovcol+nqv+ntv+nqtv+k11; + * Tposprod[k]=k11; + * Tprod[k11]=k; + * Tvardk[k][1] =m; + * Needs FixedV[Tvardk[k][1]] == 0 + */ + for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model line */ cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' cutl from left to right modelsav==V2+V1+V5*age+V4+V3*age strb=V3*age stra=V2+V1V5*age+V4 */ /* "V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */ @@ -13956,26 +11787,26 @@ int decodemodel( char model[], int lasto /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ /*scanf("%d",i);*/ if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V5*age+ V4+V3*age strb=V3*age OR double product with age strb=age*V6*V2 or V6*V2*age or V6*age*V2 */ - cutl(strc,strd,strb,'*'); /**< k=1 strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 OR strb=age*V6*V2 strc=V6*V2 strd=age OR strb=V6*age*V2 c=age*V2 d=V6 OR b=V6*V2*age c=V2*age d=V6 */ - if(strchr(strc,'*')) { /**< Model with age and DOUBLE product: allowed since 0.99r44, strc=V6*V2 or V2*age or age*V2, strd=age or V6 or V6 OR (strb=age*V6*V2 or V6*V2*age or V6*age*V2) */ + cutl(strc,strd,strb,'*'); /**< k=1 strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 OR strb=age*V6*V2 strc=V6*V2 strd=age OR c=V2*age OR c=age*V2 */ + if(strchr(strc,'*')) { /**< Model with age and DOUBLE product: allowed since 0.99r44, strc=V6*V2 or V2*age or age*V2, strd=age or V6 or V6 */ Typevar[k]=3; /* 3 for age and double product age*Vn*Vm varying of fixed */ - if(strstr(strc,"age")!=0) { /* It means that strc=V2*age=Vm*age or age*V2 (not V6*V2) and thus that strd=Vn and strb=V6*V2*age or V6*age*V2 (but not age*V6*V2) */ - cutl(stre,strf,strc,'*') ; /* if strc=age*Vm then stre=Vm and strf=age, if strc=Vm*age then stre=age and strf=Vm. */ - strcpy(strc,strb); /* save strb(=age*Vn*Vm) into strc , strd=Vn */ + if(strstr(strc,"age")!=0) { /* It means that strc=V2*age or age*V2 and thus that strd=Vn */ + cutl(stre,strf,strc,'*') ; /* strf=age or Vm, stre=Vm or age. If strc=V6*V2 then strf=V6 and stre=V2 */ + strcpy(strc,strb); /* save strb(=age*Vn*Vm) into strc */ /* We want strb=Vn*Vm */ - if(strcmp(strf,"age")==0){ /* strf is "age" so stre=Vm =V2 (strc=age*Vm and strb Vn*age*Vm) . */ - strcpy(strb,strd); /* strd=Vn */ - strcat(strb,"*"); - strcat(strb,stre);/* strb=Vn*Vm */ - }else{ /* strf=Vm so stre=age. strd=Vn If strf=V6 then stre=V2 */ + if(strcmp(strf,"age")==0){ /* strf is "age" so that stre=Vm =V2 . */ + strcpy(strb,strd); + strcat(strb,"*"); + strcat(strb,stre); + }else{ /* strf=Vm If strf=V6 then stre=V2 */ strcpy(strb,strf); strcat(strb,"*"); - strcat(strb,strd); /* strb=Vm*Vn */ - strcpy(strd,strb); /* in order for strd to not be "age" for next test (will be strd=Vn*Vm */ + strcat(strb,stre); + strcpy(strd,strb); /* in order for strd to not be "age" for next test (will be Vn*Vm */ } - printf("DEBUG FIXED k=%d, Tage[k]=%d, Tvar[Tage[k]=%d,FixedV[Tvar[Tage[k]]]=%d\n",k,Tage[k],Tvar[Tage[k]],FixedV[Tvar[Tage[k]]]); - /* FixedV[Tvar[Tage[k]]]=0;*/ /* HERY not sure */ - }else{ /* strb=age*Vn*Vm strc=Vn*Vm (and strd=age) and should be strb=Vn*Vm but want to keep original strb double product */ + /* printf("DEBUG FIXED k=%d, Tage[k]=%d, Tvar[Tage[k]=%d,FixedV[Tvar[Tage[k]]]=%d\n",k,Tage[k],Tvar[Tage[k]],FixedV[Tvar[Tage[k]]]); */ + /* FixedV[Tvar[Tage[k]]]=0; /\* HERY not sure if V7*V4*age Fixed might not exist yet*\/ */ + }else{ /* strc=Vn*Vm (and strd=age) and should be strb=Vn*Vm but want to keep original strb double product */ strcpy(stre,strb); /* save full b in stre */ strcpy(strb,strc); /* save short c in new short b for next block strb=Vn*Vm*/ strcpy(strf,strc); /* save short c in new short f */ @@ -14007,7 +11838,7 @@ int decodemodel( char model[], int lasto Tvardk[k][1] =m; /* m 1 for V1*/ Tvard[k1][2] =n; /* n 4 for V4*/ Tvardk[k][2] =n; /* n 4 for V4*/ -/* Tvar[Tage[cptcovage]]=k1;*/ /* Tvar[6=age*V3*V2]=9 (new fixed covariate) */ +/* Tvar[Tage[cptcovage]]=k1;*/ /* Tvar[6=age*V3*V2]=9 (new fixed covariate) */ /* We don't know about Fixed yet HERE */ if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */ for (i=1; i<=lastobs;i++){/* For fixed product */ /* Computes the new covariate which is a product of @@ -14043,9 +11874,9 @@ int decodemodel( char model[], int lasto /*Tage[cptcovage]=k;*/ /* For age*V3*V2 Tage[1]=V3*V3=9 HERY too*/ /* Tvar[Tage[cptcovage]]=k1; */ cptcovprodvage++; - k12=2*k11-1; FixedV[ncovcolt+k12]=1; /* We expand Vn*Vm */ - FixedV[ncovcolt+k12+1]=1; + k12++; + FixedV[ncovcolt+k12]=1; } } /* Tage[cptcovage]=k; /\* V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 *\/ */ @@ -14157,6 +11988,8 @@ int decodemodel( char model[], int lasto /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); scanf("%d",i);*/ } /* end of loop + on total covariates */ + + } /* end if strlen(modelsave == 0) age*age might exist */ } /* end if strlen(model == 0) */ cptcovs=cptcovt - cptcovdageprod - cptcovprod;/**< Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age + age*v4*V3=> V1 + V3 =4+1-3=2 */ @@ -14195,6 +12028,11 @@ Fixed[k] 0=fixed (product or simple), 1 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); for(k=-1;k<=NCOVMAX; k++){ Fixed[k]=0; Dummy[k]=0;} for(k=1;k<=NCOVMAX; k++){TvarFind[k]=0; TvarVind[k]=0;} + + + /* Second loop for calculating Fixed[k], Dummy[k]*/ + + for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0,ncovva=0,ncovvta=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt loop on k from model */ if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ Fixed[k]= 0; @@ -15448,7 +13286,7 @@ int main(int argc, char *argv[]) /* double ***mobaverage; */ double wald; - char line[MAXLINE]; + char line[MAXLINE], linetmp[MAXLINE]; char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE]; char modeltemp[MAXLINE]; @@ -15550,7 +13388,7 @@ int main(int argc, char *argv[]) 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)){ @@ -15781,7 +13619,18 @@ int main(int argc, char *argv[]) }else break; } - if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ + if((num_filled=sscanf(line,"model=%[^.\n]", model)) !=EOF){ /* Every character after model but dot and return */ + if (num_filled != 1){ + printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); + fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); + model[0]='\0'; + goto end; + }else{ + trimbtab(linetmp,line); /* Trims multiple blanks in line */ + strcpy(line, linetmp); + } + } + if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ /* Every character after 1+age but dot and return */ if (num_filled != 1){ printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); @@ -16154,7 +14003,7 @@ Please run with mle=-1 to get a correct Tvard=imatrix(1,NCOVMAX,1,2); /* n=Tvard[k1][1] and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm * For V3*V2 (in V2+V1+V1*V4+age*V3+V3*V2), V3*V2 position is 2nd. * Tvard[k1=2][1]=3 (V3) Tvard[k1=2][2]=2(V2) */ - Tvardk=imatrix(-1,NCOVMAX,1,2); + Tvardk=imatrix(0,NCOVMAX,1,2); Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age 4 covariates (3 plus signs) Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 @@ -17259,7 +15108,7 @@ Please run with mle=-1 to get a correct date2dmy(datebackf,&jbackf, &mbackf, &anbackf); } - printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, prevbcast, pathc,p, (int)anprojd-bage, (int)anbackd-fage); + printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, prevbcast, pathc,p, (int)anprojd-bage, (int)anbackd-fage);/* HERE valgrind Tvard*/ } printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \ model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,prevbcast, estepm, \ @@ -17409,18 +15258,21 @@ Please run with mle=-1 to get a correct pstamp(ficreseij); - i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */ - if (cptcovn < 1){i1=1;} + /* i1=pow(2,cptcoveff); /\* Number of combination of dummy covariates *\/ */ + /* if (cptcovn < 1){i1=1;} */ - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ - if(i1 != 1 && TKresult[nres]!= k) - continue; + for(nres=1; nres <= nresult; nres++){ /* For each resultline */ + /* for(k=1; k<=i1;k++){ /\* For any combination of dummy covariates, fixed and varying *\/ */ + /* if(i1 != 1 && TKresult[nres]!= k) */ + /* continue; */ fprintf(ficreseij,"\n#****** "); printf("\n#****** "); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); + for(j=1;j<=cptcovs;j++){ + /* for(j=1;j<=cptcoveff;j++) { */ + /* fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + fprintf(ficreseij," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); + printf(" V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); + /* printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ printf(" V%d=%lg ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); /* TvarsQ[j] gives the name of the jth quantitative (fixed or time v) */ @@ -17489,7 +15341,7 @@ Please run with mle=-1 to get a correct /* */ 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); + 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); /* It might not be a good idea to mix dummies and quantitative */ @@ -17664,7 +15516,7 @@ Please run with mle=-1 to get a correct free_vector(weight,firstobs,lastobs); - free_imatrix(Tvardk,-1,NCOVMAX,1,2); + free_imatrix(Tvardk,0,NCOVMAX,1,2); free_imatrix(Tvard,1,NCOVMAX,1,2); free_imatrix(s,1,maxwav+1,firstobs,lastobs); free_matrix(anint,1,maxwav,firstobs,lastobs);