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

Diff for /imach/src/imach.c between versions 1.361 and 1.366

version 1.361, 2024/05/12 20:29:32 version 1.366, 2024/07/02 09:42:10
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.366  2024/07/02 09:42:10  brouard
     Summary: trying clang on Linux
   
     Revision 1.365  2024/06/28 13:53:38  brouard
     * imach.c (Module): fixing some bugs in gnuplot and quantitative variables, but not completely solved
   
     Revision 1.364  2024/06/28 12:27:05  brouard
     * imach.c (Module): fixing some bugs in gnuplot and quantitative variables, but not completely solved
   
     Revision 1.363  2024/06/28 09:31:55  brouard
     Summary: Adding log lines too
   
     Revision 1.362  2024/06/28 08:00:31  brouard
     Summary: 0.99s6
   
     * imach.c (Module): s6 errors with age*age (harmless).
   
   Revision 1.361  2024/05/12 20:29:32  brouard    Revision 1.361  2024/05/12 20:29:32  brouard
   Summary: Version 0.99s5    Summary: Version 0.99s5
   
Line 1240  Important routines Line 1257  Important routines
 - Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities)  - Tricode which tests the modality of dummy variables (in order to warn with wrong or empty modalities)
   and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually.    and returns the number of efficient covariates cptcoveff and modalities nbcode[Tvar[k]][1]= 0 and nbcode[Tvar[k]][2]= 1 usually.
 - printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables  - printinghtml which outputs results like life expectancy in and from a state for a combination of modalities of dummy variables
   o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, eĢliminating 1 1 if    o There are 2**cptcoveff combinations of (0,1) for cptcoveff variables. Outputting only combinations with people, eliminating 1 1 if
     race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless.      race White (0 0), Black vs White (1 0), Hispanic (0 1) and 1 1 being meaningless.
   
   
Line 1469  int countcallfunc=0;  /* Count the numbe Line 1486  int countcallfunc=0;  /* Count the numbe
 int selected(int kvar); /* Is covariate kvar selected for printing results */  int selected(int kvar); /* Is covariate kvar selected for printing results */
   
 double jmean=1; /* Mean space between 2 waves */  double jmean=1; /* Mean space between 2 waves */
 double **matprod2(); /* test */  double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b); /* test */
   /* double **matprod2();  *//* test */
 double **oldm, **newm, **savm; /* Working pointers to matrices */  double **oldm, **newm, **savm; /* Working pointers to matrices */
 double **oldms, **newms, **savms; /* Fixed working pointers to matrices */  double **oldms, **newms, **savms; /* Fixed working pointers to matrices */
 double   **ddnewms, **ddoldms, **ddsavms; /* for freeing later */  double   **ddnewms, **ddoldms, **ddsavms; /* for freeing later */
Line 1516  char optionfilegnuplot[FILENAMELENGTH], Line 1534  char optionfilegnuplot[FILENAMELENGTH],
 /* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */  /* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */
 /* struct timezone tzp; */  /* struct timezone tzp; */
 /* extern int gettimeofday(); */  /* extern int gettimeofday(); */
 struct tm tml, *gmtime(), *localtime();  
   
 extern time_t time();  /* extern time_t time(); */ /* Commented out for clang */
   /* struct tm tml, *gmtime(), *localtime(); */
   
   
 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 */  time_t   rlast_btime; /* raw time */
 struct tm tm;  /* struct tm tm; */
   struct tm tm, tml;
   
 char strcurr[80], strfor[80];  char strcurr[80], strfor[80];
   
Line 1531  char *endptr; Line 1551  char *endptr;
 long lval;  long lval;
 double dval;  double dval;
   
   /* This for praxis gegen */
     /* int prin=1; */
     double h0=0.25;
     double macheps;
     double ffmin;
   
 #define NR_END 1  #define NR_END 1
 #define FREE_ARG char*  #define FREE_ARG char*
 #define FTOL 1.0e-10  #define FTOL 1.0e-10
