[BACK]Return to imachprax.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imachprax.c between versions 1.1 and 1.2

version 1.1, 2023/01/31 09:24:19 version 1.2, 2023/06/22 11:22:40
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
   Revision 1.1  2023/01/31 09:24:19  brouard    Revision 1.2  2023/06/22 11:22:40  brouard
   Summary: version s1 with praxis instead of Powell for large models with age and difficulties to converge    Summary: with svd but not working yet
   
     Revision 1.353  2023/05/08 18:48:22  brouard
     *** empty log message ***
   
     Revision 1.352  2023/04/29 10:46:21  brouard
     *** empty log message ***
   
     Revision 1.351  2023/04/29 10:43:47  brouard
     Summary: 099r45
   
     Revision 1.350  2023/04/24 11:38:06  brouard
     *** empty log message ***
   
     Revision 1.349  2023/01/31 09:19:37  brouard
     Summary: Improvements in models with age*Vn*Vm
   
   Revision 1.347  2022/09/18 14:36:44  brouard    Revision 1.347  2022/09/18 14:36:44  brouard
   Summary: version 0.99r42    Summary: version 0.99r42
Line 1362  double gnuplotversion=GNUPLOTVERSION; Line 1377  double gnuplotversion=GNUPLOTVERSION;
 /* $State$ */  /* $State$ */
 #include "version.h"  #include "version.h"
 char version[]=__IMACH_VERSION__;  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 fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 1474  extern time_t time(); Line 1489  extern time_t time();
   
 struct tm start_time, end_time, curr_time, last_time, forecast_time;  struct tm start_time, end_time, curr_time, last_time, forecast_time;
 time_t  rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */  time_t  rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */
   time_t   rlast_btime; /* raw time */
 struct tm tm;  struct tm tm;
   
 char strcurr[80], strfor[80];  char strcurr[80], strfor[80];
Line 1597  int **nbcode, *Tvar; /**< model=V2 => Tv Line 1613  int **nbcode, *Tvar; /**< model=V2 => Tv
 /* Tprod[i]=k             1               2             */ /* Position in model of the ith prod without age */  /* 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 */  /* 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 */  /* 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 */  /* 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*/  /* 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 */  /* 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 */
Line 1795  char *trimbb(char *out, char *in) Line 1816  char *trimbb(char *out, char *in)
   return s;    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) */  /* char *substrchaine(char *out, char *in, char *chain) */
 /* { */  /* { */
 /*   /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */  /*   /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */
Line 2578  void linmin(double p[], double xi[], int Line 2613  void linmin(double p[], double xi[], int
 }   } 
   
 /**** praxis ****/  /**** praxis ****/
 # include <float.h>  double pythag( double x, double y )
 /* # 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 ) */  
   
 /******************************************************************************/  /******************************************************************************/
 /*  /*
   Purpose:    Purpose:
       R8_HYPOT returns the value of sqrt ( X^2 + Y^2 ).
     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.  
   
   Licensing:    Licensing:
   
     This code is distributed under the GNU LGPL license.      This code is distributed under the GNU LGPL license.
   
   Modified:    Modified:
       26 March 2012
     28 July 2016  
   
   Author:    Author:
       John Burkardt
     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:    Parameters:
       Input, double X, Y, the arguments.
     Input, int N, the number of variables.      Output, double R8_HYPOT, the value of sqrt ( X^2 + Y^2 ).
   
     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.  
 */  */
 {  {
   int i;    double a;
   double *t;    double b;
   double value;    double value;
   
   t = ( double * ) malloc ( n * sizeof ( double ) );    if ( fabs ( x ) < fabs ( y ) )  {
 /*      a = fabs ( y );
   The search is linear.      b = fabs ( x );
 */    }  else  {
   if ( 0 <= jsearch )      a = fabs ( x );
   {      b = fabs ( y );
     for ( i = 0; i < n; i++ )  
     {  
       t[i] = x[i] + l * v[i+jsearch*n];  
     }  
   }    }
 /*  /*
   The search is along a parabolic space curve.    A contains the larger value.
 */  */
   else    if ( a == 0.0 )  {
   {      value = 0.0;
     *qa =                  l * ( l - *qd1 ) /        ( *qd0 + *qd1 ) / *qd0;    }  else  {
     *qb = - ( l + *qd0 ) *     ( l - *qd1 ) / *qd1                   / *qd0;      value = a * sqrt ( 1.0 + ( b / a ) * ( b / a ) );
     *qc =   ( l + *qd0 ) * l                / *qd1 / ( *qd0 + *qd1 );  
   
     for ( i = 0; i < n; i++ )  
     {  
       t[i] = *qa * q0[i] + *qb * x[i] + *qc * q1[i];  
     }  
   }    }
 /*  
   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;    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  /* void svdcmp(double **a, int m, int n, double w[], double **v)  */
     singular value decomposition is desired.  On output, the  void svdcmp(double **a, int m, int n, double w[]) 
     SVD orthogonal matrix factor V.  
   
     Input/output, double Q[N], the singular values.  
 */  
 {  {
   double c;    /* Golub 1970 Algol60 */
   double *e;  /*   Computation of the singular values and complete orthogonal decom- */
   double eps;  /* position of a real rectangular matrix A, */
   double f;  /* A = U diag (q) V^T, U^T U = V^T V =I , */
   double g;  /*     where the arrays a[1:m, 1:n], u[1:m, 1 :n], v[1:n, 1:n], q[1:n] re- */
   double h;  /* present A, U, V, q respectively. The actual parameters corresponding */
   int i;  /* to a, u, v may all be identical unless withu=withv = true . In this */
   int ii;  /* case, the actual parameters corresponding to u and v must differ. */
   int j;  /* m >= n is assumed; */
   int jj;    
   int k;    /* Simplified (as in praxis) in order to output only V (replacing A), w is the diagonal matrix q  */
   int kt;  
   const int kt_max = 30;    
   int l;    double pythag(double a, double b); 
   int l2;    int flag,i,its,j,jj,k,l,nm; 
   double s;    double anorm,c,f,g,h,s,scale,x,y,z,*rv1; 
   int skip;    
   double temp;    rv1=vector(1,n);
   double x;    /* Householder's reduction to bidiagonal form; */
   double y;    g=scale=anorm=0.0; 
   double z;    for (i=1;i<=n;i++) { 
 /*      l=i+1; 
   Householder's reduction to bidiagonal form.      rv1[i]=scale*g; 
 */      g=s=scale=0.0; 
   if ( n == 1 )      if (i <= m) { 
   {        for (k=i;k<=m;k++) scale += fabs(a[k][i]); 
     q[0] = a[0+0*n];        if (scale) { 
     a[0+0*n] = 1.0;          for (k=i;k<=m;k++) { 
     return;            a[k][i] /= scale; 
   }            s += a[k][i]*a[k][i]; 
           } 
   e = ( double * ) malloc ( n * sizeof ( double ) );          f=a[i][i]; 
           g = -SIGN(sqrt(s),f); 
   eps = DBL_EPSILON;          h=f*g-s; 
   g = 0.0;          a[i][i]=f-g; 
   x = 0.0;          for (j=l;j<=n;j++) { 
             for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j]; 
   for ( i = 1; i <= n; i++ )            f=s/h; 
   {            for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; 
     e[i-1] = g;          } 
     l = i + 1;          for (k=i;k<=m;k++) a[k][i] *= scale; 
         } 
     s = 0.0;      } 
     for ( ii = i; ii <= n; ii++ )      w[i]=scale *g; /* p 411*/
     {      g=s=scale=0.0; 
       s = s + a[ii-1+(i-1)*n] * a[ii-1+(i-1)*n];      if (i <= m && i != n) { 
     }        for (k=l;k<=n;k++) scale += fabs(a[i][k]); 
         if (scale) { 
     g = 0.0;          for (k=l;k<=n;k++) { 
             a[i][k] /= scale; 
     if ( tol <= s )            s += a[i][k]*a[i][k]; 
     {          } 
       f = a[i-1+(i-1)*n];          f=a[i][l]; 
           g = -SIGN(sqrt(s),f); 
       g = sqrt ( s );          h=f*g-s; 
           a[i][l]=f-g; 
       if ( 0.0 <= f )          for (k=l;k<=n;k++) rv1[k]=a[i][k]/h; 
       {          for (j=l;j<=m;j++) { 
         g = - g;            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]; 
           } 
       h = f * g - s;          for (k=l;k<=n;k++) a[i][k] *= scale; 
       a[i-1+(i-1)*n] = f - g;        } 
       }
       for ( j = l; j <= n; j++ )      /* y : = abs(q[i]) +abs(e[i]); if y > x then x : = y */
       {      /* anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));  */
         f = 0.0;      anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i]))); 
         for ( ii = i; ii <= n; ii++ )    }
         {    /* Comment accumulation of right-hand transformations */
           f = f + a[ii-1+(i-1)*n] * a[ii-1+(j-1)*n];    for (i=n;i>=1;i--) { 
         }      if (i < n) { 
         f = f / h;        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 ( ii = i; ii <= n; ii++ )          /* for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;  */
         {          for (j=l;j<=n;j++) { 
           a[ii-1+(j-1)*n] = a[ii-1+(j-1)*n] + f * a[ii-1+(i-1)*n];            /* 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]; 
           } 
     q[i-1] = g;        } 
         /* for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;  */
     s = 0.0;        for (j=l;j<=n;j++) a[i][j]=a[j][i]=0.0; 
     for ( j = l; j <= n; j++ )      } 
     {      /* v[i][i]=1.0;  */
       s = s + a[i-1+(j-1)*n] * a[i-1+(j-1)*n];      a[i][i]=1.0; 
     }      g=rv1[i]; 
       l=i; 
     g = 0.0;    } 
    /* Comment accumulation of left-hand transformations; */
     if ( tol <= s )   /* for (i=IMIN(m,n);i>=1;i--) {  */
     {   for (i=FMIN(m,n);i>=1;i--) { 
       if ( i < n )      l=i+1; 
       {      g=w[i]; 
         f = a[i-1+i*n];      for (j=l;j<=n;j++) a[i][j]=0.0; 
       }      if (g) { 
         g=1.0/g; 
       g = sqrt ( s );        for (j=l;j<=n;j++) { 
           for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j]; 
       if ( 0.0 <= f )          f=(s/a[i][i])*g; 
       {          for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; 
         g = - g;        } 
       }        for (j=i;j<=m;j++) a[j][i] *= g; 
       } else for (j=i;j<=m;j++) a[j][i]=0.0; 
       h = f * g - s;      ++a[i][i]; 
     } 
       if ( i < n )  /* comment diagonalization of the bidiagonal form; */
       {   for (k=n;k>=1;k--) { 
         a[i-1+i*n] = f - g;      for (its=1;its<=30;its++) { /* Loop over singular values, and over allowed iterations. */
         for ( jj = l; jj <= n; jj++ )        flag=1; 
         {        for (l=k;l>=1;l--) { /* Test for splitting. */
           e[jj-1] = a[i-1+(jj-1)*n] / h;          nm=l-1;            /* Note that rv1[1] is always zero. */
         }          if ((double)(fabs(rv1[l])+anorm) == anorm) { 
             flag=0; 
         for ( j = l; j <= n; j++ )            break; 
         {          } 
           s = 0.0;          if ((double)(fabs(w[nm])+anorm) == anorm) break; 
           for ( jj = l; jj <= n; jj++ )        } 
           {        if (flag) { 
             s = s + a[j-1+(jj-1)*n] * a[i-1+(jj-1)*n];          c=0.0;    /* Cancellation of rv1[l], if l > 1. */
           }          s=1.0; 
           for ( jj = l; jj <= n; jj++ )          for (i=l;i<=k;i++) { 
           {            f=s*rv1[i]; 
             a[j-1+(jj-1)*n] = a[j-1+(jj-1)*n] + s * e[jj-1];            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; 
     y = fabs ( q[i-1] ) + fabs ( e[i-1] );            c=g*h; 
             s = -f*h; 
     x = fmax ( x, y );            for (j=1;j<=m;j++) { 
   }              y=a[j][nm]; 
 /*              z=a[j][i]; 
   Accumulation of right-hand transformations.              a[j][nm]=y*c+z*s; 
 */              a[j][i]=z*c-y*s; 
   a[n-1+(n-1)*n] = 1.0;            } 
   g = e[n-1];          } 
   l = n;        } 
         z=w[k]; 
   for ( i = n - 1; 1 <= i; i-- )        if (l == k) { /* Convergence */
   {          if (z < 0.0) { /* Singular value is made nonnegative. */   
     if ( g != 0.0 )            w[k] = -z;   
     {            /* for (j=1;j<=n;j++) v[j][k] = -v[j][k];  */
       h = a[i-1+i*n] * g;            for (j=1;j<=n;j++) a[j][k] = -a[j][k]; 
           } 
       for ( ii = l; ii <= n; ii++ )          break; 
       {        } 
         a[ii-1+(i-1)*n] = a[i-1+(ii-1)*n] / h;        if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations"); 
       }        x=w[l];  /* shift from bottom 2 x 2 minor; */
         nm=k-1; 
       for ( j = l; j <= n; j++ )        y=w[nm]; 
       {        g=rv1[nm]; 
         s = 0.0;        h=rv1[k]; 
         for ( jj = l; jj <= n; jj++ )        f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); 
         {        g=pythag(f,1.0); 
           s = s + a[i-1+(jj-1)*n] * a[jj-1+(j-1)*n];        /* 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 ( ii = l; ii <= n; ii++ )        for (j=l;j<=nm;j++) { 
         {          i=j+1; 
           a[ii-1+(j-1)*n] = a[ii-1+(j-1)*n] + s * a[ii-1+(i-1)*n];          g=rv1[i]; 
         }          y=w[i]; 
       }          h=s*g; 
     }          g=c*g; 
           /* z=dpythag(f,h);  */
     for ( jj = l; jj <= n; jj++ )          z=pythag(f,h); 
     {          rv1[j]=z; 
       a[i-1+(jj-1)*n] = 0.0;          c=f/z; 
     }          s=h/z; 
           f=x*c+g*s; 
     for ( ii = l; ii <= n; ii++ )          g = g*c-x*s; 
     {          h=y*s; 
       a[ii-1+(i-1)*n] = 0.0;          y *= c; 
     }  /* if withv then for j:= 1 step 1 until n do */
           for (jj=1;jj<=n;jj++) { 
     a[i-1+(i-1)*n] = 1.0;            /* x=v[jj][j];  */
             /* z=v[jj][i];  */
     g = e[i-1];            x=a[jj][j]; 
             z=a[jj][i]; 
     l = i;            /* v[jj][j]=x*c+z*s;  */
   }            /* v[jj][i]=z*c-x*s;  */
 /*            a[jj][j]=x*c+z*s; 
   Diagonalization of the bidiagonal form.            a[jj][i]=z*c-x*s; 
 */          } 
   eps = eps * x;          /* z=dpythag(f,h);  */
           z=pythag(f,h); 
   for ( k = n; 1 <= k; k-- )          w[j]=z; 
   {          if (z) { 
     kt = 0;            z=1.0/z; 
             c=f*z; 
     for ( ; ; )            s=h*z; 
     {          } 
       kt = kt + 1;          f=c*g+s*y; 
           x=c*y-s*g;
       if ( kt_max < kt )          /* if withu then for j:=1 step 1 until m do */
       {          /* for (jj=1;jj<=m;jj++) {  */
         e[k-1] = 0.0;          /*   y=a[jj][j];  */
         fprintf ( stderr, "\n" );          /*   z=a[jj][i];  */
         fprintf ( stderr, "MINFIT - Fatal error!\n" );          /*   a[jj][j]=y*c+z*s;  */
         fprintf ( stderr, "  The QR algorithm failed to converge.\n" );          /*   a[jj][i]=z*c-y*s;  */
         exit ( 1 );          /* }  */
       }        } 
         rv1[l]=0.0; 
       skip = 0;        rv1[k]=f; 
         w[k]=x; 
       for ( l2 = k; 1 <= l2; l2-- )      } 
       {    } 
         l = l2;    free_vector(rv1,1,n); 
   } 
         if ( fabs ( e[l-1] ) <= eps )  
         {  
           skip = 1;  
           break;  
         }  
   
         if ( 1 < l )  
         {  
           if ( fabs ( q[l-2] ) <= eps )  
           {  
             break;  
           }  
         }  
       }  
 /*  
   Cancellation of E(L) if 1 < L.  
 */  
       if ( ! skip )  
       {  
         c = 0.0;  
         s = 1.0;  
   
         for ( i = l; i <= k; i++ )  
         {  
           f = s * e[i-1];  
           e[i-1] = c * e[i-1];  
           if ( fabs ( f ) <= eps )  
           {  
             break;  
           }  
           g = q[i-1];  
 /*  
   q(i) = h = sqrt(g*g + f*f).  
 */  
           h = r8_hypot ( f, g );  
     
           q[i-1] = h;  
   
           if ( h == 0.0 )  
           {  
             g = 1.0;  
             h = 1.0;  
           }  
   
           c =   g / h;  
           s = - f / h;  
         }  
       }  
 /*  
   Test for convergence for this index K.  
 */  
       z = q[k-1];  
   
       if ( l == k )  
       {  
         if ( z < 0.0 )  
         {  
           q[k-1] = - z;  
           for ( i = 1; i <= n; i++ )  
           {  
             a[i-1+(k-1)*n] = - a[i-1+(k-1)*n];  
           }  
         }  
         break;  
       }  
 /*  
   Shift from bottom 2*2 minor.  
 */  
       x = q[l-1];  
       y = q[k-2];  
       g = e[k-2];  
       h = e[k-1];  
       f = ( ( y - z ) * ( y + z ) + ( g - h ) * ( g + h ) ) / ( 2.0 * h * y );  
   
       g = r8_hypot ( f, 1.0 );  
   
       if ( f < 0.0 )  
       {  
         temp = f - g;  
       }  
       else  
       {  
         temp = f + g;  
       }  
   
       f = ( ( x - z ) * ( x + z ) + h * ( y / temp - h ) ) / x;  
 /*  
   Next QR transformation.  
 */  
       c = 1.0;  
       s = 1.0;  
   
       for ( i = l + 1; i <= k; i++ )  
       {  
         g = e[i-1];  
         y = q[i-1];  
         h = s * g;  
         g = g * c;  
   
         z = r8_hypot ( f, h );  
   
         e[i-2] = z;  
   
         if ( z == 0.0 )  
         {  
           f = 1.0;  
           z = 1.0;  
         }  
   
         c = f / z;  
         s = h / z;  
         f =   x * c + g * s;  
         g = - x * s + g * c;  
         h = y * s;  
         y = y * c;  
   
         for ( j = 1; j <= n; j++ )  
         {  
           x = a[j-1+(i-2)*n];  
           z = a[j-1+(i-1)*n];  
           a[j-1+(i-2)*n] =   x * c + z * s;  
           a[j-1+(i-1)*n] = - x * s + z * c;  
         }  
   
         z = r8_hypot ( f, h );  
   
         q[i-2] = z;  
   
         if ( z == 0.0 )  
         {  
           f = 1.0;  
           z = 1.0;  
         }  
   
         c = f / z;  
         s = h / z;  
         f =   c * g + s * y;  
         x = - s * g + c * y;  
       }  
   
       e[l-1] = 0.0;  
       e[k-1] = f;  
       q[k-1] = x;  
     }  
   }  
   
   free ( e );  
   
   return;  
 }  
 /******************************************************************************/  
   
 void minny ( int n, int jsearch, int nits, double *d2, double *x1, double *f1,   
              int fk, double (*func) ( double []), double x[], double t, double h,   
   double v[], double q0[], double q1[], int *nl, int *nf, double dmin,   
   double ldt, double *fx, double *qa, double *qb, double *qc, double *qd0,   
   double *qd1 )  
 /* void minny ( int n, int jsearch, int nits, double *d2, double *x1, double *f1,  */  
 /*   int fk, double f ( double x[], int n ), double x[], double t, double h,  */  
 /*   double v[], double q0[], double q1[], int *nl, int *nf, double dmin,  */  
 /*   double ldt, double *fx, double *qa, double *qb, double *qc, double *qd0,  */  
 /*   double *qd1 ) */  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     MINNY minimizes a scalar function of N variables along a line.  
   
   Discussion:  
   
     MINNY minimizes F along the line from X in the direction V(*,JSEARCH)   
     or else using a quadratic search in the plane defined by Q0, Q1 and X.  
   
     If FK = true, then F1 is FLIN(X1).  Otherwise X1 and F1 are ignored  
     on entry unless final FX is greater than F1.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     03 August 2016  
   
   Author:  
   
     Original FORTRAN77 version by Richard Brent.  
     C version by John Burkardt.  
   
   Reference:  
   
     Richard Brent,  
     Algorithms for Minimization with Derivatives,  
     Prentice Hall, 1973,  
     Reprinted by Dover, 2002.  
   
   Parameters:  
   
     Input, int N, the number of variables.  
   
     Input, int JSEARCH, indicates the kind of search.  
     If J is a legal columnindex, linear search in the direction of V(*,JSEARCH).  
     Otherwise, the search is parabolic, based on X, Q0 and Q1.  
   
     Input, int NITS, the maximum number of times the interval   
     may be halved to retry the calculation.  
   
     Input/output, double *D2, is either zero, or an approximation to   
     the value of (1/2) times the second derivative of F.  
   
     Input/output, double *X1, on entry, an estimate of the   
     distance from X to the minimum along V(*,JSEARCH), or a curve.    
     On output, the distance between X and the minimizer that was found.  
   
     Input/output, double *F1, ?  
   
     Input, int FK; if FK is TRUE, then on input F1 contains   
     the value FLIN(X1).  
   
     Input, double F ( double X[], int N ), is the name of the function to   
     be minimized.  
   
     Input/output, double X[N], ?  
   
     Input, double T, ?  
   
     Input, double H, ?  
   
     Input, double V[N,N], a matrix whose columns are direction  
     vectors along which the function may be minimized.  
   
     ?, double Q0[N], ?  
   
     ?, double Q1[N], ?  
   
     Input/output, int *NL, the number of linear searches.  
   
     Input/output, int *NF, the number of function evaluations.  
   
     Input, double DMIN, an estimate for the smallest eigenvalue.  
   
     Input, double LDT, the length of the step.  
   
     Input/output, double *FX, the value of F(X,N).  
   
     Input/output, double *QA, *QB, *QC;  
   
     Input/output, double *QD0, *QD1, ?.  
 */  
 {  
   double d1;  
   int dz;  
   double f0;  
   double f2;  
   double fm;  
   int i;  
   int k;  
   double m2;  
   double m4;  
   double machep;  
   int ok;  
   double s;  
   double sf1;  
   double small;  
   double sx1;  
   double t2;  
   double temp;  
   double x2;  
   double xm;  
   
   machep = DBL_EPSILON;  
   small = machep * machep;  
   m2 = sqrt ( machep );  
   m4 = sqrt ( m2 );  
   sf1 = *f1;  
   sx1 = *x1;  
   k = 0;  
   xm = 0.0;  
   fm = *fx;  
   f0 = *fx;  
   dz = ( *d2 < machep );  
 /*  
   Find the step size.  
 */  
   s = r8vec_norm ( n, x );  
   
   if ( dz )  
   {  
     temp = dmin;  
   }  
   else  
   {  
     temp = *d2;  
   }  
   
   t2 = m4 * sqrt ( fabs ( *fx ) / temp + s * ldt ) + m2 * ldt;  
   s = m4 * s + t;  
   if ( dz && s < t2 )  
   {  
     t2 = s;  
   }  
   
   t2 = fmax ( t2, small );  
   t2 = fmin ( t2, 0.01 * h );  
   
   if ( fk && *f1 <= fm )  
   {  
     xm = *x1;  
     fm = *f1;  
   }  
   
   if ( ( ! fk ) || fabs ( *x1 ) < t2 )  
   {  
     if ( 0.0 <= *x1 )  
     {  
       temp = 1.0;  
     }  
     else  
     {  
       temp = - 1.0;  
     }  
   
     *x1 = temp * t2;  
     *f1 = flin ( n, jsearch, *x1, func, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc );  
     /* *f1 = flin ( n, jsearch, *x1, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); */  
   }  
   
   if ( *f1 <= fm )  
   {  
     xm = *x1;  
     fm = *f1;  
   }  
 /*  
   Evaluate FLIN at another point and estimate the second derivative.  
 */  
   for ( ; ; )  
   {  
     if ( dz )  
     {  
       if ( *f1 <= f0 )  
       {  
         x2 = 2.0 * *x1;  
       }  
       else  
       {  
         x2 = - *x1;  
       }  
   
       f2 = flin ( n, jsearch, x2, func, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc );  
       /* f2 = flin ( n, jsearch, x2, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); */  
   
       if ( f2 <= fm )  
       {  
         xm = x2;  
         fm = f2;  
       }  
   
       *d2 = ( x2 * ( *f1 - f0 ) - *x1 * ( f2 - f0 ) )  
         / ( ( *x1 * x2 ) * ( *x1 - x2 ) );  
     }  
 /*  
   Estimate the first derivative at 0.  
 */  
     d1 = ( *f1 - f0 ) / *x1 - *x1 * *d2;  
     dz = 1;  
 /*  
   Predict the minimum.  
 */  
     if ( *d2 <= small )  
     {  
       if ( 0.0 <= d1 )  
       {  
         x2 = - h;  
       }  
       else  
       {  
         x2 = h;  
       }  
     }  
     else  
     {  
       x2 = ( - 0.5 * d1 ) / *d2;  
     }  
   
     if ( h < fabs ( x2 ) )  
     {  
       if ( x2 <= 0.0 )  
       {  
         x2 = - h;  
       }  
       else  
       {  
         x2 = h;  
       }  
     }  
 /*  
   Evaluate F at the predicted minimum.  
 */  
     ok = 1;  
   
     for ( ; ; )  
     {  
       f2 = flin ( n, jsearch, x2, func, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc );  
       /* f2 = flin ( n, jsearch, x2, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc ); */  
   
       if ( nits <= k || f2 <= f0 )  
       {  
         break;  
       }  
   
       k = k + 1;  
   
       if ( f0 < *f1 && 0.0 < *x1 * x2 )  
       {  
         ok = 0;  
         break;  
       }  
       x2 = 0.5 * x2;  
     }  
   
     if ( ok )  
     {  
       break;  
     }  
   }  
 /*  
   Increment the one-dimensional search counter.  
 */  
   *nl = *nl + 1;  
   
   if ( fm < f2 )  
   {  
     x2 = xm;  
   }  
   else  
   {  
     fm = f2;  
   }  
 /*  
   Get a new estimate of the second derivative.  
 */  
   if ( small < fabs ( x2 * ( x2 - *x1 ) ) )  
   {  
     *d2 = ( x2 * ( *f1 - f0 ) - *x1 * ( fm - f0 ) )   
       / ( ( *x1 * x2 ) * ( *x1 - x2 ) );  
   }  
   else  
   {  
     if ( 0 < k )  
     {  
       *d2 = 0.0;  
     }  
   }  
   
   *d2 = fmax ( *d2, small );  
   
   *x1 = x2;  
   *fx = fm;  
   
   if ( sf1 < *fx )  
   {  
     *fx = sf1;  
     *x1 = sx1;  
   }  
 /*  
   Update X for linear search.  
 */  
   if ( 0 <= jsearch )  
   {  
     for ( i = 0; i < n; i++ )  
     {  
       x[i] = x[i] + *x1 * v[i+jsearch*n];  
     }  
   }  
   
   return;  
 }  
 /******************************************************************************/  
   
 /* double praxis ( double t0, double h0, int n, int prin, double x[],  */  
 /*   double f ( double x[], int n ) ) */  
 double praxis ( double t0, double h0, int n, int prin, double x[],   
                 double (*func) ( double [] ))  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     PRAXIS seeks an N-dimensional minimizer X of a scalar function F(X).  
   
   Discussion:  
   
     PRAXIS returns the minimum of the function F(X,N) of N variables  
     using the principal axis method.  The gradient of the function is  
     not required.  
   
     The approximating quadratic form is  
   
       Q(x") = F(x,n) + (1/2) * (x"-x)" * A * (x"-x)  
   
     where X is the best estimate of the minimum and   
   
       A = inverse(V") * D * inverse(V)  
   
     V(*,*) is the matrix of search directions;   
     D(*) is the array of second differences.    
   
     If F(X) has continuous second derivatives near X0, then A will tend   
     to the hessian of F at X0 as X approaches X0.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     03 August 2016  
   
   Author:  
   
     Original FORTRAN77 version by Richard Brent.  
     C version by John Burkardt.  
   
   Reference:  
   
     Richard Brent,  
     Algorithms for Minimization with Derivatives,  
     Prentice Hall, 1973,  
     Reprinted by Dover, 2002.  
   
   Parameters:  
   
     Input, double T0, is a tolerance.  PRAXIS attempts to return   
     praxis = f(x) such that if X0 is the true local minimum near X, then  
     norm ( x - x0 ) < T0 + sqrt ( EPSILON ) * norm ( X ),  
     where EPSILON  is the machine precision.  
   
     Input, double H0, is the maximum step size.  H0 should be   
     set to about the maximum distance from the initial guess to the minimum.  
     If H0 is set too large or too small, the initial rate of  
     convergence may be slow.  
   
     Input, int N, the number of variables.  
   
     Input, int PRIN, controls printing intermediate results.  
     0, nothing is printed.  
     1, F is printed after every n+1 or n+2 linear minimizations.    
        final X is printed, but intermediate X is printed only   
        if N is at most 4.  
     2, the scale factors and the principal values of the approximating   
        quadratic form are also printed.  
     3, X is also printed after every few linear minimizations.  
     4, the principal vectors of the approximating quadratic form are   
        also printed.  
   
     Input/output, double X[N], is an array containing on entry a  
     guess of the point of minimum, on return the estimated point of minimum.  
   
     Input, double F ( double X[], int N ), is the name of the function to be  
     minimized.  
   
     Output, double PRAXIS, the function value at the minimizer.  
   
   Local parameters:  
   
     Local, double DMIN, an estimate for the smallest eigenvalue.  
   
     Local, double FX, the value of F(X,N).  
   
     Local, int ILLC, is TRUE if the system is ill-conditioned.  
   
     Local, double LDT, the length of the step.  
   
     Local, int NF, the number of function evaluations.  
   
     Local, int NL, the number of linear searches.  
 */  
 {  
   int biter=0; /* Added to count the loops */  
   double *d;  
   double d2;  
   double df;  
   double dmin;  
   double dn;  
   double dni;  
   double f1;  
   int fk;  
   double fx;  
   double h;  
   int i;  
   int illc;  
   int j;  
   int jsearch;  
   int k;  
   int k2;  
   int kl;  
   int kt;  
   int ktm;  
   double large;  
   double ldfac;  
   double lds;  
   double ldt;  
   double m2;  
   double m4;  
   double machep;  
   int nits;  
   int nl;  
   int nf;  
   double *q0;  
   double *q1;  
   double qa;  
   double qb;  
   double qc;  
   double qd0;  
   double qd1;  
   double qf1;  
   double r;  
   double s;  
   double scbd;  
   int seed;  
   double sf;  
   double sl;  
   double small;  
   double t;  
   double temp;  
   double t2;  
   double *v;  
   double value;  
   double vlarge;  
   double vsmall;  
   double *y;  
   double *z;  
 /*  
   Allocation.  
 */  
   d = ( double * ) malloc ( n * sizeof ( double ) );  
   q0 = ( double * ) malloc ( n * sizeof ( double ) );  
   q1 = ( double * ) malloc ( n * sizeof ( double ) );  
   v = ( double * ) malloc ( n * n * sizeof ( double ) );  
   y = ( double * ) malloc ( n * sizeof ( double ) );  
   z = ( double * ) malloc ( n * sizeof ( double ) );  
 /*  
   Initialization.  
 */  
   machep = DBL_EPSILON;  
   small = machep * machep;  
   vsmall = small * small;  
   large = 1.0 / small;  
   vlarge = 1.0 / vsmall;  
   m2 = sqrt ( machep );  
   m4 = sqrt ( m2 );  
   seed = 123456789;  
 /*  
   Heuristic numbers:  
   
   If the axes may be badly scaled (which is to be avoided if  
   possible), then set SCBD = 10.  Otherwise set SCBD = 1.  
   
   If the problem is known to be ill-conditioned, initialize ILLC = true.  
   
   KTM is the number of iterations without improvement before the  
   algorithm terminates.  KTM = 4 is very cautious; usually KTM = 1  
   is satisfactory.  
 */  
   scbd = 1.0;  
   illc = 0;  
   ktm = 1;  
   
   if ( illc )  
   {  
     ldfac = 0.1;  
   }  
   else  
   {  
     ldfac = 0.01;  
   }  
   
   kt = 0;  
   nl = 0;  
   nf = 1;  
   /* fx = f ( x, n ); */  
   fx = (*func) ( (x-1) );/* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/  
   qf1 = fx;  
   t = small + fabs ( t0 );  
   t2 = t;  
   dmin = small;  
   h = h0;  
   h = fmax ( h, 100.0 * t );  
   ldt = h;  
 /*  
   The initial set of search directions V is the identity matrix.  
 */  
   for ( j = 0; j < n; j++ )  
   {  
     for ( i = 0; i < n; i++ )  
     {  
       v[i+j*n] = 0.0;  
     }  
     v[j+j*n] = 1.0;  
   }  
   
   for ( i = 0; i < n; i++ )  
   {  
     d[i] = 0.0;  
   }  
   qa = 0.0;  
   qb = 0.0;  
   qc = 0.0;  
   qd0 = 0.0;  
   qd1 = 0.0;  
   r8vec_copy ( n, x, q0 );  
   r8vec_copy ( n, x, q1 );  
   
   if ( 0 < prin )  
   {  
     print2 ( n, x, prin, fx, nf, nl );  
   }  
 /*  
   The main loop starts here.  
 */  
   for ( ; ; )  
   {  
     biter++;  /* Added to count the loops */  
     printf("\n Big iteration %d \n",biter);  
     sf = d[0];  
     d[0] = 0.0;  
 /*  
   Minimize along the first direction V(*,1).  
 */  
     jsearch = 0;  
     nits = 2;  
     d2 = d[0];  
     s = 0.0;  
     value = fx;  
     fk = 0;  
   
     minny ( n, jsearch, nits, &d2, &s, &value, fk, func, x, t,   
       h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 );  
     /* minny ( n, jsearch, nits, &d2, &s, &value, fk, func, x, t,  */  
     /*   h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); */  
   
     d[0] = d2;  
   
     if ( s <= 0.0 )  
     {  
       for ( i = 0; i < n; i++ )  
       {  
         v[i+0*n] = - v[i+0*n];  
       }  
     }  
   
     if ( sf <= 0.9 * d[0] || d[0] <= 0.9 * sf )  
     {  
       for ( i = 1; i < n; i++ )  
       {  
         d[i] = 0.0;  
       }  
     }  
 /*  
   The inner loop starts here.  
 */  
     for ( k = 2; k <= n; k++ )  
     {  
       r8vec_copy ( n, x, y );  
   
       sf = fx;  
   
       if ( 0 < kt )  
       {  
         illc = 1;  
       }  
   
       for ( ; ; )  
       {  
         kl = k;  
         df = 0.0;  
 /*  
   A random step follows, to avoid resolution valleys.  
 */  
         if ( illc )  
         {  
           for ( j = 0; j < n; j++ )  
           {  
             r = r8_uniform_01 ( &seed );  
             s = ( 0.1 * ldt + t2 * pow ( 10.0, kt ) ) * ( r - 0.5 );  
             z[j] = s;  
             for ( i = 0; i < n; i++ )  
             {  
               x[i] = x[i] + s * v[i+j*n];  
             }  
           }  
   
           fx = (*func) ( (x-1) );/* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/  
           /* fx = f ( x, n ); */  
           nf = nf + 1;  
         }  
 /*  
   Minimize along the "non-conjugate" directions V(*,K),...,V(*,N).  
 */  
         for ( k2 = k; k2 <= n; k2++ )  
         {  
           sl = fx;  
   
           jsearch = k2 - 1;  
           nits = 2;  
           d2 = d[k2-1];  
           s = 0.0;  
           value = fx;  
           fk = 0;  
   
           minny ( n, jsearch, nits, &d2, &s, &value, fk, func, x, t,   
             h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 );  
           /* minny ( n, jsearch, nits, &d2, &s, &value, fk, f, x, t,  */  
           /*   h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); */  
   
           d[k2-1] = d2;  
   
           if ( illc )  
           {  
             s = d[k2-1] * pow ( s + z[k2-1], 2 );  
           }  
           else  
           {  
             s = sl - fx;  
           }  
   
           if ( df <= s )  
           {  
             df = s;  
             kl = k2;  
           }  
         }  
 /*  
   If there was not much improvement on the first try, set  
   ILLC = true and start the inner loop again.  
 */  
         if ( illc )  
         {  
           break;  
         }  
         printf("\n fabs(  100.0 * machep(=%.12lf) * fx(=%.12lf) ) <=? df(=%.12lf)\n", machep, fx, df);  
         if ( fabs ( 100.0 * machep * fx ) <= df )  
         {  
           break;  
         }  
         illc = 1;  
       }  
   
       if ( k == 2 && 1 < prin )  
       {  
         r8vec_print ( n, d, "  The second difference array:" );  
       }  
 /*  
   Minimize along the "conjugate" directions V(*,1),...,V(*,K-1).  
 */  
       for ( k2 = 1; k2 < k; k2++ )  
       {  
         jsearch = k2 - 1;  
         nits = 2;  
         d2 = d[k2-1];  
         s = 0.0;  
         value = fx;  
         fk = 0;  
   
         minny ( n, jsearch, nits, &d2, &s, &value, fk, func, x, t,   
           h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 );  
         /* minny ( n, jsearch, nits, &d2, &s, &value, fk, f, x, t,  */  
         /*   h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); */  
   
         d[k2-1] = d2;  
       }  
    
       f1 = fx;  
       fx = sf;  
   
       for ( i = 0; i < n; i++ )  
       {  
         temp = x[i];  
         x[i] = y[i];  
         y[i] = temp - y[i];  
       }  
         
       lds = r8vec_norm ( n, y );  
 /*  
   Discard direction V(*,kl).  
   
   If no random step was taken, V(*,KL) is the "non-conjugate"  
   direction along which the greatest improvement was made.  
 */  
       if ( small < lds )  
       {  
         for ( j = kl - 1; k <= j; j-- )  
         {  
           for ( i = 1; i <= n; i++ )  
           {  
             v[i-1+j*n] = v[i-1+(j-1)*n];  
           }  
           d[j] = d[j-1];  
         }  
   
         d[k-1] = 0.0;  
   
         for ( i = 1; i <= n; i++ )  
         {  
           v[i-1+(k-1)*n] = y[i-1] / lds;  
         }  
 /*  
   Minimize along the new "conjugate" direction V(*,k), which is  
   the normalized vector:  (new x) - (old x).  
 */  
         jsearch = k - 1;  
         nits = 4;  
         d2 = d[k-1];  
         value = f1;  
         fk = 1;  
   
         minny ( n, jsearch, nits, &d2, &lds, &value, fk, func, x, t,   
           h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 );  
         /* minny ( n, jsearch, nits, &d2, &lds, &value, fk, f, x, t,  */  
         /*   h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qa, &qb, &qc, &qd0, &qd1 ); */  
   
         d[k-1] = d2;  
   
         if ( lds <= 0.0 )  
         {  
           lds = - lds;  
           for ( i = 1; i <= n; i++ )  
           {  
             v[i-1+(k-1)*n] = - v[i-1+(k-1)*n];  
           }  
         }  
       }  
   
       ldt = ldfac * ldt;  
       ldt = fmax ( ldt, lds );  
   
       if ( 0 < prin )  
       {  
         printf(" k=%d",k);   
         print2 ( n, x, prin, fx, nf, nl );  
       }  
   
       t2 = r8vec_norm ( n, x );  
   
       t2 = m2 * t2 + t;  
 /*  
   See whether the length of the step taken since starting the  
   inner loop exceeds half the tolerance.  
 */  
       if ( 0.5 * t2 < ldt )  
       {  
         kt = - 1;  
       }  
   
       kt = kt + 1;  
   
       if ( ktm < kt )  
       {  
         if ( 0 < prin )  
         {  
           r8vec_print ( n, x, "  X:" );  
         }  
   
         free ( d );  
         free ( q0 );  
         free ( q1 );  
         free ( v );  
         free ( y );  
         free ( z );  
   
         return fx;  
       }  
     }  
 /*  
   The inner loop ends here.  
   
   Try quadratic extrapolation in case we are in a curved valley.  
 */  
     quad ( n, func, x, t, h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qf1,   
       &qa, &qb, &qc, &qd0, &qd1 );  
     /* quad ( n, f, x, t, h, v, q0, q1, &nl, &nf, dmin, ldt, &fx, &qf1,  */  
     /*   &qa, &qb, &qc, &qd0, &qd1 ); */  
   
     for ( j = 0; j < n; j++ )  
     {  
       d[j] = 1.0 / sqrt ( d[j] );  
     }  
    
     dn = r8vec_max ( n, d );  
   
     if ( 3 < prin )  
     {  
       r8mat_print ( n, n, v, "  The new direction vectors:" );  
     }  
   
     for ( j = 0; j < n; j++ )  
     {  
       for ( i = 0; i < n; i++ )  
       {  
         v[i+j*n] = ( d[j] / dn ) * v[i+j*n];  
       }  
     }  
 /*  
   Scale the axes to try to reduce the condition number.  
 */  
     if ( 1.0 < scbd )  
     {  
       for ( i = 0; i < n; i++ )  
       {  
         s = 0.0;  
         for ( j = 0; j < n; j++ )  
         {  
           s = s + v[i+j*n] * v[i+j*n];  
         }  
         s = sqrt ( s );  
         z[i] = fmax ( m4, s );  
       }  
   
       s = r8vec_min ( n, z );  
   
       for ( i = 0; i < n; i++ )  
       {  
         sl = s / z[i];  
         z[i] = 1.0 / sl;  
   
         if ( scbd < z[i] )  
         {  
           sl = 1.0 / scbd;  
           z[i] = scbd;  
         }  
         for ( j = 0; j < n; j++ )  
         {  
           v[i+j*n] = sl * v[i+j*n];  
         }  
       }  
     }  
 /*  
   Calculate a new set of orthogonal directions before repeating  
   the main loop.  
   
   Transpose V for MINFIT:  
 */  
     printf(" Calculate a new set of orthogonal directions before repeating  the main loop.\n  Transpose V for MINFIT:...\n");  
     r8mat_transpose_in_place ( n, v );  
 /*  
   MINFIT finds the singular value decomposition of V.  
   
   This gives the principal values and principal directions of the  
   approximating quadratic form without squaring the condition number.  
 */  
     printf(" MINFIT finds the singular value decomposition of V. \n This gives the principal values and principal directions of the\n  approximating quadratic form without squaring the condition number...\n");  
     minfit ( n, vsmall, v, d );  
 /*  
   Unscale the axes.  
 */  
     printf(" Unscale the axes.\n");  
     if ( 1.0 < scbd )  
     {  
       for ( i = 0; i < n; i++ )  
       {  
         for ( j = 0; j < n; j++ )  
         {  
           v[i+j*n] = z[i] * v[i+j*n];  
         }  
       }  
   
       for ( j = 0; j < n; j++ )  
       {  
         s = 0.0;  
         for ( i = 0; i < n; i++ )  
         {  
           s = s + v[i+j*n] * v[i+j*n];  
         }  
         s = sqrt ( s );  
   
         d[j] = s * d[j];  
         for ( i = 0; i < n; i++ )  
         {  
           v[i+j*n] = v[i+j*n] / s;  
         }  
       }  
     }  
   
     for ( i = 0; i < n; i++ )  
     {  
       dni = dn * d[i];  
   
       if ( large < dni )  
       {  
         d[i] = vsmall;  
       }  
       else if ( dni < small )  
       {  
         d[i] = vlarge;  
       }  
       else  
       {  
         d[i] = 1.0 / dni / dni;  
       }  
     }  
 /*  
   Sort the eigenvalues and eigenvectors.  
 */  
     printf(" Sort the eigenvalues and eigenvectors....\n");  
     svsort ( n, d, v );  
 /*  
   Determine the smallest eigenvalue.  
 */  
     printf("  Determine the smallest eigenvalue.\n");  
     dmin = fmax ( d[n-1], small );  
 /*  
   The ratio of the smallest to largest eigenvalue determines whether  
   the system is ill conditioned.  
 */  
       
     if ( dmin < m2 * d[0] )  
     {  
       illc = 1;  
     }  
     else  
     {  
       illc = 0;  
     }  
     printf("  The ratio of the smallest to largest eigenvalue determines whether\n  the system is ill conditioned=%d . dmin=%.12lf < m2=%.12lf * d[0]=%.12lf \n",illc, dmin,m2, d[0]);  
   
     if ( 1 < prin )  
     {  
       if ( 1.0 < scbd )  
       {  
         r8vec_print ( n, z, "  The scale factors:" );  
       }   
       r8vec_print ( n, d, "  Principal values of the quadratic form:" );  
     }  
   
     if ( 3 < prin )  
     {  
       r8mat_print ( n, n, v, "  The principal axes:" );  
     }  
 /*  
   The main loop ends here.  
 */  
   }  
   
   if ( 0 < prin )  
   {  
     r8vec_print ( n, x, "  X:" );  
   }  
 /*  
   Free memory.  
 */  
   free ( d );  
   free ( q0 );  
   free ( q1 );  
   free ( v );  
   free ( y );  
   free ( z );  
   
   return fx;  
 }  
 /******************************************************************************/  
   
 void print2 ( int n, double x[], int prin, double fx, int nf, int nl )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     PRINT2 prints certain data about the progress of the iteration.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     28 July 2016  
   
   Author:  
   
     Original FORTRAN77 version by Richard Brent.  
     C version by John Burkardt.  
   
   Reference:  
   
     Richard Brent,  
     Algorithms for Minimization with Derivatives,  
     Prentice Hall, 1973,  
     Reprinted by Dover, 2002.  
   
   Parameters:  
   
     Input, int N, the number of variables.  
   
     Input, double X[N], the current estimate of the minimizer.  
   
     Input, int PRIN, the user-specifed print level.  
     0, nothing is printed.  
     1, F is printed after every n+1 or n+2 linear minimizations.    
        final X is printed, but intermediate X is printed only   
        if N is at most 4.  
     2, the scale factors and the principal values of the approximating   
        quadratic form are also printed.  
     3, X is also printed after every few linear minimizations.  
     4, the principal vectors of the approximating quadratic form are   
        also printed.  
   
     Input, double FX, the smallest value of F(X) found so far.  
   
     Input, int NF, the number of function evaluations.  
   
     Input, int NL, the number of linear searches.  
 */  
 {  
   printf ( "\n" );  
   printf ( "  Linear searches      %d", nl );  
   /* printf ( "  Linear searches      %d\n", nl ); */  
   /* printf ( "  Function evaluations %d\n", nf ); */  
   /* printf ( "  Function value FX = %g\n", fx ); */  
   printf ( "  Function evaluations %d", nf );  
   printf ( "  Function value FX = %.12lf\n", fx );  
   
   if ( n <= 4 || 2 < prin )  
   {  
     r8vec_print ( n, x, "  X:" );  
   }  
   
   return;  
 }  
 /******************************************************************************/  
   
 void quad ( int n, double (*func) ( double [] ), double x[], double t,   
   double h, double v[], double q0[], double q1[], int *nl, int *nf, double dmin,   
   double ldt, double *fx, double *qf1, double *qa, double *qb, double *qc,   
   double *qd0, double *qd1 )  
 /* void quad ( int n, double f ( double x[], int n ), double x[], double t,  */  
 /*   double h, double v[], double q0[], double q1[], int *nl, int *nf, double dmin,  */  
 /*   double ldt, double *fx, double *qf1, double *qa, double *qb, double *qc,  */  
 /*   double *qd0, double *qd1 ) */  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     QUAD seeks to minimize the scalar function F along a particular curve.  
   
   Discussion:  
   
     The minimizer to be sought is required to lie on a curve defined  
     by Q0, Q1 and X.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     30 July 2016  
   
   Author:  
   
     Original FORTRAN77 version by Richard Brent.  
     C version by John Burkardt.  
   
   Reference:  
   
     Richard Brent,  
     Algorithms for Minimization with Derivatives,  
     Prentice Hall, 1973,  
     Reprinted by Dover, 2002.  
   
   Parameters:  
   
     Input, int N, the number of variables.  
   
     Input, double F ( double X[], int N ), the name of the function to   
     be minimized.  
   
     Input/output, double X[N], ?  
   
     Input, double T, ?  
   
     Input, double H, ?  
   
     Input, double V[N,N], the matrix of search directions.  
   
     Input/output, double Q0[N], Q1[N], two auxiliary points used to define  
     a curve through X.  
   
     Input/output, int *NL, the number of linear searches.  
   
     Input/output, int *NF, the number of function evaluations.  
   
     Input, double DMIN, an estimate for the smallest eigenvalue.  
   
     Input, double LDT, the length of the step.  
   
     Input/output, double *FX, the value of F(X,N).  
   
     Input/output, double *QF1, *QA, *QB, *QC, *QD0, *QD1 ?  
 */  
 {  
   int fk;  
   int i;  
   int jsearch;  
   double l;  
   int nits;  
   double s;  
   double temp;  
   double value;  
   
   temp = *fx;  
   *fx   = *qf1;  
   *qf1  = temp;  
   
   for ( i = 0; i < n; i++ )  
   {  
     temp  = x[i];  
     x[i]  = q1[i];  
     q1[i] = temp;  
   }  
   
   *qd1 = 0.0;  
   for ( i = 0; i < n; i++ )  
   {  
     *qd1 = *qd1 + ( x[i] - q1[i] ) * ( x[i] - q1[i] );  
   }  
   *qd1 = sqrt ( *qd1 );  
   
   if ( *qd0 <= 0.0 || *qd1 <= 0.0 || *nl < 3 * n * n )  
   {  
     *fx = *qf1;  
     *qa = 0.0;  
     *qb = 0.0;  
     *qc = 1.0;  
     s = 0.0;  
   }  
   else  
   {  
     jsearch = - 1;  
     nits = 2;  
     s = 0.0;  
     l = *qd1;  
     value = *qf1;  
     fk = 1;  
   
     minny ( n, jsearch, nits, &s, &l, &value, fk, func, x, t,   
       h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 );  
     /* minny ( n, jsearch, nits, &s, &l, &value, fk, f, x, t,  */  
     /*   h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 ); */  
   
     *qa =                  l * ( l - *qd1 )        / ( *qd0 + *qd1 ) / *qd0;  
     *qb = - ( l + *qd0 )     * ( l - *qd1 ) / *qd1                   / *qd0;  
     *qc =   ( l + *qd0 ) * l                / *qd1 / ( *qd0 + *qd1 );  
   }  
   
   *qd0 = *qd1;  
   
   for ( i = 0; i < n; i++ )  
   {  
     s = q0[i];  
     q0[i] = x[i];  
     x[i] = *qa * s + *qb * x[i] + *qc * q1[i];  
   }  
   
   return;  
 }  
 /******************************************************************************/  
   
 double r8_hypot ( double x, double y )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8_HYPOT returns the value of sqrt ( X^2 + Y^2 ).  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     26 March 2012  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, double X, Y, the arguments.  
   
     Output, double R8_HYPOT, the value of sqrt ( X^2 + Y^2 ).  
 */  
 {  
   double a;  
   double b;  
   double value;  
   
   if ( fabs ( x ) < fabs ( y ) )  
   {  
     a = fabs ( y );  
     b = fabs ( x );  
   }  
   else  
   {  
     a = fabs ( x );  
     b = fabs ( y );  
   }  
 /*  
   A contains the larger value.  
 */  
   if ( a == 0.0 )  
   {  
     value = 0.0;  
   }  
   else  
   {  
     value = a * sqrt ( 1.0 + ( b / a ) * ( b / a ) );  
   }  
   
   return value;  
 }  
 /******************************************************************************/  
   
 double r8_uniform_01 ( int *seed )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1].  
   
   Discussion:  
   
     This routine implements the recursion  
   
       seed = 16807 * seed mod ( 2^31 - 1 )  
       r8_uniform_01 = seed / ( 2^31 - 1 )  
   
     The integer arithmetic never requires more than 32 bits,  
     including a sign bit.  
   
     If the initial seed is 12345, then the first three computations are  
   
       Input     Output      R8_UNIFORM_01  
       SEED      SEED  
   
          12345   207482415  0.096616  
      207482415  1790989824  0.833995  
     1790989824  2035175616  0.947702  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     11 August 2004  
   
   Author:  
   
     John Burkardt  
   
   Reference:  
   
     Paul Bratley, Bennett Fox, Linus Schrage,  
     A Guide to Simulation,  
     Springer Verlag, pages 201-202, 1983.  
   
     Pierre L'Ecuyer,  
     Random Number Generation,  
     in Handbook of Simulation  
     edited by Jerry Banks,  
     Wiley Interscience, page 95, 1998.  
   
     Bennett Fox,  
     Algorithm 647:  
     Implementation and Relative Efficiency of Quasirandom  
     Sequence Generators,  
     ACM Transactions on Mathematical Software,  
     Volume 12, Number 4, pages 362-376, 1986.  
   
     P A Lewis, A S Goodman, J M Miller,  
     A Pseudo-Random Number Generator for the System/360,  
     IBM Systems Journal,  
     Volume 8, pages 136-143, 1969.  
   
   Parameters:  
   
     Input/output, int *SEED, the "seed" value.  Normally, this  
     value should not be 0.  On output, SEED has been updated.  
   
     Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between  
     0 and 1.  
 */  
 {  
   const int i4_huge = 2147483647;  
   int k;  
   double r;  
   
   if ( *seed == 0 )  
   {  
     fprintf ( stderr, "\n" );  
     fprintf ( stderr, "R8_UNIFORM_01 - Fatal error!\n" );  
     fprintf ( stderr, "  Input value of SEED = 0\n" );  
     exit ( 1 );  
   }  
   
   k = *seed / 127773;  
   
   *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;  
   
   if ( *seed < 0 )  
   {  
     *seed = *seed + i4_huge;  
   }  
   
   r = ( ( double ) ( *seed ) ) * 4.656612875E-10;  
   
   return r;  
 }  
 /******************************************************************************/  
   
 void r8mat_print ( int m, int n, double a[], char *title )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8MAT_PRINT prints an R8MAT.  
   
   Discussion:  
   
     An R8MAT is a doubly dimensioned array of R8 values, stored as a vector  
     in column-major order.  
   
     Entry A(I,J) is stored as A[I+J*M]  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     28 May 2008  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, int M, the number of rows in A.  
   
     Input, int N, the number of columns in A.  
   
     Input, double A[M*N], the M by N matrix.  
   
     Input, char *TITLE, a title.  
 */  
 {  
   r8mat_print_some ( m, n, a, 1, 1, m, n, title );  
   
   return;  
 }  
 /******************************************************************************/  
   
 void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,  
   int jhi, char *title )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8MAT_PRINT_SOME prints some of an R8MAT.  
   
   Discussion:  
   
     An R8MAT is a doubly dimensioned array of R8 values, stored as a vector  
     in column-major order.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     26 June 2013  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, int M, the number of rows of the matrix.  
     M must be positive.  
   
     Input, int N, the number of columns of the matrix.  
     N must be positive.  
   
     Input, double A[M*N], the matrix.  
   
     Input, int ILO, JLO, IHI, JHI, designate the first row and  
     column, and the last row and column to be printed.  
   
     Input, char *TITLE, a title.  
 */  
 {  
 # define INCX 5  
   
   int i;  
   int i2hi;  
   int i2lo;  
   int j;  
   int j2hi;  
   int j2lo;  
   
   fprintf ( stdout, "\n" );  
   fprintf ( stdout, "%s\n", title );  
   
   if ( m <= 0 || n <= 0 )  
   {  
     fprintf ( stdout, "\n" );  
     fprintf ( stdout, "  (None)\n" );  
     return;  
   }  
 /*  
   Print the columns of the matrix, in strips of 5.  
 */  
   for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )  
   {  
     j2hi = j2lo + INCX - 1;  
     if ( n < j2hi )  
     {  
       j2hi = n;  
     }  
     if ( jhi < j2hi )  
     {  
       j2hi = jhi;  
     }  
   
     fprintf ( stdout, "\n" );  
 /*  
   For each column J in the current range...  
   
   Write the header.  
 */  
     fprintf ( stdout, "  Col:  ");  
     for ( j = j2lo; j <= j2hi; j++ )  
     {  
       fprintf ( stdout, "  %7d     ", j - 1 );  
     }  
     fprintf ( stdout, "\n" );  
     fprintf ( stdout, "  Row\n" );  
     fprintf ( stdout, "\n" );  
 /*  
   Determine the range of the rows in this strip.  
 */  
     if ( 1 < ilo )  
     {  
       i2lo = ilo;  
     }  
     else  
     {  
       i2lo = 1;  
     }  
     if ( m < ihi )  
     {  
       i2hi = m;  
     }  
     else  
     {  
       i2hi = ihi;  
     }  
   
     for ( i = i2lo; i <= i2hi; i++ )  
     {  
 /*  
   Print out (up to) 5 entries in row I, that lie in the current strip.  
 */  
       fprintf ( stdout, "%5d:", i - 1 );  
       for ( j = j2lo; j <= j2hi; j++ )  
       {  
         fprintf ( stdout, "  %14g", a[i-1+(j-1)*m] );  
       }  
       fprintf ( stdout, "\n" );  
     }  
   }  
   
   return;  
 # undef INCX  
 }  
 /******************************************************************************/  
   
 void r8mat_transpose_in_place ( int n, double a[] )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8MAT_TRANSPOSE_IN_PLACE transposes a square matrix in place.  
   
   Discussion:  
   
     An R8MAT is a doubly dimensioned array of R8 values, stored as a vector  
     in column-major order.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     26 June 2008  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, int N, the number of rows and columns of the matrix A.  
   
     Input/output, double A[N*N], the matrix to be transposed.  
 */  
 {  
   int i;  
   int j;  
   double t;  
   
   for ( j = 0; j < n; j++ )  
   {  
     for ( i = 0; i < j; i++ )  
     {  
       t        = a[i+j*n];  
       a[i+j*n] = a[j+i*n];  
       a[j+i*n] = t;  
     }  
   }  
   return;  
 }  
 /******************************************************************************/  
   
 void r8vec_copy ( int n, double a1[], double a2[] )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8VEC_COPY copies an R8VEC.  
   
   Discussion:  
   
     An R8VEC is a vector of R8's.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     03 July 2005  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, int N, the number of entries in the vectors.  
   
     Input, double A1[N], the vector to be copied.  
   
     Input, double A2[N], the copy of A1.  
 */  
 {  
   int i;  
   
   for ( i = 0; i < n; i++ )  
   {  
     a2[i] = a1[i];  
   }  
   return;  
 }  
 /******************************************************************************/  
   
 double r8vec_max ( int n, double r8vec[] )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8VEC_MAX returns the value of the maximum element in a R8VEC.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     05 May 2006  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, int N, the number of entries in the array.  
   
     Input, double R8VEC[N], a pointer to the first entry of the array.  
   
     Output, double R8VEC_MAX, the value of the maximum element.  This  
     is set to 0.0 if N <= 0.  
 */  
 {  
   int i;  
   double value;  
   
   if ( n <= 0 )  
   {  
     value = 0.0;  
     return value;  
   }  
   
   value = r8vec[0];  
   
   for ( i = 1; i < n; i++ )  
   {  
     if ( value < r8vec[i] )  
     {  
       value = r8vec[i];  
     }  
   }  
   return value;  
 }  
 /******************************************************************************/  
   
 double r8vec_min ( int n, double r8vec[] )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8VEC_MIN returns the value of the minimum element in a R8VEC.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     05 May 2006  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, int N, the number of entries in the array.  
   
     Input, double R8VEC[N], the array to be checked.  
   
     Output, double R8VEC_MIN, the value of the minimum element.  
 */  
 {  
   int i;  
   double value;  
   
   value = r8vec[0];  
   
   for ( i = 1; i < n; i++ )  
   {  
     if ( r8vec[i] < value )  
     {  
       value = r8vec[i];  
     }  
   }  
   return value;  
 }  
 /******************************************************************************/  
   
 double r8vec_norm ( int n, double a[] )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8VEC_NORM returns the L2 norm of an R8VEC.  
   
   Discussion:  
   
     The vector L2 norm is defined as:  
   
       R8VEC_NORM = sqrt ( sum ( 1 <= I <= N ) A(I)^2 ).  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     01 March 2003  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, int N, the number of entries in A.  
   
     Input, double A[N], the vector whose L2 norm is desired.  
   
     Output, double R8VEC_NORM, the L2 norm of A.  
 */  
 {  
   int i;  
   double v;  
   
   v = 0.0;  
   
   for ( i = 0; i < n; i++ )  
   {  
     v = v + a[i] * a[i];  
   }  
   v = sqrt ( v );  
   
   return v;  
 }  
 /******************************************************************************/  
   
 void r8vec_print ( int n, double a[], char *title )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     R8VEC_PRINT prints an R8VEC.  
   
   Discussion:  
   
     An R8VEC is a vector of R8's.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     08 April 2009  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     Input, int N, the number of components of the vector.  
   
     Input, double A[N], the vector to be printed.  
   
     Input, char *TITLE, a title.  
 */  
 {  
   int i,j, jk, k;  
   
   double *p;  
   
   p=(a-1); /* So that a[0]=p[1]  */  
   /* for (i=1;i<=n;i++) { */  
   /*   fprintf(ficrespow," %.12lf", p[i]); */  
   /* } */  
   /* fprintf(ficrespow,"\n");fflush(ficrespow); */  
   printf("\n#model=  1      +     age ");  
   fprintf(ficlog,"\n#model=  1      +     age ");  
   if(nagesqr==1){  
     printf("  + age*age  ");  
     fprintf(ficlog,"  + age*age  ");  
   }  
   for(j=1;j <=ncovmodel-2;j++){  
     if(Typevar[j]==0) {  
       printf("  +      V%d  ",Tvar[j]);  
       fprintf(ficlog,"  +      V%d  ",Tvar[j]);  
     }else if(Typevar[j]==1) {  
       printf("  +    V%d*age ",Tvar[j]);  
       fprintf(ficlog,"  +    V%d*age ",Tvar[j]);  
     }else if(Typevar[j]==2) {  
       printf("  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);  
       fprintf(ficlog,"  +    V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);  
     }else if(Typevar[j]==3) {  
       printf("  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);  
       fprintf(ficlog,"  +    V%d*V%d*age ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);  
     }  
   }  
   printf("\n");  
 /*     printf("12   47.0114589    0.0154322   33.2424412    0.3279905    2.3731903  */  
 /* 13  -21.5392400    0.1118147    1.2680506    1.2973408   -1.0663662  */  
     fprintf(ficlog,"\n");  
     for(i=1,jk=1; i <=nlstate; i++){  
       for(k=1; k <=(nlstate+ndeath); k++){  
         if (k != i) {  
           printf("%d%d ",i,k);  
           fprintf(ficlog,"%d%d ",i,k);  
           for(j=1; j <=ncovmodel; j++){  
             printf("%12.7f ",p[jk]);  
             fprintf(ficlog,"%12.7f ",p[jk]);  
             jk++;   
           }  
           printf("\n");  
           fprintf(ficlog,"\n");  
         }  
       }  
     }  
  /* fprintf ( stdout, "\n" ); */  
   /* fprintf ( stdout, " %s\n", title ); */  
   fprintf ( stdout, " %s", title );  
   /* fprintf ( stdout, "\n" ); */  
   for ( i = 0; i < n; i++ )  
   {  
     /* fprintf ( stdout, "  %8d: %14g", i+1, a[i] ); */  
     fprintf ( stdout, "  %.12lf",  a[i] );  
   }  
   fprintf ( stdout, "\n" );  
   /* for ( i = 0; i < n; i++ ) */  
   /* { */  
   /*   fprintf ( stdout, "  %8d: %14g\n", i, a[i] ); */  
   /* } */  
   
   return;  
 }  
 /******************************************************************************/  
   
 void svsort ( int n, double d[], double v[] )   
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     SVSORT descending sorts D and adjusts the corresponding columns of V.  
   
   Discussion:  
   
     A simple bubble sort is used on D.  
   
     In our application, D contains singular values, and the columns of V are  
     the corresponding right singular vectors.  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     28 July 2016  
   
   Author:  
   
     Original FORTRAN77 version by Richard Brent.  
     C version by John Burkardt.  
   
   Reference:  
   
     Richard Brent,  
     Algorithms for Minimization with Derivatives,  
     Prentice Hall, 1973,  
     Reprinted by Dover, 2002.  
   
   Parameters:  
   
     Input, int N, the length of D, and the order of V.  
   
     Input/output, double D[N], the vector to be sorted.    
     On output, the entries of D are in descending order.  
   
     Input/output, double V[N,N], an N by N array to be adjusted   
     as D is sorted.  In particular, if the value that was in D(I) on input is  
     moved to D(J) on output, then the input column V(*,I) is moved to  
     the output column V(*,J).  
 */  
 {  
   int i;  
   int j1;  
   int j2;  
   int j3;  
   double t;  
   
   for ( j1 = 0; j1 < n - 1; j1++ )  
   {  
 /*  
   Find J3, the index of the largest entry in D[J1:N-1].  
   MAXLOC apparently requires its output to be an array.  
 */  
     j3 = j1;  
     for ( j2 = j1 + 1; j2 < n; j2++ )  
     {  
       if ( d[j3] < d[j2] )  
       {  
         j3 = j2;  
       }  
     }  
 /*  
   If J1 != J3, swap D[J1] and D[J3], and columns J1 and J3 of V.  
 */  
     if ( j1 != j3 )  
     {  
       t     = d[j1];  
       d[j1] = d[j3];  
       d[j3] = t;  
       for ( i = 0; i < n; i++ )  
       {  
         t         = v[i+j1*n];  
         v[i+j1*n] = v[i+j3*n];  
         v[i+j3*n] = t;  
       }  
     }  
   }  
   
   return;  
 }  
 /******************************************************************************/  
   
 void timestamp ( )  
   
 /******************************************************************************/  
 /*  
   Purpose:  
   
     TIMESTAMP prints the current YMDHMS date as a time stamp.  
   
   Example:  
   
     31 May 2001 09:45:54 AM  
   
   Licensing:  
   
     This code is distributed under the GNU LGPL license.  
   
   Modified:  
   
     24 September 2003  
   
   Author:  
   
     John Burkardt  
   
   Parameters:  
   
     None  
 */  
 {  
 # define TIME_SIZE 40  
   
   static char time_buffer[TIME_SIZE];  
   const struct tm *tm;  
   time_t now;  
   
   now = time ( NULL );  
   tm = localtime ( &now );  
   
   strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );  
   
   fprintf ( stdout, "%s\n", time_buffer );  
   
   return;  
 # undef TIME_SIZE  
 }  
 /* end praxis */  /* end praxis */
   
 /*************** powell ************************/  /*************** powell ************************/
Line 5187  void powell(double p[], double **xi, int Line 2898  void powell(double p[], double **xi, int
   double fp,fptt;    double fp,fptt;
   double *xits;    double *xits;
   int niterf, itmp;    int niterf, itmp;
     int Bigter=0, nBigterf=1;
     
   pt=vector(1,n);     pt=vector(1,n); 
   ptt=vector(1,n);     ptt=vector(1,n); 
   xit=vector(1,n);     xit=vector(1,n); 
Line 5200  void powell(double p[], double **xi, int Line 2912  void powell(double p[], double **xi, int
     ibig=0;       ibig=0; 
     del=0.0;       del=0.0; 
     rlast_time=rcurr_time;      rlast_time=rcurr_time;
       rlast_btime=rcurr_time;
     /* (void) gettimeofday(&curr_time,&tzp); */      /* (void) gettimeofday(&curr_time,&tzp); */
     rcurr_time = time(NULL);        rcurr_time = time(NULL);  
     curr_time = *localtime(&rcurr_time);      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); */      /* 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); */      /* 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);      Bigter=(*iter - *iter % ncovmodel)/ncovmodel +1; /* Big iteration, i.e on ncovmodel cycle */
     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);      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(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */      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 */      fp=(*fret); /* From former iteration or initial value */
     for (i=1;i<=n;i++) {      for (i=1;i<=n;i++) {
       fprintf(ficrespow," %.12lf", p[i]);        fprintf(ficrespow," %.12lf", p[i]);
Line 5262  void powell(double p[], double **xi, int Line 2976  void powell(double p[], double **xi, int
         strcurr[itmp-1]='\0';          strcurr[itmp-1]='\0';
       printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);        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);        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);          rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
         forecast_time = *localtime(&rforecast_time);          forecast_time = *localtime(&rforecast_time);
         strcpy(strfor,asctime(&forecast_time));          strcpy(strfor,asctime(&forecast_time));
         itmp = strlen(strfor);          itmp = strlen(strfor);
         if(strfor[itmp-1]=='\n')          if(strfor[itmp-1]=='\n')
           strfor[itmp-1]='\0';            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);          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 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 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 (i=1;i<=n;i++) { /* For each direction i */
Line 6568  double func( double *x) Line 4284  double func( double *x)
   
   for(k=1; k<=nlstate; k++) ll[k]=0.;    for(k=1; k<=nlstate; k++) ll[k]=0.;
   ioffset=0;    ioffset=0;
   for (i=1,ipmx=0, sw=0.; i<=imx; i++){    if(mle==1){
     /* Computes the values of the ncovmodel covariates of the model      for (i=1,ipmx=0, sw=0.; i<=imx; i++){
        depending if the covariates are fixed or varying (age dependent) and stores them in cov[]        /* Computes the values of the ncovmodel covariates of the model
        Then computes with function pmij which return a matrix p[i][j] giving the elementary probability           depending if the covariates are fixed or varying (age dependent) and stores them in cov[]
        to be observed in j being in i according to the model.           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 */        ioffset=2+nagesqr ;
     for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummy or quant or prod */     /* Fixed */
       /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */        for (kf=1; kf<=ncovf;kf++){ /* For each fixed covariate dummy or quant or prod */
       /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */          /* # V1=sex, V2=raedyrs Quant Fixed, State=livarnb4..livarnb11, V3=iadl4..iald11, V4=adlw4..adlw11, V5=r4bmi..r11bmi */
       /*  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 */          /*             V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */
       /* TvarFind;  TvarFind[1]=6,  TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod)  */          /*  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 */
       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)*/          /* TvarFind;  TvarFind[1]=6,  TvarFind[2]=7, TvarFind[3]=9 *//* Inverse V2(6) is first fixed (single or prod)  */
       /* V1*V2 (7)  TvarFind[2]=7, TvarFind[3]=9 */          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         /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] 
        has been calculated etc */           is 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6 
     /* For an individual i, wav[i] gives the number of effective waves */           has been calculated etc */
     /* We compute the contribution to Likelihood of each effective transition        /* For an individual i, wav[i] gives the number of effective waves */
        mw[mi][i] is real wave of the mi th effectve wave */        /* We compute the contribution to Likelihood of each effective transition
     /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];           mw[mi][i] is real wave of the mi th effectve wave */
        s2=s[mw[mi+1][i]][i];        /* Then statuses are computed at each begin and end of an effective wave s1=s[ mw[mi][i] ][i];
        And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] because now is moved after nvocol+nqv            s2=s[mw[mi+1][i]][i];
        But if the variable is not in the model TTvar[iv] is the real variable effective in the model:           And the iv th varying covariate is the cotvar[mw[mi+1][i]][iv][i] because now is moved after nvocol+nqv 
        meaning that decodemodel should be used cotvar[mw[mi+1][i]][TTvar[iv]][i]           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 */        */
         for(mi=1; mi<= wav[i]-1; mi++){  /* Varying with waves */
       /* Wave varying (but not age varying) */        /* 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*\/ */          /* 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]]][i]; but where is the crossproduct? *\/ */
       /*   cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i]; */          /*   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 )*/          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 */            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*/            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 */            if(FixedV[itv]!=0){ /* Not a fixed covariate */
           cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i];  /* cotvar[wav][ncovcol+nqv+iv][i] */              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) */   
           }else{ /* fixed covariate */            }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=covar[itv][i];  /* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model */
           }            }
Line 6665  double func( double *x) Line 4333  double func( double *x)
             cotvarv=cotvarv*cotvarvold;              cotvarv=cotvarv*cotvarvold;
           }            }
           iposold=ipos;            iposold=ipos;
           cov[ioffset+ipos]=cotvarv*agexact;            cov[ioffset+ipos]=cotvarv;
           /* For products */  
         }          }
           /* 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 */          /*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.          /* 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            * 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            * (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           * 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           * 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           * (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           * 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                                   * from savm to out if bh is negative or even beyond if bh is positive. bh varies
        * -stepm/2 to stepm/2 .                                   * -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 the same as for previous versions of Imach.
        * For stepm > 1 the results are less biased than in previous versions.                                    * For stepm > 1 the results are less biased than in previous versions. 
        */                                   */
       s1=s[mw[mi][i]][i];          s1=s[mw[mi][i]][i];
       s2=s[mw[mi+1][i]][i];          s2=s[mw[mi+1][i]][i];
       bbh=(double)bh[mi][i]/(double)stepm;           bbh=(double)bh[mi][i]/(double)stepm; 
       /* bias bh is positive if real duration          /* bias bh is positive if real duration
        * is higher than the multiple of stepm and negative otherwise.           * 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]));*/          /* 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){           if( s2 > nlstate){ 
         /* i.e. if s2 is a death state and if the date of death is known             /* 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                then the contribution to the likelihood is the probability to 
            die between last step unit time and current  step unit time,                die between last step unit time and current  step unit time, 
            which is also equal to probability to die before dh                which is also equal to probability to die before dh 
            minus probability to die before dh-stepm .                minus probability to die before dh-stepm . 
            In version up to 0.92 likelihood was computed               In version up to 0.92 likelihood was computed
            as if date of death was unknown. Death was treated as any other               as if date of death was unknown. Death was treated as any other
            health state: the date of the interview describes the actual state               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               and not the date of a change in health state. The former idea was
            to consider that at each interview the state was recorded               to consider that at each interview the state was recorded
            (healthy, disable or death) and IMaCh was corrected; but when we               (healthy, disable or death) and IMaCh was corrected; but when we
            introduced the exact date of death then we should have modified               introduced the exact date of death then we should have modified
            the contribution of an exact death to the likelihood. This new               the contribution of an exact death to the likelihood. This new
            contribution is smaller and very dependent of the step unit               contribution is smaller and very dependent of the step unit
            stepm. It is no more the probability to die between last interview               stepm. It is no more the probability to die between last interview
            and month of death but the probability to survive from last               and month of death but the probability to survive from last
            interview up to one month before death multiplied by the               interview up to one month before death multiplied by the
            probability to die within a month. Thanks to Chris               probability to die within a month. Thanks to Chris
            Jackson for correcting this bug.  Former versions increased               Jackson for correcting this bug.  Former versions increased
            mortality artificially. The bad side is that we add another loop               mortality artificially. The bad side is that we add another loop
            which slows down the processing. The difference can be up to 10%               which slows down the processing. The difference can be up to 10%
            lower mortality.               lower mortality.
         */            */
         /* If, at the beginning of the maximization mostly, the            /* If, at the beginning of the maximization mostly, the
            cumulative probability or probability to be dead is               cumulative probability or probability to be dead is
            constant (ie = 1) over time d, the difference is equal to               constant (ie = 1) over time d, the difference is equal to
            0.  out[s1][3] = savm[s1][3]: probability, being at state               0.  out[s1][3] = savm[s1][3]: probability, being at state
            s1 at precedent wave, to be dead a month before current               s1 at precedent wave, to be dead a month before current
            wave is equal to probability, being at state s1 at               wave is equal to probability, being at state s1 at
            precedent wave, to be dead at mont of the current               precedent wave, to be dead at mont of the current
            wave. Then the observed probability (that this person died)               wave. Then the observed probability (that this person died)
            is null according to current estimated parameter. In fact,               is null according to current estimated parameter. In fact,
            it should be very low but not zero otherwise the log go to               it should be very low but not zero otherwise the log go to
            infinity.               infinity.
         */            */
 /* #ifdef INFINITYORIGINAL */  /* #ifdef INFINITYORIGINAL */
 /*          lli=log(out[s1][s2] - savm[s1][s2]); */  /*          lli=log(out[s1][s2] - savm[s1][s2]); */
 /* #else */  /* #else */
Line 6739  double func( double *x) Line 4444  double func( double *x)
 /*        else */  /*        else */
 /*          lli=log(out[s1][s2] - savm[s1][s2]); */  /*          lli=log(out[s1][s2] - savm[s1][s2]); */
 /* #endif */  /* #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 */            out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
         for (j=1,survp=0. ; j<=nlstate; j++)                          1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];            savm=oldm;
         /*survp += out[s1][j]; */            oldm=newm;
         lli= log(survp);          } /* end mult */
       }        
       /* else if  (s2==-4) {  */          s1=s[mw[mi][i]][i];
       /*   for (j=3,survp=0. ; j<=nlstate; j++)   */          s2=s[mw[mi+1][i]][i];
       /*     survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */          if( s2 > nlstate){ 
       /*   lli= log(survp);  */            lli=log(out[s1][s2] - savm[s1][s2]);
       /* }  */          } else if  ( s2==-1 ) { /* alive */
       /* else if  (s2==-5) {  */            for (j=1,survp=0. ; j<=nlstate; j++) 
       /*   for (j=1,survp=0. ; j<=2; j++)   */              survp += out[s1][j];
       /*     survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; */            lli= log(survp);
       /*   lli= log(survp);  */          }else{
       /* }  */            lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
       else if (mle==1){          }
         lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */          ipmx +=1;
         /*  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 */          sw += weight[i];
       } else if(mle==2){          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
         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 */          /* 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])); */
       } else if(mle==3){  /* exponential inter-extrapolation */        } /* end of wave */
         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 */      } /* end of individual */
       } else if (mle==4){  /* mle=4 no inter-extrapolation */    }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */
         lli=log(out[s1][s2]); /* Original formula */      for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       } else{  /* mle=0 back to 1 */        for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
         lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */        for(mi=1; mi<= wav[i]-1; mi++){
         /*lli=log(out[s1][s2]); */ /* Original formula */          for (ii=1;ii<=nlstate+ndeath;ii++)
       }            for (j=1;j<=nlstate+ndeath;j++){
       /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/              oldm[ii][j]=(ii==j ? 1.0 : 0.0);
       /*if(lli ==000.0)*/              savm[ii][j]=(ii==j ? 1.0 : 0.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;          for(d=0; d<dh[mi][i]; d++){
       sw += weight[i];            newm=savm;
       ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;            agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
       /* if (lli < log(mytinydouble)){ */            cov[2]=agexact;
       /*   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); */            if(nagesqr==1)
       /*   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]); */              cov[3]= agexact*agexact;
       /* } */            for (kk=1; kk<=cptcovage;kk++) {
     } /* end of wave */              if(!FixedV[Tvar[Tage[kk]]])
   } /* end of individual */                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];    for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */    /* 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 */    l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
Line 6796  double funcone( double *x) Line 4646  double funcone( double *x)
   int ioffset=0;    int ioffset=0;
   int ipos=0,iposold=0,ncovv=0;    int ipos=0,iposold=0,ncovv=0;
   
   char labficresilk[NCOVMAX+1];  
   double cotvarv, cotvarvold;    double cotvarv, cotvarvold;
   double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];    double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
   double **out;    double **out;
Line 6829  double funcone( double *x) Line 4678  double funcone( double *x)
     /* Fixed */      /* Fixed */
     /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */      /* 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 *\/ */      /* 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 */      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("Debug3 TvarFind[%d]=%d",kf, TvarFind[kf]); */
       /* printf(" Tvar[TvarFind[kf]]=%d", Tvar[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]); */        /* 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)*/        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[ioffset+TvarFind[1]]=covar[Tvar[TvarFind[1]]][i];  */
 /*    cov[2+6]=covar[Tvar[6]][i];  */  /*    cov[2+6]=covar[Tvar[6]][i];  */
 /*    cov[2+6]=covar[2][i]; V2  */  /*    cov[2+6]=covar[2][i]; V2  */
Line 6848  double funcone( double *x) Line 4693  double funcone( double *x)
 /*    cov[2+9]=covar[Tvar[9]][i];  */  /*    cov[2+9]=covar[Tvar[9]][i];  */
 /*    cov[2+9]=covar[1][i]; V1  */  /*    cov[2+9]=covar[1][i]; V1  */
     }      }
     /* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4]         /* 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            is 5, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2]=6 
        has been calculated etc */           has been calculated etc */
     /* For an individual i, wav[i] gives the number of effective waves */        /* For an individual i, wav[i] gives the number of effective waves */
     /* We compute the contribution to Likelihood of each effective transition        /* We compute the contribution to Likelihood of each effective transition
        mw[mi][i] is real wave of the mi th effectve wave */           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];        /* 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];           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]           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 */      /* 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 *\/ */      /* 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?)*\/ */      /*   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?)*\/ */
