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.268 and 1.282

version 1.268, 2017/05/18 20:09:32	version 1.282, 2018/02/27 22:50:02
Line 1	Line 1
/* $Id$	/* $Id$
$State$	$State$
$Log$	$Log$
	Revision 1.282 2018/02/27 22:50:02 brouard
	* empty log message *

	Revision 1.281 2018/02/27 19:25:23 brouard
	Summary: Adding second argument for quitting

	Revision 1.280 2018/02/21 07:58:13 brouard
	Summary: 0.99r15

	New Makefile with recent VirtualBox 5.26. Bug in sqrt negatve in imach.c

	Revision 1.279 2017/07/20 13:35:01 brouard
	Summary: temporary working

	Revision 1.278 2017/07/19 14:09:02 brouard
	Summary: Bug for mobil_average=0 and prevforecast fixed(?)

	Revision 1.277 2017/07/17 08:53:49 brouard
	Summary: BOM files can be read now

	Revision 1.276 2017/06/30 15:48:31 brouard
	Summary: Graphs improvements

	Revision 1.275 2017/06/30 13:39:33 brouard
	Summary: Saito's color

	Revision 1.274 2017/06/29 09:47:08 brouard
	Summary: Version 0.99r14

	Revision 1.273 2017/06/27 11:06:02 brouard
	Summary: More documentation on projections

	Revision 1.272 2017/06/27 10:22:40 brouard
	Summary: Color of backprojection changed from 6 to 5(yellow)

	Revision 1.271 2017/06/27 10:17:50 brouard
	Summary: Some bug with rint

	Revision 1.270 2017/05/24 05:45:29 brouard
	* empty log message *

	Revision 1.269 2017/05/23 08:39:25 brouard
	Summary: Code into subroutine, cleanings

Revision 1.268 2017/05/18 20:09:32 brouard	Revision 1.268 2017/05/18 20:09:32 brouard
Summary: backprojection and confidence intervals of backprevalence	Summary: backprojection and confidence intervals of backprevalence

Line 1084 FILE *ficrescveij;	Line 1128 FILE *ficrescveij;
char filerescve[FILENAMELENGTH];	char filerescve[FILENAMELENGTH];
FILE *ficresvij;	FILE *ficresvij;
char fileresv[FILENAMELENGTH];	char fileresv[FILENAMELENGTH];
FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];
FILE *ficresvbl;
char fileresvbl[FILENAMELENGTH];
char title[MAXLINE];	char title[MAXLINE];
char model[MAXLINE]; /*< The model line /	char model[MAXLINE]; /*< The model line /
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH];	char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH], fileresplb[FILENAMELENGTH];
Line 2489 void powell(double p[], double **xi, int	Line 2530 void powell(double p[], double **xi, int

double prevalim(double prlim, int nlstate, double x[], double age, double oldm, double savm, double ftolpl, int *ncvyear, int ij, int nres)	double prevalim(double prlim, int nlstate, double x[], double age, double oldm, double savm, double ftolpl, int *ncvyear, int ij, int nres)
{	{
/* Computes the prevalence limit in each live state at age x and for covariate combination ij	/**< Computes the prevalence limit in each live state at age x and for covariate combination ij
(and selected quantitative values in nres)	* (and selected quantitative values in nres)
by left multiplying the unit	* by left multiplying the unit
matrix by transitions matrix until convergence is reached with precision ftolpl */	* matrix by transitions matrix until convergence is reached with precision ftolpl
/* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I */	* Wx= Wx-1 Px-1= Wx-2 Px-2 Px-1 = Wx-n Px-n ... Px-2 Px-1 I
/* Wx is row vector: population in state 1, population in state 2, population dead */	* Wx is row vector: population in state 1, population in state 2, population dead
/* or prevalence in state 1, prevalence in state 2, 0 */	* or prevalence in state 1, prevalence in state 2, 0
/* newm is the matrix after multiplications, its rows are identical at a factor */	* newm is the matrix after multiplications, its rows are identical at a factor.
/* Initial matrix pimij */	* Inputs are the parameter, age, a tolerance for the prevalence limit ftolpl.
	* Output is prlim.
	* Initial matrix pimij
	*/
/* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */	/* {0.85204250825084937, 0.13044499163996345, 0.017512500109187184, */
/* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */	/* 0.090851990222114765, 0.88271245433047185, 0.026435555447413338, */
/* 0, 0 , 1} */	/* 0, 0 , 1} */
Line 2932 double pmij(double ps, double *cov,	Line 2976 double pmij(double ps, double *cov,
/* Diag(w_x) */	/* Diag(w_x) */
/* Problem with prevacurrent which can be zero */	/* Problem with prevacurrent which can be zero */
sumnew=0.;	sumnew=0.;
/* for (ii=1;ii<=nlstate+ndeath;ii++){ */	/for (ii=1;ii<=nlstate+ndeath;ii++){/
for (ii=1;ii<=nlstate;ii++){ /* Only on live states */	for (ii=1;ii<=nlstate;ii++){ /* Only on live states */
	/* printf(" agefin=%d, ii=%d, ij=%d, prev=%f\n",(int)agefin,ii, ij, prevacurrent[(int)agefin][ii][ij]); */
sumnew+=prevacurrent[(int)agefin][ii][ij];	sumnew+=prevacurrent[(int)agefin][ii][ij];
}	}
if(sumnew >0.01){ /* At least some value in the prevalence */	if(sumnew >0.01){ /* At least some value in the prevalence */
for (ii=1;ii<=nlstate+ndeath;ii++){	for (ii=1;ii<=nlstate+ndeath;ii++){
for (j=1;j<=nlstate+ndeath;j++)	for (j=1;j<=nlstate+ndeath;j++)
doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0);	doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0);
}	}
}else{	}else{
for (ii=1;ii<=nlstate+ndeath;ii++){	for (ii=1;ii<=nlstate+ndeath;ii++){
Line 3239 double *hbxij(double *po, int nhstep	Line 3284 double *hbxij(double *po, int nhstep
newm=savm;	newm=savm;
/* Covariates have to be included here again */	/* Covariates have to be included here again */
cov[1]=1.;	cov[1]=1.;
agexact=age-((h-1)hstepm + (d))stepm/YEARM; /* age just before transition, d or d-1? */	agexact=age-( (h-1)hstepm + (d) )stepm/YEARM; /* age just before transition, d or d-1? */
/* agexact=age+((h-1)hstepm + (d-1))stepm/YEARM; /\* age just before transition \/ /	/* agexact=age+((h-1)hstepm + (d-1))stepm/YEARM; /\* age just before transition \/ /
cov[2]=agexact;	cov[2]=agexact;
if(nagesqr==1)	if(nagesqr==1)
Line 3841 void likelione(FILE *ficres,double p[],	Line 3886 void likelione(FILE *ficres,double p[],
else if(mle >=1)	else if(mle >=1)
fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);	fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);
fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));	fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
	fprintf(fichtm,"\n<br>Equation of the model: <b>model=1+age+%s</b><br>\n",model);

for (k=1; k<= nlstate ; k++) {	for (k=1; k<= nlstate ; k++) {
fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \	fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
Line 5463 void concatwav(int wav[], int **dh, int	Line 5508 void concatwav(int wav[], int **dh, int
/* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.	/* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.
nhstepm is the number of hstepm from age to agelim	nhstepm is the number of hstepm from age to agelim
nstepm is the number of stepm from age to agelin.	nstepm is the number of stepm from age to agelin.
Look at hpijx to understand the reason of that which relies in memory size	Look at hpijx to understand the reason which relies in memory size consideration
and note for a fixed period like estepm months */	and note for a fixed period like estepm months */
/* 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
Line 5694 void concatwav(int wav[], int **dh, int	Line 5739 void concatwav(int wav[], int **dh, int
/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/	/* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/

}	}

	/* Standard deviation of expectancies ij */
fprintf(ficresstdeij,"%3.0f",age );	fprintf(ficresstdeij,"%3.0f",age );
for(i=1; i<=nlstate;i++){	for(i=1; i<=nlstate;i++){
eip=0.;	eip=0.;
Line 5709 void concatwav(int wav[], int **dh, int	Line 5755 void concatwav(int wav[], int **dh, int
}	}
fprintf(ficresstdeij,"\n");	fprintf(ficresstdeij,"\n");

	/* Variance of expectancies ij */
fprintf(ficrescveij,"%3.0f",age );	fprintf(ficrescveij,"%3.0f",age );
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
for(j=1; j<=nlstate;j++){	for(j=1; j<=nlstate;j++){
Line 5742 void concatwav(int wav[], int **dh, int	Line 5789 void concatwav(int wav[], int **dh, int
/********** Variance ****************/	/********** Variance ****************/
void varevsij(char optionfilefiname[], double *vareij, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double *prlim, double ftolpl, int ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres)	void varevsij(char optionfilefiname[], double *vareij, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double *prlim, double ftolpl, int ncvyearp, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[], int nres)
{	{
/* Variance of health expectancies */	/** Variance of health expectancies
/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double * savm,double ftolpl);*/	* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double * savm,double ftolpl);
/* double *newm;/	* double **newm;
/* int movingaverage(double *probs, double bage,double fage, double mobaverage, int mobilav)/	* int movingaverage(double *probs, double bage,double fage, double *mobaverage, int mobilav)
	*/

/* int movingaverage(); */	/* int movingaverage(); */
double dnewm,doldm;	double dnewm,doldm;
Line 5753 void concatwav(int wav[], int **dh, int	Line 5801 void concatwav(int wav[], int **dh, int
int i, j, nhstepm, hstepm, h, nstepm ;	int i, j, nhstepm, hstepm, h, nstepm ;
int k;	int k;
double *xp;	double *xp;
double gp, gm; /* for var eij */	double gp, gm; /*< for var eij /
double *gradg, trgradg; /for var eij */	double *gradg, trgradg; /< for var eij /
double gradgp, trgradgp; /* for var p point j */	double gradgp, trgradgp; /*< for var p point j /
double gpp, gmp; /* for var p point j */	double gpp, gmp; /*< for var p point j /
double *varppt; / for var p point j nlstate to nlstate+ndeath */	double varppt; /< for var p point j nlstate to nlstate+ndeath */
double ***p3mat;	double ***p3mat;
double age,agelim, hf;	double age,agelim, hf;
/* double **mobaverage; /	/* double **mobaverage; /
Line 5818 void concatwav(int wav[], int **dh, int	Line 5866 void concatwav(int wav[], int **dh, int
/* fprintf(fichtm, "#Local time at start: %s", strstart);*/	/* fprintf(fichtm, "#Local time at start: %s", strstart);*/
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<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 <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);
pstamp(ficresvij);	pstamp(ficresvij);
fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are ");	fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are ");
Line 5873 void concatwav(int wav[], int **dh, int	Line 5921 void concatwav(int wav[], int **dh, int
for(i=1; i<=npar; i++){ /* Computes gradient x + delta*/	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);
}	}
	/**< Computes the prevalence limit with parameter theta shifted of delta up to ftolpl precision and
	* returns into prlim .
	*/
prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres);	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres);

	/* If popbased = 1 we use crossection prevalences. Previous step is useless but prlim is created */
if (popbased==1) {	if (popbased==1) {
if(mobilav ==0){	if(mobilav ==0){
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
Line 5885 void concatwav(int wav[], int **dh, int	Line 5936 void concatwav(int wav[], int **dh, int
prlim[i][i]=mobaverage[(int)age][i][ij];	prlim[i][i]=mobaverage[(int)age][i][ij];
}	}
}	}
	/**< Computes the shifted transition matrix \f$ {}{h}_p^{ij}_x\f$ at horizon h.
hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=1 to nhstepm */	*/
	hpxij(p3mat,nhstepm,age,hstepm,xp,nlstate,stepm,oldm,savm, ij,nres); /* Returns p3mat[i][j][h] for h=0 to nhstepm */
	/**< And for each alive state j, sums over i \f$ w^i_x {}{h}_p^{ij}_x\f$, which are the probability
	* at horizon h in state j including mortality.
	*/
for(j=1; j<= nlstate; j++){	for(j=1; j<= nlstate; j++){
for(h=0; h<=nhstepm; h++){	for(h=0; h<=nhstepm; h++){
for(i=1, gp[h][j]=0.;i<=nlstate;i++)	for(i=1, gp[h][j]=0.;i<=nlstate;i++)
gp[h][j] += prlim[i][i]*p3mat[i][j][h];	gp[h][j] += prlim[i][i]*p3mat[i][j][h];
}	}
}	}
/* Next for computing probability of death (h=1 means	/* Next for computing shifted+ probability of death (h=1 means
computed over hstepm matrices product = hstepm*stepm months)	computed over hstepm matrices product = hstepm*stepm months)
as a weighted average of prlim.	as a weighted average of prlim(i) * p(i,j) p.3=w1p13 + w2p23 .
*/	*/
for(j=nlstate+1;j<=nlstate+ndeath;j++){	for(j=nlstate+1;j<=nlstate+ndeath;j++){
for(i=1,gpp[j]=0.; i<= nlstate; i++)	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 probability of death */
	/* Again with minus shift */

for(i=1; i<=npar; i++) /* Computes gradient x - delta */	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);
Line 5934 void concatwav(int wav[], int **dh, int	Line 5990 void concatwav(int wav[], int **dh, int
for(i=1,gmp[j]=0.; i<= nlstate; i++)	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 probability of death */	/* end shifting computations */

	/**< Computing gradient matrix at horizon h
	*/
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];
}	}
	/**< Gradient of overall mortality p.3 (or p.j)
for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu */	*/
	for(j=nlstate+1; j<= nlstate+ndeath; j++){ /* var mu mortality from j */
gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];	gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
}	}

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

	/* We got the gradient matrix for each theta and state j */
trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */	trgradg =ma3x(0,nhstepm,1,nlstate,1,npar); /* veij */

for(h=0; h<=nhstepm; h++) /* veij */	for(h=0; h<=nhstepm; h++) /* veij */
Line 5957 void concatwav(int wav[], int **dh, int	Line 6017 void concatwav(int wav[], int **dh, int
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];
	/**< as well as its transposed matrix
	*/

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++)
for(j=1;j<=nlstate;j++)	for(j=1;j<=nlstate;j++)
vareij[i][j][(int)age] =0.;	vareij[i][j][(int)age] =0.;

	/* Computing trgradg by matcov by gradg at age and summing over h
	* and k (nhstepm) formula 15 of article
	* Lievre-Brouard-Heathcote
	*/

for(h=0;h<=nhstepm;h++){	for(h=0;h<=nhstepm;h++){
for(k=0;k<=nhstepm;k++){	for(k=0;k<=nhstepm;k++){
matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);	matprod2(dnewm,trgradg[h],1,nlstate,1,npar,1,npar,matcov);
Line 5974 void concatwav(int wav[], int **dh, int	Line 6040 void concatwav(int wav[], int **dh, int
}	}
}	}

/* pptj */	/* pptj is p.3 or p.j = trgradgp by cov by gradgp, variance of
	* p.j overall mortality formula 49 but computed directly because
	* we compute the grad (wix pijx) instead of grad (pijx),even if
	* wix is independent of theta.
	*/
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);
for(j=nlstate+1;j<=nlstate+ndeath;j++)	for(j=nlstate+1;j<=nlstate+ndeath;j++)
Line 6062 void concatwav(int wav[], int **dh, int	Line 6132 void concatwav(int wav[], int **dh, int
} /* end varevsij */	} /* end varevsij */

/********** 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 ncvyearp, int ij, char strstart[], int nres)	void varprevlim(char fileresvpl[], FILE ficresvpl, 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 ncvyearp, int ij, char strstart[], int nres)
{	{
/* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/	/* Variance of prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double savm,double ftolpl);/	/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double savm,double ftolpl);/
Line 6187 void concatwav(int wav[], int **dh, int	Line 6257 void concatwav(int wav[], int **dh, int


/********** Variance of backprevalence limit ****************/	/********** Variance of backprevalence limit ****************/
void varbrevlim(char fileres[], double varbpl, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double *bprlim, double ftolpl, int mobilavproj, int ncvyearp, int ij, char strstart[], int nres)	void varbrevlim(char fileresvbl[], FILE ficresvbl, double varbpl, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double bprlim, double ftolpl, int mobilavproj, int ncvyearp, int ij, char strstart[], int nres)
{	{
/* Variance of backward prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/	/* Variance of backward prevalence limit for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double savm,double ftolpl);/	/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double savm,double ftolpl);/
Line 6583 To be simple, these graphs help to under	Line 6653 To be simple, these graphs help to under
}	}

/* Eigen vectors */	/* Eigen vectors */
v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));	if(1+(v1-lc1)*(v1-lc1)/cv12/cv12 <1.e-5){
	printf(" Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12);
	fprintf(ficlog," Error sqrt of a negative number: %lf\n",1+(v1-lc1)*(v1-lc1)/cv12/cv12);
	v11=(1./sqrt(fabs(1+(v1-lc1)*(v1-lc1)/cv12/cv12)));
	}else
	v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
/v21=sqrt(1.-v11v11); // error */	/v21=sqrt(1.-v11v11); // error */
v21=(lc1-v1)/cv12*v11;	v21=(lc1-v1)/cv12*v11;
v12=-v21;	v12=-v21;
Line 6614 To be simple, these graphs help to under	Line 6689 To be simple, these graphs help to under
fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);	fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu1,mu2);
fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);	fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k1,l1,k2,l2);
fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f(%11.3e%11.3ecos(t)+%11.3e%11.3esin(t)), %11.3e +%.3f(%11.3e%11.3ecos(t)+%11.3e%11.3esin(t)) not", \	fprintf(ficgp,"\nplot [-pi:pi] %11.3e+ %.3f(%11.3e%11.3ecos(t)+%11.3e%11.3esin(t)), %11.3e +%.3f(%11.3e%11.3ecos(t)+%11.3e%11.3esin(t)) not", \
mu1,std,v11,sqrt(lc1),v12,sqrt(fabs(lc2)), \	mu1,std,v11,sqrt(fabs(lc1)),v12,sqrt(fabs(lc2)), \
mu2,std,v21,sqrt(lc1),v22,sqrt(fabs(lc2))); /* For gnuplot only */	mu2,std,v21,sqrt(fabs(lc1)),v22,sqrt(fabs(lc2))); /* For gnuplot only */
}else{	}else{
first=0;	first=0;
fprintf(fichtmcov," %d (%.3f),",(int) age, c12);	fprintf(fichtmcov," %d (%.3f),",(int) age, c12);
Line 6652 void printinghtml(char fileresu[], char	Line 6727 void printinghtml(char fileresu[], char
int lastpass, int stepm, int weightopt, char model[],\	int lastpass, int stepm, int weightopt, char model[],\
int imx,int jmin, int jmax, double jmeanint,char rfileres[],\	int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \	int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \
double jprev1, double mprev1,double anprev1, double dateprev1, \	double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \
double jprev2, double mprev2,double anprev2, double dateprev2){	double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){
int jj1, k1, i1, cpt, k4, nres;	int jj1, k1, i1, cpt, k4, nres;

fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \	fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
Line 6807 divided by h: <sub>h</sub>P<sub>ij</sub>	Line 6882 divided by h: <sub>h</sub>P<sub>ij</sub>
if(prevfcast==1){	if(prevfcast==1){
/* Projection of prevalence up to period (stable) prevalence in each health state */	/* Projection of prevalence up to period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \	fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);	<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
}	}
}	}
if(backcast==1){	if(backcast==1){
/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */	/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to stable (mixed) back prevalence in state %d. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \	fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);	from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \
	account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
	with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
	<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
}	}
}	}

Line 6921 true period expectancies (those weighted	Line 6999 true period expectancies (those weighted
}	}

/***************** Gnuplot file ************/	/***************** Gnuplot file ************/
void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){	void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){

char dirfileres[132],optfileres[132];	char dirfileres[132],optfileres[132];
char gplotcondition[132], gplotlabel[132];	char gplotcondition[132], gplotlabel[132];
Line 6930 void printinggnuplot(char fileresu[], ch	Line 7008 void printinggnuplot(char fileresu[], ch
int ng=0;	int ng=0;
int vpopbased;	int vpopbased;
int ioffset; /* variable offset for columns */	int ioffset; /* variable offset for columns */
	int iyearc=1; /* variable column for year of projection */
	int iagec=1; /* variable column for age of projection */
int nres=0; /* Index of resultline */	int nres=0; /* Index of resultline */
int istart=1; /* For starting graphs in projections */	int istart=1; /* For starting graphs in projections */

Line 6943 void printinggnuplot(char fileresu[], ch	Line 7023 void printinggnuplot(char fileresu[], ch
/#endif /	/#endif /
m=pow(2,cptcoveff);	m=pow(2,cptcoveff);

	/* diagram of the model */
	fprintf(ficgp,"\n#Diagram of the model \n");
	fprintf(ficgp,"\ndelta=0.03;delta2=0.07;unset arrow;\n");
	fprintf(ficgp,"yoff=(%d > 2? 0:1);\n",nlstate);
	fprintf(ficgp,"\n#Peripheral arrows\nset for [i=1:%d] for [j=1:%d] arrow i10+j from cos(pi((1-(%d/2)2./%d)/2+(i-1)2./%d))-(i!=j?(i-j)/abs(i-j)delta:0), yoff +sin(pi((1-(%d/2)2./%d)/2+(i-1)2./%d)) + (i!=j?(i-j)/abs(i-j)delta:0) rto -0.95(cos(pi((1-(%d/2)2./%d)/2+(i-1)2./%d))+(i!=j?(i-j)/abs(i-j)delta:0) - cos(pi((1-(%d/2)2./%d)/2+(j-1)2./%d)) + (i!=j?(i-j)/abs(i-j)delta2:0)), -0.95(sin(pi((1-(%d/2)2./%d)/2+(i-1)2./%d)) + (i!=j?(i-j)/abs(i-j)delta:0) - sin(pi((1-(%d/2)2./%d)/2+(j-1)2./%d))+( i!=j?(i-j)/abs(i-j)*delta2:0)) ls (i < j? 1:2)\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);

	fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] arrow (%d+1)10+i from cos(pi((1-(%d/2)2./%d)/2+(1-i)2./%d))-(i!=j?(i-j)/abs(i-j)delta:0), yoff +sin(pi((1-(%d/2)2./%d)/2+(1-i)2./%d)) + (i!=j?(i-j)/abs(i-j)delta:0) rto -0.80(cos(pi((1-(%d/2)2./%d)/2+(1-i)2./%d))+(i!=j?(i-j)/abs(i-j)delta:0) ), -0.80(sin(pi((1-(%d/2)2./%d)/2+(1-i)2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
	fprintf(ficgp,"\n#show arrow\nunset label\n");
	fprintf(ficgp,"\n#States labels, starting from 2 (2-i) instead of (1-i), was (i-1)\nset for [i=1:%d] label i sprintf(\"State %%d\",i) center at cos(pi((1-(%d/2)2./%d)/2+(2-i)2./%d)), yoff+sin(pi((1-(%d/2)2./%d)/2+(2-i)2./%d)) font \"helvetica, 16\" tc rgbcolor \"blue\"\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
	fprintf(ficgp,"\nset label %d+1 sprintf(\"State %%d\",%d+1) center at 0.,0. font \"helvetica, 16\" tc rgbcolor \"red\"\n",nlstate,nlstate);
	fprintf(ficgp,"\n#show label\nunset border;unset xtics; unset ytics;\n");
	fprintf(ficgp,"\n\nset ter svg size 640, 480;set out \"%s_.svg\" \n",subdirf2(optionfilefiname,"D_"));
	fprintf(ficgp,"unset log y; plot [-1.2:1.2][yoff-1.2:1.2] 1/0 not; set out;reset;\n");

/* Contribution to likelihood */	/* Contribution to likelihood */
/* Plot the probability implied in the likelihood */	/* Plot the probability implied in the likelihood */
fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");	fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
Line 7012 void printinggnuplot(char fileresu[], ch	Line 7106 void printinggnuplot(char fileresu[], ch

fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);	fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel);	/* fprintf(ficgp,"set label \"Alive state %d %s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",cpt,gplotlabel); */
	fprintf(ficgp,"set title \"Alive state %d %s\" font \"Helvetica,12\"\n",cpt,gplotlabel);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),nres-1,nres-1,nres);
/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */	/* fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres); */
/* k1-1 error should be nres-1*/	/* k1-1 error should be nres-1*/
Line 7034 void printinggnuplot(char fileresu[], ch	Line 7129 void printinggnuplot(char fileresu[], ch

fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_"));	fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_"));
if(cptcoveff ==0){	if(cptcoveff ==0){
fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3", 2+(cpt-1), cpt );	fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3", 2+3*(cpt-1), cpt );
}else{	}else{
kl=0;	kl=0;
for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */	for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */
Line 7092 void printinggnuplot(char fileresu[], ch	Line 7187 void printinggnuplot(char fileresu[], ch
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\"Backward (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);	fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 6 dt 3,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
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 lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);	fprintf(ficgp,"\" t\"95%% CI\" w l lt 4,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
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\"\" w l lt 1");	fprintf(ficgp,"\" t\"\" w l lt 4");
} /* end if backprojcast */	} /* end if backprojcast */
} /* end if backcast */	} /* end if backcast */
fprintf(ficgp,"\nset out ;unset label;\n");	/* fprintf(ficgp,"\nset out ;unset label;\n"); */
	fprintf(ficgp,"\nset out ;unset title;\n");
} /* nres */	} /* nres */
} /* k1 */	} /* k1 */
} /* cpt */	} /* cpt */
Line 7503 set ter svg size 640, 480\nunset log y\n	Line 7599 set ter svg size 640, 480\nunset log y\n
/# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 /	/# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 /
fprintf(ficgp," u %d:(", ioffset);	fprintf(ficgp," u %d:(", ioffset);
if(i==nlstate+1){	if(i==nlstate+1){
fprintf(ficgp," $%d/(1.-$%d)):5 t 'pw.%d' with line lc variable ", \	fprintf(ficgp," $%d/(1.-$%d)):1 t 'pw.%d' with line lc variable ", \
ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,cpt );	ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,cpt );
fprintf(ficgp,",\\\n '' ");	fprintf(ficgp,",\\\n '' ");
fprintf(ficgp," u %d:(",ioffset);	fprintf(ficgp," u %d:(",ioffset);
fprintf(ficgp," (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", \	fprintf(ficgp," (($1-$2) == %d ) ? $%d/(1.-$%d) : 1/0):1 with labels center not ", \
offyear, \	offyear, \
ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate );	ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate );
}else	}else
fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \	fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ", \
ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,i,cpt );	ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,i,cpt );
}else{ /* more than 2 covariates */	}else{ /* more than 2 covariates */
if(cptcoveff ==1){	ioffset=2cptcoveff+2; / Age is in 4 or 6 or etc.*/
ioffset=4; /* Age is in 4 */	/# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3/
}else{	/# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 /
ioffset=6; /* Age is in 6 */	iyearc=ioffset-1;
/# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3/	iagec=ioffset;
/# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 /
}
fprintf(ficgp," u %d:(",ioffset);	fprintf(ficgp," u %d:(",ioffset);
kl=0;	kl=0;
strcpy(gplotcondition,"(");	strcpy(gplotcondition,"(");
Line 7542 set ter svg size 640, 480\nunset log y\n	Line 7636 set ter svg size 640, 480\nunset log y\n
/6+1+(i-1)+(nlstate+1)nlstate; 6+1+(1-1) +(2+1)2=13 /	/6+1+(i-1)+(nlstate+1)nlstate; 6+1+(1-1) +(2+1)2=13 /
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/	/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
if(i==nlstate+1){	if(i==nlstate+1){
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):5 t 'p.%d' with line lc variable", gplotcondition, \	fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \
ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,cpt );	ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,iyearc, cpt );
fprintf(ficgp,",\\\n '' ");	fprintf(ficgp,",\\\n '' ");
fprintf(ficgp," u %d:(",ioffset);	fprintf(ficgp," u %d:(",iagec);
fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \	fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \
offyear, \	iyearc, iagec, offyear, \
ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate );	ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate, iyearc );
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/	/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
}else{	}else{
fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \	fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
Line 7618 set ter svg size 640, 480\nunset log y\n	Line 7712 set ter svg size 640, 480\nunset log y\n
/# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 /	/# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 /
fprintf(ficgp," u %d:(", ioffset);	fprintf(ficgp," u %d:(", ioffset);
if(i==nlstate+1){	if(i==nlstate+1){
fprintf(ficgp," $%d/(1.-$%d)):5 t 'bw%d' with line lc variable ", \	fprintf(ficgp," $%d/(1.-$%d)):1 t 'bw%d' with line lc variable ", \
ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,cpt );	ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,cpt );
fprintf(ficgp,",\\\n '' ");	fprintf(ficgp,",\\\n '' ");
fprintf(ficgp," u %d:(",ioffset);	fprintf(ficgp," u %d:(",ioffset);
fprintf(ficgp," (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", \	fprintf(ficgp," (($1-$2) == %d ) ? $%d : 1/0):1 with labels center not ", \
offbyear, \	offbyear, \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );	ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );
}else	}else
fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ", \	fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ", \
ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,cpt,i );	ioffset+(cpt-1)(nlstate+1)+1+(i-1), ioffset+1+(i-1)+(nlstate+1)nlstate,cpt,i );
}else{ /* more than 2 covariates */	}else{ /* more than 2 covariates */
if(cptcoveff ==1){	ioffset=2cptcoveff+2; / Age is in 4 or 6 or etc.*/
ioffset=4; /* Age is in 4 */	/# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3/
}else{	/# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 /
ioffset=6; /* Age is in 6 */	iyearc=ioffset-1;
/# V1 = 1 V2 = 0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3/	iagec=ioffset;
/# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 /
}
fprintf(ficgp," u %d:(",ioffset);	fprintf(ficgp," u %d:(",ioffset);
kl=0;	kl=0;
strcpy(gplotcondition,"(");	strcpy(gplotcondition,"(");
Line 7657 set ter svg size 640, 480\nunset log y\n	Line 7749 set ter svg size 640, 480\nunset log y\n
/6+1+(i-1)+(nlstate+1)nlstate; 6+1+(1-1) +(2+1)2=13 /	/6+1+(i-1)+(nlstate+1)nlstate; 6+1+(1-1) +(2+1)2=13 /
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/	/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
if(i==nlstate+1){	if(i==nlstate+1){
fprintf(ficgp,"%s ? $%d : 1/0):5 t 'bw%d' with line lc variable", gplotcondition, \	fprintf(ficgp,"%s ? $%d : 1/0):%d t 'bw%d' with line lc variable", gplotcondition, \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1),cpt );	ioffset+(cpt-1)*(nlstate+1)+1+(i-1),iyearc,cpt );
fprintf(ficgp,",\\\n '' ");	fprintf(ficgp,",\\\n '' ");
fprintf(ficgp," u %d:(",ioffset);	fprintf(ficgp," u %d:(",iagec);
/* fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */	/* fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */
fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", gplotcondition, \	fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d : 1/0):%d with labels center not ", gplotcondition, \
offbyear, \	iyearc,iagec,offbyear, \
ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );	ioffset+(cpt-1)*(nlstate+1)+1+(i-1), iyearc );
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/	/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
}else{	}else{
/* fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ */	/* fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ */
Line 7723 set ter svg size 640, 480\nunset log y\n	Line 7815 set ter svg size 640, 480\nunset log y\n
continue;	continue;
fprintf(ficgp,"\n\n# Combination of dummy k1=%d which is ",k1);	fprintf(ficgp,"\n\n# Combination of dummy k1=%d which is ",k1);
strcpy(gplotlabel,"(");	strcpy(gplotlabel,"(");
sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);	/sprintf(gplotlabel+strlen(gplotlabel)," Dummy combination %d ",k1);/
for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */	for (k=1; k<=cptcoveff; k++){ /* For each correspondig covariate value */
lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */	lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate value corresponding to k1 combination and kth covariate */
/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */	/* decodtabm(1,1,4) = 1 because h=1 k= (1) 1 1 1 */
Line 7740 set ter svg size 640, 480\nunset log y\n	Line 7832 set ter svg size 640, 480\nunset log y\n
strcpy(gplotlabel+strlen(gplotlabel),")");	strcpy(gplotlabel+strlen(gplotlabel),")");
fprintf(ficgp,"\n#\n");	fprintf(ficgp,"\n#\n");
fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),k1,ng,nres);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),k1,ng,nres);
fprintf(ficgp,"\nset label \"%s\" at graph 0.98,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel);	fprintf(ficgp,"\nset key outside ");
	/* fprintf(ficgp,"\nset label \"%s\" at graph 1.2,0.5 center rotate font \"Helvetica,12\"\n",gplotlabel); */
	fprintf(ficgp,"\nset title \"%s\" font \"Helvetica,12\"\n",gplotlabel);
fprintf(ficgp,"\nset ter svg size 640, 480 ");	fprintf(ficgp,"\nset ter svg size 640, 480 ");
if (ng==1){	if (ng==1){
fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12x) could be nice /	fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12x) could be nice /
Line 7860 set ter svg size 640, 480\nunset log y\n	Line 7954 set ter svg size 640, 480\nunset log y\n
}	}
fprintf(ficgp,")");	fprintf(ficgp,")");
if(ng ==2)	if(ng ==2)
fprintf(ficgp," t \"p%d%d\" ", k2,k);	fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
else /* ng= 3 */	else /* ng= 3 */
fprintf(ficgp," t \"i%d%d\" ", k2,k);	fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
}else{ /* end ng <> 1 */	}else{ /* end ng <> 1 */
if( k !=k2) /* logit p11 is hard to draw */	if( k !=k2) /* logit p11 is hard to draw */
fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);	fprintf(ficgp," w l lw 2 lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
}	}
if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)	if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
fprintf(ficgp,",");	fprintf(ficgp,",");
Line 7874 set ter svg size 640, 480\nunset log y\n	Line 7968 set ter svg size 640, 480\nunset log y\n
i=i+ncovmodel;	i=i+ncovmodel;
} /* end k */	} /* end k */
} /* end k2 */	} /* end k2 */
fprintf(ficgp,"\n set out; unset label;\n");	/* fprintf(ficgp,"\n set out; unset label;set key default;\n"); */
	fprintf(ficgp,"\n set out; unset title;set key default;\n");
} /* end k1 */	} /* end k1 */
} /* end ng */	} /* end ng */
/* avoid: */	/* avoid: */
Line 7898 set ter svg size 640, 480\nunset log y\n	Line 7993 set ter svg size 640, 480\nunset log y\n
double agemingoodr, agemaxgoodr;	double agemingoodr, agemaxgoodr;


