version 1.185, 2015/03/11 13:26:42
|
version 1.199, 2015/09/07 14:09:23
|
Line 1
|
Line 1
|
/* $Id$ |
/* $Id$ |
$State$ |
$State$ |
$Log$ |
$Log$ |
|
Revision 1.199 2015/09/07 14:09:23 brouard |
|
Summary: 0.98q6 changing default small png format for graph to vectorized svg. |
|
|
|
Revision 1.198 2015/09/03 07:14:39 brouard |
|
Summary: 0.98q5 Flavia |
|
|
|
Revision 1.197 2015/09/01 18:24:39 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.196 2015/08/18 23:17:52 brouard |
|
Summary: 0.98q5 |
|
|
|
Revision 1.195 2015/08/18 16:28:39 brouard |
|
Summary: Adding a hack for testing purpose |
|
|
|
After reading the title, ftol and model lines, if the comment line has |
|
a q, starting with #q, the answer at the end of the run is quit. It |
|
permits to run test files in batch with ctest. The former workaround was |
|
$ echo q | imach foo.imach |
|
|
|
Revision 1.194 2015/08/18 13:32:00 brouard |
|
Summary: Adding error when the covariance matrix doesn't contain the exact number of lines required by the model line. |
|
|
|
Revision 1.193 2015/08/04 07:17:42 brouard |
|
Summary: 0.98q4 |
|
|
|
Revision 1.192 2015/07/16 16:49:02 brouard |
|
Summary: Fixing some outputs |
|
|
|
Revision 1.191 2015/07/14 10:00:33 brouard |
|
Summary: Some fixes |
|
|
|
Revision 1.190 2015/05/05 08:51:13 brouard |
|
Summary: Adding digits in output parameters (7 digits instead of 6) |
|
|
|
Fix 1+age+. |
|
|
|
Revision 1.189 2015/04/30 14:45:16 brouard |
|
Summary: 0.98q2 |
|
|
|
Revision 1.188 2015/04/30 08:27:53 brouard |
|
*** empty log message *** |
|
|
|
Revision 1.187 2015/04/29 09:11:15 brouard |
|
*** empty log message *** |
|
|
|
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 |
Revision 1.185 2015/03/11 13:26:42 brouard |
Summary: Inclusion of compile and links command line for Intel Compiler |
Summary: Inclusion of compile and links command line for Intel Compiler |
|
|
Line 578
|
Line 631
|
end |
end |
*/ |
*/ |
|
|
|
/* #define DEBUG */ |
|
/* #define DEBUGBRENT */ |
#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 POWELLF1F3 /* Skip test */ |
/* #define MNBRAKORIGINAL */ /* Don't use mnbrak fix */ |
/* #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 647 typedef struct {
|
Line 703 typedef struct {
|
#define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ |
#define NLSTATEMAX 8 /**< Maximum number of live states (for func) */ |
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
#define NDEATHMAX 8 /**< Maximum number of dead states (for func) */ |
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
#define NCOVMAX 20 /**< Maximum number of covariates, including generated covariates V1*V2 */ |
#define codtabm(h,k) 1 & (h-1) >> (k-1) ; |
#define codtabm(h,k) (1 & (h-1) >> (k-1))+1 |
#define MAXN 20000 |
#define MAXN 20000 |
#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 AGEOVERFLOW 1.e20 |
#define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */ |
#define AGEGOMP 10 /**< Minimal age for Gompertz adjustment */ |
#ifdef _WIN32 |
#ifdef _WIN32 |
#define DIRSEPARATOR '\\' |
#define DIRSEPARATOR '\\' |
Line 665 typedef struct {
|
Line 722 typedef struct {
|
|
|
/* $Id$ */ |
/* $Id$ */ |
/* $State$ */ |
/* $State$ */ |
|
#include "version.h" |
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 version[]=__IMACH_VERSION__; |
|
char copyright[]="September 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]; |
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ |
int erreur=0, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */ |
int nvar=0, nforce=0; /* Number of variables, number of forces */ |
int nagesqr=0, nforce=0; /* nagesqr=1 if model is including age*age, number of forces */ |
/* Number of covariates model=V2+V1+ V3*age+V2*V4 */ |
/* Number of covariates model=V2+V1+ V3*age+V2*V4 */ |
int cptcovn=0; /**< cptcovn number of covariates added in the model (excepting constant and age and age*product) */ |
int cptcovn=0; /**< cptcovn number of covariates added in the model (excepting constant and age and age*product) */ |
int cptcovt=0; /**< cptcovt number of covariates added in the model (excepting constant and age) */ |
int cptcovt=0; /**< cptcovt number of covariates added in the model (excepting constant and age) */ |
Line 801 int estepm;
|
Line 859 int estepm;
|
|
|
int m,nb; |
int m,nb; |
long *num; |
long *num; |
int firstpass=0, lastpass=4,*cod, *ncodemax, *Tage,*cens; |
int firstpass=0, lastpass=4,*cod, *cens; |
|
int *ncodemax; /* ncodemax[j]= Number of modalities of the j th |
|
covariate for which somebody answered excluding |
|
undefined. Usually 2: 0 and 1. */ |
|
int *ncodemaxwundef; /* ncodemax[j]= Number of modalities of the j th |
|
covariate for which somebody answered including |
|
undefined. Usually 3: -1, 0 and 1. */ |
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; |
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint; |
double **pmmij, ***probs; |
double **pmmij, ***probs; |
double *ageexmed,*agecens; |
double *ageexmed,*agecens; |
Line 812 int **s; /* Status */
|
Line 876 int **s; /* Status */
|
double *agedc; |
double *agedc; |
double **covar; /**< covar[j,i], value of jth covariate for individual i, |
double **covar; /**< covar[j,i], value of jth covariate for individual i, |
* covar=matrix(0,NCOVMAX,1,n); |
* covar=matrix(0,NCOVMAX,1,n); |
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; */ |
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*age; */ |
double idx; |
double idx; |
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ |
int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ |
|
int *Tage; |
int *Ndum; /** Freq of modality (tricode */ |
int *Ndum; /** Freq of modality (tricode */ |
int **codtab; /**< codtab=imatrix(1,100,1,10); */ |
int **codtab; /**< codtab=imatrix(1,100,1,10); */ |
int **Tvard, *Tprod, cptcovprod, *Tvaraff; |
int **Tvard, *Tprod, cptcovprod, *Tvaraff; |
Line 830 static int split( char *path, char *dirc
|
Line 895 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 856 static int split( char *path, char *dirc
|
Line 921 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 908 char *trimbb(char *out, char *in)
|
Line 973 char *trimbb(char *out, char *in)
|
return s; |
return s; |
} |
} |
|
|
|
/* char *substrchaine(char *out, char *in, char *chain) */ |
|
/* { */ |
|
/* /\* Substract chain 'chain' from 'in', return and output 'out' *\/ */ |
|
/* char *s, *t; */ |
|
/* t=in;s=out; */ |
|
/* while ((*in != *chain) && (*in != '\0')){ */ |
|
/* *out++ = *in++; */ |
|
/* } */ |
|
|
|
/* /\* *in matches *chain *\/ */ |
|
/* while ((*in++ == *chain++) && (*in != '\0')){ */ |
|
/* printf("*in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ |
|
/* } */ |
|
/* in--; chain--; */ |
|
/* while ( (*in != '\0')){ */ |
|
/* printf("Bef *in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ |
|
/* *out++ = *in++; */ |
|
/* printf("Aft *in = %c, *out= %c *chain= %c \n", *in, *out, *chain); */ |
|
/* } */ |
|
/* *out='\0'; */ |
|
/* out=s; */ |
|
/* return out; */ |
|
/* } */ |
|
char *substrchaine(char *out, char *in, char *chain) |
|
{ |
|
/* Substract chain 'chain' from 'in', return and output 'out' */ |
|
/* in="V1+V1*age+age*age+V2", chain="age*age" */ |
|
|
|
char *strloc; |
|
|
|
strcpy (out, in); |
|
strloc = strstr(out, chain); /* strloc points to out at age*age+V2 */ |
|
printf("Bef strloc=%s chain=%s out=%s \n", strloc, chain, out); |
|
if(strloc != NULL){ |
|
/* will affect out */ /* strloc+strlenc(chain)=+V2 */ /* Will also work in Unicode */ |
|
memmove(strloc,strloc+strlen(chain), strlen(strloc+strlen(chain))+1); |
|
/* strcpy (strloc, strloc +strlen(chain));*/ |
|
} |
|
printf("Aft strloc=%s chain=%s in=%s out=%s \n", strloc, chain, in, out); |
|
return out; |
|
} |
|
|
|
|
char *cutl(char *blocc, char *alocc, char *in, char occ) |
char *cutl(char *blocc, char *alocc, char *in, char occ) |
{ |
{ |
/* cuts string in into blocc and alocc where blocc ends before first occurence of char 'occ' |
/* cuts string in into blocc and alocc where blocc ends before FIRST occurence of char 'occ' |
and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') |
and alocc starts after first occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') |
gives blocc="abcdef2ghi" and alocc="j". |
gives blocc="abcdef" and alocc="ghi2j". |
If occ is not found blocc is null and alocc is equal to in. Returns blocc |
If occ is not found blocc is null and alocc is equal to in. Returns blocc |
*/ |
*/ |
char *s, *t; |
char *s, *t; |
Line 938 char *cutl(char *blocc, char *alocc, cha
|
Line 1046 char *cutl(char *blocc, char *alocc, cha
|
} |
} |
char *cutv(char *blocc, char *alocc, char *in, char occ) |
char *cutv(char *blocc, char *alocc, char *in, char occ) |
{ |
{ |
/* cuts string in into blocc and alocc where blocc ends before last occurence of char 'occ' |
/* cuts string in into blocc and alocc where blocc ends before LAST occurence of char 'occ' |
and alocc starts after last occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') |
and alocc starts after last occurence of char 'occ' : ex cutv(blocc,alocc,"abcdef2ghi2j",'2') |
gives blocc="abcdef2ghi" and alocc="j". |
gives blocc="abcdef2ghi" and alocc="j". |
If occ is not found blocc is null and alocc is equal to in. Returns alocc |
If occ is not found blocc is null and alocc is equal to in. Returns alocc |
Line 1249 double f1dim(double x)
|
Line 1357 double f1dim(double x)
|
|
|
/*****************brent *************************/ |
/*****************brent *************************/ |
double brent(double ax, double bx, double cx, double (*f)(double), double tol, double *xmin) |
double brent(double ax, double bx, double cx, double (*f)(double), double tol, double *xmin) |
{ |
{ |
|
/* Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is |
|
* between ax and cx, and f(bx) is less than both f(ax) and f(cx) ), this routine isolates |
|
* the minimum to a fractional precision of about tol using Brent’s method. The abscissa of |
|
* the minimum is returned as xmin, and the minimum function value is returned as brent , the |
|
* returned function value. |
|
*/ |
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; |
Line 1332 values at the three points, fa, fb , and
|
Line 1446 values at the three points, fa, fb , and
|
*/ |
*/ |
double ulim,u,r,q, dum; |
double ulim,u,r,q, dum; |
double fu; |
double fu; |
|
|
*fa=(*func)(*ax); |
double scale=10.; |
*fb=(*func)(*bx); |
int iterscale=0; |
|
|
|
*fa=(*func)(*ax); /* xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/ |
|
*fb=(*func)(*bx); /* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */ |
|
|
|
|
|
/* while(*fb != *fb){ /\* *ax should be ok, reducing distance to *ax *\/ */ |
|
/* printf("Warning mnbrak *fb = %lf, *bx=%lf *ax=%lf *fa==%lf iter=%d\n",*fb, *bx, *ax, *fa, iterscale++); */ |
|
/* *bx = *ax - (*ax - *bx)/scale; */ |
|
/* *fb=(*func)(*bx); /\* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) *\/ */ |
|
/* } */ |
|
|
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) |
Line 1367 values at the three points, fa, fb , and
|
Line 1492 values at the three points, fa, fb , and
|
#endif |
#endif |
#ifdef MNBRAKORIGINAL |
#ifdef MNBRAKORIGINAL |
#else |
#else |
if (fu > *fc) { |
/* if (fu > *fc) { */ |
#ifdef DEBUG |
/* #ifdef DEBUG */ |
printf("mnbrak4 fu > fc \n"); |
/* printf("mnbrak4 fu > fc \n"); */ |
fprintf(ficlog, "mnbrak4 fu > fc\n"); |
/* fprintf(ficlog, "mnbrak4 fu > fc\n"); */ |
#endif |
/* #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(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 *\/ */ |
/* /\* 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 */ |
/* dum=u; /\* Shifting c and u *\/ */ |
u = *cx; |
/* u = *cx; */ |
*cx = dum; |
/* *cx = dum; */ |
dum = fu; |
/* dum = fu; */ |
fu = *fc; |
/* fu = *fc; */ |
*fc =dum; |
/* *fc =dum; */ |
} else { /* end */ |
/* } 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; */ |
|
/* } */ |
#ifdef DEBUG |
#ifdef DEBUG |
printf("mnbrak3 fu < fc \n"); |
printf("mnbrak34 fu < or >= fc \n"); |
fprintf(ficlog, "mnbrak3 fu < fc\n"); |
fprintf(ficlog, "mnbrak34 fu < fc\n"); |
#endif |
#endif |
dum=u; /* Shifting c and u */ |
dum=u; /* Shifting c and u */ |
u = *cx; |
u = *cx; |
*cx = dum; |
*cx = dum; |
dum = fu; |
dum = fu; |
fu = *fc; |
fu = *fc; |
*fc =dum; |
*fc =dum; |
} |
|
#endif |
#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 |
#ifdef DEBUG |
Line 1451 void linmin(double p[], double xi[], int
|
Line 1586 void linmin(double p[], double xi[], int
|
int j; |
int j; |
double xx,xmin,bx,ax; |
double xx,xmin,bx,ax; |
double fx,fb,fa; |
double fx,fb,fa; |
|
|
|
double scale=10., axs, xxs, xxss; /* Scale added for infinity */ |
|
|
ncom=n; |
ncom=n; |
pcom=vector(1,n); |
pcom=vector(1,n); |
Line 1460 void linmin(double p[], double xi[], int
|
Line 1597 void linmin(double p[], double xi[], int
|
pcom[j]=p[j]; |
pcom[j]=p[j]; |
xicom[j]=xi[j]; |
xicom[j]=xi[j]; |
} |
} |
ax=0.0; |
|
xx=1.0; |
/* axs=0.0; */ |
mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Find a bracket a,x,b in direction n=xi ie xicom */ |
/* xxss=1; /\* 1 and using scale *\/ */ |
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */ |
xxs=1; |
|
/* do{ */ |
|
ax=0.; |
|
xx= xxs; |
|
mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Outputs: xtx[j]=pcom[j]+(*xx)*xicom[j]; fx=f(xtx[j]) */ |
|
/* brackets with inputs ax=0 and xx=1, but points, pcom=p, and directions values, xicom=xi, are sent via f1dim(x) */ |
|
/* xt[x,j]=pcom[j]+x*xicom[j] f(ax) = f(xt(a,j=1,n)) = f(p(j) + 0 * xi(j)) and f(xx) = f(xt(x, j=1,n)) = f(p(j) + 1 * xi(j)) */ |
|
/* Outputs: fa=f(p(j)) and fx=f(p(j) + xxs * xi(j) ) and f(bx)= f(p(j)+ bx* xi(j)) */ |
|
/* Given input ax=axs and xx=xxs, xx might be too far from ax to get a finite f(xx) */ |
|
/* Searches on line, outputs (ax, xx, bx) such that fx < min(fa and fb) */ |
|
/* Find a bracket a,x,b in direction n=xi ie xicom, order may change. Scale is [0:xxs*xi[j]] et non plus [0:xi[j]]*/ |
|
/* if (fx != fx){ */ |
|
/* xxs=xxs/scale; /\* Trying a smaller xx, closer to initial ax=0 *\/ */ |
|
/* printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n", axs, xxs, fx,fb, fa, xx, ax, bx); */ |
|
/* } */ |
|
/* }while(fx != fx); */ |
|
|
|
#ifdef DEBUGLINMIN |
|
printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n", ax,xx,bx,fa,fx,fb); |
|
#endif |
|
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/ |
|
/* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */ |
|
/* fmin = f(p[j] + xmin * xi[j]) */ |
|
/* P+lambda n in that direction (lambdamin), with TOL between abscisses */ |
|
/* f1dim(xmin): for (j=1;j<=ncom;j++) xt[j]=pcom[j]+xmin*xicom[j]; */ |
#ifdef DEBUG |
#ifdef DEBUG |
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); |
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); |
fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); |
fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin); |
#endif |
#endif |
|
#ifdef DEBUGLINMIN |
|
printf("linmin end "); |
|
#endif |
for (j=1;j<=n;j++) { |
for (j=1;j<=n;j++) { |
xi[j] *= xmin; |
/* printf(" before xi[%d]=%12.8f", j,xi[j]); */ |
p[j] += xi[j]; |
xi[j] *= xmin; /* xi rescaled by xmin: if xmin=-1.237 and xi=(1,0,...,0) xi=(-1.237,0,...,0) */ |
|
/* if(xxs <1.0) */ |
|
/* printf(" after xi[%d]=%12.8f, xmin=%12.8f, ax=%12.8f, xx=%12.8f, bx=%12.8f, xxs=%12.8f", j,xi[j], xmin, ax, xx, bx,xxs ); */ |
|
p[j] += xi[j]; /* Parameters values are updated accordingly */ |
} |
} |
|
/* printf("\n"); */ |
|
#ifdef DEBUGLINMIN |
|
printf("Comparing last *frec(xmin=%12.8f)=%12.8f from Brent and frec(0.)=%12.8f \n", xmin, *fret, (*func)(p)); |
|
for (j=1;j<=n;j++) { |
|
printf(" xi[%d]= %12.7f p[%d]= %12.7f",j,xi[j],j,p[j]); |
|
if(j % ncovmodel == 0) |
|
printf("\n"); |
|
} |
|
#endif |
free_vector(xicom,1,n); |
free_vector(xicom,1,n); |
free_vector(pcom,1,n); |
free_vector(pcom,1,n); |
} |
} |
Line 1506 void powell(double p[], double **xi, int
|
Line 1682 void powell(double p[], double **xi, int
|
for (j=1;j<=n;j++) pt[j]=p[j]; |
for (j=1;j<=n;j++) pt[j]=p[j]; |
rcurr_time = time(NULL); |
rcurr_time = time(NULL); |
for (*iter=1;;++(*iter)) { |
for (*iter=1;;++(*iter)) { |
fp=(*fret); |
fp=(*fret); /* From former iteration or initial value */ |
ibig=0; |
ibig=0; |
del=0.0; |
del=0.0; |
rlast_time=rcurr_time; |
rlast_time=rcurr_time; |
Line 1516 void powell(double p[], double **xi, int
|
Line 1692 void powell(double p[], double **xi, int
|
printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); |
printf("\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); |
fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); |
fprintf(ficlog,"\nPowell iter=%d -2*LL=%.12f %ld sec. %ld sec.",*iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); |
/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ |
/* fprintf(ficrespow,"%d %.12f %ld",*iter,*fret,curr_time.tm_sec-start_time.tm_sec); */ |
for (i=1;i<=n;i++) { |
for (i=1;i<=n;i++) { |
printf(" %d %.12f",i, p[i]); |
printf(" %d %.12f",i, p[i]); |
fprintf(ficlog," %d %.12lf",i, p[i]); |
fprintf(ficlog," %d %.12lf",i, p[i]); |
fprintf(ficrespow," %.12lf", p[i]); |
fprintf(ficrespow," %.12lf", p[i]); |
Line 1544 void powell(double p[], double **xi, int
|
Line 1720 void powell(double p[], double **xi, int
|
fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); |
fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr); |
} |
} |
} |
} |
for (i=1;i<=n;i++) { |
for (i=1;i<=n;i++) { /* For each direction i */ |
for (j=1;j<=n;j++) xit[j]=xi[j][i]; |
for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */ |
fptt=(*fret); |
fptt=(*fret); |
#ifdef DEBUG |
#ifdef DEBUG |
printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); |
printf("fret=%lf, %lf, %lf \n", *fret, *fret, *fret); |
fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); |
fprintf(ficlog, "fret=%lf, %lf, %lf \n", *fret, *fret, *fret); |
#endif |
#endif |
printf("%d",i);fflush(stdout); |
printf("%d",i);fflush(stdout); /* print direction (parameter) i */ |
fprintf(ficlog,"%d",i);fflush(ficlog); |
fprintf(ficlog,"%d",i);fflush(ficlog); |
linmin(p,xit,n,fret,func); /* xit[n] has been loaded for direction i */ |
linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/ |
if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions |
/* Outputs are fret(new point p) p is updated and xit rescaled */ |
because that direction will be replaced unless the gain del is small |
if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions */ |
in comparison with the 'probable' gain, mu^2, with the last average direction. |
/* because that direction will be replaced unless the gain del is small */ |
Unless the n directions are conjugate some gain in the determinant may be obtained |
/* in comparison with the 'probable' gain, mu^2, with the last average direction. */ |
with the new 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 1578 void powell(double p[], double **xi, int
|
Line 1754 void powell(double p[], double **xi, int
|
printf("\n"); |
printf("\n"); |
fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
#endif |
#endif |
} /* end i */ |
} /* end loop on each direction i */ |
|
/* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ |
|
/* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ |
|
/* New value of last point Pn is not computed, P(n-1) */ |
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */ |
if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */ |
|
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ |
|
/* By adding age*age in a model, the new -2LL should be lower and the difference follows a */ |
|
/* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */ |
|
/* decreased of more than 3.84 */ |
|
/* By adding age*age and V1*age the gain (-2LL) should be more than 5.99 (ddl=2) */ |
|
/* By using V1+V2+V3, the gain should be 7.82, compared with basic 1+age. */ |
|
/* By adding 10 parameters more the gain should be 18.31 */ |
|
|
|
/* Starting the program with initial values given by a former maximization will simply change */ |
|
/* the scales of the directions and the directions, because the are reset to canonical directions */ |
|
/* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */ |
|
/* under the tolerance value. If the tolerance is very small 1.e-9, it could last long. */ |
#ifdef DEBUG |
#ifdef DEBUG |
int k[2],l; |
int k[2],l; |
k[0]=1; |
k[0]=1; |
Line 1609 void powell(double p[], double **xi, int
|
Line 1800 void powell(double p[], double **xi, int
|
free_vector(ptt,1,n); |
free_vector(ptt,1,n); |
free_vector(pt,1,n); |
free_vector(pt,1,n); |
return; |
return; |
} |
} /* enough precision */ |
if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); |
if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); |
for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */ |
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]; |
Line 1617 void powell(double p[], double **xi, int
|
Line 1808 void powell(double p[], double **xi, int
|
pt[j]=p[j]; |
pt[j]=p[j]; |
} |
} |
fptt=(*func)(ptt); /* f_3 */ |
fptt=(*func)(ptt); /* f_3 */ |
|
#ifdef POWELLF1F3 |
|
#else |
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 */ |
|
#endif |
/* (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. */ |
Line 1646 void powell(double p[], double **xi, int
|
Line 1840 void powell(double p[], double **xi, int
|
if (t < 0.0) { /* Then we use it for new direction */ |
if (t < 0.0) { /* Then we use it for new direction */ |
#else |
#else |
if (directest*t < 0.0) { /* Contradiction between both tests */ |
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("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); |
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,"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); |
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 */ |
if (directest < 0.0) { /* Then we use it for new direction */ |
#endif |
#endif |
linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/ |
#ifdef DEBUGLINMIN |
|
printf("Before linmin in direction P%d-P0\n",n); |
|
for (j=1;j<=n;j++) { |
|
printf("Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); |
|
if(j % ncovmodel == 0) |
|
printf("\n"); |
|
} |
|
#endif |
|
linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/ |
|
#ifdef DEBUGLINMIN |
|
for (j=1;j<=n;j++) { |
|
printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]); |
|
if(j % ncovmodel == 0) |
|
printf("\n"); |
|
} |
|
#endif |
for (j=1;j<=n;j++) { |
for (j=1;j<=n;j++) { |
xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ |
xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */ |
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ |
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */ |
Line 1671 void powell(double p[], double **xi, int
|
Line 1880 void powell(double p[], double **xi, int
|
printf("\n"); |
printf("\n"); |
fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
#endif |
#endif |
} /* end of t negative */ |
} /* end of t or directest negative */ |
|
#ifdef POWELLF1F3 |
|
#else |
} /* end if (fptt < fp) */ |
} /* end if (fptt < fp) */ |
} |
#endif |
|
} /* loop iteration */ |
} |
} |
|
|
/**** Prevalence limit (stable or period prevalence) ****************/ |
/**** Prevalence limit (stable or period prevalence) ****************/ |
Line 1702 double **prevalim(double **prlim, int nl
|
Line 1914 double **prevalim(double **prlim, int nl
|
newm=savm; |
newm=savm; |
/* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
cov[2]=agefin; |
cov[2]=agefin; |
|
if(nagesqr==1) |
|
cov[3]= agefin*agefin;; |
for (k=1; k<=cptcovn;k++) { |
for (k=1; k<=cptcovn;k++) { |
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; |
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(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 codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(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+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(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+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(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 1863 double ***hpxij(double ***po, int nhstep
|
Line 2077 double ***hpxij(double ***po, int nhstep
|
int i, j, d, h, k; |
int i, j, d, h, k; |
double **out, cov[NCOVMAX+1]; |
double **out, cov[NCOVMAX+1]; |
double **newm; |
double **newm; |
|
double agexact; |
|
|
/* Hstepm could be zero and should return the unit matrix */ |
/* Hstepm could be zero and should return the unit matrix */ |
for (i=1;i<=nlstate+ndeath;i++) |
for (i=1;i<=nlstate+ndeath;i++) |
Line 1876 double ***hpxij(double ***po, int nhstep
|
Line 2091 double ***hpxij(double ***po, int nhstep
|
newm=savm; |
newm=savm; |
/* Covariates have to be included here again */ |
/* Covariates have to be included here again */ |
cov[1]=1.; |
cov[1]=1.; |
cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM; |
agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; |
|
cov[2]=agexact; |
|
if(nagesqr==1) |
|
cov[3]= agexact*agexact; |
for (k=1; k<=cptcovn;k++) |
for (k=1; k<=cptcovn;k++) |
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; |
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(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+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(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+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; |
|
|
|
|
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/ |
Line 1934 double func( double *x)
|
Line 2153 double func( double *x)
|
int s1, s2; |
int s1, s2; |
double bbh, survp; |
double bbh, survp; |
long ipmx; |
long ipmx; |
|
double agexact; |
/*extern weight */ |
/*extern weight */ |
/* We are differentiating ll according to initial status */ |
/* We are differentiating ll according to initial status */ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
Line 1955 double func( double *x)
|
Line 2175 double func( double *x)
|
to be observed in j being in i according to the model. |
to be observed in j being in i according to the model. |
*/ |
*/ |
for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */ |
for (k=1; k<=cptcovn;k++){ /* Simple and product covariates without age* products */ |
cov[2+k]=covar[Tvar[k]][i]; |
cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
} |
} |
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
/* In model V2+V1*V4+age*V3+V3*V2 Tvar[1] is V2, Tvar[2=V1*V4] |
is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] |
is 6, Tvar[3=age*V3] should not be computed because of age Tvar[4=V3*V2] |
Line 1968 double func( double *x)
|
Line 2188 double func( double *x)
|
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
|
cov[2]=agexact; |
|
if(nagesqr==1) |
|
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; /* Tage[kk] gives the data-covariate associated with age */ |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; /* Tage[kk] gives the data-covariate associated with age */ |
} |
} |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
Line 2081 double func( double *x)
|
Line 2304 double func( double *x)
|
} /* end of individual */ |
} /* end of individual */ |
} else if(mle==2){ |
} else if(mle==2){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
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++){ |
Line 2090 double func( double *x)
|
Line 2313 double func( double *x)
|
} |
} |
for(d=0; d<=dh[mi][i]; d++){ |
for(d=0; d<=dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
|
cov[2]=agexact; |
|
if(nagesqr==1) |
|
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
Line 2111 double func( double *x)
|
Line 2337 double func( double *x)
|
} /* end of individual */ |
} /* end of individual */ |
} else if(mle==3){ /* exponential inter-extrapolation */ |
} else if(mle==3){ /* exponential inter-extrapolation */ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
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++){ |
Line 2120 double func( double *x)
|
Line 2346 double func( double *x)
|
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
|
cov[2]=agexact; |
|
if(nagesqr==1) |
|
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
Line 2141 double func( double *x)
|
Line 2370 double func( double *x)
|
} /* end of individual */ |
} /* end of individual */ |
}else if (mle==4){ /* ml=4 no inter-extrapolation */ |
}else if (mle==4){ /* ml=4 no inter-extrapolation */ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
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++){ |
Line 2150 double func( double *x)
|
Line 2379 double func( double *x)
|
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
|
cov[2]=agexact; |
|
if(nagesqr==1) |
|
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
|
|
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
Line 2176 double func( double *x)
|
Line 2408 double func( double *x)
|
} /* end of individual */ |
} /* end of individual */ |
}else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ |
}else{ /* ml=5 no inter-extrapolation no jackson =0.8a */ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
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++){ |
Line 2185 double func( double *x)
|
Line 2417 double func( double *x)
|
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
|
cov[2]=agexact; |
|
if(nagesqr==1) |
|
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
|
|
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
Line 2223 double funcone( double *x)
|
Line 2458 double funcone( double *x)
|
double llt; |
double llt; |
int s1, s2; |
int s1, s2; |
double bbh, survp; |
double bbh, survp; |
|
double agexact; |
/*extern weight */ |
/*extern weight */ |
/* We are differentiating ll according to initial status */ |
/* We are differentiating ll according to initial status */ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/ |
Line 2234 double funcone( double *x)
|
Line 2470 double funcone( double *x)
|
for(k=1; k<=nlstate; k++) ll[k]=0.; |
for(k=1; k<=nlstate; k++) ll[k]=0.; |
|
|
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; |
for(mi=1; mi<= wav[i]-1; mi++){ |
for(mi=1; mi<= wav[i]-1; mi++){ |
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++){ |
Line 2243 double funcone( double *x)
|
Line 2479 double funcone( double *x)
|
} |
} |
for(d=0; d<dh[mi][i]; d++){ |
for(d=0; d<dh[mi][i]; d++){ |
newm=savm; |
newm=savm; |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM; |
|
cov[2]=agexact; |
|
if(nagesqr==1) |
|
cov[3]= agexact*agexact; |
for (kk=1; kk<=cptcovage;kk++) { |
for (kk=1; kk<=cptcovage;kk++) { |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact; |
} |
} |
|
|
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */ |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
Line 2680 void lubksb(double **a, int n, int *indx
|
Line 2920 void lubksb(double **a, int n, int *indx
|
|
|
void pstamp(FILE *fichier) |
void pstamp(FILE *fichier) |
{ |
{ |
fprintf(fichier,"# %s.%s\n#%s\n#%s\n# %s", optionfilefiname,optionfilext,version,fullversion,strstart); |
fprintf(fichier,"# %s.%s\n#IMaCh version %s, %s\n#%s\n# %s", optionfilefiname,optionfilext,version,copyright, fullversion, strstart); |
} |
} |
|
|
/************ Frequencies ********************/ |
/************ Frequencies ********************/ |
Line 2732 void freqsummary(char fileres[], int ia
|
Line 2972 void freqsummary(char fileres[], int ia
|
bool=1; |
bool=1; |
if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ |
for (z1=1; z1<=cptcoveff; z1++) |
for (z1=1; z1<=cptcoveff; z1++) |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]){ |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){ |
/* Tests if the value of each of the covariates of i is equal to filter j1 */ |
/* Tests if the value of each of the covariates of i is equal to filter j1 */ |
bool=0; |
bool=0; |
/* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtab[%d][%d]=%d, nbcode[Tvaraff][codtab[%d][%d]=%d, j1=%d\n", |
/* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", |
bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtab[j1][z1], |
bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1), |
j1,z1,nbcode[Tvaraff[z1]][codtab[j1][z1]],j1);*/ |
j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ |
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtab[7][3]=1 and nbcde[3][?]=1*/ |
/* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ |
} |
} |
} |
} |
|
|
Line 2767 void freqsummary(char fileres[], int ia
|
Line 3007 void freqsummary(char fileres[], int ia
|
pstamp(ficresp); |
pstamp(ficresp); |
if (cptcovn>0) { |
if (cptcovn>0) { |
fprintf(ficresp, "\n#********** Variable "); |
fprintf(ficresp, "\n#********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(ficresp, "**********\n#"); |
fprintf(ficresp, "**********\n#"); |
fprintf(ficlog, "\n#********** Variable "); |
fprintf(ficlog, "\n#********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(ficlog, "**********\n#"); |
fprintf(ficlog, "**********\n#"); |
} |
} |
for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
Line 2898 void prevalence(double ***probs, double
|
Line 3138 void prevalence(double ***probs, double
|
bool=1; |
bool=1; |
if (cptcovn>0) { |
if (cptcovn>0) { |
for (z1=1; z1<=cptcoveff; z1++) |
for (z1=1; z1<=cptcoveff; z1++) |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]) |
bool=0; |
bool=0; |
} |
} |
if (bool==1) { |
if (bool==1) { |
Line 3103 void tricode(int *Tvar, int **nbcode, in
|
Line 3343 void tricode(int *Tvar, int **nbcode, in
|
|
|
cptcoveff=0; |
cptcoveff=0; |
|
|
for (k=-1; k < maxncov; k++) Ndum[k]=0; |
|
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 (k=-1; k < maxncov; k++) Ndum[k]=0; |
|
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 3130 void tricode(int *Tvar, int **nbcode, in
|
Line 3370 void tricode(int *Tvar, int **nbcode, in
|
/* getting the maximum value of the modality of the covariate |
/* getting the maximum value of the modality of the covariate |
(should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and |
(should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and |
female is 1, then modmaxcovj=1.*/ |
female is 1, then modmaxcovj=1.*/ |
} |
} /* end for loop on individuals i */ |
printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); |
printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); |
|
fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", j, Tvar[j], modmincovj, modmaxcovj); |
cptcode=modmaxcovj; |
cptcode=modmaxcovj; |
/* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ |
/* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */ |
/*for (i=0; i<=cptcode; i++) {*/ |
/*for (i=0; i<=cptcode; i++) {*/ |
for (i=modmincovj; i<=modmaxcovj; i++) { /* i=-1 ? 0 and 1*//* For each value of the modality of model-cov j */ |
for (k=modmincovj; k<=modmaxcovj; k++) { /* k=-1 ? 0 and 1*//* For each value k of the modality of model-cov j */ |
printf("Frequencies of covariates %d V%d %d\n", j, Tvar[j], Ndum[i]); |
printf("Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); |
if( Ndum[i] != 0 ){ /* Counts if nobody answered, empty modality */ |
fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", j, Tvar[j], k, Ndum[k]); |
ncodemax[j]++; /* ncodemax[j]= Number of non-null modalities of the j th covariate. */ |
if( Ndum[k] != 0 ){ /* Counts if nobody answered modality k ie empty modality, we skip it and reorder */ |
|
if( k != -1){ |
|
ncodemax[j]++; /* ncodemax[j]= Number of modalities of the j th |
|
covariate for which somebody answered excluding |
|
undefined. Usually 2: 0 and 1. */ |
|
} |
|
ncodemaxwundef[j]++; /* ncodemax[j]= Number of modalities of the j th |
|
covariate for which somebody answered including |
|
undefined. Usually 3: -1, 0 and 1. */ |
} |
} |
/* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for |
/* In fact ncodemax[j]=2 (dichotom. variables only) but it could be more for |
historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ |
historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */ |
} /* 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; |
nbcode[Tvar[j]][3]=2; |
nbcode[Tvar[j]][3]=2; |
|
To be continued (not working yet). |
*/ |
*/ |
ij=1; /* ij is similar to i but can jumps over null modalities */ |
ij=0; /* ij is similar to i but can jump over null modalities */ |
for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 */ |
for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/ |
for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */ |
if (Ndum[i] == 0) { /* If nobody responded to this modality k */ |
/*recode from 0 */ |
break; |
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. |
ij++; |
k is a modality. If we have model=V1+V1*sex |
nbcode[Tvar[j]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/ |
then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */ |
cptcode = ij; /* New max modality for covar j */ |
ij++; |
} /* end of loop on modality i=-1 to 1 or more */ |
} |
|
if (ij > ncodemax[j]) break; |
/* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */ |
} /* end of loop on */ |
/* /\*recode from 0 *\/ */ |
} /* end of loop on modality */ |
/* 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; */ |
|
/* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode *\/ */ |
|
/* } */ |
|
/* /\* cptcode = ij; *\/ /\* New max modality for covar j *\/ */ |
|
/* if (ij > ncodemax[j]) { */ |
|
/* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ |
|
/* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */ |
|
/* break; */ |
|
/* } */ |
|
/* } /\* end of loop on modality k *\/ */ |
} /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ |
} /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/ |
|
|
for (k=-1; k< maxncov; k++) Ndum[k]=0; |
for (k=-1; k< maxncov; k++) Ndum[k]=0; |
|
|
for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */ |
for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually age*age */ |
/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ |
/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ |
ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ |
ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ |
Ndum[ij]++; |
Ndum[ij]++; /* Might be supersed V1 + V1*age */ |
} |
} |
|
|
ij=1; |
ij=0; |
for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ |
for (i=0; i<= maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */ |
/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ |
/*printf("Ndum[%d]=%d\n",i, Ndum[i]);*/ |
if((Ndum[i]!=0) && (i<=ncovcol)){ |
if((Ndum[i]!=0) && (i<=ncovcol)){ |
|
ij++; |
/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ |
/*printf("diff Ndum[%d]=%d\n",i, Ndum[i]);*/ |
Tvaraff[ij]=i; /*For printing (unclear) */ |
Tvaraff[ij]=i; /*For printing (unclear) */ |
ij++; |
}else{ |
}else |
/* Tvaraff[ij]=0; */ |
Tvaraff[ij]=0; |
} |
} |
} |
ij--; |
/* ij--; */ |
cptcoveff=ij; /*Number of total covariates*/ |
cptcoveff=ij; /*Number of total covariates*/ |
|
|
} |
} |
Line 3793 void varevsij(char optionfilefiname[], d
|
Line 4054 void varevsij(char optionfilefiname[], d
|
free_vector(gmp,nlstate+1,nlstate+ndeath); |
free_vector(gmp,nlstate+1,nlstate+ndeath); |
free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); |
free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath); |
free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ |
free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/ |
fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); |
/* fprintf(ficgp,"\nunset parametric;unset label; set ter png small size 320, 240"); */ |
|
fprintf(ficgp,"\nunset parametric;unset label; set ter svg size 640, 480"); |
/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ |
/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */ |
fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); |
fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";"); |
/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ |
/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */ |
Line 3803 void varevsij(char optionfilefiname[], d
|
Line 4065 void varevsij(char optionfilefiname[], d
|
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.svg\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); |
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit); |
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.svg\"> <br>\n", stepm,YEARM,digitp,digit); |
*/ |
*/ |
/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit); */ |
/* fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.svg\";replot;",digitp,optionfilefiname,digit); */ |
fprintf(ficgp,"\nset out \"%s%s.png\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); |
fprintf(ficgp,"\nset out \"%s%s.svg\";replot;\n",subdirf3(optionfilefiname,"varmuptjgr",digitp),digit); |
|
|
free_vector(xp,1,npar); |
free_vector(xp,1,npar); |
free_matrix(doldm,1,nlstate,1,nlstate); |
free_matrix(doldm,1,nlstate,1,nlstate); |
Line 3980 void varprob(char optionfilefiname[], do
|
Line 4242 void varprob(char optionfilefiname[], do
|
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n"); |
fprintf(fichtm,"\n"); |
fprintf(fichtm,"\n"); |
|
|
fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of pairs of step probabilities (drawings)</a></h4></li>\n",optionfilehtmcov); |
fprintf(fichtm,"\n<li><h4> <a href=\"%s\">Matrix of variance-covariance of one-step probabilities (drawings)</a></h4></br>this page is important in order to visualize confidence intervals and especially correlation between disability and recovery</li>\n",optionfilehtmcov); |
fprintf(fichtmcov,"\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n\ |
fprintf(fichtmcov,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Matrix of variance-covariance of pairs of step probabilities</h4>\n",optionfilehtmcov, optionfilehtmcov); |
file %s<br>\n",optionfilehtmcov); |
fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated \ |
fprintf(fichtmcov,"\nEllipsoids of confidence centered on point (p<inf>ij</inf>, p<inf>kl</inf>) are estimated\ |
|
and drawn. It helps understanding how is the covariance between two incidences.\ |
and drawn. It helps understanding how is the covariance between two incidences.\ |
They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n"); |
They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n"); |
fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \ |
fprintf(fichtmcov,"\n<br> Contour plot corresponding to x'cov<sup>-1</sup>x = 4 (where x is the column vector (pij,pkl)) are drawn. \ |
Line 4004 To be simple, these graphs help to under
|
Line 4265 To be simple, these graphs help to under
|
/*j1++;*/ |
/*j1++;*/ |
if (cptcovn>0) { |
if (cptcovn>0) { |
fprintf(ficresprob, "\n#********** Variable "); |
fprintf(ficresprob, "\n#********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(ficresprob, "**********\n#\n"); |
fprintf(ficresprob, "**********\n#\n"); |
fprintf(ficresprobcov, "\n#********** Variable "); |
fprintf(ficresprobcov, "\n#********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(ficresprobcov, "**********\n#\n"); |
fprintf(ficresprobcov, "**********\n#\n"); |
|
|
fprintf(ficgp, "\n#********** Variable "); |
fprintf(ficgp, "\n#********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(ficgp, "**********\n#\n"); |
fprintf(ficgp, "**********\n#\n"); |
|
|
|
|
fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable "); |
fprintf(fichtmcov, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtmcov, "**********\n<hr size=\"2\" color=\"#EC5E5E\">"); |
|
|
fprintf(ficresprobcor, "\n#********** Variable "); |
fprintf(ficresprobcor, "\n#********** Variable "); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]); |
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]); |
fprintf(ficresprobcor, "**********\n#"); |
fprintf(ficresprobcor, "**********\n#"); |
} |
} |
|
|
Line 4030 To be simple, these graphs help to under
|
Line 4291 To be simple, these graphs help to under
|
gm=vector(1,(nlstate)*(nlstate+ndeath)); |
gm=vector(1,(nlstate)*(nlstate+ndeath)); |
for (age=bage; age<=fage; age ++){ |
for (age=bage; age<=fage; age ++){ |
cov[2]=age; |
cov[2]=age; |
|
if(nagesqr==1) |
|
cov[3]= age*age; |
for (k=1; k<=cptcovn;k++) { |
for (k=1; k<=cptcovn;k++) { |
cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];/* j1 1 2 3 4 |
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(j1,Tvar[k])];/* j1 1 2 3 4 |
* 1 1 1 1 1 |
* 1 1 1 1 1 |
* 2 2 1 1 1 |
* 2 2 1 1 1 |
* 3 1 2 1 1 |
* 3 1 2 1 1 |
*/ |
*/ |
/* 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]]][codtabm(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+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])]*nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; |
|
|
|
|
for(theta=1; theta <=npar; theta++){ |
for(theta=1; theta <=npar; theta++){ |
Line 4190 To be simple, these graphs help to under
|
Line 4454 To be simple, these graphs help to under
|
first=0; |
first=0; |
fprintf(ficgp,"\nset parametric;unset label"); |
fprintf(ficgp,"\nset parametric;unset label"); |
fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); |
fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k1,l1,k2,l2); |
fprintf(ficgp,"\nset ter png small size 320, 240"); |
fprintf(ficgp,"\nset ter svg size 640, 480"); |
fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\ |
fprintf(fichtmcov,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup>\ |
:<a href=\"%s%d%1d%1d-%1d%1d.png\">\ |
:<a href=\"%s%d%1d%1d-%1d%1d.svg\">\ |
%s%d%1d%1d-%1d%1d.png</A>, ",k1,l1,k2,l2,\ |
%s%d%1d%1d-%1d%1d.svg</A>, ",k1,l1,k2,l2,\ |
subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\ |
subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2,\ |
subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); |
subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); |
fprintf(fichtmcov,"\n<br><img src=\"%s%d%1d%1d-%1d%1d.png\"> ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); |
fprintf(fichtmcov,"\n<br><img src=\"%s%d%1d%1d-%1d%1d.svg\"> ",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); |
fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12); |
fprintf(fichtmcov,"\n<br> Correlation at age %d (%.3f),",(int) age, c12); |
fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); |
fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.svg\"",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); |
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2); |
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2); |
fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\ |
fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)), %11.3e +%.3f*(%11.3e*%11.3e*cos(t)+%11.3e*%11.3e*sin(t)) not",\ |
Line 4215 To be simple, these graphs help to under
|
Line 4479 To be simple, these graphs help to under
|
}/* if first */ |
}/* if first */ |
} /* age mod 5 */ |
} /* age mod 5 */ |
} /* end loop age */ |
} /* end loop age */ |
fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.png\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); |
fprintf(ficgp,"\nset out \"%s%d%1d%1d-%1d%1d.svg\";replot;",subdirf2(optionfilefiname,"varpijgr"), j1,k1,l1,k2,l2); |
first=1; |
first=1; |
} /*l12 */ |
} /*l12 */ |
} /* k12 */ |
} /* k12 */ |
Line 4272 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 4536 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
|
|
jj1=0; |
jj1=0; |
for(k1=1; k1<=m;k1++){ |
for(k1=1; k1<=m;k1++){ |
for(i1=1; i1<=ncodemax[k1];i1++){ |
/* for(i1=1; i1<=ncodemax[k1];i1++){ */ |
jj1++; |
jj1++; |
if (cptcovn > 0) { |
if (cptcovn > 0) { |
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
for (cpt=1; cpt<=cptcoveff;cpt++) |
for (cpt=1; cpt<=cptcoveff;cpt++){ |
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); |
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); |
|
printf(" V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]);fflush(stdout); |
|
} |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
} |
} |
/* Pij */ |
/* Pij */ |
fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d_1.png\">%s%d_1.png</a><br> \ |
fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s%d_1.svg\">%s%d_1.svg</a><br> \ |
<img src=\"%s%d_1.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); |
<img src=\"%s%d_1.svg\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); |
/* Quasi-incidences */ |
/* Quasi-incidences */ |
fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ |
fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months\ |
before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d_2.png\">%s%d_2.png</a><br> \ |
before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: <a href=\"%s%d_2.svg\">%s%d_2.svg</a><br> \ |
<img src=\"%s%d_2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1); |
<img src=\"%s%d_2.svg\">",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 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> \ |
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.svg\">%s%d_%d.svg</a><br> \ |
<img src=\"%s%d_%d.png\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1); |
<img src=\"%s%d_%d.svg\">", 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.svg\">%s%d%d.svg</a> <br> \ |
<img src=\"%s%d%d.png\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); |
<img src=\"%s%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1,subdirf2(optionfilefiname,"exp"),cpt,jj1); |
} |
} |
} /* end i1 */ |
/* } /\* end i1 *\/ */ |
}/* End k1 */ |
}/* End k1 */ |
fprintf(fichtm,"</ul>"); |
fprintf(fichtm,"</ul>"); |
|
|
|
|
fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
\n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\ |
\n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\ |
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres); |
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br> \ |
|
- 95%% confidence intervals and Wald tests of the estimated parameters are in the log file.<br> \ |
|
But because parameters are usually highly correlated (a higher incidence of disability \ |
|
and a higher incidence of recovery can give very close observed transition) it might \ |
|
be very useful to look not only at linear confidence intervals estimated from the \ |
|
variances but at the covariance matrix. And instead of looking at the estimated coefficients \ |
|
(parameters) of the logistic regression, it might be more meaningful to visualize the \ |
|
covariance matrix of the one-step probabilities. \ |
|
See page 'Matrix of variance-covariance of one-step probabilities' below. \n", rfileres,rfileres); |
|
|
fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
fprintf(fichtm," - Standard deviation of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
subdirf2(fileres,"prob"),subdirf2(fileres,"prob")); |
subdirf2(fileres,"prob"),subdirf2(fileres,"prob")); |
fprintf(fichtm,"\ |
fprintf(fichtm,"\ |
- Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
- Variance-covariance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n", |
Line 4346 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
Line 4619 fprintf(fichtm," \n<ul><li><b>Graphs</b>
|
|
|
jj1=0; |
jj1=0; |
for(k1=1; k1<=m;k1++){ |
for(k1=1; k1<=m;k1++){ |
for(i1=1; i1<=ncodemax[k1];i1++){ |
/* for(i1=1; i1<=ncodemax[k1];i1++){ */ |
jj1++; |
jj1++; |
if (cptcovn > 0) { |
if (cptcovn > 0) { |
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates"); |
for (cpt=1; cpt<=cptcoveff;cpt++) |
for (cpt=1; cpt<=cptcoveff;cpt++) |
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]); |
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtabm(jj1,cpt)]); |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">"); |
} |
} |
for(cpt=1; cpt<=nlstate;cpt++) { |
for(cpt=1; cpt<=nlstate;cpt++) { |
fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \ |
fprintf(fichtm,"<br>- Observed (cross-sectional) and period (incidence based) \ |
prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.png <br>\ |
prevalence (with 95%% confidence interval) in state (%d): %s%d_%d.svg <br>\ |
<img src=\"%s%d_%d.png\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1); |
<img src=\"%s%d_%d.svg\">",cpt,subdirf2(optionfilefiname,"v"),cpt,jj1,subdirf2(optionfilefiname,"v"),cpt,jj1); |
} |
} |
fprintf(fichtm,"\n<br>- Total life expectancy by age and \ |
fprintf(fichtm,"\n<br>- Total life expectancy by age and \ |
health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \ |
health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \ |
true period expectancies (those weighted with period prevalences are also\ |
true period expectancies (those weighted with period prevalences are also\ |
drawn in addition to the population based expectancies computed using\ |
drawn in addition to the population based expectancies computed using\ |
observed and cahotic prevalences: %s%d.png<br>\ |
observed and cahotic prevalences: %s%d.svg<br>\ |
<img src=\"%s%d.png\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1); |
<img src=\"%s%d.svg\">",subdirf2(optionfilefiname,"e"),jj1,subdirf2(optionfilefiname,"e"),jj1); |
} /* end i1 */ |
/* } /\* end i1 *\/ */ |
}/* End k1 */ |
}/* End k1 */ |
fprintf(fichtm,"</ul>"); |
fprintf(fichtm,"</ul>"); |
fflush(fichtm); |
fflush(fichtm); |
Line 4393 void printinggnuplot(char fileres[], cha
|
Line 4666 void printinggnuplot(char fileres[], cha
|
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n"); |
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'vpl' files\n"); |
for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
for (k1=1; k1<= m ; k1 ++) { /* plot [100000000000000000000:-100000000000000000000] "mysbiaspar/vplrmysbiaspar.txt to check */ |
fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"v"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s%d_%d.svg\" \n",subdirf2(optionfilefiname,"v"),cpt,k1); |
fprintf(ficgp,"\n#set out \"v%s%d_%d.png\" \n",optionfilefiname,cpt,k1); |
fprintf(ficgp,"\n#set out \"v%s%d_%d.svg\" \n",optionfilefiname,cpt,k1); |
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 svg size 640, 480\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 ++) { |
Line 4420 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 4693 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
/*2 eme*/ |
/*2 eme*/ |
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); |
fprintf(ficgp,"\n# 2nd: Total life expectancy with CI: 't' files\n"); |
for (k1=1; k1<= m ; k1 ++) { |
for (k1=1; k1<= m ; k1 ++) { |
fprintf(ficgp,"\nset out \"%s%d.png\" \n",subdirf2(optionfilefiname,"e"),k1); |
fprintf(ficgp,"\nset out \"%s%d.svg\" \n",subdirf2(optionfilefiname,"e"),k1); |
fprintf(ficgp,"set ylabel \"Years\" \nset ter png small size 320, 240\nplot [%.f:%.f] ",ageminpar,fage); |
fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage); |
|
|
for (i=1; i<= nlstate+1 ; i ++) { |
for (i=1; i<= nlstate+1 ; i ++) { |
k=2*i; |
k=2*i; |
Line 4454 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 4727 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
for (cpt=1; cpt<= nlstate ; cpt ++) { |
for (cpt=1; cpt<= nlstate ; cpt ++) { |
/* k=2+nlstate*(2*cpt-2); */ |
/* k=2+nlstate*(2*cpt-2); */ |
k=2+(nlstate+1)*(cpt-1); |
k=2+(nlstate+1)*(cpt-1); |
fprintf(ficgp,"\nset out \"%s%d%d.png\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s%d%d.svg\" \n",subdirf2(optionfilefiname,"exp"),cpt,k1); |
fprintf(ficgp,"set ter png small size 320, 240\n\ |
fprintf(ficgp,"set ter svg size 640, 480\n\ |
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileres,"e"),k1-1,k1-1,k,cpt); |
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileres,"e"),k1-1,k1-1,k,cpt); |
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); |
/*fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2*$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1); |
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); |
for (i=1; i<= nlstate*2 ; i ++) fprintf(ficgp,"\%%lf (\%%lf) "); |
Line 4479 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
Line 4752 plot [%.f:%.f] \"%s\" every :::%d::%d u
|
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */ |
k=3; |
k=3; |
fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt); |
fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, cov=%d state=%d",k1, cpt); |
fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"p"),cpt,k1); |
fprintf(ficgp,"\nset out \"%s%d_%d.svg\" \n",subdirf2(optionfilefiname,"p"),cpt,k1); |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\ |
set ter png small size 320, 240\n\ |
set ter svg size 640, 480\n\ |
unset log y\n\ |
unset log y\n\ |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
plot [%.f:%.f] ", ageminpar, agemaxpar); |
for (i=1; i<= nlstate ; i ++){ |
for (i=1; i<= nlstate ; i ++){ |
Line 4500 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
Line 4773 plot [%.f:%.f] ", ageminpar, agemaxpar)
|
} /* end covariate */ |
} /* end covariate */ |
|
|
/* proba elementaires */ |
/* proba elementaires */ |
|
fprintf(ficgp,"\n##############\n#MLE estimated parameters\n#############\n"); |
for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
|
fprintf(ficgp,"# initial state %d\n",i); |
for(k=1; k <=(nlstate+ndeath); k++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
if (k != i) { |
if (k != i) { |
|
fprintf(ficgp,"# current state %d\n",k); |
for(j=1; j <=ncovmodel; j++){ |
for(j=1; j <=ncovmodel; j++){ |
fprintf(ficgp,"p%d=%f ",jk,p[jk]); |
fprintf(ficgp,"p%d=%f; ",jk,p[jk]); |
jk++; |
jk++; |
fprintf(ficgp,"\n"); |
|
} |
} |
|
fprintf(ficgp,"\n"); |
} |
} |
} |
} |
} |
} |
|
fprintf(ficgp,"##############\n#\n"); |
|
|
/*goto avoid;*/ |
/*goto avoid;*/ |
|
fprintf(ficgp,"\n##############\n#Graphics of of probabilities or incidences\n#############\n"); |
|
fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n"); |
|
fprintf(ficgp,"# logi(p12/p11)=p1 +p2*age +p3*age*age+ p4*V1+ p5*V1*age\n"); |
|
fprintf(ficgp,"# logi(p13/p11)=a13+b13*age+c13age*age+d13*V1+e13*V1*age\n"); |
|
fprintf(ficgp,"# logi(p13/p11)=p6 +p7*age +p8*age*age+ p9*V1+ p10*V1*age\n"); |
|
fprintf(ficgp,"# p12+p13+p14+p11=1=p11(1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); |
|
fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age)+...)\n"); |
|
fprintf(ficgp,"# p11=1/(1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); |
|
fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age)+...)\n"); |
|
fprintf(ficgp,"# p12=exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)/\n"); |
|
fprintf(ficgp,"# (1+exp(a12+b12*age+c12age*age+d12*V1+e12*V1*age)\n"); |
|
fprintf(ficgp,"# +exp(a13+b13*age+c13age*age+d13*V1+e13*V1*age))\n"); |
|
fprintf(ficgp,"# +exp(a14+b14*age+c14age*age+d14*V1+e14*V1*age)+...)\n"); |
|
fprintf(ficgp,"#\n"); |
for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/ |
for(ng=1; ng<=2;ng++){ /* Number of graphics: first is probabilities second is incidence per year*/ |
|
fprintf(ficgp,"# ng=%d\n",ng); |
|
fprintf(ficgp,"# jk=1 to 2^%d=%d\n",cptcoveff,m); |
for(jk=1; jk <=m; jk++) { |
for(jk=1; jk <=m; jk++) { |
fprintf(ficgp,"\nset out \"%s%d_%d.png\" \n",subdirf2(optionfilefiname,"pe"),jk,ng); |
fprintf(ficgp,"# jk=%d\n",jk); |
|
fprintf(ficgp,"\nset out \"%s%d_%d.svg\" \n",subdirf2(optionfilefiname,"pe"),jk,ng); |
if (ng==2) |
if (ng==2) |
fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n"); |
fprintf(ficgp,"\nset ylabel \"Quasi-incidence per year\"\n"); |
else |
else |
fprintf(ficgp,"\nset title \"Probability\"\n"); |
fprintf(ficgp,"\nset title \"Probability\"\n"); |
fprintf(ficgp,"\nset ter png small size 320, 240\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar); |
fprintf(ficgp,"\nset ter svg size 640, 480\nset log y\nplot [%.f:%.f] ",ageminpar,agemaxpar); |
i=1; |
i=1; |
for(k2=1; k2<=nlstate; k2++) { |
for(k2=1; k2<=nlstate; k2++) { |
k3=i; |
k3=i; |
for(k=1; k<=(nlstate+ndeath); k++) { |
for(k=1; k<=(nlstate+ndeath); k++) { |
if (k != k2){ |
if (k != k2){ |
if(ng==2) |
if(ng==2) |
fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); |
if(nagesqr==0) |
|
fprintf(ficgp," %f*exp(p%d+p%d*x",YEARM/stepm,i,i+1); |
|
else /* nagesqr =1 */ |
|
fprintf(ficgp," %f*exp(p%d+p%d*x+p%d*x*x",YEARM/stepm,i,i+1,i+1+nagesqr); |
else |
else |
fprintf(ficgp," exp(p%d+p%d*x",i,i+1); |
if(nagesqr==0) |
|
fprintf(ficgp," exp(p%d+p%d*x",i,i+1); |
|
else /* nagesqr =1 */ |
|
fprintf(ficgp," exp(p%d+p%d*x+p%d*x*x",i,i+1,i+1+nagesqr); |
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-nagesqr; j++) { |
/* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /\* Bug valgrind *\/ */ |
/* printf("Tage[%d]=%d, j=%d\n", ij, Tage[ij], j); */ |
/* /\*fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);*\/ */ |
if(ij <=cptcovage) { /* Bug valgrind */ |
/* ij++; */ |
if((j-2)==Tage[ij]) { /* Bug valgrind */ |
/* } */ |
fprintf(ficgp,"+p%d*%d*x",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); |
/* else */ |
ij++; |
fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]); |
} |
|
} |
|
else |
|
fprintf(ficgp,"+p%d*%d",i+j+nagesqr-1,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
} |
} |
fprintf(ficgp,")/(1"); |
fprintf(ficgp,")/(1"); |
|
|
for(k1=1; k1 <=nlstate; k1++){ |
for(k1=1; k1 <=nlstate; k1++){ |
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); |
if(nagesqr==0) |
|
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1); |
|
else /* nagesqr =1 */ |
|
fprintf(ficgp,"+exp(p%d+p%d*x+p%d*x*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1,k3+(k1-1)*ncovmodel+1+nagesqr); |
|
|
ij=1; |
ij=1; |
for(j=3; j <=ncovmodel; j++){ |
for(j=3; j <=ncovmodel-nagesqr; j++){ |
/* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { */ |
if(ij <=cptcovage) { /* Bug valgrind */ |
/* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); */ |
if((j-2)==Tage[ij]) { /* Bug valgrind */ |
/* ij++; */ |
fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,Tvar[j-2])]); |
/* } */ |
ij++; |
/* else */ |
} |
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]); |
} |
|
else |
|
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2+nagesqr,nbcode[Tvar[j-2]][codtabm(jk,j-2)]); |
} |
} |
fprintf(ficgp,")"); |
fprintf(ficgp,")"); |
} |
} |
Line 4676 void prevforecast(char fileres[], double
|
Line 4986 void prevforecast(char fileres[], double
|
k=k+1; |
k=k+1; |
fprintf(ficresf,"\n#******"); |
fprintf(ficresf,"\n#******"); |
for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficresf," V%d=%d, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
} |
} |
fprintf(ficresf,"******\n"); |
fprintf(ficresf,"******\n"); |
fprintf(ficresf,"# Covariate valuofcovar yearproj age"); |
fprintf(ficresf,"# Covariate valuofcovar yearproj age"); |
Line 4700 void prevforecast(char fileres[], double
|
Line 5010 void prevforecast(char fileres[], double
|
if (h*hstepm/YEARM*stepm ==yearp) { |
if (h*hstepm/YEARM*stepm ==yearp) { |
fprintf(ficresf,"\n"); |
fprintf(ficresf,"\n"); |
for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); |
fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+h*hstepm/YEARM*stepm); |
} |
} |
for(j=1; j<=nlstate+ndeath;j++) { |
for(j=1; j<=nlstate+ndeath;j++) { |
Line 4798 void populforecast(char fileres[], doubl
|
Line 5108 void populforecast(char fileres[], doubl
|
k=k+1; |
k=k+1; |
fprintf(ficrespop,"\n#******"); |
fprintf(ficrespop,"\n#******"); |
for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficrespop," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
} |
} |
fprintf(ficrespop,"******\n"); |
fprintf(ficrespop,"******\n"); |
fprintf(ficrespop,"# Age"); |
fprintf(ficrespop,"# Age"); |
Line 5113 void printinghtmlmort(char fileres[], ch
|
Line 5423 void printinghtmlmort(char fileres[], ch
|
fprintf(fichtm," mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp); |
fprintf(fichtm," mu(age) =%lf*exp(%lf*(age-%d)) per year<br><br>",p[1],p[2],agegomp); |
for (i=1;i<=2;i++) |
for (i=1;i<=2;i++) |
fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
fprintf(fichtm,"<br><br><img src=\"graphmort.png\">"); |
fprintf(fichtm,"<br><br><img src=\"graphmort.svg\">"); |
fprintf(fichtm,"</ul>"); |
fprintf(fichtm,"</ul>"); |
|
|
fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>"); |
fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>"); |
Line 5142 void printinggnuplotmort(char fileres[],
|
Line 5452 void printinggnuplotmort(char fileres[],
|
|
|
strcpy(dirfileres,optionfilefiname); |
strcpy(dirfileres,optionfilefiname); |
strcpy(optfileres,"vpl"); |
strcpy(optfileres,"vpl"); |
fprintf(ficgp,"set out \"graphmort.png\"\n "); |
fprintf(ficgp,"set out \"graphmort.svg\"\n "); |
fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); |
fprintf(ficgp,"set xlabel \"Age\"\n set ylabel \"Force of mortality (per year)\" \n "); |
fprintf(ficgp, "set ter png small size 320, 240\n set log y\n"); |
fprintf(ficgp, "set ter svg size 640, 480\n set log y\n"); |
/* fprintf(ficgp, "set size 0.65,0.65\n"); */ |
/* fprintf(ficgp, "set size 0.65,0.65\n"); */ |
fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp); |
fprintf(ficgp,"plot [%d:100] %lf*exp(%lf*(x-%d))",agegomp,p[1],p[2],agegomp); |
|
|
Line 5165 int readdata(char datafile[], int firsto
|
Line 5475 int readdata(char datafile[], int firsto
|
|
|
|
|
if((fic=fopen(datafile,"r"))==NULL) { |
if((fic=fopen(datafile,"r"))==NULL) { |
printf("Problem while opening datafile: %s\n", datafile);return 1; |
printf("Problem while opening datafile: %s\n", datafile);fflush(stdout); |
fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);return 1; |
fprintf(ficlog,"Problem while opening datafile: %s\n", datafile);fflush(ficlog);return 1; |
} |
} |
|
|
i=1; |
i=1; |
Line 5343 void removespace(char *str) {
|
Line 5653 void removespace(char *str) {
|
} |
} |
|
|
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: |
* Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age |
* Model V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age |
* - cptcovt total number of covariates of the model nbocc(+)+1 = 8 |
* - nagesqr = 1 if age*age in the model, otherwise 0. |
* - cptcovn or number of covariates k of the models excluding age*products =6 |
* - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age |
|
* - cptcovn or number of covariates k of the models excluding age*products =6 and age*age |
* - cptcovage number of covariates with age*products =2 |
* - cptcovage number of covariates with age*products =2 |
* - cptcovs number of simple covariates |
* - cptcovs number of simple covariates |
* - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10 |
* - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10 |
Line 5360 int decodemodel ( char model[], int last
|
Line 5671 int decodemodel ( char model[], int last
|
int 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]; |
|
char *strpt; |
|
|
/*removespace(model);*/ |
/*removespace(model);*/ |
if (strlen(model) >1){ /* If there is at least 1 covariate */ |
if (strlen(model) >1){ /* If there is at least 1 covariate */ |
j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0; |
j=0, j1=0, k1=0, k2=-1, ks=0, cptcovn=0; |
j=nbocc(model,'+'); /**< j=Number of '+' */ |
|
j1=nbocc(model,'*'); /**< j1=Number of '*' */ |
|
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*/ |
|
/* including age products which are counted in cptcovage. |
|
* 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 */ |
|
cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */ |
|
strcpy(modelsav,model); |
|
if (strstr(model,"AGE") !=0){ |
if (strstr(model,"AGE") !=0){ |
printf("Error. AGE must be in lower case 'age' model=%s ",model); |
printf("Error. AGE must be in lower case 'age' model=1+age+%s. ",model); |
fprintf(ficlog,"Error. AGE must be in lower case model=%s ",model);fflush(ficlog); |
fprintf(ficlog,"Error. AGE must be in lower case model=1+age+%s. ",model);fflush(ficlog); |
return 1; |
return 1; |
} |
} |
if (strstr(model,"v") !=0){ |
if (strstr(model,"v") !=0){ |
Line 5383 int decodemodel ( char model[], int last
|
Line 5686 int decodemodel ( char model[], int last
|
fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog); |
fprintf(ficlog,"Error. 'v' must be in upper case model=%s ",model);fflush(ficlog); |
return 1; |
return 1; |
} |
} |
|
strcpy(modelsav,model); |
/* Design |
if ((strpt=strstr(model,"age*age")) !=0){ |
* V1 V2 V3 V4 V5 V6 V7 V8 V9 Weight |
printf(" strpt=%s, model=%s\n",strpt, model); |
* < ncovcol=8 > |
if(strpt != model){ |
* Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 |
printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ |
* k= 1 2 3 4 5 6 7 8 |
'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ |
* cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8 |
corresponding column of parameters.\n",model); |
* covar[k,i], value of kth covariate if not including age for individual i: |
fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \ |
* covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8) |
'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \ |
* Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[4]=3 Tvar[8]=8 |
corresponding column of parameters.\n",model); fflush(ficlog); |
* if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and |
return 1; |
* Tage[++cptcovage]=k |
} |
* if products, new covar are created after ncovcol with k1 |
|
* Tvar[k]=ncovcol+k1; # of the kth covariate product: Tvar[5]=ncovcol+1=10 Tvar[6]=ncovcol+1=11 |
|
* Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product |
|
* Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8 |
|
* Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2]; |
|
* Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted |
|
* V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 |
|
* < ncovcol=8 > |
|
* Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2 |
|
* k= 1 2 3 4 5 6 7 8 9 10 11 12 |
|
* Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8 |
|
* p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} |
|
* p Tprod[1]@2={ 6, 5} |
|
*p Tvard[1][1]@4= {7, 8, 5, 6} |
|
* covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8 |
|
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
|
*How to reorganize? |
|
* Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age |
|
* Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} |
|
* {2, 1, 4, 8, 5, 6, 3, 7} |
|
* Struct [] |
|
*/ |
|
|
|
/* This loop fills the array Tvar from the string 'model'.*/ |
nagesqr=1; |
/* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */ |
if (strstr(model,"+age*age") !=0) |
/* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 */ |
substrchaine(modelsav, model, "+age*age"); |
/* k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */ |
else if (strstr(model,"age*age+") !=0) |
/* k=3 V4 Tvar[k=3]= 4 (from V4) */ |
substrchaine(modelsav, model, "age*age+"); |
/* k=2 V1 Tvar[k=2]= 1 (from V1) */ |
else |
/* k=1 Tvar[1]=2 (from V2) */ |
substrchaine(modelsav, model, "age*age"); |
/* k=5 Tvar[5] */ |
}else |
/* for (k=1; k<=cptcovn;k++) { */ |
nagesqr=0; |
/* cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]; */ |
if (strlen(modelsav) >1){ |
/* } */ |
j=nbocc(modelsav,'+'); /**< j=Number of '+' */ |
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */ |
j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */ |
/* |
cptcovs=j+1-j1; /**< Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =2 */ |
* Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */ |
cptcovt= j+1; /* Number of total covariates in the model, not including |
for(k=cptcovt; k>=1;k--) /**< Number of covariates */ |
* cst, age and age*age |
|
* V1+V1*age+ V3 + V3*V4+age*age=> 4*/ |
|
/* including age products which are counted in cptcovage. |
|
* 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 */ |
|
cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1 */ |
|
|
|
|
|
/* Design |
|
* V1 V2 V3 V4 V5 V6 V7 V8 V9 Weight |
|
* < ncovcol=8 > |
|
* Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 |
|
* k= 1 2 3 4 5 6 7 8 |
|
* cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8 |
|
* covar[k,i], value of kth covariate if not including age for individual i: |
|
* covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8) |
|
* Tvar[k] # of the kth covariate: Tvar[1]=2 Tvar[4]=3 Tvar[8]=8 |
|
* if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and |
|
* Tage[++cptcovage]=k |
|
* if products, new covar are created after ncovcol with k1 |
|
* Tvar[k]=ncovcol+k1; # of the kth covariate product: Tvar[5]=ncovcol+1=10 Tvar[6]=ncovcol+1=11 |
|
* Tprod[k1]=k; Tprod[1]=5 Tprod[2]= 6; gives the position of the k1th product |
|
* Tvard[k1][1]=m Tvard[k1][2]=m; Tvard[1][1]=5 (V5) Tvard[1][2]=6 Tvard[2][1]=7 (V7) Tvard[2][2]=8 |
|
* Tvar[cptcovn+k2]=Tvard[k1][1];Tvar[cptcovn+k2+1]=Tvard[k1][2]; |
|
* Tvar[8+1]=5;Tvar[8+2]=6;Tvar[8+3]=7;Tvar[8+4]=8 inverted |
|
* V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 |
|
* < ncovcol=8 > |
|
* Model V2 + V1 + V3*age + V3 + V5*V6 + V7*V8 + V8*age + V8 d1 d1 d2 d2 |
|
* k= 1 2 3 4 5 6 7 8 9 10 11 12 |
|
* Tvar[k]= 2 1 3 3 10 11 8 8 5 6 7 8 |
|
* p Tvar[1]@12={2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} |
|
* p Tprod[1]@2={ 6, 5} |
|
*p Tvard[1][1]@4= {7, 8, 5, 6} |
|
* covar[k][i]= V2 V1 ? V3 V5*V6? V7*V8? ? V8 |
|
* cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
|
*How to reorganize? |
|
* Model V1 + V2 + V3 + V8 + V5*V6 + V7*V8 + V3*age + V8*age |
|
* Tvars {2, 1, 3, 3, 11, 10, 8, 8, 7, 8, 5, 6} |
|
* {2, 1, 4, 8, 5, 6, 3, 7} |
|
* Struct [] |
|
*/ |
|
|
|
/* This loop fills the array Tvar from the string 'model'.*/ |
|
/* j is the number of + signs in the model V1+V2+V3 j=2 i=3 to 1 */ |
|
/* modelsav=V2+V1+V4+age*V3 strb=age*V3 stra=V2+V1+V4 */ |
|
/* k=4 (age*V3) Tvar[k=4]= 3 (from V3) Tage[cptcovage=1]=4 */ |
|
/* k=3 V4 Tvar[k=3]= 4 (from V4) */ |
|
/* k=2 V1 Tvar[k=2]= 1 (from V1) */ |
|
/* k=1 Tvar[1]=2 (from V2) */ |
|
/* k=5 Tvar[5] */ |
|
/* for (k=1; k<=cptcovn;k++) { */ |
|
/* cov[2+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */ |
|
/* } */ |
|
/* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,Tvar[Tage[k])]]*cov[2]; */ |
|
/* |
|
* 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 */ |
Tvar[k]=0; |
Tvar[k]=0; |
cptcovage=0; |
cptcovage=0; |
for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */ |
for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */ |
cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' |
cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' |
modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ |
modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ |
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ |
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */ |
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ |
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/ |
/*scanf("%d",i);*/ |
/*scanf("%d",i);*/ |
if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */ |
if (strchr(strb,'*')) { /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */ |
cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ |
cutl(strc,strd,strb,'*'); /**< strd*strc Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */ |
if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ |
if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */ |
/* covar is not filled and then is empty */ |
/* covar is not filled and then is empty */ |
cptcovprod--; |
cptcovprod--; |
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; V1+V1*age Tvar[2]=1 */ |
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--; |
cutl(stre,strb,strc,'V'); |
cutl(stre,strb,strc,'V'); |
Tvar[k]=atoi(stre); |
Tvar[k]=atoi(stre); |
cptcovage++; |
cptcovage++; |
Tage[cptcovage]=k; |
Tage[cptcovage]=k; |
} else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ |
} else { /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2 strb=V3*V2*/ |
/* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ |
/* loops on k1=1 (V3*V2) and k1=2 V4*V3 */ |
|
cptcovn++; |
|
cptcovprodnoage++;k1++; |
|
cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ |
|
Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but |
|
because this model-covariate is a construction we invent a new column |
|
ncovcol + k1 |
|
If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 |
|
Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ |
|
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ |
|
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ |
|
Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ |
|
Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ |
|
k2=k2+2; |
|
Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */ |
|
Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */ |
|
for (i=1; i<=lastobs;i++){ |
|
/* Computes the new covariate which is a product of |
|
covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ |
|
covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; |
|
} |
|
} /* End age is not in the model */ |
|
} /* End if model includes a product */ |
|
else { /* no more sum */ |
|
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ |
|
/* scanf("%d",i);*/ |
|
cutl(strd,strc,strb,'V'); |
|
ks++; /**< Number of simple covariates */ |
cptcovn++; |
cptcovn++; |
cptcovprodnoage++;k1++; |
Tvar[k]=atoi(strd); |
cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/ |
} |
Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but |
strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ |
because this model-covariate is a construction we invent a new column |
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); |
ncovcol + k1 |
scanf("%d",i);*/ |
If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2 |
} /* end of loop + on total covariates */ |
Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */ |
} /* end if strlen(modelsave == 0) age*age might exist */ |
cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */ |
} /* end if strlen(model == 0) */ |
Tprod[k1]=k; /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2 */ |
|
Tvard[k1][1] =atoi(strc); /* m 1 for V1*/ |
|
Tvard[k1][2] =atoi(stre); /* n 4 for V4*/ |
|
k2=k2+2; |
|
Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */ |
|
Tvar[cptcovt+k2+1]=Tvard[k1][2]; /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */ |
|
for (i=1; i<=lastobs;i++){ |
|
/* Computes the new covariate which is a product of |
|
covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */ |
|
covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i]; |
|
} |
|
} /* End age is not in the model */ |
|
} /* End if model includes a product */ |
|
else { /* no more sum */ |
|
/*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/ |
|
/* scanf("%d",i);*/ |
|
cutl(strd,strc,strb,'V'); |
|
ks++; /**< Number of simple covariates */ |
|
cptcovn++; |
|
Tvar[k]=atoi(strd); |
|
} |
|
strcpy(modelsav,stra); /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ |
|
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav); |
|
scanf("%d",i);*/ |
|
} /* end of loop + */ |
|
} /* end model */ |
|
|
|
/*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products. |
/*The number n of Vn is stored in Tvar. cptcovage =number of age covariate. Tage gives the position of age. cptcovprod= number of products. |
If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ |
If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ |
Line 5646 BOOL IsWow64()
|
Line 5985 BOOL IsWow64()
|
} |
} |
#endif |
#endif |
|
|
void syscompilerinfo() |
void syscompilerinfo(int logged) |
{ |
{ |
/* #include "syscompilerinfo.h"*/ |
/* #include "syscompilerinfo.h"*/ |
/* command line Intel compiler 32bit windows, XP compatible:*/ |
/* command line Intel compiler 32bit windows, XP compatible:*/ |
Line 5654 void syscompilerinfo()
|
Line 5993 void syscompilerinfo()
|
/Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D |
/Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D |
"_CONSOLE" /D "_LIB" /D "_USING_V110_SDK71_" /D "_UNICODE" /D |
"_CONSOLE" /D "_LIB" /D "_USING_V110_SDK71_" /D "_UNICODE" /D |
"UNICODE" /Qipo /Zc:forScope /Gd /Oi /MT /Fa"Release\" /EHsc /nologo |
"UNICODE" /Qipo /Zc:forScope /Gd /Oi /MT /Fa"Release\" /EHsc /nologo |
/Fo"Release\" /Qprof-dir "Release\" /Fp"Release\IMaCh.pch" |
/Fo"Release\" /Qprof-dir "Release\" /Fp"Release\IMaCh.pch" |
/* 64 bits */ |
*/ |
|
/* 64 bits */ |
/* |
/* |
/GS /W3 /Gy |
/GS /W3 /Gy |
/Zc:wchar_t /Zi /O2 /Fd"x64\Release\vc120.pdb" /D "WIN32" /D "NDEBUG" |
/Zc:wchar_t /Zi /O2 /Fd"x64\Release\vc120.pdb" /D "WIN32" /D "NDEBUG" |
Line 5669 void syscompilerinfo()
|
Line 6009 void syscompilerinfo()
|
/Zc:forScope /Oi /MD /Fa"x64\Release\" /EHsc /nologo /Qparallel |
/Zc:forScope /Oi /MD /Fa"x64\Release\" /EHsc /nologo /Qparallel |
/Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\IMaCh.pch" |
/Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\IMaCh.pch" |
*/ |
*/ |
/* Link is $/ /* /OUT:"visual studio |
/* Link is */ /* /OUT:"visual studio |
2013\Projects\IMaCh\Release\IMaCh.exe" /MANIFEST /NXCOMPAT |
2013\Projects\IMaCh\Release\IMaCh.exe" /MANIFEST /NXCOMPAT |
/PDB:"visual studio |
/PDB:"visual studio |
2013\Projects\IMaCh\Release\IMaCh.pdb" /DYNAMICBASE |
2013\Projects\IMaCh\Release\IMaCh.pdb" /DYNAMICBASE |
Line 5694 void syscompilerinfo()
|
Line 6034 void syscompilerinfo()
|
int cross = CROSS; |
int cross = CROSS; |
if (cross){ |
if (cross){ |
printf("Cross-"); |
printf("Cross-"); |
fprintf(ficlog, "Cross-"); |
if(logged) fprintf(ficlog, "Cross-"); |
} |
} |
#endif |
#endif |
|
|
#include <stdint.h> |
#include <stdint.h> |
|
|
printf("Compiled with:");fprintf(ficlog,"Compiled with:"); |
printf("Compiled with:");if(logged)fprintf(ficlog,"Compiled with:"); |
#if defined(__clang__) |
#if defined(__clang__) |
printf(" Clang/LLVM");fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */ |
printf(" Clang/LLVM");if(logged)fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */ |
#endif |
#endif |
#if defined(__ICC) || defined(__INTEL_COMPILER) |
#if defined(__ICC) || defined(__INTEL_COMPILER) |
printf(" Intel ICC/ICPC");fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */ |
printf(" Intel ICC/ICPC");if(logged)fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */ |
#endif |
#endif |
#if defined(__GNUC__) || defined(__GNUG__) |
#if defined(__GNUC__) || defined(__GNUG__) |
printf(" GNU GCC/G++");fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */ |
printf(" GNU GCC/G++");if(logged)fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */ |
#endif |
#endif |
#if defined(__HP_cc) || defined(__HP_aCC) |
#if defined(__HP_cc) || defined(__HP_aCC) |
printf(" Hewlett-Packard C/aC++");fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */ |
printf(" Hewlett-Packard C/aC++");if(logged)fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */ |
#endif |
#endif |
#if defined(__IBMC__) || defined(__IBMCPP__) |
#if defined(__IBMC__) || defined(__IBMCPP__) |
printf(" IBM XL C/C++"); fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */ |
printf(" IBM XL C/C++"); if(logged) fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */ |
#endif |
#endif |
#if defined(_MSC_VER) |
#if defined(_MSC_VER) |
printf(" Microsoft Visual Studio");fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */ |
printf(" Microsoft Visual Studio");if(logged)fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */ |
#endif |
#endif |
#if defined(__PGI) |
#if defined(__PGI) |
printf(" Portland Group PGCC/PGCPP");fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */ |
printf(" Portland Group PGCC/PGCPP");if(logged) fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */ |
#endif |
#endif |
#if defined(__SUNPRO_C) || defined(__SUNPRO_CC) |
#if defined(__SUNPRO_C) || defined(__SUNPRO_CC) |
printf(" Oracle Solaris Studio");fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */ |
printf(" Oracle Solaris Studio");if(logged)fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */ |
#endif |
#endif |
printf(" for ");fprintf(ficlog," for "); |
printf(" for "); if (logged) fprintf(ficlog, " for "); |
|
|
// http://stackoverflow.com/questions/4605842/how-to-identify-platform-compiler-from-preprocessor-macros |
// 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! |
#ifdef _WIN32 // note the underscore: without it, it's not msdn official! |
// Windows (x64 and x86) |
// Windows (x64 and x86) |
printf("Windows (x64 and x86) ");fprintf(ficlog,"Windows (x64 and x86) "); |
printf("Windows (x64 and x86) ");if(logged) fprintf(ficlog,"Windows (x64 and x86) "); |
#elif __unix__ // all unices, not all compilers |
#elif __unix__ // all unices, not all compilers |
// Unix |
// Unix |
printf("Unix ");fprintf(ficlog,"Unix "); |
printf("Unix ");if(logged) fprintf(ficlog,"Unix "); |
#elif __linux__ |
#elif __linux__ |
// linux |
// linux |
printf("linux ");fprintf(ficlog,"linux "); |
printf("linux ");if(logged) fprintf(ficlog,"linux "); |
#elif __APPLE__ |
#elif __APPLE__ |
// Mac OS, not sure if this is covered by __posix__ and/or __unix__ though.. |
// Mac OS, not sure if this is covered by __posix__ and/or __unix__ though.. |
printf("Mac OS ");fprintf(ficlog,"Mac OS "); |
printf("Mac OS ");if(logged) fprintf(ficlog,"Mac OS "); |
#endif |
#endif |
|
|
/* __MINGW32__ */ |
/* __MINGW32__ */ |
Line 5753 void syscompilerinfo()
|
Line 6093 void syscompilerinfo()
|
/* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */ |
/* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */ |
|
|
#if UINTPTR_MAX == 0xffffffff |
#if UINTPTR_MAX == 0xffffffff |
printf(" 32-bit"); fprintf(ficlog," 32-bit");/* 32-bit */ |
printf(" 32-bit"); if(logged) fprintf(ficlog," 32-bit");/* 32-bit */ |
#elif UINTPTR_MAX == 0xffffffffffffffff |
#elif UINTPTR_MAX == 0xffffffffffffffff |
printf(" 64-bit"); fprintf(ficlog," 64-bit");/* 64-bit */ |
printf(" 64-bit"); if(logged) fprintf(ficlog," 64-bit");/* 64-bit */ |
#else |
#else |
printf(" wtf-bit"); fprintf(ficlog," wtf-bit");/* wtf */ |
printf(" wtf-bit"); if(logged) fprintf(ficlog," wtf-bit");/* wtf */ |
#endif |
#endif |
|
|
#if defined(__GNUC__) |
#if defined(__GNUC__) |
Line 5770 void syscompilerinfo()
|
Line 6110 void syscompilerinfo()
|
+ __GNUC_MINOR__ * 100) |
+ __GNUC_MINOR__ * 100) |
# endif |
# endif |
printf(" using GNU C version %d.\n", __GNUC_VERSION__); |
printf(" using GNU C version %d.\n", __GNUC_VERSION__); |
fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__); |
if(logged) fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__); |
|
|
if (uname(&sysInfo) != -1) { |
if (uname(&sysInfo) != -1) { |
printf("Running on: %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine); |
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); |
if(logged) fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine); |
} |
} |
else |
else |
perror("uname() error"); |
perror("uname() error"); |
//#ifndef __INTEL_COMPILER |
//#ifndef __INTEL_COMPILER |
#if !defined (__INTEL_COMPILER) && !defined(__APPLE__) |
#if !defined (__INTEL_COMPILER) && !defined(__APPLE__) |
printf("GNU libc version: %s\n", gnu_get_libc_version()); |
printf("GNU libc version: %s\n", gnu_get_libc_version()); |
fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version()); |
if(logged) fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version()); |
#endif |
#endif |
#endif |
#endif |
|
|
Line 5789 void syscompilerinfo()
|
Line 6129 void syscompilerinfo()
|
// { |
// { |
#if defined(_MSC_VER) |
#if defined(_MSC_VER) |
if (IsWow64()){ |
if (IsWow64()){ |
printf("The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n"); |
printf("\nThe 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"); |
if (logged) fprintf(ficlog, "\nThe program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n"); |
} |
} |
else{ |
else{ |
printf("The process is not running under WOW64 (i.e probably on a 64bit Windows).\n"); |
printf("\nThe program 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"); |
if (logged) fprintf(ficlog, "\nThe programm is not running under WOW64 (i.e probably on a 64bit Windows).\n"); |
} |
} |
// printf("\nPress Enter to continue..."); |
// printf("\nPress Enter to continue..."); |
// getchar(); |
// getchar(); |
Line 5838 int prevalence_limit(double *p, double *
|
Line 6178 int prevalence_limit(double *p, double *
|
//for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
//for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
k=k+1; |
k=k+1; |
/* to clean */ |
/* to clean */ |
//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]); |
//printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtabm(cptcod,cptcov)); |
fprintf(ficrespl,"\n#******"); |
fprintf(ficrespl,"\n#******"); |
printf("\n#******"); |
printf("\n#******"); |
fprintf(ficlog,"\n#******"); |
fprintf(ficlog,"\n#******"); |
for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
} |
} |
fprintf(ficrespl,"******\n"); |
fprintf(ficrespl,"******\n"); |
printf("******\n"); |
printf("******\n"); |
Line 5853 int prevalence_limit(double *p, double *
|
Line 6193 int prevalence_limit(double *p, double *
|
|
|
fprintf(ficrespl,"#Age "); |
fprintf(ficrespl,"#Age "); |
for(j=1;j<=cptcoveff;j++) { |
for(j=1;j<=cptcoveff;j++) { |
fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
} |
} |
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i); |
fprintf(ficrespl,"\n"); |
fprintf(ficrespl,"\n"); |
Line 5863 int prevalence_limit(double *p, double *
|
Line 6203 int prevalence_limit(double *p, double *
|
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); |
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); |
fprintf(ficrespl,"%.0f ",age ); |
fprintf(ficrespl,"%.0f ",age ); |
for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
fprintf(ficrespl," %.5f", prlim[i][i]); |
fprintf(ficrespl," %.5f", prlim[i][i]); |
fprintf(ficrespl,"\n"); |
fprintf(ficrespl,"\n"); |
Line 5910 int hPijx(double *p, int bage, int fage)
|
Line 6250 int hPijx(double *p, int bage, int fage)
|
for (k=1; k <= (int) pow(2,cptcoveff); k++){ |
for (k=1; k <= (int) pow(2,cptcoveff); k++){ |
fprintf(ficrespij,"\n#****** "); |
fprintf(ficrespij,"\n#****** "); |
for(j=1;j<=cptcoveff;j++) |
for(j=1;j<=cptcoveff;j++) |
fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]); |
fprintf(ficrespij,"******\n"); |
fprintf(ficrespij,"******\n"); |
|
|
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ |
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */ |
Line 5962 int main(int argc, char *argv[])
|
Line 6302 int main(int argc, char *argv[])
|
|
|
int jj, ll, li, lj, lk; |
int jj, ll, li, lj, lk; |
int numlinepar=0; /* Current linenumber of parameter file */ |
int numlinepar=0; /* Current linenumber of parameter file */ |
|
int num_filled; |
int itimes; |
int itimes; |
int NDIM=2; |
int NDIM=2; |
int vpopbased=0; |
int vpopbased=0; |
Line 5970 int main(int argc, char *argv[])
|
Line 6311 int main(int argc, char *argv[])
|
/* FILE *fichtm; *//* Html File */ |
/* FILE *fichtm; *//* Html File */ |
/* FILE *ficgp;*/ /*Gnuplot File */ |
/* FILE *ficgp;*/ /*Gnuplot File */ |
struct stat info; |
struct stat info; |
double agedeb; |
double agedeb=0.; |
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; |
|
|
double ageminpar=AGEOVERFLOW,agemin=AGEOVERFLOW, agemaxpar=-AGEOVERFLOW, agemax=-AGEOVERFLOW; |
|
|
double fret; |
double fret; |
double dum; /* Dummy variable */ |
double dum=0.; /* Dummy variable */ |
double ***p3mat; |
double ***p3mat; |
double ***mobaverage; |
double ***mobaverage; |
|
|
char line[MAXLINE]; |
char line[MAXLINE]; |
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE]; |
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE]; |
|
|
|
char model[MAXLINE], modeltemp[MAXLINE]; |
char pathr[MAXLINE], pathimach[MAXLINE]; |
char pathr[MAXLINE], pathimach[MAXLINE]; |
char *tok, *val; /* pathtot */ |
char *tok, *val; /* pathtot */ |
int firstobs=1, lastobs=10; |
int firstobs=1, lastobs=10; |
int c, h , cpt; |
int c, h , cpt, c2; |
int jl; |
int jl=0; |
int i1, j1, jk, stepsize; |
int i1, j1, jk, stepsize=0; |
|
int count=0; |
|
|
int *tab; |
int *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=0, nhstepm=0; |
int agemortsup; |
int agemortsup; |
float sumlpop=0.; |
float sumlpop=0.; |
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=0, fage=110, age, agelim, agebase; |
double bage=0, fage=110., age, agelim=0., agebase=0.; |
double ftolpl=FTOL; |
double ftolpl=FTOL; |
double **prlim; |
double **prlim; |
double ***param; /* Matrix of parameters */ |
double ***param; /* Matrix of parameters */ |
Line 6056 int main(int argc, char *argv[])
|
Line 6402 int main(int argc, char *argv[])
|
#else |
#else |
getcwd(pathcd, size); |
getcwd(pathcd, size); |
#endif |
#endif |
|
syscompilerinfo(0); |
printf("\n%s\n%s",version,fullversion); |
printf("\nIMaCh version %s, %s\n%s",version, copyright, fullversion); |
if(argc <=1){ |
if(argc <=1){ |
printf("\nEnter the parameter file name: "); |
printf("\nEnter the parameter file name: "); |
fgets(pathr,FILENAMELENGTH,stdin); |
fgets(pathr,FILENAMELENGTH,stdin); |
Line 6113 int main(int argc, char *argv[])
|
Line 6459 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 6121 int main(int argc, char *argv[])
|
Line 6467 int main(int argc, char *argv[])
|
goto end; |
goto end; |
} |
} |
fprintf(ficlog,"Log filename:%s\n",filelog); |
fprintf(ficlog,"Log filename:%s\n",filelog); |
fprintf(ficlog,"\n%s\n%s",version,fullversion); |
fprintf(ficlog,"Version %s %s",version,fullversion); |
fprintf(ficlog,"\nEnter the parameter file name: \n"); |
fprintf(ficlog,"\nEnter the parameter file name: \n"); |
fprintf(ficlog,"pathimach=%s\npathtot=%s\n\ |
fprintf(ficlog,"pathimach=%s\npathtot=%s\n\ |
path=%s \n\ |
path=%s \n\ |
Line 6129 int main(int argc, char *argv[])
|
Line 6475 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(); |
syscompilerinfo(1); |
|
|
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); |
Line 6142 int main(int argc, char *argv[])
|
Line 6488 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 6165 int main(int argc, char *argv[])
|
Line 6511 int main(int argc, char *argv[])
|
|
|
/* Reads comments: lines beginning with '#' */ |
/* Reads comments: lines beginning with '#' */ |
numlinepar=0; |
numlinepar=0; |
while((c=getc(ficpar))=='#' && c!= EOF){ |
|
ungetc(c,ficpar); |
/* First parameter line */ |
fgets(line, MAXLINE, ficpar); |
while(fgets(line, MAXLINE, ficpar)) { |
|
/* If line starts with a # it is a comment */ |
|
if (line[0] == '#') { |
|
numlinepar++; |
|
fputs(line,stdout); |
|
fputs(line,ficparo); |
|
fputs(line,ficlog); |
|
continue; |
|
}else |
|
break; |
|
} |
|
if((num_filled=sscanf(line,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", \ |
|
title, datafile, &lastobs, &firstpass,&lastpass)) !=EOF){ |
|
if (num_filled != 5) { |
|
printf("Should be 5 parameters\n"); |
|
} |
numlinepar++; |
numlinepar++; |
fputs(line,stdout); |
printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n", title, datafile, lastobs, firstpass,lastpass); |
fputs(line,ficparo); |
} |
fputs(line,ficlog); |
/* Second parameter line */ |
|
while(fgets(line, MAXLINE, ficpar)) { |
|
/* If line starts with a # it is a comment */ |
|
if (line[0] == '#') { |
|
numlinepar++; |
|
fputs(line,stdout); |
|
fputs(line,ficparo); |
|
fputs(line,ficlog); |
|
continue; |
|
}else |
|
break; |
|
} |
|
if((num_filled=sscanf(line,"ftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n", \ |
|
&ftol, &stepm, &ncovcol, &nlstate, &ndeath, &maxwav, &mle, &weightopt)) !=EOF){ |
|
if (num_filled != 8) { |
|
printf("Not 8\n"); |
|
} |
|
printf("ftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\n",ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt); |
} |
} |
ungetc(c,ficpar); |
|
|
|
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); |
/* Third parameter line */ |
numlinepar++; |
while(fgets(line, MAXLINE, ficpar)) { |
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); |
/* If line starts with a # it is a comment */ |
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); |
if (line[0] == '#') { |
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); |
numlinepar++; |
|
fputs(line,stdout); |
|
fputs(line,ficparo); |
|
fputs(line,ficlog); |
|
continue; |
|
}else |
|
break; |
|
} |
|
if((num_filled=sscanf(line,"model=1+age%[^.\n]\n", model)) !=EOF){ |
|
if (num_filled != 1) { |
|
printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line); |
|
fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line); |
|
model[0]='\0'; |
|
goto end; |
|
} |
|
else{ |
|
if (model[0]=='+'){ |
|
for(i=1; i<=strlen(model);i++) |
|
modeltemp[i-1]=model[i]; |
|
} |
|
strcpy(model,modeltemp); |
|
} |
|
/* printf(" model=1+age%s modeltemp= %s, model=%s\n",model, modeltemp, model);fflush(stdout); */ |
|
} |
|
/* 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=1+age+%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); */ |
|
/* numlinepar=numlinepar+3; /\* In general *\/ */ |
|
/* 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=1+age+%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */ |
|
if(model[strlen(model)-1]=='.') /* Suppressing leading dot in the model */ |
|
model[strlen(model)-1]='\0'; |
|
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=1+age%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=1+age%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
fflush(ficlog); |
fflush(ficlog); |
|
/* if(model[0]=='#'|| model[0]== '\0'){ */ |
|
if(model[0]=='#'){ |
|
printf("Error in 'model' line: model should start with 'model=1+age+' and end with '.' \n \ |
|
'model=1+age+.' or 'model=1+age+V1.' or 'model=1+age+age*age+V1+V1*age.' or \n \ |
|
'model=1+age+V1+V2.' or 'model=1+age+V1+V2+V1*V2.' etc. \n"); \ |
|
if(mle != -1){ |
|
printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter file.\n"); |
|
exit(1); |
|
} |
|
} |
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); |
numlinepar++; |
numlinepar++; |
|
if(line[1]=='q'){ /* This #q will quit imach (the answer is q) */ |
|
z[0]=line[1]; |
|
} |
|
/* printf("****line [1] = %c \n",line[1]); */ |
fputs(line, stdout); |
fputs(line, stdout); |
//puts(line); |
//puts(line); |
fputs(line,ficparo); |
fputs(line,ficparo); |
Line 6199 int main(int argc, char *argv[])
|
Line 6620 int main(int argc, char *argv[])
|
v1+v2*age+v2*v3 makes cptcovn = 3 |
v1+v2*age+v2*v3 makes cptcovn = 3 |
*/ |
*/ |
if (strlen(model)>1) |
if (strlen(model)>1) |
ncovmodel=2+nbocc(model,'+')+1; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7*/ |
ncovmodel=2+nbocc(model,'+')+1; /*Number of variables including intercept and age = cptcovn + intercept + age : v1+v2+v3+v2*v4+v5*age makes 5+2=7,age*age makes 3*/ |
else |
else |
ncovmodel=2; |
ncovmodel=2; /* Constant and age */ |
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ |
|
nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ |
nforce= (nlstate+ndeath-1)*nlstate; /* Number of forces ij from state i to j */ |
npar= nforce*ncovmodel; /* Number of parameters like aij*/ |
npar= nforce*ncovmodel; /* Number of parameters like aij*/ |
if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){ |
if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){ |
Line 6217 int main(int argc, char *argv[])
|
Line 6637 int main(int argc, char *argv[])
|
/*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ |
/*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/ |
if(mle==-1){ /* Print a wizard for help writing covariance matrix */ |
if(mle==-1){ /* Print a wizard for help writing covariance matrix */ |
prwizard(ncovmodel, nlstate, ndeath, model, ficparo); |
prwizard(ncovmodel, nlstate, ndeath, model, ficparo); |
printf(" You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
fprintf(ficlog," You choose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso); |
free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
fclose (ficparo); |
fclose (ficparo); |
fclose (ficlog); |
fclose (ficlog); |
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 chose 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 chose mle=-3, look at file %s for a template of covariance matrix \n",filereso); |
param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
matcov=matrix(1,npar,1,npar); |
matcov=matrix(1,npar,1,npar); |
} |
} |
Line 6252 int main(int argc, char *argv[])
|
Line 6672 int main(int argc, char *argv[])
|
if(jj==i) continue; |
if(jj==i) continue; |
j++; |
j++; |
fscanf(ficpar,"%1d%1d",&i1,&j1); |
fscanf(ficpar,"%1d%1d",&i1,&j1); |
if ((i1 != i) && (j1 != j)){ |
if ((i1 != i) || (j1 != jj)){ |
printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \ |
printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n \ |
It might be a problem of design; if ncovcol and the model are correct\n \ |
It might be a problem of design; if ncovcol and the model are correct\n \ |
run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1); |
run imach with mle=-1 to get a correct template of the parameter file.\n",numlinepar, i,j, i1, j1); |
Line 6260 run imach with mle=-1 to get a correct t
|
Line 6680 run imach with mle=-1 to get a correct t
|
} |
} |
fprintf(ficparo,"%1d%1d",i1,j1); |
fprintf(ficparo,"%1d%1d",i1,j1); |
if(mle==1) |
if(mle==1) |
printf("%1d%1d",i,j); |
printf("%1d%1d",i,jj); |
fprintf(ficlog,"%1d%1d",i,j); |
fprintf(ficlog,"%1d%1d",i,jj); |
for(k=1; k<=ncovmodel;k++){ |
for(k=1; k<=ncovmodel;k++){ |
fscanf(ficpar," %lf",¶m[i][j][k]); |
fscanf(ficpar," %lf",¶m[i][j][k]); |
if(mle==1){ |
if(mle==1){ |
Line 6342 run imach with mle=-1 to get a correct t
|
Line 6762 run imach with mle=-1 to get a correct t
|
for(i=1; i <=npar; i++) |
for(i=1; i <=npar; i++) |
for(j=1; j <=npar; j++) matcov[i][j]=0.; |
for(j=1; j <=npar; j++) matcov[i][j]=0.; |
|
|
|
/* Scans npar lines */ |
for(i=1; i <=npar; i++){ |
for(i=1; i <=npar; i++){ |
fscanf(ficpar,"%s",str); |
count=fscanf(ficpar,"%1d%1d%1d",&i1,&j1,&jk); |
|
if(count != 3){ |
|
printf("Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ |
|
This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ |
|
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); |
|
fprintf(ficlog,"Error! Error in parameter file %s at line %d after line starting with %1d%1d%1d\n\ |
|
This is probably because your covariance matrix doesn't \n contain exactly %d lines corresponding to your model line '1+age+%s'.\n\ |
|
Please run with mle=-1 to get a correct covariance matrix.\n",optionfile,numlinepar, i1,j1,jk, npar, model); |
|
exit(1); |
|
}else |
if(mle==1) |
if(mle==1) |
printf("%s",str); |
printf("%1d%1d%1d",i1,j1,jk); |
fprintf(ficlog,"%s",str); |
fprintf(ficlog,"%1d%1d%1d",i1,j1,jk); |
fprintf(ficparo,"%s",str); |
fprintf(ficparo,"%1d%1d%1d",i1,j1,jk); |
for(j=1; j <=i; j++){ |
for(j=1; j <=i; j++){ |
fscanf(ficpar," %le",&matcov[i][j]); |
fscanf(ficpar," %le",&matcov[i][j]); |
if(mle==1){ |
if(mle==1){ |
Line 6363 run imach with mle=-1 to get a correct t
|
Line 6793 run imach with mle=-1 to get a correct t
|
fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
fprintf(ficparo,"\n"); |
fprintf(ficparo,"\n"); |
} |
} |
|
/* End of read covariance matrix npar lines */ |
for(i=1; i <=npar; i++) |
for(i=1; i <=npar; i++) |
for(j=i+1;j<=npar;j++) |
for(j=i+1;j<=npar;j++) |
matcov[i][j]=matcov[j][i]; |
matcov[i][j]=matcov[j][i]; |
Line 6385 run imach with mle=-1 to get a correct t
|
Line 6816 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 6401 run imach with mle=-1 to get a correct t
|
Line 6833 run imach with mle=-1 to get a correct t
|
s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ |
s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */ |
tab=ivector(1,NCOVMAX); |
tab=ivector(1,NCOVMAX); |
ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
ncodemax=ivector(1,NCOVMAX); /* Number of code per covariate; if O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
|
ncodemaxwundef=ivector(1,NCOVMAX); /* Number of code per covariate; if - 1 O and 1 only, 2**ncov; V1+V2+V3+V4=>16 */ |
|
|
/* Reads data from file datafile */ |
/* Reads data from file datafile */ |
if (readdata(datafile, firstobs, lastobs, &imx)==1) |
if (readdata(datafile, firstobs, lastobs, &imx)==1) |
Line 6436 run imach with mle=-1 to get a correct t
|
Line 6869 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 6479 run imach with mle=-1 to get a correct t
|
Line 6915 run imach with mle=-1 to get a correct t
|
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
ncodemax[1]=1; |
ncodemax[1]=1; |
Ndum =ivector(-1,NCOVMAX); |
Ndum =ivector(-1,NCOVMAX); |
if (ncovmodel > 2) |
if (ncovmodel-nagesqr > 2 ) /* That is if covariate other than cst, age and age*age */ |
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],codtabm(100,10));*/ |
|
/* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/ |
h=0; |
h=0; |
|
|
|
|
Line 6492 run imach with mle=-1 to get a correct t
|
Line 6932 run imach with mle=-1 to get a correct t
|
|
|
m=pow(2,cptcoveff); |
m=pow(2,cptcoveff); |
|
|
for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */ |
|
for(i=1; i <=pow(2,cptcoveff-k);i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */ |
|
for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate ncodemax=2*/ |
|
for(cpt=1; cpt <=pow(2,k-1); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */ |
|
h++; |
|
if (h>m) |
|
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 6510 run imach with mle=-1 to get a correct t
|
Line 6945 run imach with mle=-1 to get a correct t
|
* 6 2 1 2 1 |
* 6 2 1 2 1 |
* 7 i=4 1 2 2 1 |
* 7 i=4 1 2 2 1 |
* 8 2 2 2 1 |
* 8 2 2 2 1 |
* 9 i=5 1 i=3 1 i=2 1 1 |
* 9 i=5 1 i=3 1 i=2 1 2 |
* 10 2 1 1 1 |
* 10 2 1 1 2 |
* 11 i=6 1 2 1 1 |
* 11 i=6 1 2 1 2 |
* 12 2 2 1 1 |
* 12 2 2 1 2 |
* 13 i=7 1 i=4 1 2 1 |
* 13 i=7 1 i=4 1 2 2 |
* 14 2 1 2 1 |
* 14 2 1 2 2 |
* 15 i=8 1 2 2 1 |
* 15 i=8 1 2 2 2 |
* 16 2 2 2 1 |
* 16 2 2 2 2 |
*/ |
*/ |
codtab[h][k]=j; |
for(h=1; h <=100 ;h++){ |
/*codtab[h][Tvar[k]]=j;*/ |
/* printf("h=%2d ", h); */ |
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]]); |
for(k=1; k <=10; k++){ |
} |
/* printf("k=%d %d ",k,codtabm(h,k)); */ |
} |
codtab[h][k]=codtabm(h,k); |
} |
} |
} |
/* printf("\n"); */ |
|
} |
|
/* for(k=1;k<=cptcoveff; k++){ /\* scans any effective covariate *\/ */ |
|
/* for(i=1; i <=pow(2,cptcoveff-k);i++){ /\* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 *\/ */ |
|
/* for(j=1; j <= ncodemax[k]; j++){ /\* For each modality of this covariate ncodemax=2*\/ */ |
|
/* for(cpt=1; cpt <=pow(2,k-1); cpt++){ /\* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 *\/ */ |
|
/* h++; */ |
|
/* if (h>m) */ |
|
/* h=1; */ |
|
/* codtab[h][k]=j; */ |
|
/* /\* codtab[12][3]=1; *\/ */ |
|
/* /\*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("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); |
/* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); |
codtab[1][2]=1;codtab[2][2]=2; */ |
codtab[1][2]=1;codtab[2][2]=2; */ |
/* for(i=1; i <=m ;i++){ |
/* for(i=1; i <=m ;i++){ */ |
for(k=1; k <=cptcovn; k++){ |
/* for(k=1; k <=cptcovn; k++){ */ |
printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); |
/* printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); */ |
} |
/* } */ |
printf("\n"); |
/* printf("\n"); */ |
} |
/* } */ |
scanf("%d",i);*/ |
/* scanf("%d",i);*/ |
|
|
free_ivector(Ndum,-1,NCOVMAX); |
free_ivector(Ndum,-1,NCOVMAX); |
|
|
|
|
|
|
/*------------ gnuplot -------------*/ |
/* Initialisation of ----------- gnuplot -------------*/ |
strcpy(optionfilegnuplot,optionfilefiname); |
strcpy(optionfilegnuplot,optionfilefiname); |
if(mle==-3) |
if(mle==-3) |
strcat(optionfilegnuplot,"-mort"); |
strcat(optionfilegnuplot,"-mort"); |
Line 6556 run imach with mle=-1 to get a correct t
|
Line 7007 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 6626 Interval (in months) between two waves:
|
Line 7079 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 6719 Interval (in months) between two waves:
|
Line 7172 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 6788 Interval (in months) between two waves:
|
Line 7241 Interval (in months) between two waves:
|
} |
} |
|
|
printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp); |
printf("iter=%d MLE=%f Eq=%lf*exp(%lf*(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp); |
for (i=1;i<=NDIM;i++) |
for (i=1;i<=NDIM;i++) { |
printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
printf("%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
|
fprintf(ficlog,"%f [%f ; %f]\n",p[i],p[i]-2*sqrt(matcov[i][i]),p[i]+2*sqrt(matcov[i][i])); |
|
} |
lsurv=vector(1,AGESUP); |
lsurv=vector(1,AGESUP); |
lpop=vector(1,AGESUP); |
lpop=vector(1,AGESUP); |
tpop=vector(1,AGESUP); |
tpop=vector(1,AGESUP); |
Line 6822 Interval (in months) between two waves:
|
Line 7276 Interval (in months) between two waves:
|
|
|
|
|
replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ |
replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ |
printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); |
if(ageminpar == AGEOVERFLOW ||agemaxpar == AGEOVERFLOW){ |
|
printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ |
|
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
|
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
|
fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ |
|
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
|
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
|
}else |
|
printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); |
printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ |
printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \ |
stepm, weightopt,\ |
stepm, weightopt,\ |
model,imx,p,matcov,agemortsup); |
model,imx,p,matcov,agemortsup); |
Line 6837 Interval (in months) between two waves:
|
Line 7298 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); |
} |
} |
|
|
/*--------- results files --------------*/ |
/*--------- results files --------------*/ |
fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); |
fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle= 0 weight=%d\nmodel=1+age+%s.\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); |
|
|
|
|
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
Line 6870 Interval (in months) between two waves:
|
Line 7332 Interval (in months) between two waves:
|
fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficlog,"%d%d ",i,k); |
fprintf(ficres,"%1d%1d ",i,k); |
fprintf(ficres,"%1d%1d ",i,k); |
for(j=1; j <=ncovmodel; j++){ |
for(j=1; j <=ncovmodel; j++){ |
printf("%lf ",p[jk]); |
printf("%12.7f ",p[jk]); |
fprintf(ficlog,"%lf ",p[jk]); |
fprintf(ficlog,"%12.7f ",p[jk]); |
fprintf(ficres,"%lf ",p[jk]); |
fprintf(ficres,"%12.7f ",p[jk]); |
jk++; |
jk++; |
} |
} |
printf("\n"); |
printf("\n"); |
Line 6886 Interval (in months) between two waves:
|
Line 7348 Interval (in months) between two waves:
|
ftolhess=ftol; /* Usually correct */ |
ftolhess=ftol; /* Usually correct */ |
hesscov(matcov, p, npar, delti, ftolhess, func); |
hesscov(matcov, p, npar, delti, ftolhess, func); |
} |
} |
|
printf("Parameters and 95%% confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W .\n But be careful that parameters are highly correlated because incidence of disability is highly correlated to incidence of recovery.\n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n"); |
|
fprintf(ficlog, "Parameters, Wald tests and Wald-based confidence intervals\n W is simply the result of the division of the parameter by the square root of covariance of the parameter.\n And Wald-based confidence intervals plus and minus 1.96 * W \n It might be better to visualize the covariance matrix. See the page 'Matrix of variance-covariance of one-step probabilities' and its graphs.\n"); |
|
for(i=1,jk=1; i <=nlstate; i++){ |
|
for(k=1; k <=(nlstate+ndeath); k++){ |
|
if (k != i) { |
|
printf("%d%d ",i,k); |
|
fprintf(ficlog,"%d%d ",i,k); |
|
for(j=1; j <=ncovmodel; j++){ |
|
printf("%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); |
|
fprintf(ficlog,"%12.7f W=%8.3f CI=[%12.7f ; %12.7f] ",p[jk], p[jk]/sqrt(matcov[jk][jk]), p[jk]-1.96*sqrt(matcov[jk][jk]),p[jk]+1.96*sqrt(matcov[jk][jk])); |
|
jk++; |
|
} |
|
printf("\n"); |
|
fprintf(ficlog,"\n"); |
|
} |
|
} |
|
} |
|
|
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); |
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n"); |
printf("# Scales (for hessian or gradient estimation)\n"); |
printf("# Scales (for hessian or gradient estimation)\n"); |
fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n"); |
fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n"); |
Line 7013 Interval (in months) between two waves:
|
Line 7493 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 7041 Interval (in months) between two waves:
|
Line 7522 Interval (in months) between two waves:
|
dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.; |
dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.; |
|
|
fscanf(ficpar,"pop_based=%d\n",&popbased); |
fscanf(ficpar,"pop_based=%d\n",&popbased); |
|
fprintf(ficlog,"pop_based=%d\n",popbased); |
fprintf(ficparo,"pop_based=%d\n",popbased); |
fprintf(ficparo,"pop_based=%d\n",popbased); |
fprintf(ficres,"pop_based=%d\n",popbased); |
fprintf(ficres,"pop_based=%d\n",popbased); |
|
|
Line 7065 Interval (in months) between two waves:
|
Line 7547 Interval (in months) between two waves:
|
/* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */ |
/* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */ |
|
|
replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ |
replace_back_to_slash(pathc,pathcd); /* Even gnuplot wants a / */ |
printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); |
if(ageminpar == AGEOVERFLOW ||agemaxpar == -AGEOVERFLOW){ |
|
printf("Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ |
|
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
|
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
|
fprintf(ficlog,"Warning! Error in gnuplot file with ageminpar %f or agemaxpar %f overflow\n\ |
|
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\ |
|
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar); |
|
}else |
|
printinggnuplot(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p); |
|
|
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ |
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,\ |
model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ |
model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,\ |
Line 7086 Interval (in months) between two waves:
|
Line 7576 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); |
prlim=matrix(1,nlstate,1,nlstate); |
Line 7125 Interval (in months) between two waves:
|
Line 7618 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 7248 Interval (in months) between two waves:
|
Line 7742 Interval (in months) between two waves:
|
|
|
|
|
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/ |
oldm=oldms;savm=savms; /* Segmentation fault */ |
oldm=oldms;savm=savms; /* ZZ Segmentation fault */ |
cptcod= 0; /* To be deleted */ |
cptcod= 0; /* To be deleted */ |
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ |
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */ |
fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); |
fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are "); |
Line 7259 Interval (in months) between two waves:
|
Line 7753 Interval (in months) between two waves:
|
fprintf(ficrest,"# Age e.. (std) "); |
fprintf(ficrest,"# Age e.. (std) "); |
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i); |
fprintf(ficrest,"\n"); |
fprintf(ficrest,"\n"); |
|
/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */ |
epj=vector(1,nlstate+1); |
epj=vector(1,nlstate+1); |
for(age=bage; age <=fage ;age++){ |
for(age=bage; age <=fage ;age++){ |
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); |
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k); /*ZZ Is it the correct prevalim */ |
if (vpopbased==1) { |
if (vpopbased==1) { |
if(mobilav ==0){ |
if(mobilav ==0){ |
for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
Line 7274 Interval (in months) between two waves:
|
Line 7768 Interval (in months) between two waves:
|
} |
} |
|
|
fprintf(ficrest," %4.0f",age); |
fprintf(ficrest," %4.0f",age); |
|
/* printf(" age %4.0f ",age); */ |
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){ |
for(i=1, epj[j]=0.;i <=nlstate;i++) { |
for(i=1, epj[j]=0.;i <=nlstate;i++) { |
epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
epj[j] += prlim[i][i]*eij[i][j][(int)age]; |
/* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ |
/*ZZZ printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]);*/ |
|
/* printf("%lf %lf ", prlim[i][i] ,eij[i][j][(int)age]); */ |
} |
} |
epj[nlstate+1] +=epj[j]; |
epj[nlstate+1] +=epj[j]; |
} |
} |
|
/* printf(" age %4.0f \n",age); */ |
|
|
for(i=1, vepp=0.;i <=nlstate;i++) |
for(i=1, vepp=0.;i <=nlstate;i++) |
for(j=1;j <=nlstate;j++) |
for(j=1;j <=nlstate;j++) |
Line 7356 Interval (in months) between two waves:
|
Line 7853 Interval (in months) between two waves:
|
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
|
|
free_ivector(ncodemax,1,NCOVMAX); |
free_ivector(ncodemax,1,NCOVMAX); |
|
free_ivector(ncodemaxwundef,1,NCOVMAX); |
free_ivector(Tvar,1,NCOVMAX); |
free_ivector(Tvar,1,NCOVMAX); |
free_ivector(Tprod,1,NCOVMAX); |
free_ivector(Tprod,1,NCOVMAX); |
free_ivector(Tvaraff,1,NCOVMAX); |
free_ivector(Tvaraff,1,NCOVMAX); |
Line 7463 Interval (in months) between two waves:
|
Line 7961 Interval (in months) between two waves:
|
} |
} |
end: |
end: |
while (z[0] != 'q') { |
while (z[0] != 'q') { |
printf("\nType q for exiting: "); |
printf("\nType q for exiting: "); fflush(stdout); |
scanf("%s",z); |
scanf("%s",z); |
} |
} |
} |
} |