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.238 and 1.247

version 1.238, 2016/08/26 14:23:35	version 1.247, 2016/09/02 11:11:21
Line 1	Line 1
/* $Id$	/* $Id$
$State$	$State$
$Log$	$Log$
	Revision 1.247 2016/09/02 11:11:21 brouard
	* empty log message *

	Revision 1.246 2016/09/02 08:49:22 brouard
	* empty log message *

	Revision 1.245 2016/09/02 07:25:01 brouard
	* empty log message *

	Revision 1.244 2016/09/02 07:17:34 brouard
	* empty log message *

	Revision 1.243 2016/09/02 06:45:35 brouard
	* empty log message *

	Revision 1.242 2016/08/30 15:01:20 brouard
	Summary: Fixing a lots

	Revision 1.241 2016/08/29 17:17:25 brouard
	Summary: gnuplot problem in Back projection to fix

	Revision 1.240 2016/08/29 07:53:18 brouard
	Summary: Better

	Revision 1.239 2016/08/26 15:51:03 brouard
	Summary: Improvement in Powell output in order to copy and paste

	Author:

Revision 1.238 2016/08/26 14:23:35 brouard	Revision 1.238 2016/08/26 14:23:35 brouard
Summary: Starting tests of 0.99	Summary: Starting tests of 0.99

Line 1032 double dval;	Line 1061 double dval;
#define FTOL 1.0e-10	#define FTOL 1.0e-10

#define NRANSI	#define NRANSI
#define ITMAX 200	#define ITMAX 200
	#define ITPOWMAX 20 /* This is now multiplied by the number of parameters */

#define TOL 2.0e-4	#define TOL 2.0e-4

Line 2061 void powell(double p[], double **xi, int	Line 2091 void powell(double p[], double **xi, int
void linmin(double p[], double xi[], int n, double *fret,	void linmin(double p[], double xi[], int n, double *fret,
double (*func)(double []));	double (*func)(double []));
#else	#else
void linmin(double p[], double xi[], int n, double *fret,	void linmin(double p[], double xi[], int n, double *fret,
double (func)(double []),int flat);	double (func)(double []),int flat);
#endif	#endif
int i,ibig,j;	int i,ibig,j,jk,k;
double del,t,pt,ptt,*xit;	double del,t,pt,ptt,*xit;
double directest;	double directest;
double fp,fptt;	double fp,fptt;
Line 2096 void powell(double p[], double **xi, int	Line 2126 void powell(double p[], double **xi, int
fprintf(ficlog,"\nPowell iter=%d -2LL=%.12f %ld sec. %ld sec.",iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);	fprintf(ficlog,"\nPowell iter=%d -2LL=%.12f %ld sec. %ld sec.",iter,*fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
/* fprintf(ficrespow,"%d %.12f %ld",iter,fret,curr_time.tm_sec-start_time.tm_sec); */	/* fprintf(ficrespow,"%d %.12f %ld",iter,fret,curr_time.tm_sec-start_time.tm_sec); */
for (i=1;i<=n;i++) {	for (i=1;i<=n;i++) {
printf(" %d %.12f",i, p[i]);
fprintf(ficlog," %d %.12lf",i, p[i]);
fprintf(ficrespow," %.12lf", p[i]);	fprintf(ficrespow," %.12lf", p[i]);
}	}
	fprintf(ficrespow,"\n");fflush(ficrespow);
	printf("\n#model= 1 + age ");
	fprintf(ficlog,"\n#model= 1 + age ");
	if(nagesqr==1){
	printf(" + age*age ");
	fprintf(ficlog," + age*age ");
	}
	for(j=1;j <=ncovmodel-2;j++){
	if(Typevar[j]==0) {
	printf(" + V%d ",Tvar[j]);
	fprintf(ficlog," + V%d ",Tvar[j]);
	}else if(Typevar[j]==1) {
	printf(" + V%d*age ",Tvar[j]);
	fprintf(ficlog," + V%d*age ",Tvar[j]);
	}else if(Typevar[j]==2) {
	printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
	fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
	}
	}
printf("\n");	printf("\n");
	/* printf("12 47.0114589 0.0154322 33.2424412 0.3279905 2.3731903 */
	/* 13 -21.5392400 0.1118147 1.2680506 1.2973408 -1.0663662 */
fprintf(ficlog,"\n");	fprintf(ficlog,"\n");
fprintf(ficrespow,"\n");fflush(ficrespow);	for(i=1,jk=1; i <=nlstate; i++){
if(*iter <=3){	for(k=1; k <=(nlstate+ndeath); k++){
	if (k != i) {
	printf("%d%d ",i,k);
	fprintf(ficlog,"%d%d ",i,k);
	for(j=1; j <=ncovmodel; j++){
	printf("%12.7f ",p[jk]);
	fprintf(ficlog,"%12.7f ",p[jk]);
	jk++;
	}
	printf("\n");
	fprintf(ficlog,"\n");
	}
	}
	}
	if(iter <=3 && iter >1){
tml = *localtime(&rcurr_time);	tml = *localtime(&rcurr_time);
strcpy(strcurr,asctime(&tml));	strcpy(strcurr,asctime(&tml));
rforecast_time=rcurr_time;	rforecast_time=rcurr_time;
itmp = strlen(strcurr);	itmp = strlen(strcurr);
if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */	if(strcurr[itmp-1]=='\n') /* Windows outputs with a new line */
strcurr[itmp-1]='\0';	strcurr[itmp-1]='\0';
printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);	printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);	fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
for(niterf=10;niterf<=30;niterf+=10){	for(niterf=10;niterf<=30;niterf+=10){
rforecast_time=rcurr_time+(niterf-iter)(rcurr_time-rlast_time);	rforecast_time=rcurr_time+(niterf-iter)(rcurr_time-rlast_time);
forecast_time = *localtime(&rforecast_time);	forecast_time = *localtime(&rforecast_time);
strcpy(strfor,asctime(&forecast_time));	strcpy(strfor,asctime(&forecast_time));
itmp = strlen(strfor);	itmp = strlen(strfor);
if(strfor[itmp-1]=='\n')	if(strfor[itmp-1]=='\n')
strfor[itmp-1]='\0';	strfor[itmp-1]='\0';
printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);	printf(" - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);	fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
}	}
}	}
for (i=1;i<=n;i++) { /* For each direction i */	for (i=1;i<=n;i++) { /* For each direction i */
Line 2174 void powell(double p[], double **xi, int	Line 2237 void powell(double p[], double **xi, int
/* printf("\n"); */	/* printf("\n"); */
/* fprintf(ficlog,"\n"); */	/* fprintf(ficlog,"\n"); */
}	}
if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret))) { /* Did we reach enough precision? */	/* if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret))) { /\* Did we reach enough precision? \/ /
	if (2.0fabs(fp-(fret)) <= ftol) { /* Did we reach enough precision? */
/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */	/* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
/* By adding ageage in a model, the new -2LL should be lower and the difference follows a /	/* By adding ageage in a model, the new -2LL should be lower and the difference follows a /
/* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */	/* a chisquare statistics with 1 degree. To be significant at the 95% level, it should have */
Line 2220 void powell(double p[], double **xi, int	Line 2284 void powell(double p[], double **xi, int
free_vector(pt,1,n);	free_vector(pt,1,n);
return;	return;
} /* enough precision */	} /* enough precision */
if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");	if (iter == ITMAXn) nrerror("powell exceeding maximum iterations.");
for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */	for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
ptt[j]=2.0*p[j]-pt[j];	ptt[j]=2.0*p[j]-pt[j];
xit[j]=p[j]-pt[j];	xit[j]=p[j]-pt[j];
Line 2500 Earliest age to start was %d-%d=%d, ncvl	Line 2564 Earliest age to start was %d-%d=%d, ncvl

/* double bprevalim(double bprlim, double *prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double oldm, double savm, double dnewm, double doldm, double dsavm, double ftolpl, int ncvyear, int ij) /	/* double bprevalim(double bprlim, double *prevacurrent, int nlstate, double x[], double age, double ageminpar, double agemaxpar, double oldm, double savm, double dnewm, double doldm, double dsavm, double ftolpl, int ncvyear, int ij) /
/* double bprevalim(double bprlim, double *prevacurrent, int nlstate, double x[], double age, double oldm, double savm, double dnewm, double doldm, double dsavm, double ftolpl, int ncvyear, int ij) /	/* double bprevalim(double bprlim, double *prevacurrent, int nlstate, double x[], double age, double oldm, double savm, double dnewm, double doldm, double dsavm, double ftolpl, int ncvyear, int ij) /
double bprevalim(double bprlim, double **prevacurrent, int nlstate, double x[], double age, double ftolpl, int ncvyear, int ij)	double bprevalim(double bprlim, double **prevacurrent, int nlstate, double x[], double age, double ftolpl, int ncvyear, int ij, int nres)
{	{
/* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit	/* Computes the prevalence limit in each live state at age x and covariate ij by left multiplying the unit
matrix by transitions matrix until convergence is reached with precision ftolpl */	matrix by transitions matrix until convergence is reached with precision ftolpl */
Line 2523 Earliest age to start was %d-%d=%d, ncvl	Line 2587 Earliest age to start was %d-%d=%d, ncvl
/* If we start from prlim again, prlim tends to a constant matrix */	/* If we start from prlim again, prlim tends to a constant matrix */

int i, ii,j,k;	int i, ii,j,k;
	int first=0;
double min, max, *meandiff, maxmax,sumnew=0.;	double min, max, *meandiff, maxmax,sumnew=0.;
/* double *matprod2(); / /* test */	/* double *matprod2(); / /* test */
double out, cov[NCOVMAX+1], bmij();	double out, cov[NCOVMAX+1], bmij();
Line 2559 Earliest age to start was %d-%d=%d, ncvl	Line 2624 Earliest age to start was %d-%d=%d, ncvl
cov[2]=agefin;	cov[2]=agefin;
if(nagesqr==1)	if(nagesqr==1)
cov[3]= agefin*agefin;;	cov[3]= agefin*agefin;;
for (k=1; k<=cptcovn;k++) {	for (k=1; k<=nsd;k++) { /* For single dummy covariates only */
/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; */	/* Here comes the value of the covariate 'ij' after renumbering k with single dummy covariates */
cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)];	cov[2+nagesqr+TvarsDind[k]]=nbcode[TvarsD[k]][codtabm(ij,k)];
/* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); */	/* printf("bprevalim Dummy combi=%d k=%d TvarsD[%d]=V%d TvarsDind[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, TvarsD[k],k,TvarsDind[k],nbcode[TvarsD[k]][codtabm(ij,k)],cov[2+nagesqr+TvarsDind[k]], ij, k, codtabm(ij,k)); */
	}
	/* for (k=1; k<=cptcovn;k++) { */
	/* /\* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,Tvar[k])]; \/ /
	/* cov[2+nagesqr+k]=nbcode[Tvar[k]][codtabm(ij,k)]; */
	/* /\* printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtabm(%d,Tvar[%d])=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtabm(ij,Tvar[k])],cov[2+k], ij, k, codtabm(ij,Tvar[k])]); \/ /
	/* } */
	for (k=1; k<=nsq;k++) { /* For single varying covariates only */
	/* Here comes the value of quantitative after renumbering k with single quantitative covariates */
	cov[2+nagesqr+TvarsQind[k]]=Tqresult[nres][k];
	/* printf("prevalim Quantitative k=%d TvarsQind[%d]=%d, TvarsQ[%d]=V%d,Tqresult[%d][%d]=%f\n",k,k,TvarsQind[k],k,TvarsQ[k],nres,k,Tqresult[nres][k]); */
	}
	/* for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]cov[2]; /
	/* for (k=1; k<=cptcovprod;k++) /\* Useless \/ /
	/* /\* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; \/ /
	/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)]; */
	for (k=1; k<=cptcovage;k++){ /* For product with age */
	if(Dummy[Tvar[Tage[k]]]){
	cov[2+nagesqr+Tage[k]]=nbcode[Tvar[Tage[k]]][codtabm(ij,k)]*cov[2];
	} else{
	cov[2+nagesqr+Tage[k]]=Tqresult[nres][k];
	}
	/* printf("prevalim Age combi=%d k=%d Tage[%d]=V%d Tqresult[%d][%d]=%f\n",ij,k,k,Tage[k],nres,k,Tqresult[nres][k]); */
	}
	for (k=1; k<=cptcovprod;k++){ /* For product without age */
	/* printf("prevalim Prod ij=%d k=%d Tprod[%d]=%d Tvard[%d][1]=V%d, Tvard[%d][2]=V%d\n",ij,k,k,Tprod[k], k,Tvard[k][1], k,Tvard[k][2]); */
	if(Dummy[Tvard[k][1]==0]){
	if(Dummy[Tvard[k][2]==0]){
	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];
	}else{
	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * Tqresult[nres][k];
	}
	}else{
	if(Dummy[Tvard[k][2]==0]){
	cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][2]][codtabm(ij,k)] * Tqinvresult[nres][Tvard[k][1]];
	}else{
	cov[2+nagesqr+Tprod[k]]=Tqinvresult[nres][Tvard[k][1]]* Tqinvresult[nres][Tvard[k][2]];
	}
	}
}	}
for (k=1; k<=cptcovage;k++) cov[2+nagesqr+Tage[k]]=nbcode[Tvar[k]][codtabm(ij,k)]*cov[2];
for (k=1; k<=cptcovprod;k++) /* Useless */
/* cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,Tvard[k][1])] * nbcode[Tvard[k][2]][codtabm(ij,Tvard[k][2])]; */
cov[2+nagesqr+Tprod[k]]=nbcode[Tvard[k][1]][codtabm(ij,k)] * nbcode[Tvard[k][2]][codtabm(ij,k)];

/printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);/	/printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);/
/printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);/	/printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);/
Line 2614 Earliest age to start was %d-%d=%d, ncvl	Line 2713 Earliest age to start was %d-%d=%d, ncvl
}	}
} /* age loop */	} /* age loop */
/* After some age loop it doesn't converge */	/* After some age loop it doesn't converge */
printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\	if(first){
	first=1;
	printf("Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. Others in log file only...\n\
	Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
	}
	fprintf(ficlog,"Warning: the back stable prevalence at age %d did not converge with the required precision (%g > ftolpl=%g) within %.0f years. Try to lower 'ftolpl'. \n\
Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);	Oldest age to start was %d-%d=%d, ncvloop=%d, ncvyear=%d\n", (int)age, maxmax, ftolpl, delaymax, (int)age, (int)delaymax, (int)agefin, ncvloop, *ncvyear);
/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */	/* Try to lower 'ftol', for example from 1.e-8 to 6.e-9.\n", ftolpl, (int)age, (int)delaymax, (int)agefin, ncvloop, (int)age-(int)agefin); */
free_vector(min,1,nlstate);	free_vector(min,1,nlstate);
Line 2747 double pmij(double ps, double *cov,	Line 2851 double pmij(double ps, double *cov,
/* }else */	/* }else */
doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);	doldm[ii][j]=(ii==j ? 1./sumnew : 0.0);
}else{	}else{
printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin);	;
	/* printf("ii=%d, i=%d, doldm=%lf dsavm=%lf, probs=%lf, sumnew=%lf,agefin=%d\n",ii,j,doldm[ii][j],dsavm[ii][j],prevacurrent[(int)agefin][ii][ij],sumnew, (int)agefin); */
}	}
} /End ii /	} /End ii /
} /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */	} /* End j, At the end doldm is diag[1/(w_1p1i+w_2 p2i)] */
Line 3130 double func( double *x)	Line 3235 double func( double *x)
Then computes with function pmij which return a matrix p[i][j] giving the elementary probability	Then computes with function pmij which return a matrix p[i][j] giving the elementary probability
to be observed in j being in i according to the model.	to be observed in j being in i according to the model.
*/	*/
ioffset=2+nagesqr+cptcovage;	ioffset=2+nagesqr ;
/* Fixed */	/* Fixed */
for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */	for (k=1; k<=ncovf;k++){ /* Simple and product fixed covariates without age* products */
cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1, only V1 is fixed (k=6)*/	cov[ioffset+TvarFind[k]]=covar[Tvar[TvarFind[k]]][i];/* V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1, only V1 is fixed (k=6)*/
Line 3149 double func( double *x)	Line 3254 double func( double *x)
*/	*/
for(mi=1; mi<= wav[i]-1; mi++){	for(mi=1; mi<= wav[i]-1; mi++){
for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/	for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i];	/* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */
	cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
}	}
for (ii=1;ii<=nlstate+ndeath;ii++)	for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){	for (j=1;j<=nlstate+ndeath;j++){
Line 3163 double func( double *x)	Line 3269 double func( double *x)
if(nagesqr==1)	if(nagesqr==1)
cov[3]= agexactagexact; / Should be changed here */	cov[3]= agexactagexact; / Should be changed here */
for (kk=1; kk<=cptcovage;kk++) {	for (kk=1; kk<=cptcovage;kk++) {
	if(!FixedV[Tvar[Tage[kk]]])
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]agexact; / Tage[kk] gives the data-covariate associated with age */	cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]agexact; / Tage[kk] gives the data-covariate associated with age */
	else
	cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact;
}	}
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));	1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
Line 3445 double funcone( double *x)	Line 3554 double funcone( double *x)
for(k=1; k<=nlstate; k++) ll[k]=0.;	for(k=1; k<=nlstate; k++) ll[k]=0.;
ioffset=0;	ioffset=0;
for (i=1,ipmx=0, sw=0.; i<=imx; i++){	for (i=1,ipmx=0, sw=0.; i<=imx; i++){
ioffset=2+nagesqr+cptcovage;	/* ioffset=2+nagesqr+cptcovage; */
	ioffset=2+nagesqr;
/* Fixed */	/* Fixed */
/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */	/* for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i]; */
/* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products \/ /	/* for (k=1; k<=ncoveff;k++){ /\* Simple and product fixed Dummy covariates without age* products \/ /
Line 3472 double funcone( double *x)	Line 3582 double funcone( double *x)
for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */	for(mi=1; mi<= wav[i]-1; mi++){ /* Varying with waves */
/* Wave varying (but not age varying) */	/* Wave varying (but not age varying) */
for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/	for(k=1; k <= ncovv ; k++){ /* Varying covariates (single and product but no age )*/
cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i];	/* cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]][i]; */
}	cov[ioffset+TvarVind[k]]=cotvar[mw[mi][i]][Tvar[TvarVind[k]]-ncovcol-nqv][i];
	}
/* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)\/ /	/* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates (single??)\/ /
/* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff \/ /	/* iv= Tvar[Tmodelind[ioffset-2-nagesqr-cptcovage+itv]]-ncovcol-nqv; /\* Counting the # varying covariate from 1 to ntveff \/ /
/* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */	/* cov[ioffset+iv]=cotvar[mw[mi][i]][iv][i]; */
/* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model \/ /	/* k=ioffset-2-nagesqr-cptcovage+itv; /\* position in simple model \/ /
/* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */	/* cov[ioffset+itv]=cotvar[mw[mi][i]][TmodelInvind[itv]][i]; */
/* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */	/* printf(" i=%d,mi=%d,itv=%d,TmodelInvind[itv]=%d,cotvar[mw[mi][i]][TmodelInvind[itv]][i]=%f\n", i, mi, itv, TmodelInvind[itv],cotvar[mw[mi][i]][TmodelInvind[itv]][i]); */
/* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates \/ /	/* for(iqtv=1; iqtv <= nqtveff; iqtv++){ /\* Varying quantitatives covariates \/ /
/* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff \/ /	/* iv=TmodelInvQind[iqtv]; /\* Counting the # varying covariate from 1 to ntveff \/ /
/* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); \/ /	/* /\* printf(" i=%d,mi=%d,iqtv=%d,TmodelInvQind[iqtv]=%d,cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]=%f\n", i, mi, iqtv, TmodelInvQind[iqtv],cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]); \/ /
/* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */	/* cov[ioffset+ntveff+iqtv]=cotqvar[mw[mi][i]][TmodelInvQind[iqtv]][i]; */
/* } */	/* } */
for (ii=1;ii<=nlstate+ndeath;ii++)	for (ii=1;ii<=nlstate+ndeath;ii++)
for (j=1;j<=nlstate+ndeath;j++){	for (j=1;j<=nlstate+ndeath;j++){
oldm[ii][j]=(ii==j ? 1.0 : 0.0);	oldm[ii][j]=(ii==j ? 1.0 : 0.0);
savm[ii][j]=(ii==j ? 1.0 : 0.0);	savm[ii][j]=(ii==j ? 1.0 : 0.0);
}	}

agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */	agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
ageend=agev[mw[mi][i]][i] + (dh[mi][i])stepm/YEARM; / Age at end of effective wave and at the end of transition */	ageend=agev[mw[mi][i]][i] + (dh[mi][i])stepm/YEARM; / Age at end of effective wave and at the end of transition */
for(d=0; d<dh[mi][i]; d++){ /* Delay between two effective waves */	for(d=0; d<dh[mi][i]; d++){ /* Delay between two effective waves */
/*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]	/* for(d=0; d<=0; d++){ /\* Delay between two effective waves Only one matrix to speed up\/ /
and mw[mi+1][i]. dh depends on stepm.*/	/*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
newm=savm;	and mw[mi+1][i]. dh depends on stepm.*/
agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;	newm=savm;
cov[2]=agexact;	agexact=agev[mw[mi][i]][i]+dstepm/YEARM; / Here d is needed */
if(nagesqr==1)	cov[2]=agexact;
cov[3]= agexact*agexact;	if(nagesqr==1)
for (kk=1; kk<=cptcovage;kk++) {	cov[3]= agexact*agexact;
cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;	for (kk=1; kk<=cptcovage;kk++) {
}	if(!FixedV[Tvar[Tage[kk]]])
/* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */	cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */	else
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,	cov[Tage[kk]+2+nagesqr]=cotvar[mw[mi][i]][Tvar[Tage[kk]]-ncovcol-nqv][i]*agexact;
1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));	}
/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */	/* printf("i=%d,mi=%d,d=%d,mw[mi][i]=%d\n",i, mi,d,mw[mi][i]); */
/* 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */	/* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
savm=oldm;	out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
oldm=newm;	1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
	/* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */
	/* 1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */
	savm=oldm;
	oldm=newm;
} /* end mult */	} /* end mult */

