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.155 and 1.165

version 1.155, 2014/08/25 18:32:34	version 1.165, 2014/12/16 11:20:36
Line 1	Line 1
/* $Id$	/* $Id$
$State$	$State$
$Log$	$Log$
	Revision 1.165 2014/12/16 11:20:36 brouard
	Summary: After compiling on Visual C

	* imach.c (Module): Merging 1.61 to 1.162

	Revision 1.164 2014/12/16 10:52:11 brouard
	Summary: Merging with Visual C after suppressing some warnings for unused variables. Also fixing Saito's bug 0.98Xn

	* imach.c (Module): Merging 1.61 to 1.162

	Revision 1.163 2014/12/16 10:30:11 brouard
	* imach.c (Module): Merging 1.61 to 1.162

	Revision 1.162 2014/09/25 11:43:39 brouard
	Summary: temporary backup 0.99!

	Revision 1.1 2014/09/16 11:06:58 brouard
	Summary: With some code (wrong) for nlopt

	Author:

	Revision 1.161 2014/09/15 20:41:41 brouard
	Summary: Problem with macro SQR on Intel compiler

	Revision 1.160 2014/09/02 09:24:05 brouard
	* empty log message *

	Revision 1.159 2014/09/01 10:34:10 brouard
	Summary: WIN32
	Author: Brouard

	Revision 1.158 2014/08/27 17:11:51 brouard
	* empty log message *

	Revision 1.157 2014/08/27 16:26:55 brouard
	Summary: Preparing windows Visual studio version
	Author: Brouard

	In order to compile on Visual studio, time.h is now correct and time_t
	and tm struct should be used. difftime should be used but sometimes I
	just make the differences in raw time format (time(&now).
	Trying to suppress #ifdef LINUX
	Add xdg-open for __linux in order to open default browser.

	Revision 1.156 2014/08/25 20:10:10 brouard
	* empty log message *

Revision 1.155 2014/08/25 18:32:34 brouard	Revision 1.155 2014/08/25 18:32:34 brouard
Summary: New compile, minor changes	Summary: New compile, minor changes
Author: Brouard	Author: Brouard
Line 455	Line 502
end	end
*/	*/

	#define POWELL /* Instead of NLOPT */



#include <math.h>	#include <math.h>
#include <stdio.h>	#include <stdio.h>
#include <stdlib.h>	#include <stdlib.h>
#include <string.h>	#include <string.h>

	#ifdef _WIN32
	#include <io.h>
	#else
#include <unistd.h>	#include <unistd.h>
	#endif

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

	/* #ifdef LINUX */
	/* #include <time.h> */
	/* #include "timeval.h" */
	/* #else */
	/* #include <sys/time.h> */
	/* #endif */

#ifdef LINUX
#include <time.h>	#include <time.h>
#include "timeval.h"
#else
#include <sys/time.h>
#endif

#ifdef GSL	#ifdef GSL
#include <gsl/gsl_errno.h>	#include <gsl/gsl_errno.h>
#include <gsl/gsl_multimin.h>	#include <gsl/gsl_multimin.h>
#endif	#endif

	#ifdef NLOPT
	#include <nlopt.h>
	typedef struct {
	double (* function)(double [] );
	} myfunc_data ;
	#endif

/* #include <libintl.h> */	/* #include <libintl.h> */
/* #define _(String) gettext (String) */	/* #define _(String) gettext (String) */

Line 506 extern int errno;	Line 566 extern int errno;
#define YEARM 12. /*< Number of months per year /	#define YEARM 12. /*< Number of months per year /
#define AGESUP 130	#define AGESUP 130
#define AGEBASE 40	#define AGEBASE 40
#define AGEGOMP 10. /*< Minimal age for Gompertz adjustment /	#define AGEGOMP 10 /*< Minimal age for Gompertz adjustment /
#ifdef UNIX	#ifdef _WIN32
#define DIRSEPARATOR '/'
#define CHARSEPARATOR "/"
#define ODIRSEPARATOR '\\'
#else
#define DIRSEPARATOR '\\'	#define DIRSEPARATOR '\\'
#define CHARSEPARATOR "\\"	#define CHARSEPARATOR "\\"
#define ODIRSEPARATOR '/'	#define ODIRSEPARATOR '/'
	#else
	#define DIRSEPARATOR '/'
	#define CHARSEPARATOR "/"
	#define ODIRSEPARATOR '\\'
#endif	#endif

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

char version[]="Imach version 0.98nV, August 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";	char version[]="Imach version 0.99, September 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";
char fullversion[]="$Revision$ $Date$";	char fullversion[]="$Revision$ $Date$";
char strstart[80];	char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];	char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 551 int *mw; / mw[mi][i] is number of the	Line 611 int *mw; / mw[mi][i] is number of the
int *dh; / dh[mi][i] is number of steps between mi,mi+1 for this individual */	int *dh; / dh[mi][i] is number of steps between mi,mi+1 for this individual */
int *bh; / bh[mi][i] is the bias (+ or -) for this individual if the delay between	int *bh; / bh[mi][i] is the bias (+ or -) for this individual if the delay between
* wave mi and wave mi+1 is not an exact multiple of stepm. */	* wave mi and wave mi+1 is not an exact multiple of stepm. */
	int countcallfunc=0; /* Count the number of calls to func */
double jmean=1; /* Mean space between 2 waves */	double jmean=1; /* Mean space between 2 waves */
double *matprod2(); / test */	double *matprod2(); / test */
double oldm, newm, *savm; / Working pointers to matrices */	double oldm, newm, *savm; / Working pointers to matrices */
Line 594 char popfile[FILENAMELENGTH];	Line 655 char popfile[FILENAMELENGTH];

char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilehtmcov[FILENAMELENGTH] ;	char optionfilegnuplot[FILENAMELENGTH], optionfilehtm[FILENAMELENGTH], optionfilehtmcov[FILENAMELENGTH] ;

struct timeval start_time, end_time, curr_time, last_time, forecast_time;	/* struct timeval start_time, end_time, curr_time, last_time, forecast_time; */
struct timezone tzp;	/* struct timezone tzp; */
extern int gettimeofday();	/* extern int gettimeofday(); */
struct tm tmg, tm, tmf, gmtime(), localtime();	struct tm tml, gmtime(), localtime();
long time_value;
extern long time();	extern time_t time();

	struct tm start_time, end_time, curr_time, last_time, forecast_time;
	time_t rstart_time, rend_time, rcurr_time, rlast_time, rforecast_time; /* raw time */
	struct tm tm;

char strcurr[80], strfor[80];	char strcurr[80], strfor[80];

