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.59 and 1.74

version 1.59, 2002/11/18 23:01:13	version 1.74, 2003/05/02 18:51:41
Line 32	Line 32
hPijx is the probability to be observed in state i at age x+h	hPijx is the probability to be observed in state i at age x+h
conditional to the observed state i at age x. The delay 'h' can be	conditional to the observed state i at age x. The delay 'h' can be
split into an exact number (nh*stepm) of unobserved intermediate	split into an exact number (nh*stepm) of unobserved intermediate
states. This elementary transition (by month or quarter trimester,	states. This elementary transition (by month, quarter,
semester or year) is model as a multinomial logistic. The hPx	semester or year) is modelled as a multinomial logistic. The hPx
matrix is simply the matrix product of nh*stepm elementary matrices	matrix is simply the matrix product of nh*stepm elementary matrices
and the contribution of each individual to the likelihood is simply	and the contribution of each individual to the likelihood is simply
hPijx.	hPijx.
Line 48	Line 48
It is copyrighted identically to a GNU software product, ie programme and	It is copyrighted identically to a GNU software product, ie programme and
software can be distributed freely for non commercial use. Latest version	software can be distributed freely for non commercial use. Latest version
can be accessed at http://euroreves.ined.fr/imach .	can be accessed at http://euroreves.ined.fr/imach .

	Help to debug: LD_PRELOAD=/usr/local/lib/libnjamd.so ./imach foo.imach
	or better on gdb : set env LD_PRELOAD=/usr/local/lib/libnjamd.so

**********************************************************************/	**********************************************************************/
	/*
	main
	read parameterfile
	read datafile
	concatwav
	if (mle >= 1)
	mlikeli
	print results files
	if mle==1
	computes hessian
	read end of parameter file: agemin, agemax, bage, fage, estepm
	begin-prev-date,...
	open gnuplot file
	open html file
	stable prevalence
	for age prevalim()
	h Pij x
	variance of p varprob
	forecasting if prevfcast==1 prevforecast call prevalence()
	health expectancies
	Variance-covariance of DFLE
	prevalence()
	movingaverage()
	varevsij()
	if popbased==1 varevsij(,popbased)
	total life expectancies
	Variance of stable prevalence
	end
	*/




#include <math.h>	#include <math.h>
#include <stdio.h>	#include <stdio.h>
Line 83	Line 119
#define ODIRSEPARATOR '\\'	#define ODIRSEPARATOR '\\'
#endif	#endif

