|
version 1.147, 2014/06/16 10:33:11
|
version 1.163, 2014/12/16 10:30:11
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
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 |
Revision 1.147 2014/06/16 10:33:11 brouard |
| *** empty log message *** |
*** empty log message *** |
| |
|
|
Line 431
|
Line 499
|
| #include <stdio.h> |
#include <stdio.h> |
| #include <stdlib.h> |
#include <stdlib.h> |
| #include <string.h> |
#include <string.h> |
| |
|
| |
#ifdef _WIN32 |
| |
#include <io.h> |
| |
#else |
| #include <unistd.h> |
#include <unistd.h> |
| |
#endif |
| |
|
| #include <limits.h> |
#include <limits.h> |
| #include <sys/types.h> |
#include <sys/types.h> |
| #include <sys/stat.h> |
#include <sys/stat.h> |
| #include <errno.h> |
#include <errno.h> |
| extern int errno; |
/* extern int errno; */ |
| |
|
| |
/* #ifdef LINUX */ |
| |
/* #include <time.h> */ |
| |
/* #include "timeval.h" */ |
| |
/* #else */ |
| |
/* #include <sys/time.h> */ |
| |
/* #endif */ |
| |
|
| #ifdef LINUX |
|
| #include <time.h> |
#include <time.h> |
| #include "timeval.h" |
|
| #else |
|
| #include <sys/time.h> |
|
| #endif |
|
| |
|
| #ifdef GSL |
#ifdef GSL |
| #include <gsl/gsl_errno.h> |
#include <gsl/gsl_errno.h> |
| #include <gsl/gsl_multimin.h> |
#include <gsl/gsl_multimin.h> |
| #endif |
#endif |
| |
|
| |
#ifdef NLOPT |
| |
#include <nlopt.h> |
| |
typedef struct { |
| |
double (* function)(double [] ); |
| |
} myfunc_data ; |
| |
#endif |
| |
|
| /* #include <libintl.h> */ |
/* #include <libintl.h> */ |
| /* #define _(String) gettext (String) */ |
/* #define _(String) gettext (String) */ |
| |
|
|
Line 476 extern int errno;
|
Line 558 extern int errno;
|
| #define AGESUP 130 |
#define AGESUP 130 |
| #define AGEBASE 40 |
#define AGEBASE 40 |
| #define AGEGOMP 10. /**< Minimal age for Gompertz adjustment */ |
#define AGEGOMP 10. /**< Minimal age for Gompertz adjustment */ |
| #ifdef UNIX |
#ifdef _WIN32 |
| #define DIRSEPARATOR '/' |
|
| #define CHARSEPARATOR "/" |
|
| #define ODIRSEPARATOR '\\' |
|
| #else |
|
| #define DIRSEPARATOR '\\' |
#define DIRSEPARATOR '\\' |
| #define CHARSEPARATOR "\\" |
#define CHARSEPARATOR "\\" |
| #define ODIRSEPARATOR '/' |
#define ODIRSEPARATOR '/' |
| |
#else |
| |
#define DIRSEPARATOR '/' |
| |
#define CHARSEPARATOR "/" |
| |
#define ODIRSEPARATOR '\\' |
| #endif |
#endif |
| |
|
| /* $Id$ */ |
/* $Id$ */ |
| /* $State$ */ |
/* $State$ */ |
| |
|
| 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 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$ $Date$"; |
char fullversion[]="$Revision$ $Date$"; |
| char strstart[80]; |
char strstart[80]; |
| char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
char optionfilext[10], optionfilefiname[FILENAMELENGTH]; |
|
Line 520 int **mw; /* mw[mi][i] is number of the
|
Line 602 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 **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 |
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. */ |
* 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 jmean=1; /* Mean space between 2 waves */ |
| double **matprod2(); /* test */ |
double **matprod2(); /* test */ |
| double **oldm, **newm, **savm; /* Working pointers to matrices */ |
double **oldm, **newm, **savm; /* Working pointers to matrices */ |
|
Line 563 char popfile[FILENAMELENGTH];
|
Line 646 char popfile[FILENAMELENGTH];
|
| |
|
| char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilehtmcov[FILENAMELENGTH] ; |
char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilehtmcov[FILENAMELENGTH] ; |
| |
|
| struct timeval start_time, end_time, curr_time, last_time, forecast_time; |
/* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */ |
| struct timezone tzp; |
/* struct timezone tzp; */ |
| extern int gettimeofday(); |
/* extern int gettimeofday(); */ |
| struct tm tmg, tm, tmf, *gmtime(), *localtime(); |
struct tm tml, *gmtime(), *localtime(); |
| long time_value; |
|
| extern long time(); |
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 strcurr[80], strfor[80]; |
| |
|
| char *endptr; |
char *endptr; |
|
Line 723 char *cutl(char *blocc, char *alocc, cha
|
Line 811 char *cutl(char *blocc, char *alocc, cha
|
| gives blocc="abcdef2ghi" and alocc="j". |
gives blocc="abcdef2ghi" and alocc="j". |
| If occ is not found blocc is null and alocc is equal to in. Returns blocc |
If occ is not found blocc is null and alocc is equal to in. Returns blocc |
| */ |
*/ |
| char *s, *t, *bl; |
char *s, *t; |
| t=in;s=in; |
t=in;s=in; |
| while ((*in != occ) && (*in != '\0')){ |
while ((*in != occ) && (*in != '\0')){ |
| *alocc++ = *in++; |
*alocc++ = *in++; |
|
Line 807 int nbocc(char *s, char occ)
|
Line 895 int nbocc(char *s, char occ)
|
| /* } */ |
/* } */ |
| /* } */ |
/* } */ |
| |
|
| |
#ifdef _WIN32 |
| |
char * strsep(char **pp, const char *delim) |
| |
{ |
| |
char *p, *q; |
| |
|
| |
if ((p = *pp) == NULL) |
| |
return 0; |
| |
if ((q = strpbrk (p, delim)) != NULL) |
| |
{ |
| |
*pp = q + 1; |
| |
*q = '\0'; |
| |
} |
| |
else |
| |
*pp = 0; |
| |
return p; |
| |
} |
| |
#endif |
| |
|
| /********************** nrerror ********************/ |
/********************** nrerror ********************/ |
| |
|
| void nrerror(char error_text[]) |
void nrerror(char error_text[]) |
|
Line 1006 char *subdirf3(char fileres[], char *pre
|
Line 1112 char *subdirf3(char fileres[], char *pre
|
| return tmpout; |
return tmpout; |
| } |
} |
| |
|
| |
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; |
| |
} |
| |
|
| /***************** f1dim *************************/ |
/***************** f1dim *************************/ |
| extern int ncom; |
extern int ncom; |
| extern double *pcom,*xicom; |
extern double *pcom,*xicom; |
|
Line 1029 double brent(double ax, double bx, doubl
|
Line 1148 double brent(double ax, double bx, doubl
|
| { |
{ |
| int iter; |
int iter; |
| double a,b,d,etemp; |
double a,b,d,etemp; |
| double fu,fv,fw,fx; |
double fu=0,fv,fw,fx; |
| double ftemp; |
double ftemp; |
| double p,q,r,tol1,tol2,u,v,w,x,xm; |
double p,q,r,tol1,tol2,u,v,w,x,xm; |
| double e=0.0; |
double e=0.0; |
|
Line 1044 double brent(double ax, double bx, doubl
|
Line 1163 double brent(double ax, double bx, doubl
|
| /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/ |
/* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret)))*/ |
| printf(".");fflush(stdout); |
printf(".");fflush(stdout); |
| fprintf(ficlog,".");fflush(ficlog); |
fprintf(ficlog,".");fflush(ficlog); |
| #ifdef DEBUG |
#ifdef DEBUGBRENT |
| printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol); |
printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol); |
| fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol); |
fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol); |
| /* if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */ |
/* if ((fabs(x-xm) <= (tol2-0.5*(b-a)))||(2.0*fabs(fu-ftemp) <= ftol*1.e-2*(fabs(fu)+fabs(ftemp)))) { */ |
|
Line 1114 void mnbrak(double *ax, double *bx, doub
|
Line 1233 void mnbrak(double *ax, double *bx, doub
|
| } |
} |
| *cx=(*bx)+GOLD*(*bx-*ax); |
*cx=(*bx)+GOLD*(*bx-*ax); |
| *fc=(*func)(*cx); |
*fc=(*func)(*cx); |
| while (*fb > *fc) { |
while (*fb > *fc) { /* Declining fa, fb, fc */ |
| r=(*bx-*ax)*(*fb-*fc); |
r=(*bx-*ax)*(*fb-*fc); |
| q=(*bx-*cx)*(*fb-*fa); |
q=(*bx-*cx)*(*fb-*fa); |
| u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ |
u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ |
| (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); |
(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); |
ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */ |
| if ((*bx-u)*(u-*cx) > 0.0) { |
if ((*bx-u)*(u-*cx) > 0.0) { /* if u between b and c */ |
| fu=(*func)(u); |
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); |
fu=(*func)(u); |
| if (fu < *fc) { |
if (fu < *fc) { |
| SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) |
SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) |
| SHFT(*fb,*fc,fu,(*func)(u)) |
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; |
u=ulim; |
| fu=(*func)(u); |
fu=(*func)(u); |
| } else { |
} else { |
|
Line 1141 void mnbrak(double *ax, double *bx, doub
|
Line 1268 void mnbrak(double *ax, double *bx, doub
|
| } |
} |
| |
|
| /*************** linmin ************************/ |
/*************** 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; |
int ncom; |
| double *pcom,*xicom; |
double *pcom,*xicom; |
| double (*nrfunc)(double []); |
double (*nrfunc)(double []); |
|
Line 1167 void linmin(double p[], double xi[], int
|
Line 1298 void linmin(double p[], double xi[], int
|
| } |
} |
| ax=0.0; |
ax=0.0; |
| xx=1.0; |
xx=1.0; |
| mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); |
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); |
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */ |
| #ifdef DEBUG |
#ifdef DEBUG |
| printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); |
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); |
| fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); |
fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); |
|
Line 1181 void linmin(double p[], double xi[], int
|
Line 1312 void linmin(double p[], double xi[], int
|
| free_vector(pcom,1,n); |
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 ************************/ |
/*************** 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, |
void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, |
| double (*func)(double [])) |
double (*func)(double [])) |
| { |
{ |
|
Line 1212 void powell(double p[], double **xi, int
|
Line 1339 void powell(double p[], double **xi, int
|
| xits=vector(1,n); |
xits=vector(1,n); |
| *fret=(*func)(p); |
*fret=(*func)(p); |
| for (j=1;j<=n;j++) pt[j]=p[j]; |
for (j=1;j<=n;j++) pt[j]=p[j]; |
| |
rcurr_time = time(NULL); |
| for (*iter=1;;++(*iter)) { |
for (*iter=1;;++(*iter)) { |
| fp=(*fret); |
fp=(*fret); |
| ibig=0; |
ibig=0; |
| del=0.0; |
del=0.0; |
| last_time=curr_time; |
rlast_time=rcurr_time; |
| (void) gettimeofday(&curr_time,&tzp); |
/* (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); |
rcurr_time = time(NULL); |
| 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); |
curr_time = *localtime(&rcurr_time); |
| /* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tv_sec-start_time.tv_sec); */ |
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++) { |
for (i=1;i<=n;i++) { |
| printf(" %d %.12f",i, p[i]); |
printf(" %d %.12f",i, p[i]); |
| fprintf(ficlog," %d %.12lf",i, p[i]); |
fprintf(ficlog," %d %.12lf",i, p[i]); |
|
Line 1230 void powell(double p[], double **xi, int
|
Line 1360 void powell(double p[], double **xi, int
|
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| fprintf(ficrespow,"\n");fflush(ficrespow); |
fprintf(ficrespow,"\n");fflush(ficrespow); |
| if(*iter <=3){ |
if(*iter <=3){ |
| tm = *localtime(&curr_time.tv_sec); |
tml = *localtime(&rcurr_time); |
| strcpy(strcurr,asctime(&tm)); |
strcpy(strcurr,asctime(&tml)); |
| /* asctime_r(&tm,strcurr); */ |
rforecast_time=rcurr_time; |
| forecast_time=curr_time; |
|
| itmp = strlen(strcurr); |
itmp = strlen(strcurr); |
| if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ |
if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */ |
| strcurr[itmp-1]='\0'; |
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); |
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,curr_time.tv_sec-last_time.tv_sec); |
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){ |
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); |
rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time); |
| tmf = *localtime(&forecast_time.tv_sec); |
forecast_time = *localtime(&rforecast_time); |
| /* asctime_r(&tmf,strfor); */ |
strcpy(strfor,asctime(&forecast_time)); |
| strcpy(strfor,asctime(&tmf)); |
|
| itmp = strlen(strfor); |
itmp = strlen(strfor); |
| if(strfor[itmp-1]=='\n') |
if(strfor[itmp-1]=='\n') |
| strfor[itmp-1]='\0'; |
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); |
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(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(rforecast_time-rcurr_time,tmpout),strfor,strcurr); |
| } |
} |
| } |
} |
| for (i=1;i<=n;i++) { |
for (i=1;i<=n;i++) { |
|
Line 1274 void powell(double p[], double **xi, int
|
Line 1402 void powell(double p[], double **xi, int
|
| fprintf(ficlog," x(%d)=%.12e",j,xit[j]); |
fprintf(ficlog," x(%d)=%.12e",j,xit[j]); |
| } |
} |
| for(j=1;j<=n;j++) { |
for(j=1;j<=n;j++) { |
| printf(" p=%.12e",p[j]); |
printf(" p(%d)=%.12e",j,p[j]); |
| fprintf(ficlog," p=%.12e",p[j]); |
fprintf(ficlog," p(%d)=%.12e",j,p[j]); |
| } |
} |
| printf("\n"); |
printf("\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| #endif |
#endif |
| } |
} /* end i */ |
| if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { |
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { |
| #ifdef DEBUG |
#ifdef DEBUG |
| int k[2],l; |
int k[2],l; |
|
Line 1313 void powell(double p[], double **xi, int
|
Line 1441 void powell(double p[], double **xi, int
|
| return; |
return; |
| } |
} |
| if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); |
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]; |
ptt[j]=2.0*p[j]-pt[j]; |
| xit[j]=p[j]-pt[j]; |
xit[j]=p[j]-pt[j]; |
| pt[j]=p[j]; |
pt[j]=p[j]; |
| } |
} |
| fptt=(*func)(ptt); |
fptt=(*func)(ptt); |
| if (fptt < fp) { |
if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */ |
| t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); |
/* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */ |
| if (t < 0.0) { |
/* From x1 (P0) distance of x2 is at h and x3 is 2h */ |
| linmin(p,xit,n,fret,func); |
/* 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++) { |
for (j=1;j<=n;j++) { |
| xi[j][ibig]=xi[j][n]; |
xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */ |
| xi[j][n]=xit[j]; |
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 |
#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); |
|
| for(j=1;j<=n;j++){ |
for(j=1;j<=n;j++){ |
| printf(" %.12e",xit[j]); |
printf(" %.12e",xit[j]); |
| fprintf(ficlog," %.12e",xit[j]); |
fprintf(ficlog," %.12e",xit[j]); |
|
Line 1337 void powell(double p[], double **xi, int
|
Line 1488 void powell(double p[], double **xi, int
|
| printf("\n"); |
printf("\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| #endif |
#endif |
| } |
} /* end of t negative */ |
| } |
} /* end if (fptt < fp) */ |
| } |
} |
| } |
} |
| |
|
|
Line 1568 double ***hpxij(double ***po, int nhstep
|
Line 1719 double ***hpxij(double ***po, int nhstep
|
| return po; |
return po; |
| } |
} |
| |
|
| |
#ifdef NLOPT |
| |
double myfunc(unsigned n, const double *p1, double *grad, void *pd){ |
| |
double fret; |
| |
double *xt; |
| |
int j; |
| |
myfunc_data *d2 = (myfunc_data *) pd; |
| |
/* xt = (p1-1); */ |
| |
xt=vector(1,n); |
| |
for (j=1;j<=n;j++) xt[j]=p1[j-1]; /* xt[1]=p1[0] */ |
| |
|
| |
fret=(d2->function)(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 *************/ |
/*************** log-likelihood *************/ |
| double func( double *x) |
double func( double *x) |
|
Line 1586 double func( double *x)
|
Line 1756 double func( double *x)
|
| /*for(i=1;i<imx;i++) |
/*for(i=1;i<imx;i++) |
| printf(" %d\n",s[4][i]); |
printf(" %d\n",s[4][i]); |
| */ |
*/ |
| |
|
| |
++countcallfunc; |
| |
|
| cov[1]=1.; |
cov[1]=1.; |
| |
|
| for(k=1; k<=nlstate; k++) ll[k]=0.; |
for(k=1; k<=nlstate; k++) ll[k]=0.; |
|
Line 1972 void mlikeli(FILE *ficres,double p[], in
|
Line 2145 void mlikeli(FILE *ficres,double p[], in
|
| double fret; |
double fret; |
| double fretone; /* Only one call to likelihood */ |
double fretone; /* Only one call to likelihood */ |
| /* char filerespow[FILENAMELENGTH];*/ |
/* char filerespow[FILENAMELENGTH];*/ |
| |
|
| |
#ifdef NLOPT |
| |
int creturn; |
| |
nlopt_opt opt; |
| |
/* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */ |
| |
double *lb; |
| |
double minf; /* the minimum objective value, upon return */ |
| |
double * p1; /* Shifted parameters from 0 instead of 1 */ |
| |
myfunc_data dinst, *d = &dinst; |
| |
#endif |
| |
|
| |
|
| xi=matrix(1,npar,1,npar); |
xi=matrix(1,npar,1,npar); |
| for (i=1;i<=npar;i++) |
for (i=1;i<=npar;i++) |
| for (j=1;j<=npar;j++) |
for (j=1;j<=npar;j++) |
|
Line 1988 void mlikeli(FILE *ficres,double p[], in
|
Line 2173 void mlikeli(FILE *ficres,double p[], in
|
| for(j=1;j<=nlstate+ndeath;j++) |
for(j=1;j<=nlstate+ndeath;j++) |
| if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); |
if(j!=i)fprintf(ficrespow," p%1d%1d",i,j); |
| fprintf(ficrespow,"\n"); |
fprintf(ficrespow,"\n"); |
| |
#ifdef POWELL |
| powell(p,xi,npar,ftol,&iter,&fret,func); |
powell(p,xi,npar,ftol,&iter,&fret,func); |
| |
#endif |
| |
|
| |
#ifdef NLOPT |
| |
#ifdef NEWUOA |
| |
opt = nlopt_create(NLOPT_LN_NEWUOA,npar); |
| |
#else |
| |
opt = nlopt_create(NLOPT_LN_BOBYQA,npar); |
| |
#endif |
| |
lb=vector(0,npar-1); |
| |
for (i=0;i<npar;i++) lb[i]= -HUGE_VAL; |
| |
nlopt_set_lower_bounds(opt, lb); |
| |
nlopt_set_initial_step1(opt, 0.1); |
| |
|
| |
p1= (p+1); /* p *(p+1)@8 and p *(p1)@8 are equal p1[0]=p[1] */ |
| |
d->function = 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); |
free_matrix(xi,1,npar,1,npar); |
| fclose(ficrespow); |
fclose(ficrespow); |
| printf("\n#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 = %d, -2 Log likelihood = %.12f \n",iter,func(p)); |
fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
| fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); |
fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
| |
|
| } |
} |
| |
|
|
Line 3868 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 4080 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
| before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d_2.png\">%s%d_2.png</a><br> \ |
before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d_2.png\">%s%d_2.png</a><br> \ |
| <img src=\"%s%d_2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); |
<img src=\"%s%d_2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); |
| /* Period (stable) prevalence in each health state */ |
/* Period (stable) prevalence in each health state */ |
| for(cpt=1; cpt<nlstate;cpt++){ |
for(cpt=1; cpt<=nlstate;cpt++){ |
| fprintf(fichtm,"<br>- Period (stable) prevalence in each health state : <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \ |
fprintf(fichtm,"<br>- Convergence from each state (1 to %d) to period (stable) prevalence in state %d <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \ |
| <img src=\"%s%d_%d.png\">",subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); |
<img src=\"%s%d_%d.png\">",nlstate, cpt, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); |
| } |
} |
| for(cpt=1; cpt<=nlstate;cpt++) { |
for(cpt=1; cpt<=nlstate;cpt++) { |
| fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \ |
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) : <a href=\"%s%d%d.png\">%s%d%d.png</a> <br> \ |
| <img src=\"%s%d%d.png\">",cpt,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); |
<img src=\"%s%d%d.png\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); |
| } |
} |
| } /* end i1 */ |
} /* end i1 */ |
| }/* End k1 */ |
}/* End k1 */ |
|
Line 3970 void printinggnuplot(char fileres[], cha
|
Line 4182 void printinggnuplot(char fileres[], cha
|
| strcpy(dirfileres,optionfilefiname); |
strcpy(dirfileres,optionfilefiname); |
| strcpy(optfileres,"vpl"); |
strcpy(optfileres,"vpl"); |
| /* 1eme*/ |
/* 1eme*/ |
| |
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n"); |
| for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (cpt=1; cpt<= nlstate ; cpt ++) { |
| for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
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); |
fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1); |
|
Line 3997 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 4210 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| } |
} |
| } |
} |
| /*2 eme*/ |
/*2 eme*/ |
| |
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); |
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
| fprintf(ficgp,"\nset out \"%s%d.png\" \n",subdirf2(optionfilefiname,"e"),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); |
fprintf(ficgp,"set ylabel \"Years\" \nset ter png small size 320, 240\nplot [%.f:%.f] ",ageminpar,fage); |
|
Line 4054 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 4267 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
| } |
} |
| |
|
| /* CV preval stable (period) */ |
/* CV preval stable (period) */ |
| for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { /* For each multivariate if any */ |
| for (cpt=1; cpt<=nlstate ; cpt ++) { |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
| k=3; |
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,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1); |
| fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
| set ter png small size 320, 240\n\ |
set ter png small size 320, 240\n\ |
| unset log y\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); |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
| |
for (i=1; i<= nlstate ; i ++){ |
| for (i=1; i< nlstate ; i ++) |
if(i==1) |
| fprintf(ficgp,"+$%d",k+i+1); |
fprintf(ficgp,"\"%s\"",subdirf2(fileres,"pij")); |
| fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1); |
else |
| |
fprintf(ficgp,", '' "); |
| l=3+(nlstate+ndeath)*cpt; |
l=(nlstate+ndeath)*(i-1)+1; |
| fprintf(ficgp,",\"%s\" u ($1==%d ? ($3):1/0):($%d/($%d",subdirf2(fileres,"pij"),k1,l+cpt+1,l+1); |
fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); |
| for (i=1; i< nlstate ; i ++) { |
for (j=1; j<= (nlstate-1) ; j ++) |
| l=3+(nlstate+ndeath)*cpt; |
fprintf(ficgp,"+$%d",k+l+j); |
| fprintf(ficgp,"+$%d",l+i+1); |
fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt); |
| } |
} /* nlstate */ |
| fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1); |
fprintf(ficgp,"\n"); |
| } |
} /* end cpt state*/ |
| } |
} /* end covariate */ |
| |
|
| /* proba elementaires */ |
/* proba elementaires */ |
| for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
|
Line 5146 calandcheckages(int imx, int maxwav, dou
|
Line 5360 calandcheckages(int imx, int maxwav, dou
|
| } |
} |
| else if(agev[m][i] >*agemax){ |
else if(agev[m][i] >*agemax){ |
| *agemax=agev[m][i]; |
*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]=anint[m][i]-annais[i];*/ |
| /* agev[m][i] = age[i]+2*m;*/ |
/* agev[m][i] = age[i]+2*m;*/ |
|
Line 5275 int main(int argc, char *argv[])
|
Line 5489 int main(int argc, char *argv[])
|
| /* setlocale (LC_MESSAGES, ""); */ |
/* setlocale (LC_MESSAGES, ""); */ |
| |
|
| /* gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ |
/* 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; |
curr_time=start_time; |
| tm = *localtime(&start_time.tv_sec); |
/*tml = *localtime(&start_time.tm_sec);*/ |
| tmg = *gmtime(&start_time.tv_sec); |
/* strcpy(strstart,asctime(&tml)); */ |
| strcpy(strstart,asctime(&tm)); |
strcpy(strstart,asctime(&start_time)); |
| |
|
| /* printf("Localtime (at start)=%s",strstart); */ |
/* printf("Localtime (at start)=%s",strstart); */ |
| /* tp.tv_sec = tp.tv_sec +86400; */ |
/* tp.tm_sec = tp.tm_sec +86400; */ |
| /* tm = *localtime(&start_time.tv_sec); */ |
/* tm = *localtime(&start_time.tm_sec); */ |
| /* tmg.tm_year=tmg.tm_year +dsign*dyear; */ |
/* tmg.tm_year=tmg.tm_year +dsign*dyear; */ |
| /* tmg.tm_mon=tmg.tm_mon +dsign*dmonth; */ |
/* tmg.tm_mon=tmg.tm_mon +dsign*dmonth; */ |
| /* tmg.tm_hour=tmg.tm_hour + 1; */ |
/* tmg.tm_hour=tmg.tm_hour + 1; */ |
| /* tp.tv_sec = mktime(&tmg); */ |
/* tp.tm_sec = mktime(&tmg); */ |
| /* strt=asctime(&tmg); */ |
/* strt=asctime(&tmg); */ |
| /* printf("Time(after) =%s",strstart); */ |
/* printf("Time(after) =%s",strstart); */ |
| /* (void) time (&time_value); |
/* (void) time (&time_value); |
|
Line 5308 int main(int argc, char *argv[])
|
Line 5524 int main(int argc, char *argv[])
|
| i=strlen(pathr); |
i=strlen(pathr); |
| if(pathr[i-1]=='\n') |
if(pathr[i-1]=='\n') |
| pathr[i-1]='\0'; |
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; ){ |
for (tok = pathr; tok != NULL; ){ |
| printf("Pathr |%s|\n",pathr); |
printf("Pathr |%s|\n",pathr); |
| while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); |
while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0'); |
|
Line 5363 int main(int argc, char *argv[])
|
Line 5582 int main(int argc, char *argv[])
|
| path=%s \n\ |
path=%s \n\ |
| optionfile=%s\n\ |
optionfile=%s\n\ |
| optionfilext=%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); |
printf("Local time (at start):%s",strstart); |
| fprintf(ficlog,"Local time (at start): %s",strstart); |
fprintf(ficlog,"Local time (at start): %s",strstart); |
| fflush(ficlog); |
fflush(ficlog); |
| /* (void) gettimeofday(&curr_time,&tzp); */ |
/* (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"); |
strcpy(fileres,"r"); |
|
Line 5379 int main(int argc, char *argv[])
|
Line 5598 int main(int argc, char *argv[])
|
| /*---------arguments file --------*/ |
/*---------arguments file --------*/ |
| |
|
| if((ficpar=fopen(optionfile,"r"))==NULL) { |
if((ficpar=fopen(optionfile,"r"))==NULL) { |
| printf("Problem with optionfile %s\n",optionfile); |
printf("Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno)); |
| fprintf(ficlog,"Problem with optionfile %s\n",optionfile); |
fprintf(ficlog,"Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno)); |
| fflush(ficlog); |
fflush(ficlog); |
| goto end; |
/* goto end; */ |
| |
exit(70); |
| } |
} |
| |
|
| |
|
|
Line 5901 Interval (in months) between two waves:
|
Line 6121 Interval (in months) between two waves:
|
| |
|
| #ifdef GSL |
#ifdef GSL |
| printf("GSL optimization\n"); fprintf(ficlog,"Powell\n"); |
printf("GSL optimization\n"); fprintf(ficlog,"Powell\n"); |
| #elsedef |
#else |
| printf("Powell\n"); fprintf(ficlog,"Powell\n"); |
printf("Powell\n"); fprintf(ficlog,"Powell\n"); |
| #endif |
#endif |
| strcpy(filerespow,"pow-mort"); |
strcpy(filerespow,"pow-mort"); |
|
Line 5912 Interval (in months) between two waves:
|
Line 6132 Interval (in months) between two waves:
|
| } |
} |
| #ifdef GSL |
#ifdef GSL |
| fprintf(ficrespow,"# GSL optimization\n# iter -2*LL"); |
fprintf(ficrespow,"# GSL optimization\n# iter -2*LL"); |
| #elsedef |
#else |
| fprintf(ficrespow,"# Powell\n# iter -2*LL"); |
fprintf(ficrespow,"# Powell\n# iter -2*LL"); |
| #endif |
#endif |
| /* for (i=1;i<=nlstate;i++) |
/* for (i=1;i<=nlstate;i++) |
|
Line 6474 Interval (in months) between two waves:
|
Line 6694 Interval (in months) between two waves:
|
| |
|
| for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
| oldm=oldms;savm=savms; /* Segmentation fault */ |
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 "); |
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) |
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); |
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); |
|
Line 6600 Interval (in months) between two waves:
|
Line 6821 Interval (in months) between two waves:
|
| } |
} |
| printf("See log file on %s\n",filelog); |
printf("See log file on %s\n",filelog); |
| /* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */ |
/* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */ |
| (void) gettimeofday(&end_time,&tzp); |
/*(void) gettimeofday(&end_time,&tzp);*/ |
| tm = *localtime(&end_time.tv_sec); |
rend_time = time(NULL); |
| tmg = *gmtime(&end_time.tv_sec); |
end_time = *localtime(&rend_time); |
| strcpy(strtend,asctime(&tm)); |
/* tml = *localtime(&end_time.tm_sec); */ |
| |
strcpy(strtend,asctime(&end_time)); |
| printf("Local time at start %s\nLocal time at end %s",strstart, strtend); |
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); |
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); |
printf("Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time)); |
| fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout)); |
fprintf(ficlog,"Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout)); |
| fprintf(ficlog,"Total time was %ld Sec.\n", end_time.tv_sec -start_time.tv_sec); |
fprintf(ficlog,"Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time)); |
| /* printf("Total time was %d uSec.\n", total_usecs);*/ |
/* printf("Total time was %d uSec.\n", total_usecs);*/ |
| /* if(fileappend(fichtm,optionfilehtm)){ */ |
/* if(fileappend(fichtm,optionfilehtm)){ */ |
| fprintf(fichtm,"<br>Local time at start %s<br>Local time at end %s<br>\n</body></html>",strstart, strtend); |
fprintf(fichtm,"<br>Local time at start %s<br>Local time at end %s<br>\n</body></html>",strstart, strtend); |
|
Line 6629 Interval (in months) between two waves:
|
Line 6851 Interval (in months) between two waves:
|
| printf("Current directory %s!\n",pathcd); |
printf("Current directory %s!\n",pathcd); |
| /*strcat(plotcmd,CHARSEPARATOR);*/ |
/*strcat(plotcmd,CHARSEPARATOR);*/ |
| sprintf(plotcmd,"gnuplot"); |
sprintf(plotcmd,"gnuplot"); |
| #ifndef UNIX |
#ifdef _WIN32 |
| sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); |
sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach); |
| #endif |
#endif |
| if(!stat(plotcmd,&info)){ |
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)){ |
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 |
}else |
| strcpy(pplotcmd,plotcmd); |
strcpy(pplotcmd,plotcmd); |
| #ifdef UNIX |
#ifdef __unix |
| strcpy(plotcmd,GNUPLOTPROGRAM); |
strcpy(plotcmd,GNUPLOTPROGRAM); |
| if(!stat(plotcmd,&info)){ |
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 |
}else |
| strcpy(pplotcmd,plotcmd); |
strcpy(pplotcmd,plotcmd); |
| #endif |
#endif |
|
Line 6649 Interval (in months) between two waves:
|
Line 6871 Interval (in months) between two waves:
|
| strcpy(pplotcmd,plotcmd); |
strcpy(pplotcmd,plotcmd); |
| |
|
| sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot); |
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){ |
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') { |
while (z[0] != 'q') { |
| /* chdir(path); */ |
/* 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); |
scanf("%s",z); |
| /* if (z[0] == 'c') system("./imach"); */ |
/* if (z[0] == 'c') system("./imach"); */ |
| if (z[0] == 'e') { |
if (z[0] == 'e') { |
| printf("Starting browser with: %s",optionfilehtm);fflush(stdout); |
#ifdef __APPLE__ |
| system(optionfilehtm); |
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] == 'g') system(plotcmd); |
| else if (z[0] == 'q') exit(0); |
else if (z[0] == 'q') exit(0); |
|
Line 6673 Interval (in months) between two waves:
|
Line 6906 Interval (in months) between two waves:
|
| scanf("%s",z); |
scanf("%s",z); |
| } |
} |
| } |
} |
| |
|
| |
|
| |
|
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src