--- imach/src/imach.c 2024/05/12 20:29:32 1.361 +++ imach/src/imach.c 2024/06/28 13:53:38 1.365 @@ -1,6 +1,20 @@ -/* $Id: imach.c,v 1.361 2024/05/12 20:29:32 brouard Exp $ +/* $Id: imach.c,v 1.365 2024/06/28 13:53:38 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + 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 Summary: Version 0.99s5 @@ -1240,7 +1254,7 @@ Important routines - 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. - 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. @@ -1405,12 +1419,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.365 2024/06/28 13:53:38 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.365 $ $Date: 2024/06/28 13:53:38 $"; 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 +1545,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 @@ -3007,24 +3027,31 @@ static void print2() /* print a line of /* printf("... after %u function calls ...\n", nf); */ /* printf("... including %u linear searches ...\n", nl); */ /* printf("%10d %10d%14.7g",nl, nf, fx); */ - printf ( "\n" ); + /* printf ( "\n" ); */ printf ( " Linear searches %d", nl ); + fprintf (ficlog, " 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 ); + fprintf (ficlog, " Function evaluations %d", nf ); + fprintf (ficlog, " Function value FX = %.12lf\n", fx ); #ifdef DEBUGPRAX printf("n=%d prin=%d\n",n,prin); #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 ) { /* 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:" ); */ } printf("\n"); + fprintf(ficlog,"\n"); } @@ -3231,7 +3258,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 +3734,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 +4245,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 +6620,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=4; + /* 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 +6633,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 +8629,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; */ @@ -9969,7 +9996,9 @@ void printinggnuplot(char fileresu[], ch char dirfileres[256],optfileres[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 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 vpopbased; int ioffset; /* variable offset for columns */ @@ -10767,8 +10796,9 @@ 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]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */ 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,lv, kl+1, vlv ); + sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%lg " ,kl,lv, kl+1, vlv );/* Solved but quantitative must be shifted */ kl++; if(k 1) sprintf(gplotcondition+strlen(gplotcondition)," && "); @@ -10782,7 +10812,9 @@ 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, \ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt ); 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, \ iyearc, iagec, offyear, \ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc ); @@ -12283,7 +12315,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,14 +15876,22 @@ 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 free_ivector(flatdir,1,npar); #endif /* LINMINORIGINAL*/ - +#endif /* POWELL */ hesscov(matcov, hess, p, NDIM, delti, 1e-4, gompertz); for(i=1; i <=NDIM; i++) @@ -15984,7 +16025,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 +16096,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]);