char version[80]="Imach version 0.9, November 2002, INED-EUROREVES ";	char version[80]="Imach version 0.95, February 2003, INED-EUROREVES ";
int erreur; /* Error number */	int erreur; /* Error number */
int nvar;	int nvar;
int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;	int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
Line 188 static int split( char path, char dirc	Line 224 static int split( char path, char dirc
if ( ss == NULL ) { /* no directory, so use current */	if ( ss == NULL ) { /* no directory, so use current */
/*if(strrchr(path, ODIRSEPARATOR )==NULL)	/*if(strrchr(path, ODIRSEPARATOR )==NULL)
printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/	printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
#if defined(__bsd__) /* get current working directory */	/* get current working directory */
extern char *getwd( );	/* extern char* getcwd ( char buf , int len);/

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

if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {	if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
#endif
return( GLOCK_ERROR_GETCWD );	return( GLOCK_ERROR_GETCWD );
}	}
strcpy( name, path ); /* we've got it */	strcpy( name, path ); /* we've got it */
Line 365 double **matrix(long nrl, long nrh, long	Line 395 double **matrix(long nrl, long nrh, long

for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;	for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
return m;	return m;
	/* print ((m+1)+70) ou print m[1][70]; print m+1 or print &(m[1])
	*/
}	}

/***********************free matrix **********************/	/***********************free matrix **********************/
Line 404 double ***ma3x(long nrl, long nrh, long	Line 436 double ***ma3x(long nrl, long nrh, long
for (j=ncl+1; j<=nch; j++)	for (j=ncl+1; j<=nch; j++)
m[i][j]=m[i][j-1]+nlay;	m[i][j]=m[i][j-1]+nlay;
}	}
return m;	return m;
	/* gdb: p *(m+1) <=> p m[1] and p (m+1) <=> p (m+1) <=> p &(m[1])
	&(m[i][j][k]) <=> (((m+i) + j)+k)
	*/
}	}

/***********************free ma3x **********************/	/***********************free ma3x **********************/
Line 856 double matprod2(double out, double *	Line 891 double matprod2(double out, double *

double *hpxij(double po, int nhstepm, double age, int hstepm, double x, int nlstate, int stepm, double oldm, double savm, int ij )	double *hpxij(double po, int nhstepm, double age, int hstepm, double x, int nlstate, int stepm, double oldm, double savm, int ij )
{	{
/* Computes the transition matrix starting at age 'age' over 'nhstepmhstepmstepm' month	/* Computes the transition matrix starting at age 'age' over
duration (i.e. until	'nhstepmhstepmstepm' months (i.e. until
age (in years) age+nhstepmstepm/12) by multiplying nhstepmhstepm matrices.	age (in years) age+nhstepmhstepmstepm/12) by multiplying
	nhstepm*hstepm matrices.
Output is stored in matrix po[i][j][h] for h every 'hstepm' step	Output is stored in matrix po[i][j][h] for h every 'hstepm' step
(typically every 2 years instead of every month which is too big).	(typically every 2 years instead of every month which is too big
	for the memory).
Model is determined by parameters x and covariates have to be	Model is determined by parameters x and covariates have to be
included manually here.	included manually here.

Line 917 double func( double *x)	Line 954 double func( double *x)
double sw; /* Sum of weights */	double sw; /* Sum of weights */
double lli; /* Individual log likelihood */	double lli; /* Individual log likelihood */
int s1, s2;	int s1, s2;
double bbh;	double bbh, survp;
long ipmx;	long ipmx;
/extern weight /	/extern weight /
/* We are differentiating ll according to initial status */	/* We are differentiating ll according to initial status */
Line 928 double func( double *x)	Line 965 double func( double *x)
cov[1]=1.;	cov[1]=1.;

for(k=1; k<=nlstate; k++) ll[k]=0.;	for(k=1; k<=nlstate; k++) ll[k]=0.;
for (i=1,ipmx=0, sw=0.; i<=imx; i++){
for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];	if(mle==1){
for(mi=1; mi<= wav[i]-1; mi++){	for (i=1,ipmx=0, sw=0.; i<=imx; i++){
for (ii=1;ii<=nlstate+ndeath;ii++)	for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
for (j=1;j<=nlstate+ndeath;j++){	for(mi=1; mi<= wav[i]-1; mi++){
oldm[ii][j]=(ii==j ? 1.0 : 0.0);	for (ii=1;ii<=nlstate+ndeath;ii++)
savm[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++){	savm[ii][j]=(ii==j ? 1.0 : 0.0);
newm=savm;	}
cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;	for(d=0; d<dh[mi][i]; d++){
for (kk=1; kk<=cptcovage;kk++) {	newm=savm;
cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
}	for (kk=1; kk<=cptcovage;kk++) {
	cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,	}
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
savm=oldm;	1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
oldm=newm;	savm=oldm;
	oldm=newm;
	} /* end mult */

	/lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);/ /* Original formula */
	/* But now since version 0.9 we anticipate for bias and large stepm.
	* If stepm is larger than one month (smallest stepm) and if the exact delay
	* (in months) between two waves is not a multiple of stepm, we rounded to
	* the nearest (and in case of equal distance, to the lowest) interval but now
	* we keep into memory the bias bh[mi][i] and also the previous matrix product
	* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
	* probability in order to take into account the bias as a fraction of the way
	* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
	* -stepm/2 to stepm/2 .
	* For stepm=1 the results are the same as for previous versions of Imach.
	* For stepm > 1 the results are less biased than in previous versions.
	*/
	s1=s[mw[mi][i]][i];
	s2=s[mw[mi+1][i]][i];
	bbh=(double)bh[mi][i]/(double)stepm;
	/* bias is positive if real duration
	* is higher than the multiple of stepm and negative otherwise.
	*/
	/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.+bbh)out[s1][s2]));/
	if( s2 > nlstate){
	/* i.e. if s2 is a death state and if the date of death is known then the contribution
	to the likelihood is the probability to die between last step unit time and current
	step unit time, which is also the differences between probability to die before dh
	and probability to die before dh-stepm .
	In version up to 0.92 likelihood was computed
	as if date of death was unknown. Death was treated as any other
	health state: the date of the interview describes the actual state
	and not the date of a change in health state. The former idea was
	to consider that at each interview the state was recorded
	(healthy, disable or death) and IMaCh was corrected; but when we
	introduced the exact date of death then we should have modified
	the contribution of an exact death to the likelihood. This new
	contribution is smaller and very dependent of the step unit
	stepm. It is no more the probability to die between last interview
	and month of death but the probability to survive from last
	interview up to one month before death multiplied by the
	probability to die within a month. Thanks to Chris
	Jackson for correcting this bug. Former versions increased
	mortality artificially. The bad side is that we add another loop
	which slows down the processing. The difference can be up to 10%
	lower mortality.
	*/
	lli=log(out[s1][s2] - savm[s1][s2]);
	}else{
	lli= log((1.+bbh)out[s1][s2]- bbhsavm[s1][s2]); /* linear interpolation */
	/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)out[s1][s2]- bbh(savm[s1][s2])):log((1.+bbh)out[s1][s2]));/ /* linear interpolation */
	}
	/lli=(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]);/
	/if(lli ==000.0)/
	/printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); /
	ipmx +=1;
	sw += weight[i];
	ll[s[mw[mi][i]][i]] += 2weight[i]lli;
	} /* end of wave */
	} /* end of individual */
	} else if(mle==2){
	for (i=1,ipmx=0, sw=0.; i<=imx; i++){
	for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
	for(mi=1; mi<= wav[i]-1; mi++){
	for (ii=1;ii<=nlstate+ndeath;ii++)
	for (j=1;j<=nlstate+ndeath;j++){
	oldm[ii][j]=(ii==j ? 1.0 : 0.0);
	savm[ii][j]=(ii==j ? 1.0 : 0.0);
	}
	for(d=0; d<=dh[mi][i]; d++){
	newm=savm;
	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
	for (kk=1; kk<=cptcovage;kk++) {
	cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
	}
	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
	1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
	savm=oldm;
	oldm=newm;
	} /* end mult */

	/lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);/ /* Original formula */
	/* But now since version 0.9 we anticipate for bias and large stepm.
	* If stepm is larger than one month (smallest stepm) and if the exact delay
	* (in months) between two waves is not a multiple of stepm, we rounded to
	* the nearest (and in case of equal distance, to the lowest) interval but now
	* we keep into memory the bias bh[mi][i] and also the previous matrix product
	* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
	* probability in order to take into account the bias as a fraction of the way
	* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
	* -stepm/2 to stepm/2 .
	* For stepm=1 the results are the same as for previous versions of Imach.
	* For stepm > 1 the results are less biased than in previous versions.
	*/
	s1=s[mw[mi][i]][i];
	s2=s[mw[mi+1][i]][i];
	bbh=(double)bh[mi][i]/(double)stepm;
	/* bias is positive if real duration
	* is higher than the multiple of stepm and negative otherwise.
	*/
	lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)out[s1][s2]- bbh(savm[s1][s2])):log((1.+bbh)out[s1][s2])); / linear interpolation */
	/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.+bbh)out[s1][s2]));/
	/lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.-+bh)out[s1][s2])); / / exponential interpolation */
	/lli=(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]);/
	/if(lli ==000.0)/
	/printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); /
	ipmx +=1;
	sw += weight[i];
	ll[s[mw[mi][i]][i]] += 2weight[i]lli;
	} /* end of wave */
	} /* end of individual */
	} else if(mle==3){ /* exponential inter-extrapolation */
	for (i=1,ipmx=0, sw=0.; i<=imx; i++){
	for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
	for(mi=1; mi<= wav[i]-1; mi++){
	for (ii=1;ii<=nlstate+ndeath;ii++)
	for (j=1;j<=nlstate+ndeath;j++){
	oldm[ii][j]=(ii==j ? 1.0 : 0.0);
	savm[ii][j]=(ii==j ? 1.0 : 0.0);
	}
	for(d=0; d<dh[mi][i]; d++){
	newm=savm;
	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
	for (kk=1; kk<=cptcovage;kk++) {
	cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
	}
	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
	1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
	savm=oldm;
	oldm=newm;
	} /* end mult */

	/lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);/ /* Original formula */
	/* But now since version 0.9 we anticipate for bias and large stepm.
	* If stepm is larger than one month (smallest stepm) and if the exact delay
	* (in months) between two waves is not a multiple of stepm, we rounded to
	* the nearest (and in case of equal distance, to the lowest) interval but now
	* we keep into memory the bias bh[mi][i] and also the previous matrix product
	* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the
	* probability in order to take into account the bias as a fraction of the way
	* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
	* -stepm/2 to stepm/2 .
	* For stepm=1 the results are the same as for previous versions of Imach.
	* For stepm > 1 the results are less biased than in previous versions.
	*/
	s1=s[mw[mi][i]][i];
	s2=s[mw[mi+1][i]][i];
	bbh=(double)bh[mi][i]/(double)stepm;
	/* bias is positive if real duration
	* is higher than the multiple of stepm and negative otherwise.
	*/
	/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)out[s1][s2]- bbh(savm[s1][s2])):log((1.+bbh)out[s1][s2])); / /* linear interpolation */
	lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.+bbh)out[s1][s2])); / exponential inter-extrapolation */
	/lli=(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]);/
	/if(lli ==000.0)/
	/printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); /
	ipmx +=1;
	sw += weight[i];
	ll[s[mw[mi][i]][i]] += 2weight[i]lli;
	} /* end of wave */
	} /* end of individual */
	}else{ /* ml=4 no inter-extrapolation */
	for (i=1,ipmx=0, sw=0.; i<=imx; i++){
	for (k=1; k<=cptcovn;k++) cov[2+k]=covar[Tvar[k]][i];
	for(mi=1; mi<= wav[i]-1; mi++){
	for (ii=1;ii<=nlstate+ndeath;ii++)
	for (j=1;j<=nlstate+ndeath;j++){
	oldm[ii][j]=(ii==j ? 1.0 : 0.0);
	savm[ii][j]=(ii==j ? 1.0 : 0.0);
	}
	for(d=0; d<dh[mi][i]; d++){
	newm=savm;
	cov[2]=agev[mw[mi][i]][i]+d*stepm/YEARM;
	for (kk=1; kk<=cptcovage;kk++) {
	cov[Tage[kk]+2]=covar[Tvar[Tage[kk]]][i]*cov[2];
	}

} /* end mult */	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
	1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
	savm=oldm;
	oldm=newm;
	} /* end mult */

/lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);/ /* Original formula */	lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
/* But now since version 0.9 we anticipate for bias and large stepm.	ipmx +=1;
* If stepm is larger than one month (smallest stepm) and if the exact delay	sw += weight[i];
* (in months) between two waves is not a multiple of stepm, we rounded to	ll[s[mw[mi][i]][i]] += 2weight[i]lli;
* the nearest (and in case of equal distance, to the lowest) interval but now	} /* end of wave */
* we keep into memory the bias bh[mi][i] and also the previous matrix product	} /* end of individual */
* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the	} /* End of if */
* probability in order to take into account the bias as a fraction of the way
* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies
* -stepm/2 to stepm/2 .
* For stepm=1 the results are the same as for previous versions of Imach.
* For stepm > 1 the results are less biased than in previous versions.
*/
s1=s[mw[mi][i]][i];
s2=s[mw[mi+1][i]][i];
bbh=(double)bh[mi][i]/(double)stepm;
lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.-bbh)*out[s1][s2]));
/lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.-bbh)out[s1][s2]));*/
/lli=(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]);/
/if(lli ==000.0)/
/printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); /
ipmx +=1;
sw += weight[i];
ll[s[mw[mi][i]][i]] += 2weight[i]lli;
} /* end of wave */
} /* end of individual */

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 */
Line 990 double func( double *x)	Line 1185 double func( double *x)
void mlikeli(FILE ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (func)(double []))	void mlikeli(FILE ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (func)(double []))
{	{
int i,j, iter;	int i,j, iter;
double *xi,delti;	double **xi;
double fret;	double fret;
xi=matrix(1,npar,1,npar);	xi=matrix(1,npar,1,npar);
for (i=1;i<=npar;i++)	for (i=1;i<=npar;i++)
Line 1000 void mlikeli(FILE *ficres,double p[], in	Line 1195 void mlikeli(FILE *ficres,double p[], in
powell(p,xi,npar,ftol,&iter,&fret,func);	powell(p,xi,npar,ftol,&iter,&fret,func);

printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));	printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
fprintf(ficlog,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));	fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));	fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));

}	}
Line 1263 void lubksb(double *a, int n, int indx	Line 1458 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 Tvaraff, int nbcode, int ncodemax,double mint,double anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)	void freqsummary(char fileres[], int iagemin, int iagemax, int s, double agev, int nlstate, int imx, int Tvaraff, int nbcode, int ncodemax,double mint,double anint, double dateprev1,double dateprev2,double jprev1, double mprev1,double anprev1,double jprev2, double mprev2,double anprev2)
{ /* Some frequencies */	{ /* Some frequencies */

int i, m, jk, k1,i1, j1, bool, z1,z2,j;	int i, m, jk, k1,i1, j1, bool, z1,z2,j;
int first;	int first;
double **freq; / Frequencies */	double **freq; / Frequencies */
double *pp;	double pp, *prop;
double pos, k2, dateintsum=0,k2cpt=0;	double pos,posprop, 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);	prop=matrix(1,nlstate,iagemin,iagemax+3);
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 1283 void freqsummary(char fileres[], int ag	Line 1478 void freqsummary(char fileres[], int ag
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);	fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
exit(0);	exit(0);
}	}
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);	freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);
j1=0;	j1=0;

