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.98 and 1.115

version 1.98, 2004/05/16 15:05:56	version 1.115, 2006/02/27 12:17:45
Line 1	Line 1
/* $Id$	/* $Id$
$State$	$State$
$Log$	$Log$
	Revision 1.115 2006/02/27 12:17:45 brouard
	(Module): One freematrix added in mlikeli! 0.98c

	Revision 1.114 2006/02/26 12:57:58 brouard
	(Module): Some improvements in processing parameter
	filename with strsep.

	Revision 1.113 2006/02/24 14:20:24 brouard
	(Module): Memory leaks checks with valgrind and:
	datafile was not closed, some imatrix were not freed and on matrix
	allocation too.

	Revision 1.112 2006/01/30 09:55:26 brouard
	(Module): Back to gnuplot.exe instead of wgnuplot.exe

	Revision 1.111 2006/01/25 20:38:18 brouard
	(Module): Lots of cleaning and bugs added (Gompertz)
	(Module): Comments can be added in data file. Missing date values
	can be a simple dot '.'.

	Revision 1.110 2006/01/25 00:51:50 brouard
	(Module): Lots of cleaning and bugs added (Gompertz)

	Revision 1.109 2006/01/24 19:37:15 brouard
	(Module): Comments (lines starting with a #) are allowed in data.

	Revision 1.108 2006/01/19 18:05:42 lievre
	Gnuplot problem appeared...
	To be fixed

	Revision 1.107 2006/01/19 16:20:37 brouard
	Test existence of gnuplot in imach path

	Revision 1.106 2006/01/19 13:24:36 brouard
	Some cleaning and links added in html output

	Revision 1.105 2006/01/05 20:23:19 lievre
	* empty log message *

	Revision 1.104 2005/09/30 16:11:43 lievre
	(Module): sump fixed, loop imx fixed, and simplifications.
	(Module): If the status is missing at the last wave but we know
	that the person is alive, then we can code his/her status as -2
	(instead of missing=-1 in earlier versions) and his/her
	contributions to the likelihood is 1 - Prob of dying from last
	health status (= 1-p13= p11+p12 in the easiest case of somebody in
	the healthy state at last known wave). Version is 0.98

	Revision 1.103 2005/09/30 15:54:49 lievre
	(Module): sump fixed, loop imx fixed, and simplifications.

	Revision 1.102 2004/09/15 17:31:30 brouard
	Add the possibility to read data file including tab characters.

	Revision 1.101 2004/09/15 10:38:38 brouard
	Fix on curr_time

	Revision 1.100 2004/07/12 18:29:06 brouard
	Add version for Mac OS X. Just define UNIX in Makefile

	Revision 1.99 2004/06/05 08:57:40 brouard
	* empty log message *

Revision 1.98 2004/05/16 15:05:56 brouard	Revision 1.98 2004/05/16 15:05:56 brouard
New version 0.97 . First attempt to estimate force of mortality	New version 0.97 . First attempt to estimate force of mortality
directly from the data i.e. without the need of knowing the health	directly from the data i.e. without the need of knowing the health
Line 194	Line 257
#include <math.h>	#include <math.h>
#include <stdio.h>	#include <stdio.h>
#include <stdlib.h>	#include <stdlib.h>
	#include <string.h>
#include <unistd.h>	#include <unistd.h>

#include <sys/time.h>	#include <limits.h>
	#include <sys/types.h>
	#include <sys/stat.h>
	#include <errno.h>
	extern int errno;

	/* #include <sys/time.h> */
#include <time.h>	#include <time.h>
#include "timeval.h"	#include "timeval.h"

Line 204	Line 274
/* #define _(String) gettext (String) */	/* #define _(String) gettext (String) */

#define MAXLINE 256	#define MAXLINE 256

#define GNUPLOTPROGRAM "gnuplot"	#define GNUPLOTPROGRAM "gnuplot"
/#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"/	/#define GNUPLOTPROGRAM "..\\gp37mgw\\wgnuplot"/
#define FILENAMELENGTH 132	#define FILENAMELENGTH 132
/#define DEBUG/
/#define windows/
#define GLOCK_ERROR_NOPATH -1 /* empty path */	#define GLOCK_ERROR_NOPATH -1 /* empty path */
#define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */	#define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */

Line 224	Line 294
#define AGESUP 130	#define AGESUP 130
#define AGEBASE 40	#define AGEBASE 40
#define AGEGOMP 10. /* Minimal age for Gompertz adjustment */	#define AGEGOMP 10. /* Minimal age for Gompertz adjustment */
#ifdef unix	#ifdef UNIX
#define DIRSEPARATOR '/'	#define DIRSEPARATOR '/'
	#define CHARSEPARATOR "/"
#define ODIRSEPARATOR '\\'	#define ODIRSEPARATOR '\\'
#else	#else
#define DIRSEPARATOR '\\'	#define DIRSEPARATOR '\\'
	#define CHARSEPARATOR "\\"
#define ODIRSEPARATOR '/'	#define ODIRSEPARATOR '/'
#endif	#endif

/* $Id$ */	/* $Id$ */
/* $State$ */	/* $State$ */

char version[]="Imach version 0.70, May 2004, INED-EUROREVES ";	char version[]="Imach version 0.98c, February 2006, INED-EUROREVES ";
char fullversion[]="$Revision$ $Date$";	char fullversion[]="$Revision$ $Date$";
int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */	int erreur, nberr=0, nbwarn=0; /* Error number, number of errors number of warnings */
int nvar;	int nvar;
Line 249 int popbased=0;	Line 321 int popbased=0;
int wav; / Number of waves for this individuual 0 is possible */	int wav; / Number of waves for this individuual 0 is possible */
int maxwav; /* Maxim number of waves */	int maxwav; /* Maxim number of waves */
int jmin, jmax; /* min, max spacing between 2 waves */	int jmin, jmax; /* min, max spacing between 2 waves */
	int ijmin, ijmax; /* Individuals having jmin and jmax */
int gipmx, gsw; /* Global variables on the number of contributions	int gipmx, gsw; /* Global variables on the number of contributions
to the likelihood and the sum of weights (done by funcone)*/	to the likelihood and the sum of weights (done by funcone)*/
int mle, weightopt;	int mle, weightopt;
Line 279 FILE *ficresvpl;	Line 352 FILE *ficresvpl;
char fileresvpl[FILENAMELENGTH];	char fileresvpl[FILENAMELENGTH];
char title[MAXLINE];	char title[MAXLINE];
char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];	char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH], filerespl[FILENAMELENGTH];
char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH];	char optionfilext[10], optionfilefiname[FILENAMELENGTH], plotcmd[FILENAMELENGTH], pplotcmd[FILENAMELENGTH];
char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH];	char tmpout[FILENAMELENGTH], tmpout2[FILENAMELENGTH];
char command[FILENAMELENGTH];	char command[FILENAMELENGTH];
int outcmd=0;	int outcmd=0;
Line 301 long time_value;	Line 374 long time_value;
extern long time();	extern long time();
char strcurr[80], strfor[80];	char strcurr[80], strfor[80];

	char *endptr;
	long lval;

#define NR_END 1	#define NR_END 1
#define FREE_ARG char*	#define FREE_ARG char*
#define FTOL 1.0e-10	#define FTOL 1.0e-10
Line 349 double *weight;	Line 425 double *weight;
int *s; / Status */	int *s; / Status */
double agedc, *covar, idx;	double agedc, *covar, idx;
int *nbcode, Tcode, Tvar, codtab, Tvard, Tprod, cptcovprod, *Tvaraff;	int *nbcode, Tcode, Tvar, codtab, Tvard, Tprod, cptcovprod, *Tvaraff;
	double lsurv, lpop, *tpop;

