|
version 1.179, 2015/01/04 09:57:06
|
version 1.186, 2015/04/23 12:01:52
|
|
Line 1
|
Line 1
|
| /* $Id$ |
/* $Id$ |
| $State$ |
$State$ |
| $Log$ |
$Log$ |
| |
Revision 1.186 2015/04/23 12:01:52 brouard |
| |
Summary: V1*age is working now, version 0.98q1 |
| |
|
| |
Some codes had been disabled in order to simplify and Vn*age was |
| |
working in the optimization phase, ie, giving correct MLE parameters, |
| |
but, as usual, outputs were not correct and program core dumped. |
| |
|
| |
Revision 1.185 2015/03/11 13:26:42 brouard |
| |
Summary: Inclusion of compile and links command line for Intel Compiler |
| |
|
| |
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 |
Revision 1.179 2015/01/04 09:57:06 brouard |
| Summary: back to OS/X |
Summary: back to OS/X |
| |
|
|
Line 553
|
Line 586
|
| */ |
*/ |
| |
|
| #define POWELL /* Instead of NLOPT */ |
#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> |
|
Line 638 typedef struct {
|
Line 673 typedef struct {
|
| /* $Id$ */ |
/* $Id$ */ |
| /* $State$ */ |
/* $State$ */ |
| |
|
| char version[]="Imach version 0.98p, December 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015"; |
char version[]="Imach version 0.98q1, April 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 754 static double maxarg1,maxarg2;
|
Line 789 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 */ |
/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */ |
| /* #define mytinydouble 1.0e-16 */ |
#define mytinydouble 1.0e-16 |
| /* #define DEQUAL(a,b) (fabs((a)-(b))<mytinydouble) */ |
/* #define DEQUAL(a,b) (fabs((a)-(b))<mytinydouble) */ |
| /* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/mynrutils_8h-source.html */ |
/* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/mynrutils_8h-source.html */ |
| /* static double dsqrarg; */ |
/* static double dsqrarg; */ |
|
Line 802 static int split( char *path, char *dirc
|
Line 837 static int split( char *path, char *dirc
|
| the name of the file (name), its extension only (ext) and its first part of the name (finame) |
the name of the file (name), its extension only (ext) and its first part of the name (finame) |
| */ |
*/ |
| char *ss; /* pointer */ |
char *ss; /* pointer */ |
| int l1, l2; /* length counters */ |
int l1=0, l2=0; /* length counters */ |
| |
|
| l1 = strlen(path ); /* length of path */ |
l1 = strlen(path ); /* length of path */ |
| if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); |
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH ); |
|
Line 813 static int split( char *path, char *dirc
|
Line 848 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 824 static int split( char *path, char *dirc
|
Line 863 static int split( char *path, char *dirc
|
| if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); |
if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); |
| strcpy( name, ss ); /* save file name */ |
strcpy( name, ss ); /* save file name */ |
| strncpy( dirc, path, l1 - l2 ); /* now the directory */ |
strncpy( dirc, path, l1 - l2 ); /* now the directory */ |
| dirc[l1-l2] = 0; /* add zero */ |
dirc[l1-l2] = '\0'; /* add zero */ |
| printf(" DIRC2 = %s \n",dirc); |
printf(" DIRC2 = %s \n",dirc); |
| } |
} |
| /* We add a separator at the end of dirc if not exists */ |
/* We add a separator at the end of dirc if not exists */ |
|
Line 1270 double brent(double ax, double bx, doubl
|
Line 1309 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 1293 double brent(double ax, double bx, doubl
|
Line 1332 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 1301 void mnbrak(double *ax, double *bx, doub
|
Line 1344 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 1320 void mnbrak(double *ax, double *bx, doub
|
Line 1367 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 1401 void powell(double p[], double **xi, int
|
Line 1500 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 1460 void powell(double p[], double **xi, int
|
Line 1560 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 1481 void powell(double p[], double **xi, int
|
Line 1586 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 1513 void powell(double p[], double **xi, int
|
Line 1618 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 1543 void powell(double p[], double **xi, int
|
Line 1649 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 %d :\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 1598 double **prevalim(double **prlim, int nl
|
Line 1714 double **prevalim(double **prlim, int nl
|
| cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; |
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; |
| /*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/ |
/*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, k, codtab[ij][Tvar[k]]);*/ |
| } |
} |
| /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
/*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
| /* for (k=1; k<=cptcovprod;k++) /\* Useless *\/ */ |
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2]; |
| /* cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; */ |
for (k=1; k<=cptcovprod;k++) /* Useless */ |
| |
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; |
| |
|
| /*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ |
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/ |
| /*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ |
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/ |
|
Line 1770 double ***hpxij(double ***po, int nhstep
|
Line 1887 double ***hpxij(double ***po, int nhstep
|
| cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM; |
cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM; |
| for (k=1; k<=cptcovn;k++) |
for (k=1; k<=cptcovn;k++) |
| cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; |
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; |
| for (k=1; k<=cptcovage;k++) |
for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */ |
| cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; |
/* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
| |
cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; |
| for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ |
for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */ |
| cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; |
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; |
| |
|
|
Line 1913 double func( double *x)
|
Line 2031 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 1945 double func( double *x)
|
Line 2082 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 2280 void mlikeli(FILE *ficres,double p[], in
|
Line 2421 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 2975 void tricode(int *Tvar, int **nbcode, in
|
Line 3116 void tricode(int *Tvar, int **nbcode, in
|
| for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ |
for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */ |
| |
|
| /* Loop on covariates without age and products */ |
/* Loop on covariates without age and products */ |
| for (j=1; j<=(cptcovs); j++) { /* model V1 + V2*age+ V3 + V3*V4 : V1 + V3 = 2 only */ |
for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */ |
| for (i=1; i<=imx; i++) { /* Lopp on individuals: reads the data file to get the maximum value of the |
for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the |
| modality of this covariate Vj*/ |
modality of this covariate Vj*/ |
| ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i |
ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i |
| * If product of Vn*Vm, still boolean *: |
* If product of Vn*Vm, still boolean *: |
|
Line 3013 void tricode(int *Tvar, int **nbcode, in
|
Line 3154 void tricode(int *Tvar, int **nbcode, in
|
| } /* Ndum[-1] number of undefined modalities */ |
} /* Ndum[-1] number of undefined modalities */ |
| |
|
| /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ |
/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */ |
| /* For covariate j, modalities could be 1, 2, 3, 4. If Ndum[2]=0 ncodemax[j] is not 4 but 3 */ |
/* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. |
| /* If Ndum[3}= 635; Ndum[4]=0; Ndum[5]=0; Ndum[6]=27; Ndum[7]=125; |
If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; |
| modmincovj=3; modmaxcovj = 7; |
modmincovj=3; modmaxcovj = 7; |
| There are only 3 modalities non empty (or 2 if 27 is too few) : ncodemax[j]=3; |
There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; |
| which will be coded 0, 1, 2 which in binary on 3-1 digits are 0=00 1=01, 2=10; defining two dummy |
which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; |
| variables V1_1 and V1_2. |
defining two dummy variables: variables V1_1 and V1_2. |
| nbcode[Tvar[j]][ij]=k; |
nbcode[Tvar[j]][ij]=k; |
| nbcode[Tvar[j]][1]=0; |
nbcode[Tvar[j]][1]=0; |
| nbcode[Tvar[j]][2]=1; |
nbcode[Tvar[j]][2]=1; |
|
Line 3029 void tricode(int *Tvar, int **nbcode, in
|
Line 3170 void tricode(int *Tvar, int **nbcode, in
|
| for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */ |
for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */ |
| /*recode from 0 */ |
/*recode from 0 */ |
| if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */ |
if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */ |
| nbcode[Tvar[j]][ij]=k; /* stores the modality in an array nbcode. |
nbcode[Tvar[j]][ij]=k; /* stores the modality k in an array nbcode. |
| k is a modality. If we have model=V1+V1*sex |
k is a modality. If we have model=V1+V1*sex |
| then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ |
then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ |
| ij++; |
ij++; |
|
Line 3906 To be simple, these graphs help to under
|
Line 4047 To be simple, these graphs help to under
|
| */ |
*/ |
| /* nbcode[1][1]=0 nbcode[1][2]=1;*/ |
/* nbcode[1][1]=0 nbcode[1][2]=1;*/ |
| } |
} |
| for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; |
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
| |
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; |
| for (k=1; k<=cptcovprod;k++) |
for (k=1; k<=cptcovprod;k++) |
| cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; |
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; |
| |
|
|
Line 4399 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 4541 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| fprintf(ficgp," exp(p%d+p%d*x",i,i+1); |
fprintf(ficgp," exp(p%d+p%d*x",i,i+1); |
| ij=1;/* To be checked else nbcode[0][0] wrong */ |
ij=1;/* To be checked else nbcode[0][0] wrong */ |
| for(j=3; j <=ncovmodel; j++) { |
for(j=3; j <=ncovmodel; j++) { |
| /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /\* Bug valgrind *\/ */ |
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */ |
| /* /\*fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);*\/ */ |
fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); |
| /* ij++; */ |
ij++; |
| /* } */ |
} |
| /* else */ |
else |
| fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]); |
fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]); |
| } |
} |
| fprintf(ficgp,")/(1"); |
fprintf(ficgp,")/(1"); |
|
Line 4412 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 4554 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
| fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); |
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); |
| ij=1; |
ij=1; |
| for(j=3; j <=ncovmodel; j++){ |
for(j=3; j <=ncovmodel; j++){ |
| /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { */ |
if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { |
| /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); */ |
fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); |
| /* ij++; */ |
ij++; |
| /* } */ |
} |
| /* else */ |
else |
| fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]); |
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]); |
| } |
} |
| fprintf(ficgp,")"); |
fprintf(ficgp,")"); |
|
Line 5297 int decodemodel ( char model[], int last
|
Line 5439 int decodemodel ( char model[], int last
|
| /* for (k=1; k<=cptcovn;k++) { */ |
/* for (k=1; k<=cptcovn;k++) { */ |
| /* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */ |
/* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */ |
| /* } */ |
/* } */ |
| /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */ |
| /* |
/* |
| * Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */ |
* Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */ |
| for(k=cptcovt; k>=1;k--) /**< Number of covariates */ |
for(k=cptcovt; k>=1;k--) /**< Number of covariates */ |
|
Line 5317 int decodemodel ( char model[], int last
|
Line 5459 int decodemodel ( char model[], int last
|
| cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ |
cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */ |
| Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2 */ |
Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2 */ |
| cptcovage++; /* Sums the number of covariates which include age as a product */ |
cptcovage++; /* Sums the number of covariates which include age as a product */ |
| Tage[cptcovage]=k; /* Tage[1] = 4 */ |
Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */ |
| /*printf("stre=%s ", stre);*/ |
/*printf("stre=%s ", stre);*/ |
| } else if (strcmp(strd,"age")==0) { /* or age*Vn */ |
} else if (strcmp(strd,"age")==0) { /* or age*Vn */ |
| cptcovprod--; |
cptcovprod--; |
|
Line 5517 BOOL IsWow64()
|
Line 5659 BOOL IsWow64()
|
| void syscompilerinfo() |
void syscompilerinfo() |
| { |
{ |
| /* #include "syscompilerinfo.h"*/ |
/* #include "syscompilerinfo.h"*/ |
| |
/* command line Intel compiler 32bit windows, XP compatible:*/ |
| |
/* /GS /W3 /Gy |
| |
/Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D |
| |
"_CONSOLE" /D "_LIB" /D "_USING_V110_SDK71_" /D "_UNICODE" /D |
| |
"UNICODE" /Qipo /Zc:forScope /Gd /Oi /MT /Fa"Release\" /EHsc /nologo |
| |
/Fo"Release\" /Qprof-dir "Release\" /Fp"Release\IMaCh.pch" |
| |
*/ |
| |
/* 64 bits */ |
| |
/* |
| |
/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" */ |
| |
/* Optimization are useless and O3 is slower than O2 */ |
| |
/* |
| |
/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" |
| |
*/ |
| |
/* Link is */ /* /OUT:"visual studio |
| |
2013\Projects\IMaCh\Release\IMaCh.exe" /MANIFEST /NXCOMPAT |
| |
/PDB:"visual studio |
| |
2013\Projects\IMaCh\Release\IMaCh.pdb" /DYNAMICBASE |
| |
"kernel32.lib" "user32.lib" "gdi32.lib" "winspool.lib" |
| |
"comdlg32.lib" "advapi32.lib" "shell32.lib" "ole32.lib" |
| |
"oleaut32.lib" "uuid.lib" "odbc32.lib" "odbccp32.lib" |
| |
/MACHINE:X86 /OPT:REF /SAFESEH /INCREMENTAL:NO |
| |
/SUBSYSTEM:CONSOLE",5.01" /MANIFESTUAC:"level='asInvoker' |
| |
uiAccess='false'" |
| |
/ManifestFile:"Release\IMaCh.exe.intermediate.manifest" /OPT:ICF |
| |
/NOLOGO /TLBID:1 |
| |
*/ |
| #if defined __INTEL_COMPILER |
#if defined __INTEL_COMPILER |
| #if defined(__GNUC__) |
#if defined(__GNUC__) |
| struct utsname sysInfo; /* For Intel on Linux and OS/X */ |
struct utsname sysInfo; /* For Intel on Linux and OS/X */ |
|
Line 5641 void syscompilerinfo()
|
Line 5816 void syscompilerinfo()
|
| |
|
| } |
} |
| |
|
| |
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 5748 int main(int argc, char *argv[])
|
Line 6062 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 5784 int main(int argc, char *argv[])
|
Line 6102 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){ |
|
Line 5801 int main(int argc, char *argv[])
|
Line 6124 int main(int argc, char *argv[])
|
| |
|
| /*-------- arguments in the command line --------*/ |
/*-------- arguments in the command line --------*/ |
| |
|
| /* Log file */ |
/* Main Log file */ |
| strcat(filelog, optionfilefiname); |
strcat(filelog, optionfilefiname); |
| strcat(filelog,".log"); /* */ |
strcat(filelog,".log"); /* */ |
| if((ficlog=fopen(filelog,"w"))==NULL) { |
if((ficlog=fopen(filelog,"w"))==NULL) { |
|
Line 5830 int main(int argc, char *argv[])
|
Line 6153 int main(int argc, char *argv[])
|
| strcat(fileres, optionfilefiname); |
strcat(fileres, optionfilefiname); |
| strcat(fileres,".txt"); /* Other files have txt extension */ |
strcat(fileres,".txt"); /* Other files have txt extension */ |
| |
|
| /*---------arguments file --------*/ |
/* Main ---------arguments file --------*/ |
| |
|
| if((ficpar=fopen(optionfile,"r"))==NULL) { |
if((ficpar=fopen(optionfile,"r"))==NULL) { |
| printf("Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno)); |
printf("Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno)); |
|
Line 5865 int main(int argc, char *argv[])
|
Line 6188 int main(int argc, char *argv[])
|
| |
|
| fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); |
fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); |
| numlinepar++; |
numlinepar++; |
| printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); |
/* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */ |
| |
printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n",title, datafile, lastobs, firstpass,lastpass); |
| |
/* |
| |
|
| |
|
| |
|
| |
*/ |
| |
printf("\nftol=%e \n", ftol); |
| |
printf("stepm=%d \n", stepm); |
| |
printf("ncovcol=%d nlstate=%d \n", ncovcol, nlstate); |
| |
printf("ndeath=%d maxwav=%d mle=%d weight=%d\n", ndeath, maxwav, mle, weightopt); |
| |
printf("model=%s\n",model); |
| fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
| fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
| fflush(ficlog); |
fflush(ficlog); |
|
Line 5913 int main(int argc, char *argv[])
|
Line 6247 int main(int argc, char *argv[])
|
| goto end; |
goto end; |
| exit(0); |
exit(0); |
| } |
} |
| else if(mle==-3) { |
else if(mle==-3) { /* Main Wizard */ |
| prwizard(ncovmodel, nlstate, ndeath, model, ficparo); |
prwizard(ncovmodel, nlstate, ndeath, model, ficparo); |
| printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); |
printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); |
| fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); |
fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso); |
|
Line 6073 run imach with mle=-1 to get a correct t
|
Line 6407 run imach with mle=-1 to get a correct t
|
| fprintf(ficres,"#%s\n",version); |
fprintf(ficres,"#%s\n",version); |
| } /* End of mle != -3 */ |
} /* End of mle != -3 */ |
| |
|
| |
/* Main data |
| |
*/ |
| n= lastobs; |
n= lastobs; |
| num=lvector(1,n); |
num=lvector(1,n); |
| moisnais=vector(1,n); |
moisnais=vector(1,n); |
|
Line 6124 run imach with mle=-1 to get a correct t
|
Line 6459 run imach with mle=-1 to get a correct t
|
| Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 |
Tage[1=V3*age]= 4; Tage[2=age*V4] = 3 |
| */ |
*/ |
| |
|
| |
/* Main decodemodel */ |
| |
|
| if(decodemodel(model, lastobs) == 1) |
if(decodemodel(model, lastobs) == 1) |
| goto end; |
goto end; |
| |
|
|
Line 6169 run imach with mle=-1 to get a correct t
|
Line 6506 run imach with mle=-1 to get a correct t
|
| Ndum =ivector(-1,NCOVMAX); |
Ndum =ivector(-1,NCOVMAX); |
| if (ncovmodel > 2) |
if (ncovmodel > 2) |
| tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ |
tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */ |
| |
/* Nbcode gives the value of the lth modality of jth covariate, in |
| |
V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/ |
| |
/* 1 to ncodemax[j] is the maximum value of this jth covariate */ |
| |
|
| codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */ |
codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */ |
| /*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtab[100][10]);*/ |
/*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtab[100][10]);*/ |
| |
/* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/ |
| h=0; |
h=0; |
| |
|
| |
|
|
Line 6188 run imach with mle=-1 to get a correct t
|
Line 6529 run imach with mle=-1 to get a correct t
|
| if (h>m) |
if (h>m) |
| h=1; |
h=1; |
| /**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 |
/**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1 |
| * h 1 2 3 4 |
* For k=4 covariates, h goes from 1 to 2**k |
| |
* codtabm(h,k)= 1 & (h-1) >> (k-1) ; |
| |
* h\k 1 2 3 4 |
| *______________________________ |
*______________________________ |
| * 1 i=1 1 i=1 1 i=1 1 i=1 1 |
* 1 i=1 1 i=1 1 i=1 1 i=1 1 |
| * 2 2 1 1 1 |
* 2 2 1 1 1 |
|
Line 6208 run imach with mle=-1 to get a correct t
|
Line 6551 run imach with mle=-1 to get a correct t
|
| * 16 2 2 2 1 |
* 16 2 2 2 1 |
| */ |
*/ |
| codtab[h][k]=j; |
codtab[h][k]=j; |
| |
/* codtab[12][3]=1; */ |
| /*codtab[h][Tvar[k]]=j;*/ |
/*codtab[h][Tvar[k]]=j;*/ |
| printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); |
printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]); |
| } |
} |
|
Line 6228 run imach with mle=-1 to get a correct t
|
Line 6572 run imach with mle=-1 to get a correct t
|
| |
|
| |
|
| |
|
| /*------------ gnuplot -------------*/ |
/* Initialisation of ----------- gnuplot -------------*/ |
| strcpy(optionfilegnuplot,optionfilefiname); |
strcpy(optionfilegnuplot,optionfilefiname); |
| if(mle==-3) |
if(mle==-3) |
| strcat(optionfilegnuplot,"-mort"); |
strcat(optionfilegnuplot,"-mort"); |
|
Line 6244 run imach with mle=-1 to get a correct t
|
Line 6588 run imach with mle=-1 to get a correct t
|
| fprintf(ficgp,"set datafile missing 'NaNq'\n"); |
fprintf(ficgp,"set datafile missing 'NaNq'\n"); |
| } |
} |
| /* fclose(ficgp);*/ |
/* fclose(ficgp);*/ |
| /*--------- index.htm --------*/ |
|
| |
|
| |
/* Initialisation of --------- index.htm --------*/ |
| |
|
| strcpy(optionfilehtm,optionfilefiname); /* Main html file */ |
strcpy(optionfilehtm,optionfilefiname); /* Main html file */ |
| if(mle==-3) |
if(mle==-3) |
|
Line 6286 Title=%s <br>Datafile=%s Firstpass=%d La
|
Line 6632 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 6309 Interval (in months) between two waves:
|
Line 6660 Interval (in months) between two waves:
|
| p=param[1][1]; /* *(*(*(param +1)+1)+0) */ |
p=param[1][1]; /* *(*(*(param +1)+1)+0) */ |
| |
|
| globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ |
globpr=0; /* To get the number ipmx of contributions and the sum of weights*/ |
| |
/* For mortality only */ |
| if (mle==-3){ |
if (mle==-3){ |
| ximort=matrix(1,NDIM,1,NDIM); |
ximort=matrix(1,NDIM,1,NDIM); |
| /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ |
/* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */ |
| cens=ivector(1,n); |
cens=ivector(1,n); |
| ageexmed=vector(1,n); |
ageexmed=vector(1,n); |
| agecens=vector(1,n); |
agecens=vector(1,n); |
|
Line 6402 Interval (in months) between two waves:
|
Line 6753 Interval (in months) between two waves:
|
| |
|
| /* Initialize method and iterate */ |
/* Initialize method and iterate */ |
| /* p[1]=0.0268; p[NDIM]=0.083; */ |
/* p[1]=0.0268; p[NDIM]=0.083; */ |
| /* gsl_vector_set(x, 0, 0.0268); */ |
/* gsl_vector_set(x, 0, 0.0268); */ |
| /* gsl_vector_set(x, 1, 0.083); */ |
/* gsl_vector_set(x, 1, 0.083); */ |
| gsl_vector_set(x, 0, p[1]); |
gsl_vector_set(x, 0, p[1]); |
| gsl_vector_set(x, 1, p[2]); |
gsl_vector_set(x, 1, p[2]); |
| |
|
|
Line 6520 Interval (in months) between two waves:
|
Line 6871 Interval (in months) between two waves:
|
| free_ivector(dcwave,1,n); |
free_ivector(dcwave,1,n); |
| free_matrix(ximort,1,NDIM,1,NDIM); |
free_matrix(ximort,1,NDIM,1,NDIM); |
| #endif |
#endif |
| } /* Endof if mle==-3 */ |
} /* Endof if mle==-3 mortality only */ |
| |
/* Standard maximisation */ |
| else{ /* For mle >=1 */ |
else{ /* For mle >=1 */ |
| globpr=0;/* debug */ |
globpr=0;/* debug */ |
| |
/* Computes likelihood for initial parameters */ |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
| printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
| for (k=1; k<=npar;k++) |
for (k=1; k<=npar;k++) |
| printf(" %d %8.5f",k,p[k]); |
printf(" %d %8.5f",k,p[k]); |
| printf("\n"); |
printf("\n"); |
| globpr=1; /* to print the contributions */ |
globpr=1; /* again, to print the contributions */ |
| likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */ |
| printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw); |
| for (k=1; k<=npar;k++) |
for (k=1; k<=npar;k++) |
| printf(" %d %8.5f",k,p[k]); |
printf(" %d %8.5f",k,p[k]); |
| printf("\n"); |
printf("\n"); |
| if(mle>=1){ /* Could be 1 or 2 */ |
if(mle>=1){ /* Could be 1 or 2, Real Maximisation */ |
| mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); |
mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func); |
| } |
} |
| |
|
|
Line 6696 Interval (in months) between two waves:
|
Line 7048 Interval (in months) between two waves:
|
| fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); |
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n"); |
| fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
| fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm); |
| |
|
| |
/* Other stuffs, more or less useful */ |
| while((c=getc(ficpar))=='#' && c!= EOF){ |
while((c=getc(ficpar))=='#' && c!= EOF){ |
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| fgets(line, MAXLINE, ficpar); |
fgets(line, MAXLINE, ficpar); |
|
Line 6769 Interval (in months) between two waves:
|
Line 7122 Interval (in months) between two waves:
|
| fclose(ficres); |
fclose(ficres); |
| |
|
| |
|
| |
/* Other results (useful)*/ |
| |
|
| |
|
| /*--------------- 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 6778 Interval (in months) between two waves:
|
Line 7136 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 6805 Interval (in months) between two waves:
|
Line 7164 Interval (in months) between two waves:
|
| /* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ |
/* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */ |
| /* } */ |
/* } */ |
| } |
} |
| |
|
| |
/* ------ Other prevalence ratios------------ */ |
| |
|
| /* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ |
/* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */ |
| |
|
|
Line 7080 Interval (in months) between two waves:
|
Line 7440 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"); |
![[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