--- imach/src/imach.c 2010/04/29 18:11:38 1.137 +++ imach/src/imach.c 2014/01/26 02:42:01 1.141 @@ -1,6 +1,19 @@ -/* $Id: imach.c,v 1.137 2010/04/29 18:11:38 brouard Exp $ +/* $Id: imach.c,v 1.141 2014/01/26 02:42:01 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + 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. + + Revision 1.138 2010/04/30 18:19:40 brouard + *** empty log message *** + Revision 1.137 2010/04/29 18:11:38 brouard (Module): Checking covariates for more complex models than V1+V2. A lot of change to be done. Unstable. @@ -381,9 +394,12 @@ #include extern int errno; -/* #include */ +#ifdef LINUX #include #include "timeval.h" +#else +#include +#endif #ifdef GSL #include @@ -393,7 +409,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"*/ @@ -424,11 +440,11 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.137 2010/04/29 18:11:38 brouard Exp $ */ +/* $Id: imach.c,v 1.141 2014/01/26 02:42:01 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98m, April 2010, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.137 $ $Date: 2010/04/29 18:11:38 $"; +char version[]="Imach version 0.98n, January 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Sicentific Research 25293121)"; +char fullversion[]="$Revision: 1.141 $ $Date: 2014/01/26 02:42:01 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -551,8 +567,14 @@ 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 */ @@ -1083,7 +1105,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; } @@ -1260,21 +1282,21 @@ double **prevalim(double **prlim, int nl for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){ newm=savm; /* Covariates have to be included here again */ - cov[2]=agefin; - - for (k=1; k<=cptcovn;k++) { - cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; - /* printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/ - } - for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; - for (k=1; k<=cptcovprod;k++) - cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; - - /*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]);*/ + cov[2]=agefin; + + for (k=1; k<=cptcovn;k++) { + cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; + /* printf("ij=%d k=%d Tvar[k]=%d nbcode=%d cov=%lf codtab[ij][Tvar[k]]=%d \n",ij,k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], codtab[ij][Tvar[k]]);*/ + } + for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; + for (k=1; k<=cptcovprod;k++) + cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + + /*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); - + savm=oldm; oldm=newm; maxmax=0.; @@ -1301,41 +1323,56 @@ double **prevalim(double **prlim, int nl double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) { - double s1, s2; + /* According to parameters values stored in x and the covariate's values stored in cov, + computes the probability to be observed in state j being in state i by appying the + model to the ncovmodel covariates (including constant and age). + lnpijopii=ln(pij/pii)= aij+bij*age+cij*v1+dij*v2+... = sum_nc=1^ncovmodel xij(nc)*cov[nc] + and, according on how parameters are entered, the position of the coefficient xij(nc) of the + ncth covariate in the global vector x is given by the formula: + j=i nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel + Computes ln(pij/pii) (lnpijopii), deduces pij/pii by exponentiation, + sums on j different of i to get 1-pii/pii, deduces pii, and then all pij. + Outputs ps[i][j] the probability to be observed in j being in j according to + the values of the covariates cov[nc] and corresponding parameter values x[nc+shiftij] + */ + double s1, lnpijopii; /*double t34;*/ int i,j,j1, nc, ii, jj; for(i=1; i<= nlstate; i++){ for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */ + for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){ + /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/ + lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc]; +/* printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */ } - ps[i][j]=s2; + ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */ } } - /*ps[3][2]=1;*/ for(i=1; i<= nlstate; i++){ s1=0; for(j=1; ji} pij/pii=(1-pii)/pii and thus pii is known from s1 */ ps[i][i]=1./(s1+1.); + /* Computing other pijs */ for(j=1; j1){ - 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 \ @@ -4672,7 +4714,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 \ @@ -4734,6 +4776,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 */ @@ -4874,7 +4921,7 @@ calandcheckages(int imx, int maxwav, dou } else if(agev[m][i] >*agemax){ *agemax=agev[m][i]; - /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/ + printf(" Max anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.2f\n",m,i,anint[m][i], i,annais[i], *agemax); } /*agev[m][i]=anint[m][i]-annais[i];*/ /* agev[m][i] = age[i]+2*m;*/ @@ -4908,8 +4955,8 @@ calandcheckages(int imx, int maxwav, dou }*/ - printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); - fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); + printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax); + fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax); return (0); endread: @@ -5130,7 +5177,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); } @@ -5146,7 +5193,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); } @@ -5199,7 +5247,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); } @@ -5249,7 +5297,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); } @@ -5290,7 +5338,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); } @@ -5381,12 +5429,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 @@ -5445,11 +5495,10 @@ run imach with mle=-1 to get a correct t 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 */ h++; - if (h>m) { + if (h>m) h=1; - codtab[h][k]=j; - codtab[h][Tvar[k]]=j; - } + 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]]); } } @@ -5477,7 +5526,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 --------*/ @@ -5916,7 +5966,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); @@ -5936,7 +5986,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); @@ -5950,7 +6000,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); @@ -5966,7 +6016,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); @@ -6343,10 +6393,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); @@ -6387,9 +6436,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);