j=cptcoveff;	j=cptcoveff;
Line 1298 void freqsummary(char fileres[], int ag	Line 1493 void freqsummary(char fileres[], int ag
scanf("%d", i);*/	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=iagemin; m <= iagemax+3; m++)
freq[i][jk][m]=0;	freq[i][jk][m]=0;

	for (i=1; i<=nlstate; i++)
	for(m=iagemin; m <= iagemax+3; m++)
	prop[i][m]=0;

dateintsum=0;	dateintsum=0;
k2cpt=0;	k2cpt=0;
Line 1314 void freqsummary(char fileres[], int ag	Line 1513 void freqsummary(char fileres[], int ag
for(m=firstpass; m<=lastpass; m++){	for(m=firstpass; m<=lastpass; m++){
k2=anint[m][i]+(mint[m][i]/12.);	k2=anint[m][i]+(mint[m][i]/12.);
if ((k2>=dateprev1) && (k2<=dateprev2)) {	if ((k2>=dateprev1) && (k2<=dateprev2)) {
if(agev[m][i]==0) agev[m][i]=agemax+1;	if(agev[m][i]==0) agev[m][i]=iagemax+1;
if(agev[m][i]==1) agev[m][i]=agemax+2;	if(agev[m][i]==1) agev[m][i]=iagemax+2;
	if (s[m][i]>0 && s[m][i]<=nlstate) prop[s[m][i]][(int)agev[m][i]] += weight[i];
if (m<lastpass) {	if (m<lastpass) {
freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];	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];	freq[s[m][i]][s[m+1][i]][iagemax+3] += weight[i];
}	}

if ((agev[m][i]>1) && (agev[m][i]< (agemax+3))) {	if ((agev[m][i]>1) && (agev[m][i]< (iagemax+3))) {
dateintsum=dateintsum+k2;	dateintsum=dateintsum+k2;
k2cpt++;	k2cpt++;
}	}
Line 1341 void freqsummary(char fileres[], int ag	Line 1541 void freqsummary(char fileres[], int ag
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");

for(i=(int)agemin; i <= (int)agemax+3; i++){	for(i=iagemin; i <= iagemax+3; i++){
if(i==(int)agemax+3){	if(i==iagemax+3){
fprintf(ficlog,"Total");	fprintf(ficlog,"Total");
}else{	}else{
if(first==1){	if(first==1){
Line 1373 void freqsummary(char fileres[], int ag	Line 1573 void freqsummary(char fileres[], int ag
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)	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,posprop=0; jk <=nlstate ; jk++){
for(jk=1,pos=0; jk <=nlstate ; jk++)
pos += pp[jk];	pos += pp[jk];
	posprop += prop[jk][i];
	}
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
if(pos>=1.e-5){	if(pos>=1.e-5){
if(first==1)	if(first==1)
Line 1387 void freqsummary(char fileres[], int ag	Line 1588 void freqsummary(char fileres[], int ag
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);	printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);	fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
}	}
if( i <= (int) agemax){	if( i <= iagemax){
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,prop[jk][i]/posprop, prop[jk][i],posprop);
probs[i][jk][j1]= 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]);/	/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,prop[jk][i],posprop);
}	}
}	}

Line 1405 void freqsummary(char fileres[], int ag	Line 1606 void freqsummary(char fileres[], int ag
printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);	printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);	fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][i]);
}	}
if(i <= (int) agemax)	if(i <= iagemax)
fprintf(ficresp,"\n");	fprintf(ficresp,"\n");
if(first==1)	if(first==1)
printf("Others in log...\n");	printf("Others in log...\n");
Line 1416 void freqsummary(char fileres[], int ag	Line 1617 void freqsummary(char fileres[], int ag
dateintmean=dateintsum/k2cpt;	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, iagemin, iagemax+3);
free_vector(pp,1,nlstate);	free_vector(pp,1,nlstate);
	free_matrix(prop,1,nlstate,iagemin, iagemax+3);
/* End of Freq */	/* End of Freq */
}	}

/********** Prevalence ******************/	/********** Prevalence ******************/
void prevalence(int agemin, float 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)	void prevalence(double agemin, double agemax, int s, double agev, int nlstate, int imx, int Tvar, int nbcode, int ncodemax,double mint,double anint, double dateprev1,double dateprev2, int firstpass, int lastpass)
{ /* Some frequencies */	{
	/* Compute observed prevalence between dateprev1 and dateprev2 by counting the number of people
	in each health status at the date of interview (if between dateprev1 and dateprev2).
	We still use firstpass and lastpass as another selection.
	*/

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, *prop;
double pos, k2;	double pos,posprop;
	double y2; /* in fractional years */
pp=vector(1,nlstate);	int iagemin, iagemax;

freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);	iagemin= (int) agemin;
	iagemax= (int) agemax;
	/pp=vector(1,nlstate);/
	prop=matrix(1,nlstate,iagemin,iagemax+3);
	/* freq=ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);*/
j1=0;	j1=0;

j=cptcoveff;	j=cptcoveff;
Line 1443 void prevalence(int agemin, float agemax	Line 1652 void prevalence(int agemin, float agemax
for(i1=1; i1<=ncodemax[k1];i1++){	for(i1=1; i1<=ncodemax[k1];i1++){
j1++;	j1++;

for (i=-1; i<=nlstate+ndeath; i++)	for (i=1; i<=nlstate; i++)
for (jk=-1; jk<=nlstate+ndeath; jk++)	for(m=iagemin; m <= iagemax+3; m++)
for(m=agemin; m <= agemax+3; m++)	prop[i][m]=0.0;
freq[i][jk][m]=0;

for (i=1; i<=imx; i++) {	for (i=1; i<=imx; i++) { /* Each individual */
bool=1;	bool=1;
if (cptcovn>0) {	if (cptcovn>0) {
for (z1=1; z1<=cptcoveff; z1++)	for (z1=1; z1<=cptcoveff; z1++)
Line 1456 void prevalence(int agemin, float agemax	Line 1664 void prevalence(int agemin, float agemax
bool=0;	bool=0;
}	}
if (bool==1) {	if (bool==1) {
for(m=firstpass; m<=lastpass; m++){	for(m=firstpass; m<=lastpass; m++){/* Other selection (we can limit to certain interviews*/
k2=anint[m][i]+(mint[m][i]/12.);	y2=anint[m][i]+(mint[m][i]/12.); /* Fractional date in year */
if ((k2>=dateprev1) && (k2<=dateprev2)) {	if ((y2>=dateprev1) && (y2<=dateprev2)) { /* Here is the main selection (fractional years) */
if(agev[m][i]==0) agev[m][i]=agemax+1;	if(agev[m][i]==0) agev[m][i]=iagemax+1;
if(agev[m][i]==1) agev[m][i]=agemax+2;	if(agev[m][i]==1) agev[m][i]=iagemax+2;
if (m<lastpass) {	if((int)agev[m][i] <iagemin \|\| (int)agev[m][i] >iagemax+3) printf("Error on individual =%d agev[m][i]=%f m=%d\n",i, agev[m][i],m);
if (calagedate>0)	if (s[m][i]>0 && s[m][i]<=nlstate) {
freq[s[m][i]][s[m+1][i]][(int)(agev[m][i]+1-((int)calagedate %12)/12.)] += weight[i];	/if(i>4620) printf(" i=%d m=%d s[m][i]=%d (int)agev[m][i]=%d weight[i]=%f prop=%f\n",i,m,s[m][i],(int)agev[m][m],weight[i],prop[s[m][i]][(int)agev[m][i]]);/
else	prop[s[m][i]][(int)agev[m][i]] += weight[i];
freq[s[m][i]][s[m+1][i]][(int)agev[m][i]] += weight[i];	prop[s[m][i]][iagemax+3] += weight[i];
freq[s[m][i]][s[m+1][i]][(int)(agemax+3)] += weight[i];	}
}
}	}
}	} /* end selection of waves */
}	}
}	}
for(i=(int)agemin; i <= (int)agemax+3; i++){	for(i=iagemin; i <= iagemax+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(jk=1,posprop=0; jk <=nlstate ; jk++) {
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)	posprop += prop[jk][i];
pp[jk] += freq[jk][m][i];	}
}
	for(jk=1; jk <=nlstate ; jk++){
for(jk=1,pos=0; jk <=nlstate ; jk++) pos += pp[jk];	if( i <= iagemax){
	if(posprop>=1.e-5){
for(jk=1; jk <=nlstate ; jk++){	probs[i][jk][j1]= prop[jk][i]/posprop;
if( i <= (int) agemax){	}
if(pos>=1.e-5){	}
probs[i][jk][j1]= pp[jk]/pos;	}/* end jk */
}	}/* end i */
}
}/* end jk */
}/* end i */
} /* end i1 */	} /* end i1 */
} /* end k1 */	} /* end k1 */


free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath,(int) agemin,(int) agemax+3);	/* free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);*/
free_vector(pp,1,nlstate);	/free_vector(pp,1,nlstate);/
	free_matrix(prop,1,nlstate, iagemin,iagemax+3);
} /* End of Freq */	} /* End of prevalence */

/*********** Waves Concatenation *************/	/*********** Waves Concatenation *************/