s1=s[mw[mi][i]][i];	s1=s[mw[mi][i]][i];
Line 3525 double funcone( double *x)	Line 3640 double funcone( double *x)
* is higher than the multiple of stepm and negative otherwise.	* is higher than the multiple of stepm and negative otherwise.
*/	*/
if( s2 > nlstate && (mle <5) ){ /* Jackson */	if( s2 > nlstate && (mle <5) ){ /* Jackson */
lli=log(out[s1][s2] - savm[s1][s2]);	lli=log(out[s1][s2] - savm[s1][s2]);
} else if ( s2==-1 ) { /* alive */	} else if ( s2==-1 ) { /* alive */
for (j=1,survp=0. ; j<=nlstate; j++)	for (j=1,survp=0. ; j<=nlstate; j++)
survp += (1.+bbh)out[s1][j]- bbhsavm[s1][j];	survp += (1.+bbh)out[s1][j]- bbhsavm[s1][j];
lli= log(survp);	lli= log(survp);
}else if (mle==1){	}else if (mle==1){
lli= log((1.+bbh)out[s1][s2]- bbhsavm[s1][s2]); /* linear interpolation */	lli= log((1.+bbh)out[s1][s2]- bbhsavm[s1][s2]); /* linear interpolation */
} else if(mle==2){	} else if(mle==2){
lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)out[s1][s2]- bbhsavm[s1][s2]):log((1.+bbh)out[s1][s2])); / linear interpolation */	lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)out[s1][s2]- bbhsavm[s1][s2]):log((1.+bbh)out[s1][s2])); / linear interpolation */
} else if(mle==3){ /* exponential inter-extrapolation */	} else if(mle==3){ /* exponential inter-extrapolation */
lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.+bbh)out[s1][s2])); / exponential inter-extrapolation */	lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.+bbh)out[s1][s2])); / exponential inter-extrapolation */
} else if (mle==4){ /* mle=4 no inter-extrapolation */	} else if (mle==4){ /* mle=4 no inter-extrapolation */
lli=log(out[s1][s2]); /* Original formula */	lli=log(out[s1][s2]); /* Original formula */
} else{ /* mle=0 back to 1 */	} else{ /* mle=0 back to 1 */
lli= log((1.+bbh)out[s1][s2]- bbhsavm[s1][s2]); /* linear interpolation */	lli= log((1.+bbh)out[s1][s2]- bbhsavm[s1][s2]); /* linear interpolation */
/lli=log(out[s1][s2]); / /* Original formula */	/lli=log(out[s1][s2]); / /* Original formula */
} /* End of if */	} /* End of if */
ipmx +=1;	ipmx +=1;
sw += weight[i];	sw += weight[i];
ll[s[mw[mi][i]][i]] += 2weight[i]lli;	ll[s[mw[mi][i]][i]] += 2weight[i]lli;
/printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); /	/printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); /
if(globpr){	if(globpr){
fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\	fprintf(ficresilk,"%09ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %15.6f %8.4f %8.3f\
%11.6f %11.6f %11.6f ", \	%11.6f %11.6f %11.6f ", \
num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,	num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
2weight[i]lli,out[s1][s2],savm[s1][s2]);	2weight[i]lli,out[s1][s2],savm[s1][s2]);
for(k=1,llt=0.,l=0.; k<=nlstate; k++){	for(k=1,llt=0.,l=0.; k<=nlstate; k++){
llt +=ll[k]*gipmx/gsw;	llt +=ll[k]*gipmx/gsw;
fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);	fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
}	}
fprintf(ficresilk," %10.6f\n", -llt);	fprintf(ficresilk," %10.6f\n", -llt);
}	}
} /* end of wave */	} /* end of wave */
} /* end of individual */	} /* end of individual */
Line 3905 double hessij( double x[], double **hess	Line 4020 double hessij( double x[], double **hess
kmax=kmax+10;	kmax=kmax+10;
}	}
if(kmax >=10 \|\| firstime ==1){	if(kmax >=10 \|\| firstime ==1){
printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);	printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you may increase ftol=%.2e\n",thetai,thetaj, ftol);	fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol);
printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e deltik=%.12e deltjk=%.12e, xi-dek=%.12e xj-dek=%.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 deltik=%.12e deltjk=%.12e, xi-dek=%.12e xj-dek=%.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 deltik=%.12e deltjk=%.12e, xi-dek=%.12e xj-dek=%.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 deltik=%.12e deltjk=%.12e, xi-dek=%.12e xj-dek=%.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);
}	}
Line 4094 void freqsummary(char fileres[], int ia	Line 4209 void freqsummary(char fileres[], int ia
fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);	fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
exit(0);	exit(0);
}	}

strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));	strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {	if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));	printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
Line 4104 void freqsummary(char fileres[], int ia	Line 4219 void freqsummary(char fileres[], int ia
}	}
else{	else{
fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \	fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\	<hr size=\"2\" color=\"#EC5E5E\"> \n \
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\	Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);	fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}	}
fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm);	fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition and dummy covariate value at beginning of transition</h4>\n",fileresphtm, fileresphtm);

strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));	strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {	if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));	printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));	fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
fflush(ficlog);	fflush(ficlog);
exit(70);	exit(70);
}	} else{
else{
fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \	fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
<hr size=\"2\" color=\"#EC5E5E\"> \n\	<hr size=\"2\" color=\"#EC5E5E\"> \n \
Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\	Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);	fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
}	}
fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);	fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions of the model, by age at begin of transition, and covariate value at the begin of transition (if the covariate is a varying covariate) </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);

freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);	freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
j1=0;	j1=0;

/* j=ncoveff; /\* Only fixed dummy covariates \/ /	/* j=ncoveff; /\* Only fixed dummy covariates \/ /
j=cptcoveff; /* Only dummy covariates of the model */	j=cptcoveff; /* Only dummy covariates of the model */
if (cptcovn<1) {j=1;ncodemax[1]=1;}	if (cptcovn<1) {j=1;ncodemax[1]=1;}

first=1;	first=1;

/* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:	/* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
reference=low_education V1=0,V2=0	reference=low_education V1=0,V2=0
med_educ V1=1 V2=0,	med_educ V1=1 V2=0,
high_educ V1=0 V2=1	high_educ V1=0 V2=1
Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff	Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff
*/	*/

for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */	for (j1 = 1; j1 <= (int) pow(2,j); j1++){ /* Loop on covariates combination in order of model, excluding quantitatives V4=0, V3=0 for example, fixed or varying covariates */
posproptt=0.;	posproptt=0.;
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);	/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
scanf("%d", i);*/	scanf("%d", i);*/
for (i=-5; i<=nlstate+ndeath; i++)	for (i=-5; i<=nlstate+ndeath; i++)
for (jk=-5; jk<=nlstate+ndeath; jk++)	for (jk=-5; jk<=nlstate+ndeath; jk++)
for(m=iagemin; m <= iagemax+3; m++)	for(m=iagemin; m <= iagemax+3; m++)
freq[i][jk][m]=0;	freq[i][jk][m]=0;

for (i=1; i<=nlstate; i++) {	for (i=1; i<=nlstate; i++) {
for(m=iagemin; m <= iagemax+3; m++)	for(m=iagemin; m <= iagemax+3; m++)
prop[i][m]=0;	prop[i][m]=0;
posprop[i]=0;	posprop[i]=0;
pospropt[i]=0;	pospropt[i]=0;
}	}
Line 4162 Title=%s <br>Datafile=%s Firstpass=%d La	Line 4276 Title=%s <br>Datafile=%s Firstpass=%d La
/* meanqt[m][z1]=0.; */	/* meanqt[m][z1]=0.; */
/* } */	/* } */
/* } */	/* } */

dateintsum=0;	dateintsum=0;
k2cpt=0;	k2cpt=0;
/* For that combination of covariate j1, we count and print the frequencies in one pass */	/* For that combination of covariate j1, we count and print the frequencies in one pass */
Line 4241 Title=%s <br>Datafile=%s Firstpass=%d La	Line 4355 Title=%s <br>Datafile=%s Firstpass=%d La
} /* end iind = 1 to imx */	} /* end iind = 1 to imx */
/* prop[s][age] is feeded for any initial and valid live state as well as	/* prop[s][age] is feeded for any initial and valid live state as well as
freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */	freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */


/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/	/* fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
pstamp(ficresp);	pstamp(ficresp);
/* if (ncoveff>0) { */	if (cptcoveff>0){
if (cptcoveff>0) {
fprintf(ficresp, "\n#********** Variable ");	fprintf(ficresp, "\n#********** Variable ");
fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");	fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");
fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");	fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");
	fprintf(ficlog, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++){	for (z1=1; z1<=cptcoveff; z1++){
fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);	if(DummyV[z1]){
fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);	fprintf(ficresp, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);	fprintf(ficresphtm, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
	fprintf(ficresphtmfr, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
	fprintf(ficlog, "V%d (fixed)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
	}else{
	fprintf(ficresp, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
	fprintf(ficresphtm, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
	fprintf(ficresphtmfr, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
	fprintf(ficlog, "V%d(varying)=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
	}
}	}
fprintf(ficresp, "**********\n#");	fprintf(ficresp, "**********\n#");
fprintf(ficresphtm, "**********</h3>\n");	fprintf(ficresphtm, "**********</h3>\n");
fprintf(ficresphtmfr, "**********</h3>\n");	fprintf(ficresphtmfr, "**********</h3>\n");
fprintf(ficlog, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
fprintf(ficlog, "**********\n");	fprintf(ficlog, "**********\n");
}	}
fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");	fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
for(i=1; i<=nlstate;i++) {	for(i=1; i<=nlstate;i++) {
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);	fprintf(ficresp, " Age Prev(%d) N(%d) N ",i,i);
fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);	fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
}	}
fprintf(ficresp, "\n");	fprintf(ficresp, "\n");
fprintf(ficresphtm, "\n");	fprintf(ficresphtm, "\n");

/* Header of frequency table by age */	/* Header of frequency table by age */
fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");	fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
fprintf(ficresphtmfr,"<th>Age</th> ");	fprintf(ficresphtmfr,"<th>Age</th> ");
Line 4280 Title=%s <br>Datafile=%s Firstpass=%d La	Line 4400 Title=%s <br>Datafile=%s Firstpass=%d La
}	}
}	}
fprintf(ficresphtmfr, "\n");	fprintf(ficresphtmfr, "\n");

/* For each age */	/* For each age */
for(iage=iagemin; iage <= iagemax+3; iage++){	for(iage=iagemin; iage <= iagemax+3; iage++){
fprintf(ficresphtm,"<tr>");	fprintf(ficresphtm,"<tr>");
if(iage==iagemax+1){	if(iage==iagemax+1){
fprintf(ficlog,"1");	fprintf(ficlog,"1");
fprintf(ficresphtmfr,"<tr><th>0</th> ");	fprintf(ficresphtmfr,"<tr><th>0</th> ");
}else if(iage==iagemax+2){	}else if(iage==iagemax+2){
fprintf(ficlog,"0");	fprintf(ficlog,"0");
fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");	fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
}else if(iage==iagemax+3){	}else if(iage==iagemax+3){
fprintf(ficlog,"Total");	fprintf(ficlog,"Total");
fprintf(ficresphtmfr,"<tr><th>Total</th> ");	fprintf(ficresphtmfr,"<tr><th>Total</th> ");
}else{	}else{
if(first==1){	if(first==1){
first=0;	first=0;
printf("See log file for details...\n");	printf("See log file for details...\n");
}	}
fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);	fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
fprintf(ficlog,"Age %d", iage);	fprintf(ficlog,"Age %d", iage);
}	}
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][iage];	pp[jk] += freq[jk][m][iage];
}	}
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][iage];	pos += freq[jk][m][iage];
if(pp[jk]>=1.e-10){	if(pp[jk]>=1.e-10){
if(first==1){	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]);
}	}
fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);	fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
}else{	}else{
if(first==1)	if(first==1)
printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);	printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
fprintf(ficlog," %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++){
/* posprop[jk]=0; */	/* posprop[jk]=0; */
for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */	for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
pp[jk] += freq[jk][m][iage];	pp[jk] += freq[jk][m][iage];
} /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */	} /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */

for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){	for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){
pos += pp[jk]; /* pos is the total number of transitions until this age */	pos += pp[jk]; /* pos is the total number of transitions until this age */
posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state	posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */	from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state	pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */	from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
}	}
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
if(pos>=1.e-5){	if(pos>=1.e-5){
if(first==1)	if(first==1)
printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);	printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);	fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
}else{	}else{
if(first==1)	if(first==1)
printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);	printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);	fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
}	}
if( iage <= iagemax){	if( iage <= iagemax){
if(pos>=1.e-5){	if(pos>=1.e-5){
fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);	fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);	fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
/probs[iage][jk][j1]= pp[jk]/pos;/	/probs[iage][jk][j1]= pp[jk]/pos;/
/printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);/	/printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);/
}	}
else{	else{
fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);	fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);	fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
}	}
}	}
pospropt[jk] +=posprop[jk];	pospropt[jk] +=posprop[jk];
} /* end loop jk */	} /* end loop jk */
/* pospropt=0.; */	/* pospropt=0.; */
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][iage] !=0 ) { /* minimizing output */	if(freq[jk][m][iage] !=0 ) { /* minimizing output */
if(first==1){	if(first==1){
printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);	printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
}	}
fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);	fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
}	}
if(jk!=0 && m!=0)	if(jk!=0 && m!=0)
fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);	fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
}	}
} /* end loop jk */	} /* end loop jk */
posproptt=0.;	posproptt=0.;
for(jk=1; jk <=nlstate; jk++){	for(jk=1; jk <=nlstate; jk++){
posproptt += pospropt[jk];	posproptt += pospropt[jk];
}	}
fprintf(ficresphtmfr,"</tr>\n ");	fprintf(ficresphtmfr,"</tr>\n ");
if(iage <= iagemax){	if(iage <= iagemax){
fprintf(ficresp,"\n");	fprintf(ficresp,"\n");
fprintf(ficresphtm,"</tr>\n");	fprintf(ficresphtm,"</tr>\n");
}	}
if(first==1)	if(first==1)
printf("Others in log...\n");	printf("Others in log...\n");
fprintf(ficlog,"\n");	fprintf(ficlog,"\n");
} /* end loop age iage */	} /* end loop age iage */
fprintf(ficresphtm,"<tr><th>Tot</th>");	fprintf(ficresphtm,"<tr><th>Tot</th>");
for(jk=1; jk <=nlstate ; jk++){	for(jk=1; jk <=nlstate ; jk++){
if(posproptt < 1.e-5){	if(posproptt < 1.e-5){
fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);	fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);
}else{	}else{
fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);	fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);
}	}
}	}
fprintf(ficresphtm,"</tr>\n");	fprintf(ficresphtm,"</tr>\n");
Line 4406 Title=%s <br>Datafile=%s Firstpass=%d La	Line 4526 Title=%s <br>Datafile=%s Firstpass=%d La
fprintf(ficresphtmfr,"</table>\n");	fprintf(ficresphtmfr,"</table>\n");
} /* end selected combination of covariate j1 */	} /* end selected combination of covariate j1 */
dateintmean=dateintsum/k2cpt;	dateintmean=dateintsum/k2cpt;

fclose(ficresp);	fclose(ficresp);
fclose(ficresphtm);	fclose(ficresphtm);
fclose(ficresphtmfr);	fclose(ficresphtmfr);
Line 4741 void concatwav(int wav[], int **dh, int	Line 4861 void concatwav(int wav[], int **dh, int

/********* Tricode **************************/	/********* Tricode **************************/
void tricode(int cptcov, int Tvar, int *nbcode, int imx, int Ndum)	void tricode(int cptcov, int Tvar, int *nbcode, int imx, int Ndum)
{	{
/*< Uses cptcovn+2cptcovprod as the number of covariates */	/*< Uses cptcovn+2cptcovprod as the number of covariates */
/* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1	/* Tvar[i]=atoi(stre); find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1
* Boring subroutine which should only output nbcode[Tvar[j]][k]	* Boring subroutine which should only output nbcode[Tvar[j]][k]
* Tvar[5] in V2+V1+V3age+V2V4 is 4 (V4) even it is a time varying or quantitative variable	* Tvar[5] in V2+V1+V3age+V2V4 is 4 (V4) even it is a time varying or quantitative variable
* nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually);	* nbcode[Tvar[5]][1]= nbcode[4][1]=0, nbcode[4][2]=1 (usually);
*/	*/

int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;	int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
int modmaxcovj=0; /* Modality max of covariates j */	int modmaxcovj=0; /* Modality max of covariates j */
int cptcode=0; /* Modality max of covariates j */	int cptcode=0; /* Modality max of covariates j */
int modmincovj=0; /* Modality min of covariates j */	int modmincovj=0; /* Modality min of covariates j */


/* cptcoveff=0; */	/* cptcoveff=0; */
/* cptcov=0; /	/* cptcov=0; /

for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */	for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */

/* Loop on covariates without age and products and no quantitative variable */	/* Loop on covariates without age and products and no quantitative variable */
/* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2age+ V3 + V3V4 keeps V1 + V3 = 2 only \/ /	/* for (j=1; j<=(cptcovs); j++) { /\* From model V1 + V2age+ V3 + V3V4 keeps V1 + V3 = 2 only \/ /
for (k=1; k<=cptcovt; k++) { /* From model V1 + V2age + V3 + V3V4 keeps V1 + V3 = 2 only */	for (k=1; k<=cptcovt; k++) { /* From model V1 + V2age + V3 + V3V4 keeps V1 + V3 = 2 only */
for (j=-1; (j < maxncov); j++) Ndum[j]=0;	for (j=-1; (j < maxncov); j++) Ndum[j]=0;
if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */	if(Dummy[k]==0 && Typevar[k] !=1){ /* Dummy covariate and not age product */
switch(Fixed[k]) {	switch(Fixed[k]) {
case 0: /* Testing on fixed dummy covariate, simple or product of fixed */	case 0: /* Testing on fixed dummy covariate, simple or product of fixed */
for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/	for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the modality of this covariate Vj*/
ij=(int)(covar[Tvar[k]][i]);	ij=(int)(covar[Tvar[k]][i]);
/* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i	/* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
* If product of VnVm, still boolean :	* If product of VnVm, still boolean :
* If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables	* If it was coded 1, 2, 3, 4 should be splitted into 3 boolean variables
* 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */	* 1 => 0 0 0, 2 => 0 0 1, 3 => 0 1 1, 4=1 0 0 */
/* Finds for covariate j, n=Tvar[j] of Vn . ij is the	/* Finds for covariate j, n=Tvar[j] of Vn . ij is the
modality of the nth covariate of individual i. */	modality of the nth covariate of individual i. */
if (ij > modmaxcovj)	if (ij > modmaxcovj)
modmaxcovj=ij;	modmaxcovj=ij;
else if (ij < modmincovj)	else if (ij < modmincovj)
modmincovj=ij;	modmincovj=ij;
if ((ij < -1) && (ij > NCOVMAX)){	if ((ij < -1) && (ij > NCOVMAX)){
printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );	printf( "Error: minimal is less than -1 or maximal is bigger than %d. Exiting. \n", NCOVMAX );
exit(1);	exit(1);
}else	}else
Ndum[ij]++; /counts and stores the occurence of this modality 0, 1, -1/	Ndum[ij]++; /counts and stores the occurence of this modality 0, 1, -1/
/* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */	/* If coded 1, 2, 3 , counts the number of 1 Ndum[1], number of 2, Ndum[2], etc */
/printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);/	/printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);/
/* getting the maximum value of the modality of the covariate	/* getting the maximum value of the modality of the covariate
(should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and	(should be 0 or 1 now) Tvar[j]. If V=sex and male is coded 0 and
female ies 1, then modmaxcovj=1.	female ies 1, then modmaxcovj=1.
*/	*/
} /* end for loop on individuals i */	} /* end for loop on individuals i */
printf(" Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);	printf(" Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
fprintf(ficlog," Minimal and maximal values of %d th covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);	fprintf(ficlog," Minimal and maximal values of %d th (fixed) covariate V%d: min=%d max=%d \n", k, Tvar[k], modmincovj, modmaxcovj);
cptcode=modmaxcovj;	cptcode=modmaxcovj;
/* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */	/* Ndum[0] = frequency of 0 for model-covariate j, Ndum[1] frequency of 1 etc. */
/for (i=0; i<=cptcode; i++) {/	/for (i=0; i<=cptcode; i++) {/
for (j=modmincovj; j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1// For each value j of the modality of model-cov k */	for (j=modmincovj; j<=modmaxcovj; j++) { /* j=-1 ? 0 and 1// For each value j of the modality of model-cov k */
printf("Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);	printf("Frequencies of (fixed) covariate %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
fprintf(ficlog, "Frequencies of covariates %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);	fprintf(ficlog, "Frequencies of (fixed) covariate %d ie V%d with value %d: %d\n", k, Tvar[k], j, Ndum[j]);
if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */	if( Ndum[j] != 0 ){ /* Counts if nobody answered modality j ie empty modality, we skip it and reorder */
if( j != -1){	if( j != -1){
ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th	ncodemax[k]++; /* ncodemax[k]= Number of modalities of the k th
covariate for which somebody answered excluding	covariate for which somebody answered excluding
undefined. Usually 2: 0 and 1. */	undefined. Usually 2: 0 and 1. */
}	}
ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th	ncodemaxwundef[k]++; /* ncodemax[j]= Number of modalities of the k th
covariate for which somebody answered including	covariate for which somebody answered including
undefined. Usually 3: -1, 0 and 1. */	undefined. Usually 3: -1, 0 and 1. */
} /* In fact ncodemax[k]=2 (dichotom. variables only) but it could be more for	} /* In fact ncodemax[k]=2 (dichotom. variables only) but it could be more for
* historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */	* historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
} /* Ndum[-1] number of undefined modalities */	} /* Ndum[-1] number of undefined modalities */

/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */	/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
/* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */	/* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7. */
/* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */	/* If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125; */
/* modmincovj=3; modmaxcovj = 7; */	/* modmincovj=3; modmaxcovj = 7; */
/* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */	/* There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3; */
/* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */	/* which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10; */
/* defining two dummy variables: variables V1_1 and V1_2.*/	/* defining two dummy variables: variables V1_1 and V1_2.*/
/* nbcode[Tvar[j]][ij]=k; */	/* nbcode[Tvar[j]][ij]=k; */
/* nbcode[Tvar[j]][1]=0; */	/* nbcode[Tvar[j]][1]=0; */
/* nbcode[Tvar[j]][2]=1; */	/* nbcode[Tvar[j]][2]=1; */
/* nbcode[Tvar[j]][3]=2; */	/* nbcode[Tvar[j]][3]=2; */
/* To be continued (not working yet). */	/* To be continued (not working yet). */
ij=0; /* ij is similar to i but can jump over null modalities */	ij=0; /* ij is similar to i but can jump over null modalities */
for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/	for (i=modmincovj; i<=modmaxcovj; i++) { /* i= 1 to 2 for dichotomous, or from 1 to 3 or from -1 or 0 to 1 currently*/
if (Ndum[i] == 0) { /* If nobody responded to this modality k */	if (Ndum[i] == 0) { /* If nobody responded to this modality k */
break;	break;
}	}
ij++;	ij++;
nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/	nbcode[Tvar[k]][ij]=i; /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality. nbcode[1][1]=0 nbcode[1][2]=1*/
cptcode = ij; /* New max modality for covar j */	cptcode = ij; /* New max modality for covar j */
} /* end of loop on modality i=-1 to 1 or more */	} /* end of loop on modality i=-1 to 1 or more */
break;	break;
case 1: /* Testing on varying covariate, could be simple and	case 1: /* Testing on varying covariate, could be simple and
* should look at waves or product of fixed *	* should look at waves or product of fixed *
* varying. No time to test -1, assuming 0 and 1 only */	* varying. No time to test -1, assuming 0 and 1 only */
ij=0;	ij=0;
for(i=0; i<=1;i++){	for(i=0; i<=1;i++){
nbcode[Tvar[k]][++ij]=i;	nbcode[Tvar[k]][++ij]=i;
}	}
break;	break;
default:	default:
break;	break;
} /* end switch */	} /* end switch */
} /* end dummy test */	} /* end dummy test */

/* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 \//\ Could be 1 to 4 \//\ cptcode=modmaxcovj \/ /	/* for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 \//\ Could be 1 to 4 \//\ cptcode=modmaxcovj \/ /
/* /\recode from 0 \/ */	/* /\recode from 0 \/ */
/* k is a modality. If we have model=V1+V1sex /	/* k is a modality. If we have model=V1+V1sex /
/* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */	/* then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
/* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode \/ /	/* But if some modality were not used, it is recoded from 0 to a newer modmaxcovj=cptcode \/ /
/* } */	/* } */
/* /\* cptcode = ij; \/ /\ New max modality for covar j \/ /	/* /\* cptcode = ij; \/ /\ New max modality for covar j \/ /
/* if (ij > ncodemax[j]) { */	/* if (ij > ncodemax[j]) { */
/* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */	/* printf( " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
/* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */	/* fprintf(ficlog, " Error ij=%d > ncodemax[%d]=%d\n", ij, j, ncodemax[j]); */
/* break; */	/* break; */
/* } */	/* } */
/* } /\* end of loop on modality k \/ /	/* } /\* end of loop on modality k \/ /
} /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/	} /* end of loop on model-covariate j. nbcode[Tvarj][1]=0 and nbcode[Tvarj][2]=1 sets the value of covariate j*/

for (k=-1; k< maxncov; k++) Ndum[k]=0;	for (k=-1; k< maxncov; k++) Ndum[k]=0;
/* Look at fixed dummy (single or product) covariates to check empty modalities */	/* Look at fixed dummy (single or product) covariates to check empty modalities */
for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually ageage /	for (i=1; i<=ncovmodel-2-nagesqr; i++) { /* -2, cste and age and eventually ageage /
/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1V2./	/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1V2./
ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V4age /	ij=Tvar[i]; /* Tvar 5,4,3,6,5,7,1,4 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V4age /
Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */	Ndum[ij]++; /* Count the # of 1, 2 etc: {1,1,1,2,2,1,1} because V1 once, V2 once, two V4 and V5 in above */
/* V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1, {2, 1, 1, 1, 2, 1, 1, 0, 0} */	/* V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1, {2, 1, 1, 1, 2, 1, 1, 0, 0} */
} /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */	} /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */

ij=0;	ij=0;
/* for (i=0; i<= maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(ageV3), 5(V3V2) 6(V1V4) \/ */	/* for (i=0; i<= maxncov-1; i++) { /\* modmaxcovj is unknown here. Only Ndum[2(V2),3(ageV3), 5(V3V2) 6(V1V4) \/ */
for (k=1; k<= cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(ageV3), 5(V3V2) 6(V1V4) /	for (k=1; k<= cptcovt; k++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(ageV3), 5(V3V2) 6(V1V4) /
/printf("Ndum[%d]=%d\n",i, Ndum[i]);/	/printf("Ndum[%d]=%d\n",i, Ndum[i]);/
/* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? \/ /	/* if((Ndum[i]!=0) && (i<=ncovcol)){ /\* Tvar[i] <= ncovmodel ? \/ /
if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */	if(Ndum[Tvar[k]]!=0 && Dummy[k] == 0 && Typevar[k]==0){ /* Only Dummy and non empty in the model */
/* If product not in single variable we don't print results */	/* If product not in single variable we don't print results */
/printf("diff Ndum[%d]=%d\n",i, Ndum[i]);/	/printf("diff Ndum[%d]=%d\n",i, Ndum[i]);/
++ij;/* V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1, */	++ij;/* V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1, */
Tvaraff[ij]=Tvar[k]; /* For printing combination // V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/	Tvaraff[ij]=Tvar[k]; /* For printing combination // V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1, Tvar {5, 4, 3, 6, 5, 2, 7, 1, 1} Tvaraff={4, 3, 1} V4, V3, V1*/
Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */	Tmodelind[ij]=k; /* Tmodelind: index in model of dummies Tmodelind[1]=2 V4: pos=2; V3: pos=3, V1=9 {2, 3, 9, ?, ?,} */
TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */	TmodelInvind[ij]=Tvar[k]- ncovcol-nqv; /* Inverse TmodelInvind[2=V4]=2 second dummy varying cov (V4)4-1-1 {0, 2, 1, } TmodelInvind[3]=1 */
if(Fixed[k]!=0)	if(Fixed[k]!=0)
anyvaryingduminmodel=1;	anyvaryingduminmodel=1;
/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */	/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv)){ */
/* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet \/ /	/* Tvaraff[++ij]=-10; /\* Dont'n know how to treat quantitative variables yet \/ /
/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */	/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv)){ */
/* Tvaraff[++ij]=i; /\For printing (unclear) \/ */	/* Tvaraff[++ij]=i; /\For printing (unclear) \/ */
/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */	/* }else if((Ndum[i]!=0) && (i<=ncovcol+nqv+ntv+nqtv)){ */
/* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet \/ /	/* Tvaraff[++ij]=-20; /\* Dont'n know how to treat quantitative variables yet \/ /
}	}
} /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */	} /* Tvaraff[1]@5 {3, 4, -20, 0, 0} Very strange */
/* ij--; */	/* ij--; */
/* cptcoveff=ij; /\Number of total covariates\/ */	/* cptcoveff=ij; /\Number of total covariates\/ */
cptcov=ij; /Number of total real effective covariates: effective	cptcov=ij; /Number of total real effective covariates: effective
* because they can be excluded from the model and real	* because they can be excluded from the model and real
* if in the model but excluded because missing values, but how to get k from ij?*/	* if in the model but excluded because missing values, but how to get k from ij?*/
for(j=ij+1; j<= cptcovt; j++){	for(j=ij+1; j<= cptcovt; j++){
Tvaraff[j]=0;	Tvaraff[j]=0;
Tmodelind[j]=0;	Tmodelind[j]=0;
}	}
for(j=ntveff+1; j<= cptcovt; j++){	for(j=ntveff+1; j<= cptcovt; j++){
TmodelInvind[j]=0;	TmodelInvind[j]=0;
}	}
/* To be sorted */	/* To be sorted */
;	;
}	}


