|
version 1.154, 2014/06/20 17:32:08
|
version 1.166, 2014/12/22 11:40:47
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
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 |
Revision 1.154 2014/06/20 17:32:08 brouard |
| Summary: Outputs now all graphs of convergence to period prevalence |
Summary: Outputs now all graphs of convergence to period prevalence |
| |
|
|
Line 451
|
Line 505
|
| end |
end |
| */ |
*/ |
| |
|
| |
#define POWELL /* Instead of NLOPT */ |
| |
|
| |
|
| |
|
| #include <math.h> |
#include <math.h> |
| #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 502 extern int errno;
|
Line 569 extern int errno;
|
| #define YEARM 12. /**< Number of months per year */ |
#define YEARM 12. /**< Number of months per year */ |
| #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.98nU, 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 547 int **mw; /* mw[mi][i] is number of the
|
Line 614 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 590 char popfile[FILENAMELENGTH];
|
Line 658 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 625 static double maxarg1,maxarg2;
|
Line 698 static double maxarg1,maxarg2;
|
| |
|
| #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a)) |
#define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a)) |
| #define rint(a) floor(a+0.5) |
#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))<mytinydouble) */ |
| |
/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/mynrutils_8h-source.html */ |
| |
/* static double dsqrarg; */ |
| |
/* #define DSQR(a) (DEQUAL((dsqrarg=(a)),0.0) ? 0.0 : dsqrarg*dsqrarg) */ |
| static double sqrarg; |
static double sqrarg; |
| #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg) |
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 :sqrarg*sqrarg) |
| #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} |
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;} |
|
Line 750 char *cutl(char *blocc, char *alocc, cha
|
Line 828 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 834 int nbocc(char *s, char occ)
|
Line 912 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 1033 char *subdirf3(char fileres[], char *pre
|
Line 1129 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 1056 double brent(double ax, double bx, doubl
|
Line 1165 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=0.; |
| 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 1071 double brent(double ax, double bx, doubl
|
Line 1180 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 1141 void mnbrak(double *ax, double *bx, doub
|
Line 1250 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 1168 void mnbrak(double *ax, double *bx, doub
|
Line 1285 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 1194 void linmin(double p[], double xi[], int
|
Line 1315 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 1208 void linmin(double p[], double xi[], int
|
Line 1329 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 1239 void powell(double p[], double **xi, int
|
Line 1356 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 1257 void powell(double p[], double **xi, int
|
Line 1377 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++) { |
| for (j=1;j<=n;j++) xit[j]=xi[j][i]; |
for (j=1;j<=n;j++) xit[j]=xi[j][i]; |
| fptt=(*fret); |
fptt=(*fret); |
| #ifdef DEBUG |
#ifdef DEBUG |
| printf("fret=%lf \n",*fret); |
printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); |
| fprintf(ficlog,"fret=%lf \n",*fret); |
fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); |
| #endif |
#endif |
| printf("%d",i);fflush(stdout); |
printf("%d",i);fflush(stdout); |
| fprintf(ficlog,"%d",i);fflush(ficlog); |
fprintf(ficlog,"%d",i);fflush(ficlog); |
|
Line 1301 void powell(double p[], double **xi, int
|
Line 1419 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 1340 void powell(double p[], double **xi, int
|
Line 1458 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); |
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); |
fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); |
|
Line 1364 void powell(double p[], double **xi, int
|
Line 1507 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 1453 double **pmij(double **ps, double *cov,
|
Line 1596 double **pmij(double **ps, double *cov,
|
| */ |
*/ |
| double s1, lnpijopii; |
double s1, lnpijopii; |
| /*double t34;*/ |
/*double t34;*/ |
| int i,j,j1, nc, ii, jj; |
int i,j, nc, ii, jj; |
| |
|
| for(i=1; i<= nlstate; i++){ |
for(i=1; i<= nlstate; i++){ |
| for(j=1; j<i;j++){ |
for(j=1; j<i;j++){ |
|
Line 1595 double ***hpxij(double ***po, int nhstep
|
Line 1738 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 1613 double func( double *x)
|
Line 1775 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 1994 void likelione(FILE *ficres,double p[],
|
Line 2159 void likelione(FILE *ficres,double p[],
|
| |
|
| void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double [])) |
void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (*func)(double [])) |
| { |
{ |
| int i,j, iter; |
int i,j, iter=0; |
| double **xi; |
double **xi; |
| 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 2015 void mlikeli(FILE *ficres,double p[], in
|
Line 2192 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 2031 void hesscov(double **matcov, double p[]
|
Line 2235 void hesscov(double **matcov, double p[]
|
| { |
{ |
| double **a,**y,*x,pd; |
double **a,**y,*x,pd; |
| double **hess; |
double **hess; |
| int i, j,jk; |
int i, j; |
| int *indx; |
int *indx; |
| |
|
| double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar); |
double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar); |
|
Line 2165 double hessii(double x[], double delta,
|
Line 2369 double hessii(double x[], double delta,
|
| k=kmax; |
k=kmax; |
| } |
} |
| else if((k1 >khi/nkhif) || (k2 >khi/nkhif)){ /* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. */ |
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)){ |
else if((k1 >khi/nkhi) || (k2 >khi/nkhi)){ |
| delts=delt; |
delts=delt; |
|
Line 2180 double hessii(double x[], double delta,
|
Line 2384 double hessii(double x[], double delta,
|
| double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) |
double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar) |
| { |
{ |
| int i; |
int i; |
| int l=1, l1, lmax=20; |
int l=1, lmax=20; |
| double k1,k2,k3,k4,res,fx; |
double k1,k2,k3,k4,res,fx; |
| double p2[MAXPARM+1]; |
double p2[MAXPARM+1]; |
| int k; |
int k; |
|
Line 2295 void pstamp(FILE *fichier)
|
Line 2499 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[]) |
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 */ |
{ /* Some frequencies */ |
| |
|
| int i, m, jk, k1,i1, j1, bool, z1,j; |
int i, m, jk, j1, bool, z1,j; |
| int first; |
int first; |
| double ***freq; /* Frequencies */ |
double ***freq; /* Frequencies */ |
| double *pp, **prop; |
double *pp, **prop; |
|
Line 2476 void prevalence(double ***probs, double
|
Line 2680 void prevalence(double ***probs, double
|
| We still use firstpass and lastpass as another selection. |
We still use firstpass and lastpass as another selection. |
| */ |
*/ |
| |
|
| int i, m, jk, k1, i1, j1, bool, z1,j; |
int i, m, jk, j1, bool, z1,j; |
| double ***freq; /* Frequencies */ |
|
| double *pp, **prop; |
double **prop; |
| double pos,posprop; |
double posprop; |
| double y2; /* in fractional years */ |
double y2; /* in fractional years */ |
| int iagemin, iagemax; |
int iagemin, iagemax; |
| int first; /** to stop verbosity which is redirected to log file */ |
int first; /** to stop verbosity which is redirected to log file */ |
|
Line 2570 void concatwav(int wav[], int **dh, int
|
Line 2774 void concatwav(int wav[], int **dh, int
|
| int j, k=0,jk, ju, jl; |
int j, k=0,jk, ju, jl; |
| double sum=0.; |
double sum=0.; |
| first=0; |
first=0; |
| jmin=1e+5; |
jmin=100000; |
| jmax=-1; |
jmax=-1; |
| jmean=0.; |
jmean=0.; |
| for(i=1; i<=imx; i++){ |
for(i=1; i<=imx; i++){ |
|
Line 2810 void evsij(double ***eij, double x[], in
|
Line 3014 void evsij(double ***eij, double x[], in
|
| |
|
| { |
{ |
| /* Health expectancies, no variances */ |
/* 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 */ |
int nhstepma, nstepma; /* Decreasing with age */ |
| double age, agelim, hf; |
double age, agelim, hf; |
| double ***p3mat; |
double ***p3mat; |
|
Line 3132 void varevsij(char optionfilefiname[], d
|
Line 3336 void varevsij(char optionfilefiname[], d
|
| double **dnewm,**doldm; |
double **dnewm,**doldm; |
| double **dnewmp,**doldmp; |
double **dnewmp,**doldmp; |
| int i, j, nhstepm, hstepm, h, nstepm ; |
int i, j, nhstepm, hstepm, h, nstepm ; |
| int k, cptcode; |
int k; |
| double *xp; |
double *xp; |
| double **gp, **gm; /* for var eij */ |
double **gp, **gm; /* for var eij */ |
| double ***gradg, ***trgradg; /*for var eij */ |
double ***gradg, ***trgradg; /*for var eij */ |
|
Line 3432 void varprevlim(char fileres[], double *
|
Line 3636 void varprevlim(char fileres[], double *
|
| { |
{ |
| /* Variance of prevalence limit */ |
/* Variance of prevalence limit */ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
| double **newm; |
|
| double **dnewm,**doldm; |
double **dnewm,**doldm; |
| int i, j, nhstepm, hstepm; |
int i, j, nhstepm, hstepm; |
| int k, cptcode; |
|
| double *xp; |
double *xp; |
| double *gp, *gm; |
double *gp, *gm; |
| double **gradg, **trgradg; |
double **gradg, **trgradg; |
|
Line 3514 void varprevlim(char fileres[], double *
|
Line 3717 void varprevlim(char fileres[], double *
|
| /************ Variance of one-step probabilities ******************/ |
/************ 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[]) |
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 k2, l2, j1, z1; |
| int k=0,l, cptcode; |
int k=0, l; |
| int first=1, first1, first2; |
int first=1, first1, first2; |
| double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; |
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp; |
| double **dnewm,**doldm; |
double **dnewm,**doldm; |
|
Line 3524 void varprob(char optionfilefiname[], do
|
Line 3727 void varprob(char optionfilefiname[], do
|
| double *gp, *gm; |
double *gp, *gm; |
| double **gradg, **trgradg; |
double **gradg, **trgradg; |
| double **mu; |
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 */ |
double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */ |
| int theta; |
int theta; |
| char fileresprob[FILENAMELENGTH]; |
char fileresprob[FILENAMELENGTH]; |
|
Line 3896 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 4099 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
| <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>- 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> \ |
fprintf(fichtm,"<br>- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \ |
| <img src=\"%s%d_%d.png\">",nlstate, cpt, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); |
<img src=\"%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++) { |
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 in each alive state (1 to %d) : <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> \ |
|
Line 3982 true period expectancies (those weighted
|
Line 4185 true period expectancies (those weighted
|
| void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ |
void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ |
| |
|
| char dirfileres[132],optfileres[132]; |
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; |
int ng=0; |
| /* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */ |
| /* printf("Problem with file %s",optionfilegnuplot); */ |
/* printf("Problem with file %s",optionfilegnuplot); */ |
|
Line 4168 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 4371 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| } /* end k2 */ |
} /* end k2 */ |
| } /* end jk */ |
} /* end jk */ |
| } /* end ng */ |
} /* end ng */ |
| avoid: |
/* avoid: */ |
| fflush(ficgp); |
fflush(ficgp); |
| } /* end gnuplot */ |
} /* end gnuplot */ |
| |
|
|
Line 4222 prevforecast(char fileres[], double anpr
|
Line 4425 prevforecast(char fileres[], double anpr
|
| dateprev1 dateprev2 range of dates during which prevalence is computed |
dateprev1 dateprev2 range of dates during which prevalence is computed |
| anproj2 year of en of projection (same day and month as proj1). |
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 yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1; |
| int *popage; |
|
| double agec; /* generic age */ |
double agec; /* generic age */ |
| double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; |
double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; |
| double *popeffectif,*popcount; |
double *popeffectif,*popcount; |
|
Line 4517 void prwizard(int ncovmodel, int nlstate
|
Line 4719 void prwizard(int ncovmodel, int nlstate
|
| |
|
| /* Wizard to print covariance matrix template */ |
/* Wizard to print covariance matrix template */ |
| |
|
| char ca[32], cb[32], cc[32]; |
char ca[32], cb[32]; |
| int i,j, k, l, li, lj, lk, ll, jj, npar, itimes; |
int i,j, k, li, lj, lk, ll, jj, npar, itimes; |
| int numlinepar; |
int numlinepar; |
| |
|
| printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
|
Line 4739 fprintf(fichtm,"<ul><li><h4>Life table</
|
Line 4941 fprintf(fichtm,"<ul><li><h4>Life table</
|
| void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ |
void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){ |
| |
|
| char dirfileres[132],optfileres[132]; |
char dirfileres[132],optfileres[132]; |
| int m,cpt,k1,i,k,j,jk,k2,k3,ij,l; |
|
| int ng; |
int ng; |
| |
|
| |
|
|
Line 4764 int readdata(char datafile[], int firsto
|
Line 4966 int readdata(char datafile[], int firsto
|
| /*-------- data file ----------*/ |
/*-------- data file ----------*/ |
| FILE *fic; |
FILE *fic; |
| char dummy[]=" "; |
char dummy[]=" "; |
| int i, j, n; |
int i=0, j=0, n=0; |
| int linei, month, year,iout; |
int linei, month, year,iout; |
| char line[MAXLINE], linetmp[MAXLINE]; |
char line[MAXLINE], linetmp[MAXLINE]; |
| char stra[80], strb[80]; |
char stra[MAXLINE], strb[MAXLINE]; |
| char *stratrunc; |
char *stratrunc; |
| int lstra; |
int lstra; |
| |
|
|
Line 4795 int readdata(char datafile[], int firsto
|
Line 4997 int readdata(char datafile[], int firsto
|
| continue; |
continue; |
| } |
} |
| trimbb(linetmp,line); /* Trims multiple blanks in line */ |
trimbb(linetmp,line); /* Trims multiple blanks in line */ |
| for (j=0; line[j]!='\0';j++){ |
strcpy(line, linetmp); |
| line[j]=linetmp[j]; |
|
| } |
|
| |
|
| |
|
| for (j=maxwav;j>=1;j--){ |
for (j=maxwav;j>=1;j--){ |
|
Line 4936 int readdata(char datafile[], int firsto
|
Line 5136 int readdata(char datafile[], int firsto
|
| fclose(fic); |
fclose(fic); |
| |
|
| return (0); |
return (0); |
| endread: |
/* endread: */ |
| printf("Exiting readdata: "); |
printf("Exiting readdata: "); |
| fclose(fic); |
fclose(fic); |
| return (1); |
return (1); |
|
Line 4967 int decodemodel ( char model[], int last
|
Line 5167 int decodemodel ( char model[], int last
|
| */ |
*/ |
| { |
{ |
| int i, j, k, ks; |
int i, j, k, ks; |
| int i1, j1, k1, k2; |
int j1, k1, k2; |
| char modelsav[80]; |
char modelsav[80]; |
| char stra[80], strb[80], strc[80], strd[80],stre[80]; |
char stra[80], strb[80], strc[80], strd[80],stre[80]; |
| |
|
|
Line 5116 int decodemodel ( char model[], int last
|
Line 5316 int decodemodel ( char model[], int last
|
| |
|
| |
|
| return (0); /* with covar[new additional covariate if product] and Tage if age */ |
return (0); /* with covar[new additional covariate if product] and Tage if age */ |
| endread: |
/*endread:*/ |
| printf("Exiting decodemodel: "); |
printf("Exiting decodemodel: "); |
| return (1); |
return (1); |
| } |
} |
|
Line 5175 calandcheckages(int imx, int maxwav, dou
|
Line 5375 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 5213 calandcheckages(int imx, int maxwav, dou
|
Line 5413 calandcheckages(int imx, int maxwav, dou
|
| fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax); |
fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, *agemin, *agemax); |
| |
|
| return (0); |
return (0); |
| endread: |
/* endread:*/ |
| printf("Exiting calandcheckages: "); |
printf("Exiting calandcheckages: "); |
| return (1); |
return (1); |
| } |
} |
|
Line 5233 int main(int argc, char *argv[])
|
Line 5433 int main(int argc, char *argv[])
|
| double ssval; |
double ssval; |
| #endif |
#endif |
| int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); |
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 i,j, k, n=MAXN,iter=0,m,size=100, cptcod; |
| int linei, month, year,iout; |
|
| int jj, ll, li, lj, lk, imk; |
int jj, ll, li, lj, lk; |
| int numlinepar=0; /* Current linenumber of parameter file */ |
int numlinepar=0; /* Current linenumber of parameter file */ |
| int itimes; |
int itimes; |
| int NDIM=2; |
int NDIM=2; |
| int vpopbased=0; |
int vpopbased=0; |
| |
|
| char ca[32], cb[32], cc[32]; |
char ca[32], cb[32]; |
| /* FILE *fichtm; *//* Html File */ |
/* FILE *fichtm; *//* Html File */ |
| /* FILE *ficgp;*/ /*Gnuplot File */ |
/* FILE *ficgp;*/ /*Gnuplot File */ |
| struct stat info; |
struct stat info; |
| double agedeb, agefin,hf; |
double agedeb; |
| double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; |
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; |
| |
|
| double fret; |
double fret; |
| double **xi,tmp,delta; |
|
| |
|
| double dum; /* Dummy variable */ |
double dum; /* Dummy variable */ |
| double ***p3mat; |
double ***p3mat; |
| double ***mobaverage; |
double ***mobaverage; |
| int *indx; |
|
| char line[MAXLINE], linepar[MAXLINE]; |
char line[MAXLINE]; |
| char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; |
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; |
| char pathr[MAXLINE], pathimach[MAXLINE]; |
char pathr[MAXLINE], pathimach[MAXLINE]; |
| char **bp, *tok, *val; /* pathtot */ |
char *tok, *val; /* pathtot */ |
| int firstobs=1, lastobs=10; |
int firstobs=1, lastobs=10; |
| int sdeb, sfin; /* Status at beginning and end */ |
int c, h , cpt; |
| int c, h , cpt,l; |
int jl; |
| int ju,jl, mi; |
int i1, j1, jk, stepsize; |
| int i1,j1, jk,aa,bb, stepsize, ij; |
int *tab; |
| int jnais,jdc,jint4,jint1,jint2,jint3,*tab; |
|
| int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ |
int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */ |
| int mobilav=0,popforecast=0; |
int mobilav=0,popforecast=0; |
| int hstepm, nhstepm; |
int hstepm, nhstepm; |
|
Line 5273 int main(int argc, char *argv[])
|
Line 5470 int main(int argc, char *argv[])
|
| double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000; |
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 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 ftolpl=FTOL; |
| double **prlim; |
double **prlim; |
| double ***param; /* Matrix of parameters */ |
double ***param; /* Matrix of parameters */ |
|
Line 5284 int main(int argc, char *argv[])
|
Line 5481 int main(int argc, char *argv[])
|
| double ***eij, ***vareij; |
double ***eij, ***vareij; |
| double **varpl; /* Variances of prevalence limits by age */ |
double **varpl; /* Variances of prevalence limits by age */ |
| double *epj, vepp; |
double *epj, vepp; |
| double kk1, kk2; |
|
| double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; |
double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000; |
| double **ximort; |
double **ximort; |
| char *alph[]={"a","a","b","c","d","e"}, str[4]="1234"; |
char *alph[]={"a","a","b","c","d","e"}, str[4]="1234"; |
| int *dcwave; |
int *dcwave; |
| |
|
| char z[1]="c", occ; |
char z[1]="c"; |
| |
|
| /*char *strt;*/ |
/*char *strt;*/ |
| char strtend[80]; |
char strtend[80]; |
| |
|
| long total_usecs; |
|
| |
|
| /* setlocale (LC_ALL, ""); */ |
/* setlocale (LC_ALL, ""); */ |
| /* bindtextdomain (PACKAGE, LOCALEDIR); */ |
/* bindtextdomain (PACKAGE, LOCALEDIR); */ |
| /* textdomain (PACKAGE); */ |
/* textdomain (PACKAGE); */ |
|
Line 5304 int main(int argc, char *argv[])
|
Line 5500 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 5337 int main(int argc, char *argv[])
|
Line 5535 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 5392 int main(int argc, char *argv[])
|
Line 5593 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 5408 int main(int argc, char *argv[])
|
Line 5609 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 with errno=%s\n",optionfile,strerror(errno)); |
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)); |
fprintf(ficlog,"Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno)); |
| fflush(ficlog); |
fflush(ficlog); |
| /* goto end; */ |
/* goto end; */ |
| exit(70); |
exit(70); |
|
Line 5562 run imach with mle=-1 to get a correct t
|
Line 5763 run imach with mle=-1 to get a correct t
|
| for(i=1; i <=nlstate; i++){ |
for(i=1; i <=nlstate; i++){ |
| for(j=1; j <=nlstate+ndeath-1; j++){ |
for(j=1; j <=nlstate+ndeath-1; j++){ |
| fscanf(ficpar,"%1d%1d",&i1,&j1); |
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); |
printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1); |
| exit(1); |
exit(1); |
| } |
} |
|
Line 5931 Interval (in months) between two waves:
|
Line 6132 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 5942 Interval (in months) between two waves:
|
Line 6143 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 6504 Interval (in months) between two waves:
|
Line 6705 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 6596 Interval (in months) between two waves:
|
Line 6798 Interval (in months) between two waves:
|
| if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| } /* mle==-3 arrives here for freeing */ |
} /* mle==-3 arrives here for freeing */ |
| endfree: |
/* endfree:*/ |
| free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */ |
| free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
|
Line 6630 Interval (in months) between two waves:
|
Line 6832 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 6659 Interval (in months) between two waves:
|
Line 6862 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 6679 Interval (in months) between two waves:
|
Line 6882 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("gnuplot command might not be in your path: %s, err=%d\n", plotcmd, outcmd); |
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"); |
printf("\n Trying if gnuplot resides on the same directory that IMaCh\n"); |
| sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); |
sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot); |
| if((outcmd=system(plotcmd)) != 0) |
if((outcmd=system(plotcmd)) != 0) |
| printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); |
printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd); |
| } |
} |
| printf(" Successul, please wait..."); |
printf(" Successful, please wait..."); |
| while (z[0] != 'q') { |
while (z[0] != 'q') { |
| /* chdir(path); */ |
/* chdir(path); */ |
| printf("\nType e to edit results with your browser, g to graph again and q for exit: "); |
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') { |
| #ifdef OSX |
#ifdef __APPLE__ |
| sprintf(pplotcmd, "open %s", optionfilehtm); |
sprintf(pplotcmd, "open %s", optionfilehtm); |
| |
#elif __linux |
| |
sprintf(pplotcmd, "xdg-open %s", optionfilehtm); |
| #else |
#else |
| sprintf(pplotcmd, "%s", optionfilehtm); |
sprintf(pplotcmd, "%s", optionfilehtm); |
| #endif |
#endif |
|
Line 6712 Interval (in months) between two waves:
|
Line 6917 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