version 1.14, 2002/02/20 17:05:44
|
version 1.22, 2002/02/22 17:54:20
|
Line 1
|
Line 1
|
|
/* $Id$
|
/*********************** Imach **************************************
|
Interpolate Markov Chain
|
This program computes Healthy Life Expectancies from cross-longitudinal
|
|
data. Cross-longitudinal consist in a first survey ("cross") where
|
Short summary of the programme:
|
individuals from different ages are interviewed on their health status
|
|
or degree of disability. At least a second wave of interviews
|
This program computes Healthy Life Expectancies from
|
("longitudinal") should measure each new individual health status.
|
cross-longitudinal data. Cross-longitudinal data consist in: -1- a
|
Health expectancies are computed from the transistions observed between
|
first survey ("cross") where individuals from different ages are
|
waves and are computed for each degree of severity of disability (number
|
interviewed on their health status or degree of disability (in the
|
of life states). More degrees you consider, more time is necessary to
|
case of a health survey which is our main interest) -2- at least a
|
reach the Maximum Likelihood of the parameters involved in the model.
|
second wave of interviews ("longitudinal") which measure each change
|
The simplest model is the multinomial logistic model where pij is
|
(if any) in individual health status. Health expectancies are
|
the probabibility to be observed in state j at the second wave conditional
|
computed from the time spent in each health state according to a
|
to be observed in state i at the first wave. Therefore the model is:
|
model. More health states you consider, more time is necessary to reach the
|
log(pij/pii)= aij + bij*age+ cij*sex + etc , where 'age' is age and 'sex'
|
Maximum Likelihood of the parameters involved in the model. The
|
is a covariate. If you want to have a more complex model than "constant and
|
simplest model is the multinomial logistic model where pij is the
|
age", you should modify the program where the markup
|
probabibility to be observed in state j at the second wave
|
*Covariates have to be included here again* invites you to do it.
|
conditional to be observed in state i at the first wave. Therefore
|
More covariates you add, less is the speed of the convergence.
|
the model is: log(pij/pii)= aij + bij*age+ cij*sex + etc , where
|
|
'age' is age and 'sex' is a covariate. If you want to have a more
|
The advantage that this computer programme claims, comes from that if the
|
complex model than "constant and age", you should modify the program
|
delay between waves is not identical for each individual, or if some
|
where the markup *Covariates have to be included here again* invites
|
individual missed an interview, the information is not rounded or lost, but
|
you to do it. More covariates you add, slower the
|
taken into account using an interpolation or extrapolation.
|
convergence.
|
hPijx is the probability to be
|
|
observed in state i at age x+h conditional to the observed state i at age
|
The advantage of this computer programme, compared to a simple
|
x. The delay 'h' can be split into an exact number (nh*stepm) of
|
multinomial logistic model, is clear when the delay between waves is not
|
unobserved intermediate states. This elementary transition (by month or
|
identical for each individual. Also, if a individual missed an
|
quarter trimester, semester or year) is model as a multinomial logistic.
|
intermediate interview, the information is lost, but taken into
|
The hPx matrix is simply the matrix product of nh*stepm elementary matrices
|
account using an interpolation or extrapolation.
|
and the contribution of each individual to the likelihood is simply hPijx.
|
|
|
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
|
|
split into an exact number (nh*stepm) of unobserved intermediate
|
|
states. This elementary transition (by month or quarter trimester,
|
|
semester or year) is model as a multinomial logistic. The hPx
|
|
matrix is simply the matrix product of nh*stepm elementary matrices
|
|
and the contribution of each individual to the likelihood is simply
|
|
hPijx.
|
|
|
Also this programme outputs the covariance matrix of the parameters but also
|
Also this programme outputs the covariance matrix of the parameters but also
|
of the life expectancies. It also computes the prevalence limits.
|
of the life expectancies. It also computes the prevalence limits.
|
Line 48
|
Line 56
|
#include <unistd.h>
|
#include <unistd.h>
|
|
|
#define MAXLINE 256
|
#define MAXLINE 256
|
|
#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"
|
#define FILENAMELENGTH 80
|
#define FILENAMELENGTH 80
|
/*#define DEBUG*/
|
/*#define DEBUG*/
|
#define windows
|
#define windows
|
Line 67
|
Line 76
|
#define AGEBASE 40
|
#define AGEBASE 40
|
|
|
|
|
|
int erreur; /* Error number */
|
int nvar;
|
int nvar;
|
int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
|
int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
|
int npar=NPARMAX;
|
int npar=NPARMAX;
|
int nlstate=2; /* Number of live states */
|
int nlstate=2; /* Number of live states */
|
int ndeath=1; /* Number of dead states */
|
int ndeath=1; /* Number of dead states */
|
int ncovmodel, ncov; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
|
int ncovmodel, ncov; /* Total number of covariables including constant a12*1 +b12*x ncovmodel=2 */
|
int popbased=0, fprev,lprev;
|
int popbased=0;
|
|
|
int *wav; /* Number of waves for this individuual 0 is possible */
|
int *wav; /* Number of waves for this individuual 0 is possible */
|
int maxwav; /* Maxim number of waves */
|
int maxwav; /* Maxim number of waves */
|
Line 85 double jmean; /* Mean space between 2 wa
|
Line 95 double jmean; /* Mean space between 2 wa
|
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 */
|
FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf;
|
FILE *fic,*ficpar, *ficparo,*ficres, *ficrespl, *ficrespij, *ficrest,*ficresf;
|
FILE *ficgp, *fichtm,*ficresprob;
|
FILE *ficgp, *fichtm,*ficresprob,*ficpop;
|
FILE *ficreseij;
|
FILE *ficreseij;
|
char filerese[FILENAMELENGTH];
|
char filerese[FILENAMELENGTH];
|
FILE *ficresvij;
|
FILE *ficresvij;
|
Line 129 int m,nb;
|
Line 139 int m,nb;
|
int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
|
int *num, firstpass=0, lastpass=4,*cod, *ncodemax, *Tage;
|
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
|
double **agev,*moisnais, *annais, *moisdc, *andc,**mint, **anint;
|
double **pmmij, ***probs, ***mobaverage;
|
double **pmmij, ***probs, ***mobaverage;
|
|
double dateintmean=0;
|
|
|
double *weight;
|
double *weight;
|
int **s; /* Status */
|
int **s; /* Status */
|
Line 139 double ftol=FTOL; /* Tolerance for compu
|
Line 150 double ftol=FTOL; /* Tolerance for compu
|
double ftolhess; /* Tolerance for computing hessian */
|
double ftolhess; /* Tolerance for computing hessian */
|
|
|
/**************** split *************************/
|
/**************** split *************************/
|
static int split( char *path, char *dirc, char *name )
|
static int split( char *path, char *dirc, char *name, char *ext, char *finame )
|
{
|
{
|
char *s; /* pointer */
|
char *s; /* 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 );
|
|
#ifdef windows
|
s = strrchr( path, '\\' ); /* find last / */
|
s = strrchr( path, '\\' ); /* find last / */
|
|
#else
|
|
s = strrchr( path, '/' ); /* find last / */
|
|
#endif
|
if ( s == NULL ) { /* no directory, so use current */
|
if ( s == NULL ) { /* no directory, so use current */
|
#if defined(__bsd__) /* get current working directory */
|
#if defined(__bsd__) /* get current working directory */
|
extern char *getwd( );
|
extern char *getwd( );
|
Line 169 static int split( char *path, char *dirc
|
Line 184 static int split( char *path, char *dirc
|
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
|
if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
|
if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }
|
|
#else
|
|
if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }
|
|
#endif
|
|
s = strrchr( name, '.' ); /* find last / */
|
|
s++;
|
|
strcpy(ext,s); /* save extension */
|
|
l1= strlen( name);
|
|
l2= strlen( s)+1;
|
|
strncpy( finame, name, l1-l2);
|
|
finame[l1-l2]= 0;
|
return( 0 ); /* we're done */
|
return( 0 ); /* we're done */
|
}
|
}
|
|
|
Line 717 double **pmij(double **ps, double *cov,
|
Line 743 double **pmij(double **ps, double *cov,
|
s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
|
s2 += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];
|
/*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
|
/*printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);*/
|
}
|
}
|
ps[i][j]=(s2);
|
ps[i][j]=s2;
|
}
|
}
|
}
|
}
|
/*ps[3][2]=1;*/
|
/*ps[3][2]=1;*/
|
Line 900 void mlikeli(FILE *ficres,double p[], in
|
Line 926 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(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f ",iter,func(p));
|
fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
|
|
|
}
|
}
|
|
|
Line 1150 void lubksb(double **a, int n, int *indx
|
Line 1176 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)
|
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)
|
{ /* 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;
|
double ***freq; /* Frequencies */
|
double ***freq; /* Frequencies */
|
double *pp;
|
double *pp;
|
double pos;
|
double pos, k2, dateintsum=0,k2cpt=0;
|
FILE *ficresp;
|
FILE *ficresp;
|
char fileresp[FILENAMELENGTH];
|
char fileresp[FILENAMELENGTH];
|
|
|
pp=vector(1,nlstate);
|
pp=vector(1,nlstate);
|
probs= ma3x(1,130 ,1,8, 1,8);
|
probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
strcpy(fileresp,"p");
|
strcpy(fileresp,"p");
|
strcat(fileresp,fileres);
|
strcat(fileresp,fileres);
|
if((ficresp=fopen(fileresp,"w"))==NULL) {
|
if((ficresp=fopen(fileresp,"w"))==NULL) {
|
Line 1183 void freqsummary(char fileres[], int ag
|
Line 1209 void freqsummary(char fileres[], int ag
|
for (jk=-1; jk<=nlstate+ndeath; jk++)
|
for (jk=-1; jk<=nlstate+ndeath; jk++)
|
for(m=agemin; m <= agemax+3; m++)
|
for(m=agemin; m <= agemax+3; m++)
|
freq[i][jk][m]=0;
|
freq[i][jk][m]=0;
|
|
|
|
dateintsum=0;
|
|
k2cpt=0;
|
for (i=1; i<=imx; i++) {
|
for (i=1; i<=imx; i++) {
|
bool=1;
|
bool=1;
|
if (cptcovn>0) {
|
if (cptcovn>0) {
|
Line 1191 void freqsummary(char fileres[], int ag
|
Line 1219 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=fprev; m<=lprev; m++){
|
for(m=firstpass; m<=lastpass; m++){
|
if(agev[m][i]==0) agev[m][i]=agemax+1;
|
k2=anint[m][i]+(mint[m][i]/12.);
|
if(agev[m][i]==1) agev[m][i]=agemax+2;
|
if ((k2>=dateprev1) && (k2<=dateprev2)) {
|
freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
|
if(agev[m][i]==0) agev[m][i]=agemax+1;
|
freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
|
if(agev[m][i]==1) agev[m][i]=agemax+2;
|
|
freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];
|
|
freq[s[m][i]][s[m+1][i]][(int) agemax+3] += weight[i];
|
|
if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {
|
|
dateintsum=dateintsum+k2;
|
|
k2cpt++;
|
|
}
|
|
|
|
}
|
}
|
}
|
}
|
}
|
}
|
}
|
Line 1258 void freqsummary(char fileres[], int ag
|
Line 1294 void freqsummary(char fileres[], int ag
|
}
|
}
|
}
|
}
|
}
|
}
|
|
dateintmean=dateintsum/k2cpt;
|
|
|
fclose(ficresp);
|
fclose(ficresp);
|
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
free_vector(pp,1,nlstate);
|
free_vector(pp,1,nlstate);
|
|
|
|
/* End of Freq */
|
|
}
|
|
|
|
/************ Prevalence ********************/
|
|
void prevalence(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 calagedate)
|
|
{ /* Some frequencies */
|
|
|
|
int i, m, jk, k1, i1, j1, bool, z1,z2,j;
|
|
double ***freq; /* Frequencies */
|
|
double *pp;
|
|
double pos, k2;
|
|
|
|
pp=vector(1,nlstate);
|
|
probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
|
|
freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
|
|
j1=0;
|
|
|
|
j=cptcoveff;
|
|
if (cptcovn<1) {j=1;ncodemax[1]=1;}
|
|
|
|
for(k1=1; k1<=j;k1++){
|
|
for(i1=1; i1<=ncodemax[k1];i1++){
|
|
j1++;
|
|
|
|
for (i=-1; i<=nlstate+ndeath; i++)
|
|
for (jk=-1; jk<=nlstate+ndeath; jk++)
|
|
for(m=agemin; m <= agemax+3; m++)
|
|
freq[i][jk][m]=0;
|
|
|
|
for (i=1; i<=imx; i++) {
|
|
bool=1;
|
|
if (cptcovn>0) {
|
|
for (z1=1; z1<=cptcoveff; z1++)
|
|
if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
|
|
bool=0;
|
|
}
|
|
if (bool==1) {
|
|
for(m=firstpass; m<=lastpass; m++){
|
|
k2=anint[m][i]+(mint[m][i]/12.);
|
|
if ((k2>=dateprev1) && (k2<=dateprev2)) {
|
|
if(agev[m][i]==0) agev[m][i]=agemax+1;
|
|
if(agev[m][i]==1) agev[m][i]=agemax+2;
|
|
freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];
|
|
freq[s[m][i]][s[m+1][i]][(int)(agemax+3+1)] += weight[i];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
for(i=(int)agemin; i <= (int)agemax+3; i++){
|
|
for(jk=1; jk <=nlstate ; jk++){
|
|
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
|
|
pp[jk] += freq[jk][m][i];
|
|
}
|
|
for(jk=1; jk <=nlstate ; jk++){
|
|
for(m=-1, pos=0; m <=0 ; m++)
|
|
pos += freq[jk][m][i];
|
|
}
|
|
|
|
for(jk=1; jk <=nlstate ; jk++){
|
|
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
|
|
pp[jk] += freq[jk][m][i];
|
|
}
|
|
|
|
for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];
|
|
|
|
for(jk=1; jk <=nlstate ; jk++){
|
|
if( i <= (int) agemax){
|
|
if(pos>=1.e-5){
|
|
probs[i][jk][j1]= pp[jk]/pos;
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);
|
|
free_vector(pp,1,nlstate);
|
|
|
} /* End of Freq */
|
} /* End of Freq */
|
|
|
/************* Waves Concatenation ***************/
|
/************* Waves Concatenation ***************/
|
Line 1763 fclose(ficresprob);
|
Line 1883 fclose(ficresprob);
|
/**************** Main Program *****************/
|
/**************** Main Program *****************/
|
/***********************************************/
|
/***********************************************/
|
|
|
/*int main(int argc, char *argv[])*/
|
int main(int argc, char *argv[])
|
int main()
|
|
{
|
{
|
|
|
int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
|
int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
|
Line 1780 int main()
|
Line 1899 int main()
|
char line[MAXLINE], linepar[MAXLINE];
|
char line[MAXLINE], linepar[MAXLINE];
|
char title[MAXLINE];
|
char title[MAXLINE];
|
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
|
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
|
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];
|
char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];
|
|
|
|
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], fileresf[FILENAMELENGTH];;
|
|
|
char filerest[FILENAMELENGTH];
|
char filerest[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
char fileregp[FILENAMELENGTH];
|
|
char popfile[FILENAMELENGTH];
|
char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
|
char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
|
int firstobs=1, lastobs=10;
|
int firstobs=1, lastobs=10;
|
int sdeb, sfin; /* Status at beginning and end */
|
int sdeb, sfin; /* Status at beginning and end */
|
Line 1790 int main()
|
Line 1913 int main()
|
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,**adl,*tab;
|
int mobilav=0, fprevfore=1, lprevfore=1;
|
int mobilav=0,popforecast=0;
|
int hstepm, nhstepm;
|
int hstepm, nhstepm;
|
|
int *popage;/*boolprev=0 if date and zero if wave*/
|
|
double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2;
|
|
|
double bage, fage, age, agelim, agebase;
|
double bage, fage, age, agelim, agebase;
|
double ftolpl=FTOL;
|
double ftolpl=FTOL;
|
Line 1805 int main()
|
Line 1930 int main()
|
double ***eij, ***vareij;
|
double ***eij, ***vareij;
|
double **varpl; /* Variances of prevalence limits by age */
|
double **varpl; /* Variances of prevalence limits by age */
|
double *epj, vepp;
|
double *epj, vepp;
|
double kk1;
|
double kk1, kk2;
|
|
double *popeffectif,*popcount;
|
|
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,jprojmean,mprojmean,anprojmean, calagedate;
|
|
double yp,yp1,yp2;
|
|
|
char version[80]="Imach version 64b, May 2001, INED-EUROREVES ";
|
char version[80]="Imach version 0.7, February 2002, INED-EUROREVES ";
|
char *alph[]={"a","a","b","c","d","e"}, str[4];
|
char *alph[]={"a","a","b","c","d","e"}, str[4];
|
|
|
|
|
Line 1815 int main()
|
Line 1943 int main()
|
#include <sys/time.h>
|
#include <sys/time.h>
|
#include <time.h>
|
#include <time.h>
|
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 */
|
|
|
|
|
printf("\nIMACH, Version 0.64b");
|
printf("\n%s",version);
|
printf("\nEnter the parameter file name: ");
|
if(argc <=1){
|
|
printf("\nEnter the parameter file name: ");
|
#ifdef windows
|
scanf("%s",pathtot);
|
scanf("%s",pathtot);
|
}
|
getcwd(pathcd, size);
|
else{
|
|
strcpy(pathtot,argv[1]);
|
|
}
|
|
/*if(getcwd(pathcd, 80)!= NULL)printf ("Error pathcd\n");*/
|
/*cygwin_split_path(pathtot,path,optionfile);
|
/*cygwin_split_path(pathtot,path,optionfile);
|
printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
|
printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
|
/* cutv(path,optionfile,pathtot,'\\');*/
|
/* cutv(path,optionfile,pathtot,'\\');*/
|
|
|
split(pathtot, path,optionfile);
|
split(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);
|
#endif
|
|
#ifdef unix
|
|
scanf("%s",optionfile);
|
|
#endif
|
|
|
|
/*-------- arguments in the command line --------*/
|
/*-------- arguments in the command line --------*/
|
|
|
strcpy(fileres,"r");
|
strcpy(fileres,"r");
|
strcat(fileres, optionfile);
|
strcat(fileres, optionfilefiname);
|
|
strcat(fileres,".txt"); /* Other files have txt extension */
|
|
|
/*---------arguments file --------*/
|
/*---------arguments file --------*/
|
|
|
Line 1877 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 2007 while((c=getc(ficpar))=='#' && c!= EOF){
|
}
|
}
|
ungetc(c,ficpar);
|
ungetc(c,ficpar);
|
|
|
fscanf(ficpar,"fprevalence=%d lprevalence=%d pop_based=%d\n",&fprev,&lprev,&popbased);
|
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
|
ungetc(c,ficpar);
|
|
fgets(line, MAXLINE, ficpar);
|
|
puts(line);
|
|
fputs(line,ficparo);
|
|
}
|
|
ungetc(c,ficpar);
|
|
|
|
fscanf(ficpar,"fprevalence=%d lprevalence=%d mob_average=%d\n",&fprevfore,&lprevfore,&mobilav);
|
|
|
|
covar=matrix(0,NCOVMAX,1,n);
|
covar=matrix(0,NCOVMAX,1,n);
|
cptcovn=0;
|
cptcovn=0;
|
if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
|
if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;
|
Line 2050 while((c=getc(ficpar))=='#' && c!= EOF){
|
Line 2170 while((c=getc(ficpar))=='#' && c!= EOF){
|
if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
|
if ((s[2][i]==3) && (s[3][i]==2)) s[3][i]=3;
|
if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
|
if ((s[3][i]==3) && (s[4][i]==2)) s[4][i]=3;
|
}
|
}
|
for (i=1; i<=imx; i++) 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]));*/
|
|
|
for (i=1; i<=imx; i++)
|
|
if (covar[1][i]==0) 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]));*/
|
|
|
/* Calculation of the number of parameter from char model*/
|
/* Calculation of the number of parameter from char model*/
|
Tvar=ivector(1,15);
|
Tvar=ivector(1,15);
|
Line 2196 printf("Total number of individuals= %d,
|
Line 2318 printf("Total number of individuals= %d,
|
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);
|
|
|
Line 2241 printf("Total number of individuals= %d,
|
Line 2363 printf("Total number of individuals= %d,
|
|
|
/* 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'. */
|
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax);
|
|
|
|
|
|
|
|
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 */
|
Line 2258 printf("Total number of individuals= %d,
|
Line 2381 printf("Total number of individuals= %d,
|
}
|
}
|
|
|
/*--------- results files --------------*/
|
/*--------- results files --------------*/
|
fprintf(ficres,"\ntitle=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);
|
fprintf(ficres,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);
|
|
|
|
|
jk=1;
|
jk=1;
|
fprintf(ficres,"# Parameters\n");
|
fprintf(ficres,"# Parameters\n");
|
printf("# Parameters\n");
|
printf("# Parameters\n");
|
Line 2300 printf("Total number of individuals= %d,
|
Line 2424 printf("Total number of individuals= %d,
|
fprintf(ficres,"\n");
|
fprintf(ficres,"\n");
|
}
|
}
|
}
|
}
|
}
|
}
|
|
|
k=1;
|
k=1;
|
fprintf(ficres,"# Covariance\n");
|
fprintf(ficres,"# Covariance\n");
|
Line 2335 printf("Total number of individuals= %d,
|
Line 2459 printf("Total number of individuals= %d,
|
bage = agemin;
|
bage = agemin;
|
fage = agemax;
|
fage = agemax;
|
}
|
}
|
|
|
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\n",agemin,agemax,bage,fage);
|
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
|
|
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
|
|
|
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
|
ungetc(c,ficpar);
|
|
fgets(line, MAXLINE, ficpar);
|
|
puts(line);
|
|
fputs(line,ficparo);
|
|
}
|
|
ungetc(c,ficpar);
|
|
|
|
fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mob_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
|
|
fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
|
|
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mob_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
|
|
|
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
|
ungetc(c,ficpar);
|
|
fgets(line, MAXLINE, ficpar);
|
|
puts(line);
|
|
fputs(line,ficparo);
|
|
}
|
|
ungetc(c,ficpar);
|
|
|
|
|
|
dateprev1=anprev1+mprev1/12.+jprev1/365.;
|
/*------------ gnuplot -------------*/
|
dateprev2=anprev2+mprev2/12.+jprev2/365.;
|
chdir(pathcd);
|
|
if((ficgp=fopen("graph.plt","w"))==NULL) {
|
fscanf(ficpar,"pop_based=%d\n",&popbased);
|
printf("Problem with file graph.gp");goto end;
|
fprintf(ficparo,"pop_based=%d\n",popbased);
|
|
fprintf(ficres,"pop_based=%d\n",popbased);
|
|
|
|
while((c=getc(ficpar))=='#' && c!= EOF){
|
|
ungetc(c,ficpar);
|
|
fgets(line, MAXLINE, ficpar);
|
|
puts(line);
|
|
fputs(line,ficparo);
|
}
|
}
|
|
ungetc(c,ficpar);
|
|
fscanf(ficpar,"popforecast=%d popfile=%s starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&popforecast,popfile,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2);
|
|
fprintf(ficparo,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
|
|
fprintf(ficres,"popforecast=%d popfile=%s starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",popforecast,popfile,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);
|
|
|
|
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2);
|
|
|
|
|
|
/*------------ gnuplot -------------*/
|
|
/*chdir(pathcd);*/
|
|
strcpy(optionfilegnuplot,optionfilefiname);
|
|
strcat(optionfilegnuplot,".plt");
|
|
if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
|
|
printf("Problem with file %s",optionfilegnuplot);goto end;
|
|
}
|
#ifdef windows
|
#ifdef windows
|
fprintf(ficgp,"cd \"%s\" \n",pathc);
|
fprintf(ficgp,"cd \"%s\" \n",pathc);
|
#endif
|
#endif
|
m=pow(2,cptcoveff);
|
m=pow(2,cptcoveff);
|
|
|
Line 2503 ij=1;
|
Line 2671 ij=1;
|
}
|
}
|
|
|
fclose(ficgp);
|
fclose(ficgp);
|
|
/* end gnuplot */
|
|
|
chdir(path);
|
chdir(path);
|
free_matrix(agev,1,maxwav,1,imx);
|
|
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(mw,1,lastpass-firstpass+1,1,imx);
|
|
|
free_imatrix(s,1,maxwav+1,1,n);
|
|
|
|
|
|
free_ivector(num,1,n);
|
free_ivector(num,1,n);
|
free_vector(agedc,1,n);
|
free_vector(agedc,1,n);
|
free_vector(weight,1,n);
|
|
/*free_matrix(covar,1,NCOVMAX,1,n);*/
|
/*free_matrix(covar,1,NCOVMAX,1,n);*/
|
fclose(ficparo);
|
fclose(ficparo);
|
fclose(ficres);
|
fclose(ficres);
|
Line 2546 chdir(path);
|
Line 2710 chdir(path);
|
printf("Problem with %s \n",optionfilehtm);goto end;
|
printf("Problem with %s \n",optionfilehtm);goto end;
|
}
|
}
|
|
|
fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.64b </font> <hr size=\"2\" color=\"#EC5E5E\">
|
fprintf(fichtm,"<body><ul> <font size=\"6\">Imach, Version 0.7 </font> <hr size=\"2\" color=\"#EC5E5E\">
|
Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
|
Titre=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>
|
Total number of observations=%d <br>
|
Total number of observations=%d <br>
|
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
|
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
|
Line 2561 Interval (in months) between two waves:
|
Line 2725 Interval (in months) between two waves:
|
- Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
|
- Variances of life expectancies by age and initial health status: <a href=\"v%s\">v%s</a><br>
|
- Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
|
- Health expectancies with their variances: <a href=\"t%s\">t%s</a> <br>
|
- Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
|
- Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%s</a> <br>
|
- Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>
|
- Prevalences and population forecasting: <a href=\"f%s\">f%s</a> <br>
|
<br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
|
<br>",title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);
|
|
|
fprintf(fichtm," <li>Graphs</li><p>");
|
fprintf(fichtm," <li>Graphs</li><p>");
|
Line 2700 fclose(fichtm);
|
Line 2864 fclose(fichtm);
|
|
|
fclose(ficrespij);
|
fclose(ficrespij);
|
|
|
|
if(stepm == 1) {
|
/*---------- Forecasting ------------------*/
|
/*---------- Forecasting ------------------*/
|
|
calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;
|
|
|
|
/*printf("calage= %f", calagedate);*/
|
|
|
|
prevalence(agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);
|
|
|
|
|
strcpy(fileresf,"f");
|
strcpy(fileresf,"f");
|
strcat(fileresf,fileres);
|
strcat(fileresf,fileres);
|
Line 2709 fclose(fichtm);
|
Line 2880 fclose(fichtm);
|
}
|
}
|
printf("Computing forecasting: result on file '%s' \n", fileresf);
|
printf("Computing forecasting: result on file '%s' \n", fileresf);
|
|
|
|
free_matrix(mint,1,maxwav,1,n);
|
|
free_matrix(anint,1,maxwav,1,n);
|
|
free_matrix(agev,1,maxwav,1,imx);
|
/* Mobile average */
|
/* Mobile average */
|
|
|
/* for (agedeb=bage; agedeb<=fage; agedeb++)
|
|
for (i=1; i<=nlstate;i++)
|
|
for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++)
|
|
printf("%f %d i=%d j1=%d\n", probs[(int)agedeb][i][cptcod],(int) agedeb,i,cptcod);*/
|
|
|
|
if (cptcoveff==0) ncodemax[cptcoveff]=1;
|
if (cptcoveff==0) ncodemax[cptcoveff]=1;
|
|
|
if (mobilav==1) {
|
if (mobilav==1) {
|
Line 2738 fclose(fichtm);
|
Line 2907 fclose(fichtm);
|
}
|
}
|
|
|
stepsize=(int) (stepm+YEARM-1)/YEARM;
|
stepsize=(int) (stepm+YEARM-1)/YEARM;
|
if (stepm<=24) stepsize=2;
|
if (stepm<=12) stepsize=1;
|
|
|
agelim=AGESUP;
|
agelim=AGESUP;
|
hstepm=stepsize*YEARM; /* Every year of age */
|
/*hstepm=stepsize*YEARM; *//* Every year of age */
|
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */
|
hstepm=1;
|
hstepm=12;
|
hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */
|
k=0;
|
yp1=modf(dateintmean,&yp);
|
|
anprojmean=yp;
|
|
yp2=modf((yp1*12),&yp);
|
|
mprojmean=yp;
|
|
yp1=modf((yp2*30.5),&yp);
|
|
jprojmean=yp;
|
|
if(jprojmean==0) jprojmean=1;
|
|
if(mprojmean==0) jprojmean=1;
|
|
|
|
fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean);
|
|
|
|
if (popforecast==1) {
|
|
if((ficpop=fopen(popfile,"r"))==NULL) {
|
|
printf("Problem with population file : %s\n",popfile);goto end;
|
|
}
|
|
popage=ivector(0,AGESUP);
|
|
popeffectif=vector(0,AGESUP);
|
|
popcount=vector(0,AGESUP);
|
|
|
|
i=1;
|
|
while ((c=fscanf(ficpop,"%d %lf\n",&popage[i],&popcount[i])) != EOF)
|
|
{
|
|
i=i+1;
|
|
}
|
|
imx=i;
|
|
|
|
for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
|
|
}
|
|
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcov=1;cptcov<=i1;cptcov++){
|
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
|
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
|
k=k+1;
|
k=k+1;
|
fprintf(ficresf,"\n#****** ");
|
fprintf(ficresf,"\n#******");
|
for(j=1;j<=cptcoveff;j++) {
|
for(j=1;j<=cptcoveff;j++) {
|
fprintf(ficresf,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
|
}
|
}
|
|
|
fprintf(ficresf,"******\n");
|
fprintf(ficresf,"******\n");
|
|
fprintf(ficresf,"# StartingAge FinalAge");
|
fprintf(ficresf,"# StartingAge FinalAge Horizon(in years)");
|
|
for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
|
for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
|
|
if (popforecast==1) fprintf(ficresf," [Population]");
|
for (agedeb=fage; agedeb>=bage; agedeb--){
|
|
fprintf(ficresf,"\n%d %.f %.f 0 ",k,agedeb, agedeb);
|
for (cpt=0; cpt<4;cpt++) {
|
if (mobilav==1) {
|
fprintf(ficresf,"\n");
|
for(j=1; j<=nlstate;j++)
|
fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);
|
fprintf(ficresf,"%.5f ",mobaverage[(int)agedeb][j][cptcod]);
|
|
}
|
for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */
|
else {
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
|
for(j=1; j<=nlstate;j++)
|
nhstepm = nhstepm/hstepm;
|
fprintf(ficresf,"%.5f ",probs[(int)agedeb][j][cptcod]);
|
/*printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);*/
|
}
|
|
|
|
for(j=1; j<=ndeath;j++) fprintf(ficresf,"0.");
|
|
}
|
|
for (cpt=1; cpt<=NCOVMAX;cpt++)
|
|
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
|
|
|
|
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */
|
|
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
|
|
/*printf("stepm=%d hstepm=%d nhstepm=%d \n",stepm,hstepm,nhstepm);*/
|
|
|
|
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);
|
|
|
for (h=0; h<=nhstepm; h++){
|
|
|
|
if (h*hstepm/YEARM*stepm==cpt)
|
for (h=0; h<=nhstepm; h++){
|
fprintf(ficresf,"\n%d %.f %.f %.f",k,agedeb, agedeb+ h*hstepm/YEARM*stepm, h*hstepm/YEARM*stepm);
|
if (h==(int) (calagedate+YEARM*cpt)) {
|
|
fprintf(ficresf,"\n %.f ",agedeb+h*hstepm/YEARM*stepm);
|
|
}
|
for(j=1; j<=nlstate+ndeath;j++) {
|
for(j=1; j<=nlstate+ndeath;j++) {
|
kk1=0.;
|
kk1=0.;kk2=0;
|
for(i=1; i<=nlstate;i++) {
|
for(i=1; i<=nlstate;i++) {
|
if (mobilav==1)
|
if (mobilav==1)
|
kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb][i][cptcod];
|
kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
|
else kk1=kk1+p3mat[i][j][h]*probs[(int)agedeb][i][cptcod];
|
else {
|
}
|
kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
|
if (h*hstepm/YEARM*stepm==cpt) fprintf(ficresf," %.5f ", kk1);
|
/* fprintf(ficresf," p3=%.3f p=%.3f ", p3mat[i][j][h], probs[(int)(agedeb)+1][i][cptcod]);*/
|
|
}
|
|
|
|
if (popforecast==1) kk2=kk1*popeffectif[(int)agedeb];
|
|
}
|
|
|
|
if (h==(int)(calagedate+12*cpt)){
|
|
fprintf(ficresf," %.3f", kk1);
|
|
|
|
if (popforecast==1) fprintf(ficresf," [%.f]", kk2);
|
|
}
|
}
|
}
|
}
|
}
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
|
}
|
}
|
|
}
|
}
|
}
|
}
|
}
|
if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
if (mobilav==1) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
|
|
if (popforecast==1) {
|
|
free_ivector(popage,0,AGESUP);
|
|
free_vector(popeffectif,0,AGESUP);
|
|
free_vector(popcount,0,AGESUP);
|
|
}
|
|
free_imatrix(s,1,maxwav+1,1,n);
|
|
free_vector(weight,1,n);
|
fclose(ficresf);
|
fclose(ficresf);
|
|
}/* End forecasting */
|
|
else{
|
|
erreur=108;
|
|
printf("Error %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d\n", erreur, stepm);
|
|
}
|
|
|
/*---------- Health expectancies and variances ------------*/
|
/*---------- Health expectancies and variances ------------*/
|
|
|
strcpy(filerest,"t");
|
strcpy(filerest,"t");
|
Line 2940 strcpy(fileresvpl,"vpl");
|
Line 3148 strcpy(fileresvpl,"vpl");
|
|
|
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
|
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
|
|
|
printf("End of Imach\n");
|
if(erreur >0)
|
|
printf("End of Imach with error %d\n",erreur);
|
|
else printf("End of Imach\n");
|
/* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */
|
/* gettimeofday(&end_time, (struct timezone*)0);*/ /* after time */
|
|
|
/* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/
|
/* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/
|
Line 2950 strcpy(fileresvpl,"vpl");
|
Line 3160 strcpy(fileresvpl,"vpl");
|
|
|
end:
|
end:
|
#ifdef windows
|
#ifdef windows
|
chdir(pathcd);
|
/* chdir(pathcd);*/
|
#endif
|
#endif
|
|
/*system("wgnuplot graph.plt");*/
|
system("..\\gp37mgw\\wgnuplot graph.plt");
|
/*system("../gp37mgw/wgnuplot graph.plt");*/
|
|
/*system("cd ../gp37mgw");*/
|
|
/* system("..\\gp37mgw\\wgnuplot graph.plt");*/
|
|
strcpy(plotcmd,GNUPLOTPROGRAM);
|
|
strcat(plotcmd," ");
|
|
strcat(plotcmd,optionfilegnuplot);
|
|
system(plotcmd);
|
|
|
#ifdef windows
|
#ifdef windows
|
while (z[0] != 'q') {
|
while (z[0] != 'q') {
|
chdir(pathcd);
|
chdir(path);
|
printf("\nType e to edit output files, c to start again, and q for exiting: ");
|
printf("\nType e to edit output files, c to start again, and q for exiting: ");
|
scanf("%s",z);
|
scanf("%s",z);
|
if (z[0] == 'c') system("./imach");
|
if (z[0] == 'c') system("./imach");
|