/* modcovmax=2cptcoveff;/\ Max number of modalities. We suppose */	/* modcovmax=2cptcoveff; Max number of modalities. We suppose /
/* a covariate has 2 modalities, should be equal to ncovcombmax \/ /	/* a covariate has 2 modalities, should be equal to ncovcombmax */

sumnewp = vector(1,ncovcombmax);	sumnewp = vector(1,ncovcombmax);
sumnewm = vector(1,ncovcombmax);	sumnewm = vector(1,ncovcombmax);
Line 8112 set ter svg size 640, 480\nunset log y\n	Line 8207 set ter svg size 640, 480\nunset log y\n


/************ Forecasting ****************/	/************ Forecasting ****************/
void 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){	void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
/* proj1, year, month, day of starting projection	/* proj1, year, month, day of starting projection
agemin, agemax range of age	agemin, agemax range of age
dateprev1 dateprev2 range of dates during which prevalence is computed	dateprev1 dateprev2 range of dates during which prevalence is computed
Line 8151 void prevforecast(char fileres[], double	Line 8246 void prevforecast(char fileres[], double
if(estepm < stepm){	if(estepm < stepm){
printf ("Problem %d lower than %d\n",estepm, stepm);	printf ("Problem %d lower than %d\n",estepm, stepm);
}	}
else hstepm=estepm;	else{
	hstepm=estepm;
	}
	if(estepm > stepm){ /* Yes every two year */
	stepsize=2;
	}

