--- imach/src/imach.c 2004/05/16 15:05:56 1.98 +++ imach/src/imach.c 2005/09/30 16:11:43 1.104 @@ -1,6 +1,30 @@ -/* $Id: imach.c,v 1.98 2004/05/16 15:05:56 brouard Exp $ +/* $Id: imach.c,v 1.104 2005/09/30 16:11:43 lievre Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.104 2005/09/30 16:11:43 lievre + (Module): sump fixed, loop imx fixed, and simplifications. + (Module): If the status is missing at the last wave but we know + that the person is alive, then we can code his/her status as -2 + (instead of missing=-1 in earlier versions) and his/her + contributions to the likelihood is 1 - Prob of dying from last + health status (= 1-p13= p11+p12 in the easiest case of somebody in + the healthy state at last known wave). Version is 0.98 + + Revision 1.103 2005/09/30 15:54:49 lievre + (Module): sump fixed, loop imx fixed, and simplifications. + + Revision 1.102 2004/09/15 17:31:30 brouard + Add the possibility to read data file including tab characters. + + Revision 1.101 2004/09/15 10:38:38 brouard + Fix on curr_time + + Revision 1.100 2004/07/12 18:29:06 brouard + Add version for Mac OS X. Just define UNIX in Makefile + + Revision 1.99 2004/06/05 08:57:40 brouard + *** empty log message *** + Revision 1.98 2004/05/16 15:05:56 brouard New version 0.97 . First attempt to estimate force of mortality directly from the data i.e. without the need of knowing the health @@ -196,7 +220,7 @@ #include #include -#include +/* #include */ #include #include "timeval.h" @@ -224,7 +248,7 @@ #define AGESUP 130 #define AGEBASE 40 #define AGEGOMP 10. /* Minimal age for Gompertz adjustment */ -#ifdef unix +#ifdef UNIX #define DIRSEPARATOR '/' #define ODIRSEPARATOR '\\' #else @@ -232,11 +256,11 @@ #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.98 2004/05/16 15:05:56 brouard Exp $ */ +/* $Id: imach.c,v 1.104 2005/09/30 16:11:43 lievre Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.70, May 2004, INED-EUROREVES "; -char fullversion[]="$Revision: 1.98 $ $Date: 2004/05/16 15:05:56 $"; +char version[]="Imach version 0.98, September 2005, INED-EUROREVES "; +char fullversion[]="$Revision: 1.104 $ $Date: 2005/09/30 16:11:43 $"; int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; @@ -356,6 +380,9 @@ double ftolhess; /* Tolerance for comput /**************** split *************************/ static int split( char *path, char *dirc, char *name, char *ext, char *finame ) { + /* From a file name with full path (either Unix or Windows) we extract the directory (dirc) + the name of the file (name), its extension only (ext) and its first part of the name (finame) + */ char *ss; /* pointer */ int l1, l2; /* length counters */ @@ -387,12 +414,14 @@ static int split( char *path, char *dirc #endif */ ss = strrchr( name, '.' ); /* find last / */ - ss++; - strcpy(ext,ss); /* save extension */ - l1= strlen( name); - l2= strlen(ss)+1; - strncpy( finame, name, l1-l2); - finame[l1-l2]= 0; + if (ss >0){ + ss++; + strcpy(ext,ss); /* save extension */ + l1= strlen( name); + l2= strlen(ss)+1; + strncpy( finame, name, l1-l2); + finame[l1-l2]= 0; + } return( 0 ); /* we're done */ } @@ -425,8 +454,8 @@ int nbocc(char *s, char occ) void cutv(char *u,char *v, char*t, char occ) { - /* cuts string t into u and v where u is ended by char occ excluding it - and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2) + /* cuts string t into u and v where u ends before first occurence of char 'occ' + and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2') gives u="abcedf" and v="ghi2j" */ int i,lg,j,p=0; i=0; @@ -868,11 +897,11 @@ void powell(double p[], double **xi, int fprintf(ficrespow,"\n");fflush(ficrespow); if(*iter <=3){ tm = *localtime(&curr_time.tv_sec); - strcpy(strcurr,asctime(&tmf)); + strcpy(strcurr,asctime(&tm)); /* asctime_r(&tm,strcurr); */ - forecast_time=curr_time; + forecast_time=curr_time; itmp = strlen(strcurr); - if(strcurr[itmp-1]=='\n') + if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ strcurr[itmp-1]='\0'; printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec); @@ -884,8 +913,8 @@ void powell(double p[], double **xi, int itmp = strlen(strfor); if(strfor[itmp-1]=='\n') strfor[itmp-1]='\0'; - printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); - fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); + printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); + fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr); } } for (i=1;i<=n;i++) { @@ -1050,57 +1079,57 @@ double **pmij(double **ps, double *cov, int i,j,j1, nc, ii, jj; for(i=1; i<= nlstate; i++){ - for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/ + for(j=1; ji s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */ + } + ps[i][j]=s2; } - ps[i][j]=s2; } - } /*ps[3][2]=1;*/ - - for(i=1; i<= nlstate; i++){ - s1=0; - for(j=1; j 1 the results are less biased than in previous versions. @@ -1239,15 +1268,16 @@ double func( double *x) s1=s[mw[mi][i]][i]; s2=s[mw[mi+1][i]][i]; bbh=(double)bh[mi][i]/(double)stepm; - /* bias is positive if real duration + /* bias bh is positive if real duration * is higher than the multiple of stepm and negative otherwise. */ /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ if( s2 > nlstate){ - /* i.e. if s2 is a death state and if the date of death is known then the contribution - to the likelihood is the probability to die between last step unit time and current - step unit time, which is also the differences between probability to die before dh - and probability to die before dh-stepm . + /* i.e. if s2 is a death state and if the date of death is known + then the contribution to the likelihood is the probability to + die between last step unit time and current step unit time, + which is also equal to probability to die before dh + minus probability to die before dh-stepm . In version up to 0.92 likelihood was computed as if date of death was unknown. Death was treated as any other health state: the date of the interview describes the actual state @@ -1267,7 +1297,16 @@ double func( double *x) lower mortality. */ lli=log(out[s1][s2] - savm[s1][s2]); - }else{ + + + } else if (s2==-2) { + for (j=1,survp=0. ; j<=nlstate; j++) + survp += out[s1][j]; + lli= survp; + } + + + else{ lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ /* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ } @@ -1851,7 +1890,7 @@ void freqsummary(char fileres[], int ia fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp); exit(0); } - freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3); + freq= ma3x(-2,nlstate+ndeath,-2,nlstate+ndeath,iagemin,iagemax+3); j1=0; j=cptcoveff; @@ -1864,8 +1903,8 @@ void freqsummary(char fileres[], int ia j1++; /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); scanf("%d", i);*/ - for (i=-1; i<=nlstate+ndeath; i++) - for (jk=-1; jk<=nlstate+ndeath; jk++) + for (i=-2; i<=nlstate+ndeath; i++) + for (jk=-2; jk<=nlstate+ndeath; jk++) for(m=iagemin; m <= iagemax+3; m++) freq[i][jk][m]=0; @@ -1990,7 +2029,7 @@ void freqsummary(char fileres[], int ia dateintmean=dateintsum/k2cpt; fclose(ficresp); - free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3); + free_ma3x(freq,-2,nlstate+ndeath,-2,nlstate+ndeath, iagemin, iagemax+3); free_vector(pp,1,nlstate); free_matrix(prop,1,nlstate,iagemin, iagemax+3); /* End of Freq */ @@ -2099,7 +2138,7 @@ void concatwav(int wav[], int **dh, int mi=0; m=firstpass; while(s[m][i] <= nlstate){ - if(s[m][i]>=1) + if(s[m][i]>=1 || s[m][i]==-2) mw[++mi][i]=m; if(m >=lastpass) break; @@ -2153,7 +2192,8 @@ void concatwav(int wav[], int **dh, int } else{ j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12)); - /* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/ +/* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */ + k=k+1; if (j >= jmax) jmax=j; else if (j <= jmin)jmin=j; @@ -2246,7 +2286,7 @@ void tricode(int *Tvar, int **nbcode, in for (k=0; k< maxncov; k++) Ndum[k]=0; for (i=1; i<=ncovmodel-2; i++) { - /* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ + /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ ij=Tvar[i]; Ndum[ij]++; } @@ -3904,24 +3944,22 @@ double gompertz(double x[]) 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; + /* sump=sump+1;*/ num=num+1; } - /* for (i=1; i<=imx; i++) + /* 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]);*/ - for (i=0;i<=imx-1 ; i++) + for (i=1;i<=imx ; i++) { if (cens[i]==1 & wav[i]>1) - A=-x[1]/(x[2])* - (exp(x[2]/YEARM*(agecens[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))); + A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); if (cens[i]==0 & wav[i]>1) - A=-x[1]/(x[2])* - (exp(x[2]/YEARM*(agedc[i]*12-agegomp*12))-exp(x[2]/YEARM*(ageexmed[i]*12-agegomp*12))) - +log(x[1]/YEARM)+x[2]/YEARM*(agedc[i]*12-agegomp*12)+log(YEARM); + A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) + +log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM); if (wav[i]>1 & agecens[i]>15) { L=L+A*weight[i]; @@ -4003,7 +4041,7 @@ int main(int argc, char *argv[]) int *indx; char line[MAXLINE], linepar[MAXLINE]; char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; - char pathr[MAXLINE]; + char pathr[MAXLINE], pathimach[MAXLINE]; int firstobs=1, lastobs=10; int sdeb, sfin; /* Status at beginning and end */ int c, h , cpt,l; @@ -4089,8 +4127,10 @@ int main(int argc, char *argv[]) printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/ /* cutv(path,optionfile,pathtot,'\\');*/ + split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname); + /* strcpy(pathimach,argv[0]); */ split(pathtot,path,optionfile,optionfilext,optionfilefiname); - printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); + printf("pathimach=%s, pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); chdir(path); strcpy(command,"mkdir "); strcat(command,optionfilefiname); @@ -4115,12 +4155,12 @@ int main(int argc, char *argv[]) } fprintf(ficlog,"Log filename:%s\n",filelog); fprintf(ficlog,"\n%s\n%s",version,fullversion); - fprintf(ficlog,"\nEnter the parameter file name: "); - fprintf(ficlog,"pathtot=%s\n\ + fprintf(ficlog,"\nEnter the parameter file name: \n"); + fprintf(ficlog,"pathimach=%s\npathtot=%s\n\ path=%s \n\ optionfile=%s\n\ optionfilext=%s\n\ - optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); + optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); printf("Local time (at start):%s",strstart); fprintf(ficlog,"Local time (at start): %s",strstart); @@ -4383,7 +4423,10 @@ int main(int argc, char *argv[]) i=1; while (fgets(line, MAXLINE, fic) != NULL) { if ((i >= firstobs) && (i <=lastobs)) { - + for(j=0; line[j] != '\n';j++){ /* Untabifies line */ + if(line[j] == '\t') + line[j] = ' '; + } for (j=maxwav;j>=1;j--){ cutv(stra, strb,line,' '); s[j][i]=atoi(strb); strcpy(line,stra); @@ -4546,7 +4589,7 @@ int main(int argc, char *argv[]) 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){ + if(s[m][i] >0 || s[m][i]==-2){ if (s[m][i] >= nlstate+1) { if(agedc[i]>0) if((int)moisdc[i]!=99 && (int)andc[i]!=9999) @@ -5411,7 +5454,10 @@ Interval (in months) between two waves: /*------ End -----------*/ chdir(path); - strcpy(plotcmd,GNUPLOTPROGRAM); + strcpy(plotcmd,"\""); + strcat(plotcmd,pathimach); + strcat(plotcmd,GNUPLOTPROGRAM); + strcat(plotcmd,"\""); strcat(plotcmd," "); strcat(plotcmd,optionfilegnuplot); printf("Starting graphs with: %s",plotcmd);fflush(stdout); @@ -5424,7 +5470,10 @@ Interval (in months) between two waves: printf("\nType e to edit output files, g to graph again and q for exiting: "); scanf("%s",z); /* if (z[0] == 'c') system("./imach"); */ - if (z[0] == 'e') system(optionfilehtm); + if (z[0] == 'e') { + printf("Starting browser with: %s",optionfilehtm);fflush(stdout); + system(optionfilehtm); + } else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); }