--- imach/src/imach.c 2023/04/24 11:38:06 1.350 +++ imach/src/imach.c 2023/05/22 17:03:18 1.355 @@ -1,6 +1,24 @@ -/* $Id: imach.c,v 1.350 2023/04/24 11:38:06 brouard Exp $ +/* $Id: imach.c,v 1.355 2023/05/22 17:03:18 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.355 2023/05/22 17:03:18 brouard + Summary: 0.99r46 + + * imach.c (Module): In the ILK....txt file, the number of columns + before the covariates values is dependent of the number of states (16+nlstate): 0.99r46 + + Revision 1.354 2023/05/21 05:05:17 brouard + Summary: Temporary change for imachprax + + Revision 1.353 2023/05/08 18:48:22 brouard + *** empty log message *** + + Revision 1.352 2023/04/29 10:46:21 brouard + *** empty log message *** + + Revision 1.351 2023/04/29 10:43:47 brouard + Summary: 099r45 + Revision 1.350 2023/04/24 11:38:06 brouard *** empty log message *** @@ -1361,12 +1379,12 @@ double gnuplotversion=GNUPLOTVERSION; #define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.350 2023/04/24 11:38:06 brouard Exp $ */ +/* $Id: imach.c,v 1.355 2023/05/22 17:03:18 brouard Exp $ */ /* $State: Exp $ */ #include "version.h" char version[]=__IMACH_VERSION__; -char copyright[]="January 2023,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-2022"; -char fullversion[]="$Revision: 1.350 $ $Date: 2023/04/24 11:38:06 $"; +char copyright[]="April 2023,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-2022"; +char fullversion[]="$Revision: 1.355 $ $Date: 2023/05/22 17:03:18 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -1603,7 +1621,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tv /* Tage[cptcovage]=k 5 8 10 */ /* Position in the model of ith cov*age */ /* model="V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/ /* p Tvard[1][1]@21 = {6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}*/ -/* p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 } +/* p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 } */ /* p Tvardk[1][1]@24 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0}*/ /* p Tvardk[1][1]@22 = {0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0} */ /* Tvard[1][1]@4={4,3,1,2} V4*V3 V1*V2 */ /* Position in model of the ith prod without age */ @@ -1804,6 +1822,20 @@ char *trimbb(char *out, char *in) return s; } +char *trimbtab(char *out, char *in) +{ /* Trim blanks or tabs in line but keeps first blanks if line starts with blanks */ + char *s; + s=out; + while (*in != '\0'){ + while( (*in == ' ' || *in == '\t')){ /* && *(in+1) != '\0'){*/ + in++; + } + *out++ = *in++; + } + *out='\0'; + return s; +} + /* char *substrchaine(char *out, char *in, char *chain) */ /* { */ /* /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */ @@ -4530,7 +4562,7 @@ double funcone( double *x) * 3 ncovta=15 +age*V3*V2+age*V2+agev3+ageV4 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 * 3 TvarAVVA[1]@15= itva 3 2 2 3 4 6 7 6 3 7 3 6 4 7 4 * 3 ncovta 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 - * TvarAVVAind[1]@15= V3 is in k=2 1 1 2 3 4 5 4,2 5,2, 4,3 5 3}TvarVVAind[] + *?TvarAVVAind[1]@15= V3 is in k=2 1 1 2 3 4 5 4,2 5,2, 4,3 5 3}TvarVVAind[] * TvarAVVAind[1]@15= V3 is in k=6 6 12 13 14 15 16 18 18 19,19, 20,20 21,21}TvarVVAind[] * 3 ncovvta=10 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva] @@ -4545,6 +4577,7 @@ double funcone( double *x) * 6, 8, 9, 10, 11} * TvarFind[itv] 0 0 0 * FixedV[itv] 1 1 1 0 1 0 1 0 1 0 0 + *? FixedV[itv] 1 1 1 0 1 0 1 0 1 0 1 0 1 0 * Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4 * Tvar[TvarFind[itv]] [0]=? ?ncovv 1 à ncovvt] * Not a fixed cotvar[mw][itv][i] 6 7 6 2 7, 2, 6, 3, 7, 3, 6, 4, 7, 4} @@ -4556,10 +4589,14 @@ double funcone( double *x) ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/ /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */ if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */ + /* printf("DEBUG ncovv=%d, Varying TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */ cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i]; /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */ + /* printf("DEBUG Varying cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */ }else{ /* fixed covariate */ /* cotvarv=covar[Tvar[TvarFind[itv]]][i]; /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */ + /* printf("DEBUG ncovv=%d, Fixed TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */ cotvarv=covar[itv][i]; /* Good: In V6*V3, 3 is fixed at position of the data */ + /* printf("DEBUG Fixed cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */ } if(ipos!=iposold){ /* Not a product or first of a product */ cotvarvold=cotvarv; @@ -4568,6 +4605,7 @@ double funcone( double *x) } iposold=ipos; cov[ioffset+ipos]=cotvarv; + /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */ /* For products */ } /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */ @@ -4910,14 +4948,15 @@ void mlikeli(FILE *ficres,double p[], in double fret; double fretone; /* Only one call to likelihood */ /* char filerespow[FILENAMELENGTH];*/ - + + double * p1; /* Shifted parameters from 0 instead of 1 */ #ifdef NLOPT int creturn; nlopt_opt opt; /* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */ double *lb; double minf; /* the minimum objective value, upon return */ - double * p1; /* Shifted parameters from 0 instead of 1 */ + myfunc_data dinst, *d = &dinst; #endif @@ -5252,6 +5291,7 @@ double hessij( double x[], double **hess kmax=kmax+10; } if(kmax >=10 || firstime ==1){ + /* What are the thetai and thetaj? thetai/ncovmodel thetai=(thetai-thetai%ncovmodel)/ncovmodel +thetai%ncovmodel=(line,pos) */ printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); @@ -8337,8 +8377,8 @@ true period expectancies (those weighted /******************* Gnuplot file **************/ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int prevbcast, char pathc[], double p[], int offyear, int offbyear){ - char dirfileres[132],optfileres[132]; - char gplotcondition[132], gplotlabel[132]; + 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 ng=0; @@ -8409,7 +8449,8 @@ void printinggnuplot(char fileresu[], ch kvar=Tvar[TvarFind[kf]]; /* variable name */ /* k=18+Tvar[TvarFind[kf]];/\*offset because there are 18 columns in the ILK_ file but could be placed else where *\/ */ /* k=18+kf;/\*offset because there are 18 columns in the ILK_ file *\/ */ - k=19+kf;/*offset because there are 19 columns in the ILK_ file */ + /* k=19+nlstate+kf;/\*offset because there are 19 columns in the ILK_ file *\/ */ + k=16+nlstate+kf;/*offset because there are 19 columns in the ILK_ file, first cov Vn on col 21 with 4 living states */ for (i=1; i<= nlstate ; i ++) { fprintf(ficgp,"\nset out \"%s-p%dj-%d.png\";set ylabel \"Probability for each individual/wave\";",subdirf2(optionfilefiname,"ILK_"),i,kvar); fprintf(ficgp,"unset log;\n# For each simple dummy covariate of the model \n plot \"%s\"",subdirf(fileresilk)); @@ -9095,18 +9136,22 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp," u %d:(",ioffset); kl=0; strcpy(gplotcondition,"("); - for (k=1; k<=cptcoveff; k++){ /* For each covariate writing the chain of conditions */ + /* 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]]); + for (k=1; k<=cptcovs; k++){ /* For each covariate k get corresponding value lv for combination k1 */ + /* lv=codtabm(k1,TnsdVar[Tvaraff[k]]); */ + lv=Tvresult[nres][k]; + vlv=TinvDoQresult[nres][Tvresult[nres][k]]; /* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */ /* decodtabm(1,2,4) = 1 because h=1 k= 1 (1) 1 1 */ /* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */ /* vlv= nbcode[Tvaraff[k]][lv]; /\* Value of the modality of Tvaraff[k] *\/ */ - vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; + /* vlv= nbcode[Tvaraff[k]][codtabm(k1,TnsdVar[Tvaraff[k]])]; */ kl++; - 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,Tvaraff[k], kl+1, nbcode[Tvaraff[k]][lv]); */ + sprintf(gplotcondition+strlen(gplotcondition),"$%d==%d && $%d==%d " ,kl,lv, kl+1, vlv ); kl++; - if(k 1) + if(k 1) sprintf(gplotcondition+strlen(gplotcondition)," && "); } strcpy(gplotcondition+strlen(gplotcondition),")"); @@ -9190,7 +9235,8 @@ 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 */ @@ -9302,7 +9348,8 @@ set ter svg size 640, 480\nunset log y\n fprintf(ficgp,"#Number of graphics: first is logit, 2nd is probabilities, third is incidences per year\n"); fprintf(ficgp,"#model=1+age+%s \n",model); fprintf(ficgp,"# Type of graphic ng=%d\n",ng); - fprintf(ficgp,"# k1=1 to 2^%d=%d\n",cptcoveff,m);/* to be checked */ + /* fprintf(ficgp,"# k1=1 to 2^%d=%d\n",cptcoveff,m);/\* to be checked *\/ */ + fprintf(ficgp,"# k1=1 to 2^%d=%d\n",cptcovs,m);/* to be checked */ /* for(k1=1; k1 <=m; k1++) /\* For each combination of covariate *\/ */ for(nres=1; nres <= nresult; nres++){ /* For each resultline */ /* k1=nres; */ @@ -9920,30 +9967,36 @@ void prevforecast(char fileres[], double /* date2dmy(dateintmean,&jintmean,&mintmean,&aintmean); */ /* date2dmy(dateprojd,&jprojd, &mprojd, &anprojd); */ /* date2dmy(dateprojf,&jprojf, &mprojf, &anprojf); */ - i1=pow(2,cptcoveff); - if (cptcovn < 1){i1=1;} + /* i1=pow(2,cptcoveff); */ + /* if (cptcovn < 1){i1=1;} */ fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); fprintf(ficresf,"#****** Routine prevforecast **\n"); /* if (h==(int)(YEARM*yearp)){ */ - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ /* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) */ - if(i1 != 1 && TKresult[nres]!= k) - continue; - if(invalidvarcomb[k]){ - printf("\nCombination (%d) projection ignored because no cases \n",k); - continue; - } + for(nres=1; nres <= nresult; nres++){ /* For each resultline */ + k=TKresult[nres]; + if(TKresult[nres]==0) k=1; /* To be checked for noresult */ + /* for(k=1; k<=i1;k++){ /\* We want to find the combination k corresponding to the values of the dummies given in this resut line (to be cleaned one day) *\/ */ + /* if(i1 != 1 && TKresult[nres]!= k) */ + /* continue; */ + /* if(invalidvarcomb[k]){ */ + /* printf("\nCombination (%d) projection ignored because no cases \n",k); */ + /* continue; */ + /* } */ fprintf(ficresf,"\n#****** hpijx=probability over h years, hp.jx is weighted by observed prev \n#"); - for(j=1;j<=cptcoveff;j++) { - /* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); */ - fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - } - for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ - fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); + for(j=1;j<=cptcovs;j++){ + /* for(j=1;j<=cptcoveff;j++) { */ + /* /\* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); *\/ */ + /* fprintf(ficresf," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + /* } */ + /* for (k4=1; k4<= nsq; k4++){ /\* For each selected (single) quantitative value *\/ */ + /* fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */ + /* } */ + fprintf(ficresf," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); } + fprintf(ficresf," yearproj age"); for(j=1; j<=nlstate+ndeath;j++){ for(i=1; i<=nlstate;i++) @@ -9968,9 +10021,11 @@ void prevforecast(char fileres[], double } } fprintf(ficresf,"\n"); - for(j=1;j<=cptcoveff;j++) + /* for(j=1;j<=cptcoveff;j++) */ + for(j=1;j<=cptcovs;j++) + fprintf(ficresf,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,Tvaraff[j])]); /\* Tvaraff not correct *\/ */ - fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /* TnsdVar[Tvaraff] correct */ + /* fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); /\* TnsdVar[Tvaraff] correct *\/ */ fprintf(ficresf,"%.f %.f ",anprojd+yearp,agec+h*hstepm/YEARM*stepm); for(j=1; j<=nlstate+ndeath;j++) { @@ -10062,29 +10117,35 @@ void prevforecast(char fileres[], double /* if(jintmean==0) jintmean=1; */ /* if(mintmean==0) jintmean=1; */ - i1=pow(2,cptcoveff); - if (cptcovn < 1){i1=1;} + /* i1=pow(2,cptcoveff); */ + /* if (cptcovn < 1){i1=1;} */ fprintf(ficresfb,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); printf("# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jintmean,mintmean,aintmean,dateintmean,dateprev1,dateprev2); fprintf(ficresfb,"#****** Routine prevbackforecast **\n"); - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ - if(i1 != 1 && TKresult[nres]!= k) - continue; - if(invalidvarcomb[k]){ - printf("\nCombination (%d) projection ignored because no cases \n",k); - continue; - } + for(nres=1; nres <= nresult; nres++){ /* For each resultline */ + k=TKresult[nres]; + if(TKresult[nres]==0) k=1; /* To be checked for noresult */ + /* for(k=1; k<=i1;k++){ */ + /* if(i1 != 1 && TKresult[nres]!= k) */ + /* continue; */ + /* if(invalidvarcomb[k]){ */ + /* printf("\nCombination (%d) projection ignored because no cases \n",k); */ + /* continue; */ + /* } */ fprintf(ficresfb,"\n#****** hbijx=probability over h years, hb.jx is weighted by observed prev \n#"); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - } - for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */ - fprintf(ficresf," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); + for(j=1;j<=cptcovs;j++){ + /* for(j=1;j<=cptcoveff;j++) { */ + /* fprintf(ficresfb," V%d (=) %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + /* } */ + fprintf(ficresfb," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); } + /* fprintf(ficrespij,"******\n"); */ + /* for (k4=1; k4<= nsq; k4++){ /\* For each selected (single) quantitative value *\/ */ + /* fprintf(ficresfb," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */ + /* } */ fprintf(ficresfb," yearbproj age"); for(j=1; j<=nlstate+ndeath;j++){ for(i=1; i<=nlstate;i++) @@ -10115,8 +10176,10 @@ void prevforecast(char fileres[], double } } fprintf(ficresfb,"\n"); - for(j=1;j<=cptcoveff;j++) - fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); + /* for(j=1;j<=cptcoveff;j++) */ + for(j=1;j<=cptcovs;j++) + fprintf(ficresfb,"%d %lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); + /* fprintf(ficresfb,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ fprintf(ficresfb,"%.f %.f ",anbackd+yearp,agec-h*hstepm/YEARM*stepm); for(i=1; i<=nlstate+ndeath;i++) { ppij=0.;ppi=0.; @@ -11080,14 +11143,15 @@ int decoderesult( char resultline[], int if (strlen(resultsav) >1){ j=nbocc(resultsav,'='); /**< j=Number of covariate values'=' in this resultline */ } - if(j == 0){ /* Resultline but no = */ + if(j == 0 && cptcovs== 0){ /* Resultline but no = and no covariate in the model */ TKresult[nres]=0; /* Combination for the nresult and the model */ return (0); } if( j != cptcovs ){ /* Be careful if a variable is in a product but not single */ - printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model); - fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, %s.\n",j, cptcovs, model); - /* return 1;*/ + fprintf(ficlog,"ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(ficlog); + printf("ERROR: the number of variables in the resultline which is %d, differs from the number %d of single variables used in the model line, 1+age+%s.\n",j, cptcovs, model);fflush(stdout); + if(j==0) + return 1; } for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */ if(nbocc(resultsav,'=') >1){ @@ -11375,7 +11439,7 @@ int decodemodel( char model[], int lasto if (strlen(modelsav) >1){ /* V2 +V3 +V4 +V6 +V7 +V6*V2 +V7*V2 +V6*V3 +V7*V3 +V6*V4 +V7*V4 +age*V2 +age*V3 +age*V4 +age*V6 +age*V7 +age*V6*V2 +V7*V2 +age*V6*V3 +age*V7*V3 +age*V6*V4 +age*V7*V4 */ j=nbocc(modelsav,'+'); /**< j=Number of '+' */ j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */ - cptcovs=j+1-j1; /**< Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age => V1 + V3 =4+1-3=2 */ + cptcovs=0; /**< Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age => V1 + V3 =4+1-3=2 Wrong */ cptcovt= j+1; /* Number of total covariates in the model, not including * cst, age and age*age * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=4*/ @@ -11440,6 +11504,22 @@ int decodemodel( char model[], int lasto Tvar[k]=0; Tprod[k]=0; Tposprod[k]=0; } cptcovage=0; + + /* First loop in order to calculate */ + /* for age*VN*Vm + * Provides, Typevar[k], Tage[cptcovage], existcomb[n][m], FixedV[ncovcolt+k12] + * Tprod[k1]=k Tposprod[k]=k1; Tvard[k1][1] =m; + */ + /* Needs FixedV[Tvardk[k][1]] */ + /* For others: + * Sets Typevar[k]; + * Tvar[k]=ncovcol+nqv+ntv+nqtv+k11; + * Tposprod[k]=k11; + * Tprod[k11]=k; + * Tvardk[k][1] =m; + * Needs FixedV[Tvardk[k][1]] == 0 + */ + for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model line */ cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' cutl from left to right modelsav==V2+V1+V5*age+V4+V3*age strb=V3*age stra=V2+V1V5*age+V4 */ /* "V5+V4+V3+V4*V3+V5*age+V1*age+V1" strb="V5" stra="V4+V3+V4*V3+V5*age+V1*age+V1" */ @@ -11465,8 +11545,8 @@ int decodemodel( char model[], int lasto strcat(strb,stre); strcpy(strd,strb); /* in order for strd to not be "age" for next test (will be Vn*Vm */ } - printf("DEBUG FIXED k=%d, Tage[k]=%d, Tvar[Tage[k]=%d,FixedV[Tvar[Tage[k]]]=%d\n",k,Tage[k],Tvar[Tage[k]],FixedV[Tvar[Tage[k]]]); - FixedV[Tvar[Tage[k]]]=0; /* HERY not sure */ + /* printf("DEBUG FIXED k=%d, Tage[k]=%d, Tvar[Tage[k]=%d,FixedV[Tvar[Tage[k]]]=%d\n",k,Tage[k],Tvar[Tage[k]],FixedV[Tvar[Tage[k]]]); */ + /* FixedV[Tvar[Tage[k]]]=0; /\* HERY not sure if V7*V4*age Fixed might not exist yet*\/ */ }else{ /* strc=Vn*Vm (and strd=age) and should be strb=Vn*Vm but want to keep original strb double product */ strcpy(stre,strb); /* save full b in stre */ strcpy(strb,strc); /* save short c in new short b for next block strb=Vn*Vm*/ @@ -11499,7 +11579,7 @@ int decodemodel( char model[], int lasto Tvardk[k][1] =m; /* m 1 for V1*/ Tvard[k1][2] =n; /* n 4 for V4*/ Tvardk[k][2] =n; /* n 4 for V4*/ -/* Tvar[Tage[cptcovage]]=k1;*/ /* Tvar[6=age*V3*V2]=9 (new fixed covariate) */ +/* Tvar[Tage[cptcovage]]=k1;*/ /* Tvar[6=age*V3*V2]=9 (new fixed covariate) */ /* We don't know about Fixed yet HERE */ if( FixedV[Tvardk[k][1]] == 0 && FixedV[Tvardk[k][2]] == 0){ /* If the product is a fixed covariate then we feed the new column with Vn*Vm */ for (i=1; i<=lastobs;i++){/* For fixed product */ /* Computes the new covariate which is a product of @@ -11649,6 +11729,8 @@ int decodemodel( char model[], int lasto /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); scanf("%d",i);*/ } /* end of loop + on total covariates */ + + } /* end if strlen(modelsave == 0) age*age might exist */ } /* end if strlen(model == 0) */ cptcovs=cptcovt - cptcovdageprod - cptcovprod;/**< Number of simple covariates V1 +V1*age +V3 +V3*V4 +age*age + age*v4*V3=> V1 + V3 =4+1-3=2 */ @@ -11687,6 +11769,11 @@ Fixed[k] 0=fixed (product or simple), 1 Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model); for(k=-1;k<=NCOVMAX; k++){ Fixed[k]=0; Dummy[k]=0;} for(k=1;k<=NCOVMAX; k++){TvarFind[k]=0; TvarVind[k]=0;} + + + /* Second loop for calculating Fixed[k], Dummy[k]*/ + + for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0,ncovva=0,ncovvta=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0, ncovvt=0;k<=cptcovt; k++){ /* or cptocvt loop on k from model */ if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */ Fixed[k]= 0; @@ -12940,7 +13027,7 @@ int main(int argc, char *argv[]) /* double ***mobaverage; */ double wald; - char line[MAXLINE]; + char line[MAXLINE], linetmp[MAXLINE]; char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE]; char modeltemp[MAXLINE]; @@ -13273,7 +13360,18 @@ int main(int argc, char *argv[]) }else break; } - if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ + if((num_filled=sscanf(line,"model=%[^.\n]", model)) !=EOF){ /* Every character after model but dot and return */ + if (num_filled != 1){ + printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); + fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); + model[0]='\0'; + goto end; + }else{ + trimbtab(linetmp,line); /* Trims multiple blanks in line */ + strcpy(line, linetmp); + } + } + if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){ /* Every character after 1+age but dot and return */ if (num_filled != 1){ printf("ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age+' instead of '%s'\n",num_filled, line); @@ -13646,7 +13744,7 @@ Please run with mle=-1 to get a correct Tvard=imatrix(1,NCOVMAX,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) */ - Tvardk=imatrix(-1,NCOVMAX,1,2); + Tvardk=imatrix(0,NCOVMAX,1,2); Tage=ivector(1,NCOVMAX); /* 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 @@ -14901,18 +14999,21 @@ Please run with mle=-1 to get a correct pstamp(ficreseij); - i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */ - if (cptcovn < 1){i1=1;} + /* i1=pow(2,cptcoveff); /\* Number of combination of dummy covariates *\/ */ + /* if (cptcovn < 1){i1=1;} */ - for(nres=1; nres <= nresult; nres++) /* For each resultline */ - for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */ - if(i1 != 1 && TKresult[nres]!= k) - continue; + for(nres=1; nres <= nresult; nres++){ /* For each resultline */ + /* for(k=1; k<=i1;k++){ /\* For any combination of dummy covariates, fixed and varying *\/ */ + /* if(i1 != 1 && TKresult[nres]!= k) */ + /* continue; */ fprintf(ficreseij,"\n#****** "); printf("\n#****** "); - for(j=1;j<=cptcoveff;j++) { - fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); - printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); + for(j=1;j<=cptcovs;j++){ + /* for(j=1;j<=cptcoveff;j++) { */ + /* fprintf(ficreseij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ + fprintf(ficreseij," V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); + printf(" V%d=%lg ",Tvresult[nres][j],TinvDoQresult[nres][Tvresult[nres][j]]); + /* printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,TnsdVar[Tvaraff[j]])]); */ } for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */ printf(" V%d=%lg ",TvarsQ[j], TinvDoQresult[nres][TvarsQ[j]]); /* TvarsQ[j] gives the name of the jth quantitative (fixed or time v) */ @@ -15156,7 +15257,7 @@ Please run with mle=-1 to get a correct free_vector(weight,firstobs,lastobs); - free_imatrix(Tvardk,-1,NCOVMAX,1,2); + free_imatrix(Tvardk,0,NCOVMAX,1,2); free_imatrix(Tvard,1,NCOVMAX,1,2); free_imatrix(s,1,maxwav+1,firstobs,lastobs); free_matrix(anint,1,maxwav,firstobs,lastobs);