imach/src/imach.c - diff

Return to imach.c CVS log

Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.2 and 1.22

version 1.2, 2001/03/13 18:10:26	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 + bijage+ cijsex + 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 + bijage+ cijsex + 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
	#define GLOCK_ERROR_NOPATH -1 /* empty path */
	#define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */

#define MAXPARM 30 /* Maximum number of parameters for the optimization */	#define MAXPARM 30 /* Maximum number of parameters for the optimization */
#define NPARMAX 64 /* (nlstate+ndeath-1)nlstatencovmodel */	#define NPARMAX 64 /* (nlstate+ndeath-1)nlstatencovmodel */
Line 59	Line 70
#define NLSTATEMAX 8 /* Maximum number of live states (for func) */	#define NLSTATEMAX 8 /* Maximum number of live states (for func) */
#define NDEATHMAX 8 /* Maximum number of dead states (for func) */	#define NDEATHMAX 8 /* Maximum number of dead states (for func) */
#define NCOVMAX 8 /* Maximum number of covariates */	#define NCOVMAX 8 /* Maximum number of covariates */
#define MAXN 80000	#define MAXN 20000
#define YEARM 12. /* Number of months per year */	#define YEARM 12. /* Number of months per year */
#define AGESUP 130	#define AGESUP 130
#define AGEBASE 40	#define AGEBASE 40


	int erreur; /* Error number */
int nvar;	int nvar;
static int cptcov;	int cptcovn, cptcovage=0, cptcoveff=0,cptcov;
int cptcovn;
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 a121 +b12x ncovmodel=2 */	int ncovmodel, ncov; /* Total number of covariables including constant a121 +b12x ncovmodel=2 */
	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 */
	int jmin, jmax; /* min, max spacing between 2 waves */
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 */
	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 */
FILE fic,ficpar, ficparo,ficres, ficrespl, ficrespij, *ficrest;	FILE fic,ficpar, ficparo,ficres, ficrespl, ficrespij, ficrest,ficresf;
FILE ficgp, fichtm;	FILE ficgp, fichtm,ficresprob,ficpop;
FILE *ficreseij;	FILE *ficreseij;
char filerese[FILENAMELENGTH];	char filerese[FILENAMELENGTH];
FILE *ficresvij;	FILE *ficresvij;
Line 89 FILE *ficreseij;	Line 103 FILE *ficreseij;
FILE *ficresvpl;	FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];	char fileresvpl[FILENAMELENGTH];




#define NR_END 1	#define NR_END 1
#define FREE_ARG char*	#define FREE_ARG char*
#define FTOL 1.0e-10	#define FTOL 1.0e-10
Line 125 int stepm;	Line 136 int stepm;
/* Stepm, step in month: minimum step interpolation*/	/* Stepm, step in month: minimum step interpolation*/

int m,nb;	int m,nb;
int num, firstpass=0, lastpass=4,cod, *ncodemax;	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;	double pmmij, probs, **mobaverage;
	double dateintmean=0;

double *weight;	double *weight;
int *s; / Status */	int *s; / Status */
double agedc, *covar, idx;	double agedc, *covar, idx;
int *nbcode, Tcode, Tvar, *codtab;	int *nbcode, Tcode, Tvar, codtab, Tvard, Tprod, cptcovprod, *Tvaraff;

double ftol=FTOL; /* Tolerance for computing Max Likelihood */	double ftol=FTOL; /* Tolerance for computing Max Likelihood */
double ftolhess; /* Tolerance for computing hessian */	double ftolhess; /* Tolerance for computing hessian */

	/************** split ***********************/
	static int split( char path, char dirc, char name, char ext, char *finame )
	{
	char s; / pointer */
	int l1, l2; /* length counters */

	l1 = strlen( path ); /* length of path */
	if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
	#ifdef windows
	s = strrchr( path, '\\' ); /* find last / */
	#else
	s = strrchr( path, '/' ); /* find last / */
	#endif
	if ( s == NULL ) { /* no directory, so use current */
	#if defined(__bsd__) /* get current working directory */
	extern char *getwd( );

	if ( getwd( dirc ) == NULL ) {
	#else
	extern char *getcwd( );

	if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
	#endif
	return( GLOCK_ERROR_GETCWD );
	}
	strcpy( name, path ); /* we've got it */
	} else { /* strip direcotry from path */
	s++; /* after this, the filename */
	l2 = strlen( s ); /* length of filename */
	if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
	strcpy( name, s ); /* save file name */
	strncpy( dirc, path, l1 - l2 ); /* now the directory */
	dirc[l1-l2] = 0; /* add zero */
	}
	l1 = strlen( dirc ); /* length of directory */
	#ifdef windows
	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 */
	}


/******************************************/	/******************************************/

Line 166 int nbocc(char *s, char occ)	Line 228 int nbocc(char *s, char occ)

void cutv(char u,char v, char*t, char occ)	void cutv(char u,char v, char*t, char occ)
{	{
int i,lg,j,p;	int i,lg,j,p=0;
i=0;	i=0;
if (t[0]== occ) p=0;
for(j=0; j<=strlen(t)-1; j++) {	for(j=0; j<=strlen(t)-1; j++) {
if((t[j]!= occ) && (t[j+1]==occ)) p=j+1;	if((t[j]!= occ) && (t[j+1]== occ)) p=j+1;
}	}

lg=strlen(t);	lg=strlen(t);
for(j=0; j<p; j++) {	for(j=0; j<p; j++) {
(u[j] = t[j]);	(u[j] = t[j]);
u[p]='\0';
}	}
	u[p]='\0';

for(j=0; j<= lg; j++) {	for(j=0; j<= lg; j++) {
if (j>=(p+1))(v[j-p-1] = t[j]);	if (j>=(p+1))(v[j-p-1] = t[j]);
}	}
}	}


/******************** nrerror ******************/	/******************** nrerror ******************/

void nrerror(char error_text[])	void nrerror(char error_text[])
Line 615 double prevalim(double prlim, int nl	Line 675 double prevalim(double prlim, int nl
for (j=1;j<=nlstate+ndeath;j++){	for (j=1;j<=nlstate+ndeath;j++){
oldm[ii][j]=(ii==j ? 1.0 : 0.0);	oldm[ii][j]=(ii==j ? 1.0 : 0.0);
}	}
/* Even if hstepm = 1, at least one multiplication by the unit matrix */
	cov[1]=1.;

	/* Even if hstepm = 1, at least one multiplication by the unit matrix */
for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){	for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
newm=savm;	newm=savm;
/* Covariates have to be included here again */	/* Covariates have to be included here again */
cov[1]=1.;	cov[2]=agefin;
cov[2]=agefin;
if (cptcovn>0){	for (k=1; k<=cptcovn;k++) {
for (k=1; k<=cptcovn;k++) {cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];/printf("Tcode[ij]=%d nbcode=%d\n",Tcode[ij],nbcode[k][Tcode[ij]]);/}	cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
}	/printf("ij=%d Tvar[k]=%d nbcode=%d cov=%lf\n",ij, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k]);/
	}
	for (k=1; k<=cptcovage;k++)
	cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
	for (k=1; k<=cptcovprod;k++)
	cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];

	/printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);/
	/printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);/

