version 1.164, 2014/12/16 10:52:11
|
version 1.184, 2015/03/11 11:52:39
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.184 2015/03/11 11:52:39 brouard |
|
Summary: Back from Windows 8. Intel Compiler |
|
|
|
Revision 1.183 2015/03/10 20:34:32 brouard |
|
Summary: 0.98q0, trying with directest, mnbrak fixed |
|
|
|
We use directest instead of original Powell test; probably no |
|
incidence on the results, but better justifications; |
|
We fixed Numerical Recipes mnbrak routine which was wrong and gave |
|
wrong results. |
|
|
|
Revision 1.182 2015/02/12 08:19:57 brouard |
|
Summary: Trying to keep directest which seems simpler and more general |
|
Author: Nicolas Brouard |
|
|
|
Revision 1.181 2015/02/11 23:22:24 brouard |
|
Summary: Comments on Powell added |
|
|
|
Author: |
|
|
|
Revision 1.180 2015/02/11 17:33:45 brouard |
|
Summary: Finishing move from main to function (hpijx and prevalence_limit) |
|
|
|
Revision 1.179 2015/01/04 09:57:06 brouard |
|
Summary: back to OS/X |
|
|
|
Revision 1.178 2015/01/04 09:35:48 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.177 2015/01/03 18:40:56 brouard |
|
Summary: Still testing ilc32 on OSX |
|
|
|
Revision 1.176 2015/01/03 16:45:04 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.175 2015/01/03 16:33:42 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.174 2015/01/03 16:15:49 brouard |
|
Summary: Still in cross-compilation |
|
|
|
Revision 1.173 2015/01/03 12:06:26 brouard |
|
Summary: trying to detect cross-compilation |
|
|
|
Revision 1.172 2014/12/27 12:07:47 brouard |
|
Summary: Back from Visual Studio and Intel, options for compiling for Windows XP |
|
|
|
Revision 1.171 2014/12/23 13:26:59 brouard |
|
Summary: Back from Visual C |
|
|
|
Still problem with utsname.h on Windows |
|
|
|
Revision 1.170 2014/12/23 11:17:12 brouard |
|
Summary: Cleaning some \%% back to %% |
|
|
|
The escape was mandatory for a specific compiler (which one?), but too many warnings. |
|
|
|
Revision 1.169 2014/12/22 23:08:31 brouard |
|
Summary: 0.98p |
|
|
|
Outputs some informations on compiler used, OS etc. Testing on different platforms. |
|
|
|
Revision 1.168 2014/12/22 15:17:42 brouard |
|
Summary: update |
|
|
|
Revision 1.167 2014/12/22 13:50:56 brouard |
|
Summary: Testing uname and compiler version and if compiled 32 or 64 |
|
|
|
Testing on Linux 64 |
|
|
|
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 |
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 |
Summary: Merging with Visual C after suppressing some warnings for unused variables. Also fixing Saito's bug 0.98Xn |
|
|
Line 497
|
Line 575
|
end |
end |
*/ |
*/ |
|
|
|
#define POWELL /* Instead of NLOPT */ |
|
/* #define POWELLORIGINAL */ /* Don't use Directest to decide new direction but original Powell test */ |
|
/* #define MNBRAKORIGINAL */ /* Don't use mnbrak fix */ |
|
|
|
|
|
|
#include <math.h> |
#include <math.h> |
#include <stdio.h> |
#include <stdio.h> |
#include <stdlib.h> |
#include <stdlib.h> |
Line 507
|
Line 586
|
|
|
#ifdef _WIN32 |
#ifdef _WIN32 |
#include <io.h> |
#include <io.h> |
|
#include <windows.h> |
|
#include <tchar.h> |
#else |
#else |
#include <unistd.h> |
#include <unistd.h> |
#endif |
#endif |
|
|
#include <limits.h> |
#include <limits.h> |
#include <sys/types.h> |
#include <sys/types.h> |
|
|
|
#if defined(__GNUC__) |
|
#include <sys/utsname.h> /* Doesn't work on Windows */ |
|
#endif |
|
|
#include <sys/stat.h> |
#include <sys/stat.h> |
#include <errno.h> |
#include <errno.h> |
/* extern int errno; */ |
/* extern int errno; */ |
Line 531
|
Line 617
|
#include <gsl/gsl_multimin.h> |
#include <gsl/gsl_multimin.h> |
#endif |
#endif |
|
|
|
|
#ifdef NLOPT |
#ifdef NLOPT |
#include <nlopt.h> |
#include <nlopt.h> |
typedef struct { |
typedef struct { |
Line 576 typedef struct {
|
Line 663 typedef struct {
|
/* $Id$ */ |
/* $Id$ */ |
/* $State$ */ |
/* $State$ */ |
|
|
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 version[]="Imach version 0.98q0, March 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; |
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 691 static double maxarg1,maxarg2;
|
Line 778 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 746 static int split( char *path, char *dirc
|
Line 838 static int split( char *path, char *dirc
|
printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ |
printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ |
/* get current working directory */ |
/* get current working directory */ |
/* extern char* getcwd ( char *buf , int len);*/ |
/* extern char* getcwd ( char *buf , int len);*/ |
if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { |
#ifdef WIN32 |
|
if (_getcwd( dirc, FILENAME_MAX ) == NULL ) { |
|
#else |
|
if (getcwd(dirc, FILENAME_MAX) == NULL) { |
|
#endif |
return( GLOCK_ERROR_GETCWD ); |
return( GLOCK_ERROR_GETCWD ); |
} |
} |
/* got dirc from getcwd*/ |
/* got dirc from getcwd*/ |
Line 1203 double brent(double ax, double bx, doubl
|
Line 1299 double brent(double ax, double bx, doubl
|
if (fu <= fx) { |
if (fu <= fx) { |
if (u >= x) a=x; else b=x; |
if (u >= x) a=x; else b=x; |
SHFT(v,w,x,u) |
SHFT(v,w,x,u) |
SHFT(fv,fw,fx,fu) |
SHFT(fv,fw,fx,fu) |
} else { |
} else { |
if (u < x) a=u; else b=u; |
if (u < x) a=u; else b=u; |
if (fu <= fw || w == x) { |
if (fu <= fw || w == x) { |
v=w; |
v=w; |
w=u; |
w=u; |
fv=fw; |
fv=fw; |
fw=fu; |
fw=fu; |
} else if (fu <= fv || v == x || v == w) { |
} else if (fu <= fv || v == x || v == w) { |
v=u; |
v=u; |
fv=fu; |
fv=fu; |
} |
} |
} |
} |
} |
} |
nrerror("Too many iterations in brent"); |
nrerror("Too many iterations in brent"); |
*xmin=x; |
*xmin=x; |
Line 1226 double brent(double ax, double bx, doubl
|
Line 1322 double brent(double ax, double bx, doubl
|
|
|
void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc, |
void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc, |
double (*func)(double)) |
double (*func)(double)) |
{ |
{ /* Given a function func , and given distinct initial points ax and bx , this routine searches in |
|
the downhill direction (defined by the function as evaluated at the initial points) and returns |
|
new points ax , bx , cx that bracket a minimum of the function. Also returned are the function |
|
values at the three points, fa, fb , and fc such that fa > fb and fb < fc. |
|
*/ |
double ulim,u,r,q, dum; |
double ulim,u,r,q, dum; |
double fu; |
double fu; |
|
|
Line 1234 void mnbrak(double *ax, double *bx, doub
|
Line 1334 void mnbrak(double *ax, double *bx, doub
|
*fb=(*func)(*bx); |
*fb=(*func)(*bx); |
if (*fb > *fa) { |
if (*fb > *fa) { |
SHFT(dum,*ax,*bx,dum) |
SHFT(dum,*ax,*bx,dum) |
SHFT(dum,*fb,*fa,dum) |
SHFT(dum,*fb,*fa,dum) |
} |
} |
*cx=(*bx)+GOLD*(*bx-*ax); |
*cx=(*bx)+GOLD*(*bx-*ax); |
*fc=(*func)(*cx); |
*fc=(*func)(*cx); |
while (*fb > *fc) { /* Declining fa, fb, fc */ |
#ifdef DEBUG |
|
printf("mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc); |
|
fprintf(ficlog,"mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc); |
|
#endif |
|
while (*fb > *fc) { /* Declining a,b,c with 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)); /* Minimum abscisse of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */ |
(2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscissa of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */ |
ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */ |
ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscissa where function should be evaluated */ |
if ((*bx-u)*(u-*cx) > 0.0) { /* if u between b and c */ |
if ((*bx-u)*(u-*cx) > 0.0) { /* if u_p is between b and c */ |
fu=(*func)(u); |
fu=(*func)(u); |
#ifdef DEBUG |
#ifdef DEBUG |
/* f(x)=A(x-u)**2+f(u) */ |
/* f(x)=A(x-u)**2+f(u) */ |
Line 1253 void mnbrak(double *ax, double *bx, doub
|
Line 1357 void mnbrak(double *ax, double *bx, doub
|
fparabu= *fa - A*(*ax-u)*(*ax-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); |
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); |
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); |
|
/* And thus,it can be that fu > *fc even if fparabu < *fc */ |
|
/* mnbrak (*ax=7.666299858533, *fa=299039.693133272231), (*bx=8.595447774979, *fb=298976.598289369489), |
|
(*cx=10.098840694817, *fc=298946.631474258087), (*u=9.852501168332, fu=298948.773013752128, fparabu=298945.434711494134) */ |
|
/* In that case, there is no bracket in the output! Routine is wrong with many consequences.*/ |
#endif |
#endif |
|
#ifdef MNBRAKORIGINAL |
|
#else |
|
if (fu > *fc) { |
|
#ifdef DEBUG |
|
printf("mnbrak4 fu > fc \n"); |
|
fprintf(ficlog, "mnbrak4 fu > fc\n"); |
|
#endif |
|
/* SHFT(u,*cx,*cx,u) /\* ie a=c, c=u and u=c; in that case, next SHFT(a,b,c,u) will give a=b=b, b=c=u, c=u=c and *\/ */ |
|
/* SHFT(*fa,*fc,fu,*fc) /\* (b, u, c) is a bracket while test fb > fc will be fu > fc will exit *\/ */ |
|
dum=u; /* Shifting c and u */ |
|
u = *cx; |
|
*cx = dum; |
|
dum = fu; |
|
fu = *fc; |
|
*fc =dum; |
|
} else { /* end */ |
|
#ifdef DEBUG |
|
printf("mnbrak3 fu < fc \n"); |
|
fprintf(ficlog, "mnbrak3 fu < fc\n"); |
|
#endif |
|
dum=u; /* Shifting c and u */ |
|
u = *cx; |
|
*cx = dum; |
|
dum = fu; |
|
fu = *fc; |
|
*fc =dum; |
|
} |
|
#endif |
} else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */ |
} else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */ |
|
#ifdef DEBUG |
|
printf("mnbrak2 u after c but before ulim\n"); |
|
fprintf(ficlog, "mnbrak2 u after c but before ulim\n"); |
|
#endif |
fu=(*func)(u); |
fu=(*func)(u); |
if (fu < *fc) { |
if (fu < *fc) { |
|
#ifdef DEBUG |
|
printf("mnbrak2 u after c but before ulim AND fu < fc\n"); |
|
fprintf(ficlog, "mnbrak2 u after c but before ulim AND fu <fc \n"); |
|
#endif |
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) { /* u outside ulim (verifying that ulim is beyond c) */ |
} else if ((u-ulim)*(ulim-*cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */ |
|
#ifdef DEBUG |
|
printf("mnbrak2 u outside ulim (verifying that ulim is beyond c)\n"); |
|
fprintf(ficlog, "mnbrak2 u outside ulim (verifying that ulim is beyond c)\n"); |
|
#endif |
u=ulim; |
u=ulim; |
fu=(*func)(u); |
fu=(*func)(u); |
} else { |
} else { /* u could be left to b (if r > q parabola has a maximum) */ |
|
#ifdef DEBUG |
|
printf("mnbrak2 u could be left to b (if r > q parabola has a maximum)\n"); |
|
fprintf(ficlog, "mnbrak2 u could be left to b (if r > q parabola has a maximum)\n"); |
|
#endif |
u=(*cx)+GOLD*(*cx-*bx); |
u=(*cx)+GOLD*(*cx-*bx); |
fu=(*func)(u); |
fu=(*func)(u); |
} |
} /* end tests */ |
SHFT(*ax,*bx,*cx,u) |
SHFT(*ax,*bx,*cx,u) |
SHFT(*fa,*fb,*fc,fu) |
SHFT(*fa,*fb,*fc,fu) |
} |
#ifdef DEBUG |
|
printf("mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu); |
|
fprintf(ficlog, "mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf), (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu); |
|
#endif |
|
} /* end while; ie return (a, b, c, fa, fb, fc) such that a < b < c with f(a) > f(b) and fb < f(c) */ |
} |
} |
|
|
/*************** linmin ************************/ |
/*************** linmin ************************/ |
Line 1334 void powell(double p[], double **xi, int
|
Line 1490 void powell(double p[], double **xi, int
|
double (*func)(double [])); |
double (*func)(double [])); |
int i,ibig,j; |
int i,ibig,j; |
double del,t,*pt,*ptt,*xit; |
double del,t,*pt,*ptt,*xit; |
|
double directest; |
double fp,fptt; |
double fp,fptt; |
double *xits; |
double *xits; |
int niterf, itmp; |
int niterf, itmp; |
Line 1393 void powell(double p[], double **xi, int
|
Line 1550 void powell(double p[], double **xi, int
|
#endif |
#endif |
printf("%d",i);fflush(stdout); |
printf("%d",i);fflush(stdout); |
fprintf(ficlog,"%d",i);fflush(ficlog); |
fprintf(ficlog,"%d",i);fflush(ficlog); |
linmin(p,xit,n,fret,func); |
linmin(p,xit,n,fret,func); /* xit[n] has been loaded for direction i */ |
if (fabs(fptt-(*fret)) > del) { |
if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions |
|
because that direction will be replaced unless the gain del is small |
|
in comparison with the 'probable' gain, mu^2, with the last average direction. |
|
Unless the n directions are conjugate some gain in the determinant may be obtained |
|
with the new direction. |
|
*/ |
del=fabs(fptt-(*fret)); |
del=fabs(fptt-(*fret)); |
ibig=i; |
ibig=i; |
} |
} |
Line 1414 void powell(double p[], double **xi, int
|
Line 1576 void powell(double p[], double **xi, int
|
fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
#endif |
#endif |
} /* end i */ |
} /* end i */ |
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { |
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */ |
#ifdef DEBUG |
#ifdef DEBUG |
int k[2],l; |
int k[2],l; |
k[0]=1; |
k[0]=1; |
Line 1446 void powell(double p[], double **xi, int
|
Line 1608 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++) { /* Computes an extrapolated point */ |
for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ |
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); /* f_3 */ |
if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */ |
if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */ |
/* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */ |
/* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */ |
/* From x1 (P0) distance of x2 is at h and x3 is 2h */ |
/* From x1 (P0) distance of x2 is at h and x3 is 2h */ |
/* Let f"(x2) be the 2nd derivative equal everywhere. */ |
/* Let f"(x2) be the 2nd derivative equal everywhere. */ |
/* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */ |
/* 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 */ |
/* 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) */ |
/* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del */ |
/* 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)-del*SQR(fp-fptt); */ |
|
#ifdef NRCORIGINAL |
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); |
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/ |
|
#else |
|
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */ |
t= t- del*SQR(fp-fptt); |
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); |
#endif |
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); |
directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */ |
#ifdef DEBUG |
#ifdef DEBUG |
|
printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); |
|
fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest); |
printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), |
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)); |
(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), |
fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt), |
Line 1476 void powell(double p[], double **xi, int
|
Line 1639 void powell(double p[], double **xi, int
|
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); |
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); |
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 |
#endif |
if (t < 0.0) { /* Then we use it for last direction */ |
#ifdef POWELLORIGINAL |
linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/ |
if (t < 0.0) { /* Then we use it for new direction */ |
|
#else |
|
if (directest*t < 0.0) { /* Contradiction between both tests */ |
|
printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del); |
|
printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); |
|
fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del); |
|
fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt); |
|
} |
|
if (directest < 0.0) { /* Then we use it for new direction */ |
|
#endif |
|
linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/ |
for (j=1;j<=n;j++) { |
for (j=1;j<=n;j++) { |
xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */ |
xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ |
xi[j][n]=xit[j]; /* and nth direction by the extrapolated */ |
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ |
} |
} |
printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig); |
printf("Gaining to use new 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); |
fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\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); |
Line 1506 double **prevalim(double **prlim, int nl
|
Line 1679 double **prevalim(double **prlim, int nl
|
{ |
{ |
/* Computes the prevalence limit in each live state at age x by left multiplying the unit |
/* Computes the prevalence limit in each live state at age x by left multiplying the unit |
matrix by transitions matrix until convergence is reached */ |
matrix by transitions matrix until convergence is reached */ |
|
|
int i, ii,j,k; |
int i, ii,j,k; |
double min, max, maxmin, maxmax,sumnew=0.; |
double min, max, maxmin, maxmax,sumnew=0.; |
/* double **matprod2(); */ /* test */ |
/* double **matprod2(); */ /* test */ |
double **out, cov[NCOVMAX+1], **pmij(); |
double **out, cov[NCOVMAX+1], **pmij(); |
double **newm; |
double **newm; |
double agefin, delaymax=50 ; /* Max number of years to converge */ |
double agefin, delaymax=50 ; /* Max number of years to converge */ |
|
|
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (ii=1;ii<=nlstate+ndeath;ii++) |
for (j=1;j<=nlstate+ndeath;j++){ |
for (j=1;j<=nlstate+ndeath;j++){ |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
} |
} |
|
|
cov[1]=1.; |
cov[1]=1.; |
|
|
/* Even if hstepm = 1, at least one multiplication by the unit matrix */ |
/* Even if hstepm = 1, at least one multiplication by the unit matrix */ |
for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){ |
for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){ |
newm=savm; |
newm=savm; |
/* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
Line 1558 double **prevalim(double **prlim, int nl
|
Line 1731 double **prevalim(double **prlim, int nl
|
} |
} |
maxmin=max-min; |
maxmin=max-min; |
maxmax=FMAX(maxmax,maxmin); |
maxmax=FMAX(maxmax,maxmin); |
} |
} /* j loop */ |
if(maxmax < ftolpl){ |
if(maxmax < ftolpl){ |
return prlim; |
return prlim; |
} |
} |
} |
} /* age loop */ |
|
return prlim; /* should not reach here */ |
} |
} |
|
|
/*************** transition probabilities ***************/ |
/*************** transition probabilities ***************/ |
Line 1845 double func( double *x)
|
Line 2019 double func( double *x)
|
which slows down the processing. The difference can be up to 10% |
which slows down the processing. The difference can be up to 10% |
lower mortality. |
lower mortality. |
*/ |
*/ |
lli=log(out[s1][s2] - savm[s1][s2]); |
/* If, at the beginning of the maximization mostly, the |
|
cumulative probability or probability to be dead is |
|
constant (ie = 1) over time d, the difference is equal to |
|
0. out[s1][3] = savm[s1][3]: probability, being at state |
|
s1 at precedent wave, to be dead a month before current |
|
wave is equal to probability, being at state s1 at |
|
precedent wave, to be dead at mont of the current |
|
wave. Then the observed probability (that this person died) |
|
is null according to current estimated parameter. In fact, |
|
it should be very low but not zero otherwise the log go to |
|
infinity. |
|
*/ |
|
/* #ifdef INFINITYORIGINAL */ |
|
/* lli=log(out[s1][s2] - savm[s1][s2]); */ |
|
/* #else */ |
|
/* if ((out[s1][s2] - savm[s1][s2]) < mytinydouble) */ |
|
/* lli=log(mytinydouble); */ |
|
/* else */ |
|
/* lli=log(out[s1][s2] - savm[s1][s2]); */ |
|
/* #endif */ |
|
lli=log(out[s1][s2] - savm[s1][s2]); |
|
|
} else if (s2==-2) { |
} else if (s2==-2) { |
for (j=1,survp=0. ; j<=nlstate; j++) |
for (j=1,survp=0. ; j<=nlstate; j++) |
Line 1877 double func( double *x)
|
Line 2070 double func( double *x)
|
ipmx +=1; |
ipmx +=1; |
sw += weight[i]; |
sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
|
/* if (lli < log(mytinydouble)){ */ |
|
/* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */ |
|
/* fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */ |
|
/* } */ |
} /* end of wave */ |
} /* end of wave */ |
} /* end of individual */ |
} /* end of individual */ |
} else if(mle==2){ |
} else if(mle==2){ |
Line 2147 void likelione(FILE *ficres,double p[],
|
Line 2344 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 */ |
Line 2212 void mlikeli(FILE *ficres,double p[], in
|
Line 2409 void mlikeli(FILE *ficres,double p[], in
|
#endif |
#endif |
free_matrix(xi,1,npar,1,npar); |
free_matrix(xi,1,npar,1,npar); |
fclose(ficrespow); |
fclose(ficrespow); |
printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
printf("#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
fprintf(ficlog,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p)); |
|
|
} |
} |
|
|
Line 2511 void freqsummary(char fileres[], int ia
|
Line 2708 void freqsummary(char fileres[], int ia
|
|
|
first=1; |
first=1; |
|
|
/* for(k1=1; k1<=j ; k1++){ /* Loop on covariates */ |
/* for(k1=1; k1<=j ; k1++){ */ /* Loop on covariates */ |
/* for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */ |
/* for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */ |
/* j1++; |
/* j1++; */ |
*/ |
|
for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ |
for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]); |
scanf("%d", i);*/ |
scanf("%d", i);*/ |
Line 2891 void tricode(int *Tvar, int **nbcode, in
|
Line 3087 void tricode(int *Tvar, int **nbcode, in
|
{ |
{ |
/**< Uses cptcovn+2*cptcovprod as the number of covariates */ |
/**< Uses cptcovn+2*cptcovprod as the number of covariates */ |
/* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 |
/* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 |
/* Boring subroutine which should only output nbcode[Tvar[j]][k] |
* Boring subroutine which should only output nbcode[Tvar[j]][k] |
* Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2) |
* Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2) |
/* nbcode[Tvar[j]][1]= |
* nbcode[Tvar[j]][1]= |
*/ |
*/ |
|
|
int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX; |
int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX; |
Line 3321 void varevsij(char optionfilefiname[], d
|
Line 3517 void varevsij(char optionfilefiname[], d
|
/* Variance of health expectancies */ |
/* Variance of health expectancies */ |
/* 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 **newm;*/ |
|
/* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/ |
|
|
|
int movingaverage(); |
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 ; |
Line 3598 void varevsij(char optionfilefiname[], d
|
Line 3797 void varevsij(char optionfilefiname[], d
|
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */ |
fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l lt 2 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev)); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev)); |
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); |
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev)); |
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); |
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); |
Line 4019 To be simple, these graphs help to under
|
Line 4218 To be simple, these graphs help to under
|
} /* k12 */ |
} /* k12 */ |
} /*l1 */ |
} /*l1 */ |
}/* k1 */ |
}/* k1 */ |
/* } /* loop covariates */ |
/* } */ /* loop covariates */ |
} |
} |
free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); |
free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); |
free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); |
free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); |
Line 4087 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 4286 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 cross-sectional prevalence in each state (1 to %d) to period (stable) prevalence in specific 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 4196 void printinggnuplot(char fileres[], cha
|
Line 4395 void printinggnuplot(char fileres[], cha
|
fprintf(ficgp,"set xlabel \"Age\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \n\ |
set ylabel \"Probability\" \n\ |
set ylabel \"Probability\" \n\ |
set ter png small size 320, 240\n\ |
set ter png small size 320, 240\n\ |
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1); |
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1); |
|
|
for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
} |
} |
fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); |
fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); |
for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
} |
} |
fprintf(ficgp,"\" t\"95\%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); |
fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); |
for (i=1; i<= nlstate ; i ++) { |
for (i=1; i<= nlstate ; i ++) { |
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)"); |
if (i==cpt) fprintf(ficgp," %%lf (%%lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
} |
} |
fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1)); |
fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1)); |
} |
} |
Line 4223 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 4422 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
|
|
for (i=1; i<= nlstate+1 ; i ++) { |
for (i=1; i<= nlstate+1 ; i ++) { |
k=2*i; |
k=2*i; |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1); |
for (j=1; j<= nlstate+1 ; j ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
} |
} |
if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,"); |
if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,"); |
else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1); |
else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1); |
for (j=1; j<= nlstate+1 ; j ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
} |
} |
fprintf(ficgp,"\" t\"\" w l lt 0,"); |
fprintf(ficgp,"\" t\"\" w l lt 0,"); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1); |
fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1); |
for (j=1; j<= nlstate+1 ; j ++) { |
for (j=1; j<= nlstate+1 ; j ++) { |
if (j==i) fprintf(ficgp," \%%lf (\%%lf)"); |
if (j==i) fprintf(ficgp," %%lf (%%lf)"); |
else fprintf(ficgp," \%%*lf (\%%*lf)"); |
else fprintf(ficgp," %%*lf (%%*lf)"); |
} |
} |
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0"); |
if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0"); |
else fprintf(ficgp,"\" t\"\" w l lt 0,"); |
else fprintf(ficgp,"\" t\"\" w l lt 0,"); |
Line 4407 int movingaverage(double ***probs, doubl
|
Line 4606 int movingaverage(double ***probs, doubl
|
|
|
|
|
/************** Forecasting ******************/ |
/************** Forecasting ******************/ |
prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){ |
void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){ |
/* proj1, year, month, day of starting projection |
/* proj1, year, month, day of starting projection |
agemin, agemax range of age |
agemin, agemax range of age |
dateprev1 dateprev2 range of dates during which prevalence is computed |
dateprev1 dateprev2 range of dates during which prevalence is computed |
Line 4530 prevforecast(char fileres[], double anpr
|
Line 4729 prevforecast(char fileres[], double anpr
|
} |
} |
|
|
/************** Forecasting *****not tested NB*************/ |
/************** Forecasting *****not tested NB*************/ |
populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ |
void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ |
|
|
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; |
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; |
int *popage; |
int *popage; |
Line 5006 int readdata(char datafile[], int firsto
|
Line 5205 int readdata(char datafile[], int firsto
|
|
|
strcpy(line,stra); |
strcpy(line,stra); |
cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ |
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
} |
} |
else if(iout=sscanf(strb,"%s.",dummy) != 0){ |
else if( (iout=sscanf(strb,"%s.",dummy)) != 0){ |
month=99; |
month=99; |
year=9999; |
year=9999; |
}else{ |
}else{ |
Line 5022 int readdata(char datafile[], int firsto
|
Line 5221 int readdata(char datafile[], int firsto
|
} /* ENd Waves */ |
} /* ENd Waves */ |
|
|
cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ |
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
} |
} |
else if(iout=sscanf(strb,"%s.",dummy) != 0){ |
else if( (iout=sscanf(strb,"%s.",dummy)) != 0){ |
month=99; |
month=99; |
year=9999; |
year=9999; |
}else{ |
}else{ |
Line 5037 int readdata(char datafile[], int firsto
|
Line 5236 int readdata(char datafile[], int firsto
|
strcpy(line,stra); |
strcpy(line,stra); |
|
|
cutv(stra, strb,line,' '); |
cutv(stra, strb,line,' '); |
if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){ |
if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){ |
} |
} |
else if(iout=sscanf(strb,"%s.", dummy) != 0){ |
else if( (iout=sscanf(strb,"%s.", dummy)) != 0){ |
month=99; |
month=99; |
year=9999; |
year=9999; |
}else{ |
}else{ |
Line 5137 void removespace(char *str) {
|
Line 5336 void removespace(char *str) {
|
do |
do |
while (*p2 == ' ') |
while (*p2 == ' ') |
p2++; |
p2++; |
while (*p1++ = *p2++); |
while (*p1++ == *p2++); |
} |
} |
|
|
int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns: |
int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns: |
Line 5167 int decodemodel ( char model[], int last
|
Line 5366 int decodemodel ( char model[], int last
|
cptcovs=j+1-j1; /**< Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2 */ |
cptcovs=j+1-j1; /**< Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2 */ |
cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/ |
cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/ |
/* including age products which are counted in cptcovage. |
/* including age products which are counted in cptcovage. |
/* but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */ |
* but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */ |
cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */ |
cptcovprod=j1; /**< Number of products V1*V2 +v3*age = 2 */ |
cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */ |
cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */ |
strcpy(modelsav,model); |
strcpy(modelsav,model); |
Line 5309 int decodemodel ( char model[], int last
|
Line 5508 int decodemodel ( char model[], int last
|
return (1); |
return (1); |
} |
} |
|
|
calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn ) |
int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn ) |
{ |
{ |
int i, m; |
int i, m; |
|
|
Line 5320 calandcheckages(int imx, int maxwav, dou
|
Line 5519 calandcheckages(int imx, int maxwav, dou
|
s[m][i]=-1; |
s[m][i]=-1; |
} |
} |
if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){ |
if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){ |
*nberr++; |
*nberr = *nberr + 1; |
printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i); |
printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr); |
fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i); |
fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr); |
s[m][i]=-1; |
s[m][i]=-1; |
} |
} |
if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){ |
if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){ |
*nberr++; |
(*nberr)++; |
printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); |
printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); |
fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); |
fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); |
s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */ |
s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */ |
Line 5339 calandcheckages(int imx, int maxwav, dou
|
Line 5538 calandcheckages(int imx, int maxwav, dou
|
for(m=firstpass; (m<= lastpass); m++){ |
for(m=firstpass; (m<= lastpass); m++){ |
if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ |
if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){ |
if (s[m][i] >= nlstate+1) { |
if (s[m][i] >= nlstate+1) { |
if(agedc[i]>0) |
if(agedc[i]>0){ |
if((int)moisdc[i]!=99 && (int)andc[i]!=9999) |
if((int)moisdc[i]!=99 && (int)andc[i]!=9999){ |
agev[m][i]=agedc[i]; |
agev[m][i]=agedc[i]; |
/*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ |
/*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ |
else { |
}else { |
if ((int)andc[i]!=9999){ |
if ((int)andc[i]!=9999){ |
nbwarn++; |
nbwarn++; |
printf("Warning negative age at death: %ld line:%d\n",num[i],i); |
printf("Warning negative age at death: %ld line:%d\n",num[i],i); |
Line 5351 calandcheckages(int imx, int maxwav, dou
|
Line 5550 calandcheckages(int imx, int maxwav, dou
|
agev[m][i]=-1; |
agev[m][i]=-1; |
} |
} |
} |
} |
|
} /* agedc > 0 */ |
} |
} |
else if(s[m][i] !=9){ /* Standard case, age in fractional |
else if(s[m][i] !=9){ /* Standard case, age in fractional |
years but with the precision of a month */ |
years but with the precision of a month */ |
Line 5381 calandcheckages(int imx, int maxwav, dou
|
Line 5581 calandcheckages(int imx, int maxwav, dou
|
for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { |
for(m=firstpass; (m<=lastpass); m++){ |
for(m=firstpass; (m<=lastpass); m++){ |
if (s[m][i] > (nlstate+ndeath)) { |
if (s[m][i] > (nlstate+ndeath)) { |
*nberr++; |
(*nberr)++; |
printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); |
printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); |
fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); |
fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); |
return 1; |
return 1; |
Line 5406 calandcheckages(int imx, int maxwav, dou
|
Line 5606 calandcheckages(int imx, int maxwav, dou
|
return (1); |
return (1); |
} |
} |
|
|
|
#if defined(_MSC_VER) |
|
/*printf("Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/ |
|
/*fprintf(ficlog, "Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/ |
|
//#include "stdafx.h" |
|
//#include <stdio.h> |
|
//#include <tchar.h> |
|
//#include <windows.h> |
|
//#include <iostream> |
|
typedef BOOL(WINAPI *LPFN_ISWOW64PROCESS) (HANDLE, PBOOL); |
|
|
|
LPFN_ISWOW64PROCESS fnIsWow64Process; |
|
|
|
BOOL IsWow64() |
|
{ |
|
BOOL bIsWow64 = FALSE; |
|
|
|
//typedef BOOL (APIENTRY *LPFN_ISWOW64PROCESS) |
|
// (HANDLE, PBOOL); |
|
|
|
//LPFN_ISWOW64PROCESS fnIsWow64Process; |
|
|
|
HMODULE module = GetModuleHandle(_T("kernel32")); |
|
const char funcName[] = "IsWow64Process"; |
|
fnIsWow64Process = (LPFN_ISWOW64PROCESS) |
|
GetProcAddress(module, funcName); |
|
|
|
if (NULL != fnIsWow64Process) |
|
{ |
|
if (!fnIsWow64Process(GetCurrentProcess(), |
|
&bIsWow64)) |
|
//throw std::exception("Unknown error"); |
|
printf("Unknown error\n"); |
|
} |
|
return bIsWow64 != FALSE; |
|
} |
|
#endif |
|
|
|
void syscompilerinfo() |
|
{ |
|
/* #include "syscompilerinfo.h"*/ |
|
/* command line Intel compiler 64bit windows: |
|
/GS /W3 /Gy /Zc:wchar_t /Zi /O2 /Fd"x64\Release\vc120.pdb" /D "WIN32" |
|
/D "NDEBUG" /D "_CONSOLE" /D "_LIB" /D "_UNICODE" /D "UNICODE" /Qipo |
|
/Zc:forScope /Oi /MD /Fa"x64\Release\" /EHsc /nologo |
|
/Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\IMaCh.pch" */ |
|
/* |
|
/GS /W3 /Gy /Zc:wchar_t /Zi /O3 /Fd"x64\Release\vc120.pdb" /D "WIN32" |
|
/D "NDEBUG" /D "_CONSOLE" /D "_LIB" /D "_UNICODE" /D "UNICODE" /Qipo |
|
/Zc:forScope /Oi /MD /Fa"x64\Release\" /EHsc /nologo /Qparallel |
|
/Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\IMaCh.pch" */ |
|
#if defined __INTEL_COMPILER |
|
#if defined(__GNUC__) |
|
struct utsname sysInfo; /* For Intel on Linux and OS/X */ |
|
#endif |
|
#elif defined(__GNUC__) |
|
#ifndef __APPLE__ |
|
#include <gnu/libc-version.h> /* Only on gnu */ |
|
#endif |
|
struct utsname sysInfo; |
|
int cross = CROSS; |
|
if (cross){ |
|
printf("Cross-"); |
|
fprintf(ficlog, "Cross-"); |
|
} |
|
#endif |
|
|
|
#include <stdint.h> |
|
|
|
printf("Compiled with:");fprintf(ficlog,"Compiled with:"); |
|
#if defined(__clang__) |
|
printf(" Clang/LLVM");fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */ |
|
#endif |
|
#if defined(__ICC) || defined(__INTEL_COMPILER) |
|
printf(" Intel ICC/ICPC");fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */ |
|
#endif |
|
#if defined(__GNUC__) || defined(__GNUG__) |
|
printf(" GNU GCC/G++");fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */ |
|
#endif |
|
#if defined(__HP_cc) || defined(__HP_aCC) |
|
printf(" Hewlett-Packard C/aC++");fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */ |
|
#endif |
|
#if defined(__IBMC__) || defined(__IBMCPP__) |
|
printf(" IBM XL C/C++"); fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */ |
|
#endif |
|
#if defined(_MSC_VER) |
|
printf(" Microsoft Visual Studio");fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */ |
|
#endif |
|
#if defined(__PGI) |
|
printf(" Portland Group PGCC/PGCPP");fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */ |
|
#endif |
|
#if defined(__SUNPRO_C) || defined(__SUNPRO_CC) |
|
printf(" Oracle Solaris Studio");fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */ |
|
#endif |
|
printf(" for ");fprintf(ficlog," for "); |
|
|
|
// http://stackoverflow.com/questions/4605842/how-to-identify-platform-compiler-from-preprocessor-macros |
|
#ifdef _WIN32 // note the underscore: without it, it's not msdn official! |
|
// Windows (x64 and x86) |
|
printf("Windows (x64 and x86) ");fprintf(ficlog,"Windows (x64 and x86) "); |
|
#elif __unix__ // all unices, not all compilers |
|
// Unix |
|
printf("Unix ");fprintf(ficlog,"Unix "); |
|
#elif __linux__ |
|
// linux |
|
printf("linux ");fprintf(ficlog,"linux "); |
|
#elif __APPLE__ |
|
// Mac OS, not sure if this is covered by __posix__ and/or __unix__ though.. |
|
printf("Mac OS ");fprintf(ficlog,"Mac OS "); |
|
#endif |
|
|
|
/* __MINGW32__ */ |
|
/* __CYGWIN__ */ |
|
/* __MINGW64__ */ |
|
// http://msdn.microsoft.com/en-us/library/b0084kay.aspx |
|
/* _MSC_VER //the Visual C++ compiler is 17.00.51106.1, the _MSC_VER macro evaluates to 1700. Type cl /? */ |
|
/* _MSC_FULL_VER //the Visual C++ compiler is 15.00.20706.01, the _MSC_FULL_VER macro evaluates to 150020706 */ |
|
/* _WIN64 // Defined for applications for Win64. */ |
|
/* _M_X64 // Defined for compilations that target x64 processors. */ |
|
/* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */ |
|
|
|
#if UINTPTR_MAX == 0xffffffff |
|
printf(" 32-bit"); fprintf(ficlog," 32-bit");/* 32-bit */ |
|
#elif UINTPTR_MAX == 0xffffffffffffffff |
|
printf(" 64-bit"); fprintf(ficlog," 64-bit");/* 64-bit */ |
|
#else |
|
printf(" wtf-bit"); fprintf(ficlog," wtf-bit");/* wtf */ |
|
#endif |
|
|
|
#if defined(__GNUC__) |
|
# if defined(__GNUC_PATCHLEVEL__) |
|
# define __GNUC_VERSION__ (__GNUC__ * 10000 \ |
|
+ __GNUC_MINOR__ * 100 \ |
|
+ __GNUC_PATCHLEVEL__) |
|
# else |
|
# define __GNUC_VERSION__ (__GNUC__ * 10000 \ |
|
+ __GNUC_MINOR__ * 100) |
|
# endif |
|
printf(" using GNU C version %d.\n", __GNUC_VERSION__); |
|
fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__); |
|
|
|
if (uname(&sysInfo) != -1) { |
|
printf("Running on: %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine); |
|
fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine); |
|
} |
|
else |
|
perror("uname() error"); |
|
//#ifndef __INTEL_COMPILER |
|
#if !defined (__INTEL_COMPILER) && !defined(__APPLE__) |
|
printf("GNU libc version: %s\n", gnu_get_libc_version()); |
|
fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version()); |
|
#endif |
|
#endif |
|
|
|
// void main() |
|
// { |
|
#if defined(_MSC_VER) |
|
if (IsWow64()){ |
|
printf("The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n"); |
|
fprintf(ficlog, "The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n"); |
|
} |
|
else{ |
|
printf("The process is not running under WOW64 (i.e probably on a 64bit Windows).\n"); |
|
fprintf(ficlog,"The programm is not running under WOW64 (i.e probably on a 64bit Windows).\n"); |
|
} |
|
// printf("\nPress Enter to continue..."); |
|
// getchar(); |
|
// } |
|
|
|
#endif |
|
|
|
|
|
} |
|
|
|
int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){ |
|
/*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
|
int i, j, k, i1 ; |
|
double ftolpl = 1.e-10; |
|
double age, agebase, agelim; |
|
|
|
strcpy(filerespl,"pl"); |
|
strcat(filerespl,fileres); |
|
if((ficrespl=fopen(filerespl,"w"))==NULL) { |
|
printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1; |
|
fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1; |
|
} |
|
printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl); |
|
fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl); |
|
pstamp(ficrespl); |
|
fprintf(ficrespl,"# Period (stable) prevalence \n"); |
|
fprintf(ficrespl,"#Age "); |
|
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
|
fprintf(ficrespl,"\n"); |
|
|
|
/* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */ |
|
|
|
agebase=ageminpar; |
|
agelim=agemaxpar; |
|
|
|
i1=pow(2,cptcoveff); |
|
if (cptcovn < 1){i1=1;} |
|
|
|
for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
|
/* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */ |
|
//for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
|
k=k+1; |
|
/* to clean */ |
|
//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]); |
|
fprintf(ficrespl,"\n#******"); |
|
printf("\n#******"); |
|
fprintf(ficlog,"\n#******"); |
|
for(j=1;j<=cptcoveff;j++) { |
|
fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
|
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
|
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
|
} |
|
fprintf(ficrespl,"******\n"); |
|
printf("******\n"); |
|
fprintf(ficlog,"******\n"); |
|
|
|
fprintf(ficrespl,"#Age "); |
|
for(j=1;j<=cptcoveff;j++) { |
|
fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
|
} |
|
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
|
fprintf(ficrespl,"\n"); |
|
|
|
for (age=agebase; age<=agelim; age++){ |
|
/* for (age=agebase; age<=agebase; age++){ */ |
|
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); |
|
fprintf(ficrespl,"%.0f ",age ); |
|
for(j=1;j<=cptcoveff;j++) |
|
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
|
for(i=1; i<=nlstate;i++) |
|
fprintf(ficrespl," %.5f", prlim[i][i]); |
|
fprintf(ficrespl,"\n"); |
|
} /* Age */ |
|
/* was end of cptcod */ |
|
} /* cptcov */ |
|
return 0; |
|
} |
|
|
|
int hPijx(double *p, int bage, int fage){ |
|
/*------------- h Pij x at various ages ------------*/ |
|
|
|
int stepsize; |
|
int agelim; |
|
int hstepm; |
|
int nhstepm; |
|
int h, i, i1, j, k; |
|
|
|
double agedeb; |
|
double ***p3mat; |
|
|
|
strcpy(filerespij,"pij"); strcat(filerespij,fileres); |
|
if((ficrespij=fopen(filerespij,"w"))==NULL) { |
|
printf("Problem with Pij resultfile: %s\n", filerespij); return 1; |
|
fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij); return 1; |
|
} |
|
printf("Computing pij: result on file '%s' \n", filerespij); |
|
fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij); |
|
|
|
stepsize=(int) (stepm+YEARM-1)/YEARM; |
|
/*if (stepm<=24) stepsize=2;*/ |
|
|
|
agelim=AGESUP; |
|
hstepm=stepsize*YEARM; /* Every year of age */ |
|
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ |
|
|
|
/* hstepm=1; aff par mois*/ |
|
pstamp(ficrespij); |
|
fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x "); |
|
i1= pow(2,cptcoveff); |
|
/* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */ |
|
/* /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */ |
|
/* k=k+1; */ |
|
for (k=1; k <= (int) pow(2,cptcoveff); k++){ |
|
fprintf(ficrespij,"\n#****** "); |
|
for(j=1;j<=cptcoveff;j++) |
|
fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
|
fprintf(ficrespij,"******\n"); |
|
|
|
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ |
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ |
|
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ |
|
|
|
/* nhstepm=nhstepm*YEARM; aff par mois*/ |
|
|
|
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
|
oldm=oldms;savm=savms; |
|
hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); |
|
fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j="); |
|
for(i=1; i<=nlstate;i++) |
|
for(j=1; j<=nlstate+ndeath;j++) |
|
fprintf(ficrespij," %1d-%1d",i,j); |
|
fprintf(ficrespij,"\n"); |
|
for (h=0; h<=nhstepm; h++){ |
|
/*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/ |
|
fprintf(ficrespij,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm ); |
|
for(i=1; i<=nlstate;i++) |
|
for(j=1; j<=nlstate+ndeath;j++) |
|
fprintf(ficrespij," %.5f", p3mat[i][j][h]); |
|
fprintf(ficrespij,"\n"); |
|
} |
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
|
fprintf(ficrespij,"\n"); |
|
} |
|
/*}*/ |
|
} |
|
return 0; |
|
} |
|
|
|
|
/***********************************************/ |
/***********************************************/ |
/**************** Main Program *****************/ |
/**************** Main Program *****************/ |
Line 5436 int main(int argc, char *argv[])
|
Line 5947 int main(int argc, char *argv[])
|
double agedeb; |
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 dum; /* Dummy variable */ |
double dum; /* Dummy variable */ |
double ***p3mat; |
double ***p3mat; |
double ***mobaverage; |
double ***mobaverage; |
Line 5514 int main(int argc, char *argv[])
|
Line 6025 int main(int argc, char *argv[])
|
|
|
nberr=0; /* Number of errors and warnings */ |
nberr=0; /* Number of errors and warnings */ |
nbwarn=0; |
nbwarn=0; |
|
#ifdef WIN32 |
|
_getcwd(pathcd, size); |
|
#else |
getcwd(pathcd, size); |
getcwd(pathcd, size); |
|
#endif |
|
|
printf("\n%s\n%s",version,fullversion); |
printf("\n%s\n%s",version,fullversion); |
if(argc <=1){ |
if(argc <=1){ |
Line 5550 int main(int argc, char *argv[])
|
Line 6065 int main(int argc, char *argv[])
|
/* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */ |
/* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */ |
split(pathtot,path,optionfile,optionfilext,optionfilefiname); |
split(pathtot,path,optionfile,optionfilext,optionfilefiname); |
printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); |
printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); |
|
#ifdef WIN32 |
|
_chdir(path); /* Can be a relative path */ |
|
if(_getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */ |
|
#else |
chdir(path); /* Can be a relative path */ |
chdir(path); /* Can be a relative path */ |
if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */ |
if (getcwd(pathcd, MAXLINE) > 0) /* So pathcd is the full path */ |
printf("Current directory %s!\n",pathcd); |
#endif |
|
printf("Current directory %s!\n",pathcd); |
strcpy(command,"mkdir "); |
strcpy(command,"mkdir "); |
strcat(command,optionfilefiname); |
strcat(command,optionfilefiname); |
if((outcmd=system(command)) != 0){ |
if((outcmd=system(command)) != 0){ |
printf("Problem creating directory or it already exists %s%s, err=%d\n",path,optionfilefiname,outcmd); |
printf("Directory already exists (or can't create it) %s%s, err=%d\n",path,optionfilefiname,outcmd); |
/* fprintf(ficlog,"Problem creating directory %s%s\n",path,optionfilefiname); */ |
/* fprintf(ficlog,"Problem creating directory %s%s\n",path,optionfilefiname); */ |
/* fclose(ficlog); */ |
/* fclose(ficlog); */ |
/* exit(1); */ |
/* exit(1); */ |
Line 5583 int main(int argc, char *argv[])
|
Line 6103 int main(int argc, char *argv[])
|
optionfilext=%s\n\ |
optionfilext=%s\n\ |
optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); |
optionfilefiname='%s'\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname); |
|
|
|
syscompilerinfo(); |
|
|
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); |
Line 6050 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 6572 Title=%s <br>Datafile=%s Firstpass=%d La
|
|
|
strcpy(pathr,path); |
strcpy(pathr,path); |
strcat(pathr,optionfilefiname); |
strcat(pathr,optionfilefiname); |
|
#ifdef WIN32 |
|
_chdir(optionfilefiname); /* Move to directory named optionfile */ |
|
#else |
chdir(optionfilefiname); /* Move to directory named optionfile */ |
chdir(optionfilefiname); /* Move to directory named optionfile */ |
|
#endif |
|
|
|
|
/* Calculates basic frequencies. Computes observed prevalence at single age |
/* Calculates basic frequencies. Computes observed prevalence at single age |
and prints on file fileres'p'. */ |
and prints on file fileres'p'. */ |
Line 6534 Interval (in months) between two waves:
|
Line 7061 Interval (in months) between two waves:
|
|
|
|
|
/*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
/*--------------- Prevalence limit (period or stable prevalence) --------------*/ |
#include "prevlim.h" /* Use ficrespl, ficlog */ |
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */ |
|
prlim=matrix(1,nlstate,1,nlstate); |
|
prevalence_limit(p, prlim, ageminpar, agemaxpar); |
fclose(ficrespl); |
fclose(ficrespl); |
|
|
#ifdef FREEEXIT2 |
#ifdef FREEEXIT2 |
Line 6542 Interval (in months) between two waves:
|
Line 7071 Interval (in months) between two waves:
|
#endif |
#endif |
|
|
/*------------- h Pij x at various ages ------------*/ |
/*------------- h Pij x at various ages ------------*/ |
#include "hpijx.h" |
/*#include "hpijx.h"*/ |
|
hPijx(p, bage, fage); |
fclose(ficrespij); |
fclose(ficrespij); |
|
|
/*-------------- Variance of one-step probabilities---*/ |
/*-------------- Variance of one-step probabilities---*/ |
Line 6844 Interval (in months) between two waves:
|
Line 7374 Interval (in months) between two waves:
|
|
|
|
|
printf("Before Current directory %s!\n",pathcd); |
printf("Before Current directory %s!\n",pathcd); |
|
#ifdef WIN32 |
|
if (_chdir(pathcd) != 0) |
|
printf("Can't move to directory %s!\n",path); |
|
if(_getcwd(pathcd,MAXLINE) > 0) |
|
#else |
if(chdir(pathcd) != 0) |
if(chdir(pathcd) != 0) |
printf("Can't move to directory %s!\n",path); |
printf("Can't move to directory %s!\n", path); |
if(getcwd(pathcd,MAXLINE) > 0) |
if (getcwd(pathcd, MAXLINE) > 0) |
|
#endif |
printf("Current directory %s!\n",pathcd); |
printf("Current directory %s!\n",pathcd); |
/*strcat(plotcmd,CHARSEPARATOR);*/ |
/*strcat(plotcmd,CHARSEPARATOR);*/ |
sprintf(plotcmd,"gnuplot"); |
sprintf(plotcmd,"gnuplot"); |