--- imach/src/imach.c 2009/06/20 16:22:47 1.131 +++ imach/src/imach.c 2009/07/06 08:22:05 1.132 @@ -1,6 +1,9 @@ -/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $ +/* $Id: imach.c,v 1.132 2009/07/06 08:22:05 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.132 2009/07/06 08:22:05 brouard + Many tings + Revision 1.131 2009/06/20 16:22:47 brouard Some dimensions resccaled @@ -397,11 +400,11 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $ */ +/* $Id: imach.c,v 1.132 2009/07/06 08:22:05 brouard Exp $ */ /* $State: Exp $ */ 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 fullversion[]="$Revision: 1.132 $ $Date: 2009/07/06 08:22:05 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -597,6 +600,20 @@ void replace_back_to_slash(char *s, char } } +char *trimbb(char *out, char *in) +{ /* Trim multiple blanks in line */ + char *s; + s=out; + while (*in != '\0'){ + while( *in == ' ' && *(in+1) == ' ' && *(in+1) != '\0'){ + in++; + } + *out++ = *in++; + } + *out='\0'; + return s; +} + int nbocc(char *s, char occ) { int i,j=0; @@ -1692,7 +1709,7 @@ double funcone( double *x) ipmx +=1; sw += weight[i]; 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]); + /*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\ %11.6f %11.6f %11.6f ", \ @@ -1897,7 +1914,7 @@ double hessii(double x[], double delta, int i; int l=1, lmax=20; double k1,k2; - double p2[NPARMAX+1]; + double p2[MAXPARM+1]; /* identical to x */ double res; double delt=0.0001, delts, nkhi=10.,nkhif=1., khi=1.e-4; double fx; @@ -1918,7 +1935,7 @@ double hessii(double x[], double delta, /*res= (k1-2.0*fx+k2)/delt/delt; */ res= (k1+k2)/delt/delt/2.; /* Divided by because L and not 2*L */ -#ifdef DEBUG +#ifdef DEBUGHESS printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx); #endif @@ -1944,7 +1961,7 @@ double hessij( double x[], double delti[ int i; int l=1, l1, lmax=20; double k1,k2,k3,k4,res,fx; - double p2[NPARMAX+1]; + double p2[MAXPARM+1]; int k; fx=func(x); @@ -2155,7 +2172,7 @@ void freqsummary(char fileres[], int ia pos += freq[jk][m][i]; if(pp[jk]>=1.e-10){ if(first==1){ - printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); + printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); } fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]); }else{ @@ -2464,7 +2481,11 @@ void tricode(int *Tvar, int **nbcode, in } 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 */ + 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; @@ -4431,7 +4452,8 @@ int main(int argc, char *argv[]) double ***mobaverage; int *indx; char line[MAXLINE], linepar[MAXLINE]; - char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; + char linetmp[MAXLINE]; + char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; char pathr[MAXLINE], pathimach[MAXLINE]; char **bp, *tok, *val; /* pathtot */ int firstobs=1, lastobs=10; @@ -4635,10 +4657,11 @@ int main(int argc, char *argv[]) covar=matrix(0,NCOVMAX,1,n); cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; - - ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ + /* where is ncovprod ?*/ + ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ - npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ + nforces= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ + npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters like aij*/ 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); @@ -4858,6 +4881,11 @@ run imach with mle=-1 to get a correct t printf("Comment line\n%s\n",line); continue; } + trimbb(linetmp,line); /* Trims multiple blanks in line */ + for (j=0; line[j]!='\0';j++){ + line[j]=linetmp[j]; + } + for (j=maxwav;j>=1;j--){ cutv(stra, strb,line,' '); @@ -5040,7 +5068,7 @@ run imach with mle=-1 to get a correct t cptcovprod--; cutv(strb,stre,strd,'V'); Tvar[i]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3 */ - cptcovage++; /* Sum the number of covariates including ages as a product */ + cptcovage++; /* Sums the number of covariates including age as a product */ Tage[cptcovage]=i; /* Tage[1] =2 */ /*printf("stre=%s ", stre);*/ } @@ -5436,8 +5464,8 @@ Interval (in months) between two waves: } /* Endof if mle==-3 */ else{ /* For mle >=1 */ - globpr=1;/* debug */ - /* likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone);*/ /* Prints the contributions to the likelihood */ + globpr=0;/* 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]);