/********* Health Expectancies **************/	/********* Health Expectancies **************/
Line 5364 void concatwav(int wav[], int **dh, int	Line 5484 void concatwav(int wav[], int **dh, int
xp[i] = x[i] + (i==theta ?delti[theta]:0);	xp[i] = x[i] + (i==theta ?delti[theta]:0);
}	}

prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nresult);	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp,ij, nres);

if (popbased==1) {	if (popbased==1) {
if(mobilav ==0){	if(mobilav ==0){
Line 5396 void concatwav(int wav[], int **dh, int	Line 5516 void concatwav(int wav[], int **dh, int
for(i=1; i<=npar; i++) /* Computes gradient x - delta */	for(i=1; i<=npar; i++) /* Computes gradient x - delta */
xp[i] = x[i] - (i==theta ?delti[theta]:0);	xp[i] = x[i] - (i==theta ?delti[theta]:0);

prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nresult);	prevalim(prlim,nlstate,xp,age,oldm,savm,ftolpl,ncvyearp, ij, nres);

if (popbased==1) {	if (popbased==1) {
if(mobilav ==0){	if(mobilav ==0){
Line 5473 void concatwav(int wav[], int **dh, int	Line 5593 void concatwav(int wav[], int **dh, int
/* end ppptj */	/* end ppptj */
/* x centered again */	/* x centered again */

prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nresult);	prevalim(prlim,nlstate,x,age,oldm,savm,ftolpl,ncvyearp,ij, nres);

if (popbased==1) {	if (popbased==1) {
if(mobilav ==0){	if(mobilav ==0){
Line 5568 void concatwav(int wav[], int **dh, int	Line 5688 void concatwav(int wav[], int **dh, int

pstamp(ficresvpl);	pstamp(ficresvpl);
fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");	fprintf(ficresvpl,"# Standard deviation of period (stable) prevalences \n");
fprintf(ficresvpl,"# Age");	fprintf(ficresvpl,"# Age ");
	if(nresult >=1)
	fprintf(ficresvpl," Result# ");
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
fprintf(ficresvpl," %1d-%1d",i,i);	fprintf(ficresvpl," %1d-%1d",i,i);
fprintf(ficresvpl,"\n");	fprintf(ficresvpl,"\n");
Line 5654 void concatwav(int wav[], int **dh, int	Line 5776 void concatwav(int wav[], int **dh, int
varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */	varpl[i][(int)age] = doldm[i][i]; /* Covariances are useless */

fprintf(ficresvpl,"%.0f ",age );	fprintf(ficresvpl,"%.0f ",age );
	if(nresult >=1)
	fprintf(ficresvpl,"%d ",nres );
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));	fprintf(ficresvpl," %.5f (%.5f)",prlim[i][i],sqrt(varpl[i][(int)age]));
fprintf(ficresvpl,"\n");	fprintf(ficresvpl,"\n");
Line 6058 void printinghtml(char fileresu[], char	Line 6182 void printinghtml(char fileresu[], char
jj1=0;	jj1=0;

for(nres=1; nres <= nresult; nres++) /* For each resultline */	for(nres=1; nres <= nresult; nres++) /* For each resultline */
for(k1=1; k1<=m;k1++){	for(k1=1; k1<=m;k1++){ /* For each combination of covariate */
if(TKresult[nres]!= k1)	if(TKresult[nres]!= k1)
continue;	continue;

Line 6086 void printinghtml(char fileresu[], char	Line 6210 void printinghtml(char fileresu[], char
}	}
}	}
/* aij, bij */	/* aij, bij */
fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1.svg\">%s_%d-1.svg</a><br> \	fprintf(fichtm,"<br>- Logit model (yours is: 1+age+%s), for example: logit(pij)=log(pij/pii)= aij+ bij age + V1 age + etc. as a function of age: <a href=\"%s_%d-1-%d.svg\">%s_%d-1-%d.svg</a><br> \
<img src=\"%s_%d-1.svg\">",model,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);	<img src=\"%s_%d-1-%d.svg\">",model,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Pij */	/* Pij */
fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2.svg\">%s_%d-2.svg</a><br> \	fprintf(fichtm,"<br>\n- P<sub>ij</sub> or conditional probabilities to be observed in state j being in state i, %d (stepm) months before: <a href=\"%s_%d-2-%d.svg\">%s_%d-2-%d.svg</a><br> \
<img src=\"%s_%d-2.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);	<img src=\"%s_%d-2-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Quasi-incidences */	/* Quasi-incidences */
fprintf(fichtm,"<br>\n- I<sub>ij</sub> or Conditional probabilities to be observed in state j being in state i %d (stepm) months\	fprintf(fichtm,"<br>\n- I<sub>ij</sub> 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, \	before but expressed in per year i.e. quasi incidences if stepm is small and probabilities too, \
incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \	incidence (rates) are the limit when h tends to zero of the ratio of the probability <sub>h</sub>P<sub>ij</sub> \
divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3.svg\">%s_%d-3.svg</a><br> \	divided by h: <sub>h</sub>P<sub>ij</sub>/h : <a href=\"%s_%d-3-%d.svg\">%s_%d-3-%d.svg</a><br> \
<img src=\"%s_%d-3.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1,subdirf2(optionfilefiname,"PE_"),jj1);	<img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Survival functions (period) in state j */	/* Survival functions (period) in state j */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \	fprintf(fichtm,"<br>\n- Survival functions in state %d. Or probability to survive in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1,subdirf2(optionfilefiname,"LIJ_"),cpt,jj1);	<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
}	}
/* State specific survival functions (period) */	/* State specific survival functions (period) */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\	fprintf(fichtm,"<br>\n- Survival functions from state %d in each live state and total.\
Or probability to survive in various states (1 to %d) being in state %d at different ages. \	Or probability to survive in various states (1 to %d) being in state %d at different ages. \
<a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> <img src=\"%s_%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1,subdirf2(optionfilefiname,"LIJT_"),cpt,jj1);	<a href=\"%s_%d-%d-%d.svg\">%s_%d%d-%d.svg</a><br> <img src=\"%s_%d-%d-%d.svg\">", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
}	}
/* Period (stable) prevalence in each health state */	/* Period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \	fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1,subdirf2(optionfilefiname,"P_"),cpt,jj1);	<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
}	}
if(backcast==1){	if(backcast==1){
/* Period (stable) back prevalence in each health state */	/* Period (stable) back prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a><br> \	fprintf(fichtm,"<br>\n- Convergence to period (stable) back prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1,subdirf2(optionfilefiname,"PB_"),cpt,jj1);	<img src=\"%s_%d-%d-%d.svg\">", cpt, cpt, nlstate, subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
}	}
}	}
if(prevfcast==1){	if(prevfcast==1){
/* Projection of prevalence up to period (stable) prevalence in each health state */	/* Projection of prevalence up to period (stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s%d_%d.svg\">%s%d_%d.svg</a><br> \	fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f) up to period (stable) prevalence in state %d. Or probability to be in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
<img src=\"%s_%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1,subdirf2(optionfilefiname,"PROJ_"),cpt,jj1);	<img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
}	}
}	}

for(cpt=1; cpt<=nlstate;cpt++) {	for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d%d.svg\">%s_%d%d.svg</a> <br> \	fprintf(fichtm,"\n<br>- Life expectancy by health state (%d) at initial age and its decomposition into health expectancies in each alive state (1 to %d) (or area under each survival functions): <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a> <br> \
<img src=\"%s_%d%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1,subdirf2(optionfilefiname,"EXP_"),cpt,jj1);	<img src=\"%s_%d-%d-%d.svg\">",cpt,nlstate,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres,subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
}	}
/* } /\* end i1 \/ /	/* } /\* end i1 \/ /
}/* End k1 */	}/* End k1 */
Line 6189 See page 'Matrix of variance-covariance	Line 6313 See page 'Matrix of variance-covariance

jj1=0;	jj1=0;

for(nres=1; nres <= nresult; nres++) /* For each resultline */	for(nres=1; nres <= nresult; nres++){ /* For each resultline */
for(k1=1; k1<=m;k1++){	for(k1=1; k1<=m;k1++){
if(TKresult[nres]!= k1)	if(TKresult[nres]!= k1)
continue;	continue;
Line 6213 See page 'Matrix of variance-covariance	Line 6337 See page 'Matrix of variance-covariance
}	}
for(cpt=1; cpt<=nlstate;cpt++) {	for(cpt=1; cpt<=nlstate;cpt++) {
fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \	fprintf(fichtm,"\n<br>- Observed (cross-sectional) and period (incidence based) \
prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d.svg\"> %s_%d-%d.svg</a>\n <br>\	prevalence (with 95%% confidence interval) in state (%d): <a href=\"%s_%d-%d-%d.svg\"> %s_%d-%d-%d.svg</a>\n <br>\
<img src=\"%s_%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1,subdirf2(optionfilefiname,"V_"),cpt,jj1);	<img src=\"%s_%d-%d-%d.svg\">",cpt,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres,subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
}	}
fprintf(fichtm,"\n<br>- Total life expectancy by age and \	fprintf(fichtm,"\n<br>- Total life expectancy by age and \
health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \	health expectancies in states (1) and (2). If popbased=1 the smooth (due to the model) \
true period expectancies (those weighted with period prevalences are also\	true period expectancies (those weighted with period prevalences are also\
drawn in addition to the population based expectancies computed using\	drawn in addition to the population based expectancies computed using\
observed and cahotic prevalences: <a href=\"%s_%d.svg\">%s_%d.svg</a>\n<br>\	observed and cahotic prevalences: <a href=\"%s_%d-%d.svg\">%s_%d-%d.svg</a>\n<br>\
<img src=\"%s_%d.svg\">",subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1,subdirf2(optionfilefiname,"E_"),jj1);	<img src=\"%s_%d-%d.svg\">",subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres,subdirf2(optionfilefiname,"E_"),k1,nres);
/* } /\* end i1 \/ /	/* } /\* end i1 \/ /
}/* End k1 */	}/* End k1 */
	}/* End nres */
fprintf(fichtm,"</ul>");	fprintf(fichtm,"</ul>");
fflush(fichtm);	fflush(fichtm);
}	}
Line 6289 void printinggnuplot(char fileresu[], ch	Line 6414 void printinggnuplot(char fileresu[], ch
if(TKresult[nres]!= k1)	if(TKresult[nres]!= k1)
continue;	continue;
/* We are interested in selected combination by the resultline */	/* We are interested in selected combination by the resultline */
printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);	/* printf("\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt); */
fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);	fprintf(ficgp,"\n# 1st: Period (stable) prevalence with CI: 'VPL_' files and live state =%d ", cpt);
for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */	for (k=1; k<=cptcoveff; k++){ /* For each covariate k get corresponding value lv for combination k1 */
lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */	lv= decodtabm(k1,k,cptcoveff); /* Should be the value of the covariate corresponding to k1 combination */
Line 6298 void printinggnuplot(char fileresu[], ch	Line 6423 void printinggnuplot(char fileresu[], ch
/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */	/* decodtabm(13,3,4)= 2 because h=13 k= 1 1 (2) 2 */
vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */	vlv= nbcode[Tvaraff[k]][lv]; /* vlv is the value of the covariate lv, 0 or 1 */
/* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */	/* For each combination of covariate k1 (V1=1, V3=0), we printed the current covariate k and its value vlv */
printf(" V%d=%d ",Tvaraff[k],vlv);	/* printf(" V%d=%d ",Tvaraff[k],vlv); */
fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);	fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
}	}
for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */	for (k4=1; k4<= nsq; k4++){ /* For each selected (single) quantitative value */
printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);	/* printf(" V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]); */
fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
}	}
printf("\n#\n");	/* printf("\n#\n"); */
fprintf(ficgp,"\n#\n");	fprintf(ficgp,"\n#\n");
if(invalidvarcomb[k1]){	if(invalidvarcomb[k1]){
fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);	fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);
continue;	continue;
}	}

fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"V_"),cpt,k1,nres);
fprintf(ficgp,"\n#set out \"V_%s_%d-%d.svg\" \n",optionfilefiname,cpt,k1);	fprintf(ficgp,"\n#set out \"V_%s_%d-%d-%d.svg\" \n",optionfilefiname,cpt,k1,nres);
fprintf(ficgp,"set xlabel \"Age\" \n\	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \nset ter svg size 640, 480\nplot [%.f:%.f] \"%s\" every :::%d::%d u 1:($2==%d ? $3:1/0) \"%%lf %%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
set ylabel \"Probability\" \n \
set ter svg size 640, 480\n \
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileresu,"VPL_"),k1-1,k1-1);

for (i=1; i<= nlstate ; i ++) {	for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");	if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%lf (%%lf)");	else fprintf(ficgp," %%lf (%%lf)");
}	}
fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);	fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
for (i=1; i<= nlstate ; i ++) {	for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");	if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%lf (%%lf)");	else fprintf(ficgp," %%lf (%%lf)");
}	}
fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1);	fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VPL_"),k1-1,k1-1,nres);
for (i=1; i<= nlstate ; i ++) {	for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," %%lf (%%lf)");	if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
else fprintf(ficgp," %%lf (%%lf)");	else fprintf(ficgp," %%lf (%%lf)");
Line 6336 plot [%.f:%.f] \"%s\" every :::%d::%d u	Line 6458 plot [%.f:%.f] \"%s\" every :::%d::%d u
fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));	fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence\" w l lt 2",subdirf2(fileresu,"P_"),k1-1,k1-1,2+4*(cpt-1));
if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */	if(backcast==1){ /* We need to get the corresponding values of the covariates involved in this combination k1 */
/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */	/* fprintf(ficgp,",\"%s\" every :::%d::%d u 1:($%d) t\"Backward stable prevalence\" w l lt 3",subdirf2(fileresu,"PLB_"),k1-1,k1-1,1+cpt); */
fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1 */	fprintf(ficgp,",\"%s\" u 1:((",subdirf2(fileresu,"PLB_")); /* Age is in 1, nres in 2 to be fixed */
if(cptcoveff ==0){	if(cptcoveff ==0){
fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line ", 2+(cpt-1), cpt );	fprintf(ficgp,"$%d)) t 'Backward prevalence in state %d' with line lt 3", 2+(cpt-1), cpt );
}else{	}else{
kl=0;	kl=0;
for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */	for (k=1; k<=cptcoveff; k++){ /* For each combination of covariate */
Line 6353 plot [%.f:%.f] \"%s\" every :::%d::%d u	Line 6475 plot [%.f:%.f] \"%s\" every :::%d::%d u
/6+1+(i-1)+(nlstate+1)nlstate; 6+1+(1-1) +(2+1)2=13 /	/6+1+(i-1)+(nlstate+1)nlstate; 6+1+(1-1) +(2+1)2=13 /
/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/	/* '' u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
if(k==cptcoveff){	if(k==cptcoveff){
fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \	fprintf(ficgp,"$%d==%d && $%d==%d)? $%d : 1/0) t 'Backward prevalence in state %d' w l lt 3",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv], \
4+(cpt-1), cpt ); /* 4 or 6 ?*/	2+cptcoveff2+(cpt-1), cpt ); / 4 or 6 ?*/
}else{	}else{
fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);	fprintf(ficgp,"$%d==%d && $%d==%d && ",kl+1, Tvaraff[k],kl+1+1,nbcode[Tvaraff[k]][lv]);
kl++;	kl++;
Line 6393 plot [%.f:%.f] \"%s\" every :::%d::%d u	Line 6515 plot [%.f:%.f] \"%s\" every :::%d::%d u
continue;	continue;
}	}

