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.60 and 1.80

version 1.60, 2002/11/19 00:17:38	version 1.80, 2003/06/05 15:34:14
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 ";	/* $Id$ */
	/* $Log$
	* Revision 1.80 2003/06/05 15:34:14 brouard
	* Trying to add the true revision is the program and log
	*
	/* Revision 1.79 2003/06/05 15:17:23 brouard
	/* * empty log message *
	/* */
	/* $Revision$ */
	/* $Date$ */
	/* $State$ */

	char version[80]="Imach version 0.95a1, June 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 105 double jmean; /* Mean space between 2 wa	Line 153 double jmean; /* Mean space between 2 wa
double oldm, newm, *savm; / Working pointers to matrices */	double oldm, newm, *savm; / Working pointers to matrices */
double oldms, newms, *savms; / Fixed working pointers to matrices */	double oldms, newms, *savms; / Fixed working pointers to matrices */
FILE fic,ficpar, ficparo,ficres, ficrespl, ficrespij, ficrest,ficresf,*ficrespop;	FILE fic,ficpar, ficparo,ficres, ficrespl, ficrespij, ficrest,ficresf,*ficrespop;
FILE *ficlog;	FILE ficlog, ficrespow;
FILE ficgp,ficresprob,ficpop, ficresprobcov, *ficresprobcor;	FILE ficgp,ficresprob,ficpop, ficresprobcov, *ficresprobcor;
FILE *ficresprobmorprev;	FILE *ficresprobmorprev;
FILE fichtm; / Html File */	FILE fichtm; / Html File */
Line 188 static int split( char path, char dirc	Line 236 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 297 void free_vector(double*v, int nl, int n	Line 339 void free_vector(double*v, int nl, int n
}	}

/**********************ivector *****************************/	/**********************ivector *****************************/
	char *cvector(long nl,long nh)
	{
	char *v;
	v=(char ) malloc((size_t)((nh-nl+1+NR_END)sizeof(char)));
	if (!v) nrerror("allocation failure in cvector");
	return v-nl+NR_END;
	}

	/****************free ivector ************************/
	void free_cvector(char *v, long nl, long nh)
	{
	free((FREE_ARG)(v+nl-NR_END));
	}

	/**********************ivector *****************************/
int *ivector(long nl,long nh)	int *ivector(long nl,long nh)
{	{
int *v;	int *v;
Line 365 double **matrix(long nrl, long nrh, long	Line 422 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 463 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 612 void powell(double p[], double **xi, int	Line 674 void powell(double p[], double **xi, int
del=0.0;	del=0.0;
printf("\nPowell iter=%d -2LL=%.12f",iter,*fret);	printf("\nPowell iter=%d -2LL=%.12f",iter,*fret);
fprintf(ficlog,"\nPowell iter=%d -2LL=%.12f",iter,*fret);	fprintf(ficlog,"\nPowell iter=%d -2LL=%.12f",iter,*fret);
for (i=1;i<=n;i++)	fprintf(ficrespow,"%d %.12f",iter,fret);
	for (i=1;i<=n;i++) {
printf(" %d %.12f",i, p[i]);	printf(" %d %.12f",i, p[i]);
fprintf(ficlog," %d %.12f",i, p[i]);	fprintf(ficlog," %d %.12lf",i, p[i]);
	fprintf(ficrespow," %.12lf", p[i]);
	}
printf("\n");	printf("\n");
fprintf(ficlog,"\n");	fprintf(ficlog,"\n");
	fprintf(ficrespow,"\n");
for (i=1;i<=n;i++) {	for (i=1;i<=n;i++) {
for (j=1;j<=n;j++) xit[j]=xi[j][i];	for (j=1;j<=n;j++) xit[j]=xi[j][i];
fptt=(*fret);	fptt=(*fret);
Line 856 double matprod2(double out, double *	Line 922 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 985 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 996 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 1216 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;
	char filerespow[FILENAMELENGTH];
xi=matrix(1,npar,1,npar);	xi=matrix(1,npar,1,npar);
for (i=1;i<=npar;i++)	for (i=1;i<=npar;i++)
for (j=1;j<=npar;j++)	for (j=1;j<=npar;j++)
xi[i][j]=(i==j ? 1.0 : 0.0);	xi[i][j]=(i==j ? 1.0 : 0.0);
printf("Powell\n"); fprintf(ficlog,"Powell\n");	printf("Powell\n"); fprintf(ficlog,"Powell\n");
	strcpy(filerespow,"pow");
	strcat(filerespow,fileres);
	if((ficrespow=fopen(filerespow,"w"))==NULL) {
	printf("Problem with resultfile: %s\n", filerespow);
	fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
	}
	fprintf(ficrespow,"# Powell\n# iter -2*LL");
	for (i=1;i<=nlstate;i++)
	for(j=1;j<=nlstate+ndeath;j++)
	if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
	fprintf(ficrespow,"\n");
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));	fclose(ficrespow);
fprintf(ficlog,"#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,"\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 1502 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 1522 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 1537 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 1557 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 1585 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 1617 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 1632 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 1650 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 1661 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 1696 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 1708 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 1548 void concatwav(int wav[], int **dh, int	Line 1787 void concatwav(int wav[], int **dh, int
wav[i]=mi;	wav[i]=mi;
if(mi==0){	if(mi==0){
if(first==0){	if(first==0){
printf("Warning, no any valid information for:%d line=%d and may be others, see log file\n",num[i],i);	printf("Warning! None valid information for:%d line=%d (skipped) and may be others, see log file\n",num[i],i);
first=1;	first=1;
}	}
if(first==1){	if(first==1){
fprintf(ficlog,"Warning, no any valid information for:%d line=%d\n",num[i],i);	fprintf(ficlog,"Warning! None valid information for:%d line=%d (skipped)\n",num[i],i);
}	}
} /* end mi==0 */	} /* end mi==0 */
}	} /* End individuals */

for(i=1; i<=imx; i++){	for(i=1; i<=imx; i++){
for(mi=1; mi<wav[i];mi++){	for(mi=1; mi<wav[i];mi++){
if (stepm <=0)	if (stepm <=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) { /* A death */
if (agedc[i] < 2*AGESUP) {	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 */
Line 1570 void concatwav(int wav[], int **dh, int	Line 1809 void concatwav(int wav[], int **dh, int
if (j >= jmax) jmax=j;	if (j >= jmax) jmax=j;
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);*/
	if(j<0)printf("Error! Negative delay (%d to death) between waves %d and %d of individual %d at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],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]);/
	if(j<0)printf("Error! Negative delay (%d) between waves %d and %d of individual %d at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],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 1936 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 1956 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 1985 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 1768 void evsij(char fileres[], double ***eij	Line 2024 void evsij(char fileres[], double ***eij
for(i=1;i<=nlstate;i++){	for(i=1;i<=nlstate;i++){
cptj=cptj+1;	cptj=cptj+1;
for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){	for(h=0, gm[h][cptj]=0.; h<=nhstepm-1; h++){

gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;	gm[h][cptj] = (p3mat[i][j][h]+p3mat[i][j][h+1])/2.;
}	}
}	}
}	}
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 2037 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 2079 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 2147 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 2169 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 2204 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 2220 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 2242 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 2273 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 2304 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 2320 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 2328 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 2341 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 2380 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 2402 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 2548 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 2591 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 2607 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 2626 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 2639 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 2651 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 2565 fprintf(fichtm," \n<ul><li><b>Graphs</b>	Line 2838 fprintf(fichtm," \n<ul><li><b>Graphs</b>
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");	fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}	}
/* Pij */	/* Pij */
fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before: pe%s%d1.png<br>	fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i, %d (stepm) months before: pe%s%d1.png<br>
<img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);	<img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
/* Quasi-incidences */	/* Quasi-incidences */
fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>	fprintf(fichtm,"<br>- Pij or Conditional probabilities to be observed in state j being in state i %d (stepm) months before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too: pe%s%d2.png<br>
Line 2596 health expectancies in states (1) and (2	Line 2869 health expectancies in states (1) and (2
- Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n	- Health expectancies with their variances (no covariance): <a href=\"t%s\">t%s</a> <br>\n
- Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);	- Standard deviation of stable prevalences: <a href=\"vpl%s\">vpl%s</a> <br>\n",rfileres,rfileres,fileres,fileres,fileres,fileres,fileres,fileres, estepm, fileres,fileres,fileres,fileres,fileres,fileres);

if(popforecast==1) fprintf(fichtm,"\n	/* if(popforecast==1) fprintf(fichtm,"\n */
- Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n	/* - Prevalences forecasting: <a href=\"f%s\">f%s</a> <br>\n */
- Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n	/* - Population forecasting (if popforecast=1): <a href=\"pop%s\">pop%s</a> <br>\n */
<br>",fileres,fileres,fileres,fileres);	/* <br>",fileres,fileres,fileres,fileres); */
else	/* else */
fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model);	/* fprintf(fichtm,"\n No population forecast: popforecast = %d (instead of 1) or stepm = %d (instead of 1) or model=%s (instead of .)<br><br></li>\n",popforecast, stepm, model); */
fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");	fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");

m=cptcoveff;	m=cptcoveff;
Line 2618 fprintf(fichtm," <ul><li><b>Graphs</b></	Line 2891 fprintf(fichtm," <ul><li><b>Graphs</b></
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");	fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}	}
for(cpt=1; cpt<=nlstate;cpt++) {	for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident	fprintf(fichtm,"<br>- Observed and period prevalence (with confident
interval) in state (%d): v%s%d%d.png <br>	interval) in state (%d): v%s%d%d.png <br>
<img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);	<img src=\"v%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
}	}
Line 2653 m=pow(2,cptcoveff);	Line 2926 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 2719 m=pow(2,cptcoveff);	Line 2992 m=pow(2,cptcoveff);
}	}
}	}

/* CV preval stat */	/* CV preval stable (period) */
for (k1=1; k1<= m ; k1 ++) {	for (k1=1; k1<= m ; k1 ++) {
for (cpt=1; cpt<nlstate ; cpt ++) {	for (cpt=1; cpt<=nlstate ; cpt ++) {
k=3;	k=3;
fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);	fprintf(ficgp,"\nset out \"p%s%d%d.png\" \n",strtok(optionfile, "."),cpt,k1);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter png small\nset size 0.65,0.65\nplot [%.f:%.f] \"pij%s\" u ($1==%d ? ($3):1/0):($%d/($%d",ageminpar,agemaxpar,fileres,k1,k+cpt+1,k+1);

for (i=1; i< nlstate ; i ++)	for (i=1; i<= nlstate ; i ++)
fprintf(ficgp,"+$%d",k+i+1);	fprintf(ficgp,"+$%d",k+i+1);
fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);	fprintf(ficgp,")) t\"prev(%d,%d)\" w l",cpt,cpt+1);

Line 2849 int movingaverage(double ***probs, doubl	Line 3122 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 3161 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 3176 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");

	/* if (h==(int)(YEARMyearp)){ /
	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 (hhstepm/YEARMstepm ==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 (hhstepm/YEARMstepm== yearp) {
}	fprintf(ficresf," %.3f", p3mat[i][j][h]);
if (h==(int)(calagedate+12*cpt)){	}
fprintf(ficresf," %.3f", kk1);	} /* end i */
	if (hhstepm/YEARMstepm==yearp) {
	fprintf(ficresf," %.3f", ppij);
}	}
}	}/* end j */
}	} /* end h */
free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);	free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
}	} /* end agec */
}	} /* end yearp */
}	} /* end cptcod */
}	} /* end cptcov */

