/* $Id$
$State$
$Log$
+ 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
/* $State$ */
#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 copyright[]="April 2023,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2022";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
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];
/* 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 <repeats 11 times>} */
+/* 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 */
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' *\/ */
}
/**** praxis ****/
-# include <float.h>
-/* # include <math.h> */
-/* # include <stdio.h> */
-/* # include <stdlib.h> */
-/* # include <string.h> */
-/* # include <time.h> */
-
-# 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];
- }
+ /* 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;
}
}
-
- 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;
- }
+ /* 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);
+}
- 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 );
-
- return;
-# undef TIME_SIZE
-}
/* end praxis */
/*************** powell ************************/
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);
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]);
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 */
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<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; /* Should be changed here */
- for (kk=1; kk<=cptcovprodage;kk++) {/* + age*V3*V2 +age*V2 +age*V3 +age*V4 For age product with simple covariates or product of fixed covariates */
- /* if(!FixedV[Tvar[Tage[kk]]]) */
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
- /* else*/
- /*cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; *//* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
- }
- for(ncovva=1, iposold=0; ncovva <= ncovvta ; ncovva++){ /* Time varying covariates with age including individual from products, product is computed dynamically */
- itv=TvarVVA[ncovva]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm */
- ipos=TvarVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
- if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
- cotvarv=cotvar[mw[mi][i]][TvarVVA[ncovva]][i]; /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+ /* 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 */
}
cotvarv=cotvarv*cotvarvold;
}
iposold=ipos;
- cov[ioffset+ipos]=cotvarv*agexact;
- /* For products */
+ cov[ioffset+ipos]=cotvarv;
}
+ /* 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<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; /* Should be changed here */
+ /* for (kk=1; kk<=cptcovage;kk++) { */
+ /* if(!FixedV[Tvar[Tage[kk]]]) */
+ /* cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /\* Tage[kk] gives the data-covariate associated with age *\/ */
+ /* else */
+ /* cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/ */
+ /* } */
+ for(ncovva=1, iposold=0; ncovva <= ncovta ; ncovva++){ /* Time varying covariates with age including individual from products, product is computed dynamically */
+ itv=TvarAVVA[ncovva]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm */
+ ipos=TvarAVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+ if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
+ cotvarv=cotvar[mw[mi][i]][TvarAVVA[ncovva]][i]; /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+ }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*agexact;
+ /* For products */
+ }
+
+ 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 */
- 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 */
-
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */
- /* 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.
- */
+ /* 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 */
/* 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; 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++) {
+ if(!FixedV[Tvar[Tage[kk]]])
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+ else
+ cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
+ }
+ 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]>1.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<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;
+ }
- } 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 if (mle==1){
- 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 */
- } 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<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++) {
+ if(!FixedV[Tvar[Tage[kk]]])
+ cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */
+ else
+ cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
+ }
+
+ 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];
+ 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("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",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/
+ } /* end of wave */
+ } /* end of individual */
+ } /* End of if */
for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */
l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
int ioffset=0;
int ipos=0,iposold=0,ncovv=0;
- char labficresilk[NCOVMAX+1];
double cotvarv, cotvarvold;
double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
double **out;
/* Fixed */
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
/* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products *\/ */
- /*strcpy(labficresilk," "); */ /* A string fed with values of covariates and print in ficresilk at the end of the line to produce graphics with covariate values and also in order to verify the process of calculating the likelihood */
for (kf=1; kf<=ncovf;kf++){ /* V2 + V3 + V4 Simple and product fixed covariates without age* products *//* Missing values are set to -1 but should be dropped */
/* printf("Debug3 TvarFind[%d]=%d",kf, TvarFind[kf]); */
/* printf(" Tvar[TvarFind[kf]]=%d", Tvar[TvarFind[kf]]); */
/* printf(" i=%d covar[Tvar[TvarFind[kf]]][i]=%f\n",i,covar[Tvar[TvarFind[kf]]][i]); */
- ipos=TvarFind[kf];
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 (k=6)*/
- /* if(globpr) */
- /* sprintf(labficresilk+strlen(labficresilk)," %g",cov[ioffset+ipos]); */
/* cov[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i]; */
/* cov[2+6]=covar[Tvar[6]][i]; */
/* cov[2+6]=covar[2][i]; V2 */
/* cov[2+9]=covar[Tvar[9]][i]; */
/* cov[2+9]=covar[1][i]; V1 */
}
- /* 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 in the DATA is the cotvar[mw[mi+1][i]][ncovcol+nqv+iv][i]
- */
+ /* 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 in the DATA is the cotvar[mw[mi+1][i]][ncovcol+nqv+iv][i]
+ */
/* This part may be useless now because everythin should be in covar */
/* for (k=1; k<=nqfveff;k++){ /\* Simple and product fixed Quantitative covariates without age* products *\/ */
/* cov[++ioffset]=coqvar[TvarFQ[k]][i];/\* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V2 and V1*V2 is fixed (k=6 and 7?)*\/ */
* Not a fixed cotvar[mw][itv][i] 6 7 6 2 7, 2, 6, 3, 7, 3, 6, 4, 7, 4}
* fixed covar[itv] [6] [7] [6][2]
*/
-
+
for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* V6 V7 V7*V2 V6*V3 V7*V3 V6*V4 V7*V4 Time varying covariates (single and extended product but no age) including individual from products, product is computed dynamically */
itv=TvarVV[ncovv]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, or fixed covariate of a varying product after exploding product Vn*Vm into Vn and then Vm */
ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
cotvarv=covar[itv][i]; /* Good: In V6*V3, 3 is fixed at position of the data */
/* printf("DEBUG Fixed cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */
}
- /* if(globpr) */
- /* sprintf(labficresilk+strlen(labficresilk)," %g",cotvarv); */
if(ipos!=iposold){ /* Not a product or first of a product */
cotvarvold=cotvarv;
}else{ /* A second product */
cotvarv=cotvarv*cotvarvold;
- /* if(globpr) */
- /* sprintf(labficresilk+strlen(labficresilk)," *"); */
- /* printf("DEBUG * \n"); */
}
iposold=ipos;
cov[ioffset+ipos]=cotvarv;
cov[2]=agexact;
if(nagesqr==1)
cov[3]= agexact*agexact;
- /* for (kk=1; kk<=cptcovage;kk++) { /\* + age*V3*V2 +age*V2 +age*V3 +age*V4 For age product with simple covariates or product of fixed covariates *\/ */
- for (kk=1; kk<=cptcovprodage;kk++) { /* + age*V3*V2 +age*V2 +age*V3 +age*V4 For age product with simple covariates or product of fixed covariates */
- ipos=Tage[kk];
- if(!FixedV[Tvar[Tage[kk]]]){ /* age*V3*V2 +age*V2 +age*V3 +age*V4 Fixed covariate with age age*V3 or age*V2*V3 Tvar[Tage[kk]] has its own already calculated column */
- /* printf("DEBUG kk=%d, Fixed Tvar[Tage[kk]]=%d agexact=%g\n",kk, Tvar[Tage[kk]], agexact); */
- cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
- /* printf("DEBUG Fixed cov[Tage[kk]+2+nagesqr=%d]=%g agexact=%g \n",Tage[kk]+2+nagesqr,cov[Tage[kk]+2+nagesqr], agexact); */
- }else{ /* +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 time varying covariates with age age*V7 or age*V7*V3 is Tvar[Tage[kk]] defined, yes*/
- /* printf("DEBUG kk=%d, Varyingd Tvar[Tage[kk]]=%d\n",kk, Tvar[Tage[kk]]); */
- cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]*agexact; /* Are you sure? because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
- /* printf("DEBUG Varying cov[Tage[kk]+2+nagesqr=%d]=%g agexact=%g \n",Tage[kk]+2+nagesqr,cov[Tage[kk]+2+nagesqr], agexact); */
- }
- /* if(globpr) */
- /* sprintf(labficresilk+strlen(labficresilk)," %g*ageF",cov[Tage[kk]+2+nagesqr]); */
- }
- /* For product with age age*Vn*Vm where Vn*Vm is time varying */
- /* for(kv=1; kv<=cptcovprodvage;kv++){ /\*HERY? +age*V6 + age*V7 +age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 Number of time varying covariates with age *\/ */
- for(ncovva=1, iposold=0; ncovva <= ncovvta ; ncovva++){ /* +age*V6 + age*V7 +age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 Time varying covariates with age including individual from products, product is computed dynamically */
- itv=TvarVVA[ncovva]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm */
- ipos=TvarVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+ for(ncovva=1, iposold=0; ncovva <= ncovta ; ncovva++){ /* Time varying covariates with age including individual from products, product is computed dynamically */
+ itv=TvarAVVA[ncovva]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm */
+ ipos=TvarAVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
/* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
- /* printf("DEBUG ncovva=%d, Varying TvarVVA[ncovva]=%d agexact=%g\n", ncovva, TvarVVA[ncovva], agexact); */
- cotvarv=cotvar[mw[mi][i]][TvarVVA[ncovva]][i]; /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+ /* printf("DEBUG ncovva=%d, Varying TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
+ cotvarv=cotvar[mw[mi][i]][TvarAVVA[ncovva]][i]; /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
}else{ /* fixed covariate */
/* cotvarv=covar[Tvar[TvarFind[itv]]][i]; /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
- /* printf("DEBUG ncovva=%d, Fixed TvarVVA[ncovva]=%d agexact=%g\n", ncovva, TvarVVA[ncovva], agexact); */
+ /* printf("DEBUG ncovva=%d, Fixed TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
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 */
cotvarv=cotvarv*cotvarvold;
}
iposold=ipos;
+ /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
cov[ioffset+ipos]=cotvarv*agexact;
- /* printf("DEBUG Product cov[ioffset+ipos=%d]=%g, agexact=%g.3 \n",ioffset+ipos,cov[ioffset+ipos], agexact); */
- /* if(globpr) */
- /* sprintf(labficresilk+strlen(labficresilk)," %g*age",cov[ioffset+ipos]); */
-
/* For products */
}
/* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */
- /* for(kk=1;kk<=ncovmodel;kk++){ */
- /* printf(" %d=%11.6f",kk,cov[kk]); */
- /* } */
- /* printf("\n"); */
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
%11.6f %11.6f %11.6f ", \
num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2]));
-
- /* printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */
+ /* printf("%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\ */
/* %11.6f %11.6f %11.6f ", \ */
/* num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw, */
/* 2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
/* printf(" %10.6f",-ll[k]*gipmx/gsw); */
}
fprintf(ficresilk," %10.6f ", -llt);
- /* fprintf(ficresilk,"%s", labficresilk); */
/* printf(" %10.6f\n", -llt); */
/* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */
/* fprintf(ficresilk,"%09ld ", num[i]); */ /* not necessary */
}
iposold=ipos;
}
- for(ncovva=1, iposold=0; ncovva <= ncovvta ; ncovva++){ /* Time varying covariates with age including individual from products, product is computed dynamically */
- itv=TvarVVA[ncovva]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm */
- ipos=TvarVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
- fprintf(ficresilk," %g*age",cov[ioffset+ipos]);
- }
- /* if(FixedV[itv]!=0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
- /* cotvarv=cotvar[mw[mi][i]][TvarVVA[ncovva]][i]; /\* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) *\/ */
- /* }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 products *\/ */
- /* fprintf(ficresilk," %g*age",cotvarv);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/ */
- /* } */
/* for (kk=1; kk<=cptcovage;kk++) { */
/* if(!FixedV[Tvar[Tage[kk]]]){ */
/* fprintf(ficresilk," %g*age",covar[Tvar[Tage[kk]]][i]); */
/* /\* printf(" %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/\\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\\/ *\/ */
/* } */
/* } */
+ for(ncovva=1, iposold=0; ncovva <= ncovta ; ncovva++){ /* Time varying covariates with age including individual from products, product is computed dynamically */
+ itv=TvarAVVA[ncovva]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm */
+ ipos=TvarAVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
+ /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
+ if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
+ /* printf("DEBUG ncovva=%d, Varying TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
+ cotvarv=cotvar[mw[mi][i]][TvarAVVA[ncovva]][i]; /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */
+ }else{ /* fixed covariate */
+ /* cotvarv=covar[Tvar[TvarFind[itv]]][i]; /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
+ /* printf("DEBUG ncovva=%d, Fixed TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
+ 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 */
+ /* printf("DEBUG * \n"); */
+ cotvarv=cotvarv*cotvarvold;
+ }
+ cotvarv=cotvarv*agexact;
+ fprintf(ficresilk," %g*age",cotvarv);
+ iposold=ipos;
+ /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
+ cov[ioffset+ipos]=cotvarv;
+ /* For products */
+ }
/* printf("\n"); */
/* } /\* End debugILK *\/ */
fprintf(ficresilk,"\n");
- /* printf("\n"); */
} /* End if globpr */
} /* end of wave */
} /* end of individual */
Plotting could be done.
*/
void pstamp(FILE *ficres);
- int k, kf, kk, kvar, kvarold, ncovv, iposold, ipos, itv;
+ int k, kf, kk, kvar, ncovv, iposold, ipos;
if(*globpri !=0){ /* Just counts and sums, no printings */
strcpy(fileresilk,"ILK_");
fprintf(ficresilk," -2*gipw/gsw*weight*ll(total) ");
/* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */
- for(kf=1;kf <= ncovf; kf++){ /* Fixed covariates */
- fprintf(ficresilk,"V%d",Tvar[TvarFind[kf]]);
- /* printf("V%d",Tvar[TvarFind[kf]]); */
- }
- for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){ /* Time varying covariates V6 V7 V7*V2 V6*V3 V7*V3 V6*V4 V7*V4 Time varying covariates (single and extended product but no age) including individual from products, product is computed dynamically */
- itv=TvarVV[ncovv]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, or fixed covariate of a varying product after exploding product Vn*Vm into Vn and then Vm */
- 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? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
- kvar=TvarVV[ncovv];
- }else{
- kvar=itv;
- }
- if(ipos!=iposold){ /* Not a product or first of a product */
- kvarold=kvar;
- /* printf(" %d",ipos); */
- fprintf(ficresilk," V%d",kvarold);
- }else{ /* a second product */
- /* printf("*"); */
- fprintf(ficresilk,"*");
- fprintf(ficresilk," V%d",kvar);
- }
- iposold=ipos;
- }
- for (kk=1; kk<=cptcovprodage;kk++) { /* Fixed Covariates with age + age*V3*V2 +age*V2 +age*V3 +age*V4 For age product with simple covariates or product of fixed covariates */
- ipos=Tage[kk];
- if(!FixedV[Tvar[Tage[kk]]]){/* age*V3*V2 +age*V2 +age*V3 +age*V4 Fixed covariate with age age*V3 or age*V2*V3 Tvar[Tage[kk]] has its own already calculated column */
- /* printf(" %d*age(Fixed)",Tvar[Tage[kk]]); */
- fprintf(ficresilk," %d*age(Fixed)",Tvar[Tage[kk]]);
- }else{
- fprintf(ficresilk," %d*age(Varying)",Tvar[Tage[kk]]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
- /* printf(" %d*age(Varying)",Tvar[Tage[kk]]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/ */
+ for(kf=1;kf <= ncovf; kf++){
+ fprintf(ficresilk,"V%d",Tvar[TvarFind[kf]]);
+ /* printf("V%d",Tvar[TvarFind[kf]]); */
}
- }
- for(ncovva=1, iposold=0; ncovva <= ncovvta ; ncovva++){ /* +age*V6 + age*V7 +age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 Time varying covariates with age including individual from products, product is computed dynamically */
- itv=TvarVVA[ncovva]; /* TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm */
- ipos=TvarVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
- /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */
- if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */
- /* printf("DEBUG ncovva=%d, Varying TvarVVA[ncovva]=%d agexact=%g\n", ncovva, TvarVVA[ncovva], agexact); */
- kvar=TvarVVA[ncovva];
- /* cotvarv=cotvar[mw[mi][i]][TvarVVA[ncovva]][i]; /\* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) *\/ */
- }else{ /* fixed covariate */
- /* cotvarv=covar[Tvar[TvarFind[itv]]][i]; /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
- /* printf("DEBUG ncovva=%d, Fixed TvarVVA[ncovva]=%d agexact=%g\n", ncovva, TvarVVA[ncovva], agexact); */
- kvar=itv;
- /* cotvarv=covar[itv][i]; /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */
+ for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){
+ ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
+ if(ipos!=iposold){ /* Not a product or first of a product */
+ /* printf(" %d",ipos); */
+ fprintf(ficresilk," V%d",TvarVV[ncovv]);
+ }else{
+ /* printf("*"); */
+ fprintf(ficresilk,"*");
+ }
+ iposold=ipos;
}
- if(ipos!=iposold){ /* Not a product or first of a product */
- kvarold=kvar;
- fprintf(ficresilk," age*V%d",kvarold);
- }else{ /* a second product */
- /* printf("*"); */
- fprintf(ficresilk," *V%d",kvar);
- /* printf("DEBUG * \n"); */
+ for (kk=1; kk<=cptcovage;kk++) {
+ if(!FixedV[Tvar[Tage[kk]]]){
+ /* printf(" %d*age(Fixed)",Tvar[Tage[kk]]); */
+ fprintf(ficresilk," %d*age(Fixed)",Tvar[Tage[kk]]);
+ }else{
+ fprintf(ficresilk," %d*age(Varying)",Tvar[Tage[kk]]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */
+ /* printf(" %d*age(Varying)",Tvar[Tage[kk]]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/ */
+ }
}
- iposold=ipos;
- /* printf("DEBUG Product cov[ioffset+ipos=%d]=%g, agexact=%g.3 \n",ioffset+ipos,cov[ioffset+ipos], agexact); */
- /* if(globpr) */
- /* sprintf(labficresilk+strlen(labficresilk)," %g*age",cov[ioffset+ipos]); */
-
- /* For products */
- }
/* } /\* End if debugILK *\/ */
/* printf("\n"); */
fprintf(ficresilk,"\n");
fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br>\n \
<img src=\"%s-p%dj.png\">\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,"<br>- Probability p<sub>%dj</sub> 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): <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
-<img src=\"%s-p%dj-%d.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,"<br>- Probability p<sub>%dj</sub> 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,"<a href=\"%s-p%dj-%d.png\">%s-p%dj-%d.png</a><br>",subdirf2(optionfilefiname,"ILK_"),k,kvar,subdirf2(optionfilefiname,"ILK_"),k,kvar);
+ fprintf(fichtm,"<img src=\"%s-p%dj-%d.png\">",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 */
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);
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));
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 <cptcoveff && cptcoveff>1)
+ if(k <cptcovs && cptcovs>1)
sprintf(gplotcondition+strlen(gplotcondition)," && ");
}
strcpy(gplotcondition+strlen(gplotcondition),")");
}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 */
/*# 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 */
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; */
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++;
/* 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]]); */
}
}
/* 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++)
}
}
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++) {
/* 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++)
}
}
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.;
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){
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]];
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
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 */ /* <model> "V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */
/* 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 */
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
/*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 *\/ */
/*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 */
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;
/* double ***mobaverage; */
double wald;
- char line[MAXLINE];
+ char line[MAXLINE], linetmp[MAXLINE];
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
char modeltemp[MAXLINE];
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)){
}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);
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
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, \
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) */
/* */
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 */
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);
git://henry.ined.fr / .git / commitdiff