Line 1981  int nboccstr(char *textin, char *chain) Line 2007  int nboccstr(char *textin, char *chain)
   /*  in="+V7*V4+age*V2+age*V3+age*V4"  chain="age" */    /*  in="+V7*V4+age*V2+age*V3+age*V4"  chain="age" */
   char *strloc;    char *strloc;
       
   int i,j=0;    int j=0;
   
   i=0;  
   
   strloc=textin; /* strloc points to "^+V7*V4+age+..." in textin */    strloc=textin; /* strloc points to "^+V7*V4+age+..." in textin */
   for(;;) {    for(;;) {
Line 2117  int **imatrix(long nrl, long nrh, long n Line 2141  int **imatrix(long nrl, long nrh, long n
 }   } 
   
 /****************** free_imatrix *************************/  /****************** free_imatrix *************************/
 void free_imatrix(m,nrl,nrh,ncl,nch)  /* void free_imatrix(m,nrl,nrh,ncl,nch); */
       int **m;  /*       int **m; */
       long nch,ncl,nrh,nrl;   /*       long nch,ncl,nrh,nrl; */
      /* free an int matrix allocated by imatrix() */   void free_imatrix(int **m,long nrl, long nrh, long ncl, long nch)
 {        /* free an int matrix allocated by imatrix() */
   free((FREE_ARG) (m[nrl]+ncl-NR_END));   {
   free((FREE_ARG) (m+nrl-NR_END));     free((FREE_ARG) (m[nrl]+ncl-NR_END));
 }     free((FREE_ARG) (m+nrl-NR_END));
   }
   
 /******************* matrix *******************************/  /******************* matrix *******************************/
 double **matrix(long nrl, long nrh, long ncl, long nch)  double **matrix(long nrl, long nrh, long ncl, long nch)
Line 2384  values at the three points, fa, fb , and Line 2409  values at the three points, fa, fb , and
   double ulim,u,r,q, dum;    double ulim,u,r,q, dum;
   double fu;     double fu; 
   
   double scale=10.;    /* double scale=10.; */
   int iterscale=0;    /* int iterscale=0; */
   
   *fa=(*func)(*ax); /*  xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/    *fa=(*func)(*ax); /*  xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/
   *fb=(*func)(*bx); /*  xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */    *fb=(*func)(*bx); /*  xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */
Line 2775  double *e; /* used in minfit, don't konw Line 2800  double *e; /* used in minfit, don't konw
 static int prin; /* added */  static int prin; /* added */
 static int n;  static int n;
 static double *x;  static double *x;
 static double (*fun)();  static double (*fun)(double *x); /* New for clang */
   /* static double (*fun)(); */
 /* static double (*fun)(double *x, int n); */  /* static double (*fun)(double *x, int n); */
   
 /* these will be set by praxis to the global control parameters */  /* these will be set by praxis to the global control parameters */
Line 3000  static void print()  /* print a line of Line 3026  static void print()  /* print a line of
 /* static void print2(int n, double *x, int prin, double fx, int nf, int nl) */ /* print a line of traces */  /* static void print2(int n, double *x, int prin, double fx, int nf, int nl) */ /* print a line of traces */
 static void print2() /* print a line of traces */  static void print2() /* print a line of traces */
 {  {
   int i; double fmin=0.;    int i; /* double fmin=0.; */
   
    /* printf("\n"); */     /* printf("\n"); */
    /* printf("... chi square reduced to ... %20.10e\n", fx); */     /* printf("... chi square reduced to ... %20.10e\n", fx); */
    /* printf("... after %u function calls ...\n", nf); */     /* printf("... after %u function calls ...\n", nf); */
    /* printf("... including %u linear searches ...\n", nl); */     /* printf("... including %u linear searches ...\n", nl); */
    /* printf("%10d    %10d%14.7g",nl, nf, fx); */     /* printf("%10d    %10d%14.7g",nl, nf, fx); */
   printf ( "\n" );    /* printf ( "\n" ); */
   printf ( "  Linear searches      %d", nl );    printf ( "  Linear searches      %d", nl );
     fprintf (ficlog, "  Linear searches      %d", nl );
   /* printf ( "  Linear searches      %d\n", nl ); */    /* printf ( "  Linear searches      %d\n", nl ); */
   /* printf ( "  Function evaluations %d\n", nf ); */    /* printf ( "  Function evaluations %d\n", nf ); */
   /* printf ( "  Function value FX = %g\n", fx ); */    /* printf ( "  Function value FX = %g\n", fx ); */
   printf ( "  Function evaluations %d", nf );    printf ( "  Function evaluations %d", nf );
   printf ( "  Function value FX = %.12lf\n", fx );    printf ( "  Function value FX = %.12lf\n", fx );
     fprintf (ficlog, "  Function evaluations %d", nf );
     fprintf (ficlog, "  Function value FX = %.12lf\n", fx );
 #ifdef DEBUGPRAX  #ifdef DEBUGPRAX
    printf("n=%d prin=%d\n",n,prin);     printf("n=%d prin=%d\n",n,prin);
 #endif  #endif
    if(fx <= fmin) printf(" UNDEFINED "); else  printf("%14.7g",log(fx-fmin));     /* if(fx <= fmin) printf(" UNDEFINED "); else  printf("%14.7g",log(fx-fmin)); */
    if ( n <= 4 || 2 < prin )     if ( n <= 4 || 2 < prin )
    {     {
      /* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */       /* for(i=1;i<=n;i++)printf("%14.7g",x[i-1]); */
      for(i=1;i<=n;i++)printf("%14.7g",x[i]);       for(i=1;i<=n;i++){
          printf(" %14.7g",x[i]);
          fprintf(ficlog," %14.7g",x[i]);
        }
      /* r8vec_print ( n, x, "  X:" ); */       /* r8vec_print ( n, x, "  X:" ); */
    }     }
    printf("\n");     printf("\n");
      fprintf(ficlog,"\n");
  }   }
   
   
Line 3231  L1:  /* L1 or try loop */ Line 3264  L1:  /* L1 or try loop */
       if (k > 0) *d2 = 0;        if (k > 0) *d2 = 0;
    }     }
 #ifdef DEBUGPRAX  #ifdef DEBUGPRAX
    printf(" bebe end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);     printf(" bebe end of min x1 might be very wrong x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
 #endif  #endif
    if (*d2 <= small_windows) *d2 = small_windows;     if (*d2 <= small_windows) *d2 = small_windows;
    *x1 = x2; fx = fm;     *x1 = x2; fx = fm;
Line 3707  mloop: Line 3740  mloop:
    printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);      printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t); 
 #endif  #endif
    /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */     /* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */
    minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global */     minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global it seems that fx doesn't correspond to f(s=*x1) */
 #ifdef DEBUGPRAX  #ifdef DEBUGPRAX
    printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx);      printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx); 
 #endif  #endif