if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);	if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);

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 3259 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 3307 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 3322 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 3331 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 3343 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 3354 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 3368 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 3376 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 3384 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 3411 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 3431 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 3450 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 3597 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 3614 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 3736 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 3535 int main(int argc, char *argv[])	Line 3832 int main(int argc, char *argv[])

for (i=1; i<=imx; i++) {	for (i=1; i<=imx; i++) {
for(m=2; (m<= maxwav); m++) {	for(m=2; (m<= maxwav); m++) {
if ((mint[m][i]== 99) && (s[m][i] <= nlstate)){	if (((int)mint[m][i]== 99) && (s[m][i] <= nlstate)){
anint[m][i]=9999;	anint[m][i]=9999;
s[m][i]=-1;	s[m][i]=-1;
}	}
if(moisdc[i]==99 && andc[i]==9999 & s[m][i]>nlstate) s[m][i]=-1;	if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
	printf("Error! Date of death (month %2d and year %4d) of individual %d on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
	fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %d on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);
	s[m][i]=-1;
	}
	if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
	printf("Error! Month of death of individual %d on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]);
	fprintf(ficlog,"Error! Month of death of individual %d on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]);
	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=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)
if(moisdc[i]!=99 && andc[i]!=9999)	if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
agev[m][i]=agedc[i];	agev[m][i]=agedc[i];
/if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;/	/if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;/
else {	else {
if (andc[i]!=9999){	if ((int)andc[i]!=9999){
printf("Warning negative age at death: %d line:%d\n",num[i],i);	printf("Warning negative age at death: %d line:%d\n",num[i],i);
fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);	fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
agev[m][i]=-1;	agev[m][i]=-1;
}	}
}	}
}	}
else if(s[m][i] !=9){ /* Should no more exist */	else if(s[m][i] !=9){ /* 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((int)mint[m][i]==99 \|\| (int)anint[m][i]==9999)
agev[m][i]=1;	agev[m][i]=1;
else if(agev[m][i] <agemin){	else if(agev[m][i] <agemin){
agemin=agev[m][i];	agemin=agev[m][i];
Line 3586 int main(int argc, char *argv[])	Line 3894 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 3903 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 3663 int main(int argc, char *argv[])	Line 3978 int main(int argc, char *argv[])
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 3779 int main(int argc, char *argv[])	Line 4094 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 3789 int main(int argc, char *argv[])	Line 4106 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 3804 int main(int argc, char *argv[])	Line 4121 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 3821 int main(int argc, char *argv[])	Line 4140 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 3848 int main(int argc, char *argv[])	Line 4168 int main(int argc, char *argv[])
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n	Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
\n	\n
Total number of observations=%d <br>\n	Total number of observations=%d <br>\n
	Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n	Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
<hr size=\"2\" color=\"#EC5E5E\">	<hr size=\"2\" color=\"#EC5E5E\">
<ul><li><h4>Parameter files</h4>\n	<ul><li><h4>Parameter files</h4>\n
- 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,agemin,agemax,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 3867 Interval (in months) between two waves:	Line 4188 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 3914 Interval (in months) between two waves:	Line 4235 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 4265 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 4283 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 4301 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 4351 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 4479 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 4496 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){
Line 4194 Interval (in months) between two waves:	Line 4526 Interval (in months) between two waves:
strcpy(plotcmd,GNUPLOTPROGRAM);	strcpy(plotcmd,GNUPLOTPROGRAM);
strcat(plotcmd," ");	strcat(plotcmd," ");
strcat(plotcmd,optionfilegnuplot);	strcat(plotcmd,optionfilegnuplot);
printf("Starting: %s\n",plotcmd);fflush(stdout);	printf("Starting graphs with: %s",plotcmd);fflush(stdout);
system(plotcmd);	system(plotcmd);
	printf(" Wait...");

/#ifdef windows/	/#ifdef windows/
while (z[0] != 'q') {	while (z[0] != 'q') {

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

Removed from v.1.60
changed lines
	Added in v.1.80