Line 6986  double funcone( double *x) Line 4831  double funcone( double *x)
         *   Not a fixed cotvar[mw][itv][i]     6       7      6  2      7, 2,     6, 3,     7, 3,     6, 4,     7, 4}          *   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]           *   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 */        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  */          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*/          ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
Line 7001  double funcone( double *x) Line 4846  double funcone( double *x)
           cotvarv=covar[itv][i];  /* Good: In V6*V3, 3 is fixed at position of the data */            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); */            /* 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 */          if(ipos!=iposold){ /* Not a product or first of a product */
           cotvarvold=cotvarv;            cotvarvold=cotvarv;
         }else{ /* A second product */          }else{ /* A second product */
           cotvarv=cotvarv*cotvarvold;            cotvarv=cotvarv*cotvarvold;
         /* if(globpr) */  
         /*   sprintf(labficresilk+strlen(labficresilk)," *"); */  
           /* printf("DEBUG * \n"); */  
         }          }
         iposold=ipos;          iposold=ipos;
         cov[ioffset+ipos]=cotvarv;          cov[ioffset+ipos]=cotvarv;
Line 7052  double funcone( double *x) Line 4892  double funcone( double *x)
         cov[2]=agexact;          cov[2]=agexact;
         if(nagesqr==1)          if(nagesqr==1)
           cov[3]= agexact*agexact;            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(ncovva=1, iposold=0; ncovva <= ncovta ; ncovva++){ /* Time varying  covariates with age including individual from products, product is computed dynamically */
         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  */            itv=TvarAVVA[ncovva]; /*  TvarVV={3, 1, 3} gives the name of each varying covariate, exploding product Vn*Vm into Vn and then Vm  */
           ipos=Tage[kk];            ipos=TvarAVVAind[ncovva]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/
           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*/  
           /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */            /* 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 */            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); */              /* printf("DEBUG  ncovva=%d, Varying TvarAVVA[ncovva]=%d\n", ncovva, TvarAVVA[ncovva]); */
             cotvarv=cotvar[mw[mi][i]][TvarVVA[ncovva]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */               cotvarv=cotvar[mw[mi][i]][TvarAVVA[ncovva]][i];  /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */ 
           }else{ /* fixed covariate */            }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 *\/ */              /* 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 */              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 */            if(ipos!=iposold){ /* Not a product or first of a product */
Line 7088  double funcone( double *x) Line 4911  double funcone( double *x)
             cotvarv=cotvarv*cotvarvold;              cotvarv=cotvarv*cotvarvold;
           }            }
           iposold=ipos;            iposold=ipos;
             /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */
           cov[ioffset+ipos]=cotvarv*agexact;            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 */            /* For products */
         }          }
   
         /* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */          /* 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); */          /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
         out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                      1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                       1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
