--- imach/src/imach.c 2016/07/22 17:45:30 1.228 +++ imach/src/imach.c 2016/08/22 06:55:53 1.230 @@ -1,6 +1,12 @@ -/* $Id: imach.c,v 1.228 2016/07/22 17:45:30 brouard Exp $ +/* $Id: imach.c,v 1.230 2016/08/22 06:55:53 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.230 2016/08/22 06:55:53 brouard + Summary: Not working + + Revision 1.229 2016/07/23 09:45:53 brouard + Summary: Completing for func too + Revision 1.228 2016/07/22 17:45:30 brouard Summary: Fixing some arrays, still debugging @@ -884,12 +890,12 @@ typedef struct { #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.228 2016/07/22 17:45:30 brouard Exp $ */ +/* $Id: imach.c,v 1.230 2016/08/22 06:55:53 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018"; -char fullversion[]="$Revision: 1.228 $ $Date: 2016/07/22 17:45:30 $"; +char fullversion[]="$Revision: 1.230 $ $Date: 2016/08/22 06:55:53 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -927,6 +933,8 @@ int **dh; /* dh[mi][i] is number of step int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between * wave mi and wave mi+1 is not an exact multiple of stepm. */ int countcallfunc=0; /* Count the number of calls to func */ +int selected(int kvar); /* Is covariate kvar selected for printing results */ + double jmean=1; /* Mean space between 2 waves */ double **matprod2(); /* test */ double **oldm, **newm, **savm; /* Working pointers to matrices */ @@ -1056,14 +1064,16 @@ double ***cotvar; /* Time varying covari double ***cotqvar; /* Time varying quantitative covariate itqv */ double idx; int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ +int *Tvarsel; /**< Selected covariates for output */ +double *Tvalsel; /**< Selected modality value of covariate for output */ int *Typevar; /**< 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product */ int *Fixed; /** Fixed[k] 0=fixed, 1 varying, 2 fixed with age product, 3 varying with age product */ int *Dummy; /** Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product */ int *Tage; int anyvaryingduminmodel=0; /**< Any varying dummy in Model=1 yes, 0 no, to avoid a loop on waves in freq */ int *Tmodelind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ -int *TmodelInvind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ -int *TmodelInvQind; /** Tmodelqind[1]=1 for V5(quantitative varying) position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ +int *TmodelInvind; /** Tmodelind[Tvaraff[3]]=9 for V1 position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ +int *TmodelInvQind; /** Tmodelqind[1]=1 for V5(quantitative varying) position,Tvaraff[1]@9={4, 3, 1, 0, 0, 0, 0, 0, 0}, model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ int *Ndum; /** Freq of modality (tricode */ /* int **codtab;*/ /**< codtab=imatrix(1,100,1,10); */ int **Tvard; @@ -2997,13 +3007,15 @@ double func( double *x) */ for(mi=1; mi<= wav[i]-1; mi++){ for(itv=1; itv <= ntveff; itv++){ /* Varying dummy covariates */ - cov[ioffset+itv]=cotvar[mw[mi][i]][Tvar[itv]][i]; /* Not sure, Tvar V4+V3+V5 Tvaraff ? */ + /* cov[ioffset+itv]=cotvar[mw[mi][i]][Tvar[itv]][i]; /\* Not sure, Tvar V4+V3+V5 Tvaraff ? *\/ */ + cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; } for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */ if(cotqvar[mw[mi][i]][iqtv][i] == -1){ printf("i=%d, mi=%d, iqtv=%d, cotqvar[mw[mi][i]][iqtv][i]=%f",i,mi,iqtv,cotqvar[mw[mi][i]][iqtv][i]); } - cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i]; + cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; + /* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][iqtv][i]; */ } /* ioffset=2+nagesqr+cptcovn+nqv+ntv+nqtv; */ for (ii=1;ii<=nlstate+ndeath;ii++) @@ -3302,11 +3314,11 @@ double funcone( double *x) for (i=1,ipmx=0, sw=0.; i<=imx; i++){ ioffset=2+nagesqr+cptcovage; /* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */ - for (k=1; k<=ncoveff+nqfveff;k++){ /* Simple and product fixed covariates without age* products */ - cov[++ioffset]=covar[Tvar[k]][i]; + for (k=1; k<=ncoveff+nqfveff;k++){ /* Simple and product fixed Dummy covariates without age* products */ + cov[++ioffset]=covar[TvarFD[k]][i];/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, only V1 is fixed (k=6)*/ } for(iqv=1; iqv <= nqfveff; iqv++){ /* Quantitative fixed covariates */ - cov[++ioffset]=coqvar[Tvar[iqv]][i]; + cov[++ioffset]=coqvar[Tvar[iqv]][i]; /* Only V1 k=9 */ } for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */ @@ -3315,11 +3327,11 @@ double funcone( double *x) /* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */ k=ioffset-2-nagesqr-cptcovage+itv; /* position in simple model */ cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; - printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); + /* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */ } for(iqtv=1; iqtv <= nqtveff; iqtv++){ /* Varying quantitatives covariates */ iv=TmodelInvQind[iqtv]; /* Counting the # varying covariate from 1 to ntveff */ - printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); + /* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); */ cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; } for (ii=1;ii<=nlstate+ndeath;ii++) @@ -4718,18 +4730,18 @@ void concatwav(int wav[], int **dh, int if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */ /* If product not in single variable we don't print results */ /*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ - ++ij; - Tvaraff[ij]=Tvar[k]; /*For printing */ - Tmodelind[ij]=k; - TmodelInvind[k]=Tvar[k]- ncovcol-nqv; + ++ij;/* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, */ + Tvaraff[ij]=Tvar[k]; /* For printing combination *//* V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/ + Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */ + TmodelInvind[k]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */ if(Fixed[k]!=0) anyvaryingduminmodel=1; - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */ - /* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */ - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */ - /* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */ - /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */ - /* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */ + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */ + /* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet *\/ */ + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */ + /* Tvaraff[++ij]=i; /\*For printing (unclear) *\/ */ + /* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */ + /* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet *\/ */ } } /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */ /* ij--; */ @@ -5895,6 +5907,7 @@ void printinghtml(char fileresu[], char fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); } + /* if(nqfveff+nqtveff 0) */ /* Test to be done */ fprintf(fichtm," ************\n
"); if(invalidvarcomb[k1]){ fprintf(fichtm,"\n

Combination (%d) ignored because no cases

\n",k1); @@ -6090,7 +6103,7 @@ void printinggnuplot(char fileresu[], ch strcpy(optfileres,"vpl"); /* 1eme*/ for (cpt=1; cpt<= nlstate ; cpt ++) { /* For each live state */ - for (k1=1; k1<= m ; k1 ++) { /* For each valid combination of covariate */ + for (k1=1; k1<= m && selected(k1) ; k1 ++) { /* For each valid combination of covariate */ /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files "); for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */ @@ -7774,16 +7787,58 @@ int readdata(char datafile[], int firsto return (1); } -void removespace(char *str) { - char *p1 = str, *p2 = str; +void removespace(char **stri){/*, char stro[]) {*/ + char *p1 = *stri, *p2 = *stri; do while (*p2 == ' ') p2++; while (*p1++ == *p2++); + *stri=p1; } -int decodemodel ( char model[], int lastobs) - /**< This routine decode the model and returns: +int decoderesult ( char resultline[]) +/**< This routine decode one result line and returns the combination # of dummy covariates only **/ +{ + int j=0, k=0; + char resultsav[MAXLINE]; + char stra[80], strb[80], strc[80], strd[80],stre[80]; + + removespace(&resultline); + printf("decoderesult=%s\n",resultline); + + if (strstr(resultline,"v") !=0){ + printf("Error. 'v' must be in upper case 'V' result: %s ",resultline); + fprintf(ficlog,"Error. 'v' must be in upper case result: %s ",resultline);fflush(ficlog); + return 1; + } + trimbb(resultsav, resultline); + if (strlen(resultsav) >1){ + j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' */ + } + + for(k=1; k<=j;k++){ /* Loop on total covariates of the model */ + cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' + resultsav= V4=1 V5=25.1 V3=0 strb=V3=0 stra= V4=1 V5=25.1 */ + cutl(strc,strd,strb,'='); /* strb:V4=1 strc=1 strd=V4 */ + Tvalsel[k]=atof(strc); /* 1 */ + + cutl(strc,stre,strd,'V'); /* strd='V4' strc=4 stre='V' */; + Tvarsel[k]=atoi(strc); + /* Typevarsel[k]=1; /\* 1 for age product *\/ */ + /* cptcovsel++; */ + if (nbocc(stra,'=') >0) + strcpy(resultsav,stra); /* and analyzes it */ + } + return (0); +} +int selected( int kvar){ /* Selected combination of covariates */ + if(Tvarsel[kvar]) + return (0); + else + return(1); +} +int decodemodel( char model[], int lastobs) + /**< This routine decodes the model and returns: * Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age * - nagesqr = 1 if age*age in the model, otherwise 0. * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age @@ -7981,8 +8036,8 @@ int decodemodel ( char model[], int last scanf("%d ",i);*/ -/* Decodemodel knows only the grammar (simple, product, age*) of the model but not what kind - of variable (dummy vs quantitative, fixed vs time varying) is behind */ +/* Until here, decodemodel knows only the grammar (simple, product, age*) of the model but not what kind + of variable (dummy vs quantitative, fixed vs time varying) is behind. But we know the # of each. */ /* ncovcol= 1, nqv=1 | ntv=2, nqtv= 1 = 5 possible variables data: 2 fixed 3, varying model= V5 + V4 +V3 + V4*V3 + V5*age + V2 + V1*V2 + V1*age + V5*age, V1 is not used saving its place k = 1 2 3 4 5 6 7 8 9 @@ -8010,14 +8065,20 @@ Dummy[k] 0=dummy (0 1), 1 quantitative ( Fixed[k]= 0; Dummy[k]= 0; ncoveff++; + TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ + TvarFDind[ncoveff]=Tvar[k]; /* TvarFDind[1]=9 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ }else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product Vn*Vm are added in k*/ Fixed[k]= 0; Dummy[k]= 1; - nqfveff++; /* Only simple fixed quantitative variable */ + nqfveff++; + TvarFQ[nqfveff]=Tvar[k]; /* TvarFQ[1]=V2 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ + TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ }else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){ Fixed[k]= 1; Dummy[k]= 0; ntveff++; /* Only simple time varying dummy variable */ + TvarVD[ntvveff]=Tvar[k]; /* TvarVD[1]=V5 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ + TvarVDind[ntveff++]=k; /* TvarVDind[1]=1 in V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 */ /* Only simple fixed quantitative variable */ printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv); printf("Quasi TmodelInvind[%d]=%d\n",k,Tvar[k]- ncovcol-nqv); }else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv && Typevar[k]==0){ @@ -8821,6 +8882,8 @@ int main(int argc, char *argv[]) char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE]; char model[MAXLINE], modeltemp[MAXLINE]; + char resultline[MAXLINE]; + char pathr[MAXLINE], pathimach[MAXLINE]; char *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; @@ -9376,6 +9439,8 @@ Please run with mle=-1 to get a correct k=1 Tvar[1]=2 (from V2) */ Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. */ + Tvarsel=ivector(1,NCOVMAX); /* */ + Tvalsel=vector(1,NCOVMAX); /* */ Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */ Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */ Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */ @@ -9402,13 +9467,15 @@ Please run with mle=-1 to get a correct 4 covariates (3 plus signs) Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 */ - Tmodelind=ivector(1,NCOVMAX);/** five the k model position of an + Tmodelind=ivector(1,NCOVMAX);/** gives the k model position of an * individual dummy, fixed or varying: * Tmodelind[Tvaraff[3]]=9,Tvaraff[1]@9={4, * 3, 1, 0, 0, 0, 0, 0, 0}, - * model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1*/ - TmodelInvind=ivector(1,NCOVMAX); - TmodelInvQind=ivector(1,NCOVMAX);/** five the k model position of an + * model=V5+V4+V3+V4*V3+V5*age+V2+V1*V2+V1*age+V1 , + * V1 df, V2 qf, V3 & V4 dv, V5 qv + * Tmodelind[1]@9={9,0,3,2,}*/ + TmodelInvind=ivector(1,NCOVMAX); /* TmodelInvind=Tvar[k]- ncovcol-nqv={5-2-1=2,*/ + TmodelInvQind=ivector(1,NCOVMAX);/** gives the k model position of an * individual quantitative, fixed or varying: * Tmodelqind[1]=1,Tvaraff[1]@9={4, * 3, 1, 0, 0, 0, 0, 0, 0}, @@ -10178,16 +10245,54 @@ Please run with mle=-1 to get a correct fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj); /* day and month of proj2 are not used but only year anproj2.*/ + /* Results */ + while(fgets(line, MAXLINE, ficpar)) { + /* If line starts with a # it is a comment */ + if (line[0] == '#') { + numlinepar++; + fputs(line,stdout); + fputs(line,ficparo); + fputs(line,ficlog); + continue; + }else + break; + } + while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){ + if (num_filled == 0) + resultline[0]='\0'; + else if (num_filled != 1){ + printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line); + } + printf("Result %d: result line should be at minimum 'line=' %s, result=%s\n",num_filled, line, resultline); + decoderesult(resultline); + while(fgets(line, MAXLINE, ficpar)) { + /* If line starts with a # it is a comment */ + if (line[0] == '#') { + numlinepar++; + fputs(line,stdout); + fputs(line,ficparo); + fputs(line,ficlog); + continue; + }else + break; + } + if (feof(ficpar)) + break; + else{ /* Processess output results for this combination of covariate values */ + } + } + + - /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */ + /* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */ /* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */ replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){ - printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ + printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); - fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ + fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); }else{ @@ -10557,6 +10662,8 @@ Please run with mle=-1 to get a correct free_ivector(Fixed,-1,NCOVMAX); free_ivector(Typevar,-1,NCOVMAX); free_ivector(Tvar,1,NCOVMAX); + free_ivector(Tvarsel,1,NCOVMAX); + free_vector(Tvalsel,1,NCOVMAX); free_ivector(Tposprod,1,NCOVMAX); free_ivector(Tprod,1,NCOVMAX); free_ivector(Tvaraff,1,NCOVMAX);