out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);	out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);

savm=oldm;	savm=oldm;
Line 648 double prevalim(double prlim, int nl	Line 720 double prevalim(double prlim, int nl
}	}
}	}

/************* transition probabilities ********/	/************* transition probabilities *************/

double pmij(double ps, double cov, int ncovmodel, double x, int nlstate )	double pmij(double ps, double cov, int ncovmodel, double x, int nlstate )
{	{
Line 674 double pmij(double ps, double *cov,	Line 746 double pmij(double ps, double *cov,
ps[i][j]=s2;	ps[i][j]=s2;
}	}
}	}
	/ps[3][2]=1;/

for(i=1; i<= nlstate; i++){	for(i=1; i<= nlstate; i++){
s1=0;	s1=0;
for(j=1; j<i; j++)	for(j=1; j<i; j++)
Line 695 double pmij(double ps, double *cov,	Line 769 double pmij(double ps, double *cov,
}	}
}	}


/* for(ii=1; ii<= nlstate+ndeath; ii++){	/* for(ii=1; ii<= nlstate+ndeath; ii++){
for(jj=1; jj<= nlstate+ndeath; jj++){	for(jj=1; jj<= nlstate+ndeath; jj++){
printf("%lf ",ps[ii][jj]);	printf("%lf ",ps[ii][jj]);
Line 712 double pmij(double ps, double *cov,	Line 787 double pmij(double ps, double *cov,

double matprod2(double out, double in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double b)	double matprod2(double out, double in,long nrl, long nrh, long ncl, long nch, long ncolol, long ncoloh, double b)
{	{
/* Computes the matric product of in(1,nrh-nrl+1)(1,nch-ncl+1) times	/* Computes the matrix product of in(1,nrh-nrl+1)(1,nch-ncl+1) times
b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */	b(1,nch-ncl+1)(1,ncoloh-ncolol+1) into out(...) */
/* in, b, out are matrice of pointers which should have been initialized	/* in, b, out are matrice of pointers which should have been initialized
before: only the contents of out is modified. The function returns	before: only the contents of out is modified. The function returns
Line 758 double *hpxij(double *po, int nhstep	Line 833 double *hpxij(double *po, int nhstep
/* Covariates have to be included here again */	/* Covariates have to be included here again */
cov[1]=1.;	cov[1]=1.;
cov[2]=age+((h-1)hstepm + (d-1))stepm/YEARM;	cov[2]=age+((h-1)hstepm + (d-1))stepm/YEARM;
if (cptcovn>0){	for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];	for (k=1; k<=cptcovage;k++)
}	cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
	for (k=1; k<=cptcovprod;k++)
	cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];


/printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);/	/printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);/
/printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);/	/printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);/
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
Line 782 double *hpxij(double *po, int nhstep	Line 861 double *hpxij(double *po, int nhstep
/************* log-likelihood ***********/	/************* log-likelihood ***********/
double func( double *x)	double func( double *x)
{	{
int i, ii, j, k, mi, d;	int i, ii, j, k, mi, d, kk;
double l, ll[NLSTATEMAX], cov[NCOVMAX];	double l, ll[NLSTATEMAX], cov[NCOVMAX];
double **out;	double **out;
double sw; /* Sum of weights */	double sw; /* Sum of weights */
Line 792 double func( double *x)	Line 871 double func( double *x)
/* We are differentiating ll according to initial status */	/* We are differentiating ll according to initial status */
/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/	/* for (i=1;i<=npar;i++) printf("%f ", x[i]);*/
/*for(i=1;i<imx;i++)	/*for(i=1;i<imx;i++)
printf(" %d\n",s[4][i]);	printf(" %d\n",s[4][i]);
*/	*/
	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 (i=1,ipmx=0, sw=0.; i<=imx; i++){
for(mi=1; mi<= wav[i]-1; mi++){	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 (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);	for (j=1;j<=nlstate+ndeath;j++) oldm[ii][j]=(ii==j ? 1.0 : 0.0);
for(d=0; d<dh[mi][i]; d++){	for(d=0; d<dh[mi][i]; d++){
newm=savm;	newm=savm;
cov[1]=1.;	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;	for (kk=1; kk<=cptcovage;kk++) {
if (cptcovn>0){	cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
for (k=1; k<=cptcovn;k++) {	}
cov[2+k]=covar[Tvar[k]][i];
/* printf("k=%d cptcovn=%d %lf\n",k,cptcovn,covar[Tvar[k]][i]);*/	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
}	1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
}	savm=oldm;
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,	oldm=newm;
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
savm=oldm;
oldm=newm;
} /* end mult */	} /* 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]]);
/* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/	/* printf(" %f ",out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);*/
ipmx +=1;	ipmx +=1;
Line 827 printf(" %d\n",s[4][i]);	Line 907 printf(" %d\n",s[4][i]);
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= lipmx/sw; / To get the same order of magnitude as if weight=1 for every body */	l= lipmx/sw; / To get the same order of magnitude as if weight=1 for every body */

return -l;	return -l;
}	}

Line 847 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 864 void hesscov(double **matcov, double p[]	Line 943 void hesscov(double **matcov, double p[]
void lubksb(double *a, int npar, int indx, double b[]) ;	void lubksb(double *a, int npar, int indx, double b[]) ;
void ludcmp(double *a, int npar, int indx, double *d) ;	void ludcmp(double *a, int npar, int indx, double *d) ;


hess=matrix(1,npar,1,npar);	hess=matrix(1,npar,1,npar);

printf("\nCalculation of the hessian matrix. Wait...\n");	printf("\nCalculation of the hessian matrix. Wait...\n");
Line 872 void hesscov(double **matcov, double p[]	Line 950 void hesscov(double **matcov, double p[]
printf("%d",i);fflush(stdout);	printf("%d",i);fflush(stdout);
hess[i][i]=hessii(p,ftolhess,i,delti);	hess[i][i]=hessii(p,ftolhess,i,delti);
/printf(" %f ",p[i]);/	/printf(" %f ",p[i]);/
	/printf(" %lf ",hess[i][i]);/
}	}

for (i=1;i<=npar;i++) {	for (i=1;i<=npar;i++) {
for (j=1;j<=npar;j++) {	for (j=1;j<=npar;j++) {
if (j>i) {	if (j>i) {
printf(".%d%d",i,j);fflush(stdout);	printf(".%d%d",i,j);fflush(stdout);
hess[i][j]=hessij(p,delti,i,j);	hess[i][j]=hessij(p,delti,i,j);
hess[j][i]=hess[i][j];	hess[j][i]=hess[i][j];
	/printf(" %lf ",hess[i][j]);/
}	}
}	}
}	}
Line 983 double hessii( double x[], double delta,	Line 1063 double hessii( double x[], double delta,
}	}
}	}
delti[theta]=delts;	delti[theta]=delts;
return res;	return res;

}	}

double hessij( double x[], double delti[], int thetai,int thetaj)	double hessij( double x[], double delti[], int thetai,int thetaj)
Line 1095 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,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 1116 void freqsummary(char fileres[], int ag	Line 1197 void freqsummary(char fileres[], int ag
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);	freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
j1=0;	j1=0;

j=cptcovn;	j=cptcoveff;
if (cptcovn<1) {j=1;ncodemax[1]=1;}	if (cptcovn<1) {j=1;ncodemax[1]=1;}

for(k1=1; k1<=j;k1++){	for(k1=1; k1<=j;k1++){
for(i1=1; i1<=ncodemax[k1];i1++){	for(i1=1; i1<=ncodemax[k1];i1++){
j1++;	j1++;
	/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
	scanf("%d", i);*/
for (i=-1; i<=nlstate+ndeath; i++)	for (i=-1; i<=nlstate+ndeath; i++)
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) {
for (z1=1; z1<=cptcovn; z1++)	for (z1=1; z1<=cptcoveff; z1++)
if (covar[Tvar[z1]][i]!= nbcode[Tvar[z1]][codtab[j1][z1]]) bool=0;	if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]])
	bool=0;
}	}
if (bool==1) {	if (bool==1) {
for(m=firstpass; m<=lastpass-1; 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++;
	}

	}
}	}
}	}
}	}
if (cptcovn>0) {	if (cptcovn>0) {
fprintf(ficresp, "\n#Variable");	fprintf(ficresp, "\n#********** Variable ");
for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvar[z1],nbcode[Tvar[z1]][codtab[j1][z1]]);	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
}	fprintf(ficresp, "**********\n#");
fprintf(ficresp, "\n#");	}
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);	fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
fprintf(ficresp, "\n");	fprintf(ficresp, "\n");
Line 1159 void freqsummary(char fileres[], int ag	Line 1252 void freqsummary(char fileres[], int ag
printf("Age %d", i);	printf("Age %d", i);
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)	for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
pp[jk] += freq[jk][m][i];	pp[jk] += freq[jk][m][i];
}	}
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
for(m=-1, pos=0; m <=0 ; m++)	for(m=-1, pos=0; m <=0 ; m++)
Line 1169 void freqsummary(char fileres[], int ag	Line 1262 void freqsummary(char fileres[], int ag
else	else
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);	printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
}	}
for(jk=1; jk <=nlstate ; jk++){
for(m=1, pp[jk]=0; m <=nlstate+ndeath; m++)	for(jk=1; jk <=nlstate ; jk++){
	for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)
pp[jk] += freq[jk][m][i];	pp[jk] += freq[jk][m][i];
}	}