hstepm=hstepm/stepm;	hstepm=hstepm/stepm;
yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and	yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
Line 8196 void prevforecast(char fileres[], double	Line 8296 void prevforecast(char fileres[], double
for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {	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+yearp);	fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);
for (agec=fage; agec>=(ageminpar-1); agec--){	/* for (agec=fage; agec>=(ageminpar-1); agec--){ */
	for (agec=fage; agec>=(bage); agec--){
nhstepm=(int) rint((agelim-agec)*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);
Line 8217 void prevforecast(char fileres[], double	Line 8318 void prevforecast(char fileres[], double
for(j=1; j<=nlstate+ndeath;j++) {	for(j=1; j<=nlstate+ndeath;j++) {
ppij=0.;	ppij=0.;
for(i=1; i<=nlstate;i++) {	for(i=1; i<=nlstate;i++) {
/* if (mobilav>=1) */	if (mobilav>=1)
ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];	ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k];
/* else { / / even if mobilav==-1 we use mobaverage */	else { /* even if mobilav==-1 we use mobaverage, probs may not sums to 1 */
/* ppij=ppij+p3mat[i][j][h]probs[(int)(agec)][i][k]; /	ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k];
/* } */	}
fprintf(ficresf," %.3f", p3mat[i][j][h]);	fprintf(ficresf," %.3f", p3mat[i][j][h]);
} /* end i */	} /* end i */
fprintf(ficresf," %.3f", ppij);	fprintf(ficresf," %.3f", ppij);
Line 8239 void prevforecast(char fileres[], double	Line 8340 void prevforecast(char fileres[], double

}	}

/* /\************ Back Forecasting ***************\/ /	/************ Back Forecasting ****************/
void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){	void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
/* back1, year, month, day of starting backection	/* back1, year, month, day of starting backection
agemin, agemax range of age	agemin, agemax range of age
dateprev1 dateprev2 range of dates during which prevalence is computed	dateprev1 dateprev2 range of dates during which prevalence is computed
anback2 year of en of backection (same day and month as back1).	anback2 year of end of backprojection (same day and month as back1).
	prevacurrent and prev are prevalences.
*/	*/
int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;	int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
double agec; /* generic age */	double agec; /* generic age */
Line 8283 void prevbackforecast(char fileres[], do	Line 8385 void prevbackforecast(char fileres[], do
if(estepm < stepm){	if(estepm < stepm){
printf ("Problem %d lower than %d\n",estepm, stepm);	printf ("Problem %d lower than %d\n",estepm, stepm);
}	}
else hstepm=estepm;	else{
	hstepm=estepm;
	}
	if(estepm >= stepm){ /* Yes every two year */
	stepsize=2;
	}

hstepm=hstepm/stepm;	hstepm=hstepm/stepm;
yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and	yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp and
Line 8304 void prevbackforecast(char fileres[], do	Line 8411 void prevbackforecast(char fileres[], do

fprintf(ficresfb,"#**** Routine prevbackforecast \n");	fprintf(ficresfb,"#**** Routine prevbackforecast \n");

/* if (h==(int)(YEARMyearp)){ /
/* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */
/* for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */
/* k=k+1; */
for(nres=1; nres <= nresult; nres++) /* For each resultline */	for(nres=1; nres <= nresult; nres++) /* For each resultline */
for(k=1; k<=i1;k++){	for(k=1; k<=i1;k++){
if(i1 != 1 && TKresult[nres]!= k)	if(i1 != 1 && TKresult[nres]!= k)
Line 8333 void prevbackforecast(char fileres[], do	Line 8436 void prevbackforecast(char fileres[], do
/* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */	/* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) { */
fprintf(ficresfb,"\n");	fprintf(ficresfb,"\n");
fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);	fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);	/* printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */
/* for (agec=fage; agec>=(ageminpar-1); agec--){ */	/* for (agec=bage; agec<=agemax-1; agec++){ /\* testing \/ /
/* nhstepm=(int) rint((agelim-agec)YEARM/stepm); /	for (agec=bage; agec<=fage; agec++){ /* testing */
for (agec=bage; agec<=agemax-1; agec++){ /* testing */
/* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/	/* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/
nhstepm=(int) rint((agec-agelim)*YEARM/stepm);	nhstepm=(int) (agec-agelim) YEARM/stepm;/ nhstepm=(int) rint((agec-agelim)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;
/* computes hbxij at age agec over 1 to nhstepm */	/* computes hbxij at age agec over 1 to nhstepm */
	/* printf("####prevbackforecast debug agec=%.2f nhstepm=%d\n",agec, nhstepm);fflush(stdout); */
hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);	hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
/* hpxij(p3mat,nhstepm,agec,hstepm,p, nlstate,stepm,oldm,savm, k,nres); */	/* hpxij(p3mat,nhstepm,agec,hstepm,p, nlstate,stepm,oldm,savm, k,nres); */
/* Then we print p3mat for h corresponding to the right agec+hstepms=yearp /	/* Then we print p3mat for h corresponding to the right agec+hstepms=yearp /
Line 8360 void prevbackforecast(char fileres[], do	Line 8463 void prevbackforecast(char fileres[], do
ppij=0.;ppi=0.;	ppij=0.;ppi=0.;
for(j=1; j<=nlstate;j++) {	for(j=1; j<=nlstate;j++) {
/* if (mobilav==1) */	/* if (mobilav==1) */
ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k];	ppij=ppij+p3mat[i][j][h]*prevacurrent[(int)agec][j][k];
ppi=ppi+mobaverage[(int)agec][j][k];	ppi=ppi+prevacurrent[(int)agec][j][k];
	/* ppij=ppij+p3mat[i][j][h]mobaverage[(int)agec][j][k]; /
	/* ppi=ppi+mobaverage[(int)agec][j][k]; */
/* else { */	/* else { */
/* ppij=ppij+p3mat[i][j][h]probs[(int)(agec)][i][k]; /	/* ppij=ppij+p3mat[i][j][h]probs[(int)(agec)][i][k]; /
/* } */	/* } */
Line 8386 void prevbackforecast(char fileres[], do	Line 8491 void prevbackforecast(char fileres[], do

}	}

	/* Variance of prevalence limit: varprlim */
	void varprlim(char fileresu[], int nresult, double *prevacurrent, int mobilavproj, double bage, double fage, double prlim, int ncvyearp, double ftolpl, double p[], double matcov, double delti, int stepm, int cptcoveff){
	/------- Variance of period (stable) prevalence------/

	char fileresvpl[FILENAMELENGTH];
	FILE *ficresvpl;
	double oldm, savm;
	double *varpl; / Variances of prevalence limits by age */
	int i1, k, nres, j ;

	strcpy(fileresvpl,"VPL_");
	strcat(fileresvpl,fileresu);
	if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
	printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl);
	exit(0);
	}
	printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
	fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);

	/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
	for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/

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

	for(nres=1; nres <= nresult; nres++) /* For each resultline */
	for(k=1; k<=i1;k++){
	if(i1 != 1 && TKresult[nres]!= k)
	continue;
	fprintf(ficresvpl,"\n#****** ");
	printf("\n#****** ");
	fprintf(ficlog,"\n#****** ");
	for(j=1;j<=cptcoveff;j++) {
	fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
	}
	for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
	fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
	}
	fprintf(ficresvpl,"******\n");
	printf("******\n");
	fprintf(ficlog,"******\n");

	varpl=matrix(1,nlstate,(int) bage, (int) fage);
	oldm=oldms;savm=savms;
	varprevlim(fileresvpl, ficresvpl, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyearp, k, strstart, nres);
	free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
	/}/
	}

	fclose(ficresvpl);
	printf("done variance-covariance of period prevalence\n");fflush(stdout);
	fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);

	}
	/* Variance of back prevalence: varbprlim */
	void varbprlim(char fileresu[], int nresult, double *prevacurrent, int mobilavproj, double bage, double fage, double bprlim, int ncvyearp, double ftolpl, double p[], double matcov, double delti, int stepm, int cptcoveff){
	/------- Variance of back (stable) prevalence------/

	char fileresvbl[FILENAMELENGTH];
	FILE *ficresvbl;

	double oldm, savm;
	double *varbpl; / Variances of back prevalence limits by age */
	int i1, k, nres, j ;

	strcpy(fileresvbl,"VBL_");
	strcat(fileresvbl,fileresu);
	if((ficresvbl=fopen(fileresvbl,"w"))==NULL) {
	printf("Problem with variance of back (stable) prevalence resultfile: %s\n", fileresvbl);
	exit(0);
	}
	printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout);
	fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog);


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

	for(nres=1; nres <= nresult; nres++) /* For each resultline */
	for(k=1; k<=i1;k++){
	if(i1 != 1 && TKresult[nres]!= k)
	continue;
	fprintf(ficresvbl,"\n#****** ");
	printf("\n#****** ");
	fprintf(ficlog,"\n#****** ");
	for(j=1;j<=cptcoveff;j++) {
	fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
	fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
	}
	for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
	printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
	fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
	fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
	}
	fprintf(ficresvbl,"******\n");
	printf("******\n");
	fprintf(ficlog,"******\n");

	varbpl=matrix(1,nlstate,(int) bage, (int) fage);
	oldm=oldms;savm=savms;

	varbrevlim(fileresvbl, ficresvbl, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, ncvyearp, k, strstart, nres);
	free_matrix(varbpl,1,nlstate,(int) bage, (int)fage);
	/}/
	}

	fclose(ficresvbl);
	printf("done variance-covariance of back prevalence\n");fflush(stdout);
	fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog);

	} /* End of varbprlim */

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