Line 1571 void concatwav(int wav[], int **dh, int	Line 1766 void concatwav(int wav[], int **dh, int
if (j <= jmin) jmin=j;	if (j <= jmin) jmin=j;
sum=sum+j;	sum=sum+j;
/if (j<0) printf("j=%d num=%d \n",j,i); /	/if (j<0) printf("j=%d num=%d \n",j,i); /
	/* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
	/printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][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("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
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); */	/* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
	/printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);/
sum=sum+j;	sum=sum+j;
}	}
jk= j/stepm;	jk= j/stepm;
jl= j -jk*stepm;	jl= j -jk*stepm;
ju= j -(jk+1)*stepm;	ju= j -(jk+1)*stepm;
if(jl <= -ju){	if(mle <=1){
dh[mi][i]=jk;	if(jl==0){
bh[mi][i]=jl;	dh[mi][i]=jk;
}	bh[mi][i]=0;
else{	}else{ /* We want a negative bias in order to only have interpolation ie
dh[mi][i]=jk+1;	* at the price of an extra matrix product in likelihood */
bh[mi][i]=ju;	dh[mi][i]=jk+1;
}	bh[mi][i]=ju;
if(dh[mi][i]==0){	}
dh[mi][i]=1; /* At least one step */	}else{
bh[mi][i]=ju; /* At least one step */	if(jl <= -ju){
printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);	dh[mi][i]=jk;
	bh[mi][i]=jl; /* bias is positive if real duration
	* is higher than the multiple of stepm and negative otherwise.
	*/
	}
	else{
	dh[mi][i]=jk+1;
	bh[mi][i]=ju;
	}
	if(dh[mi][i]==0){
	dh[mi][i]=1; /* At least one step */
	bh[mi][i]=ju; /* At least one step */
	/* printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
	}
}	}
if(i==298 \|\| i==287 \|\| i==763 \|\|i==1061)printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d",bh[mi][i],ju,jl,dh[mi][i],jk,stepm);	} /* end if mle */
}	} /* end wave */
}
}	}
jmean=sum/k;	jmean=sum/k;
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);	printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
Line 1680 void evsij(char fileres[], double ***eij	Line 1891 void evsij(char fileres[], double ***eij
double *gradg, *trgradg;	double *gradg, *trgradg;
int theta;	int theta;

varhe=ma3x(1,nlstate2,1,nlstate2,(int) bage, (int) fage);	varhe=ma3x(1,nlstatenlstate,1,nlstatenlstate,(int) bage, (int) fage);
xp=vector(1,npar);	xp=vector(1,npar);
dnewm=matrix(1,nlstate*2,1,npar);	dnewm=matrix(1,nlstate*nlstate,1,npar);
doldm=matrix(1,nlstate2,1,nlstate2);	doldm=matrix(1,nlstatenlstate,1,nlstatenlstate);

fprintf(ficreseij,"# Health expectancies\n");	fprintf(ficreseij,"# Health expectancies\n");
fprintf(ficreseij,"# Age");	fprintf(ficreseij,"# Age");
Line 1700 void evsij(char fileres[], double ***eij	Line 1911 void evsij(char fileres[], double ***eij
* This is mainly to measure the difference between two models: for example	* This is mainly to measure the difference between two models: for example
* if stepm=24 months pijx are given only every 2 years and by summing them	* if stepm=24 months pijx are given only every 2 years and by summing them
* we are calculating an estimate of the Life Expectancy assuming a linear	* we are calculating an estimate of the Life Expectancy assuming a linear
* progression inbetween and thus overestimating or underestimating according	* progression in between and thus overestimating or underestimating according
* to the curvature of the survival function. If, for the same date, we	* to the curvature of the survival function. If, for the same date, we
* estimate the model with stepm=1 month, we can keep estepm to 24 months	* estimate the model with stepm=1 month, we can keep estepm to 24 months
* to compare the new estimate of Life expectancy with the same linear	* to compare the new estimate of Life expectancy with the same linear
Line 1729 void evsij(char fileres[], double ***eij	Line 1940 void evsij(char fileres[], double ***eij
/* if (stepm >= YEARM) hstepm=1;*/	/* if (stepm >= YEARM) hstepm=1;*/
nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */	nhstepm = nstepm/hstepm;/* Expressed in hstepm, typically nhstepm=40/4=10 */
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
gradg=ma3x(0,nhstepm,1,npar,1,nlstate*2);	gradg=ma3x(0,nhstepm,1,npar,1,nlstate*nlstate);
gp=matrix(0,nhstepm,1,nlstate*2);	gp=matrix(0,nhstepm,1,nlstate*nlstate);
gm=matrix(0,nhstepm,1,nlstate*2);	gm=matrix(0,nhstepm,1,nlstate*nlstate);

/* Computed by stepm unit matrices, product of hstepm matrices, stored	/* Computed by stepm unit matrices, product of hstepm matrices, stored
in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */	in an array of nhstepm length: nhstepm=10, hstepm=4, stepm=6 months */
Line 1772 void evsij(char fileres[], double ***eij	Line 1983 void evsij(char fileres[], double ***eij
}	}
}	}
}	}
for(j=1; j<= nlstate*2; j++)	for(j=1; j<= nlstate*nlstate; j++)
for(h=0; h<=nhstepm-1; h++){	for(h=0; h<=nhstepm-1; 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 1780 void evsij(char fileres[], double ***eij	Line 1991 void evsij(char fileres[], double ***eij

/* End theta */	/* End theta */

trgradg =ma3x(0,nhstepm,1,nlstate*2,1,npar);	trgradg =ma3x(0,nhstepm,1,nlstate*nlstate,1,npar);

for(h=0; h<=nhstepm-1; h++)	for(h=0; h<=nhstepm-1; h++)
for(j=1; j<=nlstate*2;j++)	for(j=1; j<=nlstate*nlstate;j++)
for(theta=1; theta <=npar; theta++)	for(theta=1; theta <=npar; theta++)
trgradg[h][j][theta]=gradg[h][theta][j];	trgradg[h][j][theta]=gradg[h][theta][j];


for(i=1;i<=nlstate*2;i++)	for(i=1;i<=nlstate*nlstate;i++)
for(j=1;j<=nlstate*2;j++)	for(j=1;j<=nlstate*nlstate;j++)
varhe[i][j][(int)age] =0.;	varhe[i][j][(int)age] =0.;

printf("%d\|",(int)age);fflush(stdout);	printf("%d\|",(int)age);fflush(stdout);
fprintf(ficlog,"%d\|",(int)age);fflush(ficlog);	fprintf(ficlog,"%d\|",(int)age);fflush(ficlog);
for(h=0;h<=nhstepm-1;h++){	for(h=0;h<=nhstepm-1;h++){
for(k=0;k<=nhstepm-1;k++){	for(k=0;k<=nhstepm-1;k++){
matprod2(dnewm,trgradg[h],1,nlstate*2,1,npar,1,npar,matcov);	matprod2(dnewm,trgradg[h],1,nlstate*nlstate,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate2,1,npar,1,nlstate2,gradg[k]);	matprod2(doldm,dnewm,1,nlstatenlstate,1,npar,1,nlstatenlstate,gradg[k]);
for(i=1;i<=nlstate*2;i++)	for(i=1;i<=nlstate*nlstate;i++)
for(j=1;j<=nlstate*2;j++)	for(j=1;j<=nlstate*nlstate;j++)
varhe[i][j][(int)age] += doldm[i][j]hfhf;	varhe[i][j][(int)age] += doldm[i][j]hfhf;
}	}
}	}
Line 1822 void evsij(char fileres[], double ***eij	Line 2033 void evsij(char fileres[], double ***eij
}	}
fprintf(ficreseij,"\n");	fprintf(ficreseij,"\n");

free_matrix(gm,0,nhstepm,1,nlstate*2);	free_matrix(gm,0,nhstepm,1,nlstate*nlstate);
free_matrix(gp,0,nhstepm,1,nlstate*2);	free_matrix(gp,0,nhstepm,1,nlstate*nlstate);
free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*2);	free_ma3x(gradg,0,nhstepm,1,npar,1,nlstate*nlstate);
free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);	free_ma3x(trgradg,0,nhstepm,1,nlstate*nlstate,1,npar);
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
}	}
printf("\n");	printf("\n");
fprintf(ficlog,"\n");	fprintf(ficlog,"\n");

free_vector(xp,1,npar);	free_vector(xp,1,npar);
free_matrix(dnewm,1,nlstate*2,1,npar);	free_matrix(dnewm,1,nlstate*nlstate,1,npar);
free_matrix(doldm,1,nlstate2,1,nlstate2);	free_matrix(doldm,1,nlstatenlstate,1,nlstatenlstate);
free_ma3x(varhe,1,nlstate2,1,nlstate2,(int) bage, (int)fage);	free_ma3x(varhe,1,nlstatenlstate,1,nlstatenlstate,(int) bage, (int)fage);
}	}

/********** Variance ****************/	/********** Variance ****************/
Line 1890 void varevsij(char optionfilefiname[], d	Line 2101 void varevsij(char optionfilefiname[], d
}	}
printf("Computing total mortality p.j=w1p1j+w2p2j+..: result on file '%s' \n",fileresprobmorprev);	printf("Computing total mortality p.j=w1p1j+w2p2j+..: result on file '%s' \n",fileresprobmorprev);
fprintf(ficlog,"Computing total mortality p.j=w1p1j+w2p2j+..: result on file '%s' \n",fileresprobmorprev);	fprintf(ficlog,"Computing total mortality p.j=w1p1j+w2p2j+..: result on file '%s' \n",fileresprobmorprev);
fprintf(ficresprobmorprev,"# probabilities of dying during a year and weighted mean w1p1j+w2p2j+... stand dev in()\n");	fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1p1j+w2p2j+... stand dev in()\n",estepm);
fprintf(ficresprobmorprev,"# Age cov=%-d",ij);	fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
for(j=nlstate+1; j<=(nlstate+ndeath);j++){	for(j=nlstate+1; j<=(nlstate+ndeath);j++){
fprintf(ficresprobmorprev," p.%-d SE",j);	fprintf(ficresprobmorprev," p.%-d SE",j);
Line 1912 void varevsij(char optionfilefiname[], d	Line 2123 void varevsij(char optionfilefiname[], d
exit(0);	exit(0);
}	}
else{	else{
fprintf(fichtm,"\n<li><h4> Computing probabilities of dying as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");	fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
fprintf(fichtm,"\n<br>%s (� revoir) <br>\n",digitp);	fprintf(fichtm,"\n<br>%s <br>\n",digitp);
}	}
varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);	varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);

