--- imach/src/imach.c 2014/06/16 10:33:11 1.147 +++ imach/src/imach.c 2014/12/22 11:40:47 1.166 @@ -1,6 +1,87 @@ -/* $Id: imach.c,v 1.147 2014/06/16 10:33:11 brouard Exp $ +/* $Id: imach.c,v 1.166 2014/12/22 11:40:47 brouard Exp $ $State: Exp $ $Log: imach.c,v $ + Revision 1.166 2014/12/22 11:40:47 brouard + *** empty log message *** + + Revision 1.165 2014/12/16 11:20:36 brouard + Summary: After compiling on Visual C + + * imach.c (Module): Merging 1.61 to 1.162 + + Revision 1.164 2014/12/16 10:52:11 brouard + Summary: Merging with Visual C after suppressing some warnings for unused variables. Also fixing Saito's bug 0.98Xn + + * imach.c (Module): Merging 1.61 to 1.162 + + Revision 1.163 2014/12/16 10:30:11 brouard + * imach.c (Module): Merging 1.61 to 1.162 + + Revision 1.162 2014/09/25 11:43:39 brouard + Summary: temporary backup 0.99! + + Revision 1.1 2014/09/16 11:06:58 brouard + Summary: With some code (wrong) for nlopt + + Author: + + Revision 1.161 2014/09/15 20:41:41 brouard + Summary: Problem with macro SQR on Intel compiler + + Revision 1.160 2014/09/02 09:24:05 brouard + *** empty log message *** + + Revision 1.159 2014/09/01 10:34:10 brouard + Summary: WIN32 + Author: Brouard + + Revision 1.158 2014/08/27 17:11:51 brouard + *** empty log message *** + + Revision 1.157 2014/08/27 16:26:55 brouard + Summary: Preparing windows Visual studio version + Author: Brouard + + In order to compile on Visual studio, time.h is now correct and time_t + and tm struct should be used. difftime should be used but sometimes I + just make the differences in raw time format (time(&now). + Trying to suppress #ifdef LINUX + Add xdg-open for __linux in order to open default browser. + + Revision 1.156 2014/08/25 20:10:10 brouard + *** empty log message *** + + Revision 1.155 2014/08/25 18:32:34 brouard + Summary: New compile, minor changes + Author: Brouard + + Revision 1.154 2014/06/20 17:32:08 brouard + Summary: Outputs now all graphs of convergence to period prevalence + + Revision 1.153 2014/06/20 16:45:46 brouard + Summary: If 3 live state, convergence to period prevalence on same graph + Author: Brouard + + Revision 1.152 2014/06/18 17:54:09 brouard + Summary: open browser, use gnuplot on same dir than imach if not found in the path + + Revision 1.151 2014/06/18 16:43:30 brouard + *** empty log message *** + + Revision 1.150 2014/06/18 16:42:35 brouard + Summary: If gnuplot is not in the path try on same directory than imach binary (OSX) + Author: brouard + + Revision 1.149 2014/06/18 15:51:14 brouard + Summary: Some fixes in parameter files errors + Author: Nicolas Brouard + + Revision 1.148 2014/06/17 17:38:48 brouard + Summary: Nothing new + Author: Brouard + + Just a new packaging for OS/X version 0.98nS + Revision 1.147 2014/06/16 10:33:11 brouard *** empty log message *** @@ -424,33 +505,46 @@ end */ +#define POWELL /* Instead of NLOPT */ - - #include #include #include #include + +#ifdef _WIN32 +#include +#else #include +#endif #include #include #include #include -extern int errno; +/* extern int errno; */ + +/* #ifdef LINUX */ +/* #include */ +/* #include "timeval.h" */ +/* #else */ +/* #include */ +/* #endif */ -#ifdef LINUX #include -#include "timeval.h" -#else -#include -#endif #ifdef GSL #include #include #endif +#ifdef NLOPT +#include +typedef struct { + double (* function)(double [] ); +} myfunc_data ; +#endif + /* #include */ /* #define _(String) gettext (String) */ @@ -475,22 +569,22 @@ extern int errno; #define YEARM 12. /**< Number of months per year */ #define AGESUP 130 #define AGEBASE 40 -#define AGEGOMP 10. /**< Minimal age for Gompertz adjustment */ -#ifdef UNIX -#define DIRSEPARATOR '/' -#define CHARSEPARATOR "/" -#define ODIRSEPARATOR '\\' -#else +#define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */ +#ifdef _WIN32 #define DIRSEPARATOR '\\' #define CHARSEPARATOR "\\" #define ODIRSEPARATOR '/' +#else +#define DIRSEPARATOR '/' +#define CHARSEPARATOR "/" +#define ODIRSEPARATOR '\\' #endif -/* $Id: imach.c,v 1.147 2014/06/16 10:33:11 brouard Exp $ */ +/* $Id: imach.c,v 1.166 2014/12/22 11:40:47 brouard Exp $ */ /* $State: Exp $ */ -char version[]="Imach version 0.98nR2, January 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)"; -char fullversion[]="$Revision: 1.147 $ $Date: 2014/06/16 10:33:11 $"; +char version[]="Imach version 0.99, September 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)"; +char fullversion[]="$Revision: 1.166 $ $Date: 2014/12/22 11:40:47 $"; char strstart[80]; char optionfilext[10], optionfilefiname[FILENAMELENGTH]; int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ @@ -520,6 +614,7 @@ int **mw; /* mw[mi][i] is number of the int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between * wave mi and wave mi+1 is not an exact multiple of stepm. */ +int countcallfunc=0; /* Count the number of calls to func */ double jmean=1; /* Mean space between 2 waves */ double **matprod2(); /* test */ double **oldm, **newm, **savm; /* Working pointers to matrices */ @@ -563,12 +658,17 @@ char popfile[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(); +/* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */ +/* struct timezone tzp; */ +/* extern int gettimeofday(); */ +struct tm tml, *gmtime(), *localtime(); + +extern time_t time(); + +struct tm start_time, end_time, curr_time, last_time, forecast_time; +time_t rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */ +struct tm tm; + char strcurr[80], strfor[80]; char *endptr; @@ -598,7 +698,12 @@ static double maxarg1,maxarg2; #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a)) #define rint(a) floor(a+0.5) - +/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */ +/* #define mytinydouble 1.0e-16 */ +/* #define DEQUAL(a,b) (fabs((a)-(b)) *fc) { + while (*fb > *fc) { /* Declining fa, fb, fc */ r=(*bx-*ax)*(*fb-*fc); q=(*bx-*cx)*(*fb-*fa); u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ - (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); - ulim=(*bx)+GLIMIT*(*cx-*bx); - if ((*bx-u)*(u-*cx) > 0.0) { + (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscisse of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */ + ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */ + if ((*bx-u)*(u-*cx) > 0.0) { /* if u between b and c */ fu=(*func)(u); - } else if ((*cx-u)*(u-ulim) > 0.0) { +#ifdef DEBUG + /* f(x)=A(x-u)**2+f(u) */ + double A, fparabu; + A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u); + fparabu= *fa - A*(*ax-u)*(*ax-u); + printf("mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu); + fprintf(ficlog, "mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu); +#endif + } else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */ fu=(*func)(u); if (fu < *fc) { SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) SHFT(*fb,*fc,fu,(*func)(u)) } - } else if ((u-ulim)*(ulim-*cx) >= 0.0) { + } else if ((u-ulim)*(ulim-*cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */ u=ulim; fu=(*func)(u); } else { @@ -1141,7 +1285,11 @@ void mnbrak(double *ax, double *bx, doub } /*************** linmin ************************/ - +/* Given an n -dimensional point p[1..n] and an n -dimensional direction xi[1..n] , moves and +resets p to where the function func(p) takes on a minimum along the direction xi from p , +and replaces xi by the actual vector displacement that p was moved. Also returns as fret +the value of func at the returned location p . This is actually all accomplished by calling the +routines mnbrak and brent .*/ int ncom; double *pcom,*xicom; double (*nrfunc)(double []); @@ -1167,8 +1315,8 @@ void linmin(double p[], double xi[], int } ax=0.0; xx=1.0; - mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); - *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); + mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Find a bracket a,x,b in direction n=xi ie xicom */ + *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */ #ifdef DEBUG printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); @@ -1181,20 +1329,16 @@ 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,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left); - return ascdiff; -} /*************** powell ************************/ +/* +Minimization of a function func of n variables. Input consists of an initial starting point +p[1..n] ; an initial matrix xi[1..n][1..n] , whose columns contain the initial set of di- +rections (usually the n unit vectors); and ftol , the fractional tolerance in the function value +such that failure to decrease by more than this amount on one iteration signals doneness. On +output, p is set to the best point found, xi is the then-current direction set, fret is the returned +function value at p , and iter is the number of iterations taken. The routine linmin is used. + */ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, double (*func)(double [])) { @@ -1212,15 +1356,18 @@ void powell(double p[], double **xi, int xits=vector(1,n); *fret=(*func)(p); for (j=1;j<=n;j++) pt[j]=p[j]; + rcurr_time = time(NULL); for (*iter=1;;++(*iter)) { fp=(*fret); ibig=0; del=0.0; - 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); fflush(ficlog); -/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */ + rlast_time=rcurr_time; + /* (void) gettimeofday(&curr_time,&tzp); */ + rcurr_time = time(NULL); + curr_time = *localtime(&rcurr_time); + printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); + fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); +/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ for (i=1;i<=n;i++) { printf(" %d %.12f",i, p[i]); fprintf(ficlog," %d %.12lf",i, p[i]); @@ -1230,33 +1377,31 @@ void powell(double p[], double **xi, int fprintf(ficlog,"\n"); fprintf(ficrespow,"\n");fflush(ficrespow); if(*iter <=3){ - tm = *localtime(&curr_time.tv_sec); - strcpy(strcurr,asctime(&tm)); -/* asctime_r(&tm,strcurr); */ - forecast_time=curr_time; + tml = *localtime(&rcurr_time); + strcpy(strcurr,asctime(&tml)); + rforecast_time=rcurr_time; itmp = strlen(strcurr); 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); + printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); + fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time); 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)); + rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); + forecast_time = *localtime(&rforecast_time); + strcpy(strfor,asctime(&forecast_time)); 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 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); + 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(rforecast_time-rcurr_time,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(rforecast_time-rcurr_time,tmpout),strfor,strcurr); } } for (i=1;i<=n;i++) { for (j=1;j<=n;j++) xit[j]=xi[j][i]; fptt=(*fret); #ifdef DEBUG - printf("fret=%lf \n",*fret); - fprintf(ficlog,"fret=%lf \n",*fret); + printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); + fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); #endif printf("%d",i);fflush(stdout); fprintf(ficlog,"%d",i);fflush(ficlog); @@ -1274,13 +1419,13 @@ void powell(double p[], double **xi, int fprintf(ficlog," x(%d)=%.12e",j,xit[j]); } for(j=1;j<=n;j++) { - printf(" p=%.12e",p[j]); - fprintf(ficlog," p=%.12e",p[j]); + printf(" p(%d)=%.12e",j,p[j]); + fprintf(ficlog," p(%d)=%.12e",j,p[j]); } printf("\n"); fprintf(ficlog,"\n"); #endif - } + } /* end i */ if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { #ifdef DEBUG int k[2],l; @@ -1313,20 +1458,45 @@ void powell(double p[], double **xi, int return; } if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); - for (j=1;j<=n;j++) { + for (j=1;j<=n;j++) { /* Computes an extrapolated point */ ptt[j]=2.0*p[j]-pt[j]; xit[j]=p[j]-pt[j]; pt[j]=p[j]; } fptt=(*func)(ptt); - if (fptt < fp) { - t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); - if (t < 0.0) { - linmin(p,xit,n,fret,func); + if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */ + /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */ + /* From x1 (P0) distance of x2 is at h and x3 is 2h */ + /* Let f"(x2) be the 2nd derivative equal everywhere. */ + /* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */ + /* will reach at f3 = fm + h^2/2 f"m ; f" = (f1 -2f2 +f3 ) / h**2 */ + /* f1-f3 = delta(2h) = 2 h**2 f'' = 2(f1- 2f2 +f3) */ + /* Thus we compare delta(2h) with observed f1-f3 */ + /* or best gain on one ancient line 'del' with total */ + /* gain f1-f2 = f1 - f2 - 'del' with del */ + /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */ + + t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); + t= t- del*SQR(fp-fptt); + printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t); + fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t); +#ifdef DEBUG + printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), + (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt)); + fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), + (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt)); + printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t); + fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t); +#endif + if (t < 0.0) { /* Then we use it for last direction */ + linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/ for (j=1;j<=n;j++) { - xi[j][ibig]=xi[j][n]; - xi[j][n]=xit[j]; + xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */ + xi[j][n]=xit[j]; /* and nth direction by the extrapolated */ } + printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); + fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction :\n",n,ibig); + #ifdef DEBUG printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); @@ -1337,8 +1507,8 @@ void powell(double p[], double **xi, int printf("\n"); fprintf(ficlog,"\n"); #endif - } - } + } /* end of t negative */ + } /* end if (fptt < fp) */ } } @@ -1426,7 +1596,7 @@ double **pmij(double **ps, double *cov, */ double s1, lnpijopii; /*double t34;*/ - int i,j,j1, nc, ii, jj; + int i,j, nc, ii, jj; for(i=1; i<= nlstate; i++){ for(j=1; jfunction)(xt); /* p xt[1]@8 is fine */ + /* fret=(*func)(xt); /\* p xt[1]@8 is fine *\/ */ + printf("Function = %.12lf ",fret); + for (j=1;j<=n;j++) printf(" %d %.8lf", j, xt[j]); + printf("\n"); + free_vector(xt,1,n); + return fret; +} +#endif /*************** log-likelihood *************/ double func( double *x) @@ -1586,6 +1775,9 @@ double func( double *x) /*for(i=1;ifunction = func; + printf(" Func %.12lf \n",myfunc(npar,p1,NULL,d)); + nlopt_set_min_objective(opt, myfunc, d); + nlopt_set_xtol_rel(opt, ftol); + if ((creturn=nlopt_optimize(opt, p1, &minf)) < 0) { + printf("nlopt failed! %d\n",creturn); + } + else { + printf("found minimum after %d evaluations (NLOPT=%d)\n", countcallfunc ,NLOPT); + printf("found minimum at f(%g,%g) = %0.10g\n", p[0], p[1], minf); + iter=1; /* not equal */ + } + nlopt_destroy(opt); +#endif free_matrix(xi,1,npar,1,npar); fclose(ficrespow); - printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p)); - fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); - fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); + printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); + fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); } @@ -2004,7 +2235,7 @@ void hesscov(double **matcov, double p[] { double **a,**y,*x,pd; double **hess; - int i, j,jk; + int i, j; int *indx; double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar); @@ -2138,7 +2369,7 @@ double hessii(double x[], double delta, k=kmax; } else if((k1 >khi/nkhif) || (k2 >khi/nkhif)){ /* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. */ - k=kmax; l=lmax*10.; + k=kmax; l=lmax*10; } else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ delts=delt; @@ -2153,7 +2384,7 @@ double hessii(double x[], double delta, double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) { int i; - int l=1, l1, lmax=20; + int l=1, lmax=20; double k1,k2,k3,k4,res,fx; double p2[MAXPARM+1]; int k; @@ -2268,7 +2499,7 @@ void pstamp(FILE *fichier) void freqsummary(char fileres[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[]) { /* Some frequencies */ - int i, m, jk, k1,i1, j1, bool, z1,j; + int i, m, jk, j1, bool, z1,j; int first; double ***freq; /* Frequencies */ double *pp, **prop; @@ -2449,10 +2680,10 @@ void prevalence(double ***probs, double We still use firstpass and lastpass as another selection. */ - int i, m, jk, k1, i1, j1, bool, z1,j; - double ***freq; /* Frequencies */ - double *pp, **prop; - double pos,posprop; + int i, m, jk, j1, bool, z1,j; + + double **prop; + double posprop; double y2; /* in fractional years */ int iagemin, iagemax; int first; /** to stop verbosity which is redirected to log file */ @@ -2543,7 +2774,7 @@ void concatwav(int wav[], int **dh, int int j, k=0,jk, ju, jl; double sum=0.; first=0; - jmin=1e+5; + jmin=100000; jmax=-1; jmean=0.; for(i=1; i<=imx; i++){ @@ -2783,7 +3014,7 @@ void evsij(double ***eij, double x[], in { /* Health expectancies, no variances */ - int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2; + int i, j, nhstepm, hstepm, h, nstepm; int nhstepma, nstepma; /* Decreasing with age */ double age, agelim, hf; double ***p3mat; @@ -3105,7 +3336,7 @@ void varevsij(char optionfilefiname[], d double **dnewm,**doldm; double **dnewmp,**doldmp; int i, j, nhstepm, hstepm, h, nstepm ; - int k, cptcode; + int k; double *xp; double **gp, **gm; /* for var eij */ double ***gradg, ***trgradg; /*for var eij */ @@ -3405,10 +3636,9 @@ void varprevlim(char fileres[], double * { /* Variance of prevalence limit */ /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ - double **newm; + double **dnewm,**doldm; int i, j, nhstepm, hstepm; - int k, cptcode; double *xp; double *gp, *gm; double **gradg, **trgradg; @@ -3487,9 +3717,9 @@ void varprevlim(char fileres[], double * /************ Variance of one-step probabilities ******************/ void varprob(char optionfilefiname[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int *Tvar, int **nbcode, int *ncodemax, char strstart[]) { - int i, j=0, i1, k1, l1, t, tj; + int i, j=0, k1, l1, tj; int k2, l2, j1, z1; - int k=0,l, cptcode; + int k=0, l; int first=1, first1, first2; double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; double **dnewm,**doldm; @@ -3497,7 +3727,7 @@ void varprob(char optionfilefiname[], do double *gp, *gm; double **gradg, **trgradg; double **mu; - double age,agelim, cov[NCOVMAX+1]; + double age, cov[NCOVMAX+1]; double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ int theta; char fileresprob[FILENAMELENGTH]; @@ -3868,13 +4098,13 @@ fprintf(fichtm," \n
  • Graphs before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: %s%d_2.png
    \ ",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); /* Period (stable) prevalence in each health state */ - for(cpt=1; cpt- Period (stable) prevalence in each health state : %s%d_%d.png
    \ -",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); + for(cpt=1; cpt<=nlstate;cpt++){ + fprintf(fichtm,"
    - Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. %s%d_%d.png
    \ +", cpt, cpt, nlstate, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); } for(cpt=1; cpt<=nlstate;cpt++) { - fprintf(fichtm,"\n
    - Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : %s%d%d.png
    \ -",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); + fprintf(fichtm,"\n
    - Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : %s%d%d.png
    \ +",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); } } /* end i1 */ }/* End k1 */ @@ -3955,7 +4185,7 @@ true period expectancies (those weighted void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ char dirfileres[132],optfileres[132]; - int m0,cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; + int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0; int ng=0; /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ /* printf("Problem with file %s",optionfilegnuplot); */ @@ -3970,6 +4200,7 @@ void printinggnuplot(char fileres[], cha strcpy(dirfileres,optionfilefiname); strcpy(optfileres,"vpl"); /* 1eme*/ + fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n"); for (cpt=1; cpt<= nlstate ; cpt ++) { for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1); @@ -3997,7 +4228,7 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u } } /*2 eme*/ - + fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); for (k1=1; k1<= m ; k1 ++) { fprintf(ficgp,"\nset out \"%s%d.png\" \n",subdirf2(optionfilefiname,"e"),k1); fprintf(ficgp,"set ylabel \"Years\" \nset ter png small size 320, 240\nplot [%.f:%.f] ",ageminpar,fage); @@ -4054,28 +4285,29 @@ plot [%.f:%.f] \"%s\" every :::%d::%d u } /* CV preval stable (period) */ - for (k1=1; k1<= m ; k1 ++) { - for (cpt=1; cpt<=nlstate ; cpt ++) { + for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */ + for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ k=3; + fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt); fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1); fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ set ter png small size 320, 240\n\ unset log y\n\ -plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,subdirf2(fileres,"pij"),k1,k+cpt+1,k+1); - - for (i=1; i< nlstate ; i ++) - fprintf(ficgp,"+$%d",k+i+1); - fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1); - - l=3+(nlstate+ndeath)*cpt; - fprintf(ficgp,",\"%s\" u ($1==%d ? ($3):1/0):($%d/($%d",subdirf2(fileres,"pij"),k1,l+cpt+1,l+1); - for (i=1; i< nlstate ; i ++) { - l=3+(nlstate+ndeath)*cpt; - fprintf(ficgp,"+$%d",l+i+1); - } - fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1); - } - } +plot [%.f:%.f] ", ageminpar, agemaxpar); + for (i=1; i<= nlstate ; i ++){ + if(i==1) + fprintf(ficgp,"\"%s\"",subdirf2(fileres,"pij")); + else + fprintf(ficgp,", '' "); + l=(nlstate+ndeath)*(i-1)+1; + fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); + for (j=1; j<= (nlstate-1) ; j ++) + fprintf(ficgp,"+$%d",k+l+j); + fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt); + } /* nlstate */ + fprintf(ficgp,"\n"); + } /* end cpt state*/ + } /* end covariate */ /* proba elementaires */ for(i=1,jk=1; i <=nlstate; i++){ @@ -4139,7 +4371,7 @@ plot [%.f:%.f] \"%s\" u ($1==%d ? ($3):1 } /* end k2 */ } /* end jk */ } /* end ng */ - avoid: + /* avoid: */ fflush(ficgp); } /* end gnuplot */ @@ -4193,8 +4425,7 @@ prevforecast(char fileres[], double anpr dateprev1 dateprev2 range of dates during which prevalence is computed anproj2 year of en of projection (same day and month as proj1). */ - int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h, i1; - int *popage; + int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1; double agec; /* generic age */ double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; double *popeffectif,*popcount; @@ -4488,8 +4719,8 @@ void prwizard(int ncovmodel, int nlstate /* Wizard to print covariance matrix template */ - char ca[32], cb[32], cc[32]; - int i,j, k, l, li, lj, lk, ll, jj, npar, itimes; + char ca[32], cb[32]; + int i,j, k, li, lj, lk, ll, jj, npar, itimes; int numlinepar; printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); @@ -4710,7 +4941,7 @@ fprintf(fichtm,"

    • Life table=1;j--){ @@ -4907,7 +5136,7 @@ int readdata(char datafile[], int firsto fclose(fic); return (0); - endread: + /* endread: */ printf("Exiting readdata: "); fclose(fic); return (1); @@ -4938,7 +5167,7 @@ int decodemodel ( char model[], int last */ { int i, j, k, ks; - int i1, j1, k1, k2; + int j1, k1, k2; char modelsav[80]; char stra[80], strb[80], strc[80], strd[80],stre[80]; @@ -5087,7 +5316,7 @@ int decodemodel ( char model[], int last return (0); /* with covar[new additional covariate if product] and Tage if age */ - endread: + /*endread:*/ printf("Exiting decodemodel: "); return (1); } @@ -5146,7 +5375,7 @@ calandcheckages(int imx, int maxwav, dou } else if(agev[m][i] >*agemax){ *agemax=agev[m][i]; - printf(" Max anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.2f\n",m,i,anint[m][i], i,annais[i], *agemax); + /* printf(" Max anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.2f\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;*/ @@ -5184,7 +5413,7 @@ calandcheckages(int imx, int maxwav, dou fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax); return (0); - endread: + /* endread:*/ printf("Exiting calandcheckages: "); return (1); } @@ -5204,38 +5433,35 @@ int main(int argc, char *argv[]) double ssval; #endif 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 i,j, k, n=MAXN,iter=0,m,size=100, cptcod; + + int jj, ll, li, lj, lk; int numlinepar=0; /* Current linenumber of parameter file */ int itimes; int NDIM=2; int vpopbased=0; - char ca[32], cb[32], cc[32]; + char ca[32], cb[32]; /* FILE *fichtm; *//* Html File */ /* FILE *ficgp;*/ /*Gnuplot File */ struct stat info; - double agedeb, agefin,hf; + double agedeb; double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; double fret; - double **xi,tmp,delta; - double dum; /* Dummy variable */ double ***p3mat; double ***mobaverage; - int *indx; - char line[MAXLINE], linepar[MAXLINE]; + + char line[MAXLINE]; char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; char pathr[MAXLINE], pathimach[MAXLINE]; - char **bp, *tok, *val; /* pathtot */ + char *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, jk,aa,bb, stepsize, ij; - int jnais,jdc,jint4,jint1,jint2,jint3,*tab; + int c, h , cpt; + int jl; + int i1, j1, jk, stepsize; + int *tab; int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ int mobilav=0,popforecast=0; int hstepm, nhstepm; @@ -5244,7 +5470,7 @@ int main(int argc, char *argv[]) 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 bage=0, fage=110, age, agelim, agebase; double ftolpl=FTOL; double **prlim; double ***param; /* Matrix of parameters */ @@ -5255,19 +5481,18 @@ int main(int argc, char *argv[]) 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]="1234"; int *dcwave; - char z[1]="c", occ; + char z[1]="c"; /*char *strt;*/ char strtend[80]; - long total_usecs; - + /* setlocale (LC_ALL, ""); */ /* bindtextdomain (PACKAGE, LOCALEDIR); */ /* textdomain (PACKAGE); */ @@ -5275,19 +5500,21 @@ int main(int argc, char *argv[]) /* setlocale (LC_MESSAGES, ""); */ /* gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ - (void) gettimeofday(&start_time,&tzp); + rstart_time = time(NULL); + /* (void) gettimeofday(&start_time,&tzp);*/ + start_time = *localtime(&rstart_time); curr_time=start_time; - tm = *localtime(&start_time.tv_sec); - tmg = *gmtime(&start_time.tv_sec); - strcpy(strstart,asctime(&tm)); + /*tml = *localtime(&start_time.tm_sec);*/ + /* strcpy(strstart,asctime(&tml)); */ + strcpy(strstart,asctime(&start_time)); /* printf("Localtime (at start)=%s",strstart); */ -/* tp.tv_sec = tp.tv_sec +86400; */ -/* tm = *localtime(&start_time.tv_sec); */ +/* tp.tm_sec = tp.tm_sec +86400; */ +/* tm = *localtime(&start_time.tm_sec); */ /* tmg.tm_year=tmg.tm_year +dsign*dyear; */ /* tmg.tm_mon=tmg.tm_mon +dsign*dmonth; */ /* tmg.tm_hour=tmg.tm_hour + 1; */ -/* tp.tv_sec = mktime(&tmg); */ +/* tp.tm_sec = mktime(&tmg); */ /* strt=asctime(&tmg); */ /* printf("Time(after) =%s",strstart); */ /* (void) time (&time_value); @@ -5308,6 +5535,9 @@ int main(int argc, char *argv[]) i=strlen(pathr); if(pathr[i-1]=='\n') pathr[i-1]='\0'; + i=strlen(pathr); + if(pathr[i-1]==' ') /* This may happen when dragging on oS/X! */ + pathr[i-1]='\0'; for (tok = pathr; tok != NULL; ){ printf("Pathr |%s|\n",pathr); while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); @@ -5363,13 +5593,13 @@ int main(int argc, char *argv[]) path=%s \n\ optionfile=%s\n\ optionfilext=%s\n\ - optionfilefiname=%s\n",pathimach,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); fflush(ficlog); /* (void) gettimeofday(&curr_time,&tzp); */ -/* printf("Elapsed time %d\n", asc_diff_time(curr_time.tv_sec-start_time.tv_sec,tmpout)); */ +/* printf("Elapsed time %d\n", asc_diff_time(curr_time.tm_sec-start_time.tm_sec,tmpout)); */ /* */ strcpy(fileres,"r"); @@ -5379,10 +5609,11 @@ int main(int argc, char *argv[]) /*---------arguments file --------*/ if((ficpar=fopen(optionfile,"r"))==NULL) { - printf("Problem with optionfile %s\n",optionfile); - fprintf(ficlog,"Problem with optionfile %s\n",optionfile); + printf("Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno)); + fprintf(ficlog,"Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno)); fflush(ficlog); - goto end; + /* goto end; */ + exit(70); } @@ -5532,7 +5763,7 @@ run imach with mle=-1 to get a correct t for(i=1; i <=nlstate; i++){ for(j=1; j <=nlstate+ndeath-1; j++){ fscanf(ficpar,"%1d%1d",&i1,&j1); - if ((i1-i)*(j1-j)!=0){ + if ( (i1-i) * (j1-j) != 0){ printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); exit(1); } @@ -5901,7 +6132,7 @@ Interval (in months) between two waves: #ifdef GSL printf("GSL optimization\n"); fprintf(ficlog,"Powell\n"); -#elsedef +#else printf("Powell\n"); fprintf(ficlog,"Powell\n"); #endif strcpy(filerespow,"pow-mort"); @@ -5912,7 +6143,7 @@ Interval (in months) between two waves: } #ifdef GSL fprintf(ficrespow,"# GSL optimization\n# iter -2*LL"); -#elsedef +#else fprintf(ficrespow,"# Powell\n# iter -2*LL"); #endif /* for (i=1;i<=nlstate;i++) @@ -6474,7 +6705,8 @@ Interval (in months) between two waves: for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ oldm=oldms;savm=savms; /* Segmentation fault */ - varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); + cptcod= 0; /* To be deleted */ + varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); if(vpopbased==1) fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav); @@ -6566,7 +6798,7 @@ Interval (in months) between two waves: if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); } /* mle==-3 arrives here for freeing */ - endfree: + /* endfree:*/ free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */ free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); @@ -6600,17 +6832,18 @@ Interval (in months) between two waves: } printf("See log file on %s\n",filelog); /* 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)); + /*(void) gettimeofday(&end_time,&tzp);*/ + rend_time = time(NULL); + end_time = *localtime(&rend_time); + /* tml = *localtime(&end_time.tm_sec); */ + strcpy(strtend,asctime(&end_time)); printf("Local time at start %s\nLocal time 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 used %s\n", asc_diff_time(rend_time -rstart_time,tmpout)); - printf("Total time was %ld 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 %ld Sec.\n", end_time.tv_sec -start_time.tv_sec); + printf("Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time)); + fprintf(ficlog,"Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout)); + fprintf(ficlog,"Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time)); /* 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
      \n",strstart, strtend); @@ -6629,19 +6862,19 @@ Interval (in months) between two waves: printf("Current directory %s!\n",pathcd); /*strcat(plotcmd,CHARSEPARATOR);*/ sprintf(plotcmd,"gnuplot"); -#ifndef UNIX +#ifdef _WIN32 sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); #endif if(!stat(plotcmd,&info)){ - printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout); + printf("Error or gnuplot program not found: '%s'\n",plotcmd);fflush(stdout); if(!stat(getenv("GNUPLOTBIN"),&info)){ - printf("Error gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout); + printf("Error or gnuplot program not found: '%s' Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout); }else strcpy(pplotcmd,plotcmd); -#ifdef UNIX +#ifdef __unix strcpy(plotcmd,GNUPLOTPROGRAM); if(!stat(plotcmd,&info)){ - printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout); + printf("Error gnuplot program not found: '%s'\n",plotcmd);fflush(stdout); }else strcpy(pplotcmd,plotcmd); #endif @@ -6649,20 +6882,31 @@ Interval (in months) between two waves: strcpy(pplotcmd,plotcmd); sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); - printf("Starting graphs with: %s\n",plotcmd);fflush(stdout); + printf("Starting graphs with: '%s'\n",plotcmd);fflush(stdout); if((outcmd=system(plotcmd)) != 0){ - printf("\n Problem with gnuplot\n"); + printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd); + printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); + sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); + if((outcmd=system(plotcmd)) != 0) + printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); } - printf(" Wait..."); + printf(" Successful, please wait..."); while (z[0] != 'q') { /* chdir(path); */ - printf("\nType e to edit output files, g to graph again and q for exiting: "); + printf("\nType e to edit results with your browser, g to graph again and q for exit: "); scanf("%s",z); /* if (z[0] == 'c') system("./imach"); */ if (z[0] == 'e') { - printf("Starting browser with: %s",optionfilehtm);fflush(stdout); - system(optionfilehtm); +#ifdef __APPLE__ + sprintf(pplotcmd, "open %s", optionfilehtm); +#elif __linux + sprintf(pplotcmd, "xdg-open %s", optionfilehtm); +#else + sprintf(pplotcmd, "%s", optionfilehtm); +#endif + printf("Starting browser with: %s",pplotcmd);fflush(stdout); + system(pplotcmd); } else if (z[0] == 'g') system(plotcmd); else if (z[0] == 'q') exit(0); @@ -6673,6 +6917,3 @@ Interval (in months) between two waves: scanf("%s",z); } } - - -