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.43 and 1.50

version 1.43, 2002/05/24 13:00:54	version 1.50, 2002/06/26 23:25:02
Line 75	Line 75
#define YEARM 12. /* Number of months per year */	#define YEARM 12. /* Number of months per year */
#define AGESUP 130	#define AGESUP 130
#define AGEBASE 40	#define AGEBASE 40
	#ifdef windows
	#define DIRSEPARATOR '\\'
	#define ODIRSEPARATOR '/'
	#else
	#define DIRSEPARATOR '/'
	#define ODIRSEPARATOR '\\'
	#endif

	char version[80]="Imach version 0.8i, June 2002, INED-EUROREVES ";
int erreur; /* Error number */	int erreur; /* Error number */
int nvar;	int nvar;
int cptcovn, cptcovage=0, cptcoveff=0,cptcov;	int cptcovn=0, cptcovage=0, cptcoveff=0,cptcov;
int npar=NPARMAX;	int npar=NPARMAX;
int nlstate=2; /* Number of live states */	int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */	int ndeath=1; /* Number of dead states */
Line 96 double jmean; /* Mean space between 2 wa	Line 103 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 ficgp,ficresprob,*ficpop;	FILE *ficlog;
	FILE ficgp,ficresprob,ficpop, ficresprobcov, *ficresprobcor;
	FILE *ficresprobmorprev;
	FILE fichtm; / Html File */
FILE *ficreseij;	FILE *ficreseij;
char filerese[FILENAMELENGTH];	char filerese[FILENAMELENGTH];
FILE *ficresvij;	FILE *ficresvij;
char fileresv[FILENAMELENGTH];	char fileresv[FILENAMELENGTH];
FILE *ficresvpl;	FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];	char fileresvpl[FILENAMELENGTH];
	char title[MAXLINE];
	char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
	char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];

	char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];
	char filelog[FILENAMELENGTH]; /* Log file */
	char filerest[FILENAMELENGTH];
	char fileregp[FILENAMELENGTH];
	char popfile[FILENAMELENGTH];

	char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH];