Line 3772  next: Line 3805  next:
 #ifdef NR_SHIFT  #ifdef NR_SHIFT
           fx = (*fun)((x-1), n);            fx = (*fun)((x-1), n);
 #else  #else
           fx = (*fun)(x, n);            fx = (*fun)(x);
 #endif  #endif
           /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/            /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/
           nf++;            nf++;
Line 4218  void powell(double p[], double **xi, int Line 4251  void powell(double p[], double **xi, int
         printf("  + age*age  ");          printf("  + age*age  ");
         fprintf(ficlog,"  + age*age  ");          fprintf(ficlog,"  + age*age  ");
     }      }
     for(j=1;j <=ncovmodel-2;j++){      for(j=1;j <=ncovmodel-2-nagesqr;j++){
       if(Typevar[j]==0) {        if(Typevar[j]==0) {
         printf("  +      V%d  ",Tvar[j]);          printf("  +      V%d  ",Tvar[j]);
         fprintf(ficlog,"  +      V%d  ",Tvar[j]);          fprintf(ficlog,"  +      V%d  ",Tvar[j]);
Line 4566  void powell(double p[], double **xi, int Line 4599  void powell(double p[], double **xi, int
           
     int i, ii,j,k, k1;      int i, ii,j,k, k1;
   double *min, *max, *meandiff, maxmax,sumnew=0.;    double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b); /* test */ /* for clang */
   double **out, cov[NCOVMAX+1], **pmij(); /* **pmmij is a global variable feeded with oldms etc */  /* double **matprod2(); */ /* test */
     /* double **out, cov[NCOVMAX+1], **pmij(); */ /* **pmmij is a global variable feeded with oldms etc */
     double **out, cov[NCOVMAX+1], **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate); /* **pmmij is a global variable feeded with oldms etc */
   double **newm;    double **newm;
   double agefin, delaymax=200. ; /* 100 Max number of years to converge */    double agefin, delaymax=200. ; /* 100 Max number of years to converge */
   int ncvloop=0;    int ncvloop=0;
Line 4773  void powell(double p[], double **xi, int Line 4808  void powell(double p[], double **xi, int
   int first=0;    int first=0;
   double *min, *max, *meandiff, maxmax,sumnew=0.;    double *min, *max, *meandiff, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **bmij();    double **out, cov[NCOVMAX+1], **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate,  double ***prevacurrent, int ij);
     /* double **out, cov[NCOVMAX+1], **bmij(); */ /* Deprecated in clang */
   double **newm;    double **newm;
   double         **dnewm, **doldm, **dsavm;  /* for use */    double         **dnewm, **doldm, **dsavm;  /* for use */
   double         **oldm, **savm;  /* for use */    double         **oldm, **savm;  /* for use */
Line 5038  double **pmij(double **ps, double *cov, Line 5074  double **pmij(double **ps, double *cov,
    */     */
   int ii, j;    int ii, j;
       
   double  **pmij();    double  **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate);
     /* double  **pmij(); */ /* No more for clang */
   double sumnew=0.;    double sumnew=0.;
   double agefin;    double agefin;
   double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */    double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */
Line 5433  double ***hbxij(double ***po, int nhstep Line 5470  double ***hbxij(double ***po, int nhstep
   */    */
   
   int i, j, d, h, k1;    int i, j, d, h, k1;
   double **out, cov[NCOVMAX+1], **bmij();    double **out, cov[NCOVMAX+1], **bmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate,  double ***prevacurrent, int ij);
     /* double **out, cov[NCOVMAX+1], **bmij(); */ /* No more for clang */
   double **newm, ***newmm;    double **newm, ***newmm;
   double agexact;    double agexact;
   /*double agebegin, ageend;*/    /*double agebegin, ageend;*/
Line 6593  void mlikeli(FILE *ficres,double p[], in Line 6631  void mlikeli(FILE *ficres,double p[], in
 #else  /* FLATSUP */  #else  /* FLATSUP */
 /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/  /*  powell(p,xi,npar,ftol,&iter,&fret,func);*/
 /*   praxis ( t0, h0, n, prin, x, beale_f ); */  /*   praxis ( t0, h0, n, prin, x, beale_f ); */
   int prin=1;   int prin=4;
   double h0=0.25;    /* double h0=0.25; */
   double macheps;    /* double macheps; */
   double fmin;    /* double fmin; */
   macheps=pow(16.0,-13.0);    macheps=pow(16.0,-13.0);
 /* #include "praxis.h" */  /* #include "praxis.h" */
   /* Be careful that praxis start at x[0] and powell start at p[1] */    /* Be careful that praxis start at x[0] and powell start at p[1] */
Line 6606  printf("Praxis Gegenfurtner \n"); Line 6644  printf("Praxis Gegenfurtner \n");
 fprintf(ficlog, "Praxis  Gegenfurtner\n");fflush(ficlog);  fprintf(ficlog, "Praxis  Gegenfurtner\n");fflush(ficlog);
 /* praxis ( ftol, h0, npar, prin, p1, func ); */  /* praxis ( ftol, h0, npar, prin, p1, func ); */
   /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */    /* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */
   fmin = praxis(ftol,macheps, h0, npar, prin, p, func);    ffmin = praxis(ftol,macheps, h0, npar, prin, p, func);
 printf("End Praxis\n");  printf("End Praxis\n");
 #endif  /* FLATSUP */  #endif  /* FLATSUP */
   
Line 8602  void  concatwav(int wav[], int **dh, int Line 8640  void  concatwav(int wav[], int **dh, int
    double ***gradg, ***trgradg; /**< for var eij */     double ***gradg, ***trgradg; /**< for var eij */
    double **gradgp, **trgradgp; /**< for var p point j */     double **gradgp, **trgradgp; /**< for var p point j */
    double *gpp, *gmp; /**< for var p point j */     double *gpp, *gmp; /**< for var p point j */
    double **varppt; /**< for var e.. nlstate+1 to nlstate+ndeath */     double **varppt; /**< for var p.3 p.death nlstate+1 to nlstate+ndeath */
    double ***p3mat;     double ***p3mat;
    double age,agelim, hf;     double age,agelim, hf;
    /* double ***mobaverage; */     /* double ***mobaverage; */
Line 9969  void printinggnuplot(char fileresu[], ch Line 10007  void printinggnuplot(char fileresu[], ch
   char dirfileres[256],optfileres[256];    char dirfileres[256],optfileres[256];
   char gplotcondition[256], gplotlabel[256];    char gplotcondition[256], gplotlabel[256];
   int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0;    int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0;
   int lv=0, vlv=0, kl=0;    /* int lv=0, vlv=0, kl=0; */
     int lv=0, kl=0;
     double vlv=0;
   int ng=0;    int ng=0;
   int vpopbased;    int vpopbased;
   int ioffset; /* variable offset for columns */    int ioffset; /* variable offset for columns */
Line 10767  set ter svg size 640, 480\nunset log y\n Line 10807  set ter svg size 640, 480\nunset log y\n
               /* 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++;
                 /* Problem with quantitative variables TinvDoQresult[nres] */
               /* 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 );                sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%lg " ,kl,lv, kl+1, vlv );/* Solved but quantitative must be shifted */
               kl++;                kl++;
               if(k <cptcovs && cptcovs>1)                if(k <cptcovs && cptcovs>1)
                 sprintf(gplotcondition+strlen(gplotcondition)," && ");                  sprintf(gplotcondition+strlen(gplotcondition)," && ");
Line 10782  set ter svg size 640, 480\nunset log y\n Line 10823  set ter svg size 640, 480\nunset log y\n
               fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \                fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt );
               fprintf(ficgp,",\\\n '' ");                fprintf(ficgp,",\\\n '' ");
               fprintf(ficgp," u %d:(",iagec);                 fprintf(ficgp," u %d:(",iagec); /* Below iyearc should be increades if quantitative variable in the reult line */
                 /* $7==6 && $8==2.47 ) && (($9-$10) == 1953 ) ? $12/(1.-$24) : 1/0):7 with labels center not */
                 /* but was  && $7==6 && $8==2 ) && (($7-$8) == 1953 ) ? $12/(1.-$24) : 1/0):7 with labels center not */
               fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \                fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \
                       iyearc, iagec, offyear,                           \                        iyearc, iagec, offyear,                           \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc );
Line 12283  double gompertz(double x[]) Line 12326  double gompertz(double x[])
          A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));           A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
        } else if (cens[i] == 0){         } else if (cens[i] == 0){
         A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))          A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
           +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);            +log(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM);
           /* +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); */  /* To be seen */
       } else        } else
          printf("Gompertz cens[%d] neither 1 nor 0\n",i);           printf("Gompertz cens[%d] neither 1 nor 0\n",i);
       /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */        /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
Line 15843  Interval (in months) between two waves: Line 15887  Interval (in months) between two waves:
   flatdir=ivector(1,npar);     flatdir=ivector(1,npar); 
   for (j=1;j<=npar;j++) flatdir[j]=0;     for (j=1;j<=npar;j++) flatdir[j]=0; 
 #endif /*LINMINORIGINAL */  #endif /*LINMINORIGINAL */
      powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);      /* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */
 #endif      /* double h0=0.25; */
     macheps=pow(16.0,-13.0);
     printf("Praxis Gegenfurtner mle=%d\n",mle);
     fprintf(ficlog, "Praxis  Gegenfurtner mle=%d\n", mle);fflush(ficlog);
      /* ffmin = praxis(ftol,macheps, h0, npar, prin, p, gompertz); */
     /* For the Gompertz we use only two parameters */
     int _npar=2;
      ffmin = praxis(ftol,macheps, h0, _npar, 4, p, gompertz);
     printf("End Praxis\n");
     fclose(ficrespow);      fclose(ficrespow);
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
 #else  #else
       free_ivector(flatdir,1,npar);         free_ivector(flatdir,1,npar); 
 #endif  /* LINMINORIGINAL*/  #endif  /* LINMINORIGINAL*/
       #endif /* POWELL */   
     hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz);       hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz); 
   
     for(i=1; i <=NDIM; i++)      for(i=1; i <=NDIM; i++)
Line 15984  Please run with mle=-1 to get a correct Line 16036  Please run with mle=-1 to get a correct
       fprintf(ficlog,"  + age*age  ");        fprintf(ficlog,"  + age*age  ");
       fprintf(fichtm, "<th>+ age*age</th>");        fprintf(fichtm, "<th>+ age*age</th>");
     }      }
     for(j=1;j <=ncovmodel-2;j++){      for(j=1;j <=ncovmodel-2-nagesqr;j++){
       if(Typevar[j]==0) {        if(Typevar[j]==0) {
         printf("  +      V%d  ",Tvar[j]);          printf("  +      V%d  ",Tvar[j]);
         fprintf(ficres,"  +      V%d  ",Tvar[j]);          fprintf(ficres,"  +      V%d  ",Tvar[j]);
Line 16055  Please run with mle=-1 to get a correct Line 16107  Please run with mle=-1 to get a correct
         fprintf(ficlog,"  + age*age  ");          fprintf(ficlog,"  + age*age  ");
         fprintf(fichtm, "<th>+ age*age</th>");          fprintf(fichtm, "<th>+ age*age</th>");
       }        }
       for(j=1;j <=ncovmodel-2;j++){        for(j=1;j <=ncovmodel-2-nagesqr;j++){
         if(Typevar[j]==0) {          if(Typevar[j]==0) {
           printf("  +      V%d  ",Tvar[j]);            printf("  +      V%d  ",Tvar[j]);
           fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);            fprintf(fichtm, "<th>+ V%d</th>",Tvar[j]);
Removed from v.1.361  
changed lines
  Added in v.1.366

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