Line 10442 int main(int argc, char *argv[])	Line 10664 int main(int argc, char *argv[])
int vpopbased=0;	int vpopbased=0;
int nres=0;	int nres=0;
int endishere=0;	int endishere=0;
	int noffset=0;
	int ncurrv=0; /* Temporary variable */

char ca[32], cb[32];	char ca[32], cb[32];
/* FILE fichtm; //* Html File */	/* FILE fichtm; //* Html File */
/* FILE ficgp;/ /Gnuplot File /	/* FILE ficgp;/ /Gnuplot File /
Line 10494 int main(int argc, char *argv[])	Line 10718 int main(int argc, char *argv[])
double delti; / Scale */	double delti; / Scale */
double *eij, *vareij;	double *eij, *vareij;
double *varpl; / Variances of prevalence limits by age */	double *varpl; / Variances of prevalence limits by age */
double *varbpl; / Variances of back prevalence limits by age */
double *epj, vepp;	double *epj, vepp;

double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;	double dateprev1, dateprev2;
double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000;	double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0;
	double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0;

double **ximort;	double **ximort;
char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";	char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
Line 10576 int main(int argc, char *argv[])	Line 10801 int main(int argc, char *argv[])
if(pathr[0] == '\0') break; /* Dirty */	if(pathr[0] == '\0') break; /* Dirty */
}	}
}	}
	else if (argc<=2){
	strcpy(pathtot,argv[1]);
	}
