--- imach/src/imach.c 2003/06/24 12:34:15 1.90 +++ imach/src/imach.c 2003/07/08 07:54:34 1.95 @@ -1,6 +1,30 @@ -/* $Id: imach.c,v 1.90 2003/06/24 12:34:15 brouard Exp $ +/* $Id: imach.c,v 1.95 2003/07/08 07:54:34 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.95 2003/07/08 07:54:34 brouard + * imach.c (Repository): + (Repository): Using imachwizard code to output a more meaningful covariance + matrix (cov(a12,c31) instead of numbers. + + Revision 1.94 2003/06/27 13:00:02 brouard + Just cleaning + + Revision 1.93 2003/06/25 16:33:55 brouard + (Module): On windows (cygwin) function asctime_r doesn't + exist so I changed back to asctime which exists. + (Module): Version 0.96b + + Revision 1.92 2003/06/25 16:30:45 brouard + (Module): On windows (cygwin) function asctime_r doesn't + exist so I changed back to asctime which exists. + + Revision 1.91 2003/06/25 15:30:29 brouard + * imach.c (Repository): Duplicated warning errors corrected. + (Repository): Elapsed time after each iteration is now output. It + helps to forecast when convergence will be reached. Elapsed time + is stamped in powell. We created a new html file for the graphs + concerning matrix of covariance. It has extension -cov.htm. + Revision 1.90 2003/06/24 12:34:15 brouard (Module): Some bugs corrected for windows. Also, when mle=-1 a template is output in file "or"mypar.txt with the design @@ -146,6 +170,9 @@ #include #include "timeval.h" +/* #include */ +/* #define _(String) gettext (String) */ + #define MAXLINE 256 #define GNUPLOTPROGRAM "gnuplot" /*#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"*/ @@ -174,12 +201,12 @@ #define ODIRSEPARATOR '/' #endif -/* $Id: imach.c,v 1.90 2003/06/24 12:34:15 brouard Exp $ */ +/* $Id: imach.c,v 1.95 2003/07/08 07:54:34 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.96a, June 2003, INED-EUROREVES "; -char fullversion[]="$Revision: 1.90 $ $Date: 2003/06/24 12:34:15 $"; -int erreur; /* Error number */ +char version[]="Imach version 0.96b, June 2003, INED-EUROREVES "; +char fullversion[]="$Revision: 1.95 $ $Date: 2003/07/08 07:54:34 $"; +int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ int nvar; int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; int npar=NPARMAX; @@ -211,7 +238,7 @@ char fileresilk[FILENAMELENGTH]; /* File FILE *ficresilk; FILE *ficgp,*ficresprob,*ficpop, *ficresprobcov, *ficresprobcor; FILE *ficresprobmorprev; -FILE *fichtm; /* Html File */ +FILE *fichtm, *fichtmcov; /* Html File */ FILE *ficreseij; char filerese[FILENAMELENGTH]; FILE *ficresvij; @@ -226,13 +253,21 @@ char command[FILENAMELENGTH]; int outcmd=0; char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH]; -char lfileres[FILENAMELENGTH]; + char filelog[FILENAMELENGTH]; /* Log file */ char filerest[FILENAMELENGTH]; char fileregp[FILENAMELENGTH]; char popfile[FILENAMELENGTH]; -char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH]; +char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilehtmcov[FILENAMELENGTH] ; + +struct timeval start_time, end_time, curr_time, last_time, forecast_time; +struct timezone tzp; +extern int gettimeofday(); +struct tm tmg, tm, tmf, *gmtime(), *localtime(); +long time_value; +extern long time(); +char strcurr[80], strfor[80]; #define NR_END 1 #define FREE_ARG char* @@ -538,6 +573,41 @@ void free_ma3x(double ***m, long nrl, lo free((FREE_ARG)(m+nrl-NR_END)); } +/*************** function subdirf ***********/ +char *subdirf(char fileres[]) +{ + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); /* Add to the right */ + strcat(tmpout,fileres); + return tmpout; +} + +/*************** function subdirf2 ***********/ +char *subdirf2(char fileres[], char *preop) +{ + + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); + strcat(tmpout,preop); + strcat(tmpout,fileres); + return tmpout; +} + +/*************** function subdirf3 ***********/ +char *subdirf3(char fileres[], char *preop, char *preop2) +{ + + /* Caution optionfilefiname is hidden */ + strcpy(tmpout,optionfilefiname); + strcat(tmpout,"/"); + strcat(tmpout,preop); + strcat(tmpout,preop2); + strcat(tmpout,fileres); + return tmpout; +} + /***************** f1dim *************************/ extern int ncom; extern double *pcom,*xicom; @@ -713,6 +783,19 @@ void linmin(double p[], double xi[], int free_vector(pcom,1,n); } +char *asc_diff_time(long time_sec, char ascdiff[]) +{ + long sec_left, days, hours, minutes; + days = (time_sec) / (60*60*24); + sec_left = (time_sec) % (60*60*24); + hours = (sec_left) / (60*60) ; + sec_left = (sec_left) %(60*60); + minutes = (sec_left) /60; + sec_left = (sec_left) % (60); + sprintf(ascdiff,"%d day(s) %d hour(s) %d minute(s) %d second(s)",days, hours, minutes, sec_left); + return ascdiff; +} + /*************** powell ************************/ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, double (*func)(double [])) @@ -723,6 +806,8 @@ void powell(double p[], double **xi, int double del,t,*pt,*ptt,*xit; double fp,fptt; double *xits; + int niterf, itmp; + pt=vector(1,n); ptt=vector(1,n); xit=vector(1,n); @@ -733,9 +818,11 @@ void powell(double p[], double **xi, int fp=(*fret); ibig=0; del=0.0; - printf("\nPowell iter=%d -2*LL=%.12f",*iter,*fret); - fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f",*iter,*fret); - fprintf(ficrespow,"%d %.12f",*iter,*fret); + last_time=curr_time; + (void) gettimeofday(&curr_time,&tzp); + printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout); + fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); + fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); @@ -743,7 +830,29 @@ void powell(double p[], double **xi, int } printf("\n"); fprintf(ficlog,"\n"); - fprintf(ficrespow,"\n"); + fprintf(ficrespow,"\n");fflush(ficrespow); + if(*iter <=3){ + tm = *localtime(&curr_time.tv_sec); + strcpy(strcurr,asctime(&tmf)); +/* asctime_r(&tm,strcurr); */ + forecast_time=curr_time; + itmp = strlen(strcurr); + if(strcurr[itmp-1]=='\n') + 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); + for(niterf=10;niterf<=30;niterf+=10){ + forecast_time.tv_sec=curr_time.tv_sec+(niterf-*iter)*(curr_time.tv_sec-last_time.tv_sec); + tmf = *localtime(&forecast_time.tv_sec); +/* asctime_r(&tmf,strfor); */ + strcpy(strfor,asctime(&tmf)); + 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); + } + } for (i=1;i<=n;i++) { for (j=1;j<=n;j++) xit[j]=xi[j][i]; fptt=(*fret); @@ -1156,31 +1265,10 @@ double func( double *x) oldm=newm; } /* end mult */ - /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ - /* But now since version 0.9 we anticipate for bias and large stepm. - * If stepm is larger than one month (smallest stepm) and if the exact delay - * (in months) between two waves is not a multiple of stepm, we rounded to - * the nearest (and in case of equal distance, to the lowest) interval but now - * we keep into memory the bias bh[mi][i] and also the previous matrix product - * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the - * probability in order to take into account the bias as a fraction of the way - * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies - * -stepm/2 to stepm/2 . - * For stepm=1 the results are the same as for previous versions of Imach. - * For stepm > 1 the results are less biased than in previous versions. - */ 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 - * is higher than the multiple of stepm and negative otherwise. - */ 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 */ - /* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ - /*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-+bh)*out[s1][s2])); */ /* exponential interpolation */ - /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ - /*if(lli ==000.0)*/ - /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; @@ -1207,30 +1295,10 @@ double func( double *x) oldm=newm; } /* end mult */ - /*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ - /* But now since version 0.9 we anticipate for bias and large stepm. - * If stepm is larger than one month (smallest stepm) and if the exact delay - * (in months) between two waves is not a multiple of stepm, we rounded to - * the nearest (and in case of equal distance, to the lowest) interval but now - * we keep into memory the bias bh[mi][i] and also the previous matrix product - * (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the - * probability in order to take into account the bias as a fraction of the way - * from savm to out if bh is neagtive or even beyond if bh is positive. bh varies - * -stepm/2 to stepm/2 . - * For stepm=1 the results are the same as for previous versions of Imach. - * For stepm > 1 the results are less biased than in previous versions. - */ 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 - * is higher than the multiple of stepm and negative otherwise. - */ - /* 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 */ lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ - /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ - /*if(lli ==000.0)*/ - /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ ipmx +=1; sw += weight[i]; ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; @@ -1396,35 +1464,8 @@ double funcone( double *x) return -l; } -char *subdirf(char fileres[]) -{ - - strcpy(tmpout,optionfilefiname); - strcat(tmpout,"/"); /* Add to the right */ - strcat(tmpout,fileres); - return tmpout; -} - -char *subdirf2(char fileres[], char *preop) -{ - - strcpy(tmpout,optionfilefiname); - strcat(tmpout,"/"); - strcat(tmpout,preop); - strcat(tmpout,fileres); - return tmpout; -} -char *subdirf3(char fileres[], char *preop, char *preop2) -{ - - strcpy(tmpout,optionfilefiname); - strcat(tmpout,"/"); - strcat(tmpout,preop); - strcat(tmpout,preop2); - strcat(tmpout,fileres); - return tmpout; -} +/*************** function likelione ***********/ void likelione(FILE *ficres,double p[], int npar, int nlstate, int *globpri, long *ipmx, double *sw, double *fretone, double (*funcone)(double [])) { /* This routine should help understanding what is done with @@ -2036,6 +2077,7 @@ void concatwav(int wav[], int **dh, int wav[i]=mi; if(mi==0){ + nbwarn++; if(first==0){ printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i); first=1; @@ -2056,10 +2098,11 @@ void concatwav(int wav[], int **dh, int j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); if(j==0) j=1; /* Survives at least one month after exam */ else if(j<0){ + nberr++; printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); - j=1; /* Careful Patch */ + j=1; /* Temporary Dangerous patch */ printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n You MUST fix the contradiction between dates.\n",stepm); - printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); + fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n You MUST fix the contradiction between dates.\n",stepm); } k=k+1; @@ -2079,6 +2122,7 @@ void concatwav(int wav[], int **dh, int /* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */ /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/ if(j<0){ + nberr++; printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); } @@ -2255,7 +2299,7 @@ void evsij(char fileres[], double ***eij hf=hstepm*stepm/YEARM; /* Duration of hstepm expressed in year unit. */ - /* Computing Variances of health expectancies */ + /* Computing Variances of health expectancies */ for(theta=1; theta <=npar; theta++){ for(i=1; i<=npar; i++){ @@ -2809,9 +2853,19 @@ void varprob(char optionfilefiname[], do fprintf(fichtm,"\n

  • Computing and drawing one step probabilities with their confidence intervals

    • \n"); fprintf(fichtm,"\n"); - fprintf(fichtm,"\n

  • Computing matrix of variance-covariance of step probabilities

    • \n"); - fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the pij, pkl to understand the covariance between two incidences. They are expressed in year-1 in order to be less dependent of stepm.
    \n"); - fprintf(fichtm,"\n
    We have drawn x'cov-1x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis.
    When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.
    \n"); + fprintf(fichtm,"\n

  • Matrix of variance-covariance of pairs of step probabilities (drawings)

    • \n",optionfilehtmcov); + fprintf(fichtmcov,"\n

    Matrix of variance-covariance of pairs of step probabilities

    \n\ + file %s
    \n",optionfilehtmcov); + fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (pij, pkl) are estimated\ +and drawn. It helps understanding how is the covariance between two incidences.\ + They are expressed in year-1 in order to be less dependent of stepm.
    \n"); + fprintf(fichtmcov,"\n
    Contour plot corresponding to x'cov-1x = 4 (where x is the column vector (pij,pkl)) are drawn. \ +It can be understood this way: if pij and pkl where uncorrelated the (2x2) matrix of covariance \ +would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 \ +standard deviations wide on each axis.
    \ + Now, if both incidences are correlated (usual case) we diagonalised the inverse of the covariance matrix\ + and made the appropriate rotation to look at the uncorrelated principal directions.
    \ +To be simple, these graphs help to understand the significativity of each parameter in relation to a second other one.
    \n"); cov[1]=1; tj=cptcoveff; @@ -2993,13 +3047,13 @@ void varprob(char optionfilefiname[], do fprintf(ficgp,"\nset parametric;unset label"); fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); - fprintf(fichtm,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ + fprintf(fichtmcov,"\n
    Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year-1\ :\ %s%d%1d%1d-%1d%1d.png, ",k1,l1,k2,l2,\ subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\ subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); - fprintf(fichtm,"\n
    ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); - fprintf(fichtm,"\n
    Correlation at age %d (%.3f),",(int) age, c12); + fprintf(fichtmcov,"\n
    ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); + fprintf(fichtmcov,"\n
    Correlation at age %d (%.3f),",(int) age, c12); fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); @@ -3008,7 +3062,7 @@ void varprob(char optionfilefiname[], do mu2,std,v21,sqrt(lc1),v22,sqrt(lc2)); }else{ first=0; - fprintf(fichtm," %d (%.3f),",(int) age, c12); + fprintf(fichtmcov," %d (%.3f),",(int) age, c12); fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); fprintf(ficgp,"\nreplot %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\ @@ -3031,7 +3085,8 @@ void varprob(char optionfilefiname[], do fclose(ficresprob); fclose(ficresprobcov); fclose(ficresprobcor); - /* fclose(ficgp);*/ + fflush(ficgp); + fflush(fichtmcov); } @@ -3137,8 +3192,8 @@ fprintf(fichtm," \n",\ - fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt,\ - model,fileres,fileres,\ + fileres,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model,\ + fileres,fileres,\ filelog,filelog,optionfilegnuplot,optionfilegnuplot,strstart); - /*fclose(fichtm);*/ fflush(fichtm); strcpy(pathr,path); strcat(pathr,optionfilefiname); chdir(optionfilefiname); /* Move to directory named optionfile */ - strcpy(lfileres,fileres); - strcat(lfileres,"/"); - strcat(lfileres,optionfilefiname); /* Calculates basic frequencies. Computes observed prevalence at single age and prints on file fileres'p'. */ @@ -4538,27 +4619,25 @@ Interval (in months) between two waves: fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); - jk=1; fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); for(i=1,jk=1; i <=nlstate; i++){ for(k=1; k <=(nlstate+ndeath); k++){ - if (k != i) - { - printf("%d%d ",i,k); - fprintf(ficlog,"%d%d ",i,k); - fprintf(ficres,"%1d%1d ",i,k); - for(j=1; j <=ncovmodel; j++){ - printf("%f ",p[jk]); - fprintf(ficlog,"%f ",p[jk]); - fprintf(ficres,"%f ",p[jk]); - jk++; - } - printf("\n"); - fprintf(ficlog,"\n"); - fprintf(ficres,"\n"); + if (k != i) { + printf("%d%d ",i,k); + fprintf(ficlog,"%d%d ",i,k); + fprintf(ficres,"%1d%1d ",i,k); + for(j=1; j <=ncovmodel; j++){ + printf("%f ",p[jk]); + fprintf(ficlog,"%f ",p[jk]); + fprintf(ficres,"%f ",p[jk]); + jk++; } + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficres,"\n"); + } } } if(mle!=0){ @@ -4589,32 +4668,91 @@ Interval (in months) between two waves: } fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); - if(mle==1) + if(mle>=1) printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); - for(i=1,k=1;i<=npar;i++){ - /* if (k>nlstate) k=1; - i1=(i-1)/(ncovmodel*nlstate)+1; - fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]); - printf("%s%d%d",alph[k],i1,tab[i]); - */ - fprintf(ficres,"%3d",i); - if(mle==1) - printf("%3d",i); - fprintf(ficlog,"%3d",i); - for(j=1; j<=i;j++){ - fprintf(ficres," %.5e",matcov[i][j]); - if(mle==1) - printf(" %.5e",matcov[i][j]); - fprintf(ficlog," %.5e",matcov[i][j]); - } - fprintf(ficres,"\n"); - if(mle==1) - printf("\n"); - fprintf(ficlog,"\n"); - k++; - } - +/* # 121 Var(a12)\n\ */ +/* # 122 Cov(b12,a12) Var(b12)\n\ */ +/* # 131 Cov(a13,a12) Cov(a13,b12, Var(a13)\n\ */ +/* # 132 Cov(b13,a12) Cov(b13,b12, Cov(b13,a13) Var(b13)\n\ */ +/* # 212 Cov(a21,a12) Cov(a21,b12, Cov(a21,a13) Cov(a21,b13) Var(a21)\n\ */ +/* # 212 Cov(b21,a12) Cov(b21,b12, Cov(b21,a13) Cov(b21,b13) Cov(b21,a21) Var(b21)\n\ */ +/* # 232 Cov(a23,a12) Cov(a23,b12, Cov(a23,a13) Cov(a23,b13) Cov(a23,a21) Cov(a23,b21) Var(a23)\n\ */ +/* # 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n" */ + + +/* Just to have a covariance matrix which will be more understandable + even is we still don't want to manage dictionary of variables +*/ + for(itimes=1;itimes<=2;itimes++){ + jj=0; + for(i=1; i <=nlstate; i++){ + for(j=1; j <=nlstate+ndeath; j++){ + if(j==i) continue; + for(k=1; k<=ncovmodel;k++){ + jj++; + ca[0]= k+'a'-1;ca[1]='\0'; + if(itimes==1){ + if(mle>=1) + printf("#%1d%1d%d",i,j,k); + fprintf(ficlog,"#%1d%1d%d",i,j,k); + fprintf(ficres,"#%1d%1d%d",i,j,k); + }else{ + if(mle>=1) + printf("%1d%1d%d",i,j,k); + fprintf(ficlog,"%1d%1d%d",i,j,k); + fprintf(ficres,"%1d%1d%d",i,j,k); + } + ll=0; + for(li=1;li <=nlstate; li++){ + for(lj=1;lj <=nlstate+ndeath; lj++){ + if(lj==li) continue; + for(lk=1;lk<=ncovmodel;lk++){ + ll++; + if(ll<=jj){ + cb[0]= lk +'a'-1;cb[1]='\0'; + if(ll=1) + printf(" Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); + fprintf(ficlog," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); + fprintf(ficres," Cov(%s%1d%1d,%s%1d%1d)",ca,i,j,cb, li,lj); + }else{ + if(mle>=1) + printf(" %.5e",matcov[jj][ll]); + fprintf(ficlog," %.5e",matcov[jj][ll]); + fprintf(ficres," %.5e",matcov[jj][ll]); + } + }else{ + if(itimes==1){ + if(mle>=1) + printf(" Var(%s%1d%1d)",ca,i,j); + fprintf(ficlog," Var(%s%1d%1d)",ca,i,j); + fprintf(ficres," Var(%s%1d%1d)",ca,i,j); + }else{ + if(mle>=1) + printf(" %.5e",matcov[jj][ll]); + fprintf(ficlog," %.5e",matcov[jj][ll]); + fprintf(ficres," %.5e",matcov[jj][ll]); + } + } + } + } /* end lk */ + } /* end lj */ + } /* end li */ + if(mle>=1) + printf("\n"); + fprintf(ficlog,"\n"); + fprintf(ficres,"\n"); + numlinepar++; + } /* end k*/ + } /*end j */ + } /* end i */ + } /* end itimes */ + + fflush(ficlog); + fflush(ficres); + while((c=getc(ficpar))=='#' && c!= EOF){ ungetc(c,ficpar); fgets(line, MAXLINE, ficpar); @@ -5024,60 +5162,58 @@ ageminpar, agemax, s[lastpass][imx], age fflush(ficgp); - if(erreur >0){ - printf("End of Imach with error or warning %d\n",erreur); - fprintf(ficlog,"End of Imach with error or warning %d\n",erreur); + if((nberr >0) || (nbwarn>0)){ + printf("End of Imach with %d errors and/or %d warnings\n",nberr,nbwarn); + fprintf(ficlog,"End of Imach with %d errors and/or warnings %d\n",nberr,nbwarn); }else{ - printf("End of Imach\n"); - fprintf(ficlog,"End of Imach\n"); + printf("End of Imach\n"); + fprintf(ficlog,"End of Imach\n"); } printf("See log file on %s\n",filelog); - fclose(ficlog); /* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */ (void) gettimeofday(&end_time,&tzp); tm = *localtime(&end_time.tv_sec); tmg = *gmtime(&end_time.tv_sec); strcpy(strtend,asctime(&tm)); - printf("Localtime at start %s\nLocaltime at end %s",strstart, strtend); - fprintf(ficlog,"Localtime at start %s\nLocal time at end %s",strstart, strtend); - /* printf("Total time used %d Sec\n", asc_time(end_time.tv_sec -start_time.tv_sec);*/ - - printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec); - fprintf(ficlog,"Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec); + printf("Local time at start %s\nLocaltime at end %s",strstart, strtend); + fprintf(ficlog,"Local time at start %s\nLocal time at end %s\n",strstart, strtend); + printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); + + printf("Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); + fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); + fprintf(ficlog,"Total time was %d Sec.\n", end_time.tv_sec -start_time.tv_sec); /* printf("Total time was %d uSec.\n", total_usecs);*/ /* if(fileappend(fichtm,optionfilehtm)){ */ fprintf(fichtm,"
    Local time at start %s
    Local time at end %s
    ",strstart, strtend); fclose(fichtm); + fclose(fichtmcov); fclose(ficgp); + fclose(ficlog); /*------ End -----------*/ - end: -#ifdef windows - /* chdir(pathcd);*/ -#endif chdir(path); - /*system("wgnuplot graph.plt");*/ - /*system("../gp37mgw/wgnuplot graph.plt");*/ - /*system("cd ../gp37mgw");*/ - /* system("..\\gp37mgw\\wgnuplot graph.plt");*/ strcpy(plotcmd,GNUPLOTPROGRAM); strcat(plotcmd," "); strcat(plotcmd,optionfilegnuplot); printf("Starting graphs with: %s",plotcmd);fflush(stdout); - system(plotcmd); + if((outcmd=system(plotcmd)) != 0){ + printf(" Problem with gnuplot\n"); + } printf(" Wait..."); - - /*#ifdef windows*/ while (z[0] != 'q') { /* chdir(path); */ - printf("\nType e to edit output files, g to graph again, c to start again, and q for exiting: "); + printf("\nType e to edit output files, g to graph again and q for exiting: "); scanf("%s",z); - if (z[0] == 'c') system("./imach"); - else if (z[0] == 'e') system(optionfilehtm); +/* if (z[0] == 'c') system("./imach"); */ + if (z[0] == 'e') system(optionfilehtm); else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); } - /*#endif */ + end: + while (z[0] != 'q') { + printf("\nType q for exiting: "); + scanf("%s",z); + } }