double ftol=FTOL; /* Tolerance for computing Max Likelihood */	double ftol=FTOL; /* Tolerance for computing Max Likelihood */
double ftolhess; /* Tolerance for computing hessian */	double ftolhess; /* Tolerance for computing hessian */
Line 356 double ftolhess; /* Tolerance for comput	Line 433 double ftolhess; /* Tolerance for comput
/************** split ***********************/	/************** split ***********************/
static int split( char path, char dirc, char name, char ext, char *finame )	static int split( char path, char dirc, char name, char ext, char *finame )
{	{
	/* From a file name with (full) path (either Unix or Windows) we extract the directory (dirc)
	the name of the file (name), its extension only (ext) and its first part of the name (finame)
	*/
char ss; / pointer */	char ss; / pointer */
int l1, l2; /* length counters */	int l1, l2; /* length counters */

l1 = strlen(path ); /* length of path */	l1 = strlen(path ); /* length of path */
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );	if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
ss= strrchr( path, DIRSEPARATOR ); /* find last / */	ss= strrchr( path, DIRSEPARATOR ); /* find last / */
if ( ss == NULL ) { /* no directory, so use current */	if ( ss == NULL ) { /* no directory, so determine current directory */
	strcpy( name, path ); /* we got the fullname name because no directory */
/*if(strrchr(path, ODIRSEPARATOR )==NULL)	/*if(strrchr(path, ODIRSEPARATOR )==NULL)
printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/	printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
/* get current working directory */	/* get current working directory */
Line 370 static int split( char path, char dirc	Line 451 static int split( char path, char dirc
if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {	if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {
return( GLOCK_ERROR_GETCWD );	return( GLOCK_ERROR_GETCWD );
}	}
strcpy( name, path ); /* we've got it */	/* got dirc from getcwd*/
	printf(" DIRC = %s \n",dirc);
} else { /* strip direcotry from path */	} else { /* strip direcotry from path */
ss++; /* after this, the filename */	ss++; /* after this, the filename */
l2 = strlen( ss ); /* length of filename */	l2 = strlen( ss ); /* length of filename */
Line 378 static int split( char path, char dirc	Line 460 static int split( char path, char dirc
strcpy( name, ss ); /* save file name */	strcpy( name, ss ); /* save file name */
strncpy( dirc, path, l1 - l2 ); /* now the directory */	strncpy( dirc, path, l1 - l2 ); /* now the directory */
dirc[l1-l2] = 0; /* add zero */	dirc[l1-l2] = 0; /* add zero */
	printf(" DIRC2 = %s \n",dirc);
}	}
	/* We add a separator at the end of dirc if not exists */
l1 = strlen( dirc ); /* length of directory */	l1 = strlen( dirc ); /* length of directory */
/*#ifdef windows	if( dirc[l1-1] != DIRSEPARATOR ){
if ( dirc[l1-1] != '\\' ) { dirc[l1] = '\\'; dirc[l1+1] = 0; }	dirc[l1] = DIRSEPARATOR;
#else	dirc[l1+1] = 0;
if ( dirc[l1-1] != '/' ) { dirc[l1] = '/'; dirc[l1+1] = 0; }	printf(" DIRC3 = %s \n",dirc);
#endif	}
*/
ss = strrchr( name, '.' ); /* find last / */	ss = strrchr( name, '.' ); /* find last / */
ss++;	if (ss >0){
strcpy(ext,ss); /* save extension */	ss++;
l1= strlen( name);	strcpy(ext,ss); /* save extension */
l2= strlen(ss)+1;	l1= strlen( name);
strncpy( finame, name, l1-l2);	l2= strlen(ss)+1;
finame[l1-l2]= 0;	strncpy( finame, name, l1-l2);
	finame[l1-l2]= 0;
	}

return( 0 ); /* we're done */	return( 0 ); /* we're done */
}	}

Line 425 int nbocc(char *s, char occ)	Line 511 int nbocc(char *s, char occ)

void cutv(char u,char v, char*t, char occ)	void cutv(char u,char v, char*t, char occ)
{	{
/* cuts string t into u and v where u is ended by char occ excluding it	/* cuts string t into u and v where u ends before first occurence of char 'occ'
and v is after occ excluding it too : ex cutv(u,v,"abcdef2ghi2j",2)	and v starts after first occurence of char 'occ' : ex cutv(u,v,"abcdef2ghi2j",'2')
gives u="abcedf" and v="ghi2j" */	gives u="abcedf" and v="ghi2j" */
int i,lg,j,p=0;	int i,lg,j,p=0;
i=0;	i=0;
Line 868 void powell(double p[], double **xi, int	Line 954 void powell(double p[], double **xi, int
fprintf(ficrespow,"\n");fflush(ficrespow);	fprintf(ficrespow,"\n");fflush(ficrespow);
if(*iter <=3){	if(*iter <=3){
tm = *localtime(&curr_time.tv_sec);	tm = *localtime(&curr_time.tv_sec);
strcpy(strcurr,asctime(&tmf));	strcpy(strcurr,asctime(&tm));
/* asctime_r(&tm,strcurr); */	/* asctime_r(&tm,strcurr); */
forecast_time=curr_time;	forecast_time=curr_time;
itmp = strlen(strcurr);	itmp = strlen(strcurr);
if(strcurr[itmp-1]=='\n')	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 this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);	printf("\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);	fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);
Line 884 void powell(double p[], double **xi, int	Line 970 void powell(double p[], double **xi, int
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 or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,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(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
fprintf(ficlog," - if your program needs %d iterations to converge, convergence will be \n reached in %s or\n on %s (current time is %s);\n",niterf, asc_diff_time(forecast_time.tv_sec-curr_time.tv_sec,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(forecast_time.tv_sec-curr_time.tv_sec,tmpout),strfor,strcurr);
}	}
}	}
for (i=1;i<=n;i++) {	for (i=1;i<=n;i++) {
Line 1050 double pmij(double ps, double *cov,	Line 1136 double pmij(double ps, double *cov,
int i,j,j1, nc, ii, jj;	int i,j,j1, nc, ii, jj;

for(i=1; i<= nlstate; i++){	for(i=1; i<= nlstate; i++){
for(j=1; j<i;j++){	for(j=1; j<i;j++){
for (nc=1, s2=0.;nc <=ncovmodel; nc++){	for (nc=1, s2=0.;nc <=ncovmodel; nc++){
/s2 += param[i][j][nc]cov[nc];*/	/s2 += param[i][j][nc]cov[nc];*/
s2 += x[(i-1)nlstatencovmodel+(j-1)ncovmodel+nc+(i-1)(ndeath-1)ncovmodel]cov[nc];	s2 += x[(i-1)nlstatencovmodel+(j-1)ncovmodel+nc+(i-1)(ndeath-1)ncovmodel]cov[nc];
/printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2);/	/* printf("Int j<i s1=%.17e, s2=%.17e\n",s1,s2); */
}	}
ps[i][j]=s2;	ps[i][j]=s2;
/printf("s1=%.17e, s2=%.17e\n",s1,s2);/	/* printf("s1=%.17e, s2=%.17e\n",s1,s2); */
}	}
for(j=i+1; j<=nlstate+ndeath;j++){	for(j=i+1; j<=nlstate+ndeath;j++){
for (nc=1, s2=0.;nc <=ncovmodel; nc++){	for (nc=1, s2=0.;nc <=ncovmodel; nc++){
s2 += x[(i-1)nlstatencovmodel+(j-2)ncovmodel+nc+(i-1)(ndeath-1)ncovmodel]cov[nc];	s2 += x[(i-1)nlstatencovmodel+(j-2)ncovmodel+nc+(i-1)(ndeath-1)ncovmodel]cov[nc];
/printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2);/	/* printf("Int j>i s1=%.17e, s2=%.17e %lx %lx\n",s1,s2,s1,s2); */
	}
	ps[i][j]=s2;
}	}
ps[i][j]=s2;
}	}
}
/ps[3][2]=1;/	/ps[3][2]=1;/

for(i=1; i<= nlstate; i++){	for(i=1; i<= nlstate; i++){
s1=0;	s1=0;
for(j=1; j<i; j++)	for(j=1; j<i; j++)
s1+=exp(ps[i][j]);	s1+=exp(ps[i][j]);
for(j=i+1; j<=nlstate+ndeath; j++)	for(j=i+1; j<=nlstate+ndeath; j++)
s1+=exp(ps[i][j]);	s1+=exp(ps[i][j]);
ps[i][i]=1./(s1+1.);	ps[i][i]=1./(s1+1.);
for(j=1; j<i; j++)	for(j=1; j<i; j++)
ps[i][j]= exp(ps[i][j])*ps[i][i];	ps[i][j]= exp(ps[i][j])*ps[i][i];
for(j=i+1; j<=nlstate+ndeath; j++)	for(j=i+1; j<=nlstate+ndeath; j++)
ps[i][j]= exp(ps[i][j])*ps[i][i];	ps[i][j]= exp(ps[i][j])*ps[i][i];
/* ps[i][nlstate+1]=1.-s1- ps[i][i];/ / Sum should be 1 */	/* ps[i][nlstate+1]=1.-s1- ps[i][i];/ / Sum should be 1 */
} /* end i */	} /* end i */

for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){	for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
for(jj=1; jj<= nlstate+ndeath; jj++){	for(jj=1; jj<= nlstate+ndeath; jj++){
ps[ii][jj]=0;	ps[ii][jj]=0;
ps[ii][ii]=1;	ps[ii][ii]=1;
	}
}	}
}


/* for(ii=1; ii<= nlstate+ndeath; ii++){	/* for(ii=1; ii<= nlstate+ndeath; ii++){ */
for(jj=1; jj<= nlstate+ndeath; jj++){	/* for(jj=1; jj<= nlstate+ndeath; jj++){ */
printf("%lf ",ps[ii][jj]);	/* printf("ddd %lf ",ps[ii][jj]); */
}	/* } */
printf("\n ");	/* printf("\n "); */
}	/* } */
printf("\n ");printf("%lf ",cov[2]);*/	/* printf("\n ");printf("%lf ",cov[2]); */
/*	/*
for(i=1; i<= npar; i++) printf("%f ",x[i]);	for(i=1; i<= npar; i++) printf("%f ",x[i]);
goto end;*/	goto end;*/
return ps;	return ps;
}	}

Line 1224 double func( double *x)	Line 1310 double func( double *x)
} /* end mult */	} /* end mult */

/lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);/ /* Original formula */	/lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]);/ /* Original formula */
/* But now since version 0.9 we anticipate for bias and large stepm.	/* But now since version 0.9 we anticipate for bias at large stepm.
* If stepm is larger than one month (smallest stepm) and if the exact delay	* If stepm is larger than one month (smallest stepm) and if the exact delay
* (in months) between two waves is not a multiple of stepm, we rounded to	* (in months) between two waves is not a multiple of stepm, we rounded to
* the nearest (and in case of equal distance, to the lowest) interval but now	* the nearest (and in case of equal distance, to the lowest) interval but now
* we keep into memory the bias bh[mi][i] and also the previous matrix product	* we keep into memory the bias bh[mi][i] and also the previous matrix product
* (i.e to dh[mi][i]-1) saved in 'savm'. The we inter(extra)polate the	* (i.e to dh[mi][i]-1) saved in 'savm'. Then we inter(extra)polate the
* probability in order to take into account the bias as a fraction of the way	* probability in order to take into account the bias as a fraction of the way
* from savm to out if bh is neagtive or even beyond if bh is positive. bh varies	* from savm to out if bh is negative or even beyond if bh is positive. bh varies
* -stepm/2 to stepm/2 .	* -stepm/2 to stepm/2 .
* For stepm=1 the results are the same as for previous versions of Imach.	* For stepm=1 the results are the same as for previous versions of Imach.
* For stepm > 1 the results are less biased than in previous versions.	* For stepm > 1 the results are less biased than in previous versions.
Line 1239 double func( double *x)	Line 1325 double func( double *x)
s1=s[mw[mi][i]][i];	s1=s[mw[mi][i]][i];
s2=s[mw[mi+1][i]][i];	s2=s[mw[mi+1][i]][i];
bbh=(double)bh[mi][i]/(double)stepm;	bbh=(double)bh[mi][i]/(double)stepm;
/* bias is positive if real duration	/* bias bh is positive if real duration
* is higher than the multiple of stepm and negative otherwise.	* is higher than the multiple of stepm and negative otherwise.
*/	*/
/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.+bbh)out[s1][s2]));/	/* lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)log(out[s1][s2])- bbhlog(savm[s1][s2]):log((1.+bbh)out[s1][s2]));/
if( s2 > nlstate){	if( s2 > nlstate){
/* i.e. if s2 is a death state and if the date of death is known then the contribution	/* i.e. if s2 is a death state and if the date of death is known
to the likelihood is the probability to die between last step unit time and current	then the contribution to the likelihood is the probability to
step unit time, which is also the differences between probability to die before dh	die between last step unit time and current step unit time,
and probability to die before dh-stepm .	which is also equal to probability to die before dh
	minus probability to die before dh-stepm .
In version up to 0.92 likelihood was computed	In version up to 0.92 likelihood was computed
as if date of death was unknown. Death was treated as any other	as if date of death was unknown. Death was treated as any other
health state: the date of the interview describes the actual state	health state: the date of the interview describes the actual state
Line 1267 double func( double *x)	Line 1354 double func( double *x)
lower mortality.	lower mortality.
*/	*/
lli=log(out[s1][s2] - savm[s1][s2]);	lli=log(out[s1][s2] - savm[s1][s2]);
}else{

	} else if (s2==-2) {
	for (j=1,survp=0. ; j<=nlstate; j++)
	survp += out[s1][j];
	lli= survp;
	}

	else if (s2==-4) {
	for (j=3,survp=0. ; j<=nlstate; j++)
	survp += out[s1][j];
	lli= survp;
	}

	else if (s2==-5) {
	for (j=1,survp=0. ; j<=2; j++)
	survp += out[s1][j];
	lli= survp;
	}


	else{
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= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)out[s1][s2]- bbh(savm[s1][s2])):log((1.+bbh)out[s1][s2]));/ /* linear interpolation */	/* lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)out[s1][s2]- bbh(savm[s1][s2])):log((1.+bbh)out[s1][s2]));/ /* linear interpolation */
}	}
Line 1563 void mlikeli(FILE *ficres,double p[], in	Line 1671 void mlikeli(FILE *ficres,double p[], in

powell(p,xi,npar,ftol,&iter,&fret,func);	powell(p,xi,npar,ftol,&iter,&fret,func);

	free_matrix(xi,1,npar,1,npar);
fclose(ficrespow);	fclose(ficrespow);
printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));	printf("\n#Number of iterations = %d, -2 Log likelihood = %.12f\n",iter,func(p));
fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));	fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));
Line 1831 void lubksb(double *a, int n, int indx	Line 1940 void lubksb(double *a, int n, int indx
}	}