#define NR_END 1	#define NR_END 1
#define FREE_ARG char*	#define FREE_ARG char*
Line 161 static int split( char path, char dirc	Line 182 static int split( char path, char dirc

l1 = strlen( path ); /* length of path */	l1 = strlen( path ); /* length of path */
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );	if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
#ifdef windows	s= strrchr( path, DIRSEPARATOR ); /* find last / */
s = strrchr( path, '\\' ); /* find last / */
#else
s = strrchr( path, '/' ); /* find last / */
#endif
if ( s == NULL ) { /* no directory, so use current */	if ( s == NULL ) { /* no directory, so use current */
	/*if(strrchr(path, ODIRSEPARATOR )==NULL)
	printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
#if defined(__bsd__) /* get current working directory */	#if defined(__bsd__) /* get current working directory */
extern char *getwd( );	extern char *getwd( );

Line 232 int nbocc(char *s, char occ)	Line 251 int nbocc(char *s, char occ)

void cutv(char u,char v, char*t, char occ)	void cutv(char u,char v, char*t, char occ)
{	{
	/* cuts string t into u and v where u is ended by char occ excluding it
	and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2)
	gives u="abcedf" and v="ghi2j" */
int i,lg,j,p=0;	int i,lg,j,p=0;
i=0;	i=0;
for(j=0; j<=strlen(t)-1; j++) {	for(j=0; j<=strlen(t)-1; j++) {
Line 428 double brent(double ax, double bx, doubl	Line 450 double brent(double ax, double bx, doubl
tol2=2.0(tol1=tolfabs(x)+ZEPS);	tol2=2.0(tol1=tolfabs(x)+ZEPS);
/* if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret)))*/	/* if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret)))*/
printf(".");fflush(stdout);	printf(".");fflush(stdout);
	fprintf(ficlog,".");fflush(ficlog);
#ifdef DEBUG	#ifdef DEBUG
printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);	printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
	fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
/* if ((fabs(x-xm) <= (tol2-0.5(b-a)))\|\|(2.0fabs(fu-ftemp) <= ftol1.e-2(fabs(fu)+fabs(ftemp)))) { */	/* if ((fabs(x-xm) <= (tol2-0.5(b-a)))\|\|(2.0fabs(fu-ftemp) <= ftol1.e-2(fabs(fu)+fabs(ftemp)))) { */
#endif	#endif
if (fabs(x-xm) <= (tol2-0.5*(b-a))){	if (fabs(x-xm) <= (tol2-0.5*(b-a))){
Line 554 void linmin(double p[], double xi[], int	Line 578 void linmin(double p[], double xi[], int
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);	*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
#ifdef DEBUG	#ifdef DEBUG
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);	printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
	fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
#endif	#endif
for (j=1;j<=n;j++) {	for (j=1;j<=n;j++) {
xi[j] *= xmin;	xi[j] *= xmin;
Line 584 void powell(double p[], double **xi, int	Line 609 void powell(double p[], double **xi, int
ibig=0;	ibig=0;
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);
for (i=1;i<=n;i++)	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]);
printf("\n");	printf("\n");
	fprintf(ficlog,"\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);
#ifdef DEBUG	#ifdef DEBUG
printf("fret=%lf \n",*fret);	printf("fret=%lf \n",*fret);
	fprintf(ficlog,"fret=%lf \n",*fret);
#endif	#endif
printf("%d",i);fflush(stdout);	printf("%d",i);fflush(stdout);
	fprintf(ficlog,"%d",i);fflush(ficlog);
linmin(p,xit,n,fret,func);	linmin(p,xit,n,fret,func);
if (fabs(fptt-(*fret)) > del) {	if (fabs(fptt-(*fret)) > del) {
del=fabs(fptt-(*fret));	del=fabs(fptt-(*fret));
Line 601 void powell(double p[], double **xi, int	Line 631 void powell(double p[], double **xi, int
}	}
#ifdef DEBUG	#ifdef DEBUG
printf("%d %.12e",i,(*fret));	printf("%d %.12e",i,(*fret));
	fprintf(ficlog,"%d %.12e",i,(*fret));
for (j=1;j<=n;j++) {	for (j=1;j<=n;j++) {
xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);	xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
printf(" x(%d)=%.12e",j,xit[j]);	printf(" x(%d)=%.12e",j,xit[j]);
	fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
}	}
for(j=1;j<=n;j++)	for(j=1;j<=n;j++) {
printf(" p=%.12e",p[j]);	printf(" p=%.12e",p[j]);
	fprintf(ficlog," p=%.12e",p[j]);
	}
printf("\n");	printf("\n");
	fprintf(ficlog,"\n");
#endif	#endif
}	}
if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret))) {	if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret))) {
Line 616 void powell(double p[], double **xi, int	Line 651 void powell(double p[], double **xi, int
k[0]=1;	k[0]=1;
k[1]=-1;	k[1]=-1;
printf("Max: %.12e",(*func)(p));	printf("Max: %.12e",(*func)(p));
for (j=1;j<=n;j++)	fprintf(ficlog,"Max: %.12e",(*func)(p));
	for (j=1;j<=n;j++) {
printf(" %.12e",p[j]);	printf(" %.12e",p[j]);
	fprintf(ficlog," %.12e",p[j]);
	}
printf("\n");	printf("\n");
	fprintf(ficlog,"\n");
for(l=0;l<=1;l++) {	for(l=0;l<=1;l++) {
for (j=1;j<=n;j++) {	for (j=1;j<=n;j++) {
ptt[j]=p[j]+(p[j]-pt[j])*k[l];	ptt[j]=p[j]+(p[j]-pt[j])*k[l];
printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);	printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
	fprintf(ficlog,"l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
}	}
printf("func(ptt)=%.12e, deriv=%.12e\n",(func)(ptt),(ptt[j]-p[j])/((func)(ptt)-(*func)(p)));	printf("func(ptt)=%.12e, deriv=%.12e\n",(func)(ptt),(ptt[j]-p[j])/((func)(ptt)-(*func)(p)));
	fprintf(ficlog,"func(ptt)=%.12e, deriv=%.12e\n",(func)(ptt),(ptt[j]-p[j])/((func)(ptt)-(*func)(p)));
}	}
#endif	#endif

Line 652 void powell(double p[], double **xi, int	Line 693 void powell(double p[], double **xi, int
}	}
#ifdef DEBUG	#ifdef DEBUG
printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);	printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
for(j=1;j<=n;j++)	fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
	for(j=1;j<=n;j++){
printf(" %.12e",xit[j]);	printf(" %.12e",xit[j]);
	fprintf(ficlog," %.12e",xit[j]);
	}
printf("\n");	printf("\n");
	fprintf(ficlog,"\n");
#endif	#endif
}	}
}	}
Line 925 void mlikeli(FILE *ficres,double p[], in	Line 970 void mlikeli(FILE *ficres,double p[], in
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");	printf("Powell\n"); fprintf(ficlog,"Powell\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));	printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
	fprintf(ficlog,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
fprintf(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 949 void hesscov(double **matcov, double p[]	Line 995 void hesscov(double **matcov, double p[]
hess=matrix(1,npar,1,npar);	hess=matrix(1,npar,1,npar);

printf("\nCalculation of the hessian matrix. Wait...\n");	printf("\nCalculation of the hessian matrix. Wait...\n");
	fprintf(ficlog,"\nCalculation of the hessian matrix. Wait...\n");
for (i=1;i<=npar;i++){	for (i=1;i<=npar;i++){
printf("%d",i);fflush(stdout);	printf("%d",i);fflush(stdout);
	fprintf(ficlog,"%d",i);fflush(ficlog);
hess[i][i]=hessii(p,ftolhess,i,delti);	hess[i][i]=hessii(p,ftolhess,i,delti);
/printf(" %f ",p[i]);/	/printf(" %f ",p[i]);/
/printf(" %lf ",hess[i][i]);/	/printf(" %lf ",hess[i][i]);/
Line 960 void hesscov(double **matcov, double p[]	Line 1008 void hesscov(double **matcov, double p[]
for (j=1;j<=npar;j++) {	for (j=1;j<=npar;j++) {
if (j>i) {	if (j>i) {
printf(".%d%d",i,j);fflush(stdout);	printf(".%d%d",i,j);fflush(stdout);
	fprintf(ficlog,".%d%d",i,j);fflush(ficlog);
hess[i][j]=hessij(p,delti,i,j);	hess[i][j]=hessij(p,delti,i,j);
hess[j][i]=hess[i][j];	hess[j][i]=hess[i][j];
/printf(" %lf ",hess[i][j]);/	/printf(" %lf ",hess[i][j]);/
Line 967 void hesscov(double **matcov, double p[]	Line 1016 void hesscov(double **matcov, double p[]
}	}
}	}
printf("\n");	printf("\n");
	fprintf(ficlog,"\n");

printf("\nInverting the hessian to get the covariance matrix. Wait...\n");	printf("\nInverting the hessian to get the covariance matrix. Wait...\n");
	fprintf(ficlog,"\nInverting the hessian to get the covariance matrix. Wait...\n");

a=matrix(1,npar,1,npar);	a=matrix(1,npar,1,npar);
y=matrix(1,npar,1,npar);	y=matrix(1,npar,1,npar);
Line 988 void hesscov(double **matcov, double p[]	Line 1039 void hesscov(double **matcov, double p[]
}	}

printf("\n#Hessian matrix#\n");	printf("\n#Hessian matrix#\n");
	fprintf(ficlog,"\n#Hessian matrix#\n");
for (i=1;i<=npar;i++) {	for (i=1;i<=npar;i++) {
for (j=1;j<=npar;j++) {	for (j=1;j<=npar;j++) {
printf("%.3e ",hess[i][j]);	printf("%.3e ",hess[i][j]);
	fprintf(ficlog,"%.3e ",hess[i][j]);
}	}
printf("\n");	printf("\n");
	fprintf(ficlog,"\n");
}	}

/* Recompute Inverse */	/* Recompute Inverse */
Line 1009 void hesscov(double **matcov, double p[]	Line 1063 void hesscov(double **matcov, double p[]
for (i=1;i<=npar;i++){	for (i=1;i<=npar;i++){
y[i][j]=x[i];	y[i][j]=x[i];
printf("%.3e ",y[i][j]);	printf("%.3e ",y[i][j]);
	fprintf(ficlog,"%.3e ",y[i][j]);
}	}
printf("\n");	printf("\n");
	fprintf(ficlog,"\n");
}	}
*/	*/

Line 1052 double hessii( double x[], double delta,	Line 1108 double hessii( double x[], double delta,

#ifdef DEBUG	#ifdef DEBUG
printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);	printf("%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);
	fprintf(ficlog,"%d %d k1=%.12e k2=%.12e xk1=%.12e xk2=%.12e delt=%.12e res=%.12e l=%d k=%d,fx=%.12e\n",theta,theta,k1,k2,x[theta]+delt,x[theta]-delt,delt,res, l, k,fx);
#endif	#endif
/if(fabs(k1-2.0fx+k2) <1.e-13){ */	/if(fabs(k1-2.0fx+k2) <1.e-13){ */
if((k1 <khi/nkhi/2.) \|\| (k2 <khi/nkhi/2.)){	if((k1 <khi/nkhi/2.) \|\| (k2 <khi/nkhi/2.)){
Line 1099 double hessij( double x[], double delti[	Line 1156 double hessij( double x[], double delti[
res=(k1-k2-k3+k4)/4.0/delti[thetai]k/delti[thetaj]k/2.; /* Because of L not 2L /	res=(k1-k2-k3+k4)/4.0/delti[thetai]k/delti[thetaj]k/2.; /* Because of L not 2L /
#ifdef DEBUG	#ifdef DEBUG
printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);	printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
	fprintf(ficlog,"%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti/k=%.12e deltj/k=%.12e, xi-de/k=%.12e xj-de/k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4);
#endif	#endif
}	}
return res;	return res;
Line 1183 void freqsummary(char fileres[], int ag	Line 1241 void freqsummary(char fileres[], int ag
{ /* 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;
double **freq; / Frequencies */	double **freq; / Frequencies */
double *pp;	double *pp;
double pos, k2, dateintsum=0,k2cpt=0;	double pos, k2, dateintsum=0,k2cpt=0;
Line 1195 void freqsummary(char fileres[], int ag	Line 1254 void freqsummary(char fileres[], int ag
strcat(fileresp,fileres);	strcat(fileresp,fileres);
if((ficresp=fopen(fileresp,"w"))==NULL) {	if((ficresp=fopen(fileresp,"w"))==NULL) {
printf("Problem with prevalence resultfile: %s\n", fileresp);	printf("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,agemin,agemax+3);
Line 1202 void freqsummary(char fileres[], int ag	Line 1262 void freqsummary(char fileres[], int ag

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

	first=1;

for(k1=1; k1<=j;k1++){	for(k1=1; k1<=j;k1++){
for(i1=1; i1<=ncodemax[k1];i1++){	for(i1=1; i1<=ncodemax[k1];i1++){
j1++;	j1++;
Line 1254 void freqsummary(char fileres[], int ag	Line 1316 void freqsummary(char fileres[], int ag
fprintf(ficresp, "\n");	fprintf(ficresp, "\n");

for(i=(int)agemin; i <= (int)agemax+3; i++){	for(i=(int)agemin; i <= (int)agemax+3; i++){
if(i==(int)agemax+3)	if(i==(int)agemax+3){
printf("Total");	fprintf(ficlog,"Total");
else	}else{
printf("Age %d", i);	if(first==1){
	first=0;
	printf("See log file for details...\n");
	}
	fprintf(ficlog,"Age %d", i);
	}
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)	for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
pp[jk] += freq[jk][m][i];	pp[jk] += freq[jk][m][i];
Line 1265 void freqsummary(char fileres[], int ag	Line 1332 void freqsummary(char fileres[], int ag
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
for(m=-1, pos=0; m <=0 ; m++)	for(m=-1, pos=0; m <=0 ; m++)
pos += freq[jk][m][i];	pos += freq[jk][m][i];
if(pp[jk]>=1.e-10)	if(pp[jk]>=1.e-10){
	if(first==1){
printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);	printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
else	}
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);	fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
	}else{
	if(first==1)
	printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
	fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
	}
}	}

for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
Line 1279 void freqsummary(char fileres[], int ag	Line 1352 void freqsummary(char fileres[], int ag
for(jk=1,pos=0; jk <=nlstate ; jk++)	for(jk=1,pos=0; jk <=nlstate ; jk++)
pos += pp[jk];	pos += pp[jk];
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
if(pos>=1.e-5)	if(pos>=1.e-5){
printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);	if(first==1)
else	printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);	fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
	}else{
	if(first==1)
	printf(" %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 <= (int) agemax){
if(pos>=1.e-5){	if(pos>=1.e-5){
fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);	fprintf(ficresp," %d %.5f %.0f %.0f",i,pp[jk]/pos, pp[jk],pos);
Line 1296 void freqsummary(char fileres[], int ag	Line 1374 void freqsummary(char fileres[], int ag

for(jk=-1; jk <=nlstate+ndeath; jk++)	for(jk=-1; jk <=nlstate+ndeath; jk++)
for(m=-1; m <=nlstate+ndeath; m++)	for(m=-1; m <=nlstate+ndeath; m++)
if(freq[jk][m][i] !=0 ) printf(" %d%d=%.0f",jk,m,freq[jk][m][i]);	if(freq[jk][m][i] !=0 ) {
	if(first==1)
	printf(" %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 <= (int) agemax)
fprintf(ficresp,"\n");	fprintf(ficresp,"\n");
printf("\n");	if(first==1)
	printf("Others in log...\n");
	fprintf(ficlog,"\n");
}	}
}	}
}	}
Line 1386 void prevalence(int agemin, float agemax	Line 1470 void prevalence(int agemin, float agemax
probs[i][jk][j1]= pp[jk]/pos;	probs[i][jk][j1]= pp[jk]/pos;
}	}
}	}
}	}/* end jk */
	}/* end i */
}	} /* end i1 */
}	} /* 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,(int) agemin,(int) agemax+3);
Line 1412 void concatwav(int wav[], int **dh, int	Line 1495 void concatwav(int wav[], int **dh, int
int i, mi, m;	int i, mi, m;
/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;	/* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
double sum=0., jmean=0.;*/	double sum=0., jmean=0.;*/
	int first;
int j, k=0,jk, ju, jl;	int j, k=0,jk, ju, jl;
double sum=0.;	double sum=0.;
	first=0;
jmin=1e+5;	jmin=1e+5;
jmax=-1;	jmax=-1;
jmean=0.;	jmean=0.;
Line 1437 void concatwav(int wav[], int **dh, int	Line 1521 void concatwav(int wav[], int **dh, int
}	}

wav[i]=mi;	wav[i]=mi;
if(mi==0)	if(mi==0){
printf("Warning, no any valid information for:%d line=%d\n",num[i],i);	if(first==0){
	printf("Warning, no any valid information for:%d line=%d and may be others, see log file\n",num[i],i);
	first=1;
	}
	if(first==1){
	fprintf(ficlog,"Warning, no any valid information for:%d line=%d\n",num[i],i);
	}
	} /* end mi==0 */
}	}

for(i=1; i<=imx; i++){	for(i=1; i<=imx; i++){
Line 1479 void concatwav(int wav[], int **dh, int	Line 1570 void concatwav(int wav[], int **dh, int
}	}
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);
	fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);
}	}

/********* Tricode **************************/	/********* Tricode **************************/
void tricode(int Tvar, int *nbcode, int imx)	void tricode(int Tvar, int *nbcode, int imx)
{	{
Line 1519 void tricode(int Tvar, int *nbcode, in	Line 1612 void tricode(int Tvar, int *nbcode, in
for (k=0; k<19; k++) Ndum[k]=0;	for (k=0; k<19; k++) Ndum[k]=0;

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

ij=1;	ij=1;
for (i=1; i<=10; i++) {	for (i=1; i<=10; i++) {
Line 1531 void tricode(int Tvar, int *nbcode, in	Line 1624 void tricode(int Tvar, int *nbcode, in
}	}
}	}

cptcoveff=ij-1;	cptcoveff=ij-1;
}	}

/********* Health Expectancies **************/	/********* Health Expectancies **************/
Line 1641 void evsij(char fileres[], double ***eij	Line 1734 void evsij(char fileres[], double ***eij
}	}
}	}
}	}



for(j=1; j<= nlstate*2; j++)	for(j=1; j<= nlstate*2; 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];
}	}

}	}

/* End theta */	/* End theta */
Line 1658 void evsij(char fileres[], double ***eij	Line 1747 void evsij(char fileres[], double ***eij
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*2;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*2;i++)
for(j=1;j<=nlstate*2;j++)	for(j=1;j<=nlstate*2;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);
for(h=0;h<=nhstepm-1;h++){	fprintf(ficlog,"%d\|",(int)age);fflush(ficlog);
	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*2,1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,nlstate2,1,npar,1,nlstate2,gradg[k]);	matprod2(doldm,dnewm,1,nlstate2,1,npar,1,nlstate2,gradg[k]);
Line 1675 void evsij(char fileres[], double ***eij	Line 1765 void evsij(char fileres[], double ***eij
varhe[i][j][(int)age] += doldm[i][j]hfhf;	varhe[i][j][(int)age] += doldm[i][j]hfhf;
}	}
}	}


/* Computing expectancies */	/* Computing expectancies */
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 1702 void evsij(char fileres[], double ***eij	Line 1790 void evsij(char fileres[], double ***eij
free_ma3x(trgradg,0,nhstepm,1,nlstate*2,1,npar);	free_ma3x(trgradg,0,nhstepm,1,nlstate*2,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");
	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*2,1,npar);
free_matrix(doldm,1,nlstate2,1,nlstate2);	free_matrix(doldm,1,nlstate2,1,nlstate2);
Line 1709 void evsij(char fileres[], double ***eij	Line 1800 void evsij(char fileres[], double ***eij
}	}

/********** Variance ****************/	/********** Variance ****************/
void varevsij(char fileres[], 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 ij, int estepm)	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 ij, int estepm, int cptcov, int cptcod, int popbased)
{	{
/* 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;/
double dnewm,doldm;	double dnewm,doldm;
	double dnewmp,doldmp;
int i, j, nhstepm, hstepm, h, nstepm ;	int i, j, nhstepm, hstepm, h, nstepm ;
int k, cptcode;	int k, cptcode;
double *xp;	double *xp;
double gp, gm;	double gp, gm; /* for var eij */
double *gradg, *trgradg;	double *gradg, trgradg; /for var eij */
	double gradgp, trgradgp; /* 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 ***p3mat;	double ***p3mat;
double age,agelim, hf;	double age,agelim, hf;
int theta;	int theta;
	char digit[4];
	char digitp[16];

	char fileresprobmorprev[FILENAMELENGTH];

	if(popbased==1)
	strcpy(digitp,"-populbased-");
	else
	strcpy(digitp,"-stablbased-");

	strcpy(fileresprobmorprev,"prmorprev");
	sprintf(digit,"%-d",ij);
	/printf("DIGIT=%s, ij=%d ijr=%-d\|\n",digit, ij,ij);/
	strcat(fileresprobmorprev,digit); /* Tvar to be done */
	strcat(fileresprobmorprev,digitp); /* Popbased or not */
	strcat(fileresprobmorprev,fileres);
	if((ficresprobmorprev=fopen(fileresprobmorprev,"w"))==NULL) {
	printf("Problem with resultfile: %s\n", fileresprobmorprev);
	fprintf(ficlog,"Problem with resultfile: %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(ficresprobmorprev,"# probabilities of dying during a year and weighted mean w1p1j+w2p2j+... stand dev in()\n");
	fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
	for(j=nlstate+1; j<=(nlstate+ndeath);j++){
	fprintf(ficresprobmorprev," p.%-d SE",j);
	for(i=1; i<=nlstate;i++)
	fprintf(ficresprobmorprev," w%1d p%-d%-d",i,i,j);
	}
	fprintf(ficresprobmorprev,"\n");
	if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
	printf("Problem with gnuplot file: %s\n", optionfilegnuplot);
	fprintf(ficlog,"Problem with gnuplot file: %s\n", optionfilegnuplot);
	exit(0);
	}
	else{
	fprintf(ficgp,"\n# Routine varevsij");
	}
	if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
	printf("Problem with html file: %s\n", optionfilehtm);
	fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);
	exit(0);
	}
	else{
	fprintf(fichtm,"\n<li><h4> Computing step probabilities of dying and weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
	}
	varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);

fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n");	fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n");
fprintf(ficresvij,"# Age");	fprintf(ficresvij,"# Age");
Line 1734 void varevsij(char fileres[], double ***	Line 1876 void varevsij(char fileres[], double ***
xp=vector(1,npar);	xp=vector(1,npar);
dnewm=matrix(1,nlstate,1,npar);	dnewm=matrix(1,nlstate,1,npar);
doldm=matrix(1,nlstate,1,nlstate);	doldm=matrix(1,nlstate,1,nlstate);
	dnewmp= matrix(nlstate+1,nlstate+ndeath,1,npar);
	doldmp= matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);

	gradgp=matrix(1,npar,nlstate+1,nlstate+ndeath);
	gpp=vector(nlstate+1,nlstate+ndeath);
	gmp=vector(nlstate+1,nlstate+ndeath);
	trgradgp =matrix(nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/

if(estepm < stepm){	if(estepm < stepm){
printf ("Problem %d lower than %d\n",estepm, stepm);	printf ("Problem %d lower than %d\n",estepm, stepm);
Line 1761 void varevsij(char fileres[], double ***	Line 1910 void varevsij(char fileres[], double ***
gp=matrix(0,nhstepm,1,nlstate);	gp=matrix(0,nhstepm,1,nlstate);
gm=matrix(0,nhstepm,1,nlstate);	gm=matrix(0,nhstepm,1,nlstate);


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 */
xp[i] = x[i] + (i==theta ?delti[theta]:0);	xp[i] = x[i] + (i==theta ?delti[theta]:0);
Line 1779 void varevsij(char fileres[], double ***	Line 1929 void varevsij(char fileres[], double ***
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 */
	for(j=nlstate+1,gpp[j]=0.;j<=nlstate+ndeath;j++){
	for(i=1; i<= nlstate; i++)
	gpp[j] += prlim[i][i]*p3mat[i][j][1];
	}
	/* end force of mortality */

for(i=1; i<=npar; i++) /* Computes gradient */	for(i=1; i<=npar; i++) /* Computes gradient */
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 1796 void varevsij(char fileres[], double ***	Line 1952 void varevsij(char fileres[], double ***
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 */
	for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){
	for(i=1; i<= nlstate; i++)
	gmp[j] += prlim[i][i]*p3mat[i][j][1];
	}
	/* end force of mortality */

for(j=1; j<= nlstate; j++)	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 */
	gradgp[theta][j]= (gpp[j]-gmp[j])/2./delti[theta];
	}

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

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

for(h=0; h<=nhstepm; h++)	for(h=0; h<=nhstepm; h++) /* veij */
for(j=1; j<=nlstate;j++)	for(j=1; j<=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(j=nlstate+1; j<=nlstate+ndeath;j++) /* mu */
	for(theta=1; theta <=npar; theta++)
	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++)
for(j=1;j<=nlstate;j++)	for(j=1;j<=nlstate;j++)
Line 1825 void varevsij(char fileres[], double ***	Line 1995 void varevsij(char fileres[], double ***
}	}
}	}

	/* pptj */
	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);
	for(j=nlstate+1;j<=nlstate+ndeath;j++)
	for(i=nlstate+1;i<=nlstate+ndeath;i++)
	varppt[j][i]=doldmp[j][i];
	/* end ppptj */
	hpxij(p3mat,nhstepm,age,hstepm,x,nlstate,stepm,oldm,savm, ij);
	prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ij);

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

	/* This for computing force of mortality (h=1)as a weighted average */
	for(j=nlstate+1,gmp[j]=0.;j<=nlstate+ndeath;j++){
	for(i=1; i<= nlstate; i++)
	gmp[j] += prlim[i][i]*p3mat[i][j][1];
	}
	/* end force of mortality */

	fprintf(ficresprobmorprev,"%3d %d ",(int) age, ij);
	for(j=nlstate+1; j<=(nlstate+ndeath);j++){
	fprintf(ficresprobmorprev," %11.3e %11.3e",gmp[j], sqrt(varppt[j][j]));
	for(i=1; i<=nlstate;i++){
	fprintf(ficresprobmorprev," %11.3e %11.3e ",prlim[i][i],p3mat[i][j][1]);
	}
	}
	fprintf(ficresprobmorprev,"\n");

fprintf(ficresvij,"%.0f ",age );	fprintf(ficresvij,"%.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 1837 void varevsij(char fileres[], double ***	Line 2038 void varevsij(char fileres[], double ***
free_ma3x(trgradg,0,nhstepm,1,nlstate,1,npar);	free_ma3x(trgradg,0,nhstepm,1,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);
} /* End age */	} /* End age */
	free_vector(gpp,nlstate+1,nlstate+ndeath);
	free_vector(gmp,nlstate+1,nlstate+ndeath);
	free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
	free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
	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 */
	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 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(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 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);

free_vector(xp,1,npar);	free_vector(xp,1,npar);
free_matrix(doldm,1,nlstate,1,npar);	free_matrix(doldm,1,nlstate,1,nlstate);
free_matrix(dnewm,1,nlstate,1,nlstate);	free_matrix(dnewm,1,nlstate,1,npar);
	free_matrix(doldmp,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
	free_matrix(dnewmp,nlstate+1,nlstate+ndeath,1,npar);
	free_matrix(varppt,nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
	fclose(ficresprobmorprev);
	fclose(ficgp);
	fclose(fichtm);

}	}

Line 1928 void varprevlim(char fileres[], double *	Line 2148 void varprevlim(char fileres[], double *
}	}

/********** Variance of one-step probabilities ****************/	/********** Variance of one-step probabilities ****************/
void varprob(char fileres[], double *matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int Tvar, int *nbcode, int ncodemax)	void varprob(char optionfilefiname[], double *matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int Tvar, int *nbcode, int ncodemax)
{	{
int i, j, i1, k1, j1, z1;	int i, j=0, i1, k1, l1, t, tj;
int k=0, cptcode;	int k2, l2, j1, z1;
	int k=0,l, cptcode;
	int first=1, first1;
	double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2;
double dnewm,doldm;	double dnewm,doldm;
double *xp;	double *xp;
double gp, gm;	double gp, gm;
double gradg, trgradg;	double gradg, trgradg;
	double **mu;
double age,agelim, cov[NCOVMAX];	double age,agelim, cov[NCOVMAX];
	double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
int theta;	int theta;
char fileresprob[FILENAMELENGTH];	char fileresprob[FILENAMELENGTH];
	char fileresprobcov[FILENAMELENGTH];
	char fileresprobcor[FILENAMELENGTH];

	double ***varpij;

strcpy(fileresprob,"prob");	strcpy(fileresprob,"prob");
strcat(fileresprob,fileres);	strcat(fileresprob,fileres);
if((ficresprob=fopen(fileresprob,"w"))==NULL) {	if((ficresprob=fopen(fileresprob,"w"))==NULL) {
printf("Problem with resultfile: %s\n", fileresprob);	printf("Problem with resultfile: %s\n", fileresprob);
	fprintf(ficlog,"Problem with resultfile: %s\n", fileresprob);
	}
	strcpy(fileresprobcov,"probcov");
	strcat(fileresprobcov,fileres);
	if((ficresprobcov=fopen(fileresprobcov,"w"))==NULL) {
	printf("Problem with resultfile: %s\n", fileresprobcov);
	fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcov);
	}
	strcpy(fileresprobcor,"probcor");
	strcat(fileresprobcor,fileres);
	if((ficresprobcor=fopen(fileresprobcor,"w"))==NULL) {
	printf("Problem with resultfile: %s\n", fileresprobcor);
	fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobcor);
}	}
printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);	printf("Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
	fprintf(ficlog,"Computing standard deviation of one-step probabilities: result on file '%s' \n",fileresprob);
	printf("Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
	fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
	printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
	fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);

fprintf(ficresprob,"#One-step probabilities and standard deviation in parentheses\n");	fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
fprintf(ficresprob,"# Age");	fprintf(ficresprob,"# Age");
for(i=1; i<=nlstate;i++)	fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
for(j=1; j<=(nlstate+ndeath);j++)	fprintf(ficresprobcov,"# Age");
fprintf(ficresprob," p%1d-%1d (SE)",i,j);	fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
	fprintf(ficresprobcov,"# Age");


	for(i=1; i<=nlstate;i++)
	for(j=1; j<=(nlstate+ndeath);j++){
	fprintf(ficresprob," p%1d-%1d (SE)",i,j);
	fprintf(ficresprobcov," p%1d-%1d ",i,j);
	fprintf(ficresprobcor," p%1d-%1d ",i,j);
	}
fprintf(ficresprob,"\n");	fprintf(ficresprob,"\n");
	fprintf(ficresprobcov,"\n");
	fprintf(ficresprobcor,"\n");
	xp=vector(1,npar);
	dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
	doldm=matrix(1,(nlstate)(nlstate+ndeath),1,(nlstate)(nlstate+ndeath));
	mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);
	varpij=ma3x(1,nlstate(nlstate+ndeath),1,nlstate(nlstate+ndeath),(int) bage, (int) fage);
	first=1;
	if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
	printf("Problem with gnuplot file: %s\n", optionfilegnuplot);
	fprintf(ficlog,"Problem with gnuplot file: %s\n", optionfilegnuplot);
	exit(0);
	}
	else{
	fprintf(ficgp,"\n# Routine varprob");
	}
	if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
	printf("Problem with html file: %s\n", optionfilehtm);
	fprintf(ficlog,"Problem with html file: %s\n", optionfilehtm);
	exit(0);
	}
	else{
	fprintf(fichtm,"\n<li><h4> Computing and drawing one step probabilities with their confidence intervals</h4></li>\n");
	fprintf(fichtm,"\n");

	fprintf(fichtm,"\n<li><h4> Computing matrix of variance-covariance of step probabilities</h4></li>\n");
	fprintf(fichtm,"\nWe have drawn ellipsoids of confidence around the p<inf>ij</inf>, p<inf>kl</inf> to understand the covariance between two incidences. They are expressed in year<sup>-1</sup> in order to be less dependent of stepm.<br>\n");
	fprintf(fichtm,"\n<br> We have drawn x'cov<sup>-1</sup>x = 4 where x is the column vector (pij,pkl). It means that if pij and pkl where uncorrelated the (2X2) matrix would have been (1/(var pij), 0 , 0, 1/(var pkl)), and the confidence interval would be 2 standard deviations wide on each axis. <br> When both incidences are correlated we diagonalised the inverse of the covariance matrix and made the appropriate rotation.<br> \n");

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

cov[1]=1;	cov[1]=1;
j=cptcoveff;	tj=cptcoveff;
if (cptcovn<1) {j=1;ncodemax[1]=1;}	if (cptcovn<1) {tj=1;ncodemax[1]=1;}
j1=0;	j1=0;
for(k1=1; k1<=1;k1++){	for(t=1; t<=tj;t++){
for(i1=1; i1<=ncodemax[k1];i1++){	for(i1=1; i1<=ncodemax[t];i1++){
j1++;	j1++;

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#");
}	fprintf(ficresprobcov, "\n#********** Variable ");
	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprobcov, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
	fprintf(ficresprobcov, "**********\n#");

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


	fprintf(fichtm, "\n<hr size=\"2\" color=\"#EC5E5E\">********** Variable ");
	for (z1=1; z1<=cptcoveff; z1++) fprintf(fichtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
	fprintf(fichtm, "**********\n<hr size=\"2\" color=\"#EC5E5E\">");

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

for (age=bage; age<=fage; age ++){	for (age=bage; age<=fage; age ++){
cov[2]=age;	cov[2]=age;
for (k=1; k<=cptcovn;k++) {	for (k=1; k<=cptcovn;k++) {
cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];	cov[2+k]=nbcode[Tvar[k]][codtab[j1][Tvar[k]]];

}	}
for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];	for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
for (k=1; k<=cptcovprod;k++)	for (k=1; k<=cptcovprod;k++)
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];	cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];

gradg=matrix(1,npar,1,9);	gradg=matrix(1,npar,1,(nlstate)*(nlstate+ndeath));
trgradg=matrix(1,9,1,npar);	trgradg=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
gp=vector(1,(nlstate+ndeath)*(nlstate+ndeath));	gp=vector(1,(nlstate)*(nlstate+ndeath));
gm=vector(1,(nlstate+ndeath)*(nlstate+ndeath));	gm=vector(1,(nlstate)*(nlstate+ndeath));

for(theta=1; theta <=npar; theta++){	for(theta=1; theta <=npar; theta++){
for(i=1; i<=npar; i++)	for(i=1; i<=npar; i++)
Line 1997 fprintf(ficresprob,"#One-step probabilit	Line 2292 fprintf(ficresprob,"#One-step probabilit
pmij(pmmij,cov,ncovmodel,xp,nlstate);	pmij(pmmij,cov,ncovmodel,xp,nlstate);

k=0;	k=0;
for(i=1; i<= (nlstate+ndeath); i++){	for(i=1; i<= (nlstate); i++){
for(j=1; j<=(nlstate+ndeath);j++){	for(j=1; j<=(nlstate+ndeath);j++){
k=k+1;	k=k+1;
gp[k]=pmmij[i][j];	gp[k]=pmmij[i][j];
Line 2009 fprintf(ficresprob,"#One-step probabilit	Line 2304 fprintf(ficresprob,"#One-step probabilit

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

for(i=1; i<= (nlstate+ndeath)*(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])/2./delti[theta];
}	}

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

matprod2(dnewm,trgradg,1,9,1,npar,1,npar,matcov);	matprod2(dnewm,trgradg,1,(nlstate)*(nlstate+ndeath),1,npar,1,npar,matcov);
matprod2(doldm,dnewm,1,9,1,npar,1,9,gradg);	matprod2(doldm,dnewm,1,(nlstate)(nlstate+ndeath),1,npar,1,(nlstate)(nlstate+ndeath),gradg);

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

k=0;	k=0;
for(i=1; i<=(nlstate+ndeath); i++){	for(i=1; i<=(nlstate); i++){
for(j=1; j<=(nlstate+ndeath);j++){	for(j=1; j<=(nlstate+ndeath);j++){
k=k+1;	k=k+1;
gm[k]=pmmij[i][j];	mu[k][(int) age]=pmmij[i][j];
}	}
}	}
	for(i=1;i<=(nlstate)*(nlstate+ndeath);i++)
/*printf("\n%d ",(int)age);	for(j=1;j<=(nlstate)*(nlstate+ndeath);j++)
for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++){	varpij[i][j][(int)age] = doldm[i][j];

	/*printf("\n%d ",(int)age);
	for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
printf("%e [%e ;%e] ",gm[i],gm[i]-2sqrt(doldm[i][i]),gm[i]+2sqrt(doldm[i][i]));	printf("%e [%e ;%e] ",gm[i],gm[i]-2sqrt(doldm[i][i]),gm[i]+2sqrt(doldm[i][i]));
	fprintf(ficlog,"%e [%e ;%e] ",gm[i],gm[i]-2sqrt(doldm[i][i]),gm[i]+2sqrt(doldm[i][i]));
}*/	}*/

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

for (i=1; i<=(nlstate+ndeath)*(nlstate+ndeath-1);i++)	for (i=1; i<=(nlstate)*(nlstate+ndeath);i++)
fprintf(ficresprob,"%.3e (%.3e) ",gm[i],sqrt(doldm[i][i]));	fprintf(ficresprob,"%11.3e (%11.3e) ",mu[i][(int) age],sqrt(varpij[i][i][(int)age]));
	for (i=1; i<=(nlstate)*(nlstate+ndeath);i++){
}	fprintf(ficresprobcov,"%11.3e ",mu[i][(int) age]);
}	fprintf(ficresprobcor,"%11.3e ",mu[i][(int) age]);
	}
	i=0;
	for (k=1; k<=(nlstate);k++){
	for (l=1; l<=(nlstate+ndeath);l++){
	i=i++;
	fprintf(ficresprobcov,"\n%d %d-%d",(int)age,k,l);
	fprintf(ficresprobcor,"\n%d %d-%d",(int)age,k,l);
	for (j=1; j<=i;j++){
	fprintf(ficresprobcov," %11.3e",varpij[i][j][(int)age]);
	fprintf(ficresprobcor," %11.3e",varpij[i][j][(int) age]/sqrt(varpij[i][i][(int) age])/sqrt(varpij[j][j][(int)age]));
	}
	}
	}/* end of loop for state */
	} /* end of loop for age */

	/* Confidence intervalle of pij */
	/*
	fprintf(ficgp,"\nset noparametric;unset label");
	fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
	fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
	fprintf(fichtm,"\n<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
	fprintf(fichtm,"\n<br><img src=\"pijgr%s.png\"> ",optionfilefiname);
	fprintf(ficgp,"\nset out \"pijgr%s.png\"",optionfilefiname);
	fprintf(ficgp,"\nplot \"%s\" every :::%d::%d u 1:2 \"\%%lf",k1,k2,xfilevarprob);
	*/

	/* Drawing ellipsoids of confidence of two variables p(k1-l1,k2-l2)*/
	first1=1;
	for (k1=1; k1<=(nlstate);k1++){
	for (l1=1; l1<=(nlstate+ndeath);l1++){
	if(l1==k1) continue;
	i=(k1-1)*(nlstate+ndeath)+l1;
	for (k2=1; k2<=(nlstate);k2++){
	for (l2=1; l2<=(nlstate+ndeath);l2++){
	if(l2==k2) continue;
	j=(k2-1)*(nlstate+ndeath)+l2;
	if(j<=i) continue;
	for (age=bage; age<=fage; age ++){
	if ((int)age %5==0){
	v1=varpij[i][i][(int)age]/stepmYEARM/stepmYEARM;
	v2=varpij[j][j][(int)age]/stepmYEARM/stepmYEARM;
	cv12=varpij[i][j][(int)age]/stepmYEARM/stepmYEARM;
	mu1=mu[i][(int) age]/stepm*YEARM ;
	mu2=mu[j][(int) age]/stepm*YEARM;
	/* Computing eigen value of matrix of covariance */
	lc1=(v1+v2)+sqrt((v1+v2)(v1+v2) - 4(v1v2-cv12cv12));
	lc2=(v1+v2)-sqrt((v1+v2)(v1+v2) - 4(v1v2-cv12cv12));
	if(first1==1){
	first1=0;
	printf("Var %.4e %.4e cov %.4e Eigen %.3e %.3e\nOthers in log...\n",v1,v2,cv12,lc1,lc2);
	}
	fprintf(ficlog,"Var %.4e %.4e cov %.4e Eigen %.3e %.3e\n",v1,v2,cv12,lc1,lc2);
	/* Eigen vectors */
	v11=(1./sqrt(1+(v1-lc1)*(v1-lc1)/cv12/cv12));
	v21=sqrt(1.-v11*v11);
	v12=-v21;
	v22=v11;
	/printf(fignu/
	/* mu1+ v11lc1cost + v12lc2sin(t) */
	/* mu2+ v21lc1cost + v21lc2sin(t) */
	if(first==1){
	first=0;
	fprintf(ficgp,"\nset parametric;set nolabel");
	fprintf(ficgp,"\nset log y;set log x; set xlabel \"p%1d%1d (year-1)\";set ylabel \"p%1d%1d (year-1)\"",k2,l2,k1,l1);
	fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
	fprintf(fichtm,"\n<br>Ellipsoids of confidence cov(p%1d%1d,p%1d%1d) expressed in year<sup>-1</sup> :<a href=\"varpijgr%s%d%1d%1d-%1d%1d.png\">varpijgr%s%d%1d%1d-%1d%1d.png</A>, ",k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1,optionfilefiname, j1,k2,l2,k1,l1);
	fprintf(fichtm,"\n<br><img src=\"varpijgr%s%d%1d%1d-%1d%1d.png\"> ",optionfilefiname, j1,k2,l2,k1,l1);
	fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\"",optionfilefiname, j1,k2,l2,k1,l1);
	fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
	fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
	/* 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)) t \"%d\"",\
	mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
	mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
	*/
	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",\
	mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
	mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));
	}else{
	first=0;
	fprintf(ficgp,"\n# Age %d, p%1d%1d - p%1d%1d",(int) age, k2,l2,k1,l1);
	fprintf(ficgp,"\nset label \"%d\" at %11.3e,%11.3e center",(int) age, mu2,mu1);
	/*
	fprintf(ficgp,"\nreplot %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)) t \"%d\"",\
	mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
	mu1,std,v11,sqrt(lc1),v12,sqrt(lc2),(int) age);
	*/
	fprintf(ficgp,"\nreplot %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",\
	mu2,std,v21,sqrt(lc1),v21,sqrt(lc2), \
	mu1,std,v11,sqrt(lc1),v12,sqrt(lc2));
	}/* if first */
	} /* age mod 5 */
	} /* end loop age */
	fprintf(ficgp,"\nset out \"varpijgr%s%d%1d%1d-%1d%1d.png\";replot;",optionfilefiname, j1,k2,l2,k1,l1);
	first=1;
	} /l12 /
	} /* k12 */
	} /l1 /
	}/* k1 */
	} /* loop covariates */
	free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));	free_vector(gp,1,(nlstate+ndeath)*(nlstate+ndeath));
free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));	free_vector(gm,1,(nlstate+ndeath)*(nlstate+ndeath));
	free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);	free_matrix(trgradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);	free_matrix(gradg,1,(nlstate+ndeath)*(nlstate+ndeath),1,npar);
}	}
free_vector(xp,1,npar);	free_vector(xp,1,npar);
fclose(ficresprob);	fclose(ficresprob);
	fclose(ficresprobcov);
	fclose(ficresprobcor);
	fclose(ficgp);
	fclose(fichtm);
}	}


/***************** Printing html file *********/	/***************** Printing html file *********/
void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \	void printinghtml(char fileres[], char title[], char datafile[], int firstpass, \
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 optionfile[], \	int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
char optionfilehtm[],char rfileres[], char optionfilegnuplot[],\	int popforecast, int estepm ,\
char version[], int popforecast, int estepm ,/* \ */
double jprev1, double mprev1,double anprev1, \	double jprev1, double mprev1,double anprev1, \
double jprev2, double mprev2,double anprev2){	double jprev2, double mprev2,double anprev2){
int jj1, k1, i1, cpt;	int jj1, k1, i1, cpt;
FILE *fichtm;
/char optionfilehtm[FILENAMELENGTH];/	/char optionfilehtm[FILENAMELENGTH];/
	if((fichtm=fopen(optionfilehtm,"a"))==NULL) {
strcpy(optionfilehtm,optionfile);
strcat(optionfilehtm,".htm");
if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
printf("Problem with %s \n",optionfilehtm), exit(0);	printf("Problem with %s \n",optionfilehtm), exit(0);
	fprintf(ficlog,"Problem with %s \n",optionfilehtm), exit(0);
}	}

fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n	fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
\n
Total number of observations=%d <br>\n
Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
<hr size=\"2\" color=\"#EC5E5E\">
<ul><li>Parameter files<br>\n
- Copy of the parameter file: <a href=\"o%s\">o%s</a><br>\n
- Gnuplot file name: <a href=\"%s\">%s</a><br></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,optionfilegnuplot,optionfilegnuplot);

fprintf(fichtm,"<ul><li>Result files (first order: no variance)<br>\n
- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n	- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"p%s\">p%s</a> <br>\n
- Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n	- Estimated transition probabilities over %d (stepm) months: <a href=\"pij%s\">pij%s</a><br>\n
- Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n	- Stable prevalence in each health state: <a href=\"pl%s\">pl%s</a> <br>\n
Line 2095 Interval (in months) between two waves:	Line 2487 Interval (in months) between two waves:
<a href=\"e%s\">e%s</a> <br>\n</li>", \	<a href=\"e%s\">e%s</a> <br>\n</li>", \
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);	jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,fileres,fileres,stepm,fileres,fileres,fileres,fileres,estepm,fileres,fileres);

fprintf(fichtm,"\n<li> Result files (second order: variances)<br>\n	fprintf(fichtm," \n<ul><li><b>Graphs</b></li><p>");

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

	jj1=0;
	for(k1=1; k1<=m;k1++){
	for(i1=1; i1<=ncodemax[k1];i1++){
	jj1++;
	if (cptcovn > 0) {
	fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
	for (cpt=1; cpt<=cptcoveff;cpt++)
	fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
	fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
	}
	/* 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>
	<img src=\"pe%s%d1.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
	/* 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>
	<img src=\"pe%s%d2.png\">",stepm,strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
	/* Stable prevalence in each health state */
	for(cpt=1; cpt<nlstate;cpt++){
	fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>
	<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
	}
	for(cpt=1; cpt<=nlstate;cpt++) {
	fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>
	<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
	}
	fprintf(fichtm,"\n<br>- Total life expectancy by age and
	health expectancies in states (1) and (2): e%s%d.png<br>
	<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
	} /* end i1 */
	}/* End k1 */
	fprintf(fichtm,"</ul>");


	fprintf(fichtm,"\n<br><li><h4> Result files (second order: variances)</h4>\n
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n	- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n
- Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n	- Variance of one-step probabilities: <a href=\"prob%s\">prob%s</a> <br>\n
	- Variance-covariance of one-step probabilities: <a href=\"probcov%s\">probcov%s</a> <br>\n
	- Correlation matrix of one-step probabilities: <a href=\"probcor%s\">probcor%s</a> <br>\n
- Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n	- Variances and covariances of life expectancies by age and initial health status (estepm=%d months): <a href=\"v%s\">v%s</a><br>\n
- 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, 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
Line 2108 Interval (in months) between two waves:	Line 2540 Interval (in months) between two waves:
<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," <li>Graphs</li><p>");	fprintf(fichtm," <ul><li><b>Graphs</b></li><p>");

m=cptcoveff;	m=cptcoveff;
if (cptcovn < 1) {m=1;ncodemax[1]=1;}	if (cptcovn < 1) {m=1;ncodemax[1]=1;}
Line 2116 fprintf(fichtm," <li>Graphs</li><p>");	Line 2548 fprintf(fichtm," <li>Graphs</li><p>");
jj1=0;	jj1=0;
for(k1=1; k1<=m;k1++){	for(k1=1; k1<=m;k1++){
for(i1=1; i1<=ncodemax[k1];i1++){	for(i1=1; i1<=ncodemax[k1];i1++){
jj1++;	jj1++;
if (cptcovn > 0) {	if (cptcovn > 0) {
fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");	fprintf(fichtm,"<hr size=\"2\" color=\"#EC5E5E\">************ Results for covariates");
for (cpt=1; cpt<=cptcoveff;cpt++)	for (cpt=1; cpt<=cptcoveff;cpt++)
fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);	fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[jj1][cpt]]);
fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");	fprintf(fichtm," ************\n<hr size=\"2\" color=\"#EC5E5E\">");
}	}
/* Pij */	for(cpt=1; cpt<=nlstate;cpt++) {
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\">",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1);
/* 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>
<img src=\"pe%s%d2.png\">",strtok(optionfile, "."),jj1,stepm,strtok(optionfile, "."),jj1);
/* Stable prevalence in each health state */
for(cpt=1; cpt<nlstate;cpt++){
fprintf(fichtm,"<br>- Stable prevalence in each health state : p%s%d%d.png<br>
<img src=\"p%s%d%d.png\">",strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);
}
for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"<br>- Observed and stationary prevalence (with confident	fprintf(fichtm,"<br>- Observed and stationary 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);
}	}
for(cpt=1; cpt<=nlstate;cpt++) {	} /* end i1 */
fprintf(fichtm,"\n<br>- Health life expectancies by age and initial health state (%d): exp%s%d%d.png <br>	}/* End k1 */
<img src=\"exp%s%d%d.png\">",cpt,strtok(optionfile, "."),cpt,jj1,strtok(optionfile, "."),cpt,jj1);	fprintf(fichtm,"</ul>");
}
fprintf(fichtm,"\n<br>- Total life expectancy by age and
health expectancies in states (1) and (2): e%s%d.png<br>
<img src=\"e%s%d.png\">",strtok(optionfile, "."),jj1,strtok(optionfile, "."),jj1);
fprintf(fichtm,"\n</body>");
}
}
fclose(fichtm);	fclose(fichtm);
}	}

/***************** Gnuplot file ************/	/***************** Gnuplot file ************/
void printinggnuplot(char fileres[],char optionfilefiname[],char optionfile[],char optionfilegnuplot[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){	void printinggnuplot(char fileres[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){

int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;	int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
int ng;	int ng;
strcpy(optionfilegnuplot,optionfilefiname);	if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) {
strcat(optionfilegnuplot,".gp.txt");
if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
printf("Problem with file %s",optionfilegnuplot);	printf("Problem with file %s",optionfilegnuplot);
	fprintf(ficlog,"Problem with file %s",optionfilegnuplot);
}	}

#ifdef windows	#ifdef windows
Line 2280 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);	Line 2693 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
for(k=1; k <=(nlstate+ndeath); k++){	for(k=1; k <=(nlstate+ndeath); k++){
if (k != i) {	if (k != i) {
for(j=1; j <=ncovmodel; j++){	for(j=1; j <=ncovmodel; j++){

fprintf(ficgp,"p%d=%f ",jk,p[jk]);	fprintf(ficgp,"p%d=%f ",jk,p[jk]);
jk++;	jk++;
fprintf(ficgp,"\n");	fprintf(ficgp,"\n");
Line 2303 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);	Line 2715 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
for(k=1; k<=(nlstate+ndeath); k++) {	for(k=1; k<=(nlstate+ndeath); k++) {
if (k != k2){	if (k != k2){
if(ng==2)	if(ng==2)
fprintf(ficgp," %fexp(p%d+p%dx",stepm/YEARM,i,i+1);	fprintf(ficgp," %fexp(p%d+p%dx",YEARM/stepm,i,i+1);
else	else
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);	fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
ij=1;	ij=1;
Line 2334 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);	Line 2746 fprintf(ficgp,"\" t \"e%d1\" w l",cpt);
if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");	if ((k+k2)!= (nlstate*2+ndeath)) fprintf(ficgp,",");
i=i+ncovmodel;	i=i+ncovmodel;
}	}
}	} /* end k */
}	} /* end k2 */
}	} /* end jk */
}	} /* end ng */
fclose(ficgp);	fclose(ficgp);
} /* end gnuplot */	} /* end gnuplot */