else{	else{
strcpy(pathtot,argv[1]);	strcpy(pathtot,argv[1]);
	strcpy(z,argv[2]);
	printf("\nargv[2]=%s z=%c\n",argv[2],z[0]);
}	}
/if(getcwd(pathcd, MAXLINE)!= NULL)printf ("Error pathcd\n");/	/if(getcwd(pathcd, MAXLINE)!= NULL)printf ("Error pathcd\n");/
/*cygwin_split_path(pathtot,path,optionfile);	/*cygwin_split_path(pathtot,path,optionfile);
Line 10655 int main(int argc, char *argv[])	Line 10885 int main(int argc, char *argv[])
exit(70);	exit(70);
}	}



strcpy(filereso,"o");	strcpy(filereso,"o");
strcat(filereso,fileresu);	strcat(filereso,fileresu);
if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */	if((ficparo=fopen(filereso,"w"))==NULL) { /* opened on subdirectory */
Line 10665 int main(int argc, char *argv[])	Line 10893 int main(int argc, char *argv[])
fflush(ficlog);	fflush(ficlog);
goto end;	goto end;
}	}
	/-------- Rewriting parameter file ----------/
	strcpy(rfileres,"r"); /* "Rparameterfile */
	strcat(rfileres,optionfilefiname); /* Parameter file first name */
	strcat(rfileres,"."); /* */
	strcat(rfileres,optionfilext); /* Other files have txt extension */
	if((ficres =fopen(rfileres,"w"))==NULL) {
	printf("Problem writing new parameter file: %s\n", rfileres);goto end;
	fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;
	fflush(ficlog);
	goto end;
	}
	fprintf(ficres,"#IMaCh %s\n",version);


/* Reads comments: lines beginning with '#' */	/* Reads comments: lines beginning with '#' */
numlinepar=0;	numlinepar=0;
	/* Is it a BOM UTF-8 Windows file? */
/* First parameter line */	/* First parameter line */
while(fgets(line, MAXLINE, ficpar)) {	while(fgets(line, MAXLINE, ficpar)) {
	noffset=0;
	if( line[0] == (char)0xEF && line[1] == (char)0xBB) /* EF BB BF */
	{
	noffset=noffset+3;
	printf("# File is an UTF8 Bom.\n"); // 0xBF
	}
	else if( line[0] == (char)0xFE && line[1] == (char)0xFF)
	{
	noffset=noffset+2;
	printf("# File is an UTF16BE BOM file\n");
	}
	else if( line[0] == 0 && line[1] == 0)
	{
	if( line[2] == (char)0xFE && line[3] == (char)0xFF){
	noffset=noffset+4;
	printf("# File is an UTF16BE BOM file\n");
	}
	} else{
	;/printf(" Not a BOM file\n");/
	}

/* If line starts with a # it is a comment */	/* If line starts with a # it is a comment */
if (line[0] == '#') {	if (line[noffset] == '#') {
numlinepar++;	numlinepar++;
fputs(line,stdout);	fputs(line,stdout);
fputs(line,ficparo);	fputs(line,ficparo);
	fputs(line,ficres);
fputs(line,ficlog);	fputs(line,ficlog);
continue;	continue;
}else	}else
Line 10696 int main(int argc, char *argv[])	Line 10959 int main(int argc, char *argv[])
numlinepar++;	numlinepar++;
fputs(line,stdout);	fputs(line,stdout);
fputs(line,ficparo);	fputs(line,ficparo);
	fputs(line,ficres);
fputs(line,ficlog);	fputs(line,ficlog);
continue;	continue;
}else	}else
Line 10718 int main(int argc, char *argv[])	Line 10982 int main(int argc, char *argv[])
numlinepar++;	numlinepar++;
fputs(line,stdout);	fputs(line,stdout);
fputs(line,ficparo);	fputs(line,ficparo);
	fputs(line,ficres);
