--- imach/src/imach.c 2010/06/14 07:50:17 1.139 +++ imach/src/imach.c 2014/02/10 22:17:31 1.144 @@ -1,6 +1,26 @@ -/* $Id: imach.c,v 1.139 2010/06/14 07:50:17 brouard Exp $ +/* $Id: imach.c,v 1.144 2014/02/10 22:17:31 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.144 2014/02/10 22:17:31 brouard + *** empty log message *** + + Revision 1.143 2014/01/26 09:45:38 brouard + Summary: Version 0.98nR (to be improved, but gives same optimization results as 0.98k. Nice, promising + + * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... + (Module): Version 0.98nR Running ok, but output format still only works for three covariates. + + Revision 1.142 2014/01/26 03:57:36 brouard + Summary: gnuplot changed plot w l 1 has to be changed to plot w l lt 2 + + * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... + + Revision 1.141 2014/01/26 02:42:01 brouard + * imach.c (Module): Trying to merge old staffs together while being at Tokyo. Not tested... + + Revision 1.140 2011/09/02 10:37:54 brouard + Summary: times.h is ok with mingw32 now. + Revision 1.139 2010/06/14 07:50:17 brouard After the theft of my laptop, I probably lost some lines of codes which were not uploaded to the CVS tree. I remember having already fixed agemin agemax which are pointers now but not cvs saved. @@ -388,9 +408,12 @@ #include extern int errno; -/* #include */ +#ifdef LINUX #include #include "timeval.h" +#else +#include +#endif #ifdef GSL #include @@ -400,7 +423,7 @@ extern int errno; /* #include */ /* #define _(String) gettext (String) */ -#define MAXLINE 256 +#define MAXLINE 1024 /* Was 256. Overflow with 312 with 2 states and 4 covariates. Should be ok */ #define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ @@ -409,18 +432,18 @@ extern int errno; #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ -#define MAXPARM 128 /* Maximum number of parameters for the optimization */ -#define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */ +#define MAXPARM 128 /**< Maximum number of parameters for the optimization */ +#define NPARMAX 64 /**< (nlstate+ndeath-1)*nlstate*ncovmodel */ #define NINTERVMAX 8 -#define NLSTATEMAX 8 /* Maximum number of live states (for func) */ -#define NDEATHMAX 8 /* Maximum number of dead states (for func) */ -#define NCOVMAX 20 /* Maximum number of covariates */ +#define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ +#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ +#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ #define MAXN 20000 -#define YEARM 12. /* Number of months per year */ +#define YEARM 12. /**< Number of months per year */ #define AGESUP 130 #define AGEBASE 40 -#define AGEGOMP 10. /* Minimal age for Gompertz adjustment */ +#define AGEGOMP 10. /**< Minimal age for Gompertz adjustment */ #ifdef UNIX #define DIRSEPARATOR '/' #define CHARSEPARATOR "/" @@ -431,11 +454,11 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.139 2010/06/14 07:50:17 brouard Exp $ */ +/* $Id: imach.c,v 1.144 2014/02/10 22:17:31 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98m, April 2010, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.139 $ $Date: 2010/06/14 07:50:17 $"; +char version[]="Imach version 0.98nR2, January 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)"; +char fullversion[]="$Revision: 1.144 $ $Date: 2014/02/10 22:17:31 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -558,12 +581,18 @@ double dateintmean=0; double *weight; int **s; /* Status */ -double *agedc, **covar, idx; -int **nbcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; +double *agedc; +double **covar; /**< covar[i,j], value of jth covariate for individual i, + * covar=matrix(0,NCOVMAX,1,n); + * cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; */ +double idx; +int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ +int **codtab; /**< codtab=imatrix(1,100,1,10); */ +int **Tvard, *Tprod, cptcovprod, *Tvaraff; double *lsurv, *lpop, *tpop; -double ftol=FTOL; /* Tolerance for computing Max Likelihood */ -double ftolhess; /* Tolerance for computing hessian */ +double ftol=FTOL; /**< Tolerance for computing Max Likelihood */ +double ftolhess; /**< Tolerance for computing hessian */ /**************** split *************************/ static int split( char *path, char *dirc, char *name, char *ext, char *finame ) @@ -1090,7 +1119,7 @@ char *asc_diff_time(long time_sec, char sec_left = (sec_left) %(60*60); minutes = (sec_left) /60; sec_left = (sec_left) % (60); - sprintf(ascdiff,"%d day(s) %d hour(s) %d minute(s) %d second(s)",days, hours, minutes, sec_left); + sprintf(ascdiff,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left); return ascdiff; } @@ -1280,7 +1309,7 @@ double **prevalim(double **prlim, int nl /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ /*printf("ij=%d cov[3]=%lf \n",ij, cov[3]);*/ - out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); + out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm); /* Bug Valgrind */ savm=oldm; oldm=newm; @@ -1796,7 +1825,7 @@ double funcone( double *x) 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]); */ if(globpr){ - fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ + fprintf(ficresilk,"%9ld %6d %2d %2d %1d %1d %3d %11.6f %8.4f\ %11.6f %11.6f %11.6f ", \ num[i],i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i], 2*weight[i]*lli,out[s1][s2],savm[s1][s2]); @@ -2183,8 +2212,8 @@ void freqsummary(char fileres[], int ia first=1; - for(k1=1; k1<=j;k1++){ - for(i1=1; i1<=ncodemax[k1];i1++){ + for(k1=1; k1<=j ; k1++){ /* Loop on covariates */ + for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */ j1++; /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ @@ -2192,20 +2221,26 @@ void freqsummary(char fileres[], int ia for (jk=-5; jk<=nlstate+ndeath; jk++) for(m=iagemin; m <= iagemax+3; m++) freq[i][jk][m]=0; - - for (i=1; i<=nlstate; i++) - for(m=iagemin; m <= iagemax+3; m++) - prop[i][m]=0; + + for (i=1; i<=nlstate; i++) + for(m=iagemin; m <= iagemax+3; m++) + prop[i][m]=0; dateintsum=0; k2cpt=0; for (i=1; i<=imx; i++) { bool=1; - if (cptcovn>0) { - for (z1=1; z1<=cptcoveff; z1++) - if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) - bool=0; + if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ + for (z1=1; z1<=cptcoveff; z1++) + if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]){ + bool=0; + printf("bool=%d i=%d, z1=%d, i1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", + bool,i,z1, i1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1], + j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1); + /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/ + } } + if (bool==1){ for(m=firstpass; m<=lastpass; m++){ k2=anint[m][i]+(mint[m][i]/12.); @@ -2233,6 +2268,9 @@ void freqsummary(char fileres[], int ia fprintf(ficresp, "\n#********** Variable "); for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); fprintf(ficresp, "**********\n#"); + fprintf(ficlog, "\n#********** Variable "); + for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); + fprintf(ficlog, "**********\n#"); } for(i=1; i<=nlstate;i++) fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i); @@ -2538,21 +2576,24 @@ void concatwav(int wav[], int **dh, int } jmean=sum/k; printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean); - fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean); + fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %d) Max=%d (%d) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean); } /*********** Tricode ****************************/ void tricode(int *Tvar, int **nbcode, int imx) { - /* Uses cptcovn+2*cptcovprod as the number of covariates */ - /* Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */ + /**< Uses cptcovn+2*cptcovprod as the number of covariates */ + /* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 + /* Boring subroutine which should only output nbcode[Tvar[j]][k] + /* nbcode[Tvar[j][1]= + */ int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX; int modmaxcovj=0; /* Modality max of covariates j */ cptcoveff=0; - for (k=0; k File (multiple files are possible if covariates are present): %s\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); fprintf(fichtm,"\n
Probability is computed over estepm=%d months.

