--- imach/src/imach.c 2007/08/31 13:49:27 1.129 +++ imach/src/imach.c 2009/06/20 16:22:47 1.131 @@ -1,6 +1,14 @@ -/* $Id: imach.c,v 1.129 2007/08/31 13:49:27 lievre Exp $ +/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.131 2009/06/20 16:22:47 brouard + Some dimensions resccaled + + Revision 1.130 2009/05/26 06:44:34 brouard + (Module): Max Covariate is now set to 20 instead of 8. A + lot of cleaning with variables initialized to 0. Trying to make + V2+V3*age+V1+V4 strb=V3*age+V1+V4 working better. + Revision 1.129 2007/08/31 13:49:27 lievre Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting @@ -367,13 +375,13 @@ extern int errno; #define GLOCK_ERROR_NOPATH -1 /* empty path */ #define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */ -#define MAXPARM 30 /* Maximum number of parameters for the optimization */ +#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 8 /* Maximum number of covariates */ +#define NCOVMAX 20 /* Maximum number of covariates */ #define MAXN 20000 #define YEARM 12. /* Number of months per year */ #define AGESUP 130 @@ -389,41 +397,41 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.129 2007/08/31 13:49:27 lievre Exp $ */ +/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.129 $ $Date: 2007/08/31 13:49:27 $"; +char version[]="Imach version 0.98k, June 2006, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.131 $ $Date: 2009/06/20 16:22:47 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; -int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ -int nvar; -int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; +int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ +int nvar=0; +int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov=0; /* Number of covariates, of covariates with '*age' */ int npar=NPARMAX; int nlstate=2; /* Number of live states */ int ndeath=1; /* Number of dead states */ -int ncovmodel, ncovcol; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ +int ncovmodel=0, ncovcol=0; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */ int popbased=0; int *wav; /* Number of waves for this individuual 0 is possible */ -int maxwav; /* Maxim number of waves */ -int jmin, jmax; /* min, max spacing between 2 waves */ -int ijmin, ijmax; /* Individuals having jmin and jmax */ -int gipmx, gsw; /* Global variables on the number of contributions +int maxwav=0; /* Maxim number of waves */ +int jmin=0, jmax=0; /* min, max spacing between 2 waves */ +int ijmin=0, ijmax=0; /* Individuals having jmin and jmax */ +int gipmx=0, gsw=0; /* Global variables on the number of contributions to the likelihood and the sum of weights (done by funcone)*/ -int mle, weightopt; +int mle=1, weightopt=0; int **mw; /* mw[mi][i] is number of the mi wave for this individual */ int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ 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. */ -double jmean; /* Mean space between 2 waves */ +double jmean=1; /* Mean space between 2 waves */ double **oldm, **newm, **savm; /* Working pointers to matrices */ double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ FILE *fic,*ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficlog, *ficrespow; -int globpr; /* Global variable for printing or not */ +int globpr=0; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ -long ipmx; /* Number of contributions */ +long ipmx=0; /* Number of contributions */ double sw; /* Sum of weights */ char filerespow[FILENAMELENGTH]; char fileresilk[FILENAMELENGTH]; /* File of individual contributions to the likelihood */ @@ -516,7 +524,7 @@ double dateintmean=0; double *weight; int **s; /* Status */ double *agedc, **covar, idx; -int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; +int **nbcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff; double *lsurv, *lpop, *tpop; double ftol=FTOL; /* Tolerance for computing Max Likelihood */ @@ -1166,7 +1174,7 @@ double **prevalim(double **prlim, int nl int i, ii,j,k; double min, max, maxmin, maxmax,sumnew=0.; double **matprod2(); - double **out, cov[NCOVMAX], **pmij(); + double **out, cov[NCOVMAX+1], **pmij(); double **newm; double agefin, delaymax=50 ; /* Max number of years to converge */ @@ -1248,10 +1256,14 @@ double **pmij(double **ps, double *cov, for(i=1; i<= nlstate; i++){ s1=0; - for(j=1; j cptcode) cptcode=ij; /* getting the maximum of covariable + if (ij > cptcode) cptcode=ij; /* getting the maximum value of the modality of the covariate (should be 0 or 1 now) Tvar[j]. If V=sex and male is 0 and female is 1, then cptcode=1.*/ } - for (i=0; i<=cptcode; i++) { - if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */ - } + for (i=0; i<=cptcode; i++) { /* i=-1 ?*/ + if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariate. In fact ncodemax[j]=2 (dichotom. variables only) but it can be more */ + } /* Ndum[-1] number of undefined modalities */ ij=1; - for (i=1; i<=ncodemax[j]; i++) { - for (k=0; k<= maxncov; k++) { - if (Ndum[k] != 0) { - nbcode[Tvar[j]][ij]=k; - /* store the modality in an array. k is a modality. If we have model=V1+V1*sex then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ - + for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 */ + for (k=0; k<= maxncov; k++) { /* k=-1 ?*/ + if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */ + nbcode[Tvar[j]][ij]=k; /* stores the modality in an array nbcode. + k is a modality. If we have model=V1+V1*sex + then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ ij++; } if (ij > ncodemax[j]) break; @@ -2468,9 +2483,9 @@ void tricode(int *Tvar, int **nbcode, in for (k=0; k< maxncov; k++) Ndum[k]=0; - for (i=1; i<=ncovmodel-2; i++) { + for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */ /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ - ij=Tvar[i]; + ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ Ndum[ij]++; } @@ -2481,8 +2496,8 @@ void tricode(int *Tvar, int **nbcode, in ij++; } } - - cptcoveff=ij-1; /*Number of simple covariates*/ + ij--; + cptcoveff=ij; /*Number of simple covariates*/ } /*********** Health Expectancies ****************/ @@ -3080,9 +3095,9 @@ void varevsij(char optionfilefiname[], d free_vector(gmp,nlstate+1,nlstate+ndeath); free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ - fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65"); + fprintf(ficgp,"\nunset parametric;unset label; set ter png small;set size 0.65, 0.65"); /* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ - fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); + fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); /* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ /* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ @@ -3259,7 +3274,7 @@ void varprob(char optionfilefiname[], do fprintf(ficresprobcov,"\n"); fprintf(ficresprobcor,"\n"); */ - xp=vector(1,npar); + xp=vector(1,npar); dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar); doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath)); mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage); @@ -3412,7 +3427,7 @@ To be simple, these graphs help to under /* Confidence intervalle of pij */ /* - fprintf(ficgp,"\nset noparametric;unset label"); + fprintf(ficgp,"\nunset parametric;unset label"); fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\""); fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); fprintf(fichtm,"\n
Probability with confidence intervals expressed in year-1 :pijgr%s.png, ",optionfilefiname,optionfilefiname); @@ -3647,8 +3662,8 @@ true period expectancies (those weighted void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ char dirfileres[132],optfileres[132]; - int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; - int ng; + int m0,cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; + int ng=0; /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ /* printf("Problem with file %s",optionfilegnuplot); */ /* fprintf(ficlog,"Problem with file %s",optionfilegnuplot); */ @@ -3674,7 +3689,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u for (i=1; i<= nlstate ; i ++) { if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); - else 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); for (i=1; i<= nlstate ; i ++) { @@ -4624,7 +4639,13 @@ int main(int argc, char *argv[]) ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ - + if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){ + printf("Too complex model for current IMaCh: npar=(nlstate+ndeath-1)*nlstate*ncovmodel=%d >= %d(MAXPARM) or nlstate=%d >= %d(NLSTATEMAX) or ndeath=%d >= %d(NDEATHMAX) or ncovmodel=(k+age+#of+signs)=%d(NCOVMAX) >= %d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); + fprintf(ficlog,"Too complex model for current IMaCh: %d >=%d(MAXPARM) or %d >=%d(NLSTATEMAX) or %d >=%d(NDEATHMAX) or %d(NCOVMAX) >=%d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX); + fflush(stdout); + fclose (ficlog); + goto end; + } delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); delti=delti3[1][1]; /*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ @@ -4750,6 +4771,9 @@ run imach with mle=-1 to get a correct t ungetc(c,ficpar); matcov=matrix(1,npar,1,npar); + for(i=1; i <=npar; i++) + for(j=1; j <=npar; j++) matcov[i][j]=0.; + for(i=1; i <=npar; i++){ fscanf(ficpar,"%s",&str); if(mle==1) @@ -4813,7 +4837,7 @@ run imach with mle=-1 to get a correct t for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */ mint=matrix(1,maxwav,1,n); anint=matrix(1,maxwav,1,n); - s=imatrix(1,maxwav+1,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); @@ -4837,12 +4861,17 @@ run imach with mle=-1 to get a correct t for (j=maxwav;j>=1;j--){ cutv(stra, strb,line,' '); - errno=0; - lval=strtol(strb,&endptr,10); + if(strb[0]=='.') { /* Missing status */ + lval=-1; + }else{ + errno=0; + 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 '%d' at line number %d %s 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); - exit(1); + 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); + goto end; + } } s[j][i]=lval; @@ -4854,8 +4883,9 @@ run imach with mle=-1 to get a correct t month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j); - exit(1); + 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); + goto end; } anint[j][i]= (double) year; mint[j][i]= (double)month; @@ -4869,8 +4899,9 @@ run imach with mle=-1 to get a correct t month=99; year=9999; }else{ - printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line); - exit(1); + 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); + goto end; } andc[i]=(double) year; moisdc[i]=(double) month; @@ -4884,7 +4915,8 @@ run imach with mle=-1 to get a correct t year=9999; }else{ printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j); - exit(1); + fprintf(ficlog,"Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);fflush(ficlog); + goto end; } annais[i]=(double)(year); moisnais[i]=(double)(month); @@ -4895,18 +4927,25 @@ run imach with mle=-1 to get a correct t 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); - exit(1); + 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); + fflush(ficlog); + goto end; } weight[i]=dval; strcpy(line,stra); for (j=ncovcol;j>=1;j--){ cutv(stra, strb,line,' '); - errno=0; - lval=strtol(strb,&endptr,10); - if( strb[0]=='\0' || (*endptr != '\0')){ - printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line); - exit(1); + if(strb[0]=='.') { /* Missing status */ + lval=-1; + }else{ + 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); + goto end; + } } if(lval <-1 || lval >1){ printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \ @@ -4918,6 +4957,15 @@ run imach with mle=-1 to get a correct t 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 \ + 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 \ + build V1=0 V2=0 for the reference value (1),\n \ + V1=1 V2=0 for (2) \n \ + 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);fflush(ficlog); goto end; } covar[j][i]=(double)(lval); @@ -4956,7 +5004,7 @@ run imach with mle=-1 to get a correct t else weight[i]=1;*/ /* Calculation of the number of parameters from char model */ - Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */ + Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. Stores the number n of the covariates in Vm+Vn at 1 and m at 2 */ Tprod=ivector(1,15); Tvaraff=ivector(1,15); Tvard=imatrix(1,15,1,2); @@ -4966,32 +5014,35 @@ run imach with mle=-1 to get a correct t j=0, j1=0, k1=1, k2=1; j=nbocc(model,'+'); /* j=Number of '+' */ j1=nbocc(model,'*'); /* j1=Number of '*' */ - cptcovn=j+1; - cptcovprod=j1; /*Number of products */ + cptcovn=j+1; /* Number of covariates V1+V2+V3 =>2+1=3 */ + cptcovprod=j1; /*Number of products V1*V2 =1 */ strcpy(modelsav,model); if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ printf("Error. Non available option model=%s ",model); - fprintf(ficlog,"Error. Non available option model=%s ",model); + fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog); goto end; } /* 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 */ for(i=(j+1); i>=1;i--){ - cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ + cutv(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' + modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 + stra=V2 + */ if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ /* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ /*scanf("%d",i);*/ - if (strchr(strb,'*')) { /* Model includes a product */ - cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn (if not *age)*/ + if (strchr(strb,'*')) { /* Model includes a product V1+V3*age+V2 strb=V3*age*/ + cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn: V3*age strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ if (strcmp(strc,"age")==0) { /* Vn*age */ cptcovprod--; cutv(strb,stre,strd,'V'); - Tvar[i]=atoi(stre); /* computes n in Vn and stores in Tvar*/ - cptcovage++; - Tage[cptcovage]=i; - /*printf("stre=%s ", stre);*/ + Tvar[i]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3 */ + cptcovage++; /* Sum the number of covariates including ages as a product */ + Tage[cptcovage]=i; /* Tage[1] =2 */ + /*printf("stre=%s ", stre);*/ } else if (strcmp(strd,"age")==0) { /* or age*Vn */ cptcovprod--; @@ -5000,11 +5051,12 @@ run imach with mle=-1 to get a correct t cptcovage++; Tage[cptcovage]=i; } - else { /* Age is not in the model */ - cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n*/ - Tvar[i]=ncovcol+k1; + else { /* Age is not in the model V1+V3*V2+V2 strb=V3*V2*/ + cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ + Tvar[i]=ncovcol+k1; /* find 'n' in Vn and stores in Tvar. + If already ncovcol=2 and model=V2*V1 Tvar[1]=2+1 and Tvar[2]=2+2 etc */ cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */ - Tprod[k1]=i; + Tprod[k1]=i; /* Tprod[1] */ Tvard[k1][1]=atoi(strc); /* m*/ Tvard[k1][2]=atoi(stre); /* n */ Tvar[cptcovn+k2]=Tvard[k1][1]; @@ -5021,7 +5073,7 @@ run imach with mle=-1 to get a correct t cutv(strd,strc,strb,'V'); Tvar[i]=atoi(strc); } - strcpy(modelsav,stra); + strcpy(modelsav,stra); /* modelsav=V2+V3*age+V1+V4 strb=V3*age+V1+V4 */ /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); scanf("%d",i);*/ } /* end of loop + */ @@ -5151,7 +5203,6 @@ run imach with mle=-1 to get a correct t /* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ - Tcode=ivector(1,100); nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); ncodemax[1]=1; if (cptcovn > 0) tricode(Tvar,nbcode,imx); @@ -5161,13 +5212,14 @@ run imach with mle=-1 to get a correct t h=0; m=pow(2,cptcoveff); - for(k=1;k<=cptcoveff; k++){ - for(i=1; i <=(m/pow(2,k));i++){ - for(j=1; j <= ncodemax[k]; j++){ - for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ + 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 */ h++; - if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; - /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/ + if (h>m) + h=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]]); } } } @@ -5176,7 +5228,7 @@ run imach with mle=-1 to get a correct t codtab[1][2]=1;codtab[2][2]=2; */ /* for(i=1; i <=m ;i++){ for(k=1; k <=cptcovn; k++){ - printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); + printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); } printf("\n"); } @@ -5204,7 +5256,8 @@ run imach with mle=-1 to get a correct t strcat(optionfilehtm,"-mort"); strcat(optionfilehtm,".htm"); if((fichtm=fopen(optionfilehtm,"w"))==NULL) { - printf("Problem with %s \n",optionfilehtm), exit(0); + printf("Problem with %s \n",optionfilehtm); + exit(0); } strcpy(optionfilehtmcov,optionfilefiname); /* Only for matrix of covariance */ @@ -5383,8 +5436,8 @@ Interval (in months) between two waves: } /* Endof if mle==-3 */ else{ /* For mle >=1 */ - - likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ + globpr=1;/* debug */ + /* likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone);*/ /* Prints the contributions to the likelihood */ printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); for (k=1; k<=npar;k++) printf(" %d %8.5f",k,p[k]); @@ -5656,7 +5709,8 @@ Interval (in months) between two waves: for(cptcov=1,k=0;cptcov<=i1;cptcov++){ for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ k=k+1; - /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/ + /* to clean */ + printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,codtab[cptcod][cptcov],nbcode); fprintf(ficrespl,"\n#******"); printf("\n#******"); fprintf(ficlog,"\n#******"); @@ -5968,6 +6022,7 @@ Interval (in months) between two waves: 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(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); @@ -5985,7 +6040,6 @@ Interval (in months) between two waves: free_ivector(Tprod,1,15); free_ivector(Tvaraff,1,15); free_ivector(Tage,1,15); - free_ivector(Tcode,1,100); free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX); free_imatrix(codtab,1,100,1,10);