Line 2385 calagedate=(anproj1+mproj1/12.+jproj1/36	Line 2797 calagedate=(anproj1+mproj1/12.+jproj1/36
strcat(fileresf,fileres);	strcat(fileresf,fileres);
if((ficresf=fopen(fileresf,"w"))==NULL) {	if((ficresf=fopen(fileresf,"w"))==NULL) {
printf("Problem with forecast resultfile: %s\n", fileresf);	printf("Problem with forecast resultfile: %s\n", fileresf);
	fprintf(ficlog,"Problem with forecast resultfile: %s\n", fileresf);
}	}
printf("Computing forecasting: result on file '%s' \n", fileresf);	printf("Computing forecasting: result on file '%s' \n", fileresf);
	fprintf(ficlog,"Computing forecasting: result on file '%s' \n", fileresf);

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

Line 2489 populforecast(char fileres[], double anp	Line 2903 populforecast(char fileres[], double anp
strcat(filerespop,fileres);	strcat(filerespop,fileres);
if((ficrespop=fopen(filerespop,"w"))==NULL) {	if((ficrespop=fopen(filerespop,"w"))==NULL) {
printf("Problem with forecast resultfile: %s\n", filerespop);	printf("Problem with forecast resultfile: %s\n", filerespop);
	fprintf(ficlog,"Problem with forecast resultfile: %s\n", filerespop);
}	}
printf("Computing forecasting: result on file '%s' \n", filerespop);	printf("Computing forecasting: result on file '%s' \n", filerespop);
	fprintf(ficlog,"Computing forecasting: result on file '%s' \n", filerespop);

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

Line 2510 populforecast(char fileres[], double anp	Line 2926 populforecast(char fileres[], double anp
if (popforecast==1) {	if (popforecast==1) {
if((ficpop=fopen(popfile,"r"))==NULL) {	if((ficpop=fopen(popfile,"r"))==NULL) {
printf("Problem with population file : %s\n",popfile);exit(0);	printf("Problem with population file : %s\n",popfile);exit(0);
	fprintf(ficlog,"Problem with population file : %s\n",popfile);exit(0);
}	}
popage=ivector(0,AGESUP);	popage=ivector(0,AGESUP);
popeffectif=vector(0,AGESUP);	popeffectif=vector(0,AGESUP);
Line 2638 int main(int argc, char *argv[])	Line 3055 int main(int argc, char *argv[])
double ***p3mat;	double ***p3mat;
int *indx;	int *indx;
char line[MAXLINE], linepar[MAXLINE];	char line[MAXLINE], linepar[MAXLINE];
char title[MAXLINE];	char path[80],pathc[80],pathcd[80],pathtot[80],model[80];
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
char optionfilext[10], optionfilefiname[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilegnuplot[FILENAMELENGTH], plotcmd[FILENAMELENGTH];

char fileres[FILENAMELENGTH], filerespij[FILENAMELENGTH], filereso[FILENAMELENGTH], rfileres[FILENAMELENGTH];

char filerest[FILENAMELENGTH];
char fileregp[FILENAMELENGTH];
char popfile[FILENAMELENGTH];
char path[80],pathc[80],pathcd[80],pathtot[80],model[20];
int firstobs=1, lastobs=10;	int firstobs=1, lastobs=10;
int sdeb, sfin; /* Status at beginning and end */	int sdeb, sfin; /* Status at beginning and end */
int c, h , cpt,l;	int c, h , cpt,l;
Line 2674 int main(int argc, char *argv[])	Line 3082 int main(int argc, char *argv[])
double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;	double dateprev1, dateprev2,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2;


char version[80]="Imach version 0.8d, May 2002, INED-EUROREVES ";
char *alph[]={"a","a","b","c","d","e"}, str[4];	char *alph[]={"a","a","b","c","d","e"}, str[4];


Line 2709 int main(int argc, char *argv[])	Line 3116 int main(int argc, char *argv[])

/-------- arguments in the command line --------/	/-------- arguments in the command line --------/

	/* Log file */
	strcat(filelog, optionfilefiname);
	strcat(filelog,".log"); /* */
	if((ficlog=fopen(filelog,"w"))==NULL) {
	printf("Problem with logfile %s\n",filelog);
	goto end;
	}
	fprintf(ficlog,"Log filename:%s\n",filelog);
	fprintf(ficlog,"\n%s",version);
	fprintf(ficlog,"\nEnter the parameter file name: ");
	fprintf(ficlog,"pathtot=%s, path=%s, optionfile=%s optionfilext=%s optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
	fflush(ficlog);

	/* */
strcpy(fileres,"r");	strcpy(fileres,"r");
strcat(fileres, optionfilefiname);	strcat(fileres, optionfilefiname);
strcat(fileres,".txt"); /* Other files have txt extension */	strcat(fileres,".txt"); /* Other files have txt extension */
Line 2717 int main(int argc, char *argv[])	Line 3138 int main(int argc, char *argv[])

if((ficpar=fopen(optionfile,"r"))==NULL) {	if((ficpar=fopen(optionfile,"r"))==NULL) {
printf("Problem with optionfile %s\n",optionfile);	printf("Problem with optionfile %s\n",optionfile);
	fprintf(ficlog,"Problem with optionfile %s\n",optionfile);
goto end;	goto end;
}	}

strcpy(filereso,"o");	strcpy(filereso,"o");
strcat(filereso,fileres);	strcat(filereso,fileres);
if((ficparo=fopen(filereso,"w"))==NULL) {	if((ficparo=fopen(filereso,"w"))==NULL) {
printf("Problem with Output resultfile: %s\n", filereso);goto end;	printf("Problem with Output resultfile: %s\n", filereso);
	fprintf(ficlog,"Problem with Output resultfile: %s\n", filereso);
	goto end;
}	}

/* Reads comments: lines beginning with '#' */	/* Reads comments: lines beginning with '#' */
Line 2735 int main(int argc, char *argv[])	Line 3159 int main(int argc, char *argv[])
}	}
ungetc(c,ficpar);	ungetc(c,ficpar);

fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);	fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);	printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);	fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
while((c=getc(ficpar))=='#' && c!= EOF){	while((c=getc(ficpar))=='#' && c!= EOF){
Line 2769 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3193 while((c=getc(ficpar))=='#' && c!= EOF){
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);
fprintf(ficparo,"%1d%1d",i1,j1);	fprintf(ficparo,"%1d%1d",i1,j1);
printf("%1d%1d",i,j);	if(mle==1)
	printf("%1d%1d",i,j);
	fprintf(ficlog,"%1d%1d",i,j);
for(k=1; k<=ncovmodel;k++){	for(k=1; k<=ncovmodel;k++){
fscanf(ficpar," %lf",&param[i][j][k]);	fscanf(ficpar," %lf",&param[i][j][k]);
printf(" %lf",param[i][j][k]);	if(mle==1){
	printf(" %lf",param[i][j][k]);
	fprintf(ficlog," %lf",param[i][j][k]);
	}
	else
	fprintf(ficlog," %lf",param[i][j][k]);
fprintf(ficparo," %lf",param[i][j][k]);	fprintf(ficparo," %lf",param[i][j][k]);
}	}
fscanf(ficpar,"\n");	fscanf(ficpar,"\n");
printf("\n");	if(mle==1)
	printf("\n");
	fprintf(ficlog,"\n");
fprintf(ficparo,"\n");	fprintf(ficparo,"\n");
}	}

Line 2824 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3257 while((c=getc(ficpar))=='#' && c!= EOF){
matcov=matrix(1,npar,1,npar);	matcov=matrix(1,npar,1,npar);
for(i=1; i <=npar; i++){	for(i=1; i <=npar; i++){
fscanf(ficpar,"%s",&str);	fscanf(ficpar,"%s",&str);
printf("%s",str);	if(mle==1)
	printf("%s",str);
	fprintf(ficlog,"%s",str);
fprintf(ficparo,"%s",str);	fprintf(ficparo,"%s",str);
for(j=1; j <=i; j++){	for(j=1; j <=i; j++){
fscanf(ficpar," %le",&matcov[i][j]);	fscanf(ficpar," %le",&matcov[i][j]);
printf(" %.5le",matcov[i][j]);	if(mle==1){
	printf(" %.5le",matcov[i][j]);
	fprintf(ficlog," %.5le",matcov[i][j]);
	}
	else
	fprintf(ficlog," %.5le",matcov[i][j]);
fprintf(ficparo," %.5le",matcov[i][j]);	fprintf(ficparo," %.5le",matcov[i][j]);
}	}
fscanf(ficpar,"\n");	fscanf(ficpar,"\n");
printf("\n");	if(mle==1)
	printf("\n");
	fprintf(ficlog,"\n");
fprintf(ficparo,"\n");	fprintf(ficparo,"\n");
}	}
for(i=1; i <=npar; i++)	for(i=1; i <=npar; i++)
for(j=i+1;j<=npar;j++)	for(j=i+1;j<=npar;j++)
matcov[i][j]=matcov[j][i];	matcov[i][j]=matcov[j][i];

printf("\n");	if(mle==1)
	printf("\n");
	fprintf(ficlog,"\n");


/-------- Rewriting paramater file ----------/	/-------- Rewriting paramater file ----------/
Line 2849 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3293 while((c=getc(ficpar))=='#' && c!= EOF){
strcat(rfileres,optionfilext); /* Other files have txt extension */	strcat(rfileres,optionfilext); /* Other files have txt extension */
if((ficres =fopen(rfileres,"w"))==NULL) {	if((ficres =fopen(rfileres,"w"))==NULL) {
printf("Problem writing new parameter file: %s\n", fileres);goto end;	printf("Problem writing new parameter file: %s\n", fileres);goto end;
	fprintf(ficlog,"Problem writing new parameter file: %s\n", fileres);goto end;
}	}
fprintf(ficres,"#%s\n",version);	fprintf(ficres,"#%s\n",version);

/-------- data file ----------/	/-------- data file ----------/
if((fic=fopen(datafile,"r"))==NULL) {	if((fic=fopen(datafile,"r"))==NULL) {
printf("Problem with datafile: %s\n", datafile);goto end;	printf("Problem with datafile: %s\n", datafile);goto end;
	fprintf(ficlog,"Problem with datafile: %s\n", datafile);goto end;
}	}

n= lastobs;	n= lastobs;
Line 2920 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3366 while((c=getc(ficpar))=='#' && c!= EOF){


/* Calculation of the number of parameter from char model*/	/* Calculation of the number of parameter from char model*/
Tvar=ivector(1,15);	Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
Tprod=ivector(1,15);	Tprod=ivector(1,15);
Tvaraff=ivector(1,15);	Tvaraff=ivector(1,15);
Tvard=imatrix(1,15,1,2);	Tvard=imatrix(1,15,1,2);
Line 2936 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3382 while((c=getc(ficpar))=='#' && c!= EOF){
strcpy(modelsav,model);	strcpy(modelsav,model);
if ((strcmp(model,"age")==0) \|\| (strcmp(model,"age*age")==0)){	if ((strcmp(model,"age")==0) \|\| (strcmp(model,"age*age")==0)){
printf("Error. Non available option model=%s ",model);	printf("Error. Non available option model=%s ",model);
	fprintf(ficlog,"Error. Non available option model=%s ",model);
goto end;	goto end;
}	}

for(i=(j+1); i>=1;i--){	for(i=(j+1); i>=1;i--){
cutv(stra,strb,modelsav,'+');	cutv(stra,strb,modelsav,'+'); /* keeps in strb after the last + */
if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav);	if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyze it */
/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/	/* printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
/scanf("%d",i);/	/scanf("%d",i);/
if (strchr(strb,'*')) {	if (strchr(strb,'')) { / Model includes a product */
cutv(strd,strc,strb,'*');	cutv(strd,strc,strb,''); / strdstrc VmVn (if not age)/
if (strcmp(strc,"age")==0) {	if (strcmp(strc,"age")==0) { /* Vnage /
cptcovprod--;	cptcovprod--;
cutv(strb,stre,strd,'V');	cutv(strb,stre,strd,'V');
Tvar[i]=atoi(stre);	Tvar[i]=atoi(stre); /* computes n in Vn and stores in Tvar*/
cptcovage++;	cptcovage++;
Tage[cptcovage]=i;	Tage[cptcovage]=i;
/printf("stre=%s ", stre);/	/printf("stre=%s ", stre);/
}	}
else if (strcmp(strd,"age")==0) {	else if (strcmp(strd,"age")==0) { /* or ageVn /
cptcovprod--;	cptcovprod--;
cutv(strb,stre,strc,'V');	cutv(strb,stre,strc,'V');
Tvar[i]=atoi(stre);	Tvar[i]=atoi(stre);
cptcovage++;	cptcovage++;
Tage[cptcovage]=i;	Tage[cptcovage]=i;
}	}
else {	else { /* Age is not in the model */
cutv(strb,stre,strc,'V');	cutv(strb,stre,strc,'V'); /* strc= Vn, stre is n*/
Tvar[i]=ncovcol+k1;	Tvar[i]=ncovcol+k1;
cutv(strb,strc,strd,'V');	cutv(strb,strc,strd,'V'); /* strd was Vm, strc is m */
Tprod[k1]=i;	Tprod[k1]=i;
Tvard[k1][1]=atoi(strc);	Tvard[k1][1]=atoi(strc); /* m*/
Tvard[k1][2]=atoi(stre);	Tvard[k1][2]=atoi(stre); /* n */
Tvar[cptcovn+k2]=Tvard[k1][1];	Tvar[cptcovn+k2]=Tvard[k1][1];
Tvar[cptcovn+k2+1]=Tvard[k1][2];	Tvar[cptcovn+k2+1]=Tvard[k1][2];
for (k=1; k<=lastobs;k++)	for (k=1; k<=lastobs;k++)
Line 2976 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3423 while((c=getc(ficpar))=='#' && c!= EOF){
k2=k2+2;	k2=k2+2;
}	}
}	}
else {	else { /* no more sum */
/printf("d=%s c=%s b=%s\n", strd,strc,strb);/	/printf("d=%s c=%s b=%s\n", strd,strc,strb);/
/* scanf("%d",i);*/	/* scanf("%d",i);*/
cutv(strd,strc,strb,'V');	cutv(strd,strc,strb,'V');
Line 2985 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3432 while((c=getc(ficpar))=='#' && c!= EOF){
strcpy(modelsav,stra);	strcpy(modelsav,stra);
/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);	/*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
scanf("%d",i);*/	scanf("%d",i);*/
}	} /* end of loop + */
}	} /* end model */

/* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);	/* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
printf("cptcovprod=%d ", cptcovprod);	printf("cptcovprod=%d ", cptcovprod);
	fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
scanf("%d ",i);*/	scanf("%d ",i);*/
fclose(fic);	fclose(fic);

Line 3022 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3470 while((c=getc(ficpar))=='#' && c!= EOF){
else {	else {
if (andc[i]!=9999){	if (andc[i]!=9999){
printf("Warning negative age at death: %d line:%d\n",num[i],i);	printf("Warning negative age at death: %d line:%d\n",num[i],i);
	fprintf(ficlog,"Warning negative age at death: %d line:%d\n",num[i],i);
agev[m][i]=-1;	agev[m][i]=-1;
}	}
}	}
Line 3054 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3503 while((c=getc(ficpar))=='#' && c!= EOF){
for (i=1; i<=imx; i++) {	for (i=1; i<=imx; i++) {
for(m=1; (m<= maxwav); m++){	for(m=1; (m<= maxwav); m++){
if (s[m][i] > (nlstate+ndeath)) {	if (s[m][i] > (nlstate+ndeath)) {
printf("Error: Wrong value in nlstate or ndeath\n");	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);
goto end;	goto end;
}	}
}	}
}	}

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);

free_vector(severity,1,maxwav);	free_vector(severity,1,maxwav);
free_imatrix(outcome,1,maxwav+1,1,n);	free_imatrix(outcome,1,maxwav+1,1,n);
Line 3136 printf("Total number of individuals= %d,	Line 3587 printf("Total number of individuals= %d,
jk=1;	jk=1;
fprintf(ficres,"# Parameters nlstatenlstatencov a121 + b12 age + ...\n");	fprintf(ficres,"# Parameters nlstatenlstatencov a121 + b12 age + ...\n");
printf("# Parameters nlstatenlstatencov a121 + b12 age + ...\n");	printf("# Parameters nlstatenlstatencov a121 + b12 age + ...\n");
	fprintf(ficlog,"# Parameters nlstatenlstatencov a121 + b12 age + ...\n");
for(i=1,jk=1; i <=nlstate; i++){	for(i=1,jk=1; i <=nlstate; i++){
for(k=1; k <=(nlstate+ndeath); k++){	for(k=1; k <=(nlstate+ndeath); k++){
if (k != i)	if (k != i)
{	{
printf("%d%d ",i,k);	printf("%d%d ",i,k);
	fprintf(ficlog,"%d%d ",i,k);
fprintf(ficres,"%1d%1d ",i,k);	fprintf(ficres,"%1d%1d ",i,k);
for(j=1; j <=ncovmodel; j++){	for(j=1; j <=ncovmodel; j++){
printf("%f ",p[jk]);	printf("%f ",p[jk]);
	fprintf(ficlog,"%f ",p[jk]);
fprintf(ficres,"%f ",p[jk]);	fprintf(ficres,"%f ",p[jk]);
jk++;	jk++;
}	}
printf("\n");	printf("\n");
	fprintf(ficlog,"\n");
fprintf(ficres,"\n");	fprintf(ficres,"\n");
}	}
}	}
}	}
if(mle==1){	if(mle==1){
/* Computing hessian and covariance matrix */	/* Computing hessian and covariance matrix */
ftolhess=ftol; /* Usually correct */	ftolhess=ftol; /* Usually correct */
hesscov(matcov, p, npar, delti, ftolhess, func);	hesscov(matcov, p, npar, delti, ftolhess, func);
}	}
fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");	fprintf(ficres,"# Scales (for hessian or gradient estimation)\n");
printf("# Scales (for hessian or gradient estimation)\n");	printf("# Scales (for hessian or gradient estimation)\n");
for(i=1,jk=1; i <=nlstate; i++){	fprintf(ficlog,"# Scales (for hessian or gradient estimation)\n");
for(j=1; j <=nlstate+ndeath; j++){	for(i=1,jk=1; i <=nlstate; i++){
if (j!=i) {	for(j=1; j <=nlstate+ndeath; j++){
fprintf(ficres,"%1d%1d",i,j);	if (j!=i) {
printf("%1d%1d",i,j);	fprintf(ficres,"%1d%1d",i,j);
for(k=1; k<=ncovmodel;k++){	printf("%1d%1d",i,j);
printf(" %.5e",delti[jk]);	fprintf(ficlog,"%1d%1d",i,j);
fprintf(ficres," %.5e",delti[jk]);	for(k=1; k<=ncovmodel;k++){
jk++;	printf(" %.5e",delti[jk]);
}	fprintf(ficlog," %.5e",delti[jk]);
printf("\n");	fprintf(ficres," %.5e",delti[jk]);
fprintf(ficres,"\n");	jk++;
}	}
}	printf("\n");
	fprintf(ficlog,"\n");
	fprintf(ficres,"\n");
	}
}	}
	}
k=1;
fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");	k=1;
printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");	fprintf(ficres,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
for(i=1;i<=npar;i++){	if(mle==1)
/* if (k>nlstate) k=1;	printf("# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
i1=(i-1)/(ncovmodel*nlstate)+1;	fprintf(ficlog,"# Covariance matrix \n# 121 Var(a12)\n# 122 Cov(b12,a12) Var(b12)\n# ...\n# 232 Cov(b23,a12) Cov(b23,b12) ... Var (b23)\n");
fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);	for(i=1;i<=npar;i++){
printf("%s%d%d",alph[k],i1,tab[i]);*/	/* if (k>nlstate) k=1;
fprintf(ficres,"%3d",i);	i1=(i-1)/(ncovmodel*nlstate)+1;
printf("%3d",i);	fprintf(ficres,"%s%d%d",alph[k],i1,tab[i]);
for(j=1; j<=i;j++){	printf("%s%d%d",alph[k],i1,tab[i]);*/
fprintf(ficres," %.5e",matcov[i][j]);	fprintf(ficres,"%3d",i);
printf(" %.5e",matcov[i][j]);	if(mle==1)
}	printf("%3d",i);
fprintf(ficres,"\n");	fprintf(ficlog,"%3d",i);
printf("\n");	for(j=1; j<=i;j++){
k++;	fprintf(ficres," %.5e",matcov[i][j]);
}	if(mle==1)
	printf(" %.5e",matcov[i][j]);
while((c=getc(ficpar))=='#' && c!= EOF){	fprintf(ficlog," %.5e",matcov[i][j]);
ungetc(c,ficpar);	}
fgets(line, MAXLINE, ficpar);	fprintf(ficres,"\n");
puts(line);	if(mle==1)
fputs(line,ficparo);	printf("\n");
}	fprintf(ficlog,"\n");
ungetc(c,ficpar);	k++;
estepm=0;	}
fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
if (estepm==0 \|\| estepm < stepm) estepm=stepm;	while((c=getc(ficpar))=='#' && c!= EOF){
if (fage <= 2) {	ungetc(c,ficpar);
bage = ageminpar;	fgets(line, MAXLINE, ficpar);
fage = agemaxpar;	puts(line);
}	fputs(line,ficparo);
	}
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");	ungetc(c,ficpar);
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);	estepm=0;
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);	fscanf(ficpar,"agemin=%lf agemax=%lf bage=%lf fage=%lf estepm=%d\n",&ageminpar,&agemaxpar, &bage, &fage, &estepm);
	if (estepm==0 \|\| estepm < stepm) estepm=stepm;
while((c=getc(ficpar))=='#' && c!= EOF){	if (fage <= 2) {
ungetc(c,ficpar);	bage = ageminpar;
fgets(line, MAXLINE, ficpar);	fage = agemaxpar;
puts(line);	}
fputs(line,ficparo);
}	fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
ungetc(c,ficpar);	fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
	fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);

	while((c=getc(ficpar))=='#' && c!= EOF){
	ungetc(c,ficpar);
	fgets(line, MAXLINE, ficpar);
	puts(line);
	fputs(line,ficparo);
	}
	ungetc(c,ficpar);

fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);	fscanf(ficpar,"begin-prev-date=%lf/%lf/%lf end-prev-date=%lf/%lf/%lf\n",&jprev1, &mprev1,&anprev1,&jprev2, &mprev2,&anprev2);
fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);	fprintf(ficparo,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);	fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);