\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); /* fprintf(fichtm,"\n
Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year

\n", stepm,YEARM,digitp,digit); @@ -3807,17 +3848,17 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); + fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 1,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); + fprintf(ficgp,"\" t\"95\%% CI\" w l lt 2,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - fprintf(ficgp,"\" t\"\" w l 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1)); + fprintf(ficgp,"\" t\"\" w l lt 2,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 3",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1)); } } /*2 eme*/ @@ -3840,14 +3881,14 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - fprintf(ficgp,"\" t\"\" w l 0,"); + fprintf(ficgp,"\" t\"\" w l lt 1,"); fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1); for (j=1; j<= nlstate+1 ; j ++) { if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); else fprintf(ficgp," \%%*lf (\%%*lf)"); } - if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l 0"); - else fprintf(ficgp,"\" t\"\" w l 0,"); + if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 1"); + else fprintf(ficgp,"\" t\"\" w l lt 1,"); } } @@ -3932,9 +3973,9 @@ plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1 fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); else fprintf(ficgp," exp(p%d+p%d*x",i,i+1); - ij=1; + ij=1;/* To be checked else nbcode[0][0] wrong */ for(j=3; j <=ncovmodel; j++) { - if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { + if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */ fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); ij++; } @@ -4604,8 +4645,8 @@ int readdata(char datafile[], int firsto lval=strtol(strb,&endptr,10); /* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/ if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); - fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog); + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog); return 1; } } @@ -4619,8 +4660,8 @@ int readdata(char datafile[], int firsto month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); - fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog); + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog); return 1; } anint[j][i]= (double) year; @@ -4635,8 +4676,8 @@ int readdata(char datafile[], int firsto month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); - fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); return 1; } andc[i]=(double) year; @@ -4650,13 +4691,13 @@ int readdata(char datafile[], int firsto month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line); - fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); return 1; } if (year==9999) { - printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line); - fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog); + printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line); + fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog); return 1; } @@ -4668,8 +4709,8 @@ int readdata(char datafile[], int firsto errno=0; dval=strtod(strb,&endptr); if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); - fprintf(ficlog,"Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); + printf("Error reading data around '%f' at line number %d, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); + fprintf(ficlog,"Error reading data around '%f' at line number %d, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei); fflush(ficlog); return 1; } @@ -4684,13 +4725,13 @@ int readdata(char datafile[], int firsto errno=0; lval=strtol(strb,&endptr,10); if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line); - fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog); + printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line); + fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog); return 1; } } if(lval <-1 || lval >1){ - printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \ + printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ For example, for multinomial values like 1, 2 and 3,\n \ @@ -4699,7 +4740,7 @@ int readdata(char datafile[], int firsto and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \ output of IMaCh is often meaningless.\n \ Exiting.\n",lval,linei, i,line,j); - fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\n \ + fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \ For example, for multinomial values like 1, 2 and 3,\n \ @@ -4761,6 +4802,11 @@ int decodemodel ( char model[], int last fprintf(ficlog,"Error. AGE must be in lower case model=%s ",model);fflush(ficlog); return 1; } + if (strstr(model,"v") !=0){ + printf("Error. 'v' must be in upper case 'V' model=%s ",model); + fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog); + return 1; + } /* This loop fills the array Tvar from the string 'model'.*/ /* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */ @@ -5157,7 +5203,7 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); numlinepar++; - puts(line); + fputs(line,stdout); fputs(line,ficparo); fputs(line,ficlog); } @@ -5173,7 +5219,8 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); numlinepar++; - puts(line); + fputs(line, stdout); + //puts(line); fputs(line,ficparo); fputs(line,ficlog); } @@ -5226,7 +5273,7 @@ int main(int argc, char *argv[]) ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); numlinepar++; - puts(line); + fputs(line,stdout); fputs(line,ficparo); fputs(line,ficlog); } @@ -5276,7 +5323,7 @@ run imach with mle=-1 to get a correct t ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); numlinepar++; - puts(line); + fputs(line,stdout); fputs(line,ficparo); fputs(line,ficlog); } @@ -5317,7 +5364,7 @@ run imach with mle=-1 to get a correct t ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); numlinepar++; - puts(line); + fputs(line,stdout); fputs(line,ficparo); fputs(line,ficlog); } @@ -5385,7 +5432,7 @@ run imach with mle=-1 to get a correct t anint=matrix(1,maxwav,1,n); s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ tab=ivector(1,NCOVMAX); - ncodemax=ivector(1,8); /* hard coded ? */ + ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ /* Reads data from file datafile */ if (readdata(datafile, firstobs, lastobs, &imx)==1) @@ -5408,12 +5455,14 @@ run imach with mle=-1 to get a correct t ncovcol + k1 If already ncovcol=4 and model=V2+V1+V1*V4+age*V3 Tvar[3=V1*V4]=4+1 etc */ - Tprod=ivector(1,15); + Tprod=ivector(1,15); /* Gives the position of a product */ /* Tprod[k1=1]=3(=V1*V4) for V2+V1+V1*V4+age*V3 if V2+V1+V1*V4+age*V3+V3*V2 TProd[k1=2]=5 (V3*V2) */ Tvaraff=ivector(1,15); - Tvard=imatrix(1,15,1,2); /* For V3*V2 Tvard[k1=2][1]=3 (V3) Tvard[k1=2][2]=2(V2) */ + Tvard=imatrix(1,15,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. + * Tvard[k1=2][1]=3 (V3) Tvard[k1=2][2]=2(V2) */ Tage=ivector(1,15); /* Gives the covariate id of covariates associated with age: V2 + V1 + age*V4 + V3*age 4 covariates (3 plus signs) Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 @@ -5462,21 +5511,40 @@ run imach with mle=-1 to get a correct t ncodemax[1]=1; if (cptcovn > 0) tricode(Tvar,nbcode,imx); - codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of - the estimations*/ + codtab=imatrix(1,100,1,10); /**< codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + */ h=0; m=pow(2,cptcoveff); for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */ - for(i=1; i <=(m/pow(2,k));i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ - for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */ - for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ + for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ + for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/ + for(cpt=1; cpt <=pow(2,k-1); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ h++; - if (h>m) { + if (h>m) h=1; - codtab[h][k]=j; - codtab[h][Tvar[k]]=j; - } + /**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 + * h 1 2 3 4 + *______________________________ + * 1 i=1 1 i=1 1 i=1 1 i=1 1 + * 2 2 1 1 1 + * 3 i=2 1 2 1 1 + * 4 2 2 1 1 + * 5 i=3 1 i=2 1 2 1 + * 6 2 1 2 1 + * 7 i=4 1 2 2 1 + * 8 2 2 2 1 + * 9 i=5 1 i=3 1 i=2 1 1 + * 10 2 1 1 1 + * 11 i=6 1 2 1 1 + * 12 2 2 1 1 + * 13 i=7 1 i=4 1 2 1 + * 14 2 1 2 1 + * 15 i=8 1 2 2 1 + * 16 2 2 2 1 + */ + codtab[h][k]=j; + codtab[h][Tvar[k]]=j; printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); } } @@ -5504,7 +5572,8 @@ run imach with mle=-1 to get a correct t else{ fprintf(ficgp,"\n# %s\n", version); fprintf(ficgp,"# %s\n", optionfilegnuplot); - fprintf(ficgp,"set missing 'NaNq'\n"); + //fprintf(ficgp,"set missing 'NaNq'\n"); + fprintf(ficgp,"set datafile missing 'NaNq'\n"); } /* fclose(ficgp);*/ /*--------- index.htm --------*/ @@ -5943,7 +6012,7 @@ Interval (in months) between two waves: while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); - puts(line); + fputs(line,stdout); fputs(line,ficparo); } ungetc(c,ficpar); @@ -5963,7 +6032,7 @@ Interval (in months) between two waves: while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); - puts(line); + fputs(line,stdout); fputs(line,ficparo); } ungetc(c,ficpar); @@ -5977,7 +6046,7 @@ Interval (in months) between two waves: while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); - puts(line); + fputs(line,stdout); fputs(line,ficparo); } ungetc(c,ficpar); @@ -5993,7 +6062,7 @@ Interval (in months) between two waves: while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); - puts(line); + fputs(line,stdout); fputs(line,ficparo); } ungetc(c,ficpar); @@ -6053,14 +6122,14 @@ Interval (in months) between two waves: agebase=ageminpar; agelim=agemaxpar; ftolpl=1.e-10; - i1=cptcoveff; + i1=pow(2,cptcoveff); if (cptcovn < 1){i1=1;} for(cptcov=1,k=0;cptcov<=i1;cptcov++){ - for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ + //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ k=k+1; /* to clean */ - printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,codtab[cptcod][cptcov],nbcode); + //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]); fprintf(ficrespl,"\n#******"); printf("\n#******"); fprintf(ficlog,"\n#******"); @@ -6081,9 +6150,9 @@ Interval (in months) between two waves: for(i=1; i<=nlstate;i++) fprintf(ficrespl," %.5f", prlim[i][i]); fprintf(ficrespl,"\n"); - } - } - } + } /* Age */ + /* was end of cptcod */ + } /* cptcov */ fclose(ficrespl); /*------------- h Pij x at various ages ------------*/ @@ -6370,10 +6439,9 @@ Interval (in months) between two waves: /*---------- End : free ----------------*/ if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); - } /* mle==-3 arrives here for freeing */ endfree: - free_matrix(prlim,1,nlstate,1,nlstate); + free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */ free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); @@ -6414,9 +6482,9 @@ Interval (in months) between two waves: fprintf(ficlog,"Local time at start %s\nLocal time at end %s\n",strstart, strtend); printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); - printf("Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); + printf("Total time was %ld Sec.\n", end_time.tv_sec -start_time.tv_sec); fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); - fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); + fprintf(ficlog,"Total time was %ld Sec.\n", end_time.tv_sec -start_time.tv_sec); /* printf("Total time was %d uSec.\n", total_usecs);*/ /* if(fileappend(fichtm,optionfilehtm)){ */ fprintf(fichtm,"
Local time at start %s
Local time at end %s
\n",strstart, strtend);