Line 1947 void varevsij(char optionfilefiname[], d	Line 2158 void varevsij(char optionfilefiname[], d
and note for a fixed period like k years */	and note for a fixed period like k years */
/* We decided (b) to get a life expectancy respecting the most precise curvature of the	/* We decided (b) to get a life expectancy respecting the most precise curvature of the
survival function given by stepm (the optimization length). Unfortunately it	survival function given by stepm (the optimization length). Unfortunately it
means that if the survival funtion is printed only each two years of age and if	means that if the survival funtion is printed every two years of age and if
you sum them up and add 1 year (area under the trapezoids) you won't get the same	you sum them up and add 1 year (area under the trapezoids) you won't get the same
results. So we changed our mind and took the option of the best precision.	results. So we changed our mind and took the option of the best precision.
*/	*/
Line 1963 void varevsij(char optionfilefiname[], d	Line 2174 void varevsij(char optionfilefiname[], d


for(theta=1; theta <=npar; theta++){	for(theta=1; theta <=npar; theta++){
for(i=1; i<=npar; i++){ /* Computes gradient */	for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/
xp[i] = x[i] + (i==theta ?delti[theta]:0);	xp[i] = x[i] + (i==theta ?delti[theta]:0);
}	}
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);	hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
Line 1985 void varevsij(char optionfilefiname[], d	Line 2196 void varevsij(char optionfilefiname[], d
gp[h][j] += prlim[i][i]*p3mat[i][j][h];	gp[h][j] += prlim[i][i]*p3mat[i][j][h];
}	}
}	}
/* This for computing forces of mortality (h=1)as a weighted average */	/* This for computing probability of death (h=1 means
for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){	computed over hstepm matrices product = hstepm*stepm months)
for(i=1; i<= nlstate; i++)	as a weighted average of prlim.
	*/
	for(j=nlstate+1;j<=nlstate+ndeath;j++){
	for(i=1,gpp[j]=0.; i<= nlstate; i++)
gpp[j] += prlim[i][i]*p3mat[i][j][1];	gpp[j] += prlim[i][i]*p3mat[i][j][1];
}	}
/* end force of mortality */	/* end probability of death */

for(i=1; i<=npar; i++) /* Computes gradient */	for(i=1; i<=npar; i++) /* Computes gradient x - delta */
xp[i] = x[i] - (i==theta ?delti[theta]:0);	xp[i] = x[i] - (i==theta ?delti[theta]:0);
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);	hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij);
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ij);
Line 2013 void varevsij(char optionfilefiname[], d	Line 2227 void varevsij(char optionfilefiname[], d
gm[h][j] += prlim[i][i]*p3mat[i][j][h];	gm[h][j] += prlim[i][i]*p3mat[i][j][h];
}	}
}	}
/* This for computing force of mortality (h=1)as a weighted average */	/* This for computing probability of death (h=1 means
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){	computed over hstepm matrices product = hstepm*stepm months)
for(i=1; i<= nlstate; i++)	as a weighted average of prlim.
gmp[j] += prlim[i][i]*p3mat[i][j][1];	*/
	for(j=nlstate+1;j<=nlstate+ndeath;j++){
	for(i=1,gmp[j]=0.; i<= nlstate; i++)
	gmp[j] += prlim[i][i]*p3mat[i][j][1];
}	}
/* end force of mortality */	/* end probability of death */

for(j=1; j<= nlstate; j++) /* vareij */	for(j=1; j<= nlstate; j++) /* vareij */
for(h=0; h<=nhstepm; h++){	for(h=0; h<=nhstepm; h++){
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];
}	}

for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */	for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */
gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];	gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
}	}
Line 2040 void varevsij(char optionfilefiname[], d	Line 2258 void varevsij(char optionfilefiname[], d
for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */	for(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
for(theta=1; theta <=npar; theta++)	for(theta=1; theta <=npar; theta++)
trgradgp[j][theta]=gradgp[theta][j];	trgradgp[j][theta]=gradgp[theta][j];


hf=hstepmstepm/YEARM; / Duration of hstepm expressed in year unit. */	hf=hstepmstepm/YEARM; / Duration of hstepm expressed in year unit. */
for(i=1;i<=nlstate;i++)	for(i=1;i<=nlstate;i++)
Line 2055 void varevsij(char optionfilefiname[], d	Line 2274 void varevsij(char optionfilefiname[], d
vareij[i][j][(int)age] += doldm[i][j]hfhf;	vareij[i][j][(int)age] += doldm[i][j]hfhf;
}	}
}	}

/* pptj */	/* pptj */
matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);	matprod2(dnewmp,trgradgp,nlstate+1,nlstate+ndeath,1,npar,1,npar,matcov);
matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);	matprod2(doldmp,dnewmp,nlstate+1,nlstate+ndeath,1,npar,nlstate+1,nlstate+ndeath,gradgp);
Line 2063 void varevsij(char optionfilefiname[], d	Line 2282 void varevsij(char optionfilefiname[], d
for(i=nlstate+1;i<=nlstate+ndeath;i++)	for(i=nlstate+1;i<=nlstate+ndeath;i++)
varppt[j][i]=doldmp[j][i];	varppt[j][i]=doldmp[j][i];
/* end ppptj */	/* end ppptj */
	/* x centered again */
hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);	hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);	prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);

Line 2075 void varevsij(char optionfilefiname[], d	Line 2295 void varevsij(char optionfilefiname[], d
prlim[i][i]=mobaverage[(int)age][i][ij];	prlim[i][i]=mobaverage[(int)age][i][ij];
}	}
}	}

/* This for computing force of mortality (h=1)as a weighted average */	/* This for computing probability of death (h=1 means
for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){	computed over hstepm (estepm) matrices product = hstepm*stepm months)
for(i=1; i<= nlstate; i++)	as a weighted average of prlim.
	*/
	for(j=nlstate+1;j<=nlstate+ndeath;j++){
	for(i=1,gmp[j]=0.;i<= nlstate; i++)
gmp[j] += prlim[i][i]*p3mat[i][j][1];	gmp[j] += prlim[i][i]*p3mat[i][j][1];
}	}
/* end force of mortality */	/* end probability of death */

fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);	fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
for(j=nlstate+1; j<=(nlstate+ndeath);j++){	for(j=nlstate+1; j<=(nlstate+ndeath);j++){
Line 2111 void varevsij(char optionfilefiname[], d	Line 2334 void varevsij(char optionfilefiname[], d
fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65");	fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65");
/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ // Only the first actually */	/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ // Only the first actually */
fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");	fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm);	/* fprintf(ficgp,"\n plot \"%s\" u 1:($3%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); /
fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96$4)%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm);	/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96$4)%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96$4)%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm);	/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96$4)%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
	fprintf(ficgp,"\n plot \"%s\" u 1:($3) not w l 1 ",fileresprobmorprev);
	fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)) t \"95\%% interval\" w l 2 ",fileresprobmorprev);
	fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)) not w l 2 ",fileresprobmorprev);
fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev);	fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",fileresprobmorprev,fileresprobmorprev);
fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,digitp,digit);	fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"varmuptjgr%s%s%s.png\"> <br>\n", estepm,digitp,optionfilefiname,digit);
/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);	/* fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months and then divided by estepm and multiplied by %.0f in order to have the probability to die over a year <br> <img src=\"varmuptjgr%s%s.png\"> <br>\n", stepm,YEARM,digitp,digit);
*/	*/
fprintf(ficgp,"\nset out \"varmuptjgr%s%s.png\";replot;",digitp,digit);	fprintf(ficgp,"\nset out \"varmuptjgr%s%s%s.png\";replot;",digitp,optionfilefiname,digit);

free_vector(xp,1,npar);	free_vector(xp,1,npar);
free_matrix(doldm,1,nlstate,1,nlstate);	free_matrix(doldm,1,nlstate,1,nlstate);
Line 2130 void varevsij(char optionfilefiname[], d	Line 2356 void varevsij(char optionfilefiname[], d
fclose(ficresprobmorprev);	fclose(ficresprobmorprev);
fclose(ficgp);	fclose(ficgp);
fclose(fichtm);	fclose(fichtm);
}	}

/********** Variance of prevlim ****************/	/********** Variance of prevlim ****************/
void varprevlim(char fileres[], double varpl, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double **prlim, double ftolpl, int ij)	void varprevlim(char fileres[], double varpl, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double **prlim, double ftolpl, int ij)
Line 2276 void varprob(char optionfilefiname[], do	Line 2502 void varprob(char optionfilefiname[], do
fprintf(ficresprobcov," p%1d-%1d ",i,j);	fprintf(ficresprobcov," p%1d-%1d ",i,j);
fprintf(ficresprobcor," p%1d-%1d ",i,j);	fprintf(ficresprobcor," p%1d-%1d ",i,j);
}	}
fprintf(ficresprob,"\n");	/* fprintf(ficresprob,"\n");
fprintf(ficresprobcov,"\n");	fprintf(ficresprobcov,"\n");
fprintf(ficresprobcor,"\n");	fprintf(ficresprobcor,"\n");
xp=vector(1,npar);	*/
	xp=vector(1,npar);
dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);	dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
doldm=matrix(1,(nlstate)(nlstate+ndeath),1,(nlstate)(nlstate+ndeath));	doldm=matrix(1,(nlstate)(nlstate+ndeath),1,(nlstate)(nlstate+ndeath));
mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);	mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);
Line 2318 void varprob(char optionfilefiname[], do	Line 2545 void varprob(char optionfilefiname[], do
if (cptcovn>0) {	if (cptcovn>0) {
fprintf(ficresprob, "\n#********** Variable ");	fprintf(ficresprob, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
fprintf(ficresprob, "**********\n#");	fprintf(ficresprob, "**********\n#\n");
fprintf(ficresprobcov, "\n#********** Variable ");	fprintf(ficresprobcov, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
fprintf(ficresprobcov, "**********\n#");	fprintf(ficresprobcov, "**********\n#\n");

fprintf(ficgp, "\n#********** Variable ");	fprintf(ficgp, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, "# V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficgp, " V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
fprintf(ficgp, "**********\n#");	fprintf(ficgp, "**********\n#\n");


fprintf(fichtm, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable ");	fprintf(fichtm, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable ");
Line 2334 void varprob(char optionfilefiname[], do	Line 2561 void varprob(char optionfilefiname[], do

fprintf(ficresprobcor, "\n#********** Variable ");	fprintf(ficresprobcor, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcor, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
fprintf(ficgp, "**********\n#");	fprintf(ficresprobcor, "**********\n#");
}	}

for (age=bage; age<=fage; age ++){	for (age=bage; age<=fage; age ++){
Line 2353 void varprob(char optionfilefiname[], do	Line 2580 void varprob(char optionfilefiname[], do

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

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

Line 2366 void varprob(char optionfilefiname[], do	Line 2593 void varprob(char optionfilefiname[], do
}	}

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

pmij(pmmij,cov,ncovmodel,xp,nlstate);	pmij(pmmij,cov,ncovmodel,xp,nlstate);
k=0;	k=0;
Line 2378 void varprob(char optionfilefiname[], do	Line 2605 void varprob(char optionfilefiname[], do
}	}

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

for(j=1; j<=(nlstate)*(nlstate+ndeath);j++)	for(j=1; j<=(nlstate)*(nlstate+ndeath);j++)
Line 2653 m=pow(2,cptcoveff);	Line 2880 m=pow(2,cptcoveff);
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");	if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%lf (\%%lf)");	else fprintf(ficgp," \%%lf (\%%lf)");
}	}
fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+2*$3) \"\%%lf",fileres,k1-1,k1-1);	fprintf(ficgp,"\" t\"Stable prevalence\" w l 0,\"vpl%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",fileres,k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {	for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");	if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%lf (\%%lf)");	else fprintf(ficgp," \%%lf (\%%lf)");
}	}
fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-2*$3) \"\%%lf",fileres,k1-1,k1-1);	fprintf(ficgp,"\" t\"95\%% CI\" w l 1,\"vpl%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",fileres,k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {	for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");	if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
else fprintf(ficgp," \%%lf (\%%lf)");	else fprintf(ficgp," \%%lf (\%%lf)");
Line 2849 int movingaverage(double ***probs, doubl	Line 3076 int movingaverage(double ***probs, doubl


/************ Forecasting ****************/	/************ Forecasting ****************/
prevforecast(char fileres[], double anproj1,double mproj1,double jproj1,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anproj2,double p[], int i2){	prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
	/* proj1, year, month, day of starting projection
int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;	agemin, agemax range of age
	dateprev1 dateprev2 range of dates during which prevalence is computed
	anproj2 year of en of projection (same day and month as proj1).
	*/
	int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h, i1;
int *popage;	int *popage;
double calagedate, agelim, kk1, kk2, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;	double agec; /* generic age */
	double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
double popeffectif,popcount;	double popeffectif,popcount;
double ***p3mat;	double ***p3mat;
double ***mobaverage;	double ***mobaverage;
char fileresf[FILENAMELENGTH];	char fileresf[FILENAMELENGTH];

agelim=AGESUP;	agelim=AGESUP;
calagedate=(anproj1+mproj1/12.+jproj1/365.-dateintmean)*YEARM;	prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);

prevalence(ageminpar, 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 2886 prevforecast(char fileres[], double anpr	Line 3115 prevforecast(char fileres[], double anpr

stepsize=(int) (stepm+YEARM-1)/YEARM;	stepsize=(int) (stepm+YEARM-1)/YEARM;
if (stepm<=12) stepsize=1;	if (stepm<=12) stepsize=1;
	if(estepm < stepm){
agelim=AGESUP;	printf ("Problem %d lower than %d\n",estepm, stepm);
	}
hstepm=1;	else hstepm=estepm;

hstepm=hstepm/stepm;	hstepm=hstepm/stepm;
yp1=modf(dateintmean,&yp);	yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
	fractional in yp1 */
anprojmean=yp;	anprojmean=yp;
yp2=modf((yp1*12),&yp);	yp2=modf((yp1*12),&yp);
mprojmean=yp;	mprojmean=yp;
Line 2899 prevforecast(char fileres[], double anpr	Line 3130 prevforecast(char fileres[], double anpr
jprojmean=yp;	jprojmean=yp;
if(jprojmean==0) jprojmean=1;	if(jprojmean==0) jprojmean=1;
if(mprojmean==0) jprojmean=1;	if(mprojmean==0) jprojmean=1;

	i1=cptcoveff;
	if (cptcovn < 1){i1=1;}

fprintf(ficresf,"# Estimated date of observed prevalence: %.lf/%.lf/%.lf ",jprojmean,mprojmean,anprojmean);	fprintf(ficresf,"# Mean day of interviews %.lf/%.lf/%.lf (%.2f) between %.2f and %.2f \n",jprojmean,mprojmean,anprojmean,dateintmean,dateprev1,dateprev2);

for(cptcov=1;cptcov<=i2;cptcov++){	fprintf(ficresf,"#**** Routine prevforecast \n");

	for(cptcov=1, k=0;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, hpijx=probability over h years, hp.jx is weighted by observed prev ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
}	}
fprintf(ficresf,"******\n");	fprintf(ficresf,"******\n");
fprintf(ficresf,"# StartingAge FinalAge");	fprintf(ficresf,"# Covariate valuofcovar yearproj age");
for(j=1; j<=nlstate+ndeath;j++) fprintf(ficresf," P.%d",j);	for(j=1; j<=nlstate+ndeath;j++){
	for(i=1; i<=nlstate;i++)
	fprintf(ficresf," p%d%d",i,j);
for (cpt=0; cpt<=(anproj2-anproj1);cpt++) {	fprintf(ficresf," p.%d",j);
	}
	for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
fprintf(ficresf,"\n");	fprintf(ficresf,"\n");
fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+cpt);	fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);

for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){	for (agec=fage; agec>=(ageminpar-1); agec--){
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);	nhstepm=(int) rint((agelim-agec)*YEARM/stepm);
nhstepm = nhstepm/hstepm;	nhstepm = nhstepm/hstepm;

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,agec,hstepm,p,nlstate,stepm,oldm,savm, k);

for (h=0; h<=nhstepm; h++){	for (h=0; h<=nhstepm; h++){
if (h==(int) (calagedate+YEARM*cpt)) {	if (h==(int) (YEARM*yearp)) {
fprintf(ficresf,"\n %.f %.f ",anproj1+cpt,agedeb+hhstepm/YEARMstepm);	fprintf(ficresf,"\n");
	for(j=1;j<=cptcoveff;j++)
	fprintf(ficresf,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
	fprintf(ficresf,"%.f %.f ",anproj1+yearp,agec+hhstepm/YEARMstepm);
}	}
for(j=1; j<=nlstate+ndeath;j++) {	for(j=1; j<=nlstate+ndeath;j++) {
kk1=0.;kk2=0;	ppij=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+1][i][cptcod];	ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][cptcod];
else {	else {
kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];	ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][cptcod];
}	}
	if (h==(int)(YEARM*yearp))
	fprintf(ficresf," %.3f", p3mat[i][j][h]);
}	}
if (h==(int)(calagedate+12*cpt)){	if (h==(int)(YEARM*yearp)){
fprintf(ficresf," %.3f", kk1);	fprintf(ficresf," %.3f", ppij);

}	}
}	}
}	}
Line 2956 prevforecast(char fileres[], double anpr	Line 3196 prevforecast(char fileres[], double anpr

fclose(ficresf);	fclose(ficresf);
}	}
/************ Forecasting ****************/
	/************ Forecasting *not tested NB***********/
populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){	populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){

int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;	int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
int *popage;	int *popage;
double calagedate, agelim, kk1, kk2;	double calagedatem, agelim, kk1, kk2;
double popeffectif,popcount;	double popeffectif,popcount;
double *p3mat,tabpop,**tabpopprev;	double *p3mat,tabpop,**tabpopprev;
double ***mobaverage;	double ***mobaverage;
Line 2970 populforecast(char fileres[], double anp	Line 3211 populforecast(char fileres[], double anp
tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);	tabpop= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);	tabpopprev= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
agelim=AGESUP;	agelim=AGESUP;
calagedate=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;	calagedatem=(anpyram+mpyram/12.+jpyram/365.-dateintmean)*YEARM;

prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);	prevalence(ageminpar, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);


strcpy(filerespop,"pop");	strcpy(filerespop,"pop");
Line 3018 populforecast(char fileres[], double anp	Line 3259 populforecast(char fileres[], double anp
for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];	for (i=1; i<imx;i++) popeffectif[popage[i]]=popcount[i];
}	}

for(cptcov=1;cptcov<=i2;cptcov++){	for(cptcov=1,k=0;cptcov<=i2;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){	for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){
k=k+1;	k=k+1;
fprintf(ficrespop,"\n#******");	fprintf(ficrespop,"\n#******");
Line 3033 populforecast(char fileres[], double anp	Line 3274 populforecast(char fileres[], double anp
for (cpt=0; cpt<=0;cpt++) {	for (cpt=0; cpt<=0;cpt++) {
fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);	fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);

for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){	for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);	nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
nhstepm = nhstepm/hstepm;	nhstepm = nhstepm/hstepm;

Line 3042 populforecast(char fileres[], double anp	Line 3283 populforecast(char fileres[], double anp
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++){	for (h=0; h<=nhstepm; h++){
if (h==(int) (calagedate+YEARM*cpt)) {	if (h==(int) (calagedatem+YEARM*cpt)) {
fprintf(ficrespop,"\n %3.f ",agedeb+hhstepm/YEARMstepm);	fprintf(ficrespop,"\n %3.f ",agedeb+hhstepm/YEARMstepm);
}	}
for(j=1; j<=nlstate+ndeath;j++) {	for(j=1; j<=nlstate+ndeath;j++) {
Line 3054 populforecast(char fileres[], double anp	Line 3295 populforecast(char fileres[], double anp
kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];	kk1=kk1+p3mat[i][j][h]*probs[(int)(agedeb+1)][i][cptcod];
}	}
}	}
if (h==(int)(calagedate+12*cpt)){	if (h==(int)(calagedatem+12*cpt)){
tabpop[(int)(agedeb)][j][cptcod]=kk1;	tabpop[(int)(agedeb)][j][cptcod]=kk1;
/*fprintf(ficrespop," %.3f", kk1);	/*fprintf(ficrespop," %.3f", kk1);
if (popforecast==1) fprintf(ficrespop," [%.f]", kk1popeffectif[(int)agedeb+1]);/	if (popforecast==1) fprintf(ficrespop," [%.f]", kk1popeffectif[(int)agedeb+1]);/
Line 3065 populforecast(char fileres[], double anp	Line 3306 populforecast(char fileres[], double anp
for(j=1; j<=nlstate;j++){	for(j=1; j<=nlstate;j++){
kk1= kk1+tabpop[(int)(agedeb)][j][cptcod];	kk1= kk1+tabpop[(int)(agedeb)][j][cptcod];
}	}
tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1popeffectif[(int)(agedeb+(calagedate+12cpt)hstepm/YEARMstepm-1)];	tabpopprev[(int)(agedeb)][i][cptcod]=tabpop[(int)(agedeb)][i][cptcod]/kk1popeffectif[(int)(agedeb+(calagedatem+12cpt)hstepm/YEARMstepm-1)];
}	}

if (h==(int)(calagedate+12*cpt)) for(j=1; j<=nlstate;j++)	if (h==(int)(calagedatem+12*cpt)) for(j=1; j<=nlstate;j++)
fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);	fprintf(ficrespop," %15.2f",tabpopprev[(int)(agedeb+1)][j][cptcod]);
}	}
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
Line 3079 populforecast(char fileres[], double anp	Line 3320 populforecast(char fileres[], double anp

for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {	for (cpt=1; cpt<=(anpyram1-anpyram);cpt++) {
fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);	fprintf(ficrespop,"\n\n# Forecasting at date %.lf/%.lf/%.lf ",jpyram,mpyram,anpyram+cpt);
for (agedeb=(fage-((int)calagedate %12/12.)); agedeb>=(ageminpar-((int)calagedate %12)/12.); agedeb--){	for (agedeb=(fage-((int)calagedatem %12/12.)); agedeb>=(ageminpar-((int)calagedatem %12)/12.); agedeb--){
nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);	nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm);
nhstepm = nhstepm/hstepm;	nhstepm = nhstepm/hstepm;

Line 3087 populforecast(char fileres[], double anp	Line 3328 populforecast(char fileres[], double anp
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++){	for (h=0; h<=nhstepm; h++){
if (h==(int) (calagedate+YEARM*cpt)) {	if (h==(int) (calagedatem+YEARM*cpt)) {
fprintf(ficresf,"\n %3.f ",agedeb+hhstepm/YEARMstepm);	fprintf(ficresf,"\n %3.f ",agedeb+hhstepm/YEARMstepm);
}	}
for(j=1; j<=nlstate+ndeath;j++) {	for(j=1; j<=nlstate+ndeath;j++) {
Line 3095 populforecast(char fileres[], double anp	Line 3336 populforecast(char fileres[], double anp
for(i=1; i<=nlstate;i++) {	for(i=1; i<=nlstate;i++) {
kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod];	kk1=kk1+p3mat[i][j][h]*tabpopprev[(int)agedeb+1][i][cptcod];
}	}
if (h==(int)(calagedate+12*cpt)) fprintf(ficresf," %15.2f", kk1);	if (h==(int)(calagedatem+12*cpt)) fprintf(ficresf," %15.2f", kk1);
}	}
}	}
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
Line 3122 populforecast(char fileres[], double anp	Line 3363 populforecast(char fileres[], double anp

int main(int argc, char *argv[])	int main(int argc, char *argv[])
{	{
	int movingaverage(double *probs, double bage,double fage, double *mobaverage, int mobilav);
int i,j, k, n=MAXN,iter,m,size,cptcode, cptcod;	int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
double agedeb, agefin,hf;	double agedeb, agefin,hf;
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;	double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;

Line 3142 int main(int argc, char *argv[])	Line 3383 int main(int argc, char *argv[])
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,tab;	int jnais,jdc,jint4,jint1,jint2,jint3,*outcome,tab;
	int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
int mobilav=0,popforecast=0;	int mobilav=0,popforecast=0;
int hstepm, nhstepm;	int hstepm, nhstepm;
double jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,jpyram, mpyram,anpyram,jpyram1, mpyram1,anpyram1, calagedate;	double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;
	double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;

double bage, fage, age, agelim, agebase;	double bage, fage, age, agelim, agebase;
double ftolpl=FTOL;	double ftolpl=FTOL;
Line 3159 int main(int argc, char *argv[])	Line 3402 int main(int argc, char *argv[])
double *varpl; / Variances of prevalence limits by age */	double *varpl; / Variances of prevalence limits by age */
double *epj, vepp;	double *epj, vepp;
double kk1, kk2;	double kk1, kk2;
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;	double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
/int movingaverage; */

char *alph[]={"a","a","b","c","d","e"}, str[4];	char *alph[]={"a","a","b","c","d","e"}, str[4];

Line 3307 int main(int argc, char *argv[])	Line 3549 int main(int argc, char *argv[])
ungetc(c,ficpar);	ungetc(c,ficpar);

delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);	delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
delti=vector(1,npar); /* Scale of each paramater (output from hesscov) */	/* delti=vector(1,npar); // Scale of each paramater (output from hesscov) */
for(i=1; i <=nlstate; i++){	for(i=1; i <=nlstate; i++){
for(j=1; j <=nlstate+ndeath-1; j++){	for(j=1; j <=nlstate+ndeath-1; j++){
fscanf(ficpar,"%1d%1d",&i1,&j1);	fscanf(ficpar,"%1d%1d",&i1,&j1);
Line 3324 int main(int argc, char *argv[])	Line 3566 int main(int argc, char *argv[])
}	}
}	}
delti=delti3[1][1];	delti=delti3[1][1];


	/* free_ma3x(delti3,1,nlstate,1,nlstate+ndeath-1,1,ncovmodel); / / Hasn't to to freed here otherwise delti is no more allocated */

/* Reads comments: lines beginning with '#' */	/* Reads comments: lines beginning with '#' */
while((c=getc(ficpar))=='#' && c!= EOF){	while((c=getc(ficpar))=='#' && c!= EOF){
Line 3443 int main(int argc, char *argv[])	Line 3688 int main(int argc, char *argv[])
if (s[4][i]==9) s[4][i]=-1;	if (s[4][i]==9) s[4][i]=-1;
printf("%d %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}*/	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 ((s[3][i]==3) \|\| (s[4][i]==3)) weight[i]=0.08;
	else weight[i]=1;*/

/* Calculation of the number of parameter from char model*/	/* Calculation of the number of parameter from char model*/
Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */	Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
Tprod=ivector(1,15);	Tprod=ivector(1,15);
Line 3545 int main(int argc, char *argv[])	Line 3794 int main(int argc, char *argv[])

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=firstpass; (m<= lastpass); m++){
if(s[m][i] >0){	if(s[m][i] >0){
if (s[m][i] >= nlstate+1) {	if (s[m][i] >= nlstate+1) {
if(agedc[i]>0)	if(agedc[i]>0)
Line 3560 int main(int argc, char *argv[])	Line 3809 int main(int argc, char *argv[])
}	}
}	}
}	}
else if(s[m][i] !=9){ /* Should no more exist */	else if(s[m][i] !=9){ /* Standard case, age in fractional
	years but with the precision of a
	month */
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;
Line 3586 int main(int argc, char *argv[])	Line 3837 int main(int argc, char *argv[])

}	}
for (i=1; i<=imx; i++) {	for (i=1; i<=imx; i++) {
for(m=1; (m<= maxwav); m++){	for(m=firstpass; (m<=lastpass); m++){
if (s[m][i] > (nlstate+ndeath)) {	if (s[m][i] > (nlstate+ndeath)) {
printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);	printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);	fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);
Line 3595 int main(int argc, char *argv[])	Line 3846 int main(int argc, char *argv[])
}	}
}	}

	/*for (i=1; i<=imx; i++){
	for (m=firstpass; (m<lastpass); m++){
	printf("%d %d %.lf %d %d\n", num[i],(covar[1][i]),agev[m][i],s[m][i],s[m+1][i]);
	}

	}*/

printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);	printf("Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);
fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);	fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);

Line 3652 int main(int argc, char *argv[])	Line 3910 int main(int argc, char *argv[])
/* Calculates basic frequencies. Computes observed prevalence at single age	/* Calculates basic frequencies. Computes observed prevalence at single age
and prints on file fileres'p'. */	and prints on file fileres'p'. */

	pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
	oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
	newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
	savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
	oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */


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

if(mle==1){	if(mle>=1){ /* Could be 1 or 2 */
mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);	mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
}	}

Line 3774 int main(int argc, char *argv[])	Line 4037 int main(int argc, char *argv[])
fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);	fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf mov_average=%d\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2,&mobilav);
fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);	fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);	fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
	printf("begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
	fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);

while((c=getc(ficpar))=='#' && c!= EOF){	while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);	ungetc(c,ficpar);
Line 3784 int main(int argc, char *argv[])	Line 4049 int main(int argc, char *argv[])
ungetc(c,ficpar);	ungetc(c,ficpar);


dateprev1=anprev1+mprev1/12.+jprev1/365.;	dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;
dateprev2=anprev2+mprev2/12.+jprev2/365.;	dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;

fscanf(ficpar,"pop_based=%d\n",&popbased);	fscanf(ficpar,"pop_based=%d\n",&popbased);
fprintf(ficparo,"pop_based=%d\n",popbased);	fprintf(ficparo,"pop_based=%d\n",popbased);
Line 3799 int main(int argc, char *argv[])	Line 4064 int main(int argc, char *argv[])
}	}
ungetc(c,ficpar);	ungetc(c,ficpar);

fscanf(ficpar,"starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf\n",&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2);	fscanf(ficpar,"prevforecast=%d starting-proj-date=%lf/%lf/%lf final-proj-date=%lf/%lf/%lf mobil_average=%d\n",&prevfcast,&jproj1,&mproj1,&anproj1,&jproj2,&mproj2,&anproj2,&mobilavproj);
fprintf(ficparo,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);	fprintf(ficparo,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
fprintf(ficres,"starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf\n",jproj1,mproj1,anproj1,jproj2,mproj2,anproj2);	printf("prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
	fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
	fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
	/* day and month of proj2 are not used but only year anproj2.*/

while((c=getc(ficpar))=='#' && c!= EOF){	while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);	ungetc(c,ficpar);
Line 3816 int main(int argc, char *argv[])	Line 4083 int main(int argc, char *argv[])
fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);	fprintf(ficparo,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);
fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);	fprintf(ficres,"popforecast=%d popfile=%s popfiledate=%.lf/%.lf/%.lf last-popfiledate=%.lf/%.lf/%.lf\n",popforecast,popfile,jpyram,mpyram,anpyram,jpyram1,mpyram1,anpyram1);

	probs= ma3x(1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);	freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);