for(jk=1,pos=0; jk <=nlstate ; jk++)	for(jk=1,pos=0; jk <=nlstate ; jk++)
pos += pp[jk];	pos += pp[jk];
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
Line 1181 void freqsummary(char fileres[], int ag	Line 1276 void freqsummary(char fileres[], int ag
else	else
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);	printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
if( i <= (int) agemax){	if( i <= (int) agemax){
if(pos>=1.e-5)	if(pos>=1.e-5){
fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);	fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
	probs[i][jk][j1]= pp[jk]/pos;
	/printf("\ni=%d jk=%d j1=%d %.5f %.0f %.0f %f",i,jk,j1,pp[jk]/pos, pp[jk],pos,probs[i][jk][j1]);/
	}
else	else
fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);	fprintf(ficresp," %d NaNq %.0f %.0f",i,pp[jk],pos);
}	}
Line 1196 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 1215 void concatwav(int wav[], int **dh, int	Line 1397 void concatwav(int wav[], int **dh, int
*/	*/

int i, mi, m;	int i, mi, m;
int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;	/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
float sum=0.;	double sum=0., jmean=0.;*/

	int j, k=0,jk, ju, jl;
	double sum=0.;
	jmin=1e+5;
	jmax=-1;
	jmean=0.;
for(i=1; i<=imx; i++){	for(i=1; i<=imx; i++){
mi=0;	mi=0;
m=firstpass;	m=firstpass;
Line 1247 float sum=0.;	Line 1434 float sum=0.;
dh[mi][i]=1;	dh[mi][i]=1;
else{	else{
if (s[mw[mi+1][i]][i] > nlstate) {	if (s[mw[mi+1][i]][i] > nlstate) {
	if (agedc[i] < 2*AGESUP) {
j= rint(agedc[i]12-agev[mw[mi][i]][i]12);	j= rint(agedc[i]12-agev[mw[mi][i]][i]12);
if(j=0) j=1; /* Survives at least one month after exam */	if(j==0) j=1; /* Survives at least one month after exam */
	k=k+1;
	if (j >= jmax) jmax=j;
	if (j <= jmin) jmin=j;
	sum=sum+j;
	/* if (j<10) printf("j=%d num=%d ",j,i); */
	}
}	}
else{	else{
j= rint( (agev[mw[mi+1][i]][i]12 - agev[mw[mi][i]][i]12));	j= rint( (agev[mw[mi+1][i]][i]12 - agev[mw[mi][i]][i]12));
/printf("i=%d agevi+1=%lf agevi=%lf j=%d\n", i,agev[mw[mi+1][i]][i],agev[mw[mi][i]][i],j);/

k=k+1;	k=k+1;
if (j >= jmax) jmax=j;	if (j >= jmax) jmax=j;
else if (j <= jmin)jmin=j;	else if (j <= jmin)jmin=j;
	/* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
sum=sum+j;	sum=sum+j;
}	}
jk= j/stepm;	jk= j/stepm;
Line 1271 float sum=0.;	Line 1464 float sum=0.;
}	}
}	}
}	}
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,sum/k);	jmean=sum/k;
}	printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
	}
/********* Tricode **************************/	/********* Tricode **************************/
void tricode(int Tvar, int *nbcode, int imx)	void tricode(int Tvar, int *nbcode, int imx)
{	{
int Ndum[80],ij, k, j, i;	int Ndum[20],ij=1, k, j, i;
int cptcode=0;	int cptcode=0;
for (k=0; k<79; k++) Ndum[k]=0;	cptcoveff=0;

	for (k=0; k<19; 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; j++) {	for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
for (i=1; i<=imx; i++) {	for (i=1; i<=imx; i++) {
ij=(int)(covar[Tvar[j]][i]);	ij=(int)(covar[Tvar[j]][i]);
Ndum[ij]++;	Ndum[ij]++;
	/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;
}	}
/printf("cptcode=%d cptcovn=%d ",cptcode,cptcovn);/
for (i=0; i<=cptcode; i++) {	for (i=0; i<=cptcode; i++) {
if(Ndum[i]!=0) ncodemax[j]++;	if(Ndum[i]!=0) ncodemax[j]++;
}	}

ij=1;	ij=1;


for (i=1; i<=ncodemax[j]; i++) {	for (i=1; i<=ncodemax[j]; i++) {
for (k=0; k<=79; k++) {	for (k=0; k<=19; k++) {
if (Ndum[k] != 0) {	if (Ndum[k] != 0) {
nbcode[Tvar[j]][ij]=k;	nbcode[Tvar[j]][ij]=k;
ij++;	ij++;
Line 1302 void tricode(int Tvar, int *nbcode, in	Line 1500 void tricode(int Tvar, int *nbcode, in
if (ij > ncodemax[j]) break;	if (ij > ncodemax[j]) break;
}	}
}	}
}	}

}	for (k=0; k<19; k++) Ndum[k]=0;

	for (i=1; i<=ncovmodel-2; i++) {
	ij=Tvar[i];
	Ndum[ij]++;
	}

	ij=1;
	for (i=1; i<=10; i++) {
	if((Ndum[i]!=0) && (i<=ncov)){
	Tvaraff[ij]=i;
	ij++;
	}
	}

	cptcoveff=ij-1;
	}

/********* Health Expectancies **************/	/********* Health Expectancies **************/

Line 1365 void varevsij(char fileres[], double ***	Line 1579 void varevsij(char fileres[], double ***
double dnewm,doldm;	double dnewm,doldm;
int i, j, nhstepm, hstepm, h;	int i, j, nhstepm, hstepm, h;
int k, cptcode;	int k, cptcode;
double *xp;	double *xp;
double gp, gm;	double gp, gm;
double *gradg, *trgradg;	double *gradg, *trgradg;
double ***p3mat;	double ***p3mat;
Line 1401 void varevsij(char fileres[], double ***	Line 1615 void varevsij(char fileres[], double ***
}	}
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);

	if (popbased==1) {
	for(i=1; i<=nlstate;i++)
	prlim[i][i]=probs[(int)age][i][ij];
	}

for(j=1; j<= nlstate; j++){	for(j=1; j<= nlstate; j++){
for(h=0; h<=nhstepm; h++){	for(h=0; h<=nhstepm; h++){
for(i=1, gp[h][j]=0.;i<=nlstate;i++)	for(i=1, gp[h][j]=0.;i<=nlstate;i++)
Line 1412 void varevsij(char fileres[], double ***	Line 1632 void varevsij(char fileres[], double ***
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);

	if (popbased==1) {
	for(i=1; i<=nlstate;i++)
	prlim[i][i]=probs[(int)age][i][ij];
	}

for(j=1; j<= nlstate; j++){	for(j=1; j<= nlstate; j++){
for(h=0; h<=nhstepm; h++){	for(h=0; h<=nhstepm; h++){
for(i=1, gm[h][j]=0.;i<=nlstate;i++)	for(i=1, gm[h][j]=0.;i<=nlstate;i++)
gm[h][j] += prlim[i][i]*p3mat[i][j][h];	gm[h][j] += prlim[i][i]*p3mat[i][j][h];
}	}
}	}

for(j=1; j<= nlstate; j++)	for(j=1; j<= nlstate; j++)
for(h=0; h<=nhstepm; h++){	for(h=0; h<=nhstepm; h++){
gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];	gradg[h][theta][j]= (gp[h][j]-gm[h][j])/2./delti[theta];
Line 1547 void varprevlim(char fileres[], double *	Line 1774 void varprevlim(char fileres[], double *

}	}

	/********** Variance of one-step probabilities ****************/
	void varprob(char fileres[], double **matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij)
	{
	int i, j;
	int k=0, cptcode;
	double dnewm,doldm;
	double *xp;
	double gp, gm;
	double gradg, trgradg;
	double age,agelim, cov[NCOVMAX];
	int theta;
	char fileresprob[FILENAMELENGTH];

	strcpy(fileresprob,"prob");
	strcat(fileresprob,fileres);
	if((ficresprob=fopen(fileresprob,"w"))==NULL) {
	printf("Problem with resultfile: %s\n", fileresprob);
	}
	printf("Computing variance of one-step probabilities: result on file '%s' \n",fileresprob);


	xp=vector(1,npar);
	dnewm=matrix(1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
	doldm=matrix(1,(nlstate+ndeath)(nlstate+ndeath),1,(nlstate+ndeath)(nlstate+ndeath));

	cov[1]=1;
	for (age=bage; age<=fage; age ++){
	cov[2]=age;
	gradg=matrix(1,npar,1,9);
	trgradg=matrix(1,9,1,npar);
	gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));
	gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));

	for(theta=1; theta <=npar; theta++){
	for(i=1; i<=npar; i++)
	xp[i] = x[i] + (i==theta ?delti[theta]:0);

	pmij(pmmij,cov,ncovmodel,xp,nlstate);

	k=0;
	for(i=1; i<= (nlstate+ndeath); i++){
	for(j=1; j<=(nlstate+ndeath);j++){
	k=k+1;
	gp[k]=pmmij[i][j];
	}
	}

	for(i=1; i<=npar; i++)
	xp[i] = x[i] - (i==theta ?delti[theta]:0);


	pmij(pmmij,cov,ncovmodel,xp,nlstate);
	k=0;
	for(i=1; i<=(nlstate+ndeath); i++){
	for(j=1; j<=(nlstate+ndeath);j++){
	k=k+1;
	gm[k]=pmmij[i][j];
	}
	}

	for(i=1; i<= (nlstate+ndeath)*(nlstate+ndeath); i++)
	gradg[theta][i]=(gp[i]-gm[i])/2./delti[theta];
	}

	for(j=1; j<=(nlstate+ndeath)*(nlstate+ndeath);j++)
	for(theta=1; theta <=npar; theta++)
	trgradg[j][theta]=gradg[theta][j];

	matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);
	matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);

	pmij(pmmij,cov,ncovmodel,x,nlstate);

	k=0;
	for(i=1; i<=(nlstate+ndeath); i++){
	for(j=1; j<=(nlstate+ndeath);j++){
	k=k+1;
	gm[k]=pmmij[i][j];
	}
	}

	/*printf("\n%d ",(int)age);
	for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){


	printf("%e [%e ;%e] ",gm[i],gm[i]-2sqrt(doldm[i][i]),gm[i]+2sqrt(doldm[i][i]));
	}*/

	fprintf(ficresprob,"\n%d ",(int)age);

	for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){
	if (i== 2) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
	if (i== 4) fprintf(ficresprob,"%.3e %.3e ",gm[i],doldm[i][i]);
	}

	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);
	}
	free_vector(xp,1,npar);
	fclose(ficresprob);
	exit(0);
	}

/***********************************************/	/***********************************************/
/************** 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, aaa, cptcod;	int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;
double agedeb, agefin,hf;	double agedeb, agefin,hf;
double agemin=1.e20, agemax=-1.e20;	double agemin=1.e20, agemax=-1.e20;

Line 1569 int main()	Line 1898 int main()
int *indx;	int *indx;
char line[MAXLINE], linepar[MAXLINE];	char line[MAXLINE], linepar[MAXLINE];
char title[MAXLINE];	char title[MAXLINE];
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];	char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];
char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[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 */
int c, h , cpt,l;	int c, h , cpt,l;
int ju,jl, mi;	int ju,jl, mi;
int i1,j1, k1,jk,aa,bb, stepsize;	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,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;
double **prlim;	double **prlim;
Line 1594 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;
char version[80]="Imach version 62c, May 1999, INED-EUROREVES ";	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 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];


char z[1]="c", occ;	char z[1]="c", occ;
#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.63");	printf("\n%s",version);
printf("\nEnter the parameter file name: ");	if(argc <=1){
	printf("\nEnter the parameter file name: ");
	scanf("%s",pathtot);
	}
	else{
	strcpy(pathtot,argv[1]);
	}
	/if(getcwd(pathcd, 80)!= NULL)printf ("Error pathcd\n");/
	/*cygwin_split_path(pathtot,path,optionfile);
	printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
	/* cutv(path,optionfile,pathtot,'\\');*/

#ifdef windows	split(pathtot,path,optionfile,optionfilext,optionfilefiname);
scanf("%s",pathtot);	printf("pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
getcwd(pathcd, size);
cutv(path,optionfile,pathtot,'\\');
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 1650 int main()	Line 1999 int main()
fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf 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);	fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf 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);
printf("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);	printf("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);
fprintf(ficparo,"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);	fprintf(ficparo,"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);
	while((c=getc(ficpar))=='#' && c!= EOF){
covar=matrix(1,NCOVMAX,1,n);	ungetc(c,ficpar);
if (strlen(model)<=1) cptcovn=0;	fgets(line, MAXLINE, ficpar);
else {	puts(line);
j=0;	fputs(line,ficparo);
j=nbocc(model,'+');
cptcovn=j+1;
}	}
	ungetc(c,ficpar);


	covar=matrix(0,NCOVMAX,1,n);
	cptcovn=0;
	if (strlen(model)>1) cptcovn=nbocc(model,'+')+1;

ncovmodel=2+cptcovn;	ncovmodel=2+cptcovn;
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */	nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
Line 1688 int main()	Line 2041 int main()
fprintf(ficparo,"\n");	fprintf(ficparo,"\n");
}	}

npar= (nlstate+ndeath-1)nlstatencovmodel;	npar= (nlstate+ndeath-1)nlstatencovmodel;

p=param[1][1];	p=param[1][1];

/* Reads comments: lines beginning with '#' */	/* Reads comments: lines beginning with '#' */
Line 1749 int main()	Line 2103 int main()
printf("\n");	printf("\n");


if(mle==1){
/-------- data file ----------/	/-------- data file ----------/
if((ficres =fopen(fileres,"w"))==NULL) {	if((ficres =fopen(fileres,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileres);goto end;	printf("Problem with resultfile: %s\n", fileres);goto end;
Line 1777 int main()	Line 2130 int main()
s=imatrix(1,maxwav+1,1,n);	s=imatrix(1,maxwav+1,1,n);
adl=imatrix(1,maxwav+1,1,n);	adl=imatrix(1,maxwav+1,1,n);
tab=ivector(1,NCOVMAX);	tab=ivector(1,NCOVMAX);
ncodemax=ivector(1,NCOVMAX);	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)) {

Line 1801 int main()	Line 2154 int main()
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);

/* 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]));*/	/*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 %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),(covar[3][i]), (covar[4][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]));/

i=i+1;	i=i+1;
}	}
}	}
/scanf("%d",i);/	/* printf("ii=%d", ij);
	scanf("%d",i);*/
imx=i-1; /* Number of individuals */	imx=i-1; /* Number of individuals */

	/* for (i=1; i<=imx; i++){
	if ((s[1][i]==3) && (s[2][i]==2)) s[2][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;
	}

	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,8);	Tvar=ivector(1,15);
	Tprod=ivector(1,15);
	Tvaraff=ivector(1,15);
	Tvard=imatrix(1,15,1,2);
	Tage=ivector(1,15);

if (strlen(model) >1){	if (strlen(model) >1){
j=0;	j=0, j1=0, k1=1, k2=1;
j=nbocc(model,'+');	j=nbocc(model,'+');
	j1=nbocc(model,'*');
cptcovn=j+1;	cptcovn=j+1;
	cptcovprod=j1;


strcpy(modelsav,model);	strcpy(modelsav,model);
if (j==0) {	if ((strcmp(model,"age")==0) \|\| (strcmp(model,"age*age")==0)){
cutv(stra,strb,modelsav,'V'); Tvar[1]=atoi(strb);	printf("Error. Non available option model=%s ",model);
	goto end;
}	}
else {
for(i=j; i>=1;i--){	for(i=(j+1); i>=1;i--){
cutv(stra,strb,modelsav,'+');	cutv(stra,strb,modelsav,'+');
if (strchr(strb,'*')) {	if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav);
cutv(strd,strc,strb,'*');	/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
cutv(strb,stre,strc,'V');Tvar[i+1]=ncov+1;	/scanf("%d",i);/
cutv(strb,strc,strd,'V');	if (strchr(strb,'*')) {
for (k=1; k<=lastobs;k++)	cutv(strd,strc,strb,'*');
covar[ncov+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];	if (strcmp(strc,"age")==0) {
	cptcovprod--;
	cutv(strb,stre,strd,'V');
	Tvar[i]=atoi(stre);
	cptcovage++;
	Tage[cptcovage]=i;
	/printf("stre=%s ", stre);/
	}
	else if (strcmp(strd,"age")==0) {
	cptcovprod--;
	cutv(strb,stre,strc,'V');
	Tvar[i]=atoi(stre);
	cptcovage++;
	Tage[cptcovage]=i;
}	}
else {	else {
cutv(strd,strc,strb,'V');	cutv(strb,stre,strc,'V');
Tvar[i+1]=atoi(strc);	Tvar[i]=ncov+k1;
	cutv(strb,strc,strd,'V');
	Tprod[k1]=i;
	Tvard[k1][1]=atoi(strc);
	Tvard[k1][2]=atoi(stre);
	Tvar[cptcovn+k2]=Tvard[k1][1];
	Tvar[cptcovn+k2+1]=Tvard[k1][2];
	for (k=1; k<=lastobs;k++)
	covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
	k1++;
	k2=k2+2;
}	}
strcpy(modelsav,stra);
}	}
/cutv(strd,strc,stra,'V');/	else {
Tvar[1]=atoi(strc);	/printf("d=%s c=%s b=%s\n", strd,strc,strb);/
	/* scanf("%d",i);*/
	cutv(strd,strc,strb,'V');
	Tvar[i]=atoi(strc);
	}
	strcpy(modelsav,stra);
	/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
	scanf("%d",i);*/
}	}
}	}
/printf("tvar=%d ",Tvar[1]);/
/scanf("%d ",i);/	/*printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
	printf("cptcovprod=%d ", cptcovprod);
	scanf("%d ",i);*/
fclose(fic);	fclose(fic);

	/* 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(m=2; (m<= maxwav); m++)
	if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){
	anint[m][i]=9999;
	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 (mint[m][i]==99 \|\| anint[m][i]==9999) s[m][i]=-1;
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];
else{	else {
	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);
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;	agev[m][i]=1;
/* printf("i=%d m=%d agev=%lf \n",i,m, agev[m][i]); */
}
else if(agev[m][i] <agemin){	else if(agev[m][i] <agemin){
agemin=agev[m][i];	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);/	/printf(" Min anint[%d][%d]=%.2f annais[%d]=%.2f, agemin=%.2f\n",m,i,anint[m][i], i,annais[i], agemin);/
Line 1911 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 1925 printf("Total number of individuals= %d,	Line 2332 printf("Total number of individuals= %d,
concatwav(wav, dh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm);	concatwav(wav, dh, mw, s, agedc, agev, firstpass, lastpass, imx, nlstate, stepm);


Tcode=ivector(1,100);	Tcode=ivector(1,100);
nbcode=imatrix(1,nvar,1,8);	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);
h=0;	h=0;
m=pow(2,cptcovn);	m=pow(2,cptcoveff);

for(k=1;k<=cptcovn; 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,cptcovn+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;	if (h>m) h=1;codtab[h][k]=j;
}	}
Line 1945 Tcode=ivector(1,100);	Line 2352 Tcode=ivector(1,100);
}	}
}	}