fputs(line,ficlog);	fputs(line,ficlog);
continue;	continue;
}else	}else
break;	break;
}	}
if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){	if((num_filled=sscanf(line,"model=1+age%[^.\n]", model)) !=EOF){
if (num_filled == 0){	if (num_filled != 1){
printf("ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line);	printf("ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line);
fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' WITHOUT space:'%s'\n",num_filled, line);	fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age' %s\n",num_filled, line);
model[0]='\0';
goto end;
} else if (num_filled != 1){
printf("ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
fprintf(ficlog,"ERROR %d: Model should be at minimum 'model=1+age.' %s\n",num_filled, line);
model[0]='\0';	model[0]='\0';
goto end;	goto end;
}	}
Line 10753 int main(int argc, char *argv[])	Line 11013 int main(int argc, char *argv[])
fflush(ficlog);	fflush(ficlog);
/* if(model[0]=='#'\|\| model[0]== '\0'){ */	/* if(model[0]=='#'\|\| model[0]== '\0'){ */
if(model[0]=='#'){	if(model[0]=='#'){
printf("Error in 'model' line: model should start with 'model=1+age+' and end with '.' \n \	printf("Error in 'model' line: model should start with 'model=1+age+' and end without space \n \
'model=1+age+.' or 'model=1+age+V1.' or 'model=1+age+ageage+V1+V1age.' or \n \	'model=1+age+' or 'model=1+age+V1.' or 'model=1+age+ageage+V1+V1age' or \n \
'model=1+age+V1+V2.' or 'model=1+age+V1+V2+V1*V2.' etc. \n"); \	'model=1+age+V1+V2' or 'model=1+age+V1+V2+V1*V2' etc. \n"); \
if(mle != -1){	if(mle != -1){
printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter file.\n");	printf("Fix the model line and run imach with mle=-1 to get a correct template of the parameter vectors and subdiagonal covariance matrix.\n");
exit(1);	exit(1);
}	}
}	}
Line 10979 Please run with mle=-1 to get a correct	Line 11239 Please run with mle=-1 to get a correct

fflush(ficlog);	fflush(ficlog);

/-------- Rewriting parameter file ----------/
strcpy(rfileres,"r"); /* "Rparameterfile */
strcat(rfileres,optionfilefiname); /* Parameter file first name*/
strcat(rfileres,"."); /* */
strcat(rfileres,optionfilext); /* Other files have txt extension */
if((ficres =fopen(rfileres,"w"))==NULL) {
printf("Problem writing new parameter file: %s\n", rfileres);goto end;
fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;
}
fprintf(ficres,"#%s\n",version);
} /* End of mle != -3 */	} /* End of mle != -3 */

/* Main data	/* Main data
Line 11326 Title=%s <br>Datafile=%s Firstpass=%d La	Line 11576 Title=%s <br>Datafile=%s Firstpass=%d La
firstpass, lastpass, stepm, weightopt, model);	firstpass, lastpass, stepm, weightopt, model);

fprintf(fichtm,"\n");	fprintf(fichtm,"\n");
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\	fprintf(fichtm,"<h4>Parameter line 2</h4><ul><li>Tolerance for the convergence of the likelihood: ftol=%f \n<li>Interval for the elementary matrix (in month): stepm=%d",\
	ftol, stepm);
	fprintf(fichtm,"\n<li>Number of fixed dummy covariates: ncovcol=%d ", ncovcol);
	ncurrv=1;
	for(i=ncurrv; i <=ncovcol; i++) fprintf(fichtm,"V%d ", i);
	fprintf(fichtm,"\n<li> Number of fixed quantitative variables: nqv=%d ", nqv);
	ncurrv=i;
	for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", i);
	fprintf(fichtm,"\n<li> Number of time varying (wave varying) covariates: ntv=%d ", ntv);
	ncurrv=i;
	for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", i);
	fprintf(fichtm,"\n<li>Number of quantitative time varying covariates: nqtv=%d ", nqtv);
	ncurrv=i;
	for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", i);
	fprintf(fichtm,"\n<li>Weights column \n<br>Number of alive states: nlstate=%d <br>Number of death states (not really implemented): ndeath=%d \n<li>Number of waves: maxwav=%d \n<li>Parameter for maximization (1), using parameter values (0), for design of parameters and variance-covariance matrix: mle=%d \n<li>Does the weight column be taken into account (1), or not (0): weight=%d</ul>\n", \
	nlstate, ndeath, maxwav, mle, weightopt);

	fprintf(fichtm,"<h4> Diagram of states <a href=\"%s_.svg\">%s_.svg</a></h4> \n\
	<img src=\"%s_.svg\">", subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"));


	fprintf(fichtm,"\n<h4>Some descriptive statistics </h4>\n<br>Total number of observations=%d <br>\n\
Youngest age at first (selected) pass %.2f, oldest age %.2f<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",\
imx,agemin,agemax,jmin,jmax,jmean);	imx,agemin,agemax,jmin,jmax,jmean);
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */	pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */	oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */	newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
Line 11889 Please run with mle=-1 to get a correct	Line 12160 Please run with mle=-1 to get a correct
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(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);	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.*/	/* day and month of proj2 are not used but only year anproj2.*/
	dateproj1=anproj1+(mproj1-1)/12.+(jproj1-1)/365.;
	dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.;

}	}
break;	break;
case 12:	case 12:
Line 11904 Please run with mle=-1 to get a correct	Line 12178 Please run with mle=-1 to get a correct
fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);	fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);	fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
/* day and month of proj2 are not used but only year anproj2.*/	/* day and month of proj2 are not used but only year anproj2.*/
	dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.;
	dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.;
}	}
break;	break;
case 13:	case 13:
Line 11954 Please run with mle=-1 to get a correct	Line 12230 Please run with mle=-1 to get a correct
This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\	This is probably because your parameter file doesn't \n contain the exact number of lines (or columns) corresponding to your model line.\n\
Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);	Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
}else{	}else{
printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1);	/* printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1); */
	printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, backcast, pathc,p, (int)anproj1-bage, (int)anback1-fage);
}	}
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \	printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \	model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);	jprev1,mprev1,anprev1,dateprev1, dateproj1, dateback1,jprev2,mprev2,anprev2,dateprev2,dateproj2, dateback2);

/------------ free_vector -------------/	/------------ free_vector -------------/
/* chdir(path); */	/* chdir(path); */
Line 11994 Please run with mle=-1 to get a correct	Line 12271 Please run with mle=-1 to get a correct
k=1;	k=1;
varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);	varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);

/* Prevalence for each covariates in probs[age][status][cov] */	/* Prevalence for each covariate combination in probs[age][status][cov] */
probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);	probs= ma3x(AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
for(i=1;i<=AGESUP;i++)	for(i=AGEINF;i<=AGESUP;i++)
for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */	for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */
for(k=1;k<=ncovcombmax;k++)	for(k=1;k<=ncovcombmax;k++)
probs[i][j][k]=0.;	probs[i][j][k]=0.;
prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);	prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode,
	ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
if (mobilav!=0 \|\|mobilavproj !=0 ) {	if (mobilav!=0 \|\|mobilavproj !=0 ) {
mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);	mobaverages= ma3x(AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
for(i=1;i<=AGESUP;i++)	for(i=AGEINF;i<=AGESUP;i++)
for(j=1;j<=nlstate+ndeath;j++)	for(j=1;j<=nlstate+ndeath;j++)
for(k=1;k<=ncovcombmax;k++)	for(k=1;k<=ncovcombmax;k++)
mobaverages[i][j][k]=0.;	mobaverages[i][j][k]=0.;
Line 12015 Please run with mle=-1 to get a correct	Line 12293 Please run with mle=-1 to get a correct
fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);	fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
printf(" Error in movingaverage mobilav=%d\n",mobilav);	printf(" Error in movingaverage mobilav=%d\n",mobilav);
}	}
}	} else if (mobilavproj !=0) {
/* else if(mobilavproj==-1){ /\* Forcing raw observed prevalences \/ /
/* for(i=1;i<=AGESUP;i++) */
/* for(j=1;j<=nlstate;j++) */
/* for(k=1;k<=ncovcombmax;k++) */
/* mobaverages[i][j][k]=probs[i][j][k]; */
/* /\* /\\* Prevalence for each covariates in probs[age][status][cov] \\/ \/ */
/* /\* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); \/ /
/* } */
else if (mobilavproj !=0) {
printf("Movingaveraging projected observed prevalence\n");	printf("Movingaveraging projected observed prevalence\n");
fprintf(ficlog,"Movingaveraging projected observed prevalence\n");	fprintf(ficlog,"Movingaveraging projected observed prevalence\n");
if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){	if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);	fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);	printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
}	}
	}else{
	printf("Internal error moving average\n");
	fflush(stdout);
	exit(1);
}	}
}/* end if moving average */	}/* end if moving average */