while((c=getc(ficpar))=='#' && c!= EOF){	while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);	ungetc(c,ficpar);
fgets(line, MAXLINE, ficpar);	fgets(line, MAXLINE, ficpar);
puts(line);	puts(line);
fputs(line,ficparo);	fputs(line,ficparo);
}	}
ungetc(c,ficpar);	ungetc(c,ficpar);


dateprev1=anprev1+mprev1/12.+jprev1/365.;	dateprev1=anprev1+mprev1/12.+jprev1/365.;
Line 3269 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3736 while((c=getc(ficpar))=='#' && c!= EOF){
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);	freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax,mint,anint,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);

/------------ gnuplot -------------/	/------------ gnuplot -------------/
printinggnuplot(fileres,optionfilefiname,optionfile,optionfilegnuplot, ageminpar,agemaxpar,fage, pathc,p);	strcpy(optionfilegnuplot,optionfilefiname);
	strcat(optionfilegnuplot,".gp");
	if((ficgp=fopen(optionfilegnuplot,"w"))==NULL) {
	printf("Problem with file %s",optionfilegnuplot);
	}
	fclose(ficgp);
	printinggnuplot(fileres, ageminpar,agemaxpar,fage, pathc,p);
	/--------- index.htm --------/

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

	fprintf(fichtm,"<body> <font size=\"2\">%s </font> <hr size=\"2\" color=\"#EC5E5E\"> \n
	Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=%s<br>\n
	\n
	Total number of observations=%d <br>\n
	Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n
	<hr size=\"2\" color=\"#EC5E5E\">
	<ul><li><h4>Parameter files</h4>\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
	- Gnuplot file name: <a href=\"%s\">%s</a></ul>\n",version,title,datafile,firstpass,lastpass,stepm, weightopt,model,imx,jmin,jmax,jmean,fileres,fileres,filelog,filelog,optionfilegnuplot,optionfilegnuplot);
	fclose(fichtm);

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