/********** Frequencies ******************/	/********** Frequencies ******************/
void freqsummary(char fileres[], int iagemin, int iagemax, int s, double agev, int nlstate, int imx, int Tvaraff, int nbcode, int ncodemax,double mint,double anint)	void freqsummary(char fileres[], int iagemin, int iagemax, int s, double agev, int nlstate, int imx, int Tvaraff, int nbcode, int ncodemax,double mint,double anint, char strstart[])
{ /* Some frequencies */	{ /* Some frequencies */

int i, m, jk, k1,i1, j1, bool, z1,z2,j;	int i, m, jk, k1,i1, j1, bool, z1,z2,j;
Line 1851 void freqsummary(char fileres[], int ia	Line 1960 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);
}	}
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,iagemin,iagemax+3);	freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin,iagemax+3);
j1=0;	j1=0;

j=cptcoveff;	j=cptcoveff;
Line 1864 void freqsummary(char fileres[], int ia	Line 1973 void freqsummary(char fileres[], int ia
j1++;	j1++;
/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);	/*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
scanf("%d", i);*/	scanf("%d", i);*/
for (i=-1; i<=nlstate+ndeath; i++)	for (i=-5; i<=nlstate+ndeath; i++)
for (jk=-1; 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;

Line 1904 void freqsummary(char fileres[], int ia	Line 2013 void freqsummary(char fileres[], int ia
}	}

/* 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);*/
	fprintf(ficresp, "#Local time at start: %s", strstart);
if (cptcovn>0) {	if (cptcovn>0) {
fprintf(ficresp, "\n#********** Variable ");	fprintf(ficresp, "\n#********** Variable ");
for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);	for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
Line 1990 void freqsummary(char fileres[], int ia	Line 2099 void freqsummary(char fileres[], int ia
dateintmean=dateintsum/k2cpt;	dateintmean=dateintsum/k2cpt;

fclose(ficresp);	fclose(ficresp);
free_ma3x(freq,-1,nlstate+ndeath,-1,nlstate+ndeath, iagemin, iagemax+3);	free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin, iagemax+3);
free_vector(pp,1,nlstate);	free_vector(pp,1,nlstate);
free_matrix(prop,1,nlstate,iagemin, iagemax+3);	free_matrix(prop,1,nlstate,iagemin, iagemax+3);
/* End of Freq */	/* End of Freq */
Line 2099 void concatwav(int wav[], int **dh, int	Line 2208 void concatwav(int wav[], int **dh, int
mi=0;	mi=0;
m=firstpass;	m=firstpass;
while(s[m][i] <= nlstate){	while(s[m][i] <= nlstate){
if(s[m][i]>=1)	if(s[m][i]>=1 \|\| s[m][i]==-2 \|\| s[m][i]==-4 \|\| s[m][i]==-5)
mw[++mi][i]=m;	mw[++mi][i]=m;
if(m >=lastpass)	if(m >=lastpass)
break;	break;
Line 2117 void concatwav(int wav[], int **dh, int	Line 2226 void concatwav(int wav[], int **dh, int
if(mi==0){	if(mi==0){
nbwarn++;	nbwarn++;
if(first==0){	if(first==0){
printf("Warning! None valid information for:%ld line=%d (skipped) and may be others, see log file\n",num[i],i);	printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
first=1;	first=1;
}	}
if(first==1){	if(first==1){
fprintf(ficlog,"Warning! None valid information for:%ld line=%d (skipped)\n",num[i],i);	fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
}	}
} /* end mi==0 */	} /* end mi==0 */
} /* End individuals */	} /* End individuals */
Line 2139 void concatwav(int wav[], int **dh, int	Line 2248 void concatwav(int wav[], int **dh, int
nberr++;	nberr++;
printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);	printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
j=1; /* Temporary Dangerous patch */	j=1; /* Temporary Dangerous patch */
printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n You MUST fix the contradiction between dates.\n",stepm);	printf(" We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);	fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview.\n You MUST fix the contradiction between dates.\n",stepm);	fprintf(ficlog," We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
}	}
k=k+1;	k=k+1;
if (j >= jmax) jmax=j;	if (j >= jmax){
if (j <= jmin) jmin=j;	jmax=j;
	ijmax=i;
	}
	if (j <= jmin){
	jmin=j;
	ijmin=i;
	}
sum=sum+j;	sum=sum+j;
/if (j<0) printf("j=%d num=%d \n",j,i);/	/if (j<0) printf("j=%d num=%d \n",j,i);/
/* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/	/* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
Line 2153 void concatwav(int wav[], int **dh, int	Line 2268 void concatwav(int wav[], int **dh, int
}	}
else{	else{
j= rint( (agev[mw[mi+1][i]][i]12 - agev[mw[mi][i]][i]12));	j= rint( (agev[mw[mi+1][i]][i]12 - agev[mw[mi][i]][i]12));
/* printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/	/* if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */

k=k+1;	k=k+1;
if (j >= jmax) jmax=j;	if (j >= jmax) {
else if (j <= jmin)jmin=j;	jmax=j;
	ijmax=i;
	}
	else if (j <= jmin){
	jmin=j;
	ijmin=i;
	}
/* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */	/* if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
/printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);/	/printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);/
if(j<0){	if(j<0){
Line 2199 void concatwav(int wav[], int **dh, int	Line 2321 void concatwav(int wav[], int **dh, int
} /* end wave */	} /* end wave */
}	}
jmean=sum/k;	jmean=sum/k;
printf("Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);	printf("Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, num[ijmin], jmax, num[ijmax], jmean);
fprintf(ficlog,"Delay (in months) between two waves Min=%d Max=%d Mean=%f\n\n ",jmin, jmax,jmean);	fprintf(ficlog,"Delay (in months) between two waves Min=%d (for indiviudal %ld) Max=%d (%ld) Mean=%f\n\n ",jmin, ijmin, jmax, ijmax, jmean);
}	}

/********* Tricode **************************/	/********* Tricode **************************/
Line 2246 void tricode(int Tvar, int *nbcode, in	Line 2368 void tricode(int Tvar, int *nbcode, in
for (k=0; k< maxncov; k++) Ndum[k]=0;	for (k=0; k< maxncov; k++) Ndum[k]=0;

for (i=1; i<=ncovmodel-2; i++) {	for (i=1; i<=ncovmodel-2; i++) {
/* Listing of all covariables in staement 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];	ij=Tvar[i];
Ndum[ij]++;	Ndum[ij]++;
}	}
Line 2264 void tricode(int Tvar, int *nbcode, in	Line 2386 void tricode(int Tvar, int *nbcode, in

/********* Health Expectancies **************/	/********* Health Expectancies **************/

void evsij(char fileres[], double *eij, double x[], int nlstate, int stepm, int bage, int fage, double oldm, double savm, int ij, int estepm,double delti[],double matcov )	void evsij(char fileres[], double *eij, double x[], int nlstate, int stepm, int bage, int fage, double oldm, double savm, int ij, int estepm,double delti[],double matcov,char strstart[] )

{	{
/* Health expectancies */	/* Health expectancies */
Line 2282 void evsij(char fileres[], double ***eij	Line 2404 void evsij(char fileres[], double ***eij
dnewm=matrix(1,nlstate*nlstate,1,npar);	dnewm=matrix(1,nlstate*nlstate,1,npar);
doldm=matrix(1,nlstatenlstate,1,nlstatenlstate);	doldm=matrix(1,nlstatenlstate,1,nlstatenlstate);

	fprintf(ficreseij,"# Local time at start: %s", strstart);
fprintf(ficreseij,"# Health expectancies\n");	fprintf(ficreseij,"# Health expectancies\n");
fprintf(ficreseij,"# Age");	fprintf(ficreseij,"# Age");
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
Line 2436 void evsij(char fileres[], double ***eij	Line 2559 void evsij(char fileres[], double ***eij
}	}

/********** Variance ****************/	/********** Variance ****************/
void varevsij(char optionfilefiname[], double *vareij, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav)	void varevsij(char optionfilefiname[], double *vareij, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double **prlim, double ftolpl, int ij, int estepm, int cptcov, int cptcod, int popbased, int mobilav, char strstart[])
{	{
/* Variance of health expectancies */	/* Variance of health expectancies */
/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double * savm,double ftolpl);*/	/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double * savm,double ftolpl);*/
Line 2487 void varevsij(char optionfilefiname[], d	Line 2610 void varevsij(char optionfilefiname[], d
fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);	fprintf(ficlog,"Problem with resultfile: %s\n", fileresprobmorprev);
}	}
printf("Computing total mortality p.j=w1p1j+w2p2j+..: result on file '%s' \n",fileresprobmorprev);	printf("Computing total mortality p.j=w1p1j+w2p2j+..: result on file '%s' \n",fileresprobmorprev);

fprintf(ficlog,"Computing total mortality p.j=w1p1j+w2p2j+..: result on file '%s' \n",fileresprobmorprev);	fprintf(ficlog,"Computing total mortality p.j=w1p1j+w2p2j+..: result on file '%s' \n",fileresprobmorprev);
	fprintf(ficresprobmorprev, "#Local time at start: %s", strstart);
fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1p1j+w2p2j+... stand dev in()\n",estepm);	fprintf(ficresprobmorprev,"# probabilities of dying before estepm=%d months for people of exact age and weighted probabilities w1p1j+w2p2j+... stand dev in()\n",estepm);
fprintf(ficresprobmorprev,"# Age cov=%-d",ij);	fprintf(ficresprobmorprev,"# Age cov=%-d",ij);
for(j=nlstate+1; j<=(nlstate+ndeath);j++){	for(j=nlstate+1; j<=(nlstate+ndeath);j++){
Line 2497 void varevsij(char optionfilefiname[], d	Line 2622 void varevsij(char optionfilefiname[], d
}	}
fprintf(ficresprobmorprev,"\n");	fprintf(ficresprobmorprev,"\n");
fprintf(ficgp,"\n# Routine varevsij");	fprintf(ficgp,"\n# Routine varevsij");
	/* fprintf(fichtm, "#Local time at start: %s", strstart);*/
fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");	fprintf(fichtm,"\n<li><h4> Computing probabilities of dying over estepm months as a weighted average (i.e global mortality independent of initial healh state)</h4></li>\n");
fprintf(fichtm,"\n<br>%s <br>\n",digitp);	fprintf(fichtm,"\n<br>%s <br>\n",digitp);
/* } */	/* } */
varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);	varppt = matrix(nlstate+1,nlstate+ndeath,nlstate+1,nlstate+ndeath);
	fprintf(ficresvij, "#Local time at start: %s", strstart);
fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n");	fprintf(ficresvij,"# Variance and covariance of health expectancies e.j \n# (weighted average of eij where weights are the stable prevalence in health states i\n");
fprintf(ficresvij,"# Age");	fprintf(ficresvij,"# Age");
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
Line 2734 void varevsij(char optionfilefiname[], d	Line 2860 void varevsij(char optionfilefiname[], d
} /* end varevsij */	} /* end varevsij */

/********** Variance of prevlim ****************/	/********** Variance of prevlim ****************/
void varprevlim(char fileres[], double varpl, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double **prlim, double ftolpl, int ij)	void varprevlim(char fileres[], double varpl, double matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double oldm, double savm, double **prlim, double ftolpl, int ij, char strstart[])
{	{
/* Variance of prevalence limit */	/* Variance of prevalence limit */
/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double savm,double ftolpl);/	/* double prevalim(double prlim, int nlstate, double xp, double age, double oldm, double savm,double ftolpl);/
Line 2747 void varprevlim(char fileres[], double *	Line 2873 void varprevlim(char fileres[], double *
double gradg, trgradg;	double gradg, trgradg;
double age,agelim;	double age,agelim;
int theta;	int theta;
	fprintf(ficresvpl, "#Local time at start: %s", strstart);
fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");	fprintf(ficresvpl,"# Standard deviation of stable prevalences \n");
fprintf(ficresvpl,"# Age");	fprintf(ficresvpl,"# Age");
for(i=1; i<=nlstate;i++)	for(i=1; i<=nlstate;i++)
Line 2817 void varprevlim(char fileres[], double *	Line 2943 void varprevlim(char fileres[], double *
}	}

/********** Variance of one-step probabilities ****************/	/********** Variance of one-step probabilities ****************/
void varprob(char optionfilefiname[], double *matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int Tvar, int *nbcode, int ncodemax)	void varprob(char optionfilefiname[], double *matcov, double x[], double delti[], int nlstate, double bage, double fage, int ij, int Tvar, int *nbcode, int ncodemax, char strstart[])
{	{
int i, j=0, i1, k1, l1, t, tj;	int i, j=0, i1, k1, l1, t, tj;
int k2, l2, j1, z1;	int k2, l2, j1, z1;
Line 2862 void varprob(char optionfilefiname[], do	Line 2988 void varprob(char optionfilefiname[], do
fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);	fprintf(ficlog,"Computing matrix of variance covariance of one-step probabilities: result on file '%s' \n",fileresprobcov);
printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);	printf("and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);	fprintf(ficlog,"and correlation matrix of one-step probabilities: result on file '%s' \n",fileresprobcor);
	fprintf(ficresprob, "#Local time at start: %s", strstart);
fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");	fprintf(ficresprob,"#One-step probabilities and stand. devi in ()\n");
fprintf(ficresprob,"# Age");	fprintf(ficresprob,"# Age");
	fprintf(ficresprobcov, "#Local time at start: %s", strstart);
fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");	fprintf(ficresprobcov,"#One-step probabilities and covariance matrix\n");
fprintf(ficresprobcov,"# Age");	fprintf(ficresprobcov,"# Age");
	fprintf(ficresprobcor, "#Local time at start: %s", strstart);
fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");	fprintf(ficresprobcor,"#One-step probabilities and correlation matrix\n");
fprintf(ficresprobcov,"# Age");	fprintf(ficresprobcov,"# Age");

Line 3119 To be simple, these graphs help to under	Line 3247 To be simple, these graphs help to under
}	}
free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);	free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);	free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
	free_matrix(doldm,1,(nlstate)(nlstate+ndeath),1,(nlstate)(nlstate+ndeath));
	free_matrix(dnewm,1,(nlstate)*(nlstate+ndeath),1,npar);
free_vector(xp,1,npar);	free_vector(xp,1,npar);
fclose(ficresprob);	fclose(ficresprob);
fclose(ficresprobcov);	fclose(ficresprobcov);
Line 3137 void printinghtml(char fileres[], char t	Line 3267 void printinghtml(char fileres[], char t
double jprev2, double mprev2,double anprev2){	double jprev2, double mprev2,double anprev2){
int jj1, k1, i1, cpt;	int jj1, k1, i1, cpt;

fprintf(fichtm,"<ul><li><h4>Result files (first order: no variance)</h4>\n \	fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
	<li><a href='#secondorder'>Result files (second order (variance)</a>\n \
	</ul>");
	fprintf(fichtm,"<ul><li><h4><a name='firstorder'>Result files (first order: no variance)</a></h4>\n \
- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",	- Observed prevalence in each state (during the period defined between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf): <a href=\"%s\">%s</a> <br>\n ",
jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));	jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,subdirf2(fileres,"p"),subdirf2(fileres,"p"));
fprintf(fichtm,"\	fprintf(fichtm,"\
Line 3188 fprintf(fichtm," \n<ul><li><b>Graphs</b>	Line 3321 fprintf(fichtm," \n<ul><li><b>Graphs</b>


fprintf(fichtm,"\	fprintf(fichtm,"\
\n<br><li><h4> Result files (second order: variances)</h4>\n\	\n<br><li><h4> <a name='secondorder'>Result files (second order: variances)</a></h4>\n\
- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);	- Parameter file with estimated parameters and covariance matrix: <a href=\"%s\">%s</a> <br>\n", rfileres,rfileres);

fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",	fprintf(fichtm," - Variance of one-step probabilities: <a href=\"%s\">%s</a> <br>\n",
Line 3902 double gompertz(double x[])	Line 4035 double gompertz(double x[])
{	{
double A,B,L=0.0,sump=0.,num=0.;	double A,B,L=0.0,sump=0.,num=0.;
int i,n=0; /* n is the size of the sample */	int i,n=0; /* n is the size of the sample */

for (i=0;i<=imx-1 ; i++) {	for (i=0;i<=imx-1 ; i++) {
sump=sump+weight[i];	sump=sump+weight[i];
sump=sump+1;	/* sump=sump+1;*/
num=num+1;	num=num+1;
}	}


/* for (i=1; i<=imx; i++)	/* for (i=0; i<=imx; i++)
if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/	if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/

for (i=0;i<=imx-1 ; i++)	for (i=1;i<=imx ; i++)
{	{
if (cens[i]==1 & wav[i]>1)	if (cens[i] == 1 && wav[i]>1)
A=-x[1]/(x[2])*	A=-x[1]/(x[2])(exp(x[2](agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)));
(exp(x[2]/YEARM(agecens[i]12-agegomp12))-exp(x[2]/YEARM(ageexmed[i]12-agegomp12)));

if (cens[i]==0 & wav[i]>1)	if (cens[i] == 0 && wav[i]>1)
A=-x[1]/(x[2])*	A=-x[1]/(x[2])(exp(x[2](agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp)))
(exp(x[2]/YEARM(agedc[i]12-agegomp12))-exp(x[2]/YEARM(ageexmed[i]12-agegomp12)))	+log(x[1]/YEARM)+x[2]*(agedc[i]-agegomp)+log(YEARM);
+log(x[1]/YEARM)+x[2]/YEARM(agedc[i]12-agegomp*12)+log(YEARM);

if (wav[i]>1 & agecens[i]>15) {	/if (wav[i] > 1 && agecens[i] > 15) {/ /* ??? */
	if (wav[i] > 1 ) { /* ??? */
L=L+A*weight[i];	L=L+A*weight[i];
/* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]12,agecens[i]12,cens[i],agedc[i]12,weight[i]);/	/* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]12,agecens[i]12,cens[i],agedc[i]12,weight[i]);/
}	}
Line 3937 double gompertz(double x[])	Line 4070 double gompertz(double x[])
/***************** Printing html file *********/	/***************** Printing html file *********/
void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \	void printinghtmlmort(char fileres[], char title[], char datafile[], int firstpass, \
int lastpass, int stepm, int weightopt, char model[],\	int lastpass, int stepm, int weightopt, char model[],\
int imx, double p[],double **matcov){	int imx, double p[],double **matcov,double agemortsup){
int i;	int i,k;

fprintf(fichtm,"<ul><li><h4>Result files </h4>\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):<br>");	fprintf(fichtm,"<ul><li><h4>Result files </h4>\n Force of mortality. Parameters of the Gompertz fit (with confidence interval in brackets):<br>");
fprintf(fichtm," mu(age) =%lfexp(%lf(age-%d)) per year<br><br>",p[1],p[2],agegomp);	fprintf(fichtm," mu(age) =%lfexp(%lf(age-%d)) per year<br><br>",p[1],p[2],agegomp);
Line 3946 void printinghtmlmort(char fileres[], ch	Line 4079 void printinghtmlmort(char fileres[], ch
fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2sqrt(matcov[i][i]),p[i]+2sqrt(matcov[i][i]));	fprintf(fichtm," p[%d] = %lf [%f ; %f]<br>\n",i,p[i],p[i]-2sqrt(matcov[i][i]),p[i]+2sqrt(matcov[i][i]));
fprintf(fichtm,"<br><br><img src=\"graphmort.png\">");	fprintf(fichtm,"<br><br><img src=\"graphmort.png\">");
fprintf(fichtm,"</ul>");	fprintf(fichtm,"</ul>");

	fprintf(fichtm,"<ul><li><h4>Life table</h4>\n <br>");

	fprintf(fichtm,"\nAge l<inf>x</inf> q<inf>x</inf> d(x,x+1) L<inf>x</inf> T<inf>x</inf> e<infx</inf><br>");

	for (k=agegomp;k<(agemortsup-2);k++)
	fprintf(fichtm,"%d %.0lf %lf %.0lf %.0lf %.0lf %lf<br>\n",k,lsurv[k],p[1]exp(p[2](k-agegomp)),(p[1]exp(p[2](k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);


fflush(fichtm);	fflush(fichtm);
}	}

Line 3983 int main(int argc, char *argv[])	Line 4125 int main(int argc, char *argv[])
{	{
int movingaverage(double *probs, double bage,double fage, double *mobaverage, int mobilav);	int movingaverage(double *probs, double bage,double fage, double *mobaverage, int mobilav);
int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;	int i,j, k, n=MAXN,iter,m,size=100,cptcode, cptcod;
	int linei, month, year,iout;
int jj, ll, li, lj, lk, imk;	int jj, ll, li, lj, lk, imk;
int numlinepar=0; /* Current linenumber of parameter file */	int numlinepar=0; /* Current linenumber of parameter file */
int itimes;	int itimes;
int NDIM=2;	int NDIM=2;

char ca[32], cb[32], cc[32];	char ca[32], cb[32], cc[32];
	char dummy[]=" ";
/* FILE fichtm; //* Html File */	/* FILE fichtm; //* Html File */
/* FILE ficgp;/ /Gnuplot File /	/* FILE ficgp;/ /Gnuplot File /
	struct stat info;
double agedeb, agefin,hf;	double agedeb, agefin,hf;
double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;	double ageminpar=1.e20,agemin=1.e20, agemaxpar=-1.e20, agemax=-1.e20;

Line 4003 int main(int argc, char *argv[])	Line 4148 int main(int argc, char *argv[])
int *indx;	int *indx;
char line[MAXLINE], linepar[MAXLINE];	char line[MAXLINE], linepar[MAXLINE];
char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];	char path[MAXLINE],pathc[MAXLINE],pathcd[MAXLINE],pathtot[MAXLINE],model[MAXLINE];
char pathr[MAXLINE];	char pathr[MAXLINE], pathimach[MAXLINE];
	char *bp, tok, val; / pathtot */
int firstobs=1, lastobs=10;	int firstobs=1, lastobs=10;
int sdeb, sfin; /* Status at beginning and end */	int sdeb, sfin; /* Status at beginning and end */
int c, h , cpt,l;	int c, h , cpt,l;
Line 4013 int main(int argc, char *argv[])	Line 4159 int main(int argc, char *argv[])
int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */	int mobilavproj=0 , prevfcast=0 ; /* moving average of prev, If prevfcast=1 prevalence projection */
int mobilav=0,popforecast=0;	int mobilav=0,popforecast=0;
int hstepm, nhstepm;	int hstepm, nhstepm;
	int agemortsup;
	float sumlpop=0.;
double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;	double jprev1=1, mprev1=1,anprev1=2000,jprev2=1, mprev2=1,anprev2=2000;
double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;	double jpyram=1, mpyram=1,anpyram=2000,jpyram1=1, mpyram1=1,anpyram1=2000;

Line 4079 int main(int argc, char *argv[])	Line 4227 int main(int argc, char *argv[])
printf("\n%s\n%s",version,fullversion);	printf("\n%s\n%s",version,fullversion);
if(argc <=1){	if(argc <=1){
printf("\nEnter the parameter file name: ");	printf("\nEnter the parameter file name: ");
scanf("%s",pathtot);	fgets(pathr,FILENAMELENGTH,stdin);
	i=strlen(pathr);
	if(pathr[i-1]=='\n')
	pathr[i-1]='\0';
	for (tok = pathr; tok != NULL; ){
	printf("Pathr \|%s\|\n",pathr);
	while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
	printf("val= \|%s\| pathr=%s\n",val,pathr);
	strcpy (pathtot, val);
	if(pathr[0] == '\0') break; /* Un peu sale */
	}
}	}
else{	else{
strcpy(pathtot,argv[1]);	strcpy(pathtot,argv[1]);
Line 4089 int main(int argc, char *argv[])	Line 4247 int main(int argc, char *argv[])
printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/	printf("pathtot=%s, path=%s, optionfile=%s\n",pathtot,path,optionfile);*/
/* cutv(path,optionfile,pathtot,'\\');*/	/* cutv(path,optionfile,pathtot,'\\');*/

	/* Split argv[0], imach program to get pathimach */
	printf("\nargv[0]=%s argv[1]=%s, \n",argv[0],argv[1]);
	split(argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
	printf("\nargv[0]=%s pathimach=%s, \noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",argv[0],pathimach,optionfile,optionfilext,optionfilefiname);
	/* strcpy(pathimach,argv[0]); */
	/* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */
split(pathtot,path,optionfile,optionfilext,optionfilefiname);	split(pathtot,path,optionfile,optionfilext,optionfilefiname);
printf("pathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);	printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
chdir(path);	chdir(path);
strcpy(command,"mkdir ");	strcpy(command,"mkdir ");
strcat(command,optionfilefiname);	strcat(command,optionfilefiname);
Line 4115 int main(int argc, char *argv[])	Line 4279 int main(int argc, char *argv[])
}	}
fprintf(ficlog,"Log filename:%s\n",filelog);	fprintf(ficlog,"Log filename:%s\n",filelog);
fprintf(ficlog,"\n%s\n%s",version,fullversion);	fprintf(ficlog,"\n%s\n%s",version,fullversion);
fprintf(ficlog,"\nEnter the parameter file name: ");	fprintf(ficlog,"\nEnter the parameter file name: \n");
fprintf(ficlog,"pathtot=%s\n\	fprintf(ficlog,"pathimach=%s\npathtot=%s\n\
path=%s \n\	path=%s \n\
optionfile=%s\n\	optionfile=%s\n\
optionfilext=%s\n\	optionfilext=%s\n\
optionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);	optionfilefiname=%s\n",pathimach,pathtot,path,optionfile,optionfilext,optionfilefiname);

printf("Local time (at start):%s",strstart);	printf("Local time (at start):%s",strstart);
fprintf(ficlog,"Local time (at start): %s",strstart);	fprintf(ficlog,"Local time (at start): %s",strstart);
Line 4257 int main(int argc, char *argv[])	Line 4421 int main(int argc, char *argv[])
}	}
fflush(ficlog);	fflush(ficlog);


p=param[1][1];	p=param[1][1];

/* Reads comments: lines beginning with '#' */	/* Reads comments: lines beginning with '#' */
Line 4381 int main(int argc, char *argv[])	Line 4544 int main(int argc, char *argv[])
ncodemax=ivector(1,8);	ncodemax=ivector(1,8);

i=1;	i=1;
while (fgets(line, MAXLINE, fic) != NULL) {	linei=0;
if ((i >= firstobs) && (i <=lastobs)) {	while ((fgets(line, MAXLINE, fic) != NULL) &&((i >= firstobs) && (i <=lastobs))) {
	linei=linei+1;
for (j=maxwav;j>=1;j--){	for(j=strlen(line); j>=0;j--){ /* Untabifies line */
cutv(stra, strb,line,' '); s[j][i]=atoi(strb);	if(line[j] == '\t')
strcpy(line,stra);	line[j] = ' ';
cutv(stra, strb,line,'/'); anint[j][i]=(double)(atoi(strb)); strcpy(line,stra);	}
cutv(stra, strb,line,' '); mint[j][i]=(double)(atoi(strb)); strcpy(line,stra);	for(j=strlen(line)-1; (line[j]==' ')\|\|(line[j]==10)\|\|(line[j]==13);j--){
	;
	};
	line[j+1]=0; /* Trims blanks at end of line */
	if(line[0]=='#'){
	fprintf(ficlog,"Comment line\n%s\n",line);
	printf("Comment line\n%s\n",line);
	continue;
	}

	for (j=maxwav;j>=1;j--){
	cutv(stra, strb,line,' ');
	errno=0;
	lval=strtol(strb,&endptr,10);
	/* if (errno == ERANGE && (lval == LONG_MAX \|\| lval == LONG_MIN))*/
	if( strb[0]=='\0' \|\| (*endptr != '\0')){
	printf("Error reading data around '%d' at line number %d %s for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
	exit(1);
}	}
	s[j][i]=lval;
cutv(stra, strb,line,'/'); andc[i]=(double)(atoi(strb)); strcpy(line,stra);
cutv(stra, strb,line,' '); moisdc[i]=(double)(atoi(strb)); strcpy(line,stra);	strcpy(line,stra);
	cutv(stra, strb,line,' ');
cutv(stra, strb,line,'/'); annais[i]=(double)(atoi(strb)); strcpy(line,stra);	if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
cutv(stra, strb,line,' '); moisnais[i]=(double)(atoi(strb)); strcpy(line,stra);	}
	else if(iout=sscanf(strb,"%s.") != 0){
cutv(stra, strb,line,' '); weight[i]=(double)(atoi(strb)); strcpy(line,stra);	month=99;
for (j=ncovcol;j>=1;j--){	year=9999;
cutv(stra, strb,line,' '); covar[j][i]=(double)(atoi(strb)); strcpy(line,stra);	}else{
}	printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
lstra=strlen(stra);	exit(1);
if(lstra > 9){ /* More than 2*32 or max of what printf can write with %ld /	}
stratrunc = &(stra[lstra-9]);	anint[j][i]= (double) year;
num[i]=atol(stratrunc);	mint[j][i]= (double)month;
	strcpy(line,stra);
	} /* ENd Waves */

	cutv(stra, strb,line,' ');
	if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
	}
	else if(iout=sscanf(strb,"%s.",dummy) != 0){
	month=99;
	year=9999;
	}else{
	printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
	exit(1);
	}
	andc[i]=(double) year;
	moisdc[i]=(double) month;
	strcpy(line,stra);

	cutv(stra, strb,line,' ');
	if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){
	}
	else if(iout=sscanf(strb,"%s.") != 0){
	month=99;
	year=9999;
	}else{
	printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);
	exit(1);
	}
	annais[i]=(double)(year);
	moisnais[i]=(double)(month);
	strcpy(line,stra);

	cutv(stra, strb,line,' ');
	errno=0;
	lval=strtol(strb,&endptr,10);
	if( strb[0]=='\0' \|\| (*endptr != '\0')){
	printf("Error reading data around '%d' at line number %ld %s for individual %d\nShould be a weight. Exiting.\n",lval, i,line,linei);
	exit(1);
	}
	weight[i]=(double)(lval);
	strcpy(line,stra);

	for (j=ncovcol;j>=1;j--){
	cutv(stra, strb,line,' ');
	errno=0;
	lval=strtol(strb,&endptr,10);
	if( strb[0]=='\0' \|\| (*endptr != '\0')){
	printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line);
	exit(1);
	}
	if(lval <-1 \|\| lval >1){
	printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a value of the %d covar (meaning 0 for the reference or 1. IMaCh does not build design variables, do it your self). Exiting.\n",lval,linei, i,line,j);
	exit(1);
}	}
else	covar[j][i]=(double)(lval);
num[i]=atol(stra);	strcpy(line,stra);
	}
/*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){	lstra=strlen(stra);
printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/
	if(lstra > 9){ /* More than 2*32 or max of what printf can write with %ld /
i=i+1;	stratrunc = &(stra[lstra-9]);
}	num[i]=atol(stratrunc);
}	}
	else
	num[i]=atol(stra);
	/*if((s[2][i]==2) && (s[3][i]==-1)&&(s[4][i]==9)){
	printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]),weight[i], (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i])); ij=ij+1;}*/

	i=i+1;
	} /* End loop reading data */
	fclose(fic);
/* printf("ii=%d", ij);	/* printf("ii=%d", ij);
scanf("%d",i);*/	scanf("%d",i);*/
imx=i-1; /* Number of individuals */	imx=i-1; /* Number of individuals */
Line 4428 int main(int argc, char *argv[])	Line 4668 int main(int argc, char *argv[])
if (s[4][i]==9) s[4][i]=-1;	if (s[4][i]==9) s[4][i]=-1;
printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}*/	printf("%ld %.lf %.lf %.lf %.lf/%.lf %.lf/%.lf %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d %.lf/%.lf %d\n",num[i],(covar[1][i]), (covar[2][i]), (weight[i]), (moisnais[i]), (annais[i]), (moisdc[i]), (andc[i]), (mint[1][i]), (anint[1][i]), (s[1][i]), (mint[2][i]), (anint[2][i]), (s[2][i]), (mint[3][i]), (anint[3][i]), (s[3][i]), (mint[4][i]), (anint[4][i]), (s[4][i]));}*/

for (i=1; i<=imx; i++)	/* for (i=1; i<=imx; i++) */

/*if ((s[3][i]==3) \|\| (s[4][i]==3)) weight[i]=0.08;	/*if ((s[3][i]==3) \|\| (s[4][i]==3)) weight[i]=0.08;
else weight[i]=1;*/	else weight[i]=1;*/

/* Calculation of the number of parameter from char model*/	/* Calculation of the number of parameters from char model */
Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */	Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
Tprod=ivector(1,15);	Tprod=ivector(1,15);
Tvaraff=ivector(1,15);	Tvaraff=ivector(1,15);
Line 4512 int main(int argc, char *argv[])	Line 4752 int main(int argc, char *argv[])
printf("cptcovprod=%d ", cptcovprod);	printf("cptcovprod=%d ", cptcovprod);
fprintf(ficlog,"cptcovprod=%d ", cptcovprod);	fprintf(ficlog,"cptcovprod=%d ", cptcovprod);

scanf("%d ",i);	scanf("%d ",i);*/
fclose(fic);*/

/* if(mle==1){*/	/* if(mle==1){*/
if (weightopt != 1) { /* Maximisation without weights*/	if (weightopt != 1) { /* Maximisation without weights*/
Line 4546 int main(int argc, char *argv[])	Line 4785 int main(int argc, char *argv[])
for (i=1; i<=imx; i++) {	for (i=1; i<=imx; i++) {
agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);	agedc[i]=(moisdc[i]/12.+andc[i])-(moisnais[i]/12.+annais[i]);
for(m=firstpass; (m<= lastpass); m++){	for(m=firstpass; (m<= lastpass); m++){
if(s[m][i] >0){	if(s[m][i] >0 \|\| s[m][i]==-2 \|\| s[m][i]==-4 \|\| s[m][i]==-5){
if (s[m][i] >= nlstate+1) {	if (s[m][i] >= nlstate+1) {
if(agedc[i]>0)	if(agedc[i]>0)
if((int)moisdc[i]!=99 && (int)andc[i]!=9999)	if((int)moisdc[i]!=99 && (int)andc[i]!=9999)
Line 4562 int main(int argc, char *argv[])	Line 4801 int main(int argc, char *argv[])
}	}
}	}
else if(s[m][i] !=9){ /* Standard case, age in fractional	else if(s[m][i] !=9){ /* Standard case, age in fractional
years but with the precision of a	years but with the precision of a month */
month */
agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);	agev[m][i]=(mint[m][i]/12.+1./24.+anint[m][i])-(moisnais[i]/12.+1./24.+annais[i]);
if((int)mint[m][i]==99 \|\| (int)anint[m][i]==9999)	if((int)mint[m][i]==99 \|\| (int)anint[m][i]==9999)
agev[m][i]=1;	agev[m][i]=1;
Line 4720 Title=%s <br>Datafile=%s Firstpass=%d La	Line 4958 Title=%s <br>Datafile=%s Firstpass=%d La

/* Calculates basic frequencies. Computes observed prevalence at single age	/* Calculates basic frequencies. Computes observed prevalence at single age
and prints on file fileres'p'. */	and prints on file fileres'p'. */
freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint);	freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint,strstart);

fprintf(fichtm,"\n");	fprintf(fichtm,"\n");
fprintf(fichtm,"<br>Total number of observations=%d <br>\n\	fprintf(fichtm,"<br>Total number of observations=%d <br>\n\
Line 4739 Interval (in months) between two waves:	Line 4977 Interval (in months) between two waves:
p=param[1][1]; /* (((param +1)+1)+0) /	p=param[1][1]; /* (((param +1)+1)+0) /

globpr=0; /* To get the number ipmx of contributions and the sum of weights*/	globpr=0; /* To get the number ipmx of contributions and the sum of weights*/

if (mle==-3){	if (mle==-3){
ximort=matrix(1,NDIM,1,NDIM);	ximort=matrix(1,NDIM,1,NDIM);
cens=ivector(1,n);	cens=ivector(1,n);
Line 4748 Interval (in months) between two waves:	Line 4987 Interval (in months) between two waves:

for (i=1; i<=imx; i++){	for (i=1; i<=imx; i++){
dcwave[i]=-1;	dcwave[i]=-1;
for (j=1; j<=lastpass; j++)	for (m=firstpass; m<=lastpass; m++)
if (s[j][i]>nlstate) {	if (s[m][i]>nlstate) {
dcwave[i]=j;	dcwave[i]=m;
/* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/	/* printf("i=%d j=%d s=%d dcwave=%d\n",i,j, s[j][i],dcwave[i]);*/
break;	break;
}	}
Line 4759 Interval (in months) between two waves:	Line 4998 Interval (in months) between two waves:
for (i=1; i<=imx; i++) {	for (i=1; i<=imx; i++) {
if (wav[i]>0){	if (wav[i]>0){
ageexmed[i]=agev[mw[1][i]][i];	ageexmed[i]=agev[mw[1][i]][i];
j=wav[i];agecens[i]=1.;	j=wav[i];
if (ageexmed[i]>1 & wav[i]>0) agecens[i]=agev[mw[j][i]][i];	agecens[i]=1.;
cens[i]=1;
	if (ageexmed[i]> 1 && wav[i] > 0){
if (ageexmed[i]<1) cens[i]=-1;	agecens[i]=agev[mw[j][i]][i];
if (agedc[i]< AGESUP & agedc[i]>1 & dcwave[i]>firstpass & dcwave[i]<=lastpass) cens[i]=0 ;	cens[i]= 1;
	}else if (ageexmed[i]< 1)
	cens[i]= -1;
	if (agedc[i]< AGESUP && agedc[i]>1 && dcwave[i]>firstpass && dcwave[i]<=lastpass)
	cens[i]=0 ;
}	}
else cens[i]=-1;	else cens[i]=-1;
}	}
Line 4773 Interval (in months) between two waves:	Line 5016 Interval (in months) between two waves:
for (j=1;j<=NDIM;j++)	for (j=1;j<=NDIM;j++)
ximort[i][j]=(i == j ? 1.0 : 0.0);	ximort[i][j]=(i == j ? 1.0 : 0.0);
}	}

p[1]=0.1; p[2]=0.1;	p[1]=0.0268; p[NDIM]=0.083;
/printf("%lf %lf", p[1], p[2]);/	/printf("%lf %lf", p[1], p[2]);/


printf("Powell\n"); fprintf(ficlog,"Powell\n");	printf("Powell\n"); fprintf(ficlog,"Powell\n");
strcpy(filerespow,"pow-mort");	strcpy(filerespow,"pow-mort");
strcat(filerespow,fileres);	strcat(filerespow,fileres);
if((ficrespow=fopen(filerespow,"w"))==NULL) {	if((ficrespow=fopen(filerespow,"w"))==NULL) {
printf("Problem with resultfile: %s\n", filerespow);	printf("Problem with resultfile: %s\n", filerespow);
fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);	fprintf(ficlog,"Problem with resultfile: %s\n", filerespow);
}	}
fprintf(ficrespow,"# Powell\n# iter -2*LL");	fprintf(ficrespow,"# Powell\n# iter -2*LL");
/* for (i=1;i<=nlstate;i++)	/* for (i=1;i<=nlstate;i++)
for(j=1;j<=nlstate+ndeath;j++)	for(j=1;j<=nlstate+ndeath;j++)
if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);	if(j!=i)fprintf(ficrespow," p%1d%1d",i,j);
*/	*/
fprintf(ficrespow,"\n");	fprintf(ficrespow,"\n");

powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);	powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz);
fclose(ficrespow);	fclose(ficrespow);

hesscov(matcov, p, NDIM,delti, 1e-4, gompertz);	hesscov(matcov, p, NDIM, delti, 1e-4, gompertz);

for(i=1; i <=NDIM; i++)	for(i=1; i <=NDIM; i++)
for(j=i+1;j<=NDIM;j++)	for(j=i+1;j<=NDIM;j++)
Line 4812 Interval (in months) between two waves:	Line 5055 Interval (in months) between two waves:
printf("iter=%d MLE=%f Eq=%lfexp(%lf(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);	printf("iter=%d MLE=%f Eq=%lfexp(%lf(age-%d))\n",iter,-gompertz(p),p[1],p[2],agegomp);
for (i=1;i<=NDIM;i++)	for (i=1;i<=NDIM;i++)
printf("%f [%f ; %f]\n",p[i],p[i]-2sqrt(matcov[i][i]),p[i]+2sqrt(matcov[i][i]));	printf("%f [%f ; %f]\n",p[i],p[i]-2sqrt(matcov[i][i]),p[i]+2sqrt(matcov[i][i]));

	lsurv=vector(1,AGESUP);
	lpop=vector(1,AGESUP);
	tpop=vector(1,AGESUP);
	lsurv[agegomp]=100000;

	for (k=agegomp;k<=AGESUP;k++) {
	agemortsup=k;
	if (p[1]exp(p[2](k-agegomp))>1) break;
	}

	for (k=agegomp;k<agemortsup;k++)
	lsurv[k+1]=lsurv[k]-lsurv[k](p[1]exp(p[2]*(k-agegomp)));

	for (k=agegomp;k<agemortsup;k++){
	lpop[k]=(lsurv[k]+lsurv[k+1])/2.;
	sumlpop=sumlpop+lpop[k];
	}

	tpop[agegomp]=sumlpop;
	for (k=agegomp;k<(agemortsup-3);k++){
	/* tpop[k+1]=2;*/
	tpop[k+1]=tpop[k]-lpop[k];
	}


	printf("\nAge lx qx dx Lx Tx e(x)\n");
	for (k=agegomp;k<(agemortsup-2);k++)
	printf("%d %.0lf %lf %.0lf %.0lf %.0lf %lf\n",k,lsurv[k],p[1]exp(p[2](k-agegomp)),(p[1]exp(p[2](k-agegomp)))*lsurv[k],lpop[k],tpop[k],tpop[k]/lsurv[k]);


replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */	replace_back_to_slash(pathc,path); /* Even gnuplot wants a / */
printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);	printinggnuplotmort(fileres, optionfilefiname,ageminpar,agemaxpar,fage, pathc,p);

printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \	printinghtmlmort(fileres,title,datafile, firstpass, lastpass, \
stepm, weightopt,\	stepm, weightopt,\
model,imx,p,matcov);	model,imx,p,matcov,agemortsup);

	free_vector(lsurv,1,AGESUP);
	free_vector(lpop,1,AGESUP);
	free_vector(tpop,1,AGESUP);
} /* Endof if mle==-3 */	} /* Endof if mle==-3 */

else{ /* For mle >=1 */	else{ /* For mle >=1 */

likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */	likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
Line 5077 Interval (in months) between two waves:	Line 5355 Interval (in months) between two waves:
}	}
printf("Computing stable prevalence: result on file '%s' \n", filerespl);	printf("Computing stable prevalence: result on file '%s' \n", filerespl);
fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);	fprintf(ficlog,"Computing stable prevalence: result on file '%s' \n", filerespl);
	fprintf(ficrespl, "#Local time at start: %s", strstart);
fprintf(ficrespl,"#Stable prevalence \n");	fprintf(ficrespl,"#Stable prevalence \n");
fprintf(ficrespl,"#Age ");	fprintf(ficrespl,"#Age ");
for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);	for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
Line 5137 Interval (in months) between two waves:	Line 5416 Interval (in months) between two waves:
hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */	hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */

/* hstepm=1; aff par mois*/	/* hstepm=1; aff par mois*/
	fprintf(ficrespij, "#Local time at start: %s", strstart);
fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");	fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
for(cptcov=1,k=0;cptcov<=i1;cptcov++){	for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){	for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
Line 5174 Interval (in months) between two waves:	Line 5453 Interval (in months) between two waves:
}	}
}	}

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

fclose(ficrespij);	fclose(ficrespij);

Line 5263 Interval (in months) between two waves:	Line 5542 Interval (in months) between two waves:

eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);	eij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;	oldm=oldms;savm=savms;
evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov);	evsij(fileres, eij, p, nlstate, stepm, (int) bage, (int)fage, oldm, savm, k, estepm, delti, matcov, strstart);

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

	fprintf(ficrest, "#Local time at start: %s", strstart);
fprintf(ficrest,"#Total LEs with variances: e.. (std) ");	fprintf(ficrest,"#Total LEs with variances: e.. (std) ");
for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);	for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
fprintf(ficrest,"\n");	fprintf(ficrest,"\n");
Line 5345 Interval (in months) between two waves:	Line 5624 Interval (in months) between two waves:

varpl=matrix(1,nlstate,(int) bage, (int) fage);	varpl=matrix(1,nlstate,(int) bage, (int) fage);
oldm=oldms;savm=savms;	oldm=oldms;savm=savms;
varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k);	varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k,strstart);
free_matrix(varpl,1,nlstate,(int) bage, (int)fage);	free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
}	}
}	}
Line 5357 Interval (in months) between two waves:	Line 5636 Interval (in months) between two waves:
free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);	free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);

} /* mle==-3 arrives here for freeing */	} /* mle==-3 arrives here for freeing */
	free_matrix(prlim,1,nlstate,1,nlstate);
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(newms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(savms, 1,nlstate+ndeath,1,nlstate+ndeath);

free_matrix(covar,0,NCOVMAX,1,n);	free_matrix(covar,0,NCOVMAX,1,n);
free_matrix(matcov,1,npar,1,npar);	free_matrix(matcov,1,npar,1,npar);
/free_vector(delti,1,npar);/	/free_vector(delti,1,npar);/
Line 5376 Interval (in months) between two waves:	Line 5655 Interval (in months) between two waves:
free_ivector(Tage,1,15);	free_ivector(Tage,1,15);
free_ivector(Tcode,1,100);	free_ivector(Tcode,1,100);

	free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
	free_imatrix(codtab,1,100,1,10);
fflush(fichtm);	fflush(fichtm);
fflush(ficgp);	fflush(ficgp);

Line 5394 Interval (in months) between two waves:	Line 5674 Interval (in months) between two waves:
tm = *localtime(&end_time.tv_sec);	tm = *localtime(&end_time.tv_sec);
tmg = *gmtime(&end_time.tv_sec);	tmg = *gmtime(&end_time.tv_sec);
strcpy(strtend,asctime(&tm));	strcpy(strtend,asctime(&tm));
printf("Local time at start %s\nLocaltime at end %s",strstart, strtend);	printf("Local time at start %s\nLocal time at end %s",strstart, strtend);
fprintf(ficlog,"Local time at start %s\nLocal time at end %s\n",strstart, strtend);	fprintf(ficlog,"Local time at start %s\nLocal time at end %s\n",strstart, strtend);
printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));	printf("Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));

Line 5411 Interval (in months) between two waves:	Line 5691 Interval (in months) between two waves:
/------ End -----------/	/------ End -----------/

chdir(path);	chdir(path);
strcpy(plotcmd,GNUPLOTPROGRAM);	/strcat(plotcmd,CHARSEPARATOR);/
strcat(plotcmd," ");	sprintf(plotcmd,"gnuplot");
strcat(plotcmd,optionfilegnuplot);	#ifndef UNIX
printf("Starting graphs with: %s",plotcmd);fflush(stdout);	sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);
	#endif
	if(!stat(plotcmd,&info)){
	printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
	if(!stat(getenv("GNUPLOTBIN"),&info)){
	printf("Error gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout);
	}else
	strcpy(pplotcmd,plotcmd);
	#ifdef UNIX
	strcpy(plotcmd,GNUPLOTPROGRAM);
	if(!stat(plotcmd,&info)){
	printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);
	}else
	strcpy(pplotcmd,plotcmd);
	#endif
	}else
	strcpy(pplotcmd,plotcmd);

	sprintf(plotcmd,"%s %s",pplotcmd, optionfilegnuplot);
	printf("Starting graphs with: %s\n",plotcmd);fflush(stdout);

if((outcmd=system(plotcmd)) != 0){	if((outcmd=system(plotcmd)) != 0){
printf(" Problem with gnuplot\n");	printf("\n Problem with gnuplot\n");
}	}
printf(" Wait...");	printf(" Wait...");
while (z[0] != 'q') {	while (z[0] != 'q') {
Line 5424 Interval (in months) between two waves:	Line 5724 Interval (in months) between two waves:
printf("\nType e to edit output files, g to graph again and q for exiting: ");	printf("\nType e to edit output files, g to graph again and q for exiting: ");
scanf("%s",z);	scanf("%s",z);
/* if (z[0] == 'c') system("./imach"); */	/* if (z[0] == 'c') system("./imach"); */
if (z[0] == 'e') system(optionfilehtm);	if (z[0] == 'e') {
	printf("Starting browser with: %s",optionfilehtm);fflush(stdout);
	system(optionfilehtm);
	}
else if (z[0] == 'g') system(plotcmd);	else if (z[0] == 'g') system(plotcmd);
else if (z[0] == 'q') exit(0);	else if (z[0] == 'q') exit(0);
}	}

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

Removed from v.1.98
changed lines
	Added in v.1.115