fprintf(ficgp,"\nset out \"%s_%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1);	fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"E_"),k1,nres);
for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/	for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
if(vpopbased==0)	if(vpopbased==0)
fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);	fprintf(ficgp,"set ylabel \"Years\" \nset ter svg size 640, 480\nplot [%.f:%.f] ",ageminpar,fage);
Line 6423 plot [%.f:%.f] \"%s\" every :::%d::%d u	Line 6545 plot [%.f:%.f] \"%s\" every :::%d::%d u
else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");	else fprintf(ficgp,"\" t\"\" w l lt 0,\\\n");
} /* state */	} /* state */
} /* vpopbased */	} /* vpopbased */
fprintf(ficgp,"\nset out;set out \"%s_%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1); /* Buggy gnuplot */	fprintf(ficgp,"\nset out;set out \"%s_%d-%d.svg\"; replot; set out; \n",subdirf2(optionfilefiname,"E_"),k1,nres); /* Buggy gnuplot */
} /* end nres */	} /* end nres */
} /* k1 end 2 eme*/	} /* k1 end 2 eme*/

Line 6431 plot [%.f:%.f] \"%s\" every :::%d::%d u	Line 6553 plot [%.f:%.f] \"%s\" every :::%d::%d u
/3eme/	/3eme/
for (k1=1; k1<= m ; k1 ++){	for (k1=1; k1<= m ; k1 ++){
for(nres=1; nres <= nresult; nres++){ /* For each resultline */	for(nres=1; nres <= nresult; nres++){ /* For each resultline */
if(TKresult[nres]!= k)	if(TKresult[nres]!= k1)
continue;	continue;

for (cpt=1; cpt<= nlstate ; cpt ++) {	for (cpt=1; cpt<= nlstate ; cpt ++) {
Line 6455 plot [%.f:%.f] \"%s\" every :::%d::%d u	Line 6577 plot [%.f:%.f] \"%s\" every :::%d::%d u

/* k=2+nlstate(2cpt-2); */	/* k=2+nlstate(2cpt-2); */
k=2+(nlstate+1)*(cpt-1);	k=2+(nlstate+1)*(cpt-1);
fprintf(ficgp,"\nset out \"%s_%d%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"EXP_"),cpt,k1,nres);
fprintf(ficgp,"set ter svg size 640, 480\n\	fprintf(ficgp,"set ter svg size 640, 480\n\
plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);	plot [%.f:%.f] \"%s\" every :::%d::%d u 1:%d t \"e%d1\" w l",ageminpar,fage,subdirf2(fileresu,"E_"),k1-1,k1-1,k,cpt);
/fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);	/fprintf(ficgp,",\"e%s\" every :::%d::%d u 1:($%d-2$%d) \"\%%lf ",fileres,k1-1,k1-1,k,k+1);
Line 6501 plot [%.f:%.f] \"%s\" every :::%d::%d u	Line 6623 plot [%.f:%.f] \"%s\" every :::%d::%d u
continue;	continue;
}	}

fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);	set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
k=3;	k=3;
Line 6547 set ter svg size 640, 480\nunset log y\n	Line 6669 set ter svg size 640, 480\nunset log y\n
continue;	continue;
}	}

fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability to be alive\" \n\
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);	set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
k=3;	k=3;
Line 6602 set ter svg size 640, 480\nunset log y\n	Line 6724 set ter svg size 640, 480\nunset log y\n
continue;	continue;
}	}

fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);	set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
k=3; /* Offset */	k=3; /* Offset */
Line 6648 set ter svg size 640, 480\nunset log y\n	Line 6770 set ter svg size 640, 480\nunset log y\n
continue;	continue;
}	}

fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1,nres);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);	set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
k=3; /* Offset */	k=3; /* Offset */
Line 6700 set ter svg size 640, 480\nunset log y\n	Line 6822 set ter svg size 640, 480\nunset log y\n
}	}

fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");	fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\	fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);	set ter svg size 640, 480\nunset log y\nplot [%.f:%.f] ", ageminpar, agemaxpar);
for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */	for (i=1; i<= nlstate+1 ; i ++){ /* nlstate +1 p11 p21 p.1 */
Line 6816 set ter svg size 640, 480\nunset log y\n	Line 6938 set ter svg size 640, 480\nunset log y\n
fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);	fprintf(ficgp," V%d=%f ",Tvqresult[nres][k4],Tqresult[nres][k4]);
}	}
fprintf(ficgp,"\n#\n");	fprintf(ficgp,"\n#\n");
fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng);	fprintf(ficgp,"\nset out \"%s_%d-%d-%d.svg\" ",subdirf2(optionfilefiname,"PE_"),jk,ng,nres);
fprintf(ficgp,"\nset ter svg size 640, 480 ");	fprintf(ficgp,"\nset ter svg size 640, 480 ");
if (ng==1){	if (ng==1){
fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12x) could be nice /	fprintf(ficgp,"\nset ylabel \"Value of the logit of the model\"\n"); /* exp(a12+b12x) could be nice /
Line 7791 int readdata(char datafile[], int firsto	Line 7913 int readdata(char datafile[], int firsto
/-------- data file ----------/	/-------- data file ----------/
FILE *fic;	FILE *fic;
char dummy[]=" ";	char dummy[]=" ";
int i=0, j=0, n=0, iv=0;	int i=0, j=0, n=0, iv=0, v;
int lstra;	int lstra;
int linei, month, year,iout;	int linei, month, year,iout;
char line[MAXLINE], linetmp[MAXLINE];	char line[MAXLINE], linetmp[MAXLINE];
char stra[MAXLINE], strb[MAXLINE];	char stra[MAXLINE], strb[MAXLINE];
char *stratrunc;	char *stratrunc;

	DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
	FixedV=ivector(1,NCOVMAX); /* 1 to 3 */

	for(v=1; v <=ncovcol;v++){
	DummyV[v]=0;
	FixedV[v]=0;
	}
	for(v=ncovcol+1; v <=ncovcol+nqv;v++){
	DummyV[v]=1;
	FixedV[v]=0;
	}
	for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
	DummyV[v]=0;
	FixedV[v]=1;
	}
	for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
	DummyV[v]=1;
	FixedV[v]=1;
	}
	for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
	printf("Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
	fprintf(ficlog,"Covariate type in the data: V%d, DummyV(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
	}

if((fic=fopen(datafile,"r"))==NULL) {	if((fic=fopen(datafile,"r"))==NULL) {
printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);	printf("Problem while opening datafile: %s with errno='%s'\n", datafile,strerror(errno));fflush(stdout);
Line 8419 Dummy[k] 0=dummy (0 1), 1 quantitative (	Line 8563 Dummy[k] 0=dummy (0 1), 1 quantitative (
Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\	Typevar: 0 for simple covariate (dummy, quantitative, fixed or varying), 1 for age product, 2 for product \n\
Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\	Fixed[k] 0=fixed (product or simple), 1 varying, 2 fixed with age product, 3 varying with age product \n\
Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);	Dummy[k] 0=dummy (0 1), 1 quantitative (single or product without age), 2 dummy with age product, 3 quant with age product\n",model);
	for(k=1;k<=cptcovt; k++){ Fixed[k]=0; Dummy[k]=0;}
for(v=1; v <=ncovcol;v++){
DummyV[v]=0;
FixedV[v]=0;
}
for(v=ncovcol+1; v <=ncovcol+nqv;v++){
DummyV[v]=1;
FixedV[v]=0;
}
for(v=ncovcol+nqv+1; v <=ncovcol+nqv+ntv;v++){
DummyV[v]=0;
FixedV[v]=1;
}
for(v=ncovcol+nqv+ntv+1; v <=ncovcol+nqv+ntv+nqtv;v++){
DummyV[v]=1;
FixedV[v]=1;
}
for(v=1; v <=ncovcol+nqv+ntv+nqtv;v++){
printf("Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
fprintf(ficlog,"Decodemodel: V%d, Dummy(V%d)=%d, FixedV(V%d)=%d\n",v,v,DummyV[v],v,FixedV[v]);
}
for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */	for(k=1, ncovf=0, nsd=0, nsq=0, ncovv=0, ncova=0, ncoveff=0, nqfveff=0, ntveff=0, nqtveff=0;k<=cptcovt; k++){ /* or cptocvt */
if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */	if (Tvar[k] <=ncovcol && Typevar[k]==0 ){ /* Simple fixed dummy (<=ncovcol) covariates */
Fixed[k]= 0;	Fixed[k]= 0;
Line 8464 Dummy[k] 0=dummy (0 1), 1 quantitative (	Line 8588 Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarFind[ncovf]=k;	TvarFind[ncovf]=k;
TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 */	TvarFD[ncoveff]=Tvar[k]; /* TvarFD[1]=V1 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 */
TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 */	TvarFDind[ncoveff]=k; /* TvarFDind[1]=9 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 */
}else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){ /* Remind that product VnVm are added in k/ /* Only simple fixed quantitative variable */	}else if( Tvar[k] <=ncovcol+nqv && Typevar[k]==0){/* Remind that product VnVm are added in k Only simple fixed quantitative variable /
Fixed[k]= 0;	Fixed[k]= 0;
Dummy[k]= 1;	Dummy[k]= 1;
nqfveff++;	nqfveff++;
Line 8478 Dummy[k] 0=dummy (0 1), 1 quantitative (	Line 8602 Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarFind[ncovf]=k;	TvarFind[ncovf]=k;
TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple fixed quantitative variable */	TvarFQ[nqfveff]=Tvar[k]-ncovcol; /* TvarFQ[1]=V2-1=1st in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple fixed quantitative variable */
TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple fixed quantitative variable */	TvarFQind[nqfveff]=k; /* TvarFQind[1]=6 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple fixed quantitative variable */
}else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying variables */	}else if( Tvar[k] <=ncovcol+nqv+ntv && Typevar[k]==0){/* Only simple time varying dummy variables */
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 0;	Dummy[k]= 0;
ntveff++; /* Only simple time varying dummy variable */	ntveff++; /* Only simple time varying dummy variable */
Line 8489 Dummy[k] 0=dummy (0 1), 1 quantitative (	Line 8613 Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarsDind[nsd]=k;	TvarsDind[nsd]=k;
ncovv++; /* Only simple time varying variables */	ncovv++; /* Only simple time varying variables */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k; /* TvarVind[2]=2 TvarVind[3]=3 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Any time varying singele */
TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4 TvarVD[2]=V3 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple time varying dummy variable */	TvarVD[ntveff]=Tvar[k]; /* TvarVD[1]=V4 TvarVD[2]=V3 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple time varying dummy variable */
TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple time varying dummy variable */	TvarVDind[ntveff]=k; /* TvarVDind[1]=2 TvarVDind[2]=3 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple time varying dummy variable */
printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv);	printf("Quasi Tmodelind[%d]=%d,Tvar[Tmodelind[%d]]=V%d, ncovcol=%d, nqv=%d,Tvar[k]- ncovcol-nqv=%d\n",ntveff,k,ntveff,Tvar[k], ncovcol, nqv,Tvar[k]- ncovcol-nqv);
Line 8505 Dummy[k] 0=dummy (0 1), 1 quantitative (	Line 8629 Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarsQ[nsq]=Tvar[k];	TvarsQ[nsq]=Tvar[k];
TvarsQind[nsq]=k;	TvarsQind[nsq]=k;
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k; /* TvarVind[1]=1 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Any time varying singele */
TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple time varying quantitative variable */	TvarVQ[nqtveff]=Tvar[k]; /* TvarVQ[1]=V5 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple time varying quantitative variable */
TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple time varying quantitative variable */	TvarVQind[nqtveff]=k; /* TvarVQind[1]=1 in V5+V4+V3+V4V3+V5age+V2+V1V2+V1age+V1 / / Only simple time varying quantitative variable */
TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */	TmodelInvQind[nqtveff]=Tvar[k]- ncovcol-nqv-ntv;/* Only simple time varying quantitative variable */
Line 8517 Dummy[k] 0=dummy (0 1), 1 quantitative (	Line 8641 Dummy[k] 0=dummy (0 1), 1 quantitative (
TvarA[ncova]=Tvar[k];	TvarA[ncova]=Tvar[k];
TvarAind[ncova]=k;	TvarAind[ncova]=k;
if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */	if (Tvar[k] <=ncovcol ){ /* Product age with fixed dummy covariatee */
Fixed[k]= 2;	Fixed[k]= 2;
Dummy[k]= 2;	Dummy[k]= 2;
modell[k].maintype= ATYPE;	modell[k].maintype= ATYPE;
modell[k].subtype= APFD;	modell[k].subtype= APFD;
/* ncoveff++; */	/* ncoveff++; */
}else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product VnVm are added in k/	}else if( Tvar[k] <=ncovcol+nqv) { /* Remind that product VnVm are added in k/
Fixed[k]= 2;	Fixed[k]= 2;
Dummy[k]= 3;	Dummy[k]= 3;
modell[k].maintype= ATYPE;	modell[k].maintype= ATYPE;
modell[k].subtype= APFQ; /* Product age * fixed quantitative */	modell[k].subtype= APFQ; /* Product age * fixed quantitative */
/* nqfveff++; /\* Only simple fixed quantitative variable \/ /	/* nqfveff++; /\* Only simple fixed quantitative variable \/ /
}else if( Tvar[k] <=ncovcol+nqv+ntv ){	}else if( Tvar[k] <=ncovcol+nqv+ntv ){
Fixed[k]= 3;	Fixed[k]= 3;
Dummy[k]= 2;	Dummy[k]= 2;
modell[k].maintype= ATYPE;	modell[k].maintype= ATYPE;
modell[k].subtype= APVD; /* Product age * varying dummy */	modell[k].subtype= APVD; /* Product age * varying dummy */
/* ntveff++; /\* Only simple time varying dummy variable \/ /	/* ntveff++; /\* Only simple time varying dummy variable \/ /
}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){	}else if( Tvar[k] <=ncovcol+nqv+ntv+nqtv){
Fixed[k]= 3;	Fixed[k]= 3;
Dummy[k]= 3;	Dummy[k]= 3;
modell[k].maintype= ATYPE;	modell[k].maintype= ATYPE;
modell[k].subtype= APVQ; /* Product age * varying quantitative */	modell[k].subtype= APVQ; /* Product age * varying quantitative */
/* nqtveff++;/\* Only simple time varying quantitative variable \/ /	/* nqtveff++;/\* Only simple time varying quantitative variable \/ /
}	}
}else if (Typevar[k] == 2) { /* product without age */	}else if (Typevar[k] == 2) { /* product without age */
k1=Tposprod[k];	k1=Tposprod[k];
if(Tvard[k1][1] <=ncovcol){	if(Tvard[k1][1] <=ncovcol){
if(Tvard[k1][2] <=ncovcol){	if(Tvard[k1][2] <=ncovcol){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 0;	Dummy[k]= 0;
modell[k].maintype= FTYPE;	modell[k].maintype= FTYPE;
modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */	modell[k].subtype= FPDD; /* Product fixed dummy * fixed dummy */
ncovf++; /* Fixed variables without age */	ncovf++; /* Fixed variables without age */
TvarF[ncovf]=Tvar[k];	TvarF[ncovf]=Tvar[k];
TvarFind[ncovf]=k;	TvarFind[ncovf]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv){	}else if(Tvard[k1][2] <=ncovcol+nqv){
Fixed[k]= 0; /* or 2 ?*/	Fixed[k]= 0; /* or 2 ?*/
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= FTYPE;	modell[k].maintype= FTYPE;
modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */	modell[k].subtype= FPDQ; /* Product fixed dummy * fixed quantitative */
ncovf++; /* Varying variables without age */	ncovf++; /* Varying variables without age */
TvarF[ncovf]=Tvar[k];	TvarF[ncovf]=Tvar[k];
TvarFind[ncovf]=k;	TvarFind[ncovf]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 0;	Dummy[k]= 0;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */	modell[k].subtype= VPDD; /* Product fixed dummy * varying dummy */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */	modell[k].subtype= VPDQ; /* Product fixed dummy * varying quantitative */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}	}
}else if(Tvard[k1][1] <=ncovcol+nqv){	}else if(Tvard[k1][1] <=ncovcol+nqv){
if(Tvard[k1][2] <=ncovcol){	if(Tvard[k1][2] <=ncovcol){
Fixed[k]= 0; /* or 2 ?*/	Fixed[k]= 0; /* or 2 ?*/
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= FTYPE;	modell[k].maintype= FTYPE;
modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */	modell[k].subtype= FPDQ; /* Product fixed quantitative * fixed dummy */
ncovf++; /* Fixed variables without age */	ncovf++; /* Fixed variables without age */
TvarF[ncovf]=Tvar[k];	TvarF[ncovf]=Tvar[k];
TvarFind[ncovf]=k;	TvarFind[ncovf]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */	modell[k].subtype= VPDQ; /* Product fixed quantitative * varying dummy */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */	modell[k].subtype= VPQQ; /* Product fixed quantitative * varying quantitative */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}	}
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv){	}else if(Tvard[k1][1] <=ncovcol+nqv+ntv){
if(Tvard[k1][2] <=ncovcol){	if(Tvard[k1][2] <=ncovcol){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */	modell[k].subtype= VPDD; /* Product time varying dummy * fixed dummy */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv){	}else if(Tvard[k1][2] <=ncovcol+nqv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */	modell[k].subtype= VPDQ; /* Product time varying dummy * fixed quantitative */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 0;	Dummy[k]= 0;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */	modell[k].subtype= VPDD; /* Product time varying dummy * time varying dummy */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */	modell[k].subtype= VPDQ; /* Product time varying dummy * time varying quantitative */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}	}
}else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){	}else if(Tvard[k1][1] <=ncovcol+nqv+ntv+nqtv){
if(Tvard[k1][2] <=ncovcol){	if(Tvard[k1][2] <=ncovcol){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */	modell[k].subtype= VPDQ; /* Product time varying quantitative * fixed dummy */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv){	}else if(Tvard[k1][2] <=ncovcol+nqv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */	modell[k].subtype= VPQQ; /* Product time varying quantitative * fixed quantitative */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */	modell[k].subtype= VPDQ; /* Product time varying quantitative * time varying dummy */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){	}else if(Tvard[k1][2] <=ncovcol+nqv+ntv+nqtv){
Fixed[k]= 1;	Fixed[k]= 1;
Dummy[k]= 1;	Dummy[k]= 1;
modell[k].maintype= VTYPE;	modell[k].maintype= VTYPE;
modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */	modell[k].subtype= VPQQ; /* Product time varying quantitative * time varying quantitative */
ncovv++; /* Varying variables without age */	ncovv++; /* Varying variables without age */
TvarV[ncovv]=Tvar[k];	TvarV[ncovv]=Tvar[k];
TvarVind[ncovv]=k;	TvarVind[ncovv]=k;
}	}
}else{	}else{
printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);	printf("Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);	fprintf(ficlog,"Error unknown type of covariate: Tvard[%d][1]=%d,Tvard[%d][2]=%d\n",k1,Tvard[k1][1],k1,Tvard[k1][2]);
} /* end k1 */	} /end k1/
}else{	}else{
printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);	printf("Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);	fprintf(ficlog,"Error, current version can't treat for performance reasons, Tvar[%d]=%d, Typevar[%d]=%d\n", k, Tvar[k], k, Typevar[k]);
Line 9192 int back_prevalence_limit(double *p, dou	Line 9316 int back_prevalence_limit(double *p, dou
if(mobilavproj > 0){	if(mobilavproj > 0){
/* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */	/* bprevalim(bprlim, mobaverage, nlstate, p, age, ageminpar, agemaxpar, oldm, savm, doldm, dsavm, ftolpl, ncvyearp, k); */
/* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */	/* bprevalim(bprlim, mobaverage, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k);	bprevalim(bprlim, mobaverage, nlstate, p, age, ftolpl, ncvyearp, k, nres);
}else if (mobilavproj == 0){	}else if (mobilavproj == 0){
printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);	printf("There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);	fprintf(ficlog,"There is no chance to get back prevalence limit if data aren't non zero and summing to 1, please try a non null mobil_average(=%d) parameter or mobil_average=-1 if you want to try at your own risk.\n",mobilavproj);
exit(1);	exit(1);
}else{	}else{
/* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */	/* bprevalim(bprlim, probs, nlstate, p, age, oldm, savm, dnewm, doldm, dsavm, ftolpl, ncvyearp, k); */
bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k);	bprevalim(bprlim, probs, nlstate, p, age, ftolpl, ncvyearp, k,nres);
}	}
fprintf(ficresplb,"%.0f ",age );	fprintf(ficresplb,"%.0f ",age );
for(j=1;j<=cptcoveff;j++)	for(j=1;j<=cptcoveff;j++)
Line 9763 int main(int argc, char *argv[])	Line 9887 int main(int argc, char *argv[])
delti=delti3[1][1];	delti=delti3[1][1];
/delti=vector(1,npar); //* Scale of each paramater (output from hesscov)*/	/delti=vector(1,npar); //* Scale of each paramater (output from hesscov)*/
if(mle==-1){ /* Print a wizard for help writing covariance matrix */	if(mle==-1){ /* Print a wizard for help writing covariance matrix */
	/* We could also provide initial parameters values giving by simple logistic regression
	* only one way, that is without matrix product. We will have nlstate maximizations */
	/* for(i=1;i<nlstate;i++){ */
	/* /\reducing xi for 1 to npar to 1 to ncovmodel; \/ */
	/* mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
	/* } */
prwizard(ncovmodel, nlstate, ndeath, model, ficparo);	prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);	printf(" You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);	fprintf(ficlog," You chose mle=-1, look at file %s for a template of covariance matrix \n",filereso);
Line 10009 Please run with mle=-1 to get a correct	Line 10139 Please run with mle=-1 to get a correct
Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */	Typevar=ivector(-1,NCOVMAX); /* -1 to 2 */
Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */	Fixed=ivector(-1,NCOVMAX); /* -1 to 3 */
Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */	Dummy=ivector(-1,NCOVMAX); /* -1 to 3 */
DummyV=ivector(1,NCOVMAX); /* 1 to 3 */
FixedV=ivector(1,NCOVMAX); /* 1 to 3 */
/* V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs).	/* V2+V1+V4+age*V3 is a model with 4 covariates (3 plus signs).
For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4,	For each model-covariate stores the data-covariate id. Tvar[1]=2, Tvar[2]=1, Tvar[3]=4,
Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.	Tvar[4=age*V3] is 3 and 'age' is recorded in Tage.
Line 10537 Please run with mle=-1 to get a correct	Line 10665 Please run with mle=-1 to get a correct
printf("\n");	printf("\n");
if(mle>=1){ /* Could be 1 or 2, Real Maximization */	if(mle>=1){ /* Could be 1 or 2, Real Maximization */
/* mlikeli uses func not funcone */	/* mlikeli uses func not funcone */
	/* for(i=1;i<nlstate;i++){ */
	/* /\reducing xi for 1 to npar to 1 to ncovmodel; \/ */
	/* mlikeli(ficres,p, ncovmodel, ncovmodel, nlstate, ftol, funcnoprod); */
	/* } */
mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);	mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
}	}
if(mle==0) {/* No optimization, will print the likelihoods for the datafile */	if(mle==0) {/* No optimization, will print the likelihoods for the datafile */
Line 10828 Please run with mle=-1 to get a correct	Line 10960 Please run with mle=-1 to get a correct
}else	}else
break;	break;
}	}
	if (!feof(ficpar))
while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){	while((num_filled=sscanf(line,"result:%[^\n]\n",resultline)) !=EOF){
if (num_filled == 0)	if (num_filled == 0){
resultline[0]='\0';	resultline[0]='\0';
else if (num_filled != 1){	break;
	} else if (num_filled != 1){
printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);	printf("ERROR %d: result line should be at minimum 'result=' %s\n",num_filled, line);
}	}
nresult++; /* Sum of resultlines */	nresult++; /* Sum of resultlines */
Line 10861 Please run with mle=-1 to get a correct	Line 10995 Please run with mle=-1 to get a correct
break;	break;
else{ /* Processess output results for this combination of covariate values */	else{ /* Processess output results for this combination of covariate values */
}	}
}	} /* end while */



Line 11079 Please run with mle=-1 to get a correct	Line 11213 Please run with mle=-1 to get a correct
for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */	for(k=1; k<=i1;k++){ /* For any combination of dummy covariates, fixed and varying */
if(TKresult[nres]!= k)	if(TKresult[nres]!= k)
continue;	continue;
printf("\n#****** Selected:");	printf("\n#****** Result for:");
fprintf(ficrest,"\n#****** Selected:");	fprintf(ficrest,"\n#****** Result for:");
fprintf(ficlog,"\n#****** Selected:");	fprintf(ficlog,"\n#****** Result for:");
for(j=1;j<=cptcoveff;j++){	for(j=1;j<=cptcoveff;j++){
printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);	printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);	fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);

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

Removed from v.1.238
changed lines
	Added in v.1.247