--- imach/src/imach.c 2024/05/12 20:29:32 1.361
+++ imach/src/imach.c 2024/06/28 08:00:31 1.362
@@ -1,6 +1,11 @@
-/* $Id: imach.c,v 1.361 2024/05/12 20:29:32 brouard Exp $
+/* $Id: imach.c,v 1.362 2024/06/28 08:00:31 brouard Exp $
$State: Exp $
$Log: imach.c,v $
+ 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
Summary: Version 0.99s5
@@ -1405,12 +1410,12 @@ double gnuplotversion=GNUPLOTVERSION;
#define ODIRSEPARATOR '\\'
#endif
-/* $Id: imach.c,v 1.361 2024/05/12 20:29:32 brouard Exp $ */
+/* $Id: imach.c,v 1.362 2024/06/28 08:00:31 brouard Exp $ */
/* $State: Exp $ */
#include "version.h"
char version[]=__IMACH_VERSION__;
char copyright[]="April 2024,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-2024";
-char fullversion[]="$Revision: 1.361 $ $Date: 2024/05/12 20:29:32 $";
+char fullversion[]="$Revision: 1.362 $ $Date: 2024/06/28 08:00:31 $";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
@@ -1531,6 +1536,12 @@ char *endptr;
long lval;
double dval;
+/* This for praxis gegen */
+ /* int prin=1; */
+ double h0=0.25;
+ double macheps;
+ double ffmin;
+
#define NR_END 1
#define FREE_ARG char*
#define FTOL 1.0e-10
@@ -3231,7 +3242,7 @@ L1: /* L1 or try loop */
if (k > 0) *d2 = 0;
}
#ifdef DEBUGPRAX
- printf(" bebe end of min x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
+ printf(" bebe end of min x1 might be very wrong x1=%14.8e fx=%14.8e d2=%14.8e\n",*x1, fx, *d2);
#endif
if (*d2 <= small_windows) *d2 = small_windows;
*x1 = x2; fx = fm;
@@ -3707,7 +3718,7 @@ mloop:
printf("praxis4 macheps=%14g h=%14g step=%14g small=%14g t=%14g\n",macheps,h, h0,small_windows, t);
#endif
/* min(0, 2, &d[0], &s, fx, 0); /\* mac heps not global *\/ */
- minny(1, 2, &d[1], &s, fx, 0); /* mac heps not global */
+ 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
printf("praxis5 macheps=%14g h=%14g looks at sign of s=%14g fx=%14g\n",macheps,h, s,fx);
#endif
@@ -4218,7 +4229,7 @@ void powell(double p[], double **xi, int
printf(" + 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) {
printf(" + V%d ",Tvar[j]);
fprintf(ficlog," + V%d ",Tvar[j]);
@@ -6593,10 +6604,10 @@ void mlikeli(FILE *ficres,double p[], in
#else /* FLATSUP */
/* powell(p,xi,npar,ftol,&iter,&fret,func);*/
/* praxis ( t0, h0, n, prin, x, beale_f ); */
- int prin=1;
- double h0=0.25;
- double macheps;
- double fmin;
+ /* int prin=1; */
+ /* double h0=0.25; */
+ /* double macheps; */
+ /* double fmin; */
macheps=pow(16.0,-13.0);
/* #include "praxis.h" */
/* Be careful that praxis start at x[0] and powell start at p[1] */
@@ -6606,7 +6617,7 @@ printf("Praxis Gegenfurtner \n");
fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog);
/* praxis ( ftol, h0, npar, prin, p1, 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");
#endif /* FLATSUP */
@@ -8602,7 +8613,7 @@ void concatwav(int wav[], int **dh, int
double ***gradg, ***trgradg; /**< for var eij */
double **gradgp, **trgradgp; /**< 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 age,agelim, hf;
/* double ***mobaverage; */
@@ -12283,7 +12294,8 @@ double gompertz(double x[])
A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
} else if (cens[i] == 0){
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
printf("Gompertz cens[%d] neither 1 nor 0\n",i);
/*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */
@@ -15843,8 +15855,16 @@ Interval (in months) between two waves:
flatdir=ivector(1,npar);
for (j=1;j<=npar;j++) flatdir[j]=0;
#endif /*LINMINORIGINAL */
- powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
-#endif
+ /* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */
+ /* 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);
#ifdef LINMINORIGINAL
#else
@@ -15984,7 +16004,7 @@ Please run with mle=-1 to get a correct
fprintf(ficlog," + age*age ");
fprintf(fichtm, "| + age*age | ");
}
- for(j=1;j <=ncovmodel-2;j++){
+ for(j=1;j <=ncovmodel-2-nagesqr;j++){
if(Typevar[j]==0) {
printf(" + V%d ",Tvar[j]);
fprintf(ficres," + V%d ",Tvar[j]);
@@ -16055,7 +16075,7 @@ Please run with mle=-1 to get a correct
fprintf(ficlog," + age*age ");
fprintf(fichtm, "+ age*age | ");
}
- for(j=1;j <=ncovmodel-2;j++){
+ for(j=1;j <=ncovmodel-2-nagesqr;j++){
if(Typevar[j]==0) {
printf(" + V%d ",Tvar[j]);
fprintf(fichtm, "+ V%d | ",Tvar[j]);