/------------ free_vector -------------/	/------------ free_vector -------------/
chdir(path);	chdir(path);
Line 3283 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3777 while((c=getc(ficpar))=='#' && c!= EOF){
fclose(ficparo);	fclose(ficparo);
fclose(ficres);	fclose(ficres);

/--------- index.htm --------/

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


/--------------- Prevalence limit --------------/	/--------------- Prevalence limit --------------/

strcpy(filerespl,"pl");	strcpy(filerespl,"pl");
strcat(filerespl,fileres);	strcat(filerespl,fileres);
if((ficrespl=fopen(filerespl,"w"))==NULL) {	if((ficrespl=fopen(filerespl,"w"))==NULL) {
printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;	printf("Problem with Prev limit resultfile: %s\n", filerespl);goto end;
	fprintf(ficlog,"Problem with Prev limit resultfile: %s\n", filerespl);goto end;
}	}
printf("Computing prevalence limit: result on file '%s' \n", filerespl);	printf("Computing prevalence limit: result on file '%s' \n", filerespl);
	fprintf(ficlog,"Computing prevalence limit: result on file '%s' \n", filerespl);
fprintf(ficrespl,"#Prevalence limit\n");	fprintf(ficrespl,"#Prevalence limit\n");
fprintf(ficrespl,"#Age ");	fprintf(ficrespl,"#Age ");
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);	for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
Line 3319 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3811 while((c=getc(ficpar))=='#' && c!= EOF){
k=k+1;	k=k+1;
/printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);/	/printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);/
fprintf(ficrespl,"\n#******");	fprintf(ficrespl,"\n#******");
for(j=1;j<=cptcoveff;j++)	printf("\n#******");
	fprintf(ficlog,"\n#******");
	for(j=1;j<=cptcoveff;j++) {
fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);	fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
	printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
	fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
	}