/------------ gnuplot -------------/	/------------ gnuplot -------------/
Line 3849 Interval (in months) between two waves:	Line 4117 Interval (in months) between two waves:
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n	- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
- Log file of the run: <a href=\"%s\">%s</a><br>\n	- Log file of the run: <a href=\"%s\">%s</a><br>\n
- Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);	- Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);
fclose(fichtm);	fclose(fichtm);

printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);	printinghtml(fileres,title,datafile, firstpass, lastpass, stepm, weightopt,model,imx,jmin,jmax,jmean,rfileres,popforecast,estepm,jprev1,mprev1,anprev1,jprev2,mprev2,anprev2);

Line 3862 Interval (in months) between two waves:	Line 4130 Interval (in months) between two waves:
free_imatrix(mw,1,lastpass-firstpass+1,1,imx);	free_imatrix(mw,1,lastpass-firstpass+1,1,imx);
free_ivector(num,1,n);	free_ivector(num,1,n);
free_vector(agedc,1,n);	free_vector(agedc,1,n);
free_matrix(covar,0,NCOVMAX,1,n);	/free_matrix(covar,0,NCOVMAX,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 3884 Interval (in months) between two waves:	Line 4152 Interval (in months) between two waves:
fprintf(ficrespl,"\n");	fprintf(ficrespl,"\n");

prlim=matrix(1,nlstate,1,nlstate);	prlim=matrix(1,nlstate,1,nlstate);
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
savms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldm=oldms; newm=newms; savm=savms; /* Keeps fixed addresses to free */

agebase=ageminpar;	agebase=ageminpar;
agelim=agemaxpar;	agelim=agemaxpar;
Line 3914 Interval (in months) between two waves:	Line 4177 Interval (in months) between two waves:

for (age=agebase; age<=agelim; age++){	for (age=agebase; age<=agelim; age++){
prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);	prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
fprintf(ficrespl,"%.0f",age );	fprintf(ficrespl,"%.0f ",age );
	for(j=1;j<=cptcoveff;j++)
	fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
fprintf(ficrespl," %.5f", prlim[i][i]);	fprintf(ficrespl," %.5f", prlim[i][i]);
fprintf(ficrespl,"\n");	fprintf(ficrespl,"\n");
Line 3942 Interval (in months) between two waves:	Line 4207 Interval (in months) between two waves:

/* hstepm=1; aff par mois*/	/* hstepm=1; aff par mois*/

	fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
for(cptcov=1,k=0;cptcov<=i1;cptcov++){	for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){	for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;	k=k+1;
Line 3959 Interval (in months) between two waves:	Line 4225 Interval (in months) between two waves:
p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);	p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
oldm=oldms;savm=savms;	oldm=oldms;savm=savms;
hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);	hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);
fprintf(ficrespij,"# Age");	fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate+ndeath;j++)	for(j=1; j<=nlstate+ndeath;j++)
fprintf(ficrespij," %1d-%1d",i,j);	fprintf(ficrespij," %1d-%1d",i,j);
fprintf(ficrespij,"\n");	fprintf(ficrespij,"\n");
for (h=0; h<=nhstepm; h++){	for (h=0; h<=nhstepm; h++){
fprintf(ficrespij,"%d %f %f",k,agedeb, agedeb+ hhstepm/YEARMstepm );	fprintf(ficrespij,"%d %3.f %3.f",k,agedeb, agedeb+ hhstepm/YEARMstepm );
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate+ndeath;j++)	for(j=1; j<=nlstate+ndeath;j++)
fprintf(ficrespij," %.5f", p3mat[i][j][h]);	fprintf(ficrespij," %.5f", p3mat[i][j][h]);
Line 3977 Interval (in months) between two waves:	Line 4243 Interval (in months) between two waves:
}	}
}	}

