|
version 1.55, 2002/07/24 17:00:55
|
version 1.66, 2003/01/28 17:23:35
|
|
Line 32
|
Line 32
|
| hPijx is the probability to be observed in state i at age x+h |
hPijx is the probability to be observed in state i at age x+h |
| conditional to the observed state i at age x. The delay 'h' can be |
conditional to the observed state i at age x. The delay 'h' can be |
| split into an exact number (nh*stepm) of unobserved intermediate |
split into an exact number (nh*stepm) of unobserved intermediate |
| states. This elementary transition (by month or quarter trimester, |
states. This elementary transition (by month, quarter, |
| semester or year) is model as a multinomial logistic. The hPx |
semester or year) is modelled as a multinomial logistic. The hPx |
| matrix is simply the matrix product of nh*stepm elementary matrices |
matrix is simply the matrix product of nh*stepm elementary matrices |
| and the contribution of each individual to the likelihood is simply |
and the contribution of each individual to the likelihood is simply |
| hPijx. |
hPijx. |
|
Line 83
|
Line 83
|
| #define ODIRSEPARATOR '\\' |
#define ODIRSEPARATOR '\\' |
| #endif |
#endif |
| |
|
| char version[80]="Imach version 0.8k, July 2002, INED-EUROREVES "; |
char version[80]="Imach version 0.91, November 2002, INED-EUROREVES "; |
| int erreur; /* Error number */ |
int erreur; /* Error number */ |
| int nvar; |
int nvar; |
| int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; |
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov; |
|
Line 99 int jmin, jmax; /* min, max spacing betw
|
Line 99 int jmin, jmax; /* min, max spacing betw
|
| int mle, weightopt; |
int mle, weightopt; |
| int **mw; /* mw[mi][i] is number of the mi wave for this individual */ |
int **mw; /* mw[mi][i] is number of the mi wave for this individual */ |
| int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ |
int **dh; /* dh[mi][i] is number of steps between mi,mi+1 for this individual */ |
| |
int **bh; /* bh[mi][i] is the bias (+ or -) for this individual if the delay between |
| |
* wave mi and wave mi+1 is not an exact multiple of stepm. */ |
| double jmean; /* Mean space between 2 waves */ |
double jmean; /* Mean space between 2 waves */ |
| double **oldm, **newm, **savm; /* Working pointers to matrices */ |
double **oldm, **newm, **savm; /* Working pointers to matrices */ |
| double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ |
double **oldms, **newms, **savms; /* Fixed working pointers to matrices */ |
|
Line 177 double ftolhess; /* Tolerance for comput
|
Line 179 double ftolhess; /* Tolerance for comput
|
| /**************** split *************************/ |
/**************** split *************************/ |
| static int split( char *path, char *dirc, char *name, char *ext, char *finame ) |
static int split( char *path, char *dirc, char *name, char *ext, char *finame ) |
| { |
{ |
| char *s; /* pointer */ |
char *ss; /* pointer */ |
| int l1, l2; /* length counters */ |
int l1, l2; /* 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 ); |
| s= strrchr( path, DIRSEPARATOR ); /* find last / */ |
ss= strrchr( path, DIRSEPARATOR ); /* find last / */ |
| if ( s == NULL ) { /* no directory, so use current */ |
if ( ss == NULL ) { /* no directory, so use current */ |
| /*if(strrchr(path, ODIRSEPARATOR )==NULL) |
/*if(strrchr(path, ODIRSEPARATOR )==NULL) |
| printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ |
printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/ |
| #if defined(__bsd__) /* get current working directory */ |
#if defined(__bsd__) /* get current working directory */ |
| extern char *getwd( ); |
extern char *getwd( ); |
| |
|
| if ( getwd( dirc ) == NULL ) { |
if ( getwd( dirc ) == NULL ) { |
| #else |
#else |
| extern char *getcwd( ); |
extern char *getcwd( ); |
| |
|
| if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { |
if ( getcwd( dirc, FILENAME_MAX ) == NULL ) { |
| #endif |
#endif |
| return( GLOCK_ERROR_GETCWD ); |
return( GLOCK_ERROR_GETCWD ); |
| } |
} |
| strcpy( name, path ); /* we've got it */ |
strcpy( name, path ); /* we've got it */ |
| } else { /* strip direcotry from path */ |
} else { /* strip direcotry from path */ |
| s++; /* after this, the filename */ |
ss++; /* after this, the filename */ |
| l2 = strlen( s ); /* length of filename */ |
l2 = strlen( ss ); /* length of filename */ |
| if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); |
if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH ); |
| strcpy( name, s ); /* 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 */ |
| } |
} |
| l1 = strlen( dirc ); /* length of directory */ |
l1 = strlen( dirc ); /* length of directory */ |
| #ifdef windows |
#ifdef windows |
| if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; } |
if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; } |
| #else |
#else |
| if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; } |
if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; } |
| #endif |
#endif |
| s = strrchr( name, '.' ); /* find last / */ |
ss = strrchr( name, '.' ); /* find last / */ |
| s++; |
ss++; |
| strcpy(ext,s); /* save extension */ |
strcpy(ext,ss); /* save extension */ |
| l1= strlen( name); |
l1= strlen( name); |
| l2= strlen( s)+1; |
l2= strlen(ss)+1; |
| strncpy( finame, name, l1-l2); |
strncpy( finame, name, l1-l2); |
| finame[l1-l2]= 0; |
finame[l1-l2]= 0; |
| return( 0 ); /* we're done */ |
return( 0 ); /* we're done */ |
| } |
} |
| |
|
| |
|
|
Line 277 void nrerror(char error_text[])
|
Line 279 void nrerror(char error_text[])
|
| { |
{ |
| fprintf(stderr,"ERREUR ...\n"); |
fprintf(stderr,"ERREUR ...\n"); |
| fprintf(stderr,"%s\n",error_text); |
fprintf(stderr,"%s\n",error_text); |
| exit(1); |
exit(EXIT_FAILURE); |
| } |
} |
| /*********************** vector *******************/ |
/*********************** vector *******************/ |
| double *vector(int nl, int nh) |
double *vector(int nl, int nh) |
|
Line 854 double **matprod2(double **out, double *
|
Line 856 double **matprod2(double **out, double *
|
| |
|
| double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij ) |
double ***hpxij(double ***po, int nhstepm, double age, int hstepm, double *x, int nlstate, int stepm, double **oldm, double **savm, int ij ) |
| { |
{ |
| /* Computes the transition matrix starting at age 'age' over 'nhstepm*hstepm*stepm' month |
/* Computes the transition matrix starting at age 'age' over |
| duration (i.e. until |
'nhstepm*hstepm*stepm' months (i.e. until |
| age (in years) age+nhstepm*stepm/12) by multiplying nhstepm*hstepm matrices. |
age (in years) age+nhstepm*hstepm*stepm/12) by multiplying |
| |
nhstepm*hstepm matrices. |
| Output is stored in matrix po[i][j][h] for h every 'hstepm' step |
Output is stored in matrix po[i][j][h] for h every 'hstepm' step |
| (typically every 2 years instead of every month which is too big). |
(typically every 2 years instead of every month which is too big |
| |
for the memory). |
| Model is determined by parameters x and covariates have to be |
Model is determined by parameters x and covariates have to be |
| included manually here. |
included manually here. |
| |
|
|
Line 914 double func( double *x)
|
Line 918 double func( double *x)
|
| double **out; |
double **out; |
| double sw; /* Sum of weights */ |
double sw; /* Sum of weights */ |
| double lli; /* Individual log likelihood */ |
double lli; /* Individual log likelihood */ |
| |
int s1, s2; |
| |
double bbh; |
| long ipmx; |
long ipmx; |
| /*extern weight */ |
/*extern weight */ |
| /* We are differentiating ll according to initial status */ |
/* We are differentiating ll according to initial status */ |
|
Line 924 double func( double *x)
|
Line 930 double func( double *x)
|
| cov[1]=1.; |
cov[1]=1.; |
| |
|
| 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 (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
if(mle==1){ |
| for(mi=1; mi<= wav[i]-1; mi++){ |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
| for (ii=1;ii<=nlstate+ndeath;ii++) |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
| for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
for(mi=1; mi<= wav[i]-1; mi++){ |
| for(d=0; d<dh[mi][i]; d++){ |
for (ii=1;ii<=nlstate+ndeath;ii++) |
| newm=savm; |
for (j=1;j<=nlstate+ndeath;j++){ |
| cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
| for (kk=1; kk<=cptcovage;kk++) { |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
| cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
} |
| } |
for(d=0; d<dh[mi][i]; d++){ |
| |
newm=savm; |
| out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
| 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
for (kk=1; kk<=cptcovage;kk++) { |
| savm=oldm; |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
| oldm=newm; |
} |
| |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
| |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
| |
savm=oldm; |
| |
oldm=newm; |
| |
} /* end mult */ |
| |
|
| |
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ |
| |
/* But now since version 0.9 we anticipate for bias and large stepm. |
| |
* If stepm is larger than one month (smallest stepm) and if the exact delay |
| |
* (in months) between two waves is not a multiple of stepm, we rounded to |
| |
* the nearest (and in case of equal distance, to the lowest) interval but now |
| |
* we keep into memory the bias bh[mi][i] and also the previous matrix product |
| |
* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the |
| |
* probability in order to take into account the bias as a fraction of the way |
| |
* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies |
| |
* -stepm/2 to stepm/2 . |
| |
* For stepm=1 the results are the same as for previous versions of Imach. |
| |
* For stepm > 1 the results are less biased than in previous versions. |
| |
*/ |
| |
s1=s[mw[mi][i]][i]; |
| |
s2=s[mw[mi+1][i]][i]; |
| |
bbh=(double)bh[mi][i]/(double)stepm; |
| |
/* bias is positive if real duration |
| |
* is higher than the multiple of stepm and negative otherwise. |
| |
*/ |
| |
/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ |
| |
lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ |
| |
/*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ |
| |
/*if(lli ==000.0)*/ |
| |
/*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ |
| |
ipmx +=1; |
| |
sw += weight[i]; |
| |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
| |
} /* end of wave */ |
| |
} /* end of individual */ |
| |
} else if(mle==2){ |
| |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
| |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
| |
for(mi=1; mi<= wav[i]-1; mi++){ |
| |
for (ii=1;ii<=nlstate+ndeath;ii++) |
| |
for (j=1;j<=nlstate+ndeath;j++){ |
| |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
| |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
| |
} |
| |
for(d=0; d<=dh[mi][i]; d++){ |
| |
newm=savm; |
| |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
| |
for (kk=1; kk<=cptcovage;kk++) { |
| |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
| |
} |
| |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
| |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
| |
savm=oldm; |
| |
oldm=newm; |
| |
} /* end mult */ |
| |
|
| |
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ |
| |
/* But now since version 0.9 we anticipate for bias and large stepm. |
| |
* If stepm is larger than one month (smallest stepm) and if the exact delay |
| |
* (in months) between two waves is not a multiple of stepm, we rounded to |
| |
* the nearest (and in case of equal distance, to the lowest) interval but now |
| |
* we keep into memory the bias bh[mi][i] and also the previous matrix product |
| |
* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the |
| |
* probability in order to take into account the bias as a fraction of the way |
| |
* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies |
| |
* -stepm/2 to stepm/2 . |
| |
* For stepm=1 the results are the same as for previous versions of Imach. |
| |
* For stepm > 1 the results are less biased than in previous versions. |
| |
*/ |
| |
s1=s[mw[mi][i]][i]; |
| |
s2=s[mw[mi+1][i]][i]; |
| |
bbh=(double)bh[mi][i]/(double)stepm; |
| |
/* bias is positive if real duration |
| |
* is higher than the multiple of stepm and negative otherwise. |
| |
*/ |
| |
lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */ |
| |
/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2]));*/ |
| |
/*lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.-+bh)*out[s1][s2])); */ /* exponential interpolation */ |
| |
/*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ |
| |
/*if(lli ==000.0)*/ |
| |
/*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ |
| |
ipmx +=1; |
| |
sw += weight[i]; |
| |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
| |
} /* end of wave */ |
| |
} /* end of individual */ |
| |
} else if(mle==3){ /* exponential inter-extrapolation */ |
| |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
| |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
| |
for(mi=1; mi<= wav[i]-1; mi++){ |
| |
for (ii=1;ii<=nlstate+ndeath;ii++) |
| |
for (j=1;j<=nlstate+ndeath;j++){ |
| |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
| |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
| |
} |
| |
for(d=0; d<dh[mi][i]; d++){ |
| |
newm=savm; |
| |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
| |
for (kk=1; kk<=cptcovage;kk++) { |
| |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
| |
} |
| |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
| |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
| |
savm=oldm; |
| |
oldm=newm; |
| |
} /* end mult */ |
| |
|
| |
/*lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ /* Original formula */ |
| |
/* But now since version 0.9 we anticipate for bias and large stepm. |
| |
* If stepm is larger than one month (smallest stepm) and if the exact delay |
| |
* (in months) between two waves is not a multiple of stepm, we rounded to |
| |
* the nearest (and in case of equal distance, to the lowest) interval but now |
| |
* we keep into memory the bias bh[mi][i] and also the previous matrix product |
| |
* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the |
| |
* probability in order to take into account the bias as a fraction of the way |
| |
* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies |
| |
* -stepm/2 to stepm/2 . |
| |
* For stepm=1 the results are the same as for previous versions of Imach. |
| |
* For stepm > 1 the results are less biased than in previous versions. |
| |
*/ |
| |
s1=s[mw[mi][i]][i]; |
| |
s2=s[mw[mi+1][i]][i]; |
| |
bbh=(double)bh[mi][i]/(double)stepm; |
| |
/* bias is positive if real duration |
| |
* is higher than the multiple of stepm and negative otherwise. |
| |
*/ |
| |
/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); */ /* linear interpolation */ |
| |
lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */ |
| |
/*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ |
| |
/*if(lli ==000.0)*/ |
| |
/*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ |
| |
ipmx +=1; |
| |
sw += weight[i]; |
| |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
| |
} /* end of wave */ |
| |
} /* end of individual */ |
| |
}else{ /* ml=4 no inter-extrapolation */ |
| |
for (i=1,ipmx=0, sw=0.; i<=imx; i++){ |
| |
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i]; |
| |
for(mi=1; mi<= wav[i]-1; mi++){ |
| |
for (ii=1;ii<=nlstate+ndeath;ii++) |
| |
for (j=1;j<=nlstate+ndeath;j++){ |
| |
oldm[ii][j]=(ii==j ? 1.0 : 0.0); |
| |
savm[ii][j]=(ii==j ? 1.0 : 0.0); |
| |
} |
| |
for(d=0; d<dh[mi][i]; d++){ |
| |
newm=savm; |
| |
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM; |
| |
for (kk=1; kk<=cptcovage;kk++) { |
| |
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2]; |
| |
} |
| |
|
| } /* end mult */ |
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, |
| |
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); |
| |
savm=oldm; |
| |
oldm=newm; |
| |
} /* end mult */ |
| |
|
| lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); |
lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */ |
| /* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/ |
ipmx +=1; |
| ipmx +=1; |
sw += weight[i]; |
| sw += weight[i]; |
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
| ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; |
} /* end of wave */ |
| } /* end of wave */ |
} /* end of individual */ |
| } /* end of individual */ |
} /* End of if */ |
| |
|
| for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; |
for(k=1,l=0.; k<=nlstate; k++) l += ll[k]; |
| /* printf("l1=%f l2=%f ",ll[1],ll[2]); */ |
/* printf("l1=%f l2=%f ",ll[1],ll[2]); */ |
| l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ |
l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */ |
|
Line 974 void mlikeli(FILE *ficres,double p[], in
|
Line 1133 void mlikeli(FILE *ficres,double p[], in
|
| powell(p,xi,npar,ftol,&iter,&fret,func); |
powell(p,xi,npar,ftol,&iter,&fret,func); |
| |
|
| printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p)); |
printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p)); |
| fprintf(ficlog,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); |
fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); |
| fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); |
fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p)); |
| |
|
| } |
} |
|
Line 1237 void lubksb(double **a, int n, int *indx
|
Line 1396 void lubksb(double **a, int n, int *indx
|
| } |
} |
| |
|
| /************ Frequencies ********************/ |
/************ Frequencies ********************/ |
| void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvar, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2) |
void freqsummary(char fileres[], int agemin, int agemax, int **s, double **agev, int nlstate, int imx, int *Tvaraff, int **nbcode, int *ncodemax,double **mint,double **anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2) |
| { /* Some frequencies */ |
{ /* Some frequencies */ |
| |
|
| int i, m, jk, k1,i1, j1, bool, z1,z2,j; |
int i, m, jk, k1,i1, j1, bool, z1,z2,j; |
|
Line 1284 void freqsummary(char fileres[], int ag
|
Line 1443 void freqsummary(char fileres[], int ag
|
| if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) |
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) |
| bool=0; |
bool=0; |
| } |
} |
| if (bool==1) { |
if (bool==1){ |
| for(m=firstpass; m<=lastpass; m++){ |
for(m=firstpass; m<=lastpass; m++){ |
| k2=anint[m][i]+(mint[m][i]/12.); |
k2=anint[m][i]+(mint[m][i]/12.); |
| if ((k2>=dateprev1) && (k2<=dateprev2)) { |
if ((k2>=dateprev1) && (k2<=dateprev2)) { |
|
Line 1482 void prevalence(int agemin, float agemax
|
Line 1641 void prevalence(int agemin, float agemax
|
| |
|
| /************* Waves Concatenation ***************/ |
/************* Waves Concatenation ***************/ |
| |
|
| void concatwav(int wav[], int **dh, int **mw, int **s, double *agedc, double **agev, int firstpass, int lastpass, int imx, int nlstate, int stepm) |
void concatwav(int wav[], int **dh, int **bh, int **mw, int **s, double *agedc, double **agev, int firstpass, int lastpass, int imx, int nlstate, int stepm) |
| { |
{ |
| /* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i. |
/* Concatenates waves: wav[i] is the number of effective (useful waves) of individual i. |
| Death is a valid wave (if date is known). |
Death is a valid wave (if date is known). |
| mw[mi][i] is the mi (mi=1 to wav[i]) effective wave of individual i |
mw[mi][i] is the mi (mi=1 to wav[i]) effective wave of individual i |
| dh[m][i] of dh[mw[mi][i][i] is the delay between two effective waves m=mw[mi][i] |
dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i] |
| and mw[mi+1][i]. dh depends on stepm. |
and mw[mi+1][i]. dh depends on stepm. |
| */ |
*/ |
| |
|
|
Line 1558 void concatwav(int wav[], int **dh, int
|
Line 1717 void concatwav(int wav[], int **dh, int
|
| jk= j/stepm; |
jk= j/stepm; |
| jl= j -jk*stepm; |
jl= j -jk*stepm; |
| ju= j -(jk+1)*stepm; |
ju= j -(jk+1)*stepm; |
| if(jl <= -ju) |
if(mle <=1){ |
| dh[mi][i]=jk; |
if(jl==0){ |
| else |
dh[mi][i]=jk; |
| dh[mi][i]=jk+1; |
bh[mi][i]=0; |
| if(dh[mi][i]==0) |
}else{ /* We want a negative bias in order to only have interpolation ie |
| dh[mi][i]=1; /* At least one step */ |
* at the price of an extra matrix product in likelihood */ |
| } |
dh[mi][i]=jk+1; |
| } |
bh[mi][i]=ju; |
| |
} |
| |
}else{ |
| |
if(jl <= -ju){ |
| |
dh[mi][i]=jk; |
| |
bh[mi][i]=jl; /* bias is positive if real duration |
| |
* is higher than the multiple of stepm and negative otherwise. |
| |
*/ |
| |
} |
| |
else{ |
| |
dh[mi][i]=jk+1; |
| |
bh[mi][i]=ju; |
| |
} |
| |
if(dh[mi][i]==0){ |
| |
dh[mi][i]=1; /* At least one step */ |
| |
bh[mi][i]=ju; /* At least one step */ |
| |
printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i); |
| |
} |
| |
if(i==298 || i==287 || i==763 ||i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm); |
| |
} |
| |
} /* end if mle */ |
| |
} /* end wave */ |
| } |
} |
| jmean=sum/k; |
jmean=sum/k; |
| printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); |
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean); |
|
Line 1575 void concatwav(int wav[], int **dh, int
|
Line 1755 void concatwav(int wav[], int **dh, int
|
| /*********** Tricode ****************************/ |
/*********** Tricode ****************************/ |
| void tricode(int *Tvar, int **nbcode, int imx) |
void tricode(int *Tvar, int **nbcode, int imx) |
| { |
{ |
| int Ndum[20],ij=1, k, j, i; |
|
| |
int Ndum[20],ij=1, k, j, i, maxncov=19; |
| int cptcode=0; |
int cptcode=0; |
| cptcoveff=0; |
cptcoveff=0; |
| |
|
| for (k=0; k<19; k++) Ndum[k]=0; |
for (k=0; k<maxncov; k++) Ndum[k]=0; |
| for (k=1; k<=7; k++) ncodemax[k]=0; |
for (k=1; k<=7; k++) ncodemax[k]=0; |
| |
|
| for (j=1; j<=(cptcovn+2*cptcovprod); j++) { |
for (j=1; j<=(cptcovn+2*cptcovprod); j++) { |
| for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { /*reads the data file to get the maximum |
| ij=(int)(covar[Tvar[j]][i]); |
modality*/ |
| Ndum[ij]++; |
ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual*/ |
| |
Ndum[ij]++; /*store the modality */ |
| /*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ |
/*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/ |
| if (ij > cptcode) cptcode=ij; |
if (ij > cptcode) cptcode=ij; /* getting the maximum of covariable |
| |
Tvar[j]. If V=sex and male is 0 and |
| |
female is 1, then cptcode=1.*/ |
| } |
} |
| |
|
| for (i=0; i<=cptcode; i++) { |
for (i=0; i<=cptcode; i++) { |
| if(Ndum[i]!=0) ncodemax[j]++; |
if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */ |
| } |
} |
| ij=1; |
|
| |
|
| |
|
| |
ij=1; |
| for (i=1; i<=ncodemax[j]; i++) { |
for (i=1; i<=ncodemax[j]; i++) { |
| for (k=0; k<=19; k++) { |
for (k=0; k<= maxncov; k++) { |
| if (Ndum[k] != 0) { |
if (Ndum[k] != 0) { |
| nbcode[Tvar[j]][ij]=k; |
nbcode[Tvar[j]][ij]=k; |
| |
/* store the modality in an array. 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; */ |
| |
|
| ij++; |
ij++; |
| } |
} |
|
Line 1608 void tricode(int *Tvar, int **nbcode, in
|
Line 1792 void tricode(int *Tvar, int **nbcode, in
|
| } |
} |
| } |
} |
| |
|
| for (k=0; k<19; k++) Ndum[k]=0; |
for (k=0; k< maxncov; k++) Ndum[k]=0; |
| |
|
| for (i=1; i<=ncovmodel-2; i++) { |
for (i=1; i<=ncovmodel-2; i++) { |
| |
/* Listing of all covariables in staement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ |
| ij=Tvar[i]; |
ij=Tvar[i]; |
| Ndum[ij]++; |
Ndum[ij]++; |
| } |
} |
| |
|
| ij=1; |
ij=1; |
| for (i=1; i<=10; i++) { |
for (i=1; i<= maxncov; i++) { |
| if((Ndum[i]!=0) && (i<=ncovcol)){ |
if((Ndum[i]!=0) && (i<=ncovcol)){ |
| Tvaraff[ij]=i; |
Tvaraff[ij]=i; /*For printing */ |
| ij++; |
ij++; |
| } |
} |
| } |
} |
| |
|
| cptcoveff=ij-1; |
cptcoveff=ij-1; /*Number of simple covariates*/ |
| } |
} |
| |
|
| /*********** Health Expectancies ****************/ |
/*********** Health Expectancies ****************/ |
|
Line 1661 void evsij(char fileres[], double ***eij
|
Line 1846 void evsij(char fileres[], double ***eij
|
| * This is mainly to measure the difference between two models: for example |
* This is mainly to measure the difference between two models: for example |
| * if stepm=24 months pijx are given only every 2 years and by summing them |
* if stepm=24 months pijx are given only every 2 years and by summing them |
| * we are calculating an estimate of the Life Expectancy assuming a linear |
* we are calculating an estimate of the Life Expectancy assuming a linear |
| * progression inbetween and thus overestimating or underestimating according |
* progression in between and thus overestimating or underestimating according |
| * to the curvature of the survival function. If, for the same date, we |
* to the curvature of the survival function. If, for the same date, we |
| * estimate the model with stepm=1 month, we can keep estepm to 24 months |
* estimate the model with stepm=1 month, we can keep estepm to 24 months |
| * to compare the new estimate of Life expectancy with the same linear |
* to compare the new estimate of Life expectancy with the same linear |
|
Line 1824 void varevsij(char optionfilefiname[], d
|
Line 2009 void varevsij(char optionfilefiname[], d
|
| char fileresprobmorprev[FILENAMELENGTH]; |
char fileresprobmorprev[FILENAMELENGTH]; |
| |
|
| if(popbased==1){ |
if(popbased==1){ |
| if(mobilav==1) |
if(mobilav!=0) |
| strcpy(digitp,"-populbased-mobilav-"); |
strcpy(digitp,"-populbased-mobilav-"); |
| else strcpy(digitp,"-populbased-nomobil-"); |
else strcpy(digitp,"-populbased-nomobil-"); |
| } |
} |
| else |
else |
| strcpy(digitp,"-stablbased-"); |
strcpy(digitp,"-stablbased-"); |
| <<<<<<< imach.c |
|
| if (mobilav!=0) { |
|
| ======= |
|
| if(mobilav!=0) |
|
| strcat(digitp,"mobilav-"); |
|
| else |
|
| strcat(digitp,"nomobil-"); |
|
| if (mobilav!=0) { |
if (mobilav!=0) { |
| >>>>>>> 1.54 |
|
| mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ |
if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ |
| fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); |
fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav); |
|
Line 1858 void varevsij(char optionfilefiname[], d
|
Line 2036 void varevsij(char optionfilefiname[], d
|
| } |
} |
| printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
printf("Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
| fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
fprintf(ficlog,"Computing total mortality p.j=w1*p1j+w2*p2j+..: result on file '%s' \n",fileresprobmorprev); |
| fprintf(ficresprobmorprev,"# probabilities of dying during a year and weighted mean w1*p1j+w2*p2j+... stand dev in()\n"); |
fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1*p1j+w2*p2j+... stand dev in()\n",estepm); |
| fprintf(ficresprobmorprev,"# Age cov=%-d",ij); |
fprintf(ficresprobmorprev,"# Age cov=%-d",ij); |
| for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
| fprintf(ficresprobmorprev," p.%-d SE",j); |
fprintf(ficresprobmorprev," p.%-d SE",j); |
|
Line 1915 void varevsij(char optionfilefiname[], d
|
Line 2093 void varevsij(char optionfilefiname[], d
|
| and note for a fixed period like k years */ |
and note for a fixed period like k years */ |
| /* We decided (b) to get a life expectancy respecting the most precise curvature of the |
/* We decided (b) to get a life expectancy respecting the most precise curvature of the |
| survival function given by stepm (the optimization length). Unfortunately it |
survival function given by stepm (the optimization length). Unfortunately it |
| means that if the survival funtion is printed only each two years of age and if |
means that if the survival funtion is printed every two years of age and if |
| you sum them up and add 1 year (area under the trapezoids) you won't get the same |
you sum them up and add 1 year (area under the trapezoids) you won't get the same |
| results. So we changed our mind and took the option of the best precision. |
results. So we changed our mind and took the option of the best precision. |
| */ |
*/ |
|
Line 1931 void varevsij(char optionfilefiname[], d
|
Line 2109 void varevsij(char optionfilefiname[], d
|
| |
|
| |
|
| for(theta=1; theta <=npar; theta++){ |
for(theta=1; theta <=npar; theta++){ |
| for(i=1; i<=npar; i++){ /* Computes gradient */ |
for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/ |
| xp[i] = x[i] + (i==theta ?delti[theta]:0); |
xp[i] = x[i] + (i==theta ?delti[theta]:0); |
| } |
} |
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
|
Line 1953 void varevsij(char optionfilefiname[], d
|
Line 2131 void varevsij(char optionfilefiname[], d
|
| gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gp[h][j] += prlim[i][i]*p3mat[i][j][h]; |
| } |
} |
| } |
} |
| /* This for computing forces of mortality (h=1)as a weighted average */ |
/* This for computing probability of death (h=1 means |
| |
computed over hstepm matrices product = hstepm*stepm months) |
| |
as a weighted average of prlim. |
| |
*/ |
| for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){ |
| for(i=1; i<= nlstate; i++) |
for(i=1; i<= nlstate; i++) |
| gpp[j] += prlim[i][i]*p3mat[i][j][1]; |
gpp[j] += prlim[i][i]*p3mat[i][j][1]; |
| } |
} |
| /* end force of mortality */ |
/* end probability of death */ |
| |
|
| for(i=1; i<=npar; i++) /* Computes gradient */ |
for(i=1; i<=npar; i++) /* Computes gradient x - delta */ |
| xp[i] = x[i] - (i==theta ?delti[theta]:0); |
xp[i] = x[i] - (i==theta ?delti[theta]:0); |
| hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij); |
| prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); |
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij); |
|
Line 1981 void varevsij(char optionfilefiname[], d
|
Line 2162 void varevsij(char optionfilefiname[], d
|
| gm[h][j] += prlim[i][i]*p3mat[i][j][h]; |
gm[h][j] += prlim[i][i]*p3mat[i][j][h]; |
| } |
} |
| } |
} |
| /* This for computing force of mortality (h=1)as a weighted average */ |
/* This for computing probability of death (h=1 means |
| |
computed over hstepm matrices product = hstepm*stepm months) |
| |
as a weighted average of prlim. |
| |
*/ |
| for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ |
| for(i=1; i<= nlstate; i++) |
for(i=1; i<= nlstate; i++) |
| gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
| } |
} |
| /* end force of mortality */ |
/* end probability of death */ |
| |
|
| for(j=1; j<= nlstate; j++) /* vareij */ |
for(j=1; j<= nlstate; j++) /* vareij */ |
| for(h=0; h<=nhstepm; h++){ |
for(h=0; h<=nhstepm; h++){ |
|
Line 2031 void varevsij(char optionfilefiname[], d
|
Line 2215 void varevsij(char optionfilefiname[], d
|
| for(i=nlstate+1;i<=nlstate+ndeath;i++) |
for(i=nlstate+1;i<=nlstate+ndeath;i++) |
| varppt[j][i]=doldmp[j][i]; |
varppt[j][i]=doldmp[j][i]; |
| /* end ppptj */ |
/* end ppptj */ |
| |
/* x centered again */ |
| hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); |
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij); |
| prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); |
prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij); |
| |
|
|
Line 2044 void varevsij(char optionfilefiname[], d
|
Line 2229 void varevsij(char optionfilefiname[], d
|
| } |
} |
| } |
} |
| |
|
| /* This for computing force of mortality (h=1)as a weighted average */ |
/* This for computing probability of death (h=1 means |
| |
computed over hstepm (estepm) matrices product = hstepm*stepm months) |
| |
as a weighted average of prlim. |
| |
*/ |
| for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ |
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){ |
| for(i=1; i<= nlstate; i++) |
for(i=1; i<= nlstate; i++) |
| gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
gmp[j] += prlim[i][i]*p3mat[i][j][1]; |
| } |
} |
| /* end force of mortality */ |
/* end probability of death */ |
| |
|
| fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); |
fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij); |
| for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
for(j=nlstate+1; j<=(nlstate+ndeath);j++){ |
|
Line 2083 void varevsij(char optionfilefiname[], d
|
Line 2271 void varevsij(char optionfilefiname[], d
|
| fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); |
| fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); |
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); |
| fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev); |
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev); |
| fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,digitp,digit); |
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", estepm,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.png\"> <br>\n", stepm,YEARM,digitp,digit); |
| */ |
*/ |
| fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit); |
fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit); |
|
Line 2104 void varevsij(char optionfilefiname[], d
|
Line 2292 void varevsij(char optionfilefiname[], d
|
| void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) |
void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int ij) |
| { |
{ |
| /* Variance of prevalence limit */ |
/* Variance of prevalence limit */ |
| /* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/ |
/* double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double **savm,double ftolpl);*/ |
| double **newm; |
double **newm; |
| double **dnewm,**doldm; |
double **dnewm,**doldm; |
| int i, j, nhstepm, hstepm; |
int i, j, nhstepm, hstepm; |
|
Line 2276 void varprob(char optionfilefiname[], do
|
Line 2464 void varprob(char optionfilefiname[], do
|
| |
|
| } |
} |
| |
|
| |
|
| cov[1]=1; |
cov[1]=1; |
| tj=cptcoveff; |
tj=cptcoveff; |
| if (cptcovn<1) {tj=1;ncodemax[1]=1;} |
if (cptcovn<1) {tj=1;ncodemax[1]=1;} |
|
Line 2284 void varprob(char optionfilefiname[], do
|
Line 2471 void varprob(char optionfilefiname[], do
|
| for(t=1; t<=tj;t++){ |
for(t=1; t<=tj;t++){ |
| for(i1=1; i1<=ncodemax[t];i1++){ |
for(i1=1; i1<=ncodemax[t];i1++){ |
| 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]][codtab[j1][z1]]); |
|
Line 2357 void varprob(char optionfilefiname[], do
|
Line 2543 void varprob(char optionfilefiname[], do
|
| |
|
| matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); |
matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov); |
| matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg); |
matprod2(doldm,dnewm,1,(nlstate)*(nlstate+ndeath),1,npar,1,(nlstate)*(nlstate+ndeath),gradg); |
| |
free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); |
| |
free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); |
| |
free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); |
| |
free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); |
| |
|
| pmij(pmmij,cov,ncovmodel,x,nlstate); |
pmij(pmmij,cov,ncovmodel,x,nlstate); |
| |
|
| k=0; |
k=0; |
|
Line 2372 void varprob(char optionfilefiname[], do
|
Line 2562 void varprob(char optionfilefiname[], do
|
| varpij[i][j][(int)age] = doldm[i][j]; |
varpij[i][j][(int)age] = doldm[i][j]; |
| |
|
| /*printf("\n%d ",(int)age); |
/*printf("\n%d ",(int)age); |
| for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ |
for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){ |
| printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); |
printf("%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); |
| fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); |
fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2*sqrt(doldm[i][i]),gm[i]+2*sqrt(doldm[i][i])); |
| }*/ |
}*/ |
| |
|
| fprintf(ficresprob,"\n%d ",(int)age); |
fprintf(ficresprob,"\n%d ",(int)age); |
| fprintf(ficresprobcov,"\n%d ",(int)age); |
fprintf(ficresprobcov,"\n%d ",(int)age); |
|
Line 2403 void varprob(char optionfilefiname[], do
|
Line 2593 void varprob(char optionfilefiname[], do
|
| |
|
| /* Confidence intervalle of pij */ |
/* Confidence intervalle of pij */ |
| /* |
/* |
| fprintf(ficgp,"\nset noparametric;unset label"); |
fprintf(ficgp,"\nset noparametric;unset label"); |
| fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\""); |
fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\""); |
| fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); |
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65"); |
| fprintf(fichtm,"\n<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname); |
fprintf(fichtm,"\n<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname); |
| fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname); |
fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname); |
| fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname); |
fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname); |
| fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob); |
fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob); |
| */ |
*/ |
| |
|
| /* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ |
/* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/ |
|
Line 2481 void varprob(char optionfilefiname[], do
|
Line 2671 void varprob(char optionfilefiname[], do
|
| } /*l1 */ |
} /*l1 */ |
| }/* k1 */ |
}/* k1 */ |
| } /* loop covariates */ |
} /* loop covariates */ |
| free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); |
|
| free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath)); |
|
| free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath)); |
|
| free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); |
|
| free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); |
|
| free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar); |
|
| } |
} |
| |
free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage); |
| |
free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage); |
| free_vector(xp,1,npar); |
free_vector(xp,1,npar); |
| fclose(ficresprob); |
fclose(ficresprob); |
| fclose(ficresprobcov); |
fclose(ficresprobcov); |
|
Line 2780 m=pow(2,cptcoveff);
|
Line 2966 m=pow(2,cptcoveff);
|
| int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){ |
int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav){ |
| |
|
| int i, cpt, cptcod; |
int i, cpt, cptcod; |
| |
int modcovmax =1; |
| int mobilavrange, mob; |
int mobilavrange, mob; |
| double age; |
double age; |
| |
|
| |
modcovmax=2*cptcoveff;/* Max number of modalities. We suppose |
| |
a covariate has 2 modalities */ |
| |
if (cptcovn<1) modcovmax=1; /* At least 1 pass */ |
| |
|
| if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ |
if(mobilav==1||mobilav ==3 ||mobilav==5 ||mobilav== 7){ |
| if(mobilav==1) mobilavrange=5; /* default */ |
if(mobilav==1) mobilavrange=5; /* default */ |
| else mobilavrange=mobilav; |
else mobilavrange=mobilav; |
| for (age=bage; age<=fage; age++) |
for (age=bage; age<=fage; age++) |
| for (i=1; i<=nlstate;i++) |
for (i=1; i<=nlstate;i++) |
| for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++) |
for (cptcod=1;cptcod<=modcovmax;cptcod++) |
| mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; |
mobaverage[(int)age][i][cptcod]=probs[(int)age][i][cptcod]; |
| /* We keep the original values on the extreme ages bage, fage and for |
/* We keep the original values on the extreme ages bage, fage and for |
| fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 |
fage+1 and bage-1 we use a 3 terms moving average; for fage+2 bage+2 |
|
Line 2796 int movingaverage(double ***probs, doubl
|
Line 2988 int movingaverage(double ***probs, doubl
|
| for (mob=3;mob <=mobilavrange;mob=mob+2){ |
for (mob=3;mob <=mobilavrange;mob=mob+2){ |
| for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ |
for (age=bage+(mob-1)/2; age<=fage-(mob-1)/2; age++){ |
| for (i=1; i<=nlstate;i++){ |
for (i=1; i<=nlstate;i++){ |
| for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ |
for (cptcod=1;cptcod<=modcovmax;cptcod++){ |
| mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; |
mobaverage[(int)age][i][cptcod] =probs[(int)age][i][cptcod]; |
| for (cpt=1;cpt<=(mob-1)/2;cpt++){ |
for (cpt=1;cpt<=(mob-1)/2;cpt++){ |
| mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; |
mobaverage[(int)age][i][cptcod] +=probs[(int)age-cpt][i][cptcod]; |
|
Line 2824 prevforecast(char fileres[], double anpr
|
Line 3016 prevforecast(char fileres[], double anpr
|
| char fileresf[FILENAMELENGTH]; |
char fileresf[FILENAMELENGTH]; |
| |
|
| agelim=AGESUP; |
agelim=AGESUP; |
| calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; |
calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM; |
| |
|
| prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); |
prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); |
| |
|
|
Line 2925 populforecast(char fileres[], double anp
|
Line 3117 populforecast(char fileres[], double anp
|
| |
|
| int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; |
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h; |
| int *popage; |
int *popage; |
| double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean; |
double calagedate, agelim, kk1, kk2; |
| double *popeffectif,*popcount; |
double *popeffectif,*popcount; |
| double ***p3mat,***tabpop,***tabpopprev; |
double ***p3mat,***tabpop,***tabpopprev; |
| double ***mobaverage; |
double ***mobaverage; |
|
Line 3086 populforecast(char fileres[], double anp
|
Line 3278 populforecast(char fileres[], double anp
|
| |
|
| int main(int argc, char *argv[]) |
int main(int argc, char *argv[]) |
| { |
{ |
| |
int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav); |
| int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod; |
int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod; |
| double agedeb, agefin,hf; |
double agedeb, agefin,hf; |
| double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; |
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20; |
|
Line 3105 int main(int argc, char *argv[])
|
Line 3297 int main(int argc, char *argv[])
|
| int c, h , cpt,l; |
int c, h , cpt,l; |
| int ju,jl, mi; |
int ju,jl, mi; |
| int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij; |
int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij; |
| int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab; |
int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,*tab; |
| int mobilav=0,popforecast=0; |
int mobilav=0,popforecast=0; |
| int hstepm, nhstepm; |
int hstepm, nhstepm; |
| double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate; |
double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate; |
|
Line 3124 int main(int argc, char *argv[])
|
Line 3316 int main(int argc, char *argv[])
|
| double *epj, vepp; |
double *epj, vepp; |
| double kk1, kk2; |
double kk1, kk2; |
| double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2; |
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2; |
| |
|
| |
|
| char *alph[]={"a","a","b","c","d","e"}, str[4]; |
char *alph[]={"a","a","b","c","d","e"}, str[4]; |
| |
|
|
Line 3135 int main(int argc, char *argv[])
|
Line 3326 int main(int argc, char *argv[])
|
| char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80]; |
char stra[80], strb[80], strc[80], strd[80],stre[80],modelsav[80]; |
| |
|
| /* long total_usecs; |
/* long total_usecs; |
| struct timeval start_time, end_time; |
struct timeval start_time, end_time; |
| |
|
| gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ |
gettimeofday(&start_time, (struct timezone*)0); */ /* at first time */ |
| getcwd(pathcd, size); |
getcwd(pathcd, size); |
| |
|
| printf("\n%s",version); |
printf("\n%s",version); |
|
Line 3154 int main(int argc, char *argv[])
|
Line 3345 int main(int argc, char *argv[])
|
| /* cutv(path,optionfile,pathtot,'\\');*/ |
/* cutv(path,optionfile,pathtot,'\\');*/ |
| |
|
| split(pathtot,path,optionfile,optionfilext,optionfilefiname); |
split(pathtot,path,optionfile,optionfilext,optionfilefiname); |
| printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); |
printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname); |
| chdir(path); |
chdir(path); |
| replace(pathc,path); |
replace(pathc,path); |
| |
|
| /*-------- arguments in the command line --------*/ |
/*-------- arguments in the command line --------*/ |
| |
|
| /* Log file */ |
/* Log file */ |
| strcat(filelog, optionfilefiname); |
strcat(filelog, optionfilefiname); |
|
Line 3206 int main(int argc, char *argv[])
|
Line 3397 int main(int argc, char *argv[])
|
| fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); |
fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model); |
| printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); |
printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); |
| fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model); |
| 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); |
| puts(line); |
puts(line); |
|
Line 3216 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3407 while((c=getc(ficpar))=='#' && c!= EOF){
|
| |
|
| |
|
| covar=matrix(0,NCOVMAX,1,n); |
covar=matrix(0,NCOVMAX,1,n); |
| cptcovn=0; |
cptcovn=0; /*Number of covariates, i.e. number of '+' in model statement*/ |
| if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; |
if (strlen(model)>1) cptcovn=nbocc(model,'+')+1; |
| |
|
| ncovmodel=2+cptcovn; |
ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */ |
| nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ |
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */ |
| |
|
| /* Read guess parameters */ |
/* Read guess parameters */ |
|
Line 3233 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3424 while((c=getc(ficpar))=='#' && c!= EOF){
|
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| |
|
| param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
param= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); |
| for(i=1; i <=nlstate; i++) |
for(i=1; i <=nlstate; i++) |
| for(j=1; j <=nlstate+ndeath-1; j++){ |
for(j=1; j <=nlstate+ndeath-1; j++){ |
| fscanf(ficpar,"%1d%1d",&i1,&j1); |
fscanf(ficpar,"%1d%1d",&i1,&j1); |
| fprintf(ficparo,"%1d%1d",i1,j1); |
fprintf(ficparo,"%1d%1d",i1,j1); |
|
Line 3257 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3448 while((c=getc(ficpar))=='#' && c!= EOF){
|
| fprintf(ficparo,"\n"); |
fprintf(ficparo,"\n"); |
| } |
} |
| |
|
| npar= (nlstate+ndeath-1)*nlstate*ncovmodel; |
npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/ |
| |
|
| p=param[1][1]; |
p=param[1][1]; |
| |
|
|
Line 3330 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3521 while((c=getc(ficpar))=='#' && c!= EOF){
|
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| |
|
| |
|
| /*-------- Rewriting paramater file ----------*/ |
/*-------- Rewriting paramater file ----------*/ |
| strcpy(rfileres,"r"); /* "Rparameterfile */ |
strcpy(rfileres,"r"); /* "Rparameterfile */ |
| strcat(rfileres,optionfilefiname); /* Parameter file first name*/ |
strcat(rfileres,optionfilefiname); /* Parameter file first name*/ |
| strcat(rfileres,"."); /* */ |
strcat(rfileres,"."); /* */ |
| strcat(rfileres,optionfilext); /* Other files have txt extension */ |
strcat(rfileres,optionfilext); /* Other files have txt extension */ |
| if((ficres =fopen(rfileres,"w"))==NULL) { |
if((ficres =fopen(rfileres,"w"))==NULL) { |
| printf("Problem writing new parameter file: %s\n", fileres);goto end; |
printf("Problem writing new parameter file: %s\n", fileres);goto end; |
| fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end; |
fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end; |
| } |
} |
| fprintf(ficres,"#%s\n",version); |
fprintf(ficres,"#%s\n",version); |
| |
|
| /*-------- data file ----------*/ |
/*-------- data file ----------*/ |
| if((fic=fopen(datafile,"r"))==NULL) { |
if((fic=fopen(datafile,"r"))==NULL) { |
| printf("Problem with datafile: %s\n", datafile);goto end; |
printf("Problem with datafile: %s\n", datafile);goto end; |
| fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end; |
fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end; |
| } |
} |
| |
|
| n= lastobs; |
n= lastobs; |
| severity = vector(1,maxwav); |
severity = vector(1,maxwav); |
| outcome=imatrix(1,maxwav+1,1,n); |
outcome=imatrix(1,maxwav+1,1,n); |
| num=ivector(1,n); |
num=ivector(1,n); |
| moisnais=vector(1,n); |
moisnais=vector(1,n); |
| annais=vector(1,n); |
annais=vector(1,n); |
| moisdc=vector(1,n); |
moisdc=vector(1,n); |
| andc=vector(1,n); |
andc=vector(1,n); |
| agedc=vector(1,n); |
agedc=vector(1,n); |
| cod=ivector(1,n); |
cod=ivector(1,n); |
| weight=vector(1,n); |
weight=vector(1,n); |
| for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */ |
for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */ |
| mint=matrix(1,maxwav,1,n); |
mint=matrix(1,maxwav,1,n); |
| anint=matrix(1,maxwav,1,n); |
anint=matrix(1,maxwav,1,n); |
| s=imatrix(1,maxwav+1,1,n); |
s=imatrix(1,maxwav+1,1,n); |
| adl=imatrix(1,maxwav+1,1,n); |
tab=ivector(1,NCOVMAX); |
| tab=ivector(1,NCOVMAX); |
ncodemax=ivector(1,8); |
| ncodemax=ivector(1,8); |
|
| |
i=1; |
| i=1; |
while (fgets(line, MAXLINE, fic) != NULL) { |
| while (fgets(line, MAXLINE, fic) != NULL) { |
if ((i >= firstobs) && (i <=lastobs)) { |
| if ((i >= firstobs) && (i <=lastobs)) { |
|
| |
|
| for (j=maxwav;j>=1;j--){ |
for (j=maxwav;j>=1;j--){ |
| cutv(stra, strb,line,' '); s[j][i]=atoi(strb); |
cutv(stra, strb,line,' '); s[j][i]=atoi(strb); |
| strcpy(line,stra); |
strcpy(line,stra); |
| cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra); |
cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra); |
| cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra); |
cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra); |
| } |
} |
| |
|
| cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra); |
cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra); |
| cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra); |
cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra); |
| |
|
| cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra); |
cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra); |
| cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra); |
cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra); |
| |
|
| cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra); |
cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra); |
| for (j=ncovcol;j>=1;j--){ |
for (j=ncovcol;j>=1;j--){ |
| cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra); |
cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra); |
| } |
} |
| num[i]=atol(stra); |
num[i]=atol(stra); |
| |
|
| /*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ |
/*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){ |
| printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ |
printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/ |
| |
|
| i=i+1; |
i=i+1; |
| } |
} |
| } |
} |
| /* printf("ii=%d", ij); |
/* printf("ii=%d", ij); |
| scanf("%d",i);*/ |
scanf("%d",i);*/ |
| imx=i-1; /* Number of individuals */ |
imx=i-1; /* Number of individuals */ |
| |
|
| /* for (i=1; i<=imx; i++){ |
/* for (i=1; i<=imx; i++){ |
|
Line 3416 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3606 while((c=getc(ficpar))=='#' && c!= EOF){
|
| Tvard=imatrix(1,15,1,2); |
Tvard=imatrix(1,15,1,2); |
| Tage=ivector(1,15); |
Tage=ivector(1,15); |
| |
|
| if (strlen(model) >1){ |
if (strlen(model) >1){ /* If there is at least 1 covariate */ |
| j=0, j1=0, k1=1, k2=1; |
j=0, j1=0, k1=1, k2=1; |
| j=nbocc(model,'+'); |
j=nbocc(model,'+'); /* j=Number of '+' */ |
| j1=nbocc(model,'*'); |
j1=nbocc(model,'*'); /* j1=Number of '*' */ |
| cptcovn=j+1; |
cptcovn=j+1; |
| cptcovprod=j1; |
cptcovprod=j1; /*Number of products */ |
| |
|
| strcpy(modelsav,model); |
strcpy(modelsav,model); |
| if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ |
if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){ |
|
Line 3430 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3620 while((c=getc(ficpar))=='#' && c!= EOF){
|
| goto end; |
goto end; |
| } |
} |
| |
|
| |
/* This loop fills the array Tvar from the string 'model'.*/ |
| |
|
| for(i=(j+1); i>=1;i--){ |
for(i=(j+1); i>=1;i--){ |
| cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ |
cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */ |
| if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze 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 */ |
if (strchr(strb,'*')) { /* Model includes a product */ |
|
Line 3479 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3671 while((c=getc(ficpar))=='#' && c!= EOF){
|
| } /* end of loop + */ |
} /* end of loop + */ |
| } /* end model */ |
} /* 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. |
| |
If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/ |
| |
|
| /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); |
/* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]); |
| printf("cptcovprod=%d ", cptcovprod); |
printf("cptcovprod=%d ", cptcovprod); |
| fprintf(ficlog,"cptcovprod=%d ", cptcovprod); |
fprintf(ficlog,"cptcovprod=%d ", cptcovprod); |
| scanf("%d ",i);*/ |
|
| fclose(fic); |
scanf("%d ",i); |
| |
fclose(fic);*/ |
| |
|
| /* if(mle==1){*/ |
/* if(mle==1){*/ |
| if (weightopt != 1) { /* Maximisation without weights*/ |
if (weightopt != 1) { /* Maximisation without weights*/ |
| for(i=1;i<=n;i++) weight[i]=1.0; |
for(i=1;i<=n;i++) weight[i]=1.0; |
| } |
} |
| /*-calculation of age at interview from date of interview and age at death -*/ |
/*-calculation of age at interview from date of interview and age at death -*/ |
| agev=matrix(1,maxwav,1,imx); |
agev=matrix(1,maxwav,1,imx); |
| |
|
| for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { |
| for(m=2; (m<= maxwav); m++) { |
for(m=2; (m<= maxwav); m++) { |
| if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){ |
if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){ |
| anint[m][i]=9999; |
anint[m][i]=9999; |
| s[m][i]=-1; |
s[m][i]=-1; |
| } |
} |
| if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1; |
if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1; |
| } |
|
| } |
} |
| |
} |
| |
|
| for (i=1; i<=imx; i++) { |
for (i=1; i<=imx; i++) { |
| agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); |
agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]); |
| for(m=1; (m<= maxwav); m++){ |
for(m=1; (m<= maxwav); m++){ |
| if(s[m][i] >0){ |
if(s[m][i] >0){ |
| if (s[m][i] >= nlstate+1) { |
if (s[m][i] >= nlstate+1) { |
| if(agedc[i]>0) |
if(agedc[i]>0) |
| if(moisdc[i]!=99 && andc[i]!=9999) |
if(moisdc[i]!=99 && andc[i]!=9999) |
| agev[m][i]=agedc[i]; |
agev[m][i]=agedc[i]; |
| /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ |
/*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/ |
| else { |
else { |
| if (andc[i]!=9999){ |
if (andc[i]!=9999){ |
| printf("Warning negative age at death: %d line:%d\n",num[i],i); |
printf("Warning negative age at death: %d line:%d\n",num[i],i); |
| fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i); |
fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i); |
| agev[m][i]=-1; |
agev[m][i]=-1; |
| } |
} |
| } |
} |
| } |
} |
| else if(s[m][i] !=9){ /* Should no more exist */ |
else if(s[m][i] !=9){ /* Should no more exist */ |
| agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); |
agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]); |
| if(mint[m][i]==99 || anint[m][i]==9999) |
if(mint[m][i]==99 || anint[m][i]==9999) |
| agev[m][i]=1; |
|
| else if(agev[m][i] <agemin){ |
|
| agemin=agev[m][i]; |
|
| /*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/ |
|
| } |
|
| else if(agev[m][i] >agemax){ |
|
| agemax=agev[m][i]; |
|
| /* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/ |
|
| } |
|
| /*agev[m][i]=anint[m][i]-annais[i];*/ |
|
| /* agev[m][i] = age[i]+2*m;*/ |
|
| } |
|
| else { /* =9 */ |
|
| agev[m][i]=1; |
agev[m][i]=1; |
| s[m][i]=-1; |
else if(agev[m][i] <agemin){ |
| |
agemin=agev[m][i]; |
| |
/*printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);*/ |
| } |
} |
| |
else if(agev[m][i] >agemax){ |
| |
agemax=agev[m][i]; |
| |
/* printf(" anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.0f\n",m,i,anint[m][i], i,annais[i], agemax);*/ |
| |
} |
| |
/*agev[m][i]=anint[m][i]-annais[i];*/ |
| |
/* agev[m][i] = age[i]+2*m;*/ |
| } |
} |
| else /*= 0 Unknown */ |
else { /* =9 */ |
| agev[m][i]=1; |
agev[m][i]=1; |
| |
s[m][i]=-1; |
| |
} |
| } |
} |
| |
else /*= 0 Unknown */ |
| |
agev[m][i]=1; |
| } |
} |
| for (i=1; i<=imx; i++) { |
|
| for(m=1; (m<= maxwav); m++){ |
} |
| if (s[m][i] > (nlstate+ndeath)) { |
for (i=1; i<=imx; i++) { |
| printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); |
for(m=1; (m<= maxwav); m++){ |
| fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); |
if (s[m][i] > (nlstate+ndeath)) { |
| goto end; |
printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); |
| } |
fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath); |
| |
goto end; |
| } |
} |
| } |
} |
| |
} |
| |
|
| printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); |
printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); |
| fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); |
fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax); |
| |
|
| free_vector(severity,1,maxwav); |
free_vector(severity,1,maxwav); |
| free_imatrix(outcome,1,maxwav+1,1,n); |
free_imatrix(outcome,1,maxwav+1,1,n); |
| free_vector(moisnais,1,n); |
free_vector(moisnais,1,n); |
| free_vector(annais,1,n); |
free_vector(annais,1,n); |
| /* free_matrix(mint,1,maxwav,1,n); |
/* free_matrix(mint,1,maxwav,1,n); |
| free_matrix(anint,1,maxwav,1,n);*/ |
free_matrix(anint,1,maxwav,1,n);*/ |
| free_vector(moisdc,1,n); |
free_vector(moisdc,1,n); |
| free_vector(andc,1,n); |
free_vector(andc,1,n); |
| |
|
| |
|
| wav=ivector(1,imx); |
wav=ivector(1,imx); |
| dh=imatrix(1,lastpass-firstpass+1,1,imx); |
dh=imatrix(1,lastpass-firstpass+1,1,imx); |
| mw=imatrix(1,lastpass-firstpass+1,1,imx); |
bh=imatrix(1,lastpass-firstpass+1,1,imx); |
| |
mw=imatrix(1,lastpass-firstpass+1,1,imx); |
| |
|
| /* Concatenates waves */ |
/* Concatenates waves */ |
| concatwav(wav, dh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); |
concatwav(wav, dh, bh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm); |
| |
|
| |
/* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */ |
| |
|
| Tcode=ivector(1,100); |
Tcode=ivector(1,100); |
| nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX); |
| ncodemax[1]=1; |
ncodemax[1]=1; |
| if (cptcovn > 0) tricode(Tvar,nbcode,imx); |
if (cptcovn > 0) tricode(Tvar,nbcode,imx); |
| |
|
| codtab=imatrix(1,100,1,10); |
codtab=imatrix(1,100,1,10); /* Cross tabulation to get the order of |
| h=0; |
the estimations*/ |
| m=pow(2,cptcoveff); |
h=0; |
| |
m=pow(2,cptcoveff); |
| |
|
| for(k=1;k<=cptcoveff; k++){ |
for(k=1;k<=cptcoveff; k++){ |
| for(i=1; i <=(m/pow(2,k));i++){ |
for(i=1; i <=(m/pow(2,k));i++){ |
| for(j=1; j <= ncodemax[k]; j++){ |
for(j=1; j <= ncodemax[k]; j++){ |
| for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ |
for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ |
| h++; |
h++; |
| if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; |
if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j; |
| /* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/ |
/* printf("h=%d k=%d j=%d codtab[h][k]=%d tvar[k]=%d \n",h, k,j,codtab[h][k],Tvar[k]);*/ |
| } |
} |
| } |
|
| } |
|
| } |
|
| /* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); |
|
| codtab[1][2]=1;codtab[2][2]=2; */ |
|
| /* for(i=1; i <=m ;i++){ |
|
| for(k=1; k <=cptcovn; k++){ |
|
| printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); |
|
| } |
|
| printf("\n"); |
|
| } |
} |
| scanf("%d",i);*/ |
} |
| |
} |
| |
/* printf("codtab[1][2]=%d codtab[2][2]=%d",codtab[1][2],codtab[2][2]); |
| |
codtab[1][2]=1;codtab[2][2]=2; */ |
| |
/* for(i=1; i <=m ;i++){ |
| |
for(k=1; k <=cptcovn; k++){ |
| |
printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff); |
| |
} |
| |
printf("\n"); |
| |
} |
| |
scanf("%d",i);*/ |
| |
|
| /* Calculates basic frequencies. Computes observed prevalence at single age |
/* Calculates basic frequencies. Computes observed prevalence at single age |
| and prints on file fileres'p'. */ |
and prints on file fileres'p'. */ |
| |
|
| |
|
| |
|
| pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
| oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
| newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
| savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
| oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */ |
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */ |
| |
|
| /* For Powell, parameters are in a vector p[] starting at p[1] |
|
| so we point p on param[1][1] so that p[1] maps on param[1][1][1] */ |
/* For Powell, parameters are in a vector p[] starting at p[1] |
| p=param[1][1]; /* *(*(*(param +1)+1)+0) */ |
so we point p on param[1][1] so that p[1] maps on param[1][1][1] */ |
| |
p=param[1][1]; /* *(*(*(param +1)+1)+0) */ |
| |
|
| if(mle==1){ |
if(mle>=1){ /* Could be 1 or 2 */ |
| 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=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate, ndeath, maxwav, weightopt,model); |
| |
|
| |
|
| jk=1; |
jk=1; |
| fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
fprintf(ficres,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
| printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
printf("# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
| fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
fprintf(ficlog,"# Parameters nlstate*nlstate*ncov a12*1 + b12 * age + ...\n"); |
| for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
| for(k=1; k <=(nlstate+ndeath); k++){ |
for(k=1; k <=(nlstate+ndeath); k++){ |
| if (k != i) |
if (k != i) |
| { |
{ |
| printf("%d%d ",i,k); |
printf("%d%d ",i,k); |
| 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("%f ",p[jk]); |
printf("%f ",p[jk]); |
| fprintf(ficlog,"%f ",p[jk]); |
fprintf(ficlog,"%f ",p[jk]); |
| fprintf(ficres,"%f ",p[jk]); |
fprintf(ficres,"%f ",p[jk]); |
| jk++; |
jk++; |
| } |
} |
| printf("\n"); |
printf("\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| fprintf(ficres,"\n"); |
fprintf(ficres,"\n"); |
| } |
} |
| } |
} |
| } |
} |
| if(mle==1){ |
if(mle==1){ |
| /* Computing hessian and covariance matrix */ |
/* Computing hessian and covariance matrix */ |
| ftolhess=ftol; /* Usually correct */ |
ftolhess=ftol; /* Usually correct */ |
| hesscov(matcov, p, npar, delti, ftolhess, func); |
hesscov(matcov, p, npar, delti, ftolhess, func); |
| } |
} |
| 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"); |
| for(i=1,jk=1; i <=nlstate; i++){ |
for(i=1,jk=1; i <=nlstate; i++){ |
| for(j=1; j <=nlstate+ndeath; j++){ |
for(j=1; j <=nlstate+ndeath; j++){ |
| if (j!=i) { |
if (j!=i) { |
| fprintf(ficres,"%1d%1d",i,j); |
fprintf(ficres,"%1d%1d",i,j); |
| printf("%1d%1d",i,j); |
printf("%1d%1d",i,j); |
| fprintf(ficlog,"%1d%1d",i,j); |
fprintf(ficlog,"%1d%1d",i,j); |
| for(k=1; k<=ncovmodel;k++){ |
for(k=1; k<=ncovmodel;k++){ |
| printf(" %.5e",delti[jk]); |
printf(" %.5e",delti[jk]); |
| fprintf(ficlog," %.5e",delti[jk]); |
fprintf(ficlog," %.5e",delti[jk]); |
| fprintf(ficres," %.5e",delti[jk]); |
fprintf(ficres," %.5e",delti[jk]); |
| jk++; |
jk++; |
| } |
} |
| printf("\n"); |
printf("\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| fprintf(ficres,"\n"); |
fprintf(ficres,"\n"); |
| } |
} |
| } |
} |
| } |
} |
| |
|
| k=1; |
fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); |
| fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); |
if(mle==1) |
| if(mle==1) |
printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); |
| printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); |
fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); |
| fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n"); |
for(i=1,k=1;i<=npar;i++){ |
| for(i=1;i<=npar;i++){ |
/* if (k>nlstate) k=1; |
| /* if (k>nlstate) k=1; |
i1=(i-1)/(ncovmodel*nlstate)+1; |
| i1=(i-1)/(ncovmodel*nlstate)+1; |
fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]); |
| fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]); |
printf("%s%d%d",alph[k],i1,tab[i]); |
| printf("%s%d%d",alph[k],i1,tab[i]);*/ |
*/ |
| fprintf(ficres,"%3d",i); |
fprintf(ficres,"%3d",i); |
| if(mle==1) |
if(mle==1) |
| printf("%3d",i); |
printf("%3d",i); |
| fprintf(ficlog,"%3d",i); |
fprintf(ficlog,"%3d",i); |
| for(j=1; j<=i;j++){ |
for(j=1; j<=i;j++){ |
| fprintf(ficres," %.5e",matcov[i][j]); |
fprintf(ficres," %.5e",matcov[i][j]); |
| if(mle==1) |
if(mle==1) |
| printf(" %.5e",matcov[i][j]); |
printf(" %.5e",matcov[i][j]); |
| fprintf(ficlog," %.5e",matcov[i][j]); |
fprintf(ficlog," %.5e",matcov[i][j]); |
| } |
} |
| fprintf(ficres,"\n"); |
fprintf(ficres,"\n"); |
| if(mle==1) |
if(mle==1) |
| printf("\n"); |
printf("\n"); |
| fprintf(ficlog,"\n"); |
fprintf(ficlog,"\n"); |
| k++; |
k++; |
| } |
} |
| |
|
| 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); |
| puts(line); |
puts(line); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| } |
} |
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| estepm=0; |
|
| fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); |
estepm=0; |
| if (estepm==0 || estepm < stepm) estepm=stepm; |
fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm); |
| if (fage <= 2) { |
if (estepm==0 || estepm < stepm) estepm=stepm; |
| bage = ageminpar; |
if (fage <= 2) { |
| fage = agemaxpar; |
bage = ageminpar; |
| } |
fage = agemaxpar; |
| |
} |
| |
|
| 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); |
| |
|
| 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); |
| puts(line); |
puts(line); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| } |
} |
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| |
|
| fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); |
fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav); |
| fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,&mobilav); |
fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); |
| fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,&mobilav); |
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav); |
| |
|
| 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); |
| puts(line); |
puts(line); |
| fputs(line,ficparo); |
fputs(line,ficparo); |
| } |
} |
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| |
|
| |
|
| dateprev1=anprev1+mprev1/12.+jprev1/365.; |
dateprev1=anprev1+mprev1/12.+jprev1/365.; |
| dateprev2=anprev2+mprev2/12.+jprev2/365.; |
dateprev2=anprev2+mprev2/12.+jprev2/365.; |
| |
|
| fscanf(ficpar,"pop_based=%d\n",&popbased); |
fscanf(ficpar,"pop_based=%d\n",&popbased); |
| fprintf(ficparo,"pop_based=%d\n",popbased); |
fprintf(ficparo,"pop_based=%d\n",popbased); |
|
Line 3761 printf("Total number of individuals= %d,
|
Line 3960 printf("Total number of individuals= %d,
|
| ungetc(c,ficpar); |
ungetc(c,ficpar); |
| |
|
| fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2); |
fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2); |
| fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2); |
fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2); |
| fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2); |
fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2); |
| |
|
| |
|
| 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); |
| puts(line); |
puts(line); |
|
Line 3777 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 3976 while((c=getc(ficpar))=='#' && c!= EOF){
|
| fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); |
fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); |
| fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); |
fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1); |
| |
|
| freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); |
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); |
| /*------------ gnuplot -------------*/ |
|
| strcpy(optionfilegnuplot,optionfilefiname); |
/*------------ gnuplot -------------*/ |
| strcat(optionfilegnuplot,".gp"); |
strcpy(optionfilegnuplot,optionfilefiname); |
| if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { |
strcat(optionfilegnuplot,".gp"); |
| printf("Problem with file %s",optionfilegnuplot); |
if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) { |
| } |
printf("Problem with file %s",optionfilegnuplot); |
| else{ |
} |
| fprintf(ficgp,"\n# %s\n", version); |
else{ |
| fprintf(ficgp,"# %s\n", optionfilegnuplot); |
fprintf(ficgp,"\n# %s\n", version); |
| fprintf(ficgp,"set missing 'NaNq'\n"); |
fprintf(ficgp,"# %s\n", optionfilegnuplot); |
| } |
fprintf(ficgp,"set missing 'NaNq'\n"); |
| fclose(ficgp); |
} |
| printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p); |
fclose(ficgp); |
| /*--------- index.htm --------*/ |
printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p); |
| |
/*--------- index.htm --------*/ |
| |
|
| strcpy(optionfilehtm,optionfile); |
strcpy(optionfilehtm,optionfile); |
| strcat(optionfilehtm,".htm"); |
strcat(optionfilehtm,".htm"); |
|
Line 3811 Interval (in months) between two waves:
|
Line 4011 Interval (in months) between two waves:
|
| - Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot); |
- Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot); |
| fclose(fichtm); |
fclose(fichtm); |
| |
|
| printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); |
printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2); |
| |
|
| /*------------ free_vector -------------*/ |
/*------------ free_vector -------------*/ |
| chdir(path); |
chdir(path); |
| |
|
| free_ivector(wav,1,imx); |
free_ivector(wav,1,imx); |
| free_imatrix(dh,1,lastpass-firstpass+1,1,imx); |
free_imatrix(dh,1,lastpass-firstpass+1,1,imx); |
| free_imatrix(mw,1,lastpass-firstpass+1,1,imx); |
free_imatrix(bh,1,lastpass-firstpass+1,1,imx); |
| free_ivector(num,1,n); |
free_imatrix(mw,1,lastpass-firstpass+1,1,imx); |
| free_vector(agedc,1,n); |
free_ivector(num,1,n); |
| /*free_matrix(covar,1,NCOVMAX,1,n);*/ |
free_vector(agedc,1,n); |
| fclose(ficparo); |
/*free_matrix(covar,0,NCOVMAX,1,n);*/ |
| fclose(ficres); |
/*free_matrix(covar,1,NCOVMAX,1,n);*/ |
| |
fclose(ficparo); |
| |
fclose(ficres); |
| |
|
| |
|
| /*--------------- Prevalence limit (stable prevalence) --------------*/ |
/*--------------- Prevalence limit (stable prevalence) --------------*/ |
|
Line 3842 Interval (in months) between two waves:
|
Line 4044 Interval (in months) between two waves:
|
| fprintf(ficrespl,"\n"); |
fprintf(ficrespl,"\n"); |
| |
|
| prlim=matrix(1,nlstate,1,nlstate); |
prlim=matrix(1,nlstate,1,nlstate); |
| pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
|
| oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
|
| newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
|
| savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */ |
|
| oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */ |
|
| k=0; |
|
| agebase=ageminpar; |
agebase=ageminpar; |
| agelim=agemaxpar; |
agelim=agemaxpar; |
| ftolpl=1.e-10; |
ftolpl=1.e-10; |
| i1=cptcoveff; |
i1=cptcoveff; |
| if (cptcovn < 1){i1=1;} |
if (cptcovn < 1){i1=1;} |
| |
|
| for(cptcov=1;cptcov<=i1;cptcov++){ |
for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| k=k+1; |
k=k+1; |
| /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/ |
/*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[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]][codtab[k][j]]); |
| printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]); |
| } |
} |
| fprintf(ficrespl,"******\n"); |
fprintf(ficrespl,"******\n"); |
| printf("******\n"); |
printf("******\n"); |
| fprintf(ficlog,"******\n"); |
fprintf(ficlog,"******\n"); |
| |
|
| for (age=agebase; age<=agelim; age++){ |
for (age=agebase; age<=agelim; age++){ |
| 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(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"); |
| } |
|
| } |
} |
| } |
} |
| |
} |
| fclose(ficrespl); |
fclose(ficrespl); |
| |
|
| /*------------- h Pij x at various ages ------------*/ |
/*------------- h Pij x at various ages ------------*/ |
|
Line 3900 Interval (in months) between two waves:
|
Line 4097 Interval (in months) between two waves:
|
| |
|
| /* hstepm=1; aff par mois*/ |
/* hstepm=1; aff par mois*/ |
| |
|
| k=0; |
for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcov=1;cptcov<=i1;cptcov++){ |
|
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| k=k+1; |
k=k+1; |
| 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]][codtab[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 */ |
| nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ |
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ |
| nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ |
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */ |
| |
|
| /* nhstepm=nhstepm*YEARM; aff par mois*/ |
/* nhstepm=nhstepm*YEARM; aff par mois*/ |
| |
|
| p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); |
hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k); |
| fprintf(ficrespij,"# Age"); |
fprintf(ficrespij,"# Age"); |
| |
for(i=1; i<=nlstate;i++) |
| |
for(j=1; j<=nlstate+ndeath;j++) |
| |
fprintf(ficrespij," %1d-%1d",i,j); |
| |
fprintf(ficrespij,"\n"); |
| |
for (h=0; h<=nhstepm; h++){ |
| |
fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); |
| for(i=1; i<=nlstate;i++) |
for(i=1; i<=nlstate;i++) |
| for(j=1; j<=nlstate+ndeath;j++) |
for(j=1; j<=nlstate+ndeath;j++) |
| fprintf(ficrespij," %1d-%1d",i,j); |
fprintf(ficrespij," %.5f", p3mat[i][j][h]); |
| fprintf(ficrespij,"\n"); |
|
| for (h=0; h<=nhstepm; h++){ |
|
| fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm ); |
|
| for(i=1; i<=nlstate;i++) |
|
| for(j=1; j<=nlstate+ndeath;j++) |
|
| fprintf(ficrespij," %.5f", p3mat[i][j][h]); |
|
| fprintf(ficrespij,"\n"); |
|
| } |
|
| free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
|
| fprintf(ficrespij,"\n"); |
fprintf(ficrespij,"\n"); |
| } |
} |
| |
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm); |
| |
fprintf(ficrespij,"\n"); |
| |
} |
| } |
} |
| } |
} |
| |
|
|
Line 3982 Interval (in months) between two waves:
|
Line 4178 Interval (in months) between two waves:
|
| } |
} |
| printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
| fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv); |
| |
|
| calagedate=-1; |
calagedate=-1; |
| |
|
| prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); |
prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate); |
| |
|
| if (mobilav!=0) { |
if (mobilav!=0) { |
| mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ |
if (movingaverage(probs, bage, fage, mobaverage,mobilav)!=0){ |
|
Line 3992 Interval (in months) between two waves:
|
Line 4191 Interval (in months) between two waves:
|
| } |
} |
| } |
} |
| |
|
| k=0; |
for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcov=1;cptcov<=i1;cptcov++){ |
|
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| k=k+1; |
k=k+1; |
| fprintf(ficrest,"\n#****** "); |
fprintf(ficrest,"\n#****** "); |
|
Line 4020 Interval (in months) between two waves:
|
Line 4218 Interval (in months) between two waves:
|
| varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav); |
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0, mobilav); |
| if(popbased==1){ |
if(popbased==1){ |
| varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav); |
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased,mobilav); |
| } |
} |
| |
|
| |
|
| fprintf(ficrest,"#Total LEs with variances: e.. (std) "); |
fprintf(ficrest,"#Total LEs with variances: e.. (std) "); |
|
Line 4058 Interval (in months) between two waves:
|
Line 4256 Interval (in months) between two waves:
|
| } |
} |
| fprintf(ficrest,"\n"); |
fprintf(ficrest,"\n"); |
| } |
} |
| |
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| |
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
| |
free_vector(epj,1,nlstate+1); |
| } |
} |
| } |
} |
| free_matrix(mint,1,maxwav,1,n); |
free_vector(weight,1,n); |
| free_matrix(anint,1,maxwav,1,n); free_imatrix(s,1,maxwav+1,1,n); |
free_imatrix(Tvard,1,15,1,2); |
| free_vector(weight,1,n); |
free_imatrix(s,1,maxwav+1,1,n); |
| |
free_matrix(anint,1,maxwav,1,n); |
| |
free_matrix(mint,1,maxwav,1,n); |
| |
free_ivector(cod,1,n); |
| |
free_ivector(tab,1,NCOVMAX); |
| fclose(ficreseij); |
fclose(ficreseij); |
| fclose(ficresvij); |
fclose(ficresvij); |
| fclose(ficrest); |
fclose(ficrest); |
| fclose(ficpar); |
fclose(ficpar); |
| free_vector(epj,1,nlstate+1); |
|
| |
|
| /*------- Variance of stable prevalence------*/ |
/*------- Variance of stable prevalence------*/ |
| |
|
|
Line 4079 free_matrix(mint,1,maxwav,1,n);
|
Line 4283 free_matrix(mint,1,maxwav,1,n);
|
| } |
} |
| printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); |
printf("Computing Variance-covariance of stable prevalence: file '%s' \n", fileresvpl); |
| |
|
| k=0; |
for(cptcov=1,k=0;cptcov<=i1;cptcov++){ |
| for(cptcov=1;cptcov<=i1;cptcov++){ |
|
| for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){ |
| k=k+1; |
k=k+1; |
| fprintf(ficresvpl,"\n#****** "); |
fprintf(ficresvpl,"\n#****** "); |
|
Line 4090 free_matrix(mint,1,maxwav,1,n);
|
Line 4293 free_matrix(mint,1,maxwav,1,n);
|
| |
|
| varpl=matrix(1,nlstate,(int) bage, (int) fage); |
varpl=matrix(1,nlstate,(int) bage, (int) fage); |
| oldm=oldms;savm=savms; |
oldm=oldms;savm=savms; |
| varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); |
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k); |
| |
free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
| } |
} |
| } |
} |
| |
|
| fclose(ficresvpl); |
fclose(ficresvpl); |
| |
|
| /*---------- End : free ----------------*/ |
/*---------- End : free ----------------*/ |
| free_matrix(varpl,1,nlstate,(int) bage, (int)fage); |
|
| |
|
| free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
|
| free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage); |
|
| |
|
| |
|
| free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath); |
| |
|
| |
free_matrix(covar,0,NCOVMAX,1,n); |
| free_matrix(matcov,1,npar,1,npar); |
free_matrix(matcov,1,npar,1,npar); |
| free_vector(delti,1,npar); |
free_vector(delti,1,npar); |
| free_matrix(agev,1,maxwav,1,imx); |
free_matrix(agev,1,maxwav,1,imx); |
| free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); |
| if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX); |
| |
free_ivector(ncodemax,1,8); |
| |
free_ivector(Tvar,1,15); |
| |
free_ivector(Tprod,1,15); |
| |
free_ivector(Tvaraff,1,15); |
| |
free_ivector(Tage,1,15); |
| |
free_ivector(Tcode,1,100); |
| |
|
| fprintf(fichtm,"\n</body>"); |
fprintf(fichtm,"\n</body>"); |
| fclose(fichtm); |
fclose(fichtm); |
|
Line 4134 free_matrix(mint,1,maxwav,1,n);
|
Line 4339 free_matrix(mint,1,maxwav,1,n);
|
| /*printf("Total time was %d uSec.\n", total_usecs);*/ |
/*printf("Total time was %d uSec.\n", total_usecs);*/ |
| /*------ End -----------*/ |
/*------ End -----------*/ |
| |
|
| |
end: |
| end: |
|
| #ifdef windows |
#ifdef windows |
| /* chdir(pathcd);*/ |
/* chdir(pathcd);*/ |
| #endif |
#endif |
|
Line 4143 free_matrix(mint,1,maxwav,1,n);
|
Line 4347 free_matrix(mint,1,maxwav,1,n);
|
| /*system("../gp37mgw/wgnuplot graph.plt");*/ |
/*system("../gp37mgw/wgnuplot graph.plt");*/ |
| /*system("cd ../gp37mgw");*/ |
/*system("cd ../gp37mgw");*/ |
| /* system("..\\gp37mgw\\wgnuplot graph.plt");*/ |
/* system("..\\gp37mgw\\wgnuplot graph.plt");*/ |
| strcpy(plotcmd,GNUPLOTPROGRAM); |
strcpy(plotcmd,GNUPLOTPROGRAM); |
| strcat(plotcmd," "); |
strcat(plotcmd," "); |
| strcat(plotcmd,optionfilegnuplot); |
strcat(plotcmd,optionfilegnuplot); |
| printf("Starting: %s\n",plotcmd);fflush(stdout); |
printf("Starting: %s\n",plotcmd);fflush(stdout); |
| system(plotcmd); |
system(plotcmd); |
| |
|
| /*#ifdef windows*/ |
/*#ifdef windows*/ |
| while (z[0] != 'q') { |
while (z[0] != 'q') { |
![[BACK]](/icons/small/back.gif){kind=icon}
Return to imach.c CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src