char *endptr;	char *endptr;
Line 754 char cutl(char blocc, char *alocc, cha	Line 820 char cutl(char blocc, char *alocc, cha
gives blocc="abcdef2ghi" and alocc="j".	gives blocc="abcdef2ghi" and alocc="j".
If occ is not found blocc is null and alocc is equal to in. Returns blocc	If occ is not found blocc is null and alocc is equal to in. Returns blocc
*/	*/
char s, t, *bl;	char s, t;
t=in;s=in;	t=in;s=in;
while ((in != occ) && (in != '\0')){	while ((in != occ) && (in != '\0')){
alocc++ = in++;	alocc++ = in++;
Line 838 int nbocc(char *s, char occ)	Line 904 int nbocc(char *s, char occ)
/* } */	/* } */
/* } */	/* } */

	#ifdef _WIN32
	char * strsep(char *pp, const char delim)
	{
	char p, q;

	if ((p = *pp) == NULL)
	return 0;
	if ((q = strpbrk (p, delim)) != NULL)
	{
	*pp = q + 1;
	*q = '\0';
	}
	else
	*pp = 0;
	return p;
	}
	#endif

/******************** nrerror ******************/	/******************** nrerror ******************/

void nrerror(char error_text[])	void nrerror(char error_text[])
Line 1037 char subdirf3(char fileres[], char pre	Line 1121 char subdirf3(char fileres[], char pre
return tmpout;	return tmpout;
}	}

	char *asc_diff_time(long time_sec, char ascdiff[])
	{
	long sec_left, days, hours, minutes;
	days = (time_sec) / (606024);
	sec_left = (time_sec) % (606024);
	hours = (sec_left) / (60*60) ;
	sec_left = (sec_left) %(60*60);
	minutes = (sec_left) /60;
	sec_left = (sec_left) % (60);
	sprintf(ascdiff,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left);
	return ascdiff;
	}

/*************** f1dim ***********************/	/*************** f1dim ***********************/
extern int ncom;	extern int ncom;
extern double pcom,xicom;	extern double pcom,xicom;
Line 1060 double brent(double ax, double bx, doubl	Line 1157 double brent(double ax, double bx, doubl
{	{
int iter;	int iter;
double a,b,d,etemp;	double a,b,d,etemp;
double fu,fv,fw,fx;	double fu=0,fv,fw,fx;
double ftemp;	double ftemp=0.;
double p,q,r,tol1,tol2,u,v,w,x,xm;	double p,q,r,tol1,tol2,u,v,w,x,xm;
double e=0.0;	double e=0.0;

Line 1075 double brent(double ax, double bx, doubl	Line 1172 double brent(double ax, double bx, doubl
/* if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret)))*/	/* if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret)))*/
printf(".");fflush(stdout);	printf(".");fflush(stdout);
fprintf(ficlog,".");fflush(ficlog);	fprintf(ficlog,".");fflush(ficlog);
#ifdef DEBUG	#ifdef DEBUGBRENT
printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);	printf("br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);	fprintf(ficlog,"br %d,x=%.10e xm=%.10e b=%.10e a=%.10e tol=%.10e tol1=%.10e tol2=%.10e x-xm=%.10e fx=%.12e fu=%.12e,fw=%.12e,ftemp=%.12e,ftol=%.12e\n",iter,x,xm,b,a,tol,tol1,tol2,(x-xm),fx,fu,fw,ftemp,ftol);
/* if ((fabs(x-xm) <= (tol2-0.5(b-a)))\|\|(2.0fabs(fu-ftemp) <= ftol1.e-2(fabs(fu)+fabs(ftemp)))) { */	/* if ((fabs(x-xm) <= (tol2-0.5(b-a)))\|\|(2.0fabs(fu-ftemp) <= ftol1.e-2(fabs(fu)+fabs(ftemp)))) { */
Line 1145 void mnbrak(double ax, double bx, doub	Line 1242 void mnbrak(double ax, double bx, doub
}	}
cx=(bx)+GOLD(bx-*ax);	cx=(bx)+GOLD(bx-*ax);
fc=(func)(*cx);	fc=(func)(*cx);
while (fb > fc) {	while (fb > fc) { /* Declining fa, fb, fc */
r=(bx-ax)(fb-*fc);	r=(bx-ax)(fb-*fc);
q=(bx-cx)(fb-*fa);	q=(bx-cx)(fb-*fa);
u=(bx)-((bx-cx)q-(bx-ax)*r)/	u=(bx)-((bx-cx)q-(bx-ax)*r)/
(2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));	(2.0SIGN(FMAX(fabs(q-r),TINY),q-r)); / Minimum abscisse of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */
ulim=(bx)+GLIMIT(cx-bx);	ulim=(bx)+GLIMIT(cx-bx); /* Maximum abscisse where function can be evaluated */
if ((bx-u)(u-*cx) > 0.0) {	if ((bx-u)(u-cx) > 0.0) { / if u between b and c */
fu=(*func)(u);	fu=(*func)(u);
} else if ((cx-u)(u-ulim) > 0.0) {	#ifdef DEBUG
	/* f(x)=A(x-u)*2+f(u) /
	double A, fparabu;
	A= (fb - fa)/(bx-ax)/(bx+ax-2*u);
	fparabu= fa - A(ax-u)(*ax-u);
	printf("mnbrak (ax=%.12f, fa=%.12lf), (bx=%.12f, fb=%.12lf), (cx=%.12f, fc=%.12lf), (u=%.12f, fu=%.12lf, fparabu=%.12f)\n",ax,fa,bx,fb,cx,*fc,u,fu, fparabu);
	fprintf(ficlog, "mnbrak (ax=%.12f, fa=%.12lf), (bx=%.12f, fb=%.12lf), (cx=%.12f, fc=%.12lf), (u=%.12f, fu=%.12lf, fparabu=%.12f)\n",ax,fa,bx,fb,cx,*fc,u,fu, fparabu);
	#endif
	} else if ((cx-u)(u-ulim) > 0.0) { /* u is after c but before ulim */
fu=(*func)(u);	fu=(*func)(u);
if (fu < *fc) {	if (fu < *fc) {
SHFT(bx,cx,u,cx+GOLD(cx-bx))	SHFT(bx,cx,u,cx+GOLD(cx-bx))
SHFT(fb,fc,fu,(*func)(u))	SHFT(fb,fc,fu,(*func)(u))
}	}
} else if ((u-ulim)(ulim-cx) >= 0.0) {	} else if ((u-ulim)(ulim-cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */
u=ulim;	u=ulim;
fu=(*func)(u);	fu=(*func)(u);
} else {	} else {
Line 1172 void mnbrak(double ax, double bx, doub	Line 1277 void mnbrak(double ax, double bx, doub
}	}

/************* linmin **********************/	/************* linmin **********************/
	/* Given an n -dimensional point p[1..n] and an n -dimensional direction xi[1..n] , moves and
	resets p to where the function func(p) takes on a minimum along the direction xi from p ,
	and replaces xi by the actual vector displacement that p was moved. Also returns as fret
	the value of func at the returned location p . This is actually all accomplished by calling the
	routines mnbrak and brent .*/
int ncom;	int ncom;
double pcom,xicom;	double pcom,xicom;
double (*nrfunc)(double []);	double (*nrfunc)(double []);
Line 1198 void linmin(double p[], double xi[], int	Line 1307 void linmin(double p[], double xi[], int
}	}
ax=0.0;	ax=0.0;
xx=1.0;	xx=1.0;
mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);	mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); /* Find a bracket a,x,b in direction n=xi ie xicom */
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);	fret=brent(ax,xx,bx,f1dim,TOL,&xmin); / Find a minimum P+lambda n in that direction (lambdamin), with TOL between abscisses */
#ifdef DEBUG	#ifdef DEBUG
printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);	printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);	fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);
Line 1212 void linmin(double p[], double xi[], int	Line 1321 void linmin(double p[], double xi[], int
free_vector(pcom,1,n);	free_vector(pcom,1,n);
}	}

char *asc_diff_time(long time_sec, char ascdiff[])
{
long sec_left, days, hours, minutes;
days = (time_sec) / (606024);
sec_left = (time_sec) % (606024);
hours = (sec_left) / (60*60) ;
sec_left = (sec_left) %(60*60);
minutes = (sec_left) /60;
sec_left = (sec_left) % (60);
sprintf(ascdiff,"%ld day(s) %ld hour(s) %ld minute(s) %ld second(s)",days, hours, minutes, sec_left);
return ascdiff;
}

/************* powell **********************/	/************* powell **********************/
	/*
	Minimization of a function func of n variables. Input consists of an initial starting point
	p[1..n] ; an initial matrix xi[1..n][1..n] , whose columns contain the initial set of di-
	rections (usually the n unit vectors); and ftol , the fractional tolerance in the function value
	such that failure to decrease by more than this amount on one iteration signals doneness. On
	output, p is set to the best point found, xi is the then-current direction set, fret is the returned
	function value at p , and iter is the number of iterations taken. The routine linmin is used.
	*/
void powell(double p[], double *xi, int n, double ftol, int iter, double *fret,	void powell(double p[], double *xi, int n, double ftol, int iter, double *fret,
double (*func)(double []))	double (*func)(double []))
{	{
Line 1243 void powell(double p[], double **xi, int	Line 1348 void powell(double p[], double **xi, int
xits=vector(1,n);	xits=vector(1,n);
fret=(func)(p);	fret=(func)(p);
for (j=1;j<=n;j++) pt[j]=p[j];	for (j=1;j<=n;j++) pt[j]=p[j];
	rcurr_time = time(NULL);
for (iter=1;;++(iter)) {	for (iter=1;;++(iter)) {
fp=(*fret);	fp=(*fret);
ibig=0;	ibig=0;
del=0.0;	del=0.0;
last_time=curr_time;	rlast_time=rcurr_time;
(void) gettimeofday(&curr_time,&tzp);	/* (void) gettimeofday(&curr_time,&tzp); */
printf("\nPowell iter=%d -2LL=%.12f %ld sec. %ld sec.",iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec);fflush(stdout);	rcurr_time = time(NULL);
fprintf(ficlog,"\nPowell iter=%d -2LL=%.12f %ld sec. %ld sec.",iter,*fret, curr_time.tv_sec-last_time.tv_sec, curr_time.tv_sec-start_time.tv_sec); fflush(ficlog);	curr_time = *localtime(&rcurr_time);
/* fprintf(ficrespow,"%d %.12f %ld",iter,fret,curr_time.tv_sec-start_time.tv_sec); */	printf("\nPowell iter=%d -2LL=%.12f %ld sec. %ld sec.",iter,*fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
	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); */
for (i=1;i<=n;i++) {	for (i=1;i<=n;i++) {
printf(" %d %.12f",i, p[i]);	printf(" %d %.12f",i, p[i]);
fprintf(ficlog," %d %.12lf",i, p[i]);	fprintf(ficlog," %d %.12lf",i, p[i]);
Line 1261 void powell(double p[], double **xi, int	Line 1369 void powell(double p[], double **xi, int
fprintf(ficlog,"\n");	fprintf(ficlog,"\n");
fprintf(ficrespow,"\n");fflush(ficrespow);	fprintf(ficrespow,"\n");fflush(ficrespow);
if(*iter <=3){	if(*iter <=3){
tm = *localtime(&curr_time.tv_sec);	tml = *localtime(&rcurr_time);
strcpy(strcurr,asctime(&tm));	strcpy(strcurr,asctime(&tml));
/* asctime_r(&tm,strcurr); */	rforecast_time=rcurr_time;
forecast_time=curr_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 this last iteration #%d: %ld seconds,\n",*iter,curr_time.tv_sec-last_time.tv_sec);	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,curr_time.tv_sec-last_time.tv_sec);	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){
forecast_time.tv_sec=curr_time.tv_sec+(niterf-iter)(curr_time.tv_sec-last_time.tv_sec);	rforecast_time=rcurr_time+(niterf-iter)(rcurr_time-rlast_time);
tmf = *localtime(&forecast_time.tv_sec);	forecast_time = *localtime(&rforecast_time);
/* asctime_r(&tmf,strfor); */	strcpy(strfor,asctime(&forecast_time));
strcpy(strfor,asctime(&tmf));
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(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(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(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(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
}	}
}	}
for (i=1;i<=n;i++) {	for (i=1;i<=n;i++) {
for (j=1;j<=n;j++) xit[j]=xi[j][i];	for (j=1;j<=n;j++) xit[j]=xi[j][i];
fptt=(*fret);	fptt=(*fret);
#ifdef DEBUG	#ifdef DEBUG
printf("fret=%lf \n",*fret);	printf("fret=%lf, %lf, %lf \n", fret, fret, *fret);
fprintf(ficlog,"fret=%lf \n",*fret);	fprintf(ficlog, "fret=%lf, %lf, %lf \n", fret, fret, *fret);
#endif	#endif
printf("%d",i);fflush(stdout);	printf("%d",i);fflush(stdout);
fprintf(ficlog,"%d",i);fflush(ficlog);	fprintf(ficlog,"%d",i);fflush(ficlog);
Line 1305 void powell(double p[], double **xi, int	Line 1411 void powell(double p[], double **xi, int
fprintf(ficlog," x(%d)=%.12e",j,xit[j]);	fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
}	}
for(j=1;j<=n;j++) {	for(j=1;j<=n;j++) {
printf(" p=%.12e",p[j]);	printf(" p(%d)=%.12e",j,p[j]);
fprintf(ficlog," p=%.12e",p[j]);	fprintf(ficlog," p(%d)=%.12e",j,p[j]);
}	}
printf("\n");	printf("\n");
fprintf(ficlog,"\n");	fprintf(ficlog,"\n");
#endif	#endif
}	} /* end i */
if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret))) {	if (2.0fabs(fp-(fret)) <= ftol(fabs(fp)+fabs(fret))) {
#ifdef DEBUG	#ifdef DEBUG
int k[2],l;	int k[2],l;
Line 1344 void powell(double p[], double **xi, int	Line 1450 void powell(double p[], double **xi, int
return;	return;
}	}
if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");	if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");
for (j=1;j<=n;j++) {	for (j=1;j<=n;j++) { /* Computes an extrapolated point */
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];
pt[j]=p[j];	pt[j]=p[j];
}	}
fptt=(*func)(ptt);	fptt=(*func)(ptt);
if (fptt < fp) {	if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
t=2.0(fp-2.0(fret)+fptt)SQR(fp-(fret)-del)-delSQR(fp-fptt);	/* (x1 f1=fp), (x2 f2=fret), (x3 f3=fptt), (xm fm) /
if (t < 0.0) {	/* From x1 (P0) distance of x2 is at h and x3 is 2h */
linmin(p,xit,n,fret,func);	/* Let f"(x2) be the 2nd derivative equal everywhere. */
	/* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */
	/* will reach at f3 = fm + h^2/2 f"m ; f" = (f1 -2f2 +f3 ) / h*2 /
	/* f1-f3 = delta(2h) = 2 h*2 f'' = 2(f1- 2f2 +f3) /
	/* Thus we compare delta(2h) with observed f1-f3 */
	/* or best gain on one ancient line 'del' with total */
	/* gain f1-f2 = f1 - f2 - 'del' with del */
	/* t=2.0(fp-2.0(fret)+fptt)SQR(fp-(fret)-del)-delSQR(fp-fptt); */

	t=2.0(fp-2.0(fret)+fptt)SQR(fp-(*fret)-del);
	t= t- del*SQR(fp-fptt);
	printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0(fp-2.0(fret)+fptt)SQR(fp-(fret)-del),delSQR(fp-fptt),t);
	fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0(fp-2.0(fret)+fptt)SQR(fp-(fret)-del),delSQR(fp-fptt),t);
	#ifdef DEBUG
	printf("t3= %.12lf, t4= %.12lf, t3= %.12lf, t4= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
	(fp-(fret)-del)(fp-(fret)-del),(fp-fptt)(fp-fptt));
	fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3= %.12lf, t4= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
	(fp-(fret)-del)(fp-(fret)-del),(fp-fptt)(fp-fptt));
	printf("tt= %.12lf, t=%.12lf\n",2.0(fp-2.0(fret)+fptt)(fp-(fret)-del)(fp-(fret)-del)-del(fp-fptt)*(fp-fptt),t);
	fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0(fp-2.0(fret)+fptt)(fp-(fret)-del)(fp-(fret)-del)-del(fp-fptt)*(fp-fptt),t);
	#endif
	if (t < 0.0) { /* Then we use it for last direction */
	linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/
for (j=1;j<=n;j++) {	for (j=1;j<=n;j++) {
xi[j][ibig]=xi[j][n];	xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */
xi[j][n]=xit[j];	xi[j][n]=xit[j]; /* and nth direction by the extrapolated */
}	}
	printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
	fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction :\n",n,ibig);

#ifdef DEBUG	#ifdef DEBUG
printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);	printf("Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);	fprintf(ficlog,"Direction changed last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
Line 1368 void powell(double p[], double **xi, int	Line 1499 void powell(double p[], double **xi, int
printf("\n");	printf("\n");
fprintf(ficlog,"\n");	fprintf(ficlog,"\n");
#endif	#endif
}	} /* end of t negative */
}	} /* end if (fptt < fp) */
}	}
}	}

Line 1457 double pmij(double ps, double *cov,	Line 1588 double pmij(double ps, double *cov,
*/	*/
double s1, lnpijopii;	double s1, lnpijopii;
/double t34;/	/double t34;/
int i,j,j1, nc, ii, jj;	int i,j, 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++){
Line 1599 double *hpxij(double *po, int nhstep	Line 1730 double *hpxij(double *po, int nhstep
return po;	return po;
}	}

	#ifdef NLOPT
	double myfunc(unsigned n, const double p1, double grad, void *pd){
	double fret;
	double *xt;
	int j;
	myfunc_data d2 = (myfunc_data ) pd;
	/* xt = (p1-1); */
	xt=vector(1,n);
	for (j=1;j<=n;j++) xt[j]=p1[j-1]; /* xt[1]=p1[0] */

	fret=(d2->function)(xt); /* p xt[1]@8 is fine */
	/* fret=(func)(xt); /\ p xt[1]@8 is fine \/ /
	printf("Function = %.12lf ",fret);
	for (j=1;j<=n;j++) printf(" %d %.8lf", j, xt[j]);
	printf("\n");
	free_vector(xt,1,n);
	return fret;
	}
	#endif

/************* log-likelihood ***********/	/************* log-likelihood ***********/
double func( double *x)	double func( double *x)
Line 1617 double func( double *x)	Line 1767 double func( double *x)
/*for(i=1;i<imx;i++)	/*for(i=1;i<imx;i++)
printf(" %d\n",s[4][i]);	printf(" %d\n",s[4][i]);
*/	*/

	++countcallfunc;

cov[1]=1.;	cov[1]=1.;

for(k=1; k<=nlstate; k++) ll[k]=0.;	for(k=1; k<=nlstate; k++) ll[k]=0.;
Line 1998 void likelione(FILE *ficres,double p[],	Line 2151 void likelione(FILE *ficres,double p[],

void mlikeli(FILE ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (func)(double []))	void mlikeli(FILE ficres,double p[], int npar, int ncovmodel, int nlstate, double ftol, double (func)(double []))
{	{
int i,j, iter;	int i,j, iter=0;
double **xi;	double **xi;
double fret;	double fret;
double fretone; /* Only one call to likelihood */	double fretone; /* Only one call to likelihood */
/* char filerespow[FILENAMELENGTH];*/	/* char filerespow[FILENAMELENGTH];*/

	#ifdef NLOPT
	int creturn;
	nlopt_opt opt;
	/* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds \/ /
	double *lb;
	double minf; /* the minimum objective value, upon return */
	double * p1; /* Shifted parameters from 0 instead of 1 */
	myfunc_data dinst, *d = &dinst;
	#endif


xi=matrix(1,npar,1,npar);	xi=matrix(1,npar,1,npar);
for (i=1;i<=npar;i++)	for (i=1;i<=npar;i++)
for (j=1;j<=npar;j++)	for (j=1;j<=npar;j++)
Line 2019 void mlikeli(FILE *ficres,double p[], in	Line 2184 void mlikeli(FILE *ficres,double p[], in
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");
	#ifdef POWELL
powell(p,xi,npar,ftol,&iter,&fret,func);	powell(p,xi,npar,ftol,&iter,&fret,func);
	#endif

	#ifdef NLOPT
	#ifdef NEWUOA
	opt = nlopt_create(NLOPT_LN_NEWUOA,npar);
	#else
	opt = nlopt_create(NLOPT_LN_BOBYQA,npar);
	#endif
	lb=vector(0,npar-1);
	for (i=0;i<npar;i++) lb[i]= -HUGE_VAL;
	nlopt_set_lower_bounds(opt, lb);
	nlopt_set_initial_step1(opt, 0.1);

	p1= (p+1); /* p (p+1)@8 and p (p1)@8 are equal p1[0]=p[1] */
	d->function = func;
	printf(" Func %.12lf \n",myfunc(npar,p1,NULL,d));
	nlopt_set_min_objective(opt, myfunc, d);
	nlopt_set_xtol_rel(opt, ftol);
	if ((creturn=nlopt_optimize(opt, p1, &minf)) < 0) {
	printf("nlopt failed! %d\n",creturn);
	}
	else {
	printf("found minimum after %d evaluations (NLOPT=%d)\n", countcallfunc ,NLOPT);
	printf("found minimum at f(%g,%g) = %0.10g\n", p[0], p[1], minf);
	iter=1; /* not equal */
	}
	nlopt_destroy(opt);
	#endif
free_matrix(xi,1,npar,1,npar);	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 & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
fprintf(ficlog,"\n#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));	fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
fprintf(ficres,"#Number of iterations = %d, -2 Log likelihood = %.12f \n",iter,func(p));	fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));

}	}

Line 2035 void hesscov(double **matcov, double p[]	Line 2227 void hesscov(double **matcov, double p[]
{	{
double a,y,*x,pd;	double a,y,*x,pd;
double **hess;	double **hess;
int i, j,jk;	int i, j;
int *indx;	int *indx;

double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar);	double hessii(double p[], double delta, int theta, double delti[],double (*func)(double []),int npar);
Line 2169 double hessii(double x[], double delta,	Line 2361 double hessii(double x[], double delta,
k=kmax;	k=kmax;
}	}
else if((k1 >khi/nkhif) \|\| (k2 >khi/nkhif)){ /* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. */	else if((k1 >khi/nkhif) \|\| (k2 >khi/nkhif)){ /* Keeps lastvalue before 3.84/2 KHI2 5% 1d.f. */
k=kmax; l=lmax*10.;	k=kmax; l=lmax*10;
}	}
else if((k1 >khi/nkhi) \|\| (k2 >khi/nkhi)){	else if((k1 >khi/nkhi) \|\| (k2 >khi/nkhi)){
delts=delt;	delts=delt;
Line 2184 double hessii(double x[], double delta,	Line 2376 double hessii(double x[], double delta,
double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar)	double hessij( double x[], double delti[], int thetai,int thetaj,double (*func)(double []),int npar)
{	{
int i;	int i;
int l=1, l1, lmax=20;	int l=1, lmax=20;
double k1,k2,k3,k4,res,fx;	double k1,k2,k3,k4,res,fx;
double p2[MAXPARM+1];	double p2[MAXPARM+1];
int k;	int k;
Line 2299 void pstamp(FILE *fichier)	Line 2491 void pstamp(FILE *fichier)
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[])	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,j;	int i, m, jk, j1, bool, z1,j;
int first;	int first;
double **freq; / Frequencies */	double **freq; / Frequencies */
double pp, *prop;	double pp, *prop;
Line 2480 void prevalence(double ***probs, double	Line 2672 void prevalence(double ***probs, double
We still use firstpass and lastpass as another selection.	We still use firstpass and lastpass as another selection.
*/	*/

int i, m, jk, k1, i1, j1, bool, z1,j;	int i, m, jk, j1, bool, z1,j;
double **freq; / Frequencies */
double pp, *prop;	double **prop;
double pos,posprop;	double posprop;
double y2; /* in fractional years */	double y2; /* in fractional years */
int iagemin, iagemax;	int iagemin, iagemax;
int first; /** to stop verbosity which is redirected to log file */	int first; /** to stop verbosity which is redirected to log file */
Line 2574 void concatwav(int wav[], int **dh, int	Line 2766 void concatwav(int wav[], int **dh, int
int j, k=0,jk, ju, jl;	int j, k=0,jk, ju, jl;
double sum=0.;	double sum=0.;
first=0;	first=0;
jmin=1e+5;	jmin=100000;
jmax=-1;	jmax=-1;
jmean=0.;	jmean=0.;
for(i=1; i<=imx; i++){	for(i=1; i<=imx; i++){
Line 2814 void evsij(double ***eij, double x[], in	Line 3006 void evsij(double ***eij, double x[], in

{	{
/* Health expectancies, no variances */	/* Health expectancies, no variances */
int i, j, nhstepm, hstepm, h, nstepm, k, cptj, cptj2, i2, j2;	int i, j, nhstepm, hstepm, h, nstepm;
int nhstepma, nstepma; /* Decreasing with age */	int nhstepma, nstepma; /* Decreasing with age */
double age, agelim, hf;	double age, agelim, hf;
double ***p3mat;	double ***p3mat;
Line 3136 void varevsij(char optionfilefiname[], d	Line 3328 void varevsij(char optionfilefiname[], d
double dnewm,doldm;	double dnewm,doldm;
double dnewmp,doldmp;	double dnewmp,doldmp;
int i, j, nhstepm, hstepm, h, nstepm ;	int i, j, nhstepm, hstepm, h, nstepm ;
int k, cptcode;	int k;
double *xp;	double *xp;
double gp, gm; /* for var eij */	double gp, gm; /* for var eij */
double *gradg, trgradg; /for var eij */	double *gradg, trgradg; /for var eij */
Line 3436 void varprevlim(char fileres[], double *	Line 3628 void varprevlim(char fileres[], double *
{	{
/* 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);/
double **newm;
double dnewm,doldm;	double dnewm,doldm;
int i, j, nhstepm, hstepm;	int i, j, nhstepm, hstepm;
int k, cptcode;
double *xp;	double *xp;
double gp, gm;	double gp, gm;
double gradg, trgradg;	double gradg, trgradg;
Line 3518 void varprevlim(char fileres[], double *	Line 3709 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, char strstart[])	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, k1, l1, tj;
int k2, l2, j1, z1;	int k2, l2, j1, z1;
int k=0,l, cptcode;	int k=0, l;
int first=1, first1, first2;	int first=1, first1, first2;
double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;	double cv12, mu1, mu2, lc1, lc2, v12, v21, v11, v22,v1,v2, c12, tnalp;
double dnewm,doldm;	double dnewm,doldm;
Line 3528 void varprob(char optionfilefiname[], do	Line 3719 void varprob(char optionfilefiname[], do
double gp, gm;	double gp, gm;
double gradg, trgradg;	double gradg, trgradg;
double **mu;	double **mu;
double age,agelim, cov[NCOVMAX+1];	double age, cov[NCOVMAX+1];
double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */	double std=2.0; /* Number of standard deviation wide of confidence ellipsoids */
int theta;	int theta;
char fileresprob[FILENAMELENGTH];	char fileresprob[FILENAMELENGTH];
Line 3900 fprintf(fichtm," \n<ul><li><b>Graphs</b>	Line 4091 fprintf(fichtm," \n<ul><li><b>Graphs</b>
<img src=\"%s%d_2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);	<img src=\"%s%d_2.png\">",stepm,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1,subdirf2(optionfilefiname,"pe"),jj1);
/* 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>- Convergence from each state (1 to %d) to period (stable) prevalence in state %d <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \	fprintf(fichtm,"<br>- Convergence from cross-sectional prevalence in each state (1 to %d) to period (stable) prevalence in specific state %d <a href=\"%s%d_%d.png\">%s%d_%d.png</a><br> \
<img src=\"%s%d_%d.png\">",nlstate, cpt, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);	<img src=\"%s%d_%d.png\">",nlstate, cpt, subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1,subdirf2(optionfilefiname,"p"),cpt,jj1);
}	}
for(cpt=1; cpt<=nlstate;cpt++) {	for(cpt=1; cpt<=nlstate;cpt++) {
Line 3986 true period expectancies (those weighted	Line 4177 true period expectancies (those weighted
void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){	void printinggnuplot(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){

char dirfileres[132],optfileres[132];	char dirfileres[132],optfileres[132];
int m0,cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;	int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,ij=0,l=0;
int ng=0;	int ng=0;
/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */	/* if((ficgp=fopen(optionfilegnuplot,"a"))==NULL) { */
/* printf("Problem with file %s",optionfilegnuplot); */	/* printf("Problem with file %s",optionfilegnuplot); */
Line 4172 plot [%.f:%.f] ", ageminpar, agemaxpar)	Line 4363 plot [%.f:%.f] ", ageminpar, agemaxpar)
} /* end k2 */	} /* end k2 */
} /* end jk */	} /* end jk */
} /* end ng */	} /* end ng */
avoid:	/* avoid: */
fflush(ficgp);	fflush(ficgp);
} /* end gnuplot */	} /* end gnuplot */

Line 4226 prevforecast(char fileres[], double anpr	Line 4417 prevforecast(char fileres[], double anpr
dateprev1 dateprev2 range of dates during which prevalence is computed	dateprev1 dateprev2 range of dates during which prevalence is computed
anproj2 year of en of projection (same day and month as proj1).	anproj2 year of en of projection (same day and month as proj1).
*/	*/
int yearp, stepsize, hstepm, nhstepm, j, k, c, cptcod, i, h, i1;	int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1;
int *popage;
double agec; /* generic age */	double agec; /* generic age */
double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;	double agelim, ppij, yp,yp1,yp2,jprojmean,mprojmean,anprojmean;
double popeffectif,popcount;	double popeffectif,popcount;
Line 4521 void prwizard(int ncovmodel, int nlstate	Line 4711 void prwizard(int ncovmodel, int nlstate

/* Wizard to print covariance matrix template */	/* Wizard to print covariance matrix template */

char ca[32], cb[32], cc[32];	char ca[32], cb[32];
int i,j, k, l, li, lj, lk, ll, jj, npar, itimes;	int i,j, k, li, lj, lk, ll, jj, npar, itimes;
int numlinepar;	int numlinepar;

printf("# Parameters nlstatenlstatencov a121 + b12 age + ...\n");	printf("# Parameters nlstatenlstatencov a121 + b12 age + ...\n");
Line 4743 fprintf(fichtm,"<ul><li><h4>Life table</	Line 4933 fprintf(fichtm,"<ul><li><h4>Life table</
void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){	void printinggnuplotmort(char fileres[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , char pathc[], double p[]){

char dirfileres[132],optfileres[132];	char dirfileres[132],optfileres[132];
int m,cpt,k1,i,k,j,jk,k2,k3,ij,l;
int ng;	int ng;


Line 4768 int readdata(char datafile[], int firsto	Line 4958 int readdata(char datafile[], int firsto
/-------- data file ----------/	/-------- data file ----------/
FILE *fic;	FILE *fic;
char dummy[]=" ";	char dummy[]=" ";
int i, j, n;	int i=0, j=0, n=0;
int linei, month, year,iout;	int linei, month, year,iout;
char line[MAXLINE], linetmp[MAXLINE];	char line[MAXLINE], linetmp[MAXLINE];
char stra[80], strb[80];	char stra[MAXLINE], strb[MAXLINE];
char *stratrunc;	char *stratrunc;
int lstra;	int lstra;

Line 4799 int readdata(char datafile[], int firsto	Line 4989 int readdata(char datafile[], int firsto
continue;	continue;
}	}
trimbb(linetmp,line); /* Trims multiple blanks in line */	trimbb(linetmp,line); /* Trims multiple blanks in line */
for (j=0; line[j]!='\0';j++){	strcpy(line, linetmp);
line[j]=linetmp[j];
}


for (j=maxwav;j>=1;j--){	for (j=maxwav;j>=1;j--){
Line 4940 int readdata(char datafile[], int firsto	Line 5128 int readdata(char datafile[], int firsto
fclose(fic);	fclose(fic);

return (0);	return (0);
endread:	/* endread: */
printf("Exiting readdata: ");	printf("Exiting readdata: ");
fclose(fic);	fclose(fic);
return (1);	return (1);
Line 4971 int decodemodel ( char model[], int last	Line 5159 int decodemodel ( char model[], int last
*/	*/
{	{
int i, j, k, ks;	int i, j, k, ks;
int i1, j1, k1, k2;	int j1, k1, k2;
char modelsav[80];	char modelsav[80];
char stra[80], strb[80], strc[80], strd[80],stre[80];	char stra[80], strb[80], strc[80], strd[80],stre[80];

Line 5120 int decodemodel ( char model[], int last	Line 5308 int decodemodel ( char model[], int last


return (0); /* with covar[new additional covariate if product] and Tage if age */	return (0); /* with covar[new additional covariate if product] and Tage if age */
endread:	/endread:/
printf("Exiting decodemodel: ");	printf("Exiting decodemodel: ");
return (1);	return (1);
}	}
Line 5179 calandcheckages(int imx, int maxwav, dou	Line 5367 calandcheckages(int imx, int maxwav, dou
}	}
else if(agev[m][i] >*agemax){	else if(agev[m][i] >*agemax){
*agemax=agev[m][i];	*agemax=agev[m][i];
printf(" Max anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.2f\n",m,i,anint[m][i], i,annais[i], *agemax);	/* printf(" Max anint[%d][%d]=%.0f annais[%d]=%.0f, agemax=%.2f\n",m,i,anint[m][i], i,annais[i], agemax);/
}	}
/agev[m][i]=anint[m][i]-annais[i];/	/agev[m][i]=anint[m][i]-annais[i];/
/* agev[m][i] = age[i]+2m;/	/* agev[m][i] = age[i]+2m;/
Line 5217 calandcheckages(int imx, int maxwav, dou	Line 5405 calandcheckages(int imx, int maxwav, dou
fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);	fprintf(ficlog,"Total number of individuals= %d, Agemin = %.2f, Agemax= %.2f\n\n", imx, agemin, agemax);

return (0);	return (0);
endread:	/* endread:*/
printf("Exiting calandcheckages: ");	printf("Exiting calandcheckages: ");
return (1);	return (1);
}	}
Line 5237 int main(int argc, char *argv[])	Line 5425 int main(int argc, char *argv[])
double ssval;	double ssval;
#endif	#endif
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=0,m,size=100, cptcod;
int linei, month, year,iout;
int jj, ll, li, lj, lk, imk;	int jj, ll, li, lj, lk;
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;
int vpopbased=0;	int vpopbased=0;

char ca[32], cb[32], cc[32];	char ca[32], cb[32];
/* FILE fichtm; //* Html File */	/* FILE fichtm; //* Html File */
/* FILE ficgp;/ /Gnuplot File /	/* FILE ficgp;/ /Gnuplot File /
struct stat info;	struct stat info;
double agedeb, agefin,hf;	double agedeb;
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;

double fret;	double fret;
double **xi,tmp,delta;

double dum; /* Dummy variable */	double dum; /* Dummy variable */
double ***p3mat;	double ***p3mat;
double ***mobaverage;	double ***mobaverage;
int *indx;
char line[MAXLINE], linepar[MAXLINE];	char line[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], pathimach[MAXLINE];	char pathr[MAXLINE], pathimach[MAXLINE];
char *bp, tok, val; / pathtot */	char tok, val; /* pathtot */
int firstobs=1, lastobs=10;	int firstobs=1, lastobs=10;
int sdeb, sfin; /* Status at beginning and end */	int c, h , cpt;
int c, h , cpt,l;	int jl;
int ju,jl, mi;	int i1, j1, jk, stepsize;
int i1,j1, jk,aa,bb, stepsize, ij;	int *tab;
int jnais,jdc,jint4,jint1,jint2,jint3,*tab;
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;
Line 5277 int main(int argc, char *argv[])	Line 5462 int main(int argc, char *argv[])
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;

double bage, fage, age, agelim, agebase;	double bage=0, fage=110, age, agelim, agebase;
double ftolpl=FTOL;	double ftolpl=FTOL;
double **prlim;	double **prlim;
double **param; / Matrix of parameters */	double **param; / Matrix of parameters */
Line 5288 int main(int argc, char *argv[])	Line 5473 int main(int argc, char *argv[])
double *eij, *vareij;	double *eij, *vareij;
double *varpl; / Variances of prevalence limits by age */	double *varpl; / Variances of prevalence limits by age */
double *epj, vepp;	double *epj, vepp;
double kk1, kk2;
double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;	double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;
double **ximort;	double **ximort;
char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";	char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
int *dcwave;	int *dcwave;

char z[1]="c", occ;	char z[1]="c";

/char strt;*/	/char strt;*/
char strtend[80];	char strtend[80];

long total_usecs;

/* setlocale (LC_ALL, ""); */	/* setlocale (LC_ALL, ""); */
/* bindtextdomain (PACKAGE, LOCALEDIR); */	/* bindtextdomain (PACKAGE, LOCALEDIR); */
/* textdomain (PACKAGE); */	/* textdomain (PACKAGE); */
Line 5308 int main(int argc, char *argv[])	Line 5492 int main(int argc, char *argv[])
/* setlocale (LC_MESSAGES, ""); */	/* setlocale (LC_MESSAGES, ""); */

/* gettimeofday(&start_time, (struct timezone)0); / /* at first time */	/* gettimeofday(&start_time, (struct timezone)0); / /* at first time */
(void) gettimeofday(&start_time,&tzp);	rstart_time = time(NULL);
	/* (void) gettimeofday(&start_time,&tzp);*/
	start_time = *localtime(&rstart_time);
curr_time=start_time;	curr_time=start_time;
tm = *localtime(&start_time.tv_sec);	/tml = localtime(&start_time.tm_sec);*/
tmg = *gmtime(&start_time.tv_sec);	/* strcpy(strstart,asctime(&tml)); */
strcpy(strstart,asctime(&tm));	strcpy(strstart,asctime(&start_time));

/* printf("Localtime (at start)=%s",strstart); */	/* printf("Localtime (at start)=%s",strstart); */
/* tp.tv_sec = tp.tv_sec +86400; */	/* tp.tm_sec = tp.tm_sec +86400; */
/* tm = localtime(&start_time.tv_sec); /	/* tm = localtime(&start_time.tm_sec); /
/* tmg.tm_year=tmg.tm_year +dsigndyear; /	/* tmg.tm_year=tmg.tm_year +dsigndyear; /
/* tmg.tm_mon=tmg.tm_mon +dsigndmonth; /	/* tmg.tm_mon=tmg.tm_mon +dsigndmonth; /
/* tmg.tm_hour=tmg.tm_hour + 1; */	/* tmg.tm_hour=tmg.tm_hour + 1; */
/* tp.tv_sec = mktime(&tmg); */	/* tp.tm_sec = mktime(&tmg); */
/* strt=asctime(&tmg); */	/* strt=asctime(&tmg); */
/* printf("Time(after) =%s",strstart); */	/* printf("Time(after) =%s",strstart); */
/* (void) time (&time_value);	/* (void) time (&time_value);
Line 5341 int main(int argc, char *argv[])	Line 5527 int main(int argc, char *argv[])
i=strlen(pathr);	i=strlen(pathr);
if(pathr[i-1]=='\n')	if(pathr[i-1]=='\n')
pathr[i-1]='\0';	pathr[i-1]='\0';
	i=strlen(pathr);
	if(pathr[i-1]==' ') /* This may happen when dragging on oS/X! */
	pathr[i-1]='\0';
for (tok = pathr; tok != NULL; ){	for (tok = pathr; tok != NULL; ){
printf("Pathr \|%s\|\n",pathr);	printf("Pathr \|%s\|\n",pathr);
while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');	while ((val = strsep(&tok, "\"" )) != NULL && *val == '\0');
Line 5396 int main(int argc, char *argv[])	Line 5585 int main(int argc, char *argv[])
path=%s \n\	path=%s \n\
optionfile=%s\n\	optionfile=%s\n\
optionfilext=%s\n\	optionfilext=%s\n\
optionfilefiname=%s\n",pathimach,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);
fflush(ficlog);	fflush(ficlog);
/* (void) gettimeofday(&curr_time,&tzp); */	/* (void) gettimeofday(&curr_time,&tzp); */
/* printf("Elapsed time %d\n", asc_diff_time(curr_time.tv_sec-start_time.tv_sec,tmpout)); */	/* printf("Elapsed time %d\n", asc_diff_time(curr_time.tm_sec-start_time.tm_sec,tmpout)); */

/* */	/* */
strcpy(fileres,"r");	strcpy(fileres,"r");
Line 5566 run imach with mle=-1 to get a correct t	Line 5755 run imach with mle=-1 to get a correct t
for(i=1; i <=nlstate; i++){	for(i=1; i <=nlstate; i++){
for(j=1; j <=nlstate+ndeath-1; j++){	for(j=1; j <=nlstate+ndeath-1; j++){
fscanf(ficpar,"%1d%1d",&i1,&j1);	fscanf(ficpar,"%1d%1d",&i1,&j1);
if ((i1-i)*(j1-j)!=0){	if ( (i1-i) * (j1-j) != 0){
printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);	printf("Error in line parameters number %d, %1d%1d instead of %1d%1d \n",numlinepar, i,j, i1, j1);
exit(1);	exit(1);
}	}
Line 5935 Interval (in months) between two waves:	Line 6124 Interval (in months) between two waves:

#ifdef GSL	#ifdef GSL
printf("GSL optimization\n"); fprintf(ficlog,"Powell\n");	printf("GSL optimization\n"); fprintf(ficlog,"Powell\n");
#elsedef	#else
printf("Powell\n"); fprintf(ficlog,"Powell\n");	printf("Powell\n"); fprintf(ficlog,"Powell\n");
#endif	#endif
strcpy(filerespow,"pow-mort");	strcpy(filerespow,"pow-mort");
Line 5946 Interval (in months) between two waves:	Line 6135 Interval (in months) between two waves:
}	}
#ifdef GSL	#ifdef GSL
fprintf(ficrespow,"# GSL optimization\n# iter -2*LL");	fprintf(ficrespow,"# GSL optimization\n# iter -2*LL");
#elsedef	#else
fprintf(ficrespow,"# Powell\n# iter -2*LL");	fprintf(ficrespow,"# Powell\n# iter -2*LL");
#endif	#endif
/* for (i=1;i<=nlstate;i++)	/* for (i=1;i<=nlstate;i++)
Line 6508 Interval (in months) between two waves:	Line 6697 Interval (in months) between two waves:

for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/	for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
oldm=oldms;savm=savms; /* Segmentation fault */	oldm=oldms;savm=savms; /* Segmentation fault */
varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart);	cptcod= 0; /* To be deleted */
	varevsij(optionfilefiname, vareij, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl,k, estepm, cptcov,cptcod,vpopbased,mobilav, strstart); /* cptcod not initialized Intel */
fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");	fprintf(ficrest,"# Total life expectancy with std error and decomposition into time to be expected in each health state\n# (weighted average of eij where weights are ");
if(vpopbased==1)	if(vpopbased==1)
fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);	fprintf(ficrest,"the age specific prevalence observed (cross-sectionally) in the population i.e cross-sectionally\n in each health state (popbased=1) (mobilav=%d)\n",mobilav);
Line 6600 Interval (in months) between two waves:	Line 6790 Interval (in months) between two waves:
if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);	if (mobilav!=0) free_ma3x(mobaverage,1, AGESUP,1,NCOVMAX, 1,NCOVMAX);
free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);	free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
} /* mle==-3 arrives here for freeing */	} /* mle==-3 arrives here for freeing */
endfree:	/* endfree:*/
free_matrix(prlim,1,nlstate,1,nlstate); /here or after loop ? /	free_matrix(prlim,1,nlstate,1,nlstate); /here or after loop ? /
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
Line 6634 Interval (in months) between two waves:	Line 6824 Interval (in months) between two waves:
}	}
printf("See log file on %s\n",filelog);	printf("See log file on %s\n",filelog);
/* gettimeofday(&end_time, (struct timezone)0);/ /* after time */	/* gettimeofday(&end_time, (struct timezone)0);/ /* after time */
(void) gettimeofday(&end_time,&tzp);	/(void) gettimeofday(&end_time,&tzp);/
tm = *localtime(&end_time.tv_sec);	rend_time = time(NULL);
tmg = *gmtime(&end_time.tv_sec);	end_time = *localtime(&rend_time);
strcpy(strtend,asctime(&tm));	/* tml = localtime(&end_time.tm_sec); /
	strcpy(strtend,asctime(&end_time));
printf("Local time at start %s\nLocal time 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(rend_time -rstart_time,tmpout));

printf("Total time was %ld Sec.\n", end_time.tv_sec -start_time.tv_sec);	printf("Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time));
fprintf(ficlog,"Total time used %s\n", asc_diff_time(end_time.tv_sec -start_time.tv_sec,tmpout));	fprintf(ficlog,"Total time used %s\n", asc_diff_time(rend_time -rstart_time,tmpout));
fprintf(ficlog,"Total time was %ld Sec.\n", end_time.tv_sec -start_time.tv_sec);	fprintf(ficlog,"Total time was %.0lf Sec.\n", difftime(rend_time,rstart_time));
/* printf("Total time was %d uSec.\n", total_usecs);*/	/* printf("Total time was %d uSec.\n", total_usecs);*/
/* if(fileappend(fichtm,optionfilehtm)){ */	/* if(fileappend(fichtm,optionfilehtm)){ */
fprintf(fichtm,"<br>Local time at start %s<br>Local time at end %s<br>\n</body></html>",strstart, strtend);	fprintf(fichtm,"<br>Local time at start %s<br>Local time at end %s<br>\n</body></html>",strstart, strtend);
Line 6663 Interval (in months) between two waves:	Line 6854 Interval (in months) between two waves:
printf("Current directory %s!\n",pathcd);	printf("Current directory %s!\n",pathcd);
/strcat(plotcmd,CHARSEPARATOR);/	/strcat(plotcmd,CHARSEPARATOR);/
sprintf(plotcmd,"gnuplot");	sprintf(plotcmd,"gnuplot");
#ifndef UNIX	#ifdef _WIN32
sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);	sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);
#endif	#endif
if(!stat(plotcmd,&info)){	if(!stat(plotcmd,&info)){
printf("Error or gnuplot program not found: %s\n",plotcmd);fflush(stdout);	printf("Error or gnuplot program not found: '%s'\n",plotcmd);fflush(stdout);
if(!stat(getenv("GNUPLOTBIN"),&info)){	if(!stat(getenv("GNUPLOTBIN"),&info)){
printf("Error or gnuplot program not found: %s Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout);	printf("Error or gnuplot program not found: '%s' Environment GNUPLOTBIN not set.\n",plotcmd);fflush(stdout);
}else	}else
strcpy(pplotcmd,plotcmd);	strcpy(pplotcmd,plotcmd);
#ifdef UNIX	#ifdef __unix
strcpy(plotcmd,GNUPLOTPROGRAM);	strcpy(plotcmd,GNUPLOTPROGRAM);
if(!stat(plotcmd,&info)){	if(!stat(plotcmd,&info)){
printf("Error gnuplot program not found: %s\n",plotcmd);fflush(stdout);	printf("Error gnuplot program not found: '%s'\n",plotcmd);fflush(stdout);
}else	}else
strcpy(pplotcmd,plotcmd);	strcpy(pplotcmd,plotcmd);
#endif	#endif
Line 6683 Interval (in months) between two waves:	Line 6874 Interval (in months) between two waves:
strcpy(pplotcmd,plotcmd);	strcpy(pplotcmd,plotcmd);

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

if((outcmd=system(plotcmd)) != 0){	if((outcmd=system(plotcmd)) != 0){
printf("gnuplot command might not be in your path: %s, err=%d\n", plotcmd, outcmd);	printf("gnuplot command might not be in your path: '%s', err=%d\n", plotcmd, outcmd);
printf("\n Trying if gnuplot resides on the same directory that IMaCh\n");	printf("\n Trying if gnuplot resides on the same directory that IMaCh\n");
sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot);	sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot);
if((outcmd=system(plotcmd)) != 0)	if((outcmd=system(plotcmd)) != 0)
printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd);	printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd);
}	}
printf(" Successul, please wait...");	printf(" Successful, please wait...");
while (z[0] != 'q') {	while (z[0] != 'q') {
/* chdir(path); */	/* chdir(path); */
printf("\nType e to edit results with your browser, g to graph again and q for exit: ");	printf("\nType e to edit results with your browser, g to graph again and q for exit: ");
scanf("%s",z);	scanf("%s",z);
/* if (z[0] == 'c') system("./imach"); */	/* if (z[0] == 'c') system("./imach"); */
if (z[0] == 'e') {	if (z[0] == 'e') {
#ifdef OSX	#ifdef __APPLE__
sprintf(pplotcmd, "open %s", optionfilehtm);	sprintf(pplotcmd, "open %s", optionfilehtm);
	#elif __linux
	sprintf(pplotcmd, "xdg-open %s", optionfilehtm);
#else	#else
sprintf(pplotcmd, "%s", optionfilehtm);	sprintf(pplotcmd, "%s", optionfilehtm);
#endif	#endif
Line 6716 Interval (in months) between two waves:	Line 6909 Interval (in months) between two waves:
scanf("%s",z);	scanf("%s",z);
}	}
}	}

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

Removed from v.1.155
changed lines
	Added in v.1.165