varprob(optionfilefiname, matcov, p, delti, nlstate, (int) bage, (int) fage,k,Tvar,nbcode, ncodemax);	varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax);

fclose(ficrespij);	fclose(ficrespij);


/---------- Forecasting ------------------/	/---------- Forecasting ------------------/
if((stepm == 1) && (strcmp(model,".")==0)){	/if((stepm == 1) && (strcmp(model,".")==0)){/
prevforecast(fileres, anproj1,mproj1,jproj1, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anproj2,p, i1);	if(prevfcast==1){
if (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);	/* if(stepm ==1){*/
}	prevforecast(fileres, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
else{	/* (popforecast==1) populforecast(fileres, anpyram,mpyram,jpyram, agemin,agemax, dateprev1, dateprev2,mobilav, agedeb, fage, popforecast, popfile, anpyram1,p, i1);*/
erreur=108;	/* } */
printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);	/* else{ */
fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model);	/* erreur=108; */
	/* printf("Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
	/* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
	/* } */
}	}


Line 4024 Interval (in months) between two waves:	Line 4293 Interval (in months) between two waves:
printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);	printf("Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);
fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);	fprintf(ficlog,"Computing Variance-covariance of DFLEs: file '%s' \n", fileresv);

calagedate=-1;	/* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
	prevalence(agemin, agemax, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
prevalence(ageminpar, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2, calagedate);	/* printf("ageminpar=%f, agemax=%f, s[lastpass][imx]=%d, agev[lastpass][imx]=%f, nlstate=%d, imx=%d, mint[lastpass][imx]=%f, anint[lastpass][imx]=%f,dateprev1=%f, dateprev2=%f, firstpass=%d, lastpass=%d\n",\
	ageminpar, agemax, s[lastpass][imx], agev[lastpass][imx], nlstate, imx, mint[lastpass][imx],anint[lastpass][imx], dateprev1, dateprev2, firstpass, lastpass);
	*/

if (mobilav!=0) {	if (mobilav!=0) {
mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);	mobaverage= ma3x(1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
Line 4150 Interval (in months) between two waves:	Line 4421 Interval (in months) between two waves:
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);

	free_matrix(covar,0,NCOVMAX,1,n);
free_matrix(matcov,1,npar,1,npar);	free_matrix(matcov,1,npar,1,npar);
free_vector(delti,1,npar);	/free_vector(delti,1,npar);/
	free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
free_matrix(agev,1,maxwav,1,imx);	free_matrix(agev,1,maxwav,1,imx);
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);	free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);	if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
	free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);

free_ivector(ncodemax,1,8);	free_ivector(ncodemax,1,8);
free_ivector(Tvar,1,15);	free_ivector(Tvar,1,15);
free_ivector(Tprod,1,15);	free_ivector(Tprod,1,15);
Line 4163 Interval (in months) between two waves:	Line 4438 Interval (in months) between two waves:
free_ivector(Tage,1,15);	free_ivector(Tage,1,15);
free_ivector(Tcode,1,100);	free_ivector(Tcode,1,100);

fprintf(fichtm,"\n</body>");	/* fclose(fichtm);*/
fclose(fichtm);	/* fclose(ficgp);/ / ALready done */
fclose(ficgp);


if(erreur >0){	if(erreur >0){

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

Removed from v.1.59
changed lines
	Added in v.1.74