--- imach/src/imach.c 2024/06/28 08:00:31 1.362 +++ imach/src/imach.c 2024/07/08 14:26:18 1.367 @@ -1,6 +1,24 @@ -/* $Id: imach.c,v 1.362 2024/06/28 08:00:31 brouard Exp $ +/* $Id: imach.c,v 1.367 2024/07/08 14:26:18 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.367 2024/07/08 14:26:18 brouard + Summary: 0.99s7 + + * imach.c (Module): Some bug fixes: in drawings when age*age is + included in the model as well as with quantitative variables. + + 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 @@ -1245,7 +1263,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. @@ -1410,12 +1428,12 @@ double gnuplotversion=GNUPLOTVERSION; #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.362 2024/06/28 08:00:31 brouard Exp $ */ +/* $Id: imach.c,v 1.367 2024/07/08 14:26:18 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.362 $ $Date: 2024/06/28 08:00:31 $"; +char fullversion[]="$Revision: 1.367 $ $Date: 2024/07/08 14:26:18 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1474,7 +1492,8 @@ int countcallfunc=0; /* Count the numbe int selected(int kvar); /* Is covariate kvar selected for printing results */ 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 **oldms, **newms, **savms; /* Fixed working pointers to matrices */ double **ddnewms, **ddoldms, **ddsavms; /* for freeing later */ @@ -1521,14 +1540,16 @@ char optionfilegnuplot[FILENAMELENGTH], /* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */ /* struct timezone tzp; */ /* 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; 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; */ +struct tm tm, tml; char strcurr[80], strfor[80]; @@ -1992,9 +2013,7 @@ int nboccstr(char *textin, char *chain) /* in="+V7*V4+age*V2+age*V3+age*V4" chain="age" */ char *strloc; - int i,j=0; - - i=0; + int j=0; strloc=textin; /* strloc points to "^+V7*V4+age+..." in textin */ for(;;) { @@ -2128,14 +2147,15 @@ int **imatrix(long nrl, long nrh, long n } /****************** free_imatrix *************************/ -void free_imatrix(m,nrl,nrh,ncl,nch) - int **m; - long nch,ncl,nrh,nrl; - /* free an int matrix allocated by imatrix() */ -{ - free((FREE_ARG) (m[nrl]+ncl-NR_END)); - free((FREE_ARG) (m+nrl-NR_END)); -} +/* void free_imatrix(m,nrl,nrh,ncl,nch); */ +/* int **m; */ +/* long nch,ncl,nrh,nrl; */ +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)); +} /******************* matrix *******************************/ double **matrix(long nrl, long nrh, long ncl, long nch) @@ -2395,8 +2415,8 @@ values at the three points, fa, fb , and double ulim,u,r,q, dum; double fu; - double scale=10.; - int iterscale=0; + /* double scale=10.; */ + /* int iterscale=0; */ *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]) */ @@ -2786,7 +2806,8 @@ double *e; /* used in minfit, don't konw static int prin; /* added */ static int n; static double *x; -static double (*fun)(); +static double (*fun)(double *x); /* New for clang */ +/* static double (*fun)(); */ /* static double (*fun)(double *x, int n); */ /* these will be set by praxis to the global control parameters */ @@ -3011,31 +3032,38 @@ 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() /* print a line of traces */ { - int i; double fmin=0.; + int i; /* double fmin=0.; */ /* printf("\n"); */ /* printf("... chi square reduced to ... %20.10e\n", fx); */ /* 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"); } @@ -3783,7 +3811,7 @@ next: #ifdef NR_SHIFT fx = (*fun)((x-1), n); #else - fx = (*fun)(x, n); + fx = (*fun)(x); #endif /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/ nf++; @@ -4577,8 +4605,10 @@ void powell(double p[], double **xi, int int i, ii,j,k, k1; double *min, *max, *meandiff, maxmax,sumnew=0.; - /* double **matprod2(); */ /* test */ - double **out, cov[NCOVMAX+1], **pmij(); /* **pmmij is a global variable feeded with oldms etc */ + double **matprod2(double **out, double **in,int nrl, int nrh, int ncl, int nch, int ncolol, int ncoloh, double **b); /* test */ /* for clang */ +/* 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 agefin, delaymax=200. ; /* 100 Max number of years to converge */ int ncvloop=0; @@ -4784,7 +4814,8 @@ void powell(double p[], double **xi, int int first=0; double *min, *max, *meandiff, maxmax,sumnew=0.; /* 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 **dnewm, **doldm, **dsavm; /* for use */ double **oldm, **savm; /* for use */ @@ -5049,7 +5080,8 @@ double **pmij(double **ps, double *cov, */ 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 agefin; double k3=0.; /* constant of the w_x diagonal matrix (in order for B to sum to 1 even for death state) */ @@ -5444,7 +5476,8 @@ double ***hbxij(double ***po, int nhstep */ 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 agexact; /*double agebegin, ageend;*/ @@ -6604,7 +6637,7 @@ 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; */ + int prin=4; /* double h0=0.25; */ /* double macheps; */ /* double fmin; */ @@ -9980,7 +10013,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 */ @@ -10174,7 +10209,7 @@ void printinggnuplot(char fileresu[], ch fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres); /* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */ fprintf(ficgp,"set title \"Alive state %d %s model=1+age+%s\" font \"Helvetica,12\"\n",cpt,gplotlabel,model); - fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); + fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] [0:1] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres); /* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */ /* k1-1 error should be nres-1*/ for (i=1; i<= nlstate ; i ++) { @@ -10739,12 +10774,15 @@ set ter svg size 640, 480\nunset log y\n }else{ fprintf(ficgp,",\\\n '' "); } - if(cptcoveff ==0){ /* No covariate */ + /* if(cptcoveff ==0){ /\* No covariate *\/ */ + if(cptcovs ==0){ /* No covariate */ ioffset=2; /* Age is in 2 */ /*# yearproj age p11 p21 p31 p.1 p12 p22 p32 p.2 p13 p23 p33 p.3 p14 p24 p34 p.4*/ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 */ /*# V1 = 1 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 */ + /*# V1 = 1 yearproj age 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 */ fprintf(ficgp," u %d:(", ioffset); if(i==nlstate+1){ fprintf(ficgp," $%d/(1.-$%d)):1 t 'pw.%d' with line lc variable ", \ @@ -10758,9 +10796,15 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \ ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt ); }else{ /* more than 2 covariates */ - ioffset=2*cptcoveff+2; /* Age is in 4 or 6 or etc.*/ + /* ioffset=2*cptcoveff+2; */ /* Age is in 4 or 6 or etc.*/ + ioffset=2*cptcovs+2; /* Age is in 4 or 6 or etc.*/ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ + /* # Forecasting at date 3/1/2003 */ + /* V1=0 V2=1 V3=0 V6=2.47 yearproj age */ + /* # 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 */ + /* # p11 p21 p31 wp.1 p12 p22 p32 wp.2 p13 p23 p33 wp.3 p14 p24 p34 wp.4 */ + /* 1 0 2 1 3 0 6 2.47 2003 100 1.000 0.000 0.000 0.297 0.000 1.000 0.000 0.207 0.000 0.000 1.000 0.497 0.000 0.000 0.000 0.000 */ iyearc=ioffset-1; iagec=ioffset; fprintf(ficgp," u %d:(",ioffset); @@ -10778,8 +10822,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)," && "); @@ -10793,7 +10838,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 ); @@ -10874,27 +10921,38 @@ 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 16 17 18 */ fprintf(ficgp," u %d:(", ioffset); if(i==nlstate+1){ - fprintf(ficgp," $%d/(1.-$%d)):1 t 'bw%d' with line lc variable ", \ - ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); + fprintf(ficgp," $%d):1 t 'bw%d' with line lc variable ", \ + ioffset+(cpt-1)*(nlstate+1)+1+(i-1),cpt ); + /* fprintf(ficgp," $%d/(1.-$%d)):1 t 'bw%d' with line lc variable ", \ */ + /* ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt ); */ fprintf(ficgp,",\\\n '' "); fprintf(ficgp," u %d:(",ioffset); fprintf(ficgp," (($1-$2) == %d ) ? $%d : 1/0):1 with labels center not ", \ offbyear, \ ioffset+(cpt-1)*(nlstate+1)+1+(i-1) ); - }else - fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ", \ - ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt,i ); + }else /* not sure divided by 1- to be checked */ + fprintf(ficgp," $%d) t 'b%d%d' with line ", \ + ioffset+(cpt-1)*(nlstate+1)+1+(i-1),cpt,i ); + /* fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ", \ */ + /* ioffset+(cpt-1)*(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt,i ); */ }else{ /* more than 2 covariates */ - ioffset=2*cptcoveff+2; /* Age is in 4 or 6 or etc.*/ + /* ioffset=2*cptcoveff+2; /\* Age is in 4 or 6 or etc.*\/ */ + ioffset=2*cptcovs+2; /* Age is in 4 or 6 or etc.*/ /*# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/ /*# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 */ +/* #****** hbijx=probability over h years, hb.jx is weighted by observed prev */ +/* # V1=0 V2=1 V3=0 V6=2.47 */ +/* yearbproj age b11 b21 b31 b.1 b12 b22 b32 b.2 b13 b23 b33 b.3 b14 b24 b34 b.4 */ +/* # Back Forecasting at date 3/1/2003 */ +/* 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 */ +/* 1 0 2 1 3 0 6 2.47 2003 50 1.000 0.000 0.000 0.714 0.000 1.000 0.000 0.286 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 */ iyearc=ioffset-1; iagec=ioffset; fprintf(ficgp," u %d:(",ioffset); kl=0; strcpy(gplotcondition,"("); for (k=1; k<=cptcovs; k++){ /* For each covariate k of the resultline, get corresponding value lv for combination k1 */ - if(Dummy[modelresult[nres][k]]==0){ /* To be verified */ + /* if(Dummy[modelresult[nres][k]]==0){ /\* To be verified *\/ */ /* for (k=1; k<=cptcoveff; k++){ /\* For each covariate writing the chain of conditions *\/ */ /* lv= decodtabm(k1,k,cptcoveff); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */ /* lv= codtabm(k1,TnsdVar[Tvaraff[k]]); /\* Should be the covariate value corresponding to combination k1 and covariate k *\/ */ @@ -10911,7 +10969,7 @@ set ter svg size 640, 480\nunset log y\n kl++; if(k 1) sprintf(gplotcondition+strlen(gplotcondition)," && "); - } + /* } */ /* end dummy */ } strcpy(gplotcondition+strlen(gplotcondition),")"); /* kl=6+(cpt-1)*(nlstate+1)+1+(i-1); /\* 6+(1-1)*(2+1)+1+(1-1)=7, 6+(2-1)(2+1)+1+(1-1)=10 *\/ */ @@ -15870,7 +15928,7 @@ Interval (in months) between two waves: #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++)