/---------- Forecasting ------------------/	/---------- Forecasting ------------------/
/if((stepm == 1) && (strcmp(model,".")==0)){/
if(prevfcast==1){	if(prevfcast==1){
/* if(stepm ==1){*/	/* if(stepm ==1){*/
prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);	prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
}	}

	/* Backcasting */
if(backcast==1){	if(backcast==1){
ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);	ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);	ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);
Line 12048 Please run with mle=-1 to get a correct	Line 12322 Please run with mle=-1 to get a correct
/--------------- Back Prevalence limit (period or stable prevalence) --------------/	/--------------- Back Prevalence limit (period or stable prevalence) --------------/

bprlim=matrix(1,nlstate,1,nlstate);	bprlim=matrix(1,nlstate,1,nlstate);

back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);	back_prevalence_limit(p, bprlim, ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);
fclose(ficresplb);	fclose(ficresplb);

hBijx(p, bage, fage, mobaverage);	hBijx(p, bage, fage, mobaverage);
fclose(ficrespijb);	fclose(ficrespijb);
/* free_matrix(bprlim,1,nlstate,1,nlstate); /\here or after loop ? \/ */

prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,	prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2,
bage, fage, firstpass, lastpass, anback2, p, cptcoveff);	mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
	varbprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, bprlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);

/------- Variance of back (stable) prevalence------/

strcpy(fileresvbl,"VBL_");	free_matrix(bprlim,1,nlstate,1,nlstate); /here or after loop ? /
strcat(fileresvbl,fileresu);
if((ficresvbl=fopen(fileresvbl,"w"))==NULL) {
printf("Problem with variance of back (stable) prevalence resultfile: %s\n", fileresvbl);
exit(0);
}
printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout);
fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog);

/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/

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

for(nres=1; nres <= nresult; nres++) /* For each resultline */
for(k=1; k<=i1;k++){
if(i1 != 1 && TKresult[nres]!= k)
continue;
fprintf(ficresvbl,"\n#****** ");
printf("\n#****** ");
fprintf(ficlog,"\n#****** ");
for(j=1;j<=cptcoveff;j++) {
fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
}
fprintf(ficresvbl,"******\n");
printf("******\n");
fprintf(ficlog,"******\n");

varbpl=matrix(1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;

varbrevlim(fileres, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, &ncvyear, k, strstart, nres);
free_matrix(varbpl,1,nlstate,(int) bage, (int)fage);
/}/
}

fclose(ficresvbl);
printf("done variance-covariance of back prevalence\n");fflush(stdout);
fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog);

free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);	free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);	free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);	free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);
}	} /* end Backcasting */
free_matrix(bprlim,1,nlstate,1,nlstate); /here or after loop ? /


/* ------ Other prevalence ratios------------ */	/* ------ Other prevalence ratios------------ */
Line 12165 Please run with mle=-1 to get a correct	Line 12391 Please run with mle=-1 to get a correct
fclose(ficreseij);	fclose(ficreseij);
printf("done evsij\n");fflush(stdout);	printf("done evsij\n");fflush(stdout);
fprintf(ficlog,"done evsij\n");fflush(ficlog);	fprintf(ficlog,"done evsij\n");fflush(ficlog);


/---------- State-specific expectancies and variances ------------/	/---------- State-specific expectancies and variances ------------/


strcpy(filerest,"T_");	strcpy(filerest,"T_");
strcat(filerest,fileresu);	strcat(filerest,fileresu);
if((ficrest=fopen(filerest,"w"))==NULL) {	if((ficrest=fopen(filerest,"w"))==NULL) {
Line 12177 Please run with mle=-1 to get a correct	Line 12403 Please run with mle=-1 to get a correct
}	}
printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);	printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);
fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);	fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);


strcpy(fileresstde,"STDE_");	strcpy(fileresstde,"STDE_");
strcat(fileresstde,fileresu);	strcat(fileresstde,fileresu);
if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {	if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {
Line 12206 Please run with mle=-1 to get a correct	Line 12430 Please run with mle=-1 to get a correct
printf(" Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout);	printf(" Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout);
fprintf(ficlog," Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog);	fprintf(ficlog," Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog);

/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/

i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */	i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
if (cptcovn < 1){i1=1;}	if (cptcovn < 1){i1=1;}

Line 12270 Please run with mle=-1 to get a correct	Line 12491 Please run with mle=-1 to get a correct
vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);	vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
pstamp(ficrest);	pstamp(ficrest);

	epj=vector(1,nlstate+1);
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/	for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
oldm=oldms;savm=savms; /* ZZ Segmentation fault */	oldm=oldms;savm=savms; /* ZZ Segmentation fault */
cptcod= 0; /* To be deleted */	cptcod= 0; /* To be deleted */
Line 12286 Please run with mle=-1 to get a correct	Line 12507 Please run with mle=-1 to get a correct
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);	for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
fprintf(ficrest,"\n");	fprintf(ficrest,"\n");
/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */	/* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
epj=vector(1,nlstate+1);
printf("Computing age specific period (stable) prevalences in each health state \n");	printf("Computing age specific period (stable) prevalences in each health state \n");
fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");	fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
for(age=bage; age <=fage ;age++){	for(age=bage; age <=fage ;age++){
Line 12324 Please run with mle=-1 to get a correct	Line 12544 Please run with mle=-1 to get a correct
fprintf(ficrest,"\n");	fprintf(ficrest,"\n");
}	}
} /* End vpopbased */	} /* End vpopbased */
	free_vector(epj,1,nlstate+1);
free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);	free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);	free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
free_vector(epj,1,nlstate+1);
printf("done selection\n");fflush(stdout);	printf("done selection\n");fflush(stdout);
fprintf(ficlog,"done selection\n");fflush(ficlog);	fprintf(ficlog,"done selection\n");fflush(ficlog);

/}/
} /* End k selection */	} /* End k selection */

printf("done State-specific expectancies\n");fflush(stdout);	printf("done State-specific expectancies\n");fflush(stdout);
fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);	fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);

/------- Variance of period (stable) prevalence------/	/* variance-covariance of period prevalence*/
	varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
strcpy(fileresvpl,"VPL_");
strcat(fileresvpl,fileresu);
if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
printf("Problem with variance of period (stable) prevalence resultfile: %s\n", fileresvpl);
exit(0);
}
printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);

/*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/

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

for(nres=1; nres <= nresult; nres++) /* For each resultline */
for(k=1; k<=i1;k++){
if(i1 != 1 && TKresult[nres]!= k)
continue;
fprintf(ficresvpl,"\n#****** ");
printf("\n#****** ");
fprintf(ficlog,"\n#****** ");
for(j=1;j<=cptcoveff;j++) {
fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
}
for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
}
fprintf(ficresvpl,"******\n");
printf("******\n");
fprintf(ficlog,"******\n");

varpl=matrix(1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart, nres);
free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
/}/
}

fclose(ficresvpl);
printf("done variance-covariance of period prevalence\n");fflush(stdout);
fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);


free_vector(weight,1,n);	free_vector(weight,1,n);
Line 12402 Please run with mle=-1 to get a correct	Line 12575 Please run with mle=-1 to get a correct

/---------- End : free ----------------/	/---------- End : free ----------------/
if (mobilav!=0 \|\|mobilavproj !=0)	if (mobilav!=0 \|\|mobilavproj !=0)
free_ma3x(mobaverages,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */	free_ma3x(mobaverages,AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);	free_ma3x(probs,AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
free_matrix(prlim,1,nlstate,1,nlstate); /here or after loop ? /	free_matrix(prlim,1,nlstate,1,nlstate); /here or after loop ? /
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
} /* mle==-3 arrives here for freeing */	} /* mle==-3 arrives here for freeing */
Line 12420 Please run with mle=-1 to get a correct	Line 12593 Please run with mle=-1 to get a correct
/free_vector(delti,1,npar);/	/free_vector(delti,1,npar);/
free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);	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(paramstart,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);	free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);

free_ivector(ncodemax,1,NCOVMAX);	free_ivector(ncodemax,1,NCOVMAX);
Line 12496 Please run with mle=-1 to get a correct	Line 12670 Please run with mle=-1 to get a correct
fclose(ficlog);	fclose(ficlog);
/------ End -----------/	/------ End -----------/


	/* Executes gnuplot */

printf("Before Current directory %s!\n",pathcd);	printf("Before Current directory %s!\n",pathcd);
#ifdef WIN32	#ifdef WIN32
Line 12564 end:	Line 12740 end:
printf("\nType q for exiting: "); fflush(stdout);	printf("\nType q for exiting: "); fflush(stdout);
scanf("%s",z);	scanf("%s",z);
}	}
	exit(0);
}	}

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

Removed from v.1.268
changed lines
	Added in v.1.282