/*for(i=1; i <=m ;i++){	/*for(i=1; i <=m ;i++){
for(k=1; k <=cptcovn; k++){	for(k=1; k <=cptcovn; k++){
printf("i=%d k=%d %d ",i,k,codtab[i][k]);	printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff);
}	}
printf("\n");	printf("\n");
}	}
scanf("%d",i);*/	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'. */
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 */
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]	/* 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] */	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) /	p=param[1][1]; /* (((param +1)+1)+0) /
/scanf("%d",i);/
mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);

	if(mle==1){
	mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
	}

/--------- 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 1992 Tcode=ivector(1,100);	Line 2403 Tcode=ivector(1,100);
}	}
}	}
}	}
	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\n");	fprintf(ficres,"# Scales\n");
printf("# Scales\n");	printf("# Scales\n");
for(i=1,jk=1; i <=nlstate; i++){	for(i=1,jk=1; i <=nlstate; i++){
Line 2012 Tcode=ivector(1,100);	Line 2424 Tcode=ivector(1,100);
fprintf(ficres,"\n");	fprintf(ficres,"\n");
}	}
}	}
}	}

k=1;	k=1;
fprintf(ficres,"# Covariance\n");	fprintf(ficres,"# Covariance\n");
Line 2047 Tcode=ivector(1,100);	Line 2459 Tcode=ivector(1,100);
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);
/------------ gnuplot -------------/	fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f\n",agemin,agemax,bage,fage);
chdir(pathcd);
if((ficgp=fopen("graph.plt","w"))==NULL) {	while((c=getc(ficpar))=='#' && c!= EOF){
printf("Problem with file graph.gp");goto end;	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.;
	dateprev2=anprev2+mprev2/12.+jprev2/365.;

	fscanf(ficpar,"pop_based=%d\n",&popbased);
	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,cptcovn);	m=pow(2,cptcoveff);

/* 1eme*/	/* 1eme*/
for (cpt=1; cpt<= nlstate ; cpt ++) {	for (cpt=1; cpt<= nlstate ; cpt ++) {
Line 2154 fprintf(ficgp,"\nset out \"v%s%d%d.gif\"	Line 2612 fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);	fprintf(ficgp,")) t\"prev(%d,%d)\" w l\n",cpt+1,cpt+1);
fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);	fprintf(ficgp,"set out \"p%s%d%d.gif\" \nreplot\n\n",strtok(optionfile, "."),cpt,k1);
}	}
}	}

