--- imach/src/imach.c 2009/07/06 08:22:05 1.132 +++ imach/src/imach.c 2010/04/26 20:30:53 1.136 @@ -1,6 +1,21 @@ -/* $Id: imach.c,v 1.132 2009/07/06 08:22:05 brouard Exp $ +/* $Id: imach.c,v 1.136 2010/04/26 20:30:53 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.136 2010/04/26 20:30:53 brouard + (Module): merging some libgsl code. Fixing computation + of likelione (using inter/intrapolation if mle = 0) in order to + get same likelihood as if mle=1. + Some cleaning of code and comments added. + + Revision 1.135 2009/10/29 15:33:14 brouard + (Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code. + + Revision 1.134 2009/10/29 13:18:53 brouard + (Module): Now imach stops if date of birth, at least year of birth, is not given. Some cleaning of the code. + + Revision 1.133 2009/07/06 10:21:25 brouard + just nforces + Revision 1.132 2009/07/06 08:22:05 brouard Many tings @@ -171,7 +186,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 @@ -304,8 +319,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 @@ -366,6 +381,11 @@ extern int errno; #include #include "timeval.h" +#ifdef GSL +#include +#include +#endif + /* #include */ /* #define _(String) gettext (String) */ @@ -400,15 +420,15 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.132 2009/07/06 08:22:05 brouard Exp $ */ +/* $Id: imach.c,v 1.136 2010/04/26 20:30:53 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98k, June 2006, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.132 $ $Date: 2009/07/06 08:22:05 $"; +char version[]="Imach version 0.98m, April 2010, INED-EUROREVES-Institut de longevite "; +char fullversion[]="$Revision: 1.136 $ $Date: 2010/04/26 20:30:53 $"; 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 */ @@ -430,7 +450,8 @@ int **bh; /* bh[mi][i] is the bias (+ or 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 *fic ; */ /* Used in readdata only */ +FILE *ficpar, *ficparo,*ficres, *ficresp, *ficrespl, *ficrespij, *ficrest,*ficresf,*ficrespop; FILE *ficlog, *ficrespow; int globpr=0; /* Global variable for printing or not */ double fretone; /* Only one call to likelihood */ @@ -1214,7 +1235,7 @@ double **prevalim(double **prlim, int nl } for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; for (k=1; k<=cptcovprod;k++) - cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; + cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ @@ -1703,8 +1724,9 @@ double funcone( double *x) lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ } else if (mle==4){ /* mle=4 no inter-extrapolation */ lli=log(out[s1][s2]); /* Original formula */ - } else{ /* ml>=5 no inter-extrapolation no jackson =0.8a */ - lli=log(out[s1][s2]); /* Original formula */ + } else{ /* mle=0 back to 1 */ + lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ + /*lli=log(out[s1][s2]); */ /* Original formula */ } /* End of if */ ipmx +=1; sw += weight[i]; @@ -2427,7 +2449,7 @@ void concatwav(int wav[], int **dh, int dh[mi][i]=jk; bh[mi][i]=0; }else{ /* We want a negative bias in order to only have interpolation ie - * at the price of an extra matrix product in likelihood */ + * to avoid the price of an extra matrix product in likelihood */ dh[mi][i]=jk+1; bh[mi][i]=ju; } @@ -2459,38 +2481,41 @@ void concatwav(int wav[], int **dh, int /*********** Tricode ****************************/ void tricode(int *Tvar, int **nbcode, int imx) { - + /* Uses cptcovn+2*cptcovprod as the number of covariates */ /* Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */ int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX; - int cptcode=0; + int modmaxcovj=0; /* Modality max of covariates j */ cptcoveff=0; for (k=0; k 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.*/ + if (ij > modmaxcovj) modmaxcovj=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 coded 0 and + female is 1, then modmaxcovj=1.*/ } - 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 */ + for (i=0; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*/ + if( Ndum[i] != 0 ) + ncodemax[j]++; + /* Number of modalities of the j th covariate. In fact + ncodemax[j]=2 (dichotom. variables only) but it could be more for + historical reasons */ } /* Ndum[-1] number of undefined modalities */ + /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ ij=1; - for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 */ - for (k=0; k<= maxncov; k++) { /* k=-1 ?*/ + for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 for dichotomous */ + for (k=0; k<= maxncov; k++) { /* k=-1 ? NCOVMAX*/ 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 @@ -3479,6 +3504,13 @@ To be simple, these graphs help to under /* Computing eigen value of matrix of covariance */ lc1=((v1+v2)+sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; lc2=((v1+v2)-sqrt((v1+v2)*(v1+v2) - 4*(v1*v2-cv12*cv12)))/2.; + if ((lc2 <0) || (lc1 <0) ){ + printf("Error: One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e.\n It means that the matrix was not well estimated (varpij), for i=%2d, j=%2d, age=%4d .\n See files %s and %s. Continuing by making them positive: WRONG RESULTS.\n", lc1, lc2, v1, v2, cv12, i, j, (int)age,fileresprobcov, fileresprobcor); + fprintf(ficlog,"Error: One eigen value of 2x2 matrix of covariance is negative, lc1=%11.3e, lc2=%11.3e, v1=%11.3e, v2=%11.3e, cv12=%11.3e\n", lc1, lc2, v1, v2, cv12);fflush(ficlog); + lc1=fabs(lc1); + lc2=fabs(lc2); + } + /* Eigen vectors */ v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12)); /*v21=sqrt(1.-v11*v11); *//* error */ @@ -4370,6 +4402,47 @@ double gompertz(double x[]) return -2*L*num/sump; } +#ifdef GSL +/******************* Gompertz_f Likelihood ******************************/ +double gompertz_f(const gsl_vector *v, void *params) +{ + double A,B,LL=0.0,sump=0.,num=0.; + double *x= (double *) v->data; + int i,n=0; /* n is the size of the sample */ + + for (i=0;i<=imx-1 ; i++) { + sump=sump+weight[i]; + /* sump=sump+1;*/ + num=num+1; + } + + + /* for (i=0; i<=imx; i++) + if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ + printf("x[0]=%lf x[1]=%lf\n",x[0],x[1]); + for (i=1;i<=imx ; i++) + { + if (cens[i] == 1 && wav[i]>1) + A=-x[0]/(x[1])*(exp(x[1]*(agecens[i]-agegomp))-exp(x[1]*(ageexmed[i]-agegomp))); + + if (cens[i] == 0 && wav[i]>1) + A=-x[0]/(x[1])*(exp(x[1]*(agedc[i]-agegomp))-exp(x[1]*(ageexmed[i]-agegomp))) + +log(x[0]/YEARM)+x[1]*(agedc[i]-agegomp)+log(YEARM); + + /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ + if (wav[i] > 1 ) { /* ??? */ + LL=LL+A*weight[i]; + /* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/ + } + } + + /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/ + printf("x[0]=%lf x[1]=%lf -2*LL*num/sump=%lf\n",x[0],x[1],-2*LL*num/sump); + + return -2*LL*num/sump; +} +#endif + /******************* Printing html file ***********/ void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \ int lastpass, int stepm, int weightopt, char model[],\ @@ -4417,90 +4490,475 @@ void printinggnuplotmort(char fileres[], } +int readdata(char datafile[], int firstobs, int lastobs, int *imax) +{ + /*-------- data file ----------*/ + FILE *fic; + char dummy[]=" "; + int i, j, n; + int linei, month, year,iout; + char line[MAXLINE], linetmp[MAXLINE]; + char stra[80], strb[80]; + char *stratrunc; + int lstra; + if((fic=fopen(datafile,"r"))==NULL) { + printf("Problem while opening datafile: %s\n", datafile);return 1; + fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1; + } -/***********************************************/ -/**************** Main Program *****************/ -/***********************************************/ + i=1; + linei=0; + while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { + linei=linei+1; + for(j=strlen(line); j>=0;j--){ /* Untabifies line */ + if(line[j] == '\t') + line[j] = ' '; + } + for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){ + ; + }; + line[j+1]=0; /* Trims blanks at end of line */ + if(line[0]=='#'){ + fprintf(ficlog,"Comment line\n%s\n",line); + 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]; + } + -int main(int argc, char *argv[]) -{ - int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); - int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod; - int linei, month, year,iout; - int jj, ll, li, lj, lk, imk; - int numlinepar=0; /* Current linenumber of parameter file */ - int itimes; - int NDIM=2; - int vpopbased=0; + for (j=maxwav;j>=1;j--){ + cutv(stra, strb,line,' '); + 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 '%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); + return 1; + } + } + s[j][i]=lval; + + strcpy(line,stra); + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,"%s.") != 0){ + month=99; + year=9999; + }else{ + 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); + return 1; + } + anint[j][i]= (double) year; + mint[j][i]= (double)month; + strcpy(line,stra); + } /* ENd Waves */ + + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,"%s.",dummy) != 0){ + month=99; + year=9999; + }else{ + 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); + return 1; + } + andc[i]=(double) year; + moisdc[i]=(double) month; + strcpy(line,stra); + + cutv(stra, strb,line,' '); + if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ + } + else if(iout=sscanf(strb,"%s.") != 0){ + month=99; + year=9999; + }else{ + printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (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 birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog); + return 1; + } + if (year==9999) { + printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. 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 birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog); + return 1; - char ca[32], cb[32], cc[32]; - char dummy[]=" "; - /* FILE *fichtm; *//* Html File */ - /* FILE *ficgp;*/ /*Gnuplot File */ - struct stat info; - double agedeb, agefin,hf; - double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; + } + annais[i]=(double)(year); + moisnais[i]=(double)(month); + strcpy(line,stra); + + cutv(stra, strb,line,' '); + errno=0; + 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); + 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); + return 1; + } + weight[i]=dval; + strcpy(line,stra); + + for (j=ncovcol;j>=1;j--){ + cutv(stra, strb,line,' '); + 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); + return 1; + } + } + if(lval <-1 || lval >1){ + printf("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); + 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); + return 1; + } + covar[j][i]=(double)(lval); + strcpy(line,stra); + } + lstra=strlen(stra); + + if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ + stratrunc = &(stra[lstra-9]); + num[i]=atol(stratrunc); + } + else + num[i]=atol(stra); + /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ + printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ + + i=i+1; + } /* End loop reading data */ - double fret; - double **xi,tmp,delta; + *imax=i-1; /* Number of individuals */ + fclose(fic); + + return (0); + endread: + printf("Exiting readdata: "); + fclose(fic); + return (1); - double dum; /* Dummy variable */ - double ***p3mat; - double ***mobaverage; - int *indx; - char line[MAXLINE], linepar[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; - int sdeb, sfin; /* Status at beginning and end */ - int c, h , cpt,l; - int ju,jl, mi; - int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij; - int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,*tab; - int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ - int mobilav=0,popforecast=0; - int hstepm, nhstepm; - int agemortsup; - float sumlpop=0.; - double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000; - double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000; - double bage, fage, age, agelim, agebase; - double ftolpl=FTOL; - double **prlim; - double *severity; - double ***param; /* Matrix of parameters */ - double *p; - double **matcov; /* Matrix of covariance */ - double ***delti3; /* Scale */ - double *delti; /* Scale */ - double ***eij, ***vareij; - double **varpl; /* Variances of prevalence limits by age */ - double *epj, vepp; - double kk1, kk2; - double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; - double **ximort; - char *alph[]={"a","a","b","c","d","e"}, str[4]; - int *dcwave; - char z[1]="c", occ; +} - char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80]; - char *strt, strtend[80]; - char *stratrunc; - int lstra; +int decodemodel ( char model[], int lastobs) +{ + int i, j, k; + int i1, j1, k1, k2; + char modelsav[80]; + char stra[80], strb[80], strc[80], strd[80],stre[80]; - long total_usecs; - -/* setlocale (LC_ALL, ""); */ -/* bindtextdomain (PACKAGE, LOCALEDIR); */ -/* textdomain (PACKAGE); */ -/* setlocale (LC_CTYPE, ""); */ + if (strlen(model) >1){ /* If there is at least 1 covariate */ + j=0, j1=0, k1=1, k2=1; + j=nbocc(model,'+'); /* j=Number of '+' */ + j1=nbocc(model,'*'); /* j1=Number of '*' */ + cptcovn=j+1; /* Number of covariates V1+V2*age+V3 =>(2 plus signs) + 1=3 + but the covariates which are product must be computed and stored. */ + cptcovprod=j1; /*Number of products V1*V2 +v3*age = 2 */ + + 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);fflush(ficlog); + return 1; + } + + /* 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 */ + /* modelsav=V3*age+V2+V1+V4 strb=V3*age stra=V2+V1+V4 + i=1 Tvar[1]=3 Tage[1]=1 + i=2 Tvar[2]=2 + i=3 Tvar[3]=1 + i=4 Tvar[4]= 4 + i=5 Tvar[5] + for (k=1; k<=cptcovn;k++) { + cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; + */ + for(k=1; k<=(j+1);k++){ + cutv(strb,stra,modelsav,'+'); /* keeps in strb after the first '+' + modelsav=V3*age+V2+V1+V4 strb=V3*age stra=V2+V1+V4 + */ + /* 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 V3*age+V2+V1+V4 strb=V3*age */ + cutv(strd,strc,strb,'*'); /* strd*strc Vm*Vn: strb=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'); /* stre="V3" */ + Tvar[k]=atoi(stre); /* V1+V3*age+V2 Tvar[2]=3, and Tvar[3]=2 */ + cptcovage++; /* Sums the number of covariates which include age as a product */ + Tage[cptcovage]=k; /* Tage[1] =2 */ + /*printf("stre=%s ", stre);*/ + } + else if (strcmp(strd,"age")==0) { /* or age*Vn */ + cptcovprod--; + cutv(strb,stre,strc,'V'); + Tvar[k]=atoi(stre); + cptcovage++; + Tage[cptcovage]=k; + } + 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[k]=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]=k; /* Tprod[1] */ + Tvard[k1][1]=atoi(strc); /* m*/ + Tvard[k1][2]=atoi(stre); /* n */ + Tvar[cptcovn+k2]=Tvard[k1][1]; + Tvar[cptcovn+k2+1]=Tvard[k1][2]; + for (i=1; i<=lastobs;i++) /* Computes the new covariate which is a product of covar[n][i]* covar[m][i] + and is stored at ncovol+k1 */ + covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; + k1++; + k2=k2+2; + } + } + else { /* no more sum */ + /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ + /* scanf("%d",i);*/ + cutv(strd,strc,strb,'V'); + Tvar[i]=atoi(strc); + } + 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 + */ + } /* end model */ + + /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products. + If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ + + /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); + printf("cptcovprod=%d ", cptcovprod); + fprintf(ficlog,"cptcovprod=%d ", cptcovprod); + + scanf("%d ",i);*/ + + + return (0); + endread: + printf("Exiting decodemodel: "); + return (1); +} + +calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn ) +{ + int i, m; + + for (i=1; i<=imx; i++) { + for(m=2; (m<= maxwav); m++) { + if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){ + anint[m][i]=9999; + s[m][i]=-1; + } + if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){ + *nberr++; + printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i); + fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i); + s[m][i]=-1; + } + if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){ + *nberr++; + printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); + fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); + s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */ + } + } + } + + for (i=1; i<=imx; i++) { + agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); + for(m=firstpass; (m<= lastpass); m++){ + if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ + if (s[m][i] >= nlstate+1) { + if(agedc[i]>0) + if((int)moisdc[i]!=99 && (int)andc[i]!=9999) + agev[m][i]=agedc[i]; + /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ + else { + if ((int)andc[i]!=9999){ + nbwarn++; + printf("Warning negative age at death: %ld line:%d\n",num[i],i); + fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i); + agev[m][i]=-1; + } + } + } + else if(s[m][i] !=9){ /* Standard case, age in fractional + years but with the precision of a month */ + agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); + if((int)mint[m][i]==99 || (int)anint[m][i]==9999) + agev[m][i]=1; + else if(agev[m][i] < *agemin){ + *agemin=agev[m][i]; + printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], *agemin); + } + else if(agev[m][i] >*agemax){ + *agemax=agev[m][i]; + /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/ + } + /*agev[m][i]=anint[m][i]-annais[i];*/ + /* agev[m][i] = age[i]+2*m;*/ + } + else { /* =9 */ + agev[m][i]=1; + s[m][i]=-1; + } + } + else /*= 0 Unknown */ + agev[m][i]=1; + } + + } + for (i=1; i<=imx; i++) { + for(m=firstpass; (m<=lastpass); m++){ + if (s[m][i] > (nlstate+ndeath)) { + *nberr++; + printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); + fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); + return 1; + } + } + } + + /*for (i=1; i<=imx; i++){ + for (m=firstpass; (m1) cptcovn=nbocc(model,'+')+1; - /* where is ncovprod ?*/ - ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/ + cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement plus one, indepently of n in Vn*/ + /* v1+v2+v3+v2*v4+v5*age makes cptcovn = 5 + v1+v2*age+v2*v3 makes cptcovn = 3 + */ + if (strlen(model)>1) + cptcovn=nbocc(model,'+')+1; + /* ncovprod */ + ncovmodel=2+cptcovn; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/ nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ - 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*/ + 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); @@ -4840,15 +5302,8 @@ run imach with mle=-1 to get a correct t fprintf(ficres,"#%s\n",version); } /* End of mle != -3 */ - /*-------- data file ----------*/ - if((fic=fopen(datafile,"r"))==NULL) { - printf("Problem while opening datafile: %s\n", datafile);goto end; - fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);goto end; - } n= lastobs; - severity = vector(1,maxwav); - outcome=imatrix(1,maxwav+1,1,n); num=lvector(1,n); moisnais=vector(1,n); annais=vector(1,n); @@ -4864,172 +5319,9 @@ run imach with mle=-1 to get a correct t tab=ivector(1,NCOVMAX); ncodemax=ivector(1,8); - i=1; - linei=0; - while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) { - linei=linei+1; - for(j=strlen(line); j>=0;j--){ /* Untabifies line */ - if(line[j] == '\t') - line[j] = ' '; - } - for(j=strlen(line)-1; (line[j]==' ')||(line[j]==10)||(line[j]==13);j--){ - ; - }; - line[j+1]=0; /* Trims blanks at end of line */ - if(line[0]=='#'){ - fprintf(ficlog,"Comment line\n%s\n",line); - 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,' '); - 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 '%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; - - strcpy(line,stra); - cutv(stra, strb,line,' '); - if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ - } - else if(iout=sscanf(strb,"%s.") != 0){ - month=99; - year=9999; - }else{ - 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; - strcpy(line,stra); - } /* ENd Waves */ - - cutv(stra, strb,line,' '); - if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ - } - else if(iout=sscanf(strb,"%s.",dummy) != 0){ - month=99; - year=9999; - }else{ - 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; - strcpy(line,stra); - - cutv(stra, strb,line,' '); - if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ - } - else if(iout=sscanf(strb,"%s.") != 0){ - 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 birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j); - 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); - strcpy(line,stra); - - cutv(stra, strb,line,' '); - errno=0; - 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); - 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,' '); - 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 \ - 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); - 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); - strcpy(line,stra); - } - lstra=strlen(stra); - - if(lstra > 9){ /* More than 2**32 or max of what printf can write with %ld */ - stratrunc = &(stra[lstra-9]); - num[i]=atol(stratrunc); - } - else - num[i]=atol(stra); - /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ - printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ - - i=i+1; - } /* End loop reading data */ - fclose(fic); - /* printf("ii=%d", ij); - scanf("%d",i);*/ - imx=i-1; /* Number of individuals */ - - /* for (i=1; i<=imx; i++){ - if ((s[1][i]==3) && (s[2][i]==2)) s[2][i]=3; - if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3; - if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3; - }*/ - /* for (i=1; i<=imx; i++){ - if (s[4][i]==9) s[4][i]=-1; - printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}*/ - - /* for (i=1; i<=imx; i++) */ - - /*if ((s[3][i]==3) || (s[4][i]==3)) weight[i]=0.08; - else weight[i]=1;*/ + /* Reads data from file datafile */ + if (readdata(datafile, firstobs, lastobs, &imx)==1) + goto end; /* Calculation of the number of parameters from char model */ 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 */ @@ -5037,182 +5329,23 @@ run imach with mle=-1 to get a correct t Tvaraff=ivector(1,15); Tvard=imatrix(1,15,1,2); Tage=ivector(1,15); - - if (strlen(model) >1){ /* If there is at least 1 covariate */ - j=0, j1=0, k1=1, k2=1; - j=nbocc(model,'+'); /* j=Number of '+' */ - j1=nbocc(model,'*'); /* j1=Number of '*' */ - 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);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 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 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); /* V1+V3*age+V2 Tvar[2]=3 */ - cptcovage++; /* Sums the number of covariates including age as a product */ - Tage[cptcovage]=i; /* Tage[1] =2 */ - /*printf("stre=%s ", stre);*/ - } - else if (strcmp(strd,"age")==0) { /* or age*Vn */ - cptcovprod--; - cutv(strb,stre,strc,'V'); - Tvar[i]=atoi(stre); - cptcovage++; - Tage[cptcovage]=i; - } - 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[1] */ - Tvard[k1][1]=atoi(strc); /* m*/ - Tvard[k1][2]=atoi(stre); /* n */ - Tvar[cptcovn+k2]=Tvard[k1][1]; - Tvar[cptcovn+k2+1]=Tvard[k1][2]; - for (k=1; k<=lastobs;k++) - covar[ncovcol+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k]; - k1++; - k2=k2+2; - } - } - else { /* no more sum */ - /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ - /* scanf("%d",i);*/ - cutv(strd,strc,strb,'V'); - Tvar[i]=atoi(strc); - } - 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 + */ - } /* end model */ - - /*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products. - If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ - /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); - printf("cptcovprod=%d ", cptcovprod); - fprintf(ficlog,"cptcovprod=%d ", cptcovprod); - - scanf("%d ",i);*/ + if(decodemodel(model, lastobs) == 1) + goto end; /* if(mle==1){*/ if (weightopt != 1) { /* Maximisation without weights*/ for(i=1;i<=n;i++) weight[i]=1.0; } + /*-calculation of age at interview from date of interview and age at death -*/ agev=matrix(1,maxwav,1,imx); - for (i=1; i<=imx; i++) { - for(m=2; (m<= maxwav); m++) { - if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){ - anint[m][i]=9999; - s[m][i]=-1; - } - if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){ - nberr++; - printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i); - fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i); - s[m][i]=-1; - } - if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){ - nberr++; - printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); - fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); - s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */ - } - } - } - - for (i=1; i<=imx; i++) { - agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); - for(m=firstpass; (m<= lastpass); m++){ - if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ - if (s[m][i] >= nlstate+1) { - if(agedc[i]>0) - if((int)moisdc[i]!=99 && (int)andc[i]!=9999) - agev[m][i]=agedc[i]; - /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ - else { - if ((int)andc[i]!=9999){ - nbwarn++; - printf("Warning negative age at death: %ld line:%d\n",num[i],i); - fprintf(ficlog,"Warning negative age at death: %ld line:%d\n",num[i],i); - agev[m][i]=-1; - } - } - } - else if(s[m][i] !=9){ /* Standard case, age in fractional - years but with the precision of a month */ - agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); - if((int)mint[m][i]==99 || (int)anint[m][i]==9999) - agev[m][i]=1; - else if(agev[m][i] agemax){ - agemax=agev[m][i]; - /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/ - } - /*agev[m][i]=anint[m][i]-annais[i];*/ - /* agev[m][i] = age[i]+2*m;*/ - } - else { /* =9 */ - agev[m][i]=1; - s[m][i]=-1; - } - } - else /*= 0 Unknown */ - agev[m][i]=1; - } - - } - for (i=1; i<=imx; i++) { - for(m=firstpass; (m<=lastpass); m++){ - if (s[m][i] > (nlstate+ndeath)) { - nberr++; - printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); - fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); - goto end; - } - } - } - - /*for (i=1; i<=imx; i++){ - for (m=firstpass; (mm) - h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; + 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]]); } } @@ -5344,7 +5480,8 @@ Interval (in months) between two waves: globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ if (mle==-3){ - ximort=matrix(1,NDIM,1,NDIM); + ximort=matrix(1,NDIM,1,NDIM); +/* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ cens=ivector(1,n); ageexmed=vector(1,n); agecens=vector(1,n); @@ -5386,21 +5523,106 @@ Interval (in months) between two waves: /*printf("%lf %lf", p[1], p[2]);*/ +#ifdef GSL + printf("GSL optimization\n"); fprintf(ficlog,"Powell\n"); +#elsedef printf("Powell\n"); fprintf(ficlog,"Powell\n"); +#endif strcpy(filerespow,"pow-mort"); strcat(filerespow,fileres); if((ficrespow=fopen(filerespow,"w"))==NULL) { printf("Problem with resultfile: %s\n", filerespow); fprintf(ficlog,"Problem with resultfile: %s\n", filerespow); } +#ifdef GSL + fprintf(ficrespow,"# GSL optimization\n# iter -2*LL"); +#elsedef fprintf(ficrespow,"# Powell\n# iter -2*LL"); +#endif /* for (i=1;i<=nlstate;i++) for(j=1;j<=nlstate+ndeath;j++) if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); */ fprintf(ficrespow,"\n"); +#ifdef GSL + /* gsl starts here */ + T = gsl_multimin_fminimizer_nmsimplex; + gsl_multimin_fminimizer *sfm = NULL; + gsl_vector *ss, *x; + gsl_multimin_function minex_func; + + /* Initial vertex size vector */ + ss = gsl_vector_alloc (NDIM); + + if (ss == NULL){ + GSL_ERROR_VAL ("failed to allocate space for ss", GSL_ENOMEM, 0); + } + /* Set all step sizes to 1 */ + gsl_vector_set_all (ss, 0.001); + + /* Starting point */ + + x = gsl_vector_alloc (NDIM); - powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); + if (x == NULL){ + gsl_vector_free(ss); + GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0); + } + + /* Initialize method and iterate */ + /* p[1]=0.0268; p[NDIM]=0.083; */ +/* gsl_vector_set(x, 0, 0.0268); */ +/* gsl_vector_set(x, 1, 0.083); */ + gsl_vector_set(x, 0, p[1]); + gsl_vector_set(x, 1, p[2]); + + minex_func.f = &gompertz_f; + minex_func.n = NDIM; + minex_func.params = (void *)&p; /* ??? */ + + sfm = gsl_multimin_fminimizer_alloc (T, NDIM); + gsl_multimin_fminimizer_set (sfm, &minex_func, x, ss); + + printf("Iterations beginning .....\n\n"); + printf("Iter. # Intercept Slope -Log Likelihood Simplex size\n"); + + iteri=0; + while (rval == GSL_CONTINUE){ + iteri++; + status = gsl_multimin_fminimizer_iterate(sfm); + + if (status) printf("error: %s\n", gsl_strerror (status)); + fflush(0); + + if (status) + break; + + rval = gsl_multimin_test_size (gsl_multimin_fminimizer_size (sfm), 1e-6); + ssval = gsl_multimin_fminimizer_size (sfm); + + if (rval == GSL_SUCCESS) + printf ("converged to a local maximum at\n"); + + printf("%5d ", iteri); + for (it = 0; it < NDIM; it++){ + printf ("%10.5f ", gsl_vector_get (sfm->x, it)); + } + printf("f() = %-10.5f ssize = %.7f\n", sfm->fval, ssval); + } + + printf("\n\n Please note: Program should be run many times with varying starting points to detemine global maximum\n\n"); + + gsl_vector_free(x); /* initial values */ + gsl_vector_free(ss); /* inital step size */ + for (it=0; itx,it); + fprintf(ficrespow," %.12lf", p[it]); + } + gsl_multimin_fminimizer_free (sfm); /* p *(sfm.x.data) et p *(sfm.x.data+1) */ +#endif +#ifdef POWELL + powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); +#endif fclose(ficrespow); hesscov(matcov, p, NDIM, delti, 1e-4, gompertz); @@ -5461,6 +5683,12 @@ Interval (in months) between two waves: free_vector(lsurv,1,AGESUP); free_vector(lpop,1,AGESUP); free_vector(tpop,1,AGESUP); +#ifdef GSL + free_ivector(cens,1,n); + free_vector(agecens,1,n); + free_ivector(dcwave,1,n); + free_matrix(ximort,1,NDIM,1,NDIM); +#endif } /* Endof if mle==-3 */ else{ /* For mle >=1 */