--- imach/src/imach.c 2009/06/20 16:22:47 1.131 +++ imach/src/imach.c 2009/07/06 10:21:25 1.133 @@ -1,6 +1,12 @@ -/* $Id: imach.c,v 1.131 2009/06/20 16:22:47 brouard Exp $ +/* $Id: imach.c,v 1.133 2009/07/06 10:21:25 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.133 2009/07/06 10:21:25 brouard + just nforces + + 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 @@ -168,7 +174,7 @@ The same imach parameter file can be used but the option for mle should be -3. - Agnès, who wrote this part of the code, tried to keep most of the + Agnès, who wrote this part of the code, tried to keep most of the former routines in order to include the new code within the former code. The output is very simple: only an estimate of the intercept and of @@ -301,8 +307,8 @@ Also this programme outputs the covariance matrix of the parameters but also of the life expectancies. It also computes the period (stable) prevalence. - Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). - Institut national d'études démographiques, Paris. + Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr). + Institut national d'études démographiques, Paris. This software have been partly granted by Euro-REVES, a concerted action from the European Union. It is copyrighted identically to a GNU software product, ie programme and @@ -397,15 +403,15 @@ 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.133 2009/07/06 10:21:25 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 version[]="Imach version 0.98k, June 2009, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.133 $ $Date: 2009/07/06 10:21:25 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ -int nvar=0; +int nvar=0, nforce=0; /* Number of variables, number of forces */ 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 */ @@ -597,6 +603,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 +1712,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 +1917,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 +1938,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 +1964,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 +2175,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 +2484,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 +4455,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 +4660,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*/ + nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ + npar= nforce*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 +4884,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 +5071,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 +5467,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]);