/* proba elementaires */	/* proba elementaires */
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) {
/* fprintf(ficgp,"%1d%1d ",i,k);*/
for(j=1; j <=ncovmodel; j++){	for(j=1; j <=ncovmodel; j++){
fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);	/fprintf(ficgp,"%s%1d%1d=%f ",alph[j],i,k,p[jk]);/
	/fprintf(ficgp,"%s",alph[1]);/
	fprintf(ficgp,"p%d=%f ",jk,p[jk]);
jk++;	jk++;
fprintf(ficgp,"\n");	fprintf(ficgp,"\n");
}	}
}	}
}	}
}	}

for(jk=1; jk <=m; jk++) {	for(jk=1; jk <=m; jk++) {
fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot [%.f:%.f] ",agemin,agemax);	fprintf(ficgp,"\nset ter gif small size 400,300\nset log y\nplot [%.f:%.f] ",agemin,agemax);
for(i=1; i <=nlstate; i++) {	i=1;
for(k=1; k <=(nlstate+ndeath); k++){	for(k2=1; k2<=nlstate; k2++) {
if (k != i) {	k3=i;
fprintf(ficgp," exp(a%d%d+b%d%d*x",i,k,i,k);	for(k=1; k<=(nlstate+ndeath); k++) {
for(j=3; j <=ncovmodel; j++)	if (k != k2){
fprintf(ficgp,"+%s%d%d*%d",alph[j],i,k,nbcode[Tvar[j-2]][codtab[jk][j-2]]);	fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
fprintf(ficgp,")/(1");	ij=1;
for(k1=1; k1 <=(nlstate+ndeath); k1++)	for(j=3; j <=ncovmodel; j++) {
if (k1 != i) {	if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
fprintf(ficgp,"+exp(a%d%d+b%d%d*x",i,k1,i,k1);	fprintf(ficgp,"+p%d%dx",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
for(j=3; j <=ncovmodel; j++)	ij++;
fprintf(ficgp,"+%s%d%d*%d",alph[j],i,k,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
fprintf(ficgp,")");
}	}
fprintf(ficgp,") t \"p%d%d\" ", i,k);	else
if ((i+k)!= (nlstate*2+ndeath)) fprintf(ficgp,",");	fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
}	}
}	fprintf(ficgp,")/(1");
}
fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);	for(k1=1; k1 <=nlstate; k1++){
	fprintf(ficgp,"+exp(p%d+p%dx",k3+(k1-1)ncovmodel,k3+(k1-1)*ncovmodel+1);
	ij=1;
	for(j=3; j <=ncovmodel; j++){
	if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
	fprintf(ficgp,"+p%d%dx",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
	ij++;
	}
	else
	fprintf(ficgp,"+p%d%d",k3+(k1-1)ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
	}
	fprintf(ficgp,")");
	}
	fprintf(ficgp,") t \"p%d%d\" ", k2,k);
	if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
	i=i+ncovmodel;
	}
	}
	}
	fprintf(ficgp,"\nset out \"pe%s%d.gif\" \nreplot\n\n",strtok(optionfile, "."),jk);
}	}

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);
}	/* }*/

/________fin mle=1_________/	/________fin mle=1_________/

Line 2230 chdir(path);	Line 2704 chdir(path);
fprintf(ficparo,"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);
/--------- index.htm --------/	/--------- index.htm --------/

if((fichtm=fopen("index.htm","w"))==NULL) {	strcpy(optionfilehtm,optionfile);
printf("Problem with index.htm \n");goto end;	strcat(optionfilehtm,".htm");
	if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
	printf("Problem with %s \n",optionfilehtm);goto end;
}	}

fprintf(fichtm,"<body><ul> Imach, Version 0.63<hr> <li>Outputs files<br><br>\n	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>
	Total number of observations=%d <br>
	Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>
	<hr size=\"2\" color=\"#EC5E5E\">
	<li>Outputs files<br><br>\n
- Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n	- Observed prevalence in each state: <a href=\"p%s\">p%s</a> <br>\n
- Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>	- Estimated parameters and the covariance matrix: <a href=\"%s\">%s</a> <br>
- Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>	- Stationary prevalence in each state: <a href=\"pl%s\">pl%s</a> <br>
Line 2243 chdir(path);	Line 2724 chdir(path);
- Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>	- Life expectancies by age and initial health status: <a href=\"e%s\">e%s</a> <br>
- 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><br>",fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres,fileres);	- Standard deviation of stationary prevalences: <a href=\"vpl%s\">vpl%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);

fprintf(fichtm," <li>Graphs</li>\n<p>");	fprintf(fichtm," <li>Graphs</li><p>");

m=cptcovn;	m=cptcoveff;
if (cptcovn < 1) {m=1;ncodemax[1]=1;}	if (cptcovn < 1) {m=1;ncodemax[1]=1;}

j1=0;	j1=0;
Line 2255 chdir(path);	Line 2738 chdir(path);
for(i1=1; i1<=ncodemax[k1];i1++){	for(i1=1; i1<=ncodemax[k1];i1++){
j1++;	j1++;
if (cptcovn > 0) {	if (cptcovn > 0) {
fprintf(fichtm,"<hr>************ Results for covariates");	fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
for (cpt=1; cpt<=cptcovn;cpt++)	for (cpt=1; cpt<=cptcoveff;cpt++)
fprintf(fichtm," V%d=%d ",Tvar[cpt],nbcode[Tvar[cpt]][codtab[j1][cpt]]);	fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[j1][cpt]]);
fprintf(fichtm," ************\n<hr>");	fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}	}
fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>	fprintf(fichtm,"<br>- Probabilities: pe%s%d.gif<br>
<img src=\"pe%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);	<img src=\"pe%s%d.gif\">",strtok(optionfile, "."),j1,strtok(optionfile, "."),j1);
Line 2306 fclose(fichtm);	Line 2789 fclose(fichtm);
agebase=agemin;	agebase=agemin;
agelim=agemax;	agelim=agemax;
ftolpl=1.e-10;	ftolpl=1.e-10;
i1=cptcovn;	i1=cptcoveff;
if (cptcovn < 1){i1=1;}	if (cptcovn < 1){i1=1;}

for(cptcov=1;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;
/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#******");
for(j=1;j<=cptcovn;j++)	for(j=1;j<=cptcoveff;j++)
fprintf(ficrespl,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);	fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrespl,"******\n");	fprintf(ficrespl,"******\n");

for (age=agebase; age<=agelim; age++){	for (age=agebase; age<=agelim; age++){
Line 2328 fclose(fichtm);	Line 2811 fclose(fichtm);
}	}
}	}
fclose(ficrespl);	fclose(ficrespl);

/------------- h Pij x at various ages ------------/	/------------- h Pij x at various ages ------------/

strcpy(filerespij,"pij"); strcat(filerespij,fileres);	strcpy(filerespij,"pij"); strcat(filerespij,fileres);
Line 2337 fclose(fichtm);	Line 2821 fclose(fichtm);
printf("Computing pij: result on file '%s' \n", filerespij);	printf("Computing pij: result on file '%s' \n", filerespij);

stepsize=(int) (stepm+YEARM-1)/YEARM;	stepsize=(int) (stepm+YEARM-1)/YEARM;
if (stepm<=24) stepsize=2;	/if (stepm<=24) stepsize=2;/

agelim=AGESUP;	agelim=AGESUP;
hstepm=stepsizeYEARM; / Every year of age */	hstepm=stepsizeYEARM; / Every year of age */
Line 2348 fclose(fichtm);	Line 2832 fclose(fichtm);
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<=cptcovn;j++)	for(j=1;j<=cptcoveff;j++)
fprintf(ficrespij,"V%d=%d ",Tvar[j],nbcode[Tvar[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 */
Line 2376 fclose(fichtm);	Line 2860 fclose(fichtm);
}	}
}	}

	/* varprob(fileres, matcov, p, delti, nlstate, (int) bage, (int) fage,k);*/

fclose(ficrespij);	fclose(ficrespij);

	if(stepm == 1) {
	/---------- 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");
	strcat(fileresf,fileres);
	if((ficresf=fopen(fileresf,"w"))==NULL) {
	printf("Problem with forecast resultfile: %s\n", fileresf);goto end;
	}
	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 */

	if (cptcoveff==0) ncodemax[cptcoveff]=1;

	if (mobilav==1) {
	mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
	for (agedeb=bage+3; agedeb<=fage-2; agedeb++)
	for (i=1; i<=nlstate;i++)
	for (cptcod=1;cptcod<=ncodemax[cptcov];cptcod++)
	mobaverage[(int)agedeb][i][cptcod]=0.;

	for (agedeb=bage+4; agedeb<=fage; agedeb++){
	for (i=1; i<=nlstate;i++){
	for (cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
	for (cpt=0;cpt<=4;cpt++){
	mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]+probs[(int)agedeb-cpt][i][cptcod];
	}
	mobaverage[(int)agedeb-2][i][cptcod]=mobaverage[(int)agedeb-2][i][cptcod]/5;
	}
	}
	}
	}

	stepsize=(int) (stepm+YEARM-1)/YEARM;
	if (stepm<=12) stepsize=1;

	agelim=AGESUP;
	/hstepm=stepsizeYEARM; // Every year of age */
	hstepm=1;
	hstepm=hstepm/stepm; /* Typically 2 years, = 2 years/6 months = 4 */
	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(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
	k=k+1;
	fprintf(ficresf,"\n#******");
	for(j=1;j<=cptcoveff;j++) {
	fprintf(ficresf," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
	}
	fprintf(ficresf,"******\n");
	fprintf(ficresf,"# StartingAge FinalAge");
	for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);
	if (popforecast==1) fprintf(ficresf," [Population]");

	for (cpt=0; cpt<4;cpt++) {
	fprintf(ficresf,"\n");
	fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);

	for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(bage-((int)calagedate %12)/12.); agedeb--){ /* If stepm=6 months */
	nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
	nhstepm = nhstepm/hstepm;
	/printf("agedeb=%.lf stepm=%d hstepm=%d nhstepm=%d \n",agedeb,stepm,hstepm,nhstepm);/

	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
	oldm=oldms;savm=savms;
	hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);

	for (h=0; h<=nhstepm; h++){
	if (h==(int) (calagedate+YEARM*cpt)) {
	fprintf(ficresf,"\n %.f ",agedeb+hhstepm/YEARMstepm);
	}
	for(j=1; j<=nlstate+ndeath;j++) {
	kk1=0.;kk2=0;
	for(i=1; i<=nlstate;i++) {
	if (mobilav==1)
	kk1=kk1+p3mat[i][j][h]*mobaverage[(int)agedeb+1][i][cptcod];
	else {
	kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
	/* 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);
	}
	}
	}
	}
	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);
	}/* 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 2407 fclose(fichtm);	Line 3040 fclose(fichtm);
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#****** ");
for(j=1;j<=cptcovn;j++)	for(j=1;j<=cptcoveff;j++)
fprintf(ficrest,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);	fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrest,"******\n");	fprintf(ficrest,"******\n");

fprintf(ficreseij,"\n#****** ");	fprintf(ficreseij,"\n#****** ");
for(j=1;j<=cptcovn;j++)	for(j=1;j<=cptcoveff;j++)
fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);	fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
fprintf(ficreseij,"******\n");	fprintf(ficreseij,"******\n");

fprintf(ficresvij,"\n#****** ");	fprintf(ficresvij,"\n#****** ");
for(j=1;j<=cptcovn;j++)	for(j=1;j<=cptcoveff;j++)
fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);	fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
fprintf(ficresvij,"******\n");	fprintf(ficresvij,"******\n");

Line 2437 fclose(fichtm);	Line 3070 fclose(fichtm);
epj=vector(1,nlstate+1);	epj=vector(1,nlstate+1);
for(age=bage; age <=fage ;age++){	for(age=bage; age <=fage ;age++){
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);	prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
	if (popbased==1) {
	for(i=1; i<=nlstate;i++)
	prlim[i][i]=probs[(int)age][i][k];
	}

fprintf(ficrest," %.0f",age);	fprintf(ficrest," %.0f",age);
for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){	for(j=1, epj[nlstate+1]=0.;j <=nlstate;j++){
for(i=1, epj[j]=0.;i <=nlstate;i++) {	for(i=1, epj[j]=0.;i <=nlstate;i++) {
Line 2456 fclose(fichtm);	Line 3094 fclose(fichtm);
}	}
}	}




fclose(ficreseij);	fclose(ficreseij);
fclose(ficresvij);	fclose(ficresvij);
fclose(ficrest);	fclose(ficrest);
Line 2478 strcpy(fileresvpl,"vpl");	Line 3119 strcpy(fileresvpl,"vpl");
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#****** ");
for(j=1;j<=cptcovn;j++)	for(j=1;j<=cptcoveff;j++)
fprintf(ficresvpl,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);	fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficresvpl,"******\n");	fprintf(ficresvpl,"******\n");

varpl=matrix(1,nlstate,(int) bage, (int) fage);	varpl=matrix(1,nlstate,(int) bage, (int) fage);
Line 2501 strcpy(fileresvpl,"vpl");	Line 3142 strcpy(fileresvpl,"vpl");
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(matcov,1,npar,1,npar);	free_matrix(matcov,1,npar,1,npar);
free_vector(delti,1,npar);	free_vector(delti,1,npar);

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);*/
/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
system("wgnuplot ../gp37mgw/graph.plt");	/system("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");
else if (z[0] == 'e') {	else if (z[0] == 'e') {
chdir(path);	chdir(path);
system("index.htm");	system(optionfilehtm);
}	}
else if (z[0] == 'q') exit(0);	else if (z[0] == 'q') exit(0);
}	}

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>

Removed from v.1.2
changed lines
	Added in v.1.22