fprintf(ficrespl,"******\n");	fprintf(ficrespl,"******\n");
	printf("******\n");
	fprintf(ficlog,"******\n");

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);
Line 3339 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3838 while((c=getc(ficpar))=='#' && c!= EOF){
strcpy(filerespij,"pij"); strcat(filerespij,fileres);	strcpy(filerespij,"pij"); strcat(filerespij,fileres);
if((ficrespij=fopen(filerespij,"w"))==NULL) {	if((ficrespij=fopen(filerespij,"w"))==NULL) {
printf("Problem with Pij resultfile: %s\n", filerespij);goto end;	printf("Problem with Pij resultfile: %s\n", filerespij);goto end;
	fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij);goto end;
}	}
printf("Computing pij: result on file '%s' \n", filerespij);	printf("Computing pij: result on file '%s' \n", filerespij);
	fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);

stepsize=(int) (stepm+YEARM-1)/YEARM;	stepsize=(int) (stepm+YEARM-1)/YEARM;
/if (stepm<=24) stepsize=2;/	/if (stepm<=24) stepsize=2;/
Line 3348 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3849 while((c=getc(ficpar))=='#' && c!= EOF){
agelim=AGESUP;	agelim=AGESUP;
hstepm=stepsizeYEARM; / Every year of age */	hstepm=stepsizeYEARM; / Every year of age */
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */	hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */

	/* hstepm=1; aff par mois*/

k=0;	k=0;
for(cptcov=1;cptcov<=i1;cptcov++){	for(cptcov=1;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){	for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
Line 3361 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3864 while((c=getc(ficpar))=='#' && c!= EOF){
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */	for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
nhstepm=(int) rint((agelim-agedeb)YEARM/stepm); / Typically 20 years = 2012/6=40 /	nhstepm=(int) rint((agelim-agedeb)YEARM/stepm); / Typically 20 years = 2012/6=40 /
nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */	nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */

	/* nhstepm=nhstepmYEARM; aff par mois/

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);
Line 3370 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3876 while((c=getc(ficpar))=='#' && c!= EOF){
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 %.0f %.0f",k,agedeb, agedeb+ hhstepm/YEARMstepm );	fprintf(ficrespij,"%d %f %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 3382 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3888 while((c=getc(ficpar))=='#' && c!= EOF){
}	}
}	}

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

fclose(ficrespij);	fclose(ficrespij);

Line 3395 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3901 while((c=getc(ficpar))=='#' && c!= EOF){
else{	else{
erreur=108;	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);	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 3404 while((c=getc(ficpar))=='#' && c!= EOF){	Line 3911 while((c=getc(ficpar))=='#' && c!= EOF){
strcat(filerest,fileres);	strcat(filerest,fileres);
if((ficrest=fopen(filerest,"w"))==NULL) {	if((ficrest=fopen(filerest,"w"))==NULL) {
printf("Problem with total LE resultfile: %s\n", filerest);goto end;	printf("Problem with total LE resultfile: %s\n", filerest);goto end;
	fprintf(ficlog,"Problem with total LE resultfile: %s\n", filerest);goto end;
}	}
printf("Computing Total LEs with variances: file '%s' \n", filerest);	printf("Computing Total LEs with variances: file '%s' \n", filerest);
	fprintf(ficlog,"Computing Total LEs with variances: file '%s' \n", filerest);


strcpy(filerese,"e");	strcpy(filerese,"e");
strcat(filerese,fileres);	strcat(filerese,fileres);
if((ficreseij=fopen(filerese,"w"))==NULL) {	if((ficreseij=fopen(filerese,"w"))==NULL) {
printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);	printf("Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
	fprintf(ficlog,"Problem with Health Exp. resultfile: %s\n", filerese); exit(0);
}	}
printf("Computing Health Expectancies: result on file '%s' \n", filerese);	printf("Computing Health Expectancies: result on file '%s' \n", filerese);
	fprintf(ficlog,"Computing Health Expectancies: result on file '%s' \n", filerese);

strcpy(fileresv,"v");	strcpy(fileresv,"v");
strcat(fileresv,fileres);	strcat(fileresv,fileres);
if((ficresvij=fopen(fileresv,"w"))==NULL) {	if((ficresvij=fopen(fileresv,"w"))==NULL) {
printf("Problem with variance resultfile: %s\n", fileresv);exit(0);	printf("Problem with variance resultfile: %s\n", fileresv);exit(0);
	fprintf(ficlog,"Problem with variance resultfile: %s\n", fileresv);exit(0);
}	}
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);
calagedate=-1;	calagedate=-1;
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, calagedate);

k=0;	k=0;
for(cptcov=1;cptcov<=i1;cptcov++){	for(cptcov=1;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){	for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;	k=k+1;
fprintf(ficrest,"\n#****** ");	fprintf(ficrest,"\n#****** ");
for(j=1;j<=cptcoveff;j++)	for(j=1;j<=cptcoveff;j++)
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);	fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
Line 3449 prevalence(ageminpar, agemax, s, agev, n	Line 3962 prevalence(ageminpar, agemax, s, agev, n

vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);	vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;	oldm=oldms;savm=savms;
varevsij(fileres, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm);	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,0);
	if(popbased==1){
	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,popbased);
	}


fprintf(ficrest,"#Total LEs with variances: e.. (std) ");	fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
Line 3538 free_matrix(mint,1,maxwav,1,n);	Line 4053 free_matrix(mint,1,maxwav,1,n);
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(erreur >0)	fprintf(fichtm,"\n</body>");
	fclose(fichtm);
	fclose(ficgp);


	if(erreur >0){
printf("End of Imach with error or warning %d\n",erreur);	printf("End of Imach with error or warning %d\n",erreur);
else printf("End of Imach\n");	fprintf(ficlog,"End of Imach with error or warning %d\n",erreur);
	}else{
	printf("End of Imach\n");
	fprintf(ficlog,"End of Imach\n");
	}
	printf("See log file on %s\n",filelog);
	fclose(ficlog);
/* gettimeofday(&end_time, (struct timezone)0);/ /* after time */	/* gettimeofday(&end_time, (struct timezone)0);/ /* after time */

/* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/	/* printf("Total time was %d Sec. %d uSec.\n", end_time.tv_sec -start_time.tv_sec, end_time.tv_usec -start_time.tv_usec);*/

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

Removed from v.1.43
changed lines
	Added in v.1.50