Line 7160  double funcone( double *x) Line 4976  double funcone( double *x)
  %11.6f %11.6f %11.6f ", \   %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,                  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]));                  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 ", \ */          /* %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, */          /*              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])); */          /*              2*weight[i]*lli,(s2==-1? -1: out[s1][s2]),(s2==-1? -1: savm[s1][s2])); */
Line 7171  double funcone( double *x) Line 4986  double funcone( double *x)
           /* printf(" %10.6f",-ll[k]*gipmx/gsw); */            /* printf(" %10.6f",-ll[k]*gipmx/gsw); */
         }          }
         fprintf(ficresilk," %10.6f ", -llt);          fprintf(ficresilk," %10.6f ", -llt);
         /* fprintf(ficresilk,"%s", labficresilk); */  
         /* printf(" %10.6f\n", -llt); */          /* printf(" %10.6f\n", -llt); */
         /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */          /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */
         /* fprintf(ficresilk,"%09ld ", num[i]); */ /* not necessary */          /* fprintf(ficresilk,"%09ld ", num[i]); */ /* not necessary */
Line 7189  double funcone( double *x) Line 5003  double funcone( double *x)
           }            }
           iposold=ipos;            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++) { */          /* for (kk=1; kk<=cptcovage;kk++) { */
         /*   if(!FixedV[Tvar[Tage[kk]]]){ */          /*   if(!FixedV[Tvar[Tage[kk]]]){ */
         /*     fprintf(ficresilk," %g*age",covar[Tvar[Tage[kk]]][i]); */          /*     fprintf(ficresilk," %g*age",covar[Tvar[Tage[kk]]][i]); */
Line 7218  double funcone( double *x) Line 5012  double funcone( double *x)
         /*     /\* printf(" %g*age",cotvar[mw[mi][i]][Tvar[Tage[kk]]][i]);/\\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\\/  *\/ */          /*     /\* 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"); */          /* printf("\n"); */
         /* } /\*  End debugILK *\/ */          /* } /\*  End debugILK *\/ */
         fprintf(ficresilk,"\n");          fprintf(ficresilk,"\n");
         /* printf("\n"); */  
       } /* End if globpr */        } /* End if globpr */
     } /* end of wave */      } /* end of wave */
   } /* end of individual */    } /* end of individual */
Line 7245  void likelione(FILE *ficres,double p[], Line 5063  void likelione(FILE *ficres,double p[],
      Plotting could be done.       Plotting could be done.
   */    */
   void pstamp(FILE *ficres);    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 */    if(*globpri !=0){ /* Just counts and sums, no printings */
     strcpy(fileresilk,"ILK_");       strcpy(fileresilk,"ILK_"); 
Line 7263  void likelione(FILE *ficres,double p[], Line 5081  void likelione(FILE *ficres,double p[],
     fprintf(ficresilk," -2*gipw/gsw*weight*ll(total) ");      fprintf(ficresilk," -2*gipw/gsw*weight*ll(total) ");
   
     /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */      /* if(debugILK){ /\* debugILK is set by a #d in a comment line *\/ */
     for(kf=1;kf <= ncovf; kf++){  /* Fixed covariates */        for(kf=1;kf <= ncovf; kf++){
       fprintf(ficresilk,"V%d",Tvar[TvarFind[kf]]);          fprintf(ficresilk,"V%d",Tvar[TvarFind[kf]]);
       /* printf("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 */        for(ncovv=1, iposold=0; ncovv <= ncovvt ; ncovv++){
         kvarold=kvar;          ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
         /* printf(" %d",ipos); */          if(ipos!=iposold){ /* Not a product or first of a product */
         fprintf(ficresilk," V%d",kvarold);            /* printf(" %d",ipos); */
       }else{ /* a second product */            fprintf(ficresilk," V%d",TvarVV[ncovv]);
         /* printf("*"); */          }else{
         fprintf(ficresilk,"*");            /* printf("*"); */
         fprintf(ficresilk," V%d",kvar);            fprintf(ficresilk,"*");
       }          }
       iposold=ipos;          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  */        for (kk=1; kk<=cptcovage;kk++) {
       ipos=Tage[kk];          if(!FixedV[Tvar[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]]); */
         /* printf(" %d*age(Fixed)",Tvar[Tage[kk]]); */            fprintf(ficresilk," %d*age(Fixed)",Tvar[Tage[kk]]);
         fprintf(ficresilk," %d*age(Fixed)",Tvar[Tage[kk]]);          }else{
       }else{            fprintf(ficresilk," %d*age(Varying)",Tvar[Tage[kk]]);/* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) */ 
         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) *\/  */
         /* printf(" %d*age(Varying)",Tvar[Tage[kk]]);/\* because cotvar starts now at first ncovcol+nqv+ (1 to nqtv) *\/  */          }
       }        }
     }  
     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 *\/ */  
       }  
       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"); */  
       }  
       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 *\/ */      /* } /\* End if debugILK *\/ */
     /* printf("\n"); */      /* printf("\n"); */
     fprintf(ficresilk,"\n");      fprintf(ficresilk,"\n");
Line 7349  void likelione(FILE *ficres,double p[], Line 5129  void likelione(FILE *ficres,double p[],
       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 \        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);  <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 */        for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */
         /* kvar=Tvar[TvarFind[kf]]; */ /* variable */           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> \           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]]);
 <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]]);           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 */        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 */          ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate */
Line 10107  void varprob(char optionfilefiname[], do Line 7888  void varprob(char optionfilefiname[], do
    double ***varpij;     double ***varpij;
   
    strcpy(fileresprob,"PROB_");      strcpy(fileresprob,"PROB_"); 
    strcat(fileresprob,fileres);     strcat(fileresprob,fileresu);
    if((ficresprob=fopen(fileresprob,"w"))==NULL) {     if((ficresprob=fopen(fileresprob,"w"))==NULL) {
      printf("Problem with resultfile: %s\n", fileresprob);       printf("Problem with resultfile: %s\n", fileresprob);
      fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);       fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
Line 10926  void printinggnuplot(char fileresu[], ch Line 8707  void printinggnuplot(char fileresu[], ch
   for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */    for(kf=1; kf <= ncovf; kf++){ /* For each simple dummy covariate of the model */
     kvar=Tvar[TvarFind[kf]]; /* variable name */      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+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 ++) {      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,"\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,"unset log;\n# For each simple dummy covariate of the model \n plot  \"%s\"",subdirf(fileresilk));
Line 11612  set ter svg size 640, 480\nunset log y\n Line 9395  set ter svg size 640, 480\nunset log y\n
             fprintf(ficgp," u %d:(",ioffset);               fprintf(ficgp," u %d:(",ioffset); 
             kl=0;              kl=0;
             strcpy(gplotcondition,"(");              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= 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,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(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 */                /* 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]][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++;                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++;                kl++;
               if(k <cptcoveff && cptcoveff>1)                if(k <cptcovs && cptcovs>1)
                 sprintf(gplotcondition+strlen(gplotcondition)," && ");                  sprintf(gplotcondition+strlen(gplotcondition)," && ");
             }              }
             strcpy(gplotcondition+strlen(gplotcondition),")");              strcpy(gplotcondition+strlen(gplotcondition),")");
Line 11707  set ter svg size 640, 480\nunset log y\n Line 9494  set ter svg size 640, 480\nunset log y\n
           }else{            }else{
             fprintf(ficgp,",\\\n '' ");              fprintf(ficgp,",\\\n '' ");
           }            }
           if(cptcoveff ==0){ /* No covariate */            /* if(cptcoveff ==0){ /\* No covariate *\/ */
             if(cptcovs ==0){ /* No covariate */
             ioffset=2; /* Age is in 2 */              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*/              /*# 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  16  17  18 */
Line 11731  set ter svg size 640, 480\nunset log y\n Line 9519  set ter svg size 640, 480\nunset log y\n
             /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */              /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
             iyearc=ioffset-1;              iyearc=ioffset-1;
             iagec=ioffset;              iagec=ioffset;
             fprintf(ficgp," u %d:(",ioffset); /* PROBLEM HERE VERIFY */              fprintf(ficgp," u %d:(",ioffset); 
             kl=0;              kl=0;
             strcpy(gplotcondition,"(");              strcpy(gplotcondition,"(");
             for (k=1; k<=cptcovs; k++){    /* For each covariate k of the resultline, get corresponding value lv for combination k1 */              for (k=1; k<=cptcovs; k++){    /* For each covariate k of the resultline, get corresponding value lv for combination k1 */
Line 11819  set ter svg size 640, 480\nunset log y\n Line 9607  set ter svg size 640, 480\nunset log y\n
     fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n");      fprintf(ficgp,"#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,"#model=1+age+%s \n",model);
     fprintf(ficgp,"# Type of graphic ng=%d\n",ng);      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(k1=1; k1 <=m; k1++)  /\* For each combination of covariate *\/ */
     for(nres=1; nres <= nresult; nres++){ /* For each resultline */      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
      /* k1=nres; */       /* k1=nres; */
Line 11948  set ter svg size 640, 480\nunset log y\n Line 9737  set ter svg size 640, 480\nunset log y\n
                 if(cptcovdageprod >0){                  if(cptcovdageprod >0){
                   /* if(j==Tprod[ijp]) { */ /* not necessary */                     /* if(j==Tprod[ijp]) { */ /* not necessary */ 
                     /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */                      /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
                     if(ijp <=cptcovprod) { /* Product */                      if(ijp <=cptcovprod) { /* Product Vn*Vm and age*VN*Vm*/
                       if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */                        if(DummyV[Tvardk[ijp][1]]==0){/* Vn is dummy */
                         if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are 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",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]]);                            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 */                          }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",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){                          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 */                          }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++;                        ijp++;
Line 12060  set ter svg size 640, 480\nunset log y\n Line 9849  set ter svg size 640, 480\nunset log y\n
                     /* if(j==Tprod[ijp]) { /\* *\/  */                      /* if(j==Tprod[ijp]) { /\* *\/  */
                       /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */                        /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
                       if(ijp <=cptcovprod) { /* Product */                        if(ijp <=cptcovprod) { /* Product */
                         if(DummyV[Tvard[ijp][1]]==0){/* Vn is dummy */                          if(DummyV[Tvardk[ijp][1]]==0){/* Vn is dummy */
                           if(DummyV[Tvard[ijp][2]]==0){/* Vn and Vm are 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",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]]); */                              /* 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 */                            }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",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]]); */                              /* 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 */                          }else{ /* Vn*Vm Vn is quanti */
                           if(DummyV[Tvard[ijp][2]]==0){                            if(DummyV[Tvardk[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]]);                              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]]); */                              /* 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 */                            }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]]); */                              /* fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */
                           }                             } 
                         }                          }
Line 12437  void prevforecast(char fileres[], double Line 10226  void prevforecast(char fileres[], double
   /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */    /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */
   /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */    /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */
   /* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */    /* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */
   i1=pow(2,cptcoveff);    /* i1=pow(2,cptcoveff); */
   if (cptcovn < 1){i1=1;}    /* 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,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); 
       
   fprintf(ficresf,"#****** Routine prevforecast **\n");    fprintf(ficresf,"#****** Routine prevforecast **\n");
       
 /*            if (h==(int)(YEARM*yearp)){ */  /*            if (h==(int)(YEARM*yearp)){ */
   for(nres=1; nres <= nresult; nres++) /* For each resultline */    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) */      k=TKresult[nres];
     if(i1 != 1 && TKresult[nres]!= k)      if(TKresult[nres]==0) k=1; /* To be checked for noresult */
       continue;      /*  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(invalidvarcomb[k]){      /* if(i1 != 1 && TKresult[nres]!= k) */
       printf("\nCombination (%d) projection ignored because no cases \n",k);       /*   continue; */
       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#");      fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#");
     for(j=1;j<=cptcoveff;j++) {      for(j=1;j<=cptcovs;j++){
       /* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); */        /* for(j=1;j<=cptcoveff;j++) { */
       fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[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 (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");      fprintf(ficresf," yearproj age");
     for(j=1; j<=nlstate+ndeath;j++){       for(j=1; j<=nlstate+ndeath;j++){ 
       for(i=1; i<=nlstate;i++)                for(i=1; i<=nlstate;i++)        
Line 12485  void prevforecast(char fileres[], double Line 10280  void prevforecast(char fileres[], double
           }            }
         }          }
         fprintf(ficresf,"\n");          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,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);          fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm);
                   
         for(j=1; j<=nlstate+ndeath;j++) {          for(j=1; j<=nlstate+ndeath;j++) {
Line 12579  void prevforecast(char fileres[], double Line 10376  void prevforecast(char fileres[], double
   /* if(jintmean==0) jintmean=1; */    /* if(jintmean==0) jintmean=1; */
   /* if(mintmean==0) jintmean=1; */    /* if(mintmean==0) jintmean=1; */
       
   i1=pow(2,cptcoveff);    /* i1=pow(2,cptcoveff); */
   if (cptcovn < 1){i1=1;}    /* 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);    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);    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");    fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
       
   for(nres=1; nres <= nresult; nres++) /* For each resultline */    for(nres=1; nres <= nresult; nres++){ /* For each resultline */
   for(k=1; k<=i1;k++){      k=TKresult[nres];
     if(i1 != 1 && TKresult[nres]!= k)      if(TKresult[nres]==0) k=1; /* To be checked for noresult */
       continue;    /* for(k=1; k<=i1;k++){ */
     if(invalidvarcomb[k]){    /*   if(i1 != 1 && TKresult[nres]!= k) */
       printf("\nCombination (%d) projection ignored because no cases \n",k);     /*     continue; */
       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#");      fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#");
     for(j=1;j<=cptcoveff;j++) {      for(j=1;j<=cptcovs;j++){
       fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);      /* 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]);        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");      fprintf(ficresfb," yearbproj age");
     for(j=1; j<=nlstate+ndeath;j++){      for(j=1; j<=nlstate+ndeath;j++){
       for(i=1; i<=nlstate;i++)        for(i=1; i<=nlstate;i++)
Line 12632  void prevforecast(char fileres[], double Line 10435  void prevforecast(char fileres[], double
           }            }
         }          }
         fprintf(ficresfb,"\n");          fprintf(ficresfb,"\n");
         for(j=1;j<=cptcoveff;j++)          /* for(j=1;j<=cptcoveff;j++) */
           fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[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);          fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm);
         for(i=1; i<=nlstate+ndeath;i++) {          for(i=1; i<=nlstate+ndeath;i++) {
           ppij=0.;ppi=0.;            ppij=0.;ppi=0.;
Line 13597  int decoderesult( char resultline[], int Line 11402  int decoderesult( char resultline[], int
   if (strlen(resultsav) >1){    if (strlen(resultsav) >1){
     j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */      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 */      TKresult[nres]=0; /* Combination for the nresult and the model */
     return (0);      return (0);
   }    }
   if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */    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, 1+age+%s.\n",j, cptcovs, model);fflush(ficlog);
     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);      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);
     /* return 1;*/      if(j==0)
         return 1;
   }    }
   for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */    for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */
     if(nbocc(resultsav,'=') >1){      if(nbocc(resultsav,'=') >1){
Line 13787  int decoderesult( char resultline[], int Line 11593  int decoderesult( char resultline[], int
       precov[nres][k1]=Tvalsel[k3q];        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]); */        /* 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++;;        k4q++;;
     }else if( Dummy[k1]==2 ){ /* For dummy with age product */      }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 */        /* Tvar[k1]; */ /* Age variable */ /* 17 age*V6*V2 ?*/
       /* Wrong we want the value of variable name Tvar[k1] */        /* Wrong we want the value of variable name Tvar[k1] */
               if(Typevar[k1]==2 || Typevar[k1]==3 ){ /* For product quant or dummy (with or without age) */
       k3= resultmodel[nres][k1]; /* nres=1 k1=2 resultmodel[2(V4)] = 1=k3 ; k1=3 resultmodel[3(V3)] = 2=k3*/          precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];      
       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)*/        /* 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]]); */
       TinvDoQresult[nres][(int)Tvarsel[k3]]=Tvalsel[k3]; /* TinvDoQresult[nres][4]=1 */        }else{
       precov[nres][k1]=Tvalsel[k3];          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]); */        /* 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 */      }else if( Dummy[k1]==3 ){ /* For quant with age product */
       k3q= resultmodel[nres][k1]; /* resultmodel[1(V5)] = 25.1=k3q */        if(Typevar[k1]==2 || Typevar[k1]==3 ){ /* For product quant or dummy (with or without age) */
       k2q=(int)Tvarsel[k3q]; /*  Tvarsel[resultmodel[1]]= Tvarsel[1] = 4=k2 */          precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];      
       TinvDoQresult[nres][(int)Tvarsel[k3q]]=Tvalsel[k3q]; /* TinvDoQresult[nres][5]=25.1 */        /* 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]]); */
       precov[nres][k1]=Tvalsel[k3q];        }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]); */        /* 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) */      }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]];              precov[nres][k1]=TinvDoQresult[nres][Tvardk[k1][1]] * TinvDoQresult[nres][Tvardk[k1][2]];      
Line 13883  int decodemodel( char model[], int lasto Line 11698  int decodemodel( char model[], int lasto
     if (strlen(modelsav) >1){ /* V2 +V3 +V4 +V6 +V7 +V6*V2 +V7*V2 +V6*V3 +V7*V3 +V6*V4 +V7*V4 +age*V2 +age*V3 +age*V4 +age*V6 +age*V7 +age*V6*V2 +V7*V2 +age*V6*V3 +age*V7*V3 +age*V6*V4 +age*V7*V4 */      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 '+' */        j=nbocc(modelsav,'+'); /**< j=Number of '+' */
       j1=nbocc(modelsav,'*'); /**< j1=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        cptcovt= j+1; /* Number of total covariates in the model, not including
                      * cst, age and age*age                        * cst, age and age*age 
                      * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/                       * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/
       /* including age products which are counted in cptcovage.        /* including age products which are counted in cptcovage.
        * but the covariates which are products must be treated          * but the covariates which are products must be treated 
        * separately: ncovn=4- 2=2 (V1+V3). */         * separately: ncovn=4- 2=2 (V1+V3). */
       cptcovprod=0; /**< Number of products and single product with age  V1*V2 +v3*age = 2 */        cptcovprod=0; /**< Number of products  V1*V2 +v3*age = 2 */
       cptcovdageprod=0; /* Number of double products with age age*Vn*VM or Vn*age*Vm or Vn*Vm*age */        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  */        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");*/        /* cptcovprodage=nboccstr(modelsav,"age");*/
               
       /*   Design        /*   Design
Line 13948  int decodemodel( char model[], int lasto Line 11763  int decodemodel( char model[], int lasto
         Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;          Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0;
       }        }
       cptcovage=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 */        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          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" */                                           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" */
Line 13956  int decodemodel( char model[], int lasto Line 11787  int decodemodel( char model[], int lasto
         /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/          /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
         /*scanf("%d",i);*/          /*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 */          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   */            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 OR (strb=age*V6*V2 or V6*V2*age or V6*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 */              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) */              if(strstr(strc,"age")!=0) { /* It means that strc=V2*age or age*V2 and thus that strd=Vn */
               cutl(stre,strf,strc,'*') ; /* if strc=age*Vm then stre=Vm and strf=age, if strc=Vm*age then stre=age and strf=Vm. */                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 ,  strd=Vn */                strcpy(strc,strb); /* save strb(=age*Vn*Vm) into strc */
               /* We want strb=Vn*Vm */                /* 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) . */                if(strcmp(strf,"age")==0){ /* strf is "age" so that stre=Vm =V2 . */
                 strcpy(strb,strd); /* strd=Vn  */                   strcpy(strb,strd);
                 strcat(strb,"*");                   strcat(strb,"*");
                 strcat(strb,stre);/* strb=Vn*Vm */                  strcat(strb,stre);
               }else{  /* strf=Vm so stre=age. strd=Vn  If strf=V6 then stre=V2 */                }else{  /* strf=Vm  If strf=V6 then stre=V2 */
                 strcpy(strb,strf);                  strcpy(strb,strf);
                 strcat(strb,"*");                  strcat(strb,"*");
                 strcat(strb,strd); /* strb=Vm*Vn */                  strcat(strb,stre);
                 strcpy(strd,strb); /* in order for strd to not be "age"  for next test (will be strd=Vn*Vm */                  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]]]);                /* 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 */                /* FixedV[Tvar[Tage[k]]]=0; /\* HERY not sure if V7*V4*age Fixed might not exist  yet*\/ */
             }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  */              }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(stre,strb); /* save full b in stre */
               strcpy(strb,strc); /* save short c in new short b for next block strb=Vn*Vm*/                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 */                strcpy(strf,strc); /* save short c in new short f */
Line 14007  int decodemodel( char model[], int lasto Line 11838  int decodemodel( char model[], int lasto
               Tvardk[k][1] =m; /* m 1 for V1*/                Tvardk[k][1] =m; /* m 1 for V1*/
               Tvard[k1][2] =n; /* n 4 for V4*/                Tvard[k1][2] =n; /* n 4 for V4*/
               Tvardk[k][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 */                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 */                  for (i=1; i<=lastobs;i++){/* For fixed product */
                   /* Computes the new covariate which is a product of                    /* Computes the new covariate which is a product of
Line 14043  int decodemodel( char model[], int lasto Line 11874  int decodemodel( char model[], int lasto
                 /*Tage[cptcovage]=k;*/ /* For age*V3*V2 Tage[1]=V3*V3=9 HERY too*/                  /*Tage[cptcovage]=k;*/ /* For age*V3*V2 Tage[1]=V3*V3=9 HERY too*/
                 /* Tvar[Tage[cptcovage]]=k1; */                  /* Tvar[Tage[cptcovage]]=k1; */
                 cptcovprodvage++;                  cptcovprodvage++;
                 k12=2*k11-1;  
                 FixedV[ncovcolt+k12]=1; /* We expand Vn*Vm */                  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 *\/ */              /* Tage[cptcovage]=k;  /\*  V2+V1+V4+V3*age Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 *\/ */
Line 14157  int decodemodel( char model[], int lasto Line 11988  int decodemodel( char model[], int lasto
                                 /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);                                  /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
                                   scanf("%d",i);*/                                    scanf("%d",i);*/
       } /* end of loop + on total covariates */        } /* end of loop + on total covariates */
   
         
     } /* end if strlen(modelsave == 0) age*age might exist */      } /* end if strlen(modelsave == 0) age*age might exist */
   } /* end if strlen(model == 0) */    } /* 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  */    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  */
Line 14195  Fixed[k] 0=fixed (product or simple), 1 Line 12028  Fixed[k] 0=fixed (product or simple), 1
 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);  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++){ Fixed[k]=0; Dummy[k]=0;}
   for(k=1;k<=NCOVMAX; k++){TvarFind[k]=0; TvarVind[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 */    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 */      if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
       Fixed[k]= 0;        Fixed[k]= 0;
Line 15448  int main(int argc, char *argv[]) Line 13286  int main(int argc, char *argv[])
   /* double ***mobaverage; */    /* double ***mobaverage; */
   double wald;    double wald;
   
   char line[MAXLINE];    char line[MAXLINE], linetmp[MAXLINE];
   char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];    char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE];
   
   char  modeltemp[MAXLINE];    char  modeltemp[MAXLINE];
Line 15550  int main(int argc, char *argv[]) Line 13388  int main(int argc, char *argv[])
   getcwd(pathcd, size);    getcwd(pathcd, size);
 #endif  #endif
   syscompilerinfo(0);    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){    if(argc <=1){
     printf("\nEnter the parameter file name: ");      printf("\nEnter the parameter file name: ");
     if(!fgets(pathr,FILENAMELENGTH,stdin)){      if(!fgets(pathr,FILENAMELENGTH,stdin)){
Line 15781  int main(int argc, char *argv[]) Line 13619  int main(int argc, char *argv[])
     }else      }else
       break;        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){      if (num_filled != 1){
       printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);        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);        fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line);
Line 16154  Please run with mle=-1 to get a correct Line 14003  Please run with mle=-1 to get a correct
   Tvard=imatrix(1,NCOVMAX,1,2); /* n=Tvard[k1][1]  and m=Tvard[k1][2] gives the couple n,m of the k1 th product Vn*Vm    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.                               * 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) */                              * 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    Tage=ivector(1,NCOVMAX); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age
                          4 covariates (3 plus signs)                           4 covariates (3 plus signs)
                          Tage[1=V3*age]= 4; Tage[2=age*V4] = 3                           Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
Line 17259  Please run with mle=-1 to get a correct Line 15108  Please run with mle=-1 to get a correct
         date2dmy(datebackf,&jbackf, &mbackf, &anbackf);          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, \      printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
                  model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,prevbcast, estepm, \                   model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,prevbcast, estepm, \
Line 17409  Please run with mle=-1 to get a correct Line 15258  Please run with mle=-1 to get a correct
   
     pstamp(ficreseij);      pstamp(ficreseij);
                                   
     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */      /* i1=pow(2,cptcoveff); /\* Number of combination of dummy covariates *\/ */
     if (cptcovn < 1){i1=1;}      /* if (cptcovn < 1){i1=1;} */
           
     for(nres=1; nres <= nresult; nres++) /* For each resultline */      for(nres=1; nres <= nresult; nres++){ /* For each resultline */
     for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */      /* for(k=1; k<=i1;k++){ /\* For any combination of dummy covariates, fixed and varying *\/ */
       if(i1 != 1 && TKresult[nres]!= k)        /* if(i1 != 1 && TKresult[nres]!= k) */
         continue;        /*        continue; */
       fprintf(ficreseij,"\n#****** ");        fprintf(ficreseij,"\n#****** ");
       printf("\n#****** ");        printf("\n#****** ");
       for(j=1;j<=cptcoveff;j++) {        for(j=1;j<=cptcovs;j++){
         fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]);        /* for(j=1;j<=cptcoveff;j++) { */
         printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[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 */        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) */          printf(" V%d=%lg ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); /* TvarsQ[j] gives the name of the jth quantitative (fixed or time v) */
Line 17489  Please run with mle=-1 to get a correct Line 15341  Please run with mle=-1 to get a correct
       /* */        /* */
       if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */        if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
         continue;          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(ficrest,"\n# model %s \n#****** Result for:", model);
       fprintf(ficlog,"\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 */        /* It might not be a good idea to mix dummies and quantitative */
Line 17664  Please run with mle=-1 to get a correct Line 15516  Please run with mle=-1 to get a correct
   
           
     free_vector(weight,firstobs,lastobs);      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(Tvard,1,NCOVMAX,1,2);
     free_imatrix(s,1,maxwav+1,firstobs,lastobs);      free_imatrix(s,1,maxwav+1,firstobs,lastobs);
     free_matrix(anint,1,maxwav,firstobs,lastobs);       free_matrix(anint,1,maxwav,firstobs,lastobs); 
Removed from v.1.1  
changed lines
  Added in v.1.2

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>