|
version 1.162, 2014/09/25 11:43:39
|
version 1.183, 2015/03/10 20:34:32
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
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 |
| |
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 |
Revision 1.162 2014/09/25 11:43:39 brouard |
| Summary: temporary backup 0.99! |
Summary: temporary backup 0.99! |
| |
|
|
Line 489
|
Line 572
|
| end |
end |
| */ |
*/ |
| |
|
| |
#define POWELL /* Instead of NLOPT */ |
| |
/* #define POWELLORIGINAL */ /* Don't use Directest to decide new direction but original Powell test */ |
| |
|
| |
|
| |
|
| #include <math.h> |
#include <math.h> |
| #include <stdio.h> |
#include <stdio.h> |
| #include <stdlib.h> |
#include <stdlib.h> |
|
Line 499
|
Line 582
|
| |
|
| #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 523
|
Line 613
|
| #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 554 typedef struct {
|
Line 645 typedef struct {
|
| #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 _WIN32 |
#ifdef _WIN32 |
| #define DIRSEPARATOR '\\' |
#define DIRSEPARATOR '\\' |
| #define CHARSEPARATOR "\\" |
#define CHARSEPARATOR "\\" |
|
Line 568 typedef struct {
|
Line 659 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 683 static double maxarg1,maxarg2;
|
Line 774 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 1146 double brent(double ax, double bx, doubl
|
Line 1242 double brent(double ax, double bx, doubl
|
| int iter; |
int iter; |
| double a,b,d,etemp; |
double a,b,d,etemp; |
| double fu=0,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 1195 double brent(double ax, double bx, doubl
|
Line 1291 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 1218 double brent(double ax, double bx, doubl
|
Line 1314 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 1226 void mnbrak(double *ax, double *bx, doub
|
Line 1326 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 |
| |
/* 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); |
| |
/* 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 |
| |
#ifdef MNBRAKORI |
| |
#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 1318 void powell(double p[], double **xi, int
|
Line 1482 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 1372 void powell(double p[], double **xi, int
|
Line 1537 void powell(double p[], double **xi, int
|
| 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); |
| 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 1398 void powell(double p[], double **xi, int
|
Line 1568 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 1430 void powell(double p[], double **xi, int
|
Line 1600 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 1460 void powell(double p[], double **xi, int
|
Line 1631 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); |
| |
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); |
| |
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); |
| |
fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig); |
| for(j=1;j<=n;j++){ |
for(j=1;j<=n;j++){ |
| printf(" %.12e",xit[j]); |
printf(" %.12e",xit[j]); |
| fprintf(ficlog," %.12e",xit[j]); |
fprintf(ficlog," %.12e",xit[j]); |
|
Line 1488 double **prevalim(double **prlim, int nl
|
Line 1671 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 1540 double **prevalim(double **prlim, int nl
|
Line 1723 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 1566 double **pmij(double **ps, double *cov,
|
Line 1750 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 1827 double func( double *x)
|
Line 2011 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 1859 double func( double *x)
|
Line 2062 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 2129 void likelione(FILE *ficres,double p[],
|
Line 2336 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 2194 void mlikeli(FILE *ficres,double p[], in
|
Line 2401 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 2205 void hesscov(double **matcov, double p[]
|
Line 2412 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 2339 double hessii(double x[], double delta,
|
Line 2546 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 2354 double hessii(double x[], double delta,
|
Line 2561 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 2469 void pstamp(FILE *fichier)
|
Line 2676 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 2493 void freqsummary(char fileres[], int ia
|
Line 2700 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 2650 void prevalence(double ***probs, double
|
Line 2856 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 2744 void concatwav(int wav[], int **dh, int
|
Line 2950 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 2873 void tricode(int *Tvar, int **nbcode, in
|
Line 3079 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 2984 void evsij(double ***eij, double x[], in
|
Line 3190 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 3303 void varevsij(char optionfilefiname[], d
|
Line 3509 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 ; |
| 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 3580 void varevsij(char optionfilefiname[], d
|
Line 3789 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 3606 void varprevlim(char fileres[], double *
|
Line 3815 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 3688 void varprevlim(char fileres[], double *
|
Line 3896 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 3698 void varprob(char optionfilefiname[], do
|
Line 3906 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 4002 To be simple, these graphs help to under
|
Line 4210 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 4070 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 4278 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 4156 true period expectancies (those weighted
|
Line 4364 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 4179 void printinggnuplot(char fileres[], cha
|
Line 4387 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 4206 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 4414 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 4342 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 4550 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 4390 int movingaverage(double ***probs, doubl
|
Line 4598 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 |
| 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 4514 prevforecast(char fileres[], double anpr
|
Line 4721 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 4691 void prwizard(int ncovmodel, int nlstate
|
Line 4898 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 4913 fprintf(fichtm,"<ul><li><h4>Life table</
|
Line 5120 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 4938 int readdata(char datafile[], int firsto
|
Line 5145 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 4969 int readdata(char datafile[], int firsto
|
Line 5176 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 4992 int readdata(char datafile[], int firsto
|
Line 5197 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 5008 int readdata(char datafile[], int firsto
|
Line 5213 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 5023 int readdata(char datafile[], int firsto
|
Line 5228 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 5110 int readdata(char datafile[], int firsto
|
Line 5315 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 5123 void removespace(char *str) {
|
Line 5328 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 5141 int decodemodel ( char model[], int last
|
Line 5346 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 5153 int decodemodel ( char model[], int last
|
Line 5358 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 5290 int decodemodel ( char model[], int last
|
Line 5495 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); |
| } |
} |
| |
|
| 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 5306 calandcheckages(int imx, int maxwav, dou
|
Line 5511 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 5325 calandcheckages(int imx, int maxwav, dou
|
Line 5530 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 5337 calandcheckages(int imx, int maxwav, dou
|
Line 5542 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 5367 calandcheckages(int imx, int maxwav, dou
|
Line 5573 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 5387 calandcheckages(int imx, int maxwav, dou
|
Line 5593 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); |
| } |
} |
| |
|
| |
#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"*/ |
| |
|
| |
#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 */ |
| |
} |
| |
|
| |
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"); |
| |
} |
| |
/*}*/ |
| |
} |
| |
} |
| |
|
| |
|
| /***********************************************/ |
/***********************************************/ |
| /**************** Main Program *****************/ |
/**************** Main Program *****************/ |
|
Line 5407 int main(int argc, char *argv[])
|
Line 5913 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 5447 int main(int argc, char *argv[])
|
Line 5950 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 5458 int main(int argc, char *argv[])
|
Line 5961 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 5546 int main(int argc, char *argv[])
|
Line 6048 int main(int argc, char *argv[])
|
| 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 5573 int main(int argc, char *argv[])
|
Line 6075 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 5741 run imach with mle=-1 to get a correct t
|
Line 6245 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 6524 Interval (in months) between two waves:
|
Line 7028 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 6532 Interval (in months) between two waves:
|
Line 7038 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 6776 Interval (in months) between two waves:
|
Line 7283 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); |
![[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