--- imach/src/imach.c 2009/10/29 15:33:14 1.135 +++ imach/src/imach.c 2010/04/26 20:30:53 1.136 @@ -1,6 +1,12 @@ -/* $Id: imach.c,v 1.135 2009/10/29 15:33:14 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. @@ -375,6 +381,11 @@ extern int errno; #include #include "timeval.h" +#ifdef GSL +#include +#include +#endif + /* #include */ /* #define _(String) gettext (String) */ @@ -409,11 +420,11 @@ extern int errno; #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.135 2009/10/29 15:33:14 brouard Exp $ */ +/* $Id: imach.c,v 1.136 2010/04/26 20:30:53 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98l, October 2009, INED-EUROREVES-Institut de longevite "; -char fullversion[]="$Revision: 1.135 $ $Date: 2009/10/29 15:33:14 $"; +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 */ @@ -439,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 */ @@ -1223,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]);*/ @@ -1712,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]; @@ -2436,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; } @@ -2468,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 @@ -4386,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[],\ @@ -4433,89 +4490,474 @@ 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); */ + 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 */ nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ npar= nforce*ncovmodel; /* Number of parameters like aij*/ @@ -4856,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); @@ -4880,178 +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 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); - goto end; - } - 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); - 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 */ @@ -5059,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]]); } } @@ -5366,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); @@ -5408,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); + + 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; /* ??? */ - powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); + 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); @@ -5483,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 */