imach/src/imachprax.c - diff

Return to imachprax.c CVS log

Up to [Local Repository] / imach / src

Diff for /imach/src/imachprax.c between versions 1.4 and 1.5

version 1.4, 2023/06/22 12:50:51	version 1.5, 2023/10/09 09:10:01
Line 1	Line 1
/* $Id$	/* $Id$
$State$	$State$
$Log$	$Log$
	Revision 1.5 2023/10/09 09:10:01 brouard
	Summary: trying to reconsider

Revision 1.4 2023/06/22 12:50:51 brouard	Revision 1.4 2023/06/22 12:50:51 brouard
Summary: stil on going	Summary: stil on going

Line 1288 Important routines	Line 1291 Important routines
/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix \/ /	/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix \/ /
/* #define FLATSUP // Suppresses directions where likelihood is flat */	/* #define FLATSUP // Suppresses directions where likelihood is flat */
/* #define POWELLORIGINCONJUGATE /\* Don't use conjugate but biggest decrease if valuable \/ /	/* #define POWELLORIGINCONJUGATE /\* Don't use conjugate but biggest decrease if valuable \/ /
	/* #define NOTMINFIT */

#include <math.h>	#include <math.h>
#include <stdio.h>	#include <stdio.h>
Line 2620 void linmin(double p[], double xi[], int	Line 2624 void linmin(double p[], double xi[], int
}	}

/** praxis **/	/** praxis **/
	# include <float.h>

void transpose_in_place ( int n, double **a )	void transpose_in_place ( int n, double **a )

Line 2694 double pythag( double x, double y )	Line 2699 double pythag( double x, double y )
return value;	return value;
}	}

	void svsort ( int n, double d[], double **v )

	/******************************************************************************/
	/*
	Purpose:

	SVSORT descending sorts D and adjusts the corresponding columns of V.

	Discussion:
	A simple bubble sort is used on D.
	In our application, D contains singular values, and the columns of V are
	the corresponding right singular vectors.
	Author:
	Original FORTRAN77 version by Richard Brent.
	Richard Brent,
	Algorithms for Minimization with Derivatives,
	Prentice Hall, 1973,
	Reprinted by Dover, 2002.

	Parameters:
	Input, int N, the length of D, and the order of V.
	Input/output, double D[N], the vector to be sorted.
	On output, the entries of D are in descending order.

	Input/output, double V[N,N], an N by N array to be adjusted
	as D is sorted. In particular, if the value that was in D(I) on input is
	moved to D(J) on output, then the input column V(*,I) is moved to
	the output column V(*,J).
	*/
	{
	int i, j1, j2, j3;
	double t;

	for (j1 = 1; j1 < n; j1++) {
	/*
	* Find J3, the index of the largest entry in D[J1:N-1].
	* MAXLOC apparently requires its output to be an array.
	*/
	j3 = j1;
	for (j2 = j1+1; j2 < n; j2++) {
	if (d[j3] < d[j2]) {
	j3 = j2;
	}
	}
	/*
	* If J1 != J3, swap D[J1] and D[J3], and columns J1 and J3 of V.
	*/
	if (j1 != j3) {
	t = d[j1];
	d[j1] = d[j3];
	d[j3] = t;
	for (i = 1; i <= n; i++) {
	t = v[i][j1];
	v[i][j1] = v[i][j3];
	v[i][j3] = t;
	} /* end i */
	} /* end j1 != j3 */
	} /* end j1 */
	return;
	}

/* void svdcmp(double a, int m, int n, double w[], double v) */	/* void svdcmp(double a, int m, int n, double w[], double v) */
void svdminfit(double **a, int m, int n, double w[])	void svdminfit(double **a, int m, int n, double w[])
{	{
Line 2967 void powell(double p[], double **xi, int	Line 3033 void powell(double p[], double **xi, int
curr_time = *localtime(&rcurr_time);	curr_time = *localtime(&rcurr_time);
/* printf("\nPowell iter=%d -2LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",iter,fret, fp-fret,fp-fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); /	/* printf("\nPowell iter=%d -2LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",iter,fret, fp-fret,fp-fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout); /
/* fprintf(ficlog,"\nPowell iter=%d -2LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",iter,fret, fp-fret,fp-fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); /	/* fprintf(ficlog,"\nPowell iter=%d -2LL=%.12f gain=%.12f=%.3g %ld sec. %ld sec.",iter,fret, fp-fret,fp-fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog); /
Bigter=(iter - iter % ncovmodel)/ncovmodel +1; /* Big iteration, i.e on ncovmodel cycle */	/* Bigter=(iter - iter % ncovmodel)/ncovmodel +1; /\* Big iteration, i.e on ncovmodel cycle \/ /
	Bigter=(iter - (iter-1) % n)/n +1; /* Big iteration, i.e on ncovmodel cycle */
printf("\nPowell iter=%d Big Iter=%d -2LL=%.12f gain=%.3lg %ld sec. %ld sec.",iter,Bigter,fret,fp-fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);	printf("\nPowell iter=%d Big Iter=%d -2LL=%.12f gain=%.3lg %ld sec. %ld sec.",iter,Bigter,fret,fp-fret, rcurr_time-rlast_time, rcurr_time-rstart_time);fflush(stdout);
fprintf(ficlog,"\nPowell iter=%d Big Iter=%d -2LL=%.12f gain=%.3lg %ld sec. %ld sec.",iter,Bigter,fret,fp-fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);	fprintf(ficlog,"\nPowell iter=%d Big Iter=%d -2LL=%.12f gain=%.3lg %ld sec. %ld sec.",iter,Bigter,fret,fp-fret,rcurr_time-rlast_time, rcurr_time-rstart_time); fflush(ficlog);
fprintf(ficrespow,"%d %d %.12f %d",iter,Bigter, fret,curr_time.tm_sec-start_time.tm_sec);	fprintf(ficrespow,"%d %d %.12f %d",iter,Bigter, fret,curr_time.tm_sec-start_time.tm_sec);
Line 3038 void powell(double p[], double **xi, int	Line 3105 void powell(double p[], double **xi, int
fprintf(ficlog," - if your program needs %d BIG iterations (%d iterations) to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);	fprintf(ficlog," - if your program needs %d BIG iterations (%d iterations) to converge, convergence will be \n reached in %s i.e.\n on %s (current time is %s);\n",nBigterf, niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
}	}
}	}
for (i=1;i<=n;i++) { /* For each direction i */	for (i=1;i<=n;i++) { /* For each direction i, maximisation after loading directions */
for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */	for (j=1;j<=n;j++) xit[j]=xi[j][i]; /* Directions stored from previous iteration with previous scales */
fptt=(*fret);	fptt=(*fret);
#ifdef DEBUG	#ifdef DEBUG
Line 3048 void powell(double p[], double **xi, int	Line 3115 void powell(double p[], double **xi, int
printf("%d",i);fflush(stdout); /* print direction (parameter) i */	printf("%d",i);fflush(stdout); /* print direction (parameter) i */
fprintf(ficlog,"%d",i);fflush(ficlog);	fprintf(ficlog,"%d",i);fflush(ficlog);
#ifdef LINMINORIGINAL	#ifdef LINMINORIGINAL
linmin(p,xit,n,fret,func); /* New point i minimizing in direction i has coordinates p[j].*/	linmin(p,xit,n,fret,func); /* New point i minimizing in direction xit i has coordinates p[j].*/
/* xit[j] gives the n coordinates of direction i as input.*/	/* xit[j] gives the n coordinates of direction i as input.*/
/* fret gives the maximum value on direction xit /	/* fret gives the maximum value on direction xit /
#else	#else
linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/	linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
flatdir[i]=flat; /* Function is vanishing in that direction i */	flatdir[i]=flat; /* Function is vanishing in that direction i */
#endif	#endif
/* Outputs are fret(new point p) p is updated and xit rescaled */	/* Outputs are fret(new point p) p is updated and xit rescaled */
if (fabs(fptt-(fret)) > del) { / We are keeping the max gain on each of the n directions */	if (fabs(fptt-(fret)) > del) { / We are keeping the max gain on each of the n directions */
/* because that direction will be replaced unless the gain del is small */	/* because that direction will be replaced unless the gain del is small */
/* in comparison with the 'probable' gain, mu^2, with the last average direction. */	/* in comparison with the 'probable' gain, mu^2, with the last average direction. */
/* Unless the n directions are conjugate some gain in the determinant may be obtained */	/* Unless the n directions are conjugate some gain in the determinant may be obtained */
/* with the new direction. */	/* with the new direction. */
del=fabs(fptt-(*fret));	del=fabs(fptt-(*fret));
ibig=i;	ibig=i;
}	}
#ifdef DEBUG	#ifdef DEBUG
printf("%d %.12e",i,(*fret));	printf("%d %.12e",i,(*fret));
fprintf(ficlog,"%d %.12e",i,(*fret));	fprintf(ficlog,"%d %.12e",i,(*fret));
for (j=1;j<=n;j++) {	for (j=1;j<=n;j++) {
xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);	xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
printf(" x(%d)=%.12e",j,xit[j]);	printf(" x(%d)=%.12e",j,xit[j]);
fprintf(ficlog," x(%d)=%.12e",j,xit[j]);	fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
}	}
for(j=1;j<=n;j++) {	for(j=1;j<=n;j++) {
printf(" p(%d)=%.12e",j,p[j]);	printf(" p(%d)=%.12e",j,p[j]);
fprintf(ficlog," p(%d)=%.12e",j,p[j]);	fprintf(ficlog," p(%d)=%.12e",j,p[j]);
}	}
printf("\n");	printf("\n");
fprintf(ficlog,"\n");	fprintf(ficlog,"\n");
Line 3104 void powell(double p[], double **xi, int	Line 3171 void powell(double p[], double **xi, int
/* the scales of the directions and the directions, because the are reset to canonical directions */	/* the scales of the directions and the directions, because the are reset to canonical directions */
/* Thus the first calls to linmin will give new points and better maximizations until fp-(fret) is /	/* Thus the first calls to linmin will give new points and better maximizations until fp-(fret) is /
/* under the tolerance value. If the tolerance is very small 1.e-9, it could last long. */	/* under the tolerance value. If the tolerance is very small 1.e-9, it could last long. */
#ifdef DEBUG
int k[2],l;
k[0]=1;
k[1]=-1;
printf("Max: %.12e",(*func)(p));
fprintf(ficlog,"Max: %.12e",(*func)(p));
for (j=1;j<=n;j++) {
printf(" %.12e",p[j]);
fprintf(ficlog," %.12e",p[j]);
}
printf("\n");
fprintf(ficlog,"\n");
for(l=0;l<=1;l++) {
for (j=1;j<=n;j++) {
ptt[j]=p[j]+(p[j]-pt[j])*k[l];
printf("l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
fprintf(ficlog,"l=%d j=%d ptt=%.12e, xits=%.12e, p=%.12e, xit=%.12e", l,j,ptt[j],xits[j],p[j],xit[j]);
}
printf("func(ptt)=%.12e, deriv=%.12e\n",(func)(ptt),(ptt[j]-p[j])/((func)(ptt)-(*func)(p)));
fprintf(ficlog,"func(ptt)=%.12e, deriv=%.12e\n",(func)(ptt),(ptt[j]-p[j])/((func)(ptt)-(*func)(p)));
}
#endif

free_vector(xit,1,n);	free_vector(xit,1,n);
free_vector(xits,1,n);	free_vector(xits,1,n);
Line 3134 void powell(double p[], double **xi, int	Line 3179 void powell(double p[], double **xi, int
return;	return;
} /* enough precision */	} /* enough precision */
if (iter == ITMAXn) nrerror("powell exceeding maximum iterations.");	if (iter == ITMAXn) nrerror("powell exceeding maximum iterations.");
for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */	for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0)=2Pn-P0 */
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]; /* Coordinate j of last direction xi_n=P_n-_0 */
	printf("\n %d xit=%12.7g p=%12.7g pt=%12.7g ",j,xit[j],p[j],pt[j]);
pt[j]=p[j];	pt[j]=p[j];
}	}
	printf("\n");
fptt=(func)(ptt); / f_3 */	fptt=(func)(ptt); / f_3 */
#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */	#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */
if (*iter <=4) {	if (*iter <=4) {
Line 3159 void powell(double p[], double **xi, int	Line 3206 void powell(double p[], double **xi, int
/* 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)-delSQR(fp-fptt); */
/* Even if f3 <f1, directest can be negative and t >0 */	/* Even if f3 <f1, directest can be negative and t >0 */
/* mu² and del² are equal when f3=f1 */	/* mu² and del² are equal when f3=f1 */
/* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */	/* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */
/* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */	/* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */
/* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0 */	/* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0 */
/* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */	/* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */
#ifdef NRCORIGINAL	#ifdef NRCORIGINAL
t=2.0(fp-2.0(fret)+fptt)SQR(fp-(fret)-del)- delSQR(fp-fptt); /* Original Numerical Recipes in C*/	t=2.0(fp-2.0(fret)+fptt)SQR(fp-(fret)-del)- delSQR(fp-fptt); /* Original Numerical Recipes in C*/
#else	#else
Line 3170 void powell(double p[], double **xi, int	Line 3217 void powell(double p[], double **xi, int
t= t- del*SQR(fp-fptt);	t= t- del*SQR(fp-fptt);
#endif	#endif
directest = fp-2.0(fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */	directest = fp-2.0(fret)+fptt - 2.0 * del; /* If delta was big enough we change it for a new direction */
#ifdef DEBUG	printf(" t=%g, directest=%g\n",t, directest);
printf("t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0(fp-2.0(fret)+fptt)SQR(fp-(fret)-del),delSQR(fp-fptt),t,directest);	#ifdef POWELLNOTTRUECONJUGATE /* Searching for IBIG and testing for replacement */
fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0(fp-2.0(fret)+fptt)SQR(fp-(fret)-del),delSQR(fp-fptt),t,directest);
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
#ifdef POWELLORIGINAL	#ifdef POWELLORIGINAL
if (t < 0.0) { /* Then we use it for new direction */	if (t < 0.0) { /* Then we use it for new direction */
#else	#else /* Not POWELLOriginal but Brouard's */
if (directestt < 0.0) { / Contradiction between both tests */	if (directestt < 0.0) { / Contradiction between both tests */
printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);	printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0(fret)+fptt, fp -(*fret) -del, fp-fptt);	printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0(fret)+fptt, fp -(*fret) -del, fp-fptt);
fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);	fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0(fret)+fptt, fp -(*fret) -del, fp-fptt);	fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0(fret)+fptt, fp -(*fret) -del, fp-fptt);
}	}
if (directest < 0.0) { /* Then we use it for new direction */	if (directest < 0.0) { /* Then we use (P0, Pn) for new direction Xi_n or Xi_iBig */
#endif	#endif /* end POWELLOriginal */
#ifdef DEBUGLINMIN	#endif /* POWELLNOTTRUECONJUGATE else means systematic replacement by new direction P_0P_n */
printf("Before linmin in direction P%d-P0\n",n);
for (j=1;j<=n;j++) {
printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
if(j % ncovmodel == 0){
printf("\n");
fprintf(ficlog,"\n");
}
}
#endif
#ifdef LINMINORIGINAL	#ifdef LINMINORIGINAL
linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/	/* xit[j]=p[j]-pt[j] */
#else	printf("\n Computes min on P_0, P_n direction iter=%d Bigter=%d\n",*iter,Bigter);
	linmin(p,xit,n,fret,func); /* computes minimum on P_0,P_n direction: changes p and rescales xit.*/
	#else /* NOT LINMINORIGINAL but with searching for flat directions */
	printf("\n Flat Computes min on P_0, P_n direction iter=%d Bigter=%d\n",*iter,Bigter);
linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/	linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
flatdir[i]=flat; /* Function is vanishing in that direction i */	flatdir[i]=flat; /* Function is vanishing in that direction i */
#endif	#endif

#ifdef DEBUGLINMIN	#ifdef POWELLNOTTRUECONJUGATE
for (j=1;j<=n;j++) {	#else
printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
if(j % ncovmodel == 0){
printf("\n");
fprintf(ficlog,"\n");
}
}
#endif
#ifdef POWELLORIGINCONJUGATE	#ifdef POWELLORIGINCONJUGATE
for (j=1;j<=n;j++) {	for (j=1;j<=n;j++) {
xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */	xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */	xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */
}	}
#else	#else
for (j=1;j<=n-1;j++) {	for (i=1;i<=n-1;i++) {
xi[j][1]=xi[j][j+1]; /* Standard method of conjugate directions */	for (j=1;j<=n;j++) {
	xi[j][i]=xi[j][i+1]; /* Standard method of conjugate directions, not Powell who changes the nth direction by p0 pn . */
	}
	}
	for (j=1;j<=n;j++) {
xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */	xi[j][n]=xit[j]; /* and this nth direction by the by the average p_0 p_n */
}	}
/*	#endif /* POWELLORIGINCONJUGATE*/
Calculate a new set of orthogonal directions before repeating	#endif /POWELLNOTTRUECONJUGATE/
the main loop.	printf(" Standard method of conjugate directions\n");
Transpose V for SVD (minfit) (because minfit returns the right V in ULV=A):	printf("\n#A Before prax Bigter=%d model= 1 + age ", Bigter);
*/	for(j=1;j<=n;j++){
printf(" Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n");	printf("%d \n",j);
transpose_in_place ( n, xi );	for(i=1;i<=n;i++){
/*	printf(" %f",xi[j][i]);
SVD/MINFIT finds the singular value decomposition of V.	}
	}
	printf("\n");

This gives the principal values and principal directions of the	#ifdef NOTMINFIT
approximating quadratic form without squaring the condition number.	#else
*/	if(*iter >n){
printf(" SVDMINFIT finds the singular value decomposition of V. \n This gives the principal values and principal directions of the\n approximating quadratic form without squaring the condition number...\n");	/* if(Bigter >n){ */
double w; / eigenvalues of principal directions */	printf("\n#Before prax Bigter=%d model= 1 + age ", Bigter);
w=vector(1,n);	printf("\n");
svdminfit (xi, n, n, w ); /* In Brent's notation find d such that V=Q Diagonal(d) R, and Lambda=d^(-1/2) */	for(j=1;j<=n;j++){
/* minfit ( n, vsmall, v, d ); */	printf("%d \n",j);
/* v is overwritten with R. */	for(i=1;i<=n;i++){
free_vector(w,1,n);	printf(" %f",xi[j][i]);
#endif	}
	}
	printf("\n");
	/*
	* Calculate a new set of orthogonal directions before repeating
	* the main loop.
	* Transpose V for SVD (minfit) (because minfit returns the right V in ULV=A):
	*/
	printf(" Bigter=%d Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n",Bigter);
	transpose_in_place ( n, xi );
	/*
	SVD/MINFIT finds the singular value decomposition of V.

	This gives the principal values and principal directions of the
	approximating quadratic form without squaring the condition number.
	*/
	printf(" SVDMINFIT finds the singular value decomposition of V. \n This gives the principal values and principal directions of the\n approximating quadratic form without squaring the condition number...\n");
	double d; / eigenvalues of principal directions */
	d=vector(1,n);


	svdminfit (xi, n, n, d ); /* In Brent's notation find d such that V=Q Diagonal(d) R, and Lambda=d^(-1/2) */

	printf("\n#After prax model= 1 + age ");
	fprintf(ficlog,"\n#model= 1 + age ");

	if(nagesqr==1){
	printf(" + age*age ");
	fprintf(ficlog," + age*age ");
	}
	for(j=1;j <=ncovmodel-2;j++){
	if(Typevar[j]==0) {
	printf(" + V%d ",Tvar[j]);
	fprintf(ficlog," + V%d ",Tvar[j]);
	}else if(Typevar[j]==1) {
	printf(" + V%d*age ",Tvar[j]);
	fprintf(ficlog," + V%d*age ",Tvar[j]);
	}else if(Typevar[j]==2) {
	printf(" + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
	fprintf(ficlog," + V%d*V%d ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
	}else if(Typevar[j]==3) {
	printf(" + V%dV%dage ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
	fprintf(ficlog," + V%dV%dage ",Tvard[Tposprod[j]][1],Tvard[Tposprod[j]][2]);
	}
	}
	printf("\n");
	fprintf(ficlog,"\n");
	for(i=1,jk=1; i <=nlstate; i++){
	for(k=1; k <=(nlstate+ndeath); k++){
	if (k != i) {
	printf("%d%d ",i,k);
	fprintf(ficlog,"%d%d ",i,k);
	for(j=1; j <=ncovmodel; j++){
	printf("%12.7f ",p[jk]);
	fprintf(ficlog,"%12.7f ",p[jk]);
	jk++;
	}
	printf("\n");
	fprintf(ficlog,"\n");
	}
	}
	}
	/* minfit ( n, vsmall, v, d ); */
	/* v is overwritten with R. */
	/*
	Heuristic numbers:
	If the axes may be badly scaled (which is to be avoided if
	possible), then set SCBD = 10. Otherwise set SCBD = 1.

	If the problem is known to be ill-conditioned, initialize ILLC = true.
	KTM is the number of iterations without improvement before the
	algorithm terminates. KTM = 4 is very cautious; usually KTM = 1
	is satisfactory.
	*/
	double machep, small;
	double dmin;
	int illc=0; /* Local, int ILLC, is TRUE if the system is ill-conditioned. */
	machep = DBL_EPSILON;
	small = machep * machep;
	/* m2 = dsqrt(machep); */

	/*
	* Sort the eigenvalues and eigenvectors.
	*/
	printf(" Sort the eigenvalues and eigenvectors....\n");
	svsort ( n, d, xi );
	printf("Sorted Eigenvalues:\n");
	for(i=1; i<=n;i++){
	printf(" d[%d]=%g",i,d[i]);
	}
	printf("\n");
	/*
	* Determine the smallest eigenvalue.
	*/
	printf(" Determine the smallest eigenvalue.\n");
	dmin = fmax ( d[n], small );
	/*
	* The ratio of the smallest to largest eigenvalue determines whether
	* the system is ill conditioned.
	*/
	if ( dmin < sqrt(machep) * d[1] ) { /* m2d[0] /
	illc = 1;
	} else {
	illc = 0;
	}
	printf(" The ratio of the smallest to largest eigenvalue determines whether\n \
	the system is ill conditioned=%d . dmin=%.12lf < m2=%.12lf * d[1]=%.12lf \n",illc, dmin,sqrt(machep), d[1]);
	/* if ( 1.0 < scbd ) { */
	/* r8vec_print ( n, z, " The scale factors:" ); */
	/* } */
	/* r8vec_print ( n, d, " Principal values of the quadratic form:" ); */
	/* } */
	/* if ( 3 < prin ) { */
	/* r8mat_print ( n, n, v, " The principal axes:" ); */
	/* } */
	free_vector(d,1,n);
	/*
	* The main loop ends here.
	*/

	/* if ( 0 < prin ) */
	/* { */
	/* r8vec_print ( n, x, " X:" ); */
	/* } */
	}
	#endif /* NOTMINFIT */
#ifdef LINMINORIGINAL	#ifdef LINMINORIGINAL
#else	#else
for (j=1, flatd=0;j<=n;j++) {	for (j=1, flatd=0;j<=n;j++) {
Line 3274 void powell(double p[], double **xi, int	Line 3429 void powell(double p[], double **xi, int
free_vector(pt,1,n);	free_vector(pt,1,n);
return;	return;
#endif	#endif
}	} /* endif(flatd >0) */
#endif	#endif
printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);	printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);	fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
Line 3286 void powell(double p[], double **xi, int	Line 3441 void powell(double p[], double **xi, int
/* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$. As $P^{(2)}_1$ and */	/* Minimization amongst $\xi^{(1)}_2=(P_2-P_0)$ gives point $P^{(1)}_2$. As $P^{(2)}_1$ and */
/* $P^{(1)}_0$ are minimizing in the same direction $P^{(1)}_2 - P^{(1)}_1= P_2-P_0$, directions $P^{(1)}_2-P^{(1)}_0$ */	/* $P^{(1)}_0$ are minimizing in the same direction $P^{(1)}_2 - P^{(1)}_1= P_2-P_0$, directions $P^{(1)}_2-P^{(1)}_0$ */
/* and $P_2-P_0$ (parallel to $\xi$ and $\xi^c$) are conjugate. } */	/* and $P_2-P_0$ (parallel to $\xi$ and $\xi^c$) are conjugate. } */
	#ifdef POWELLNOTTRUECONJUGATE

#ifdef DEBUG
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);
for(j=1;j<=n;j++){
printf(" %lf",xit[j]);
fprintf(ficlog," %lf",xit[j]);
}
printf("\n");
fprintf(ficlog,"\n");
#endif
} /* end of t or directest negative */	} /* end of t or directest negative */
	#endif
#ifdef POWELLNOF3INFF1TEST	#ifdef POWELLNOF3INFF1TEST
#else	#else
} /* end if (fptt < fp) */	} /* end if (fptt < fp) */
#endif	#endif
	#ifdef POWELLORIGINCONJUGATE
	} /* if (t < 0.0) */
	#endif
#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */	#ifdef NODIRECTIONCHANGEDUNTILNITER /* No change in drections until some iterations are done */
} /NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done /	} /NODIRECTIONCHANGEDUNTILNITER No change in drections until some iterations are done /
#else	#else
Line 5300 void mlikeli(FILE *ficres,double p[], in	Line 5448 void mlikeli(FILE *ficres,double p[], in
for (i=1;i<=npar;i++) /* Starting with canonical directions j=1,n xi[i=1,n][j] */	for (i=1;i<=npar;i++) /* Starting with canonical directions j=1,n xi[i=1,n][j] */
for (j=1;j<=npar;j++)	for (j=1;j<=npar;j++)
xi[i][j]=(i==j ? 1.0 : 0.0);	xi[i][j]=(i==j ? 1.0 : 0.0);
printf("Powell\n"); fprintf(ficlog,"Powell\n");	printf("Powell-prax\n"); fprintf(ficlog,"Powell-prax\n");
strcpy(filerespow,"POW_");	strcpy(filerespow,"POW_");
strcat(filerespow,fileres);	strcat(filerespow,fileres);
if((ficrespow=fopen(filerespow,"w"))==NULL) {	if((ficrespow=fopen(filerespow,"w"))==NULL) {
Line 5364 void mlikeli(FILE *ficres,double p[], in	Line 5512 void mlikeli(FILE *ficres,double p[], in
}	}
powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);	powell(p,xi,npar,ftol,&iter,&fret,flatdir,func);
#else /* FLATSUP */	#else /* FLATSUP */
/* powell(p,xi,npar,ftol,&iter,&fret,func);*/	powell(p,xi,npar,ftol,&iter,&fret,func);
/* praxis ( t0, h0, n, prin, x, beale_f ); */	/* praxis ( t0, h0, n, prin, x, beale_f ); */
/* int prin=4; */	/* int prin=4; */
/* double h0=0.25; */	/* double h0=0.25; */
Line 6745 void concatwav(int wav[], int **dh, int	Line 6893 void concatwav(int wav[], int **dh, int
if(j==0) j=1; /* Survives at least one month after exam */	if(j==0) j=1; /* Survives at least one month after exam */
else if(j<0){	else if(j<0){
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 (around 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. 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 (around 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. 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;
Line 6782 void concatwav(int wav[], int **dh, int	Line 6930 void concatwav(int wav[], int **dh, int
/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){
nberr++;	nberr++;
printf("Error! Negative delay (%d) 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) between waves %d and %d of individual %ld (around 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) 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) between waves %d and %d of individual %ld (around 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]);
}	}
sum=sum+j;	sum=sum+j;
}	}
Line 8530 divided by h: <sub>h</sub>P<sub>ij</sub>	Line 8678 divided by h: <sub>h</sub>P<sub>ij</sub>
<img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);	<img src=\"%s_%d-3-%d.svg\">",stepm,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres,subdirf2(optionfilefiname,"PE_"),k1,nres);
/* Survival functions (period) in state j */	/* Survival functions (period) in state j */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);	fprintf(fichtm,"<br>\n- Survival functions in state %d. And probability to be observed in state %d being in state (1 to %d) at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, cpt, nlstate, subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));	fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);	fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJ_"),cpt,k1,nres);
}	}
/* State specific survival functions (period) */	/* State specific survival functions (period) */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\	fprintf(fichtm,"<br>\n- Survival functions in state %d and in any other live state (total).\
And probability to be observed in various states (up to %d) being in state %d at different ages. \	And probability to be observed in various states (up to %d) being in state %d at different ages. Mean times spent in state (or Life Expectancy or Health Expectancy etc.) are the areas under each curve. \
<a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);	<a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> ", cpt, nlstate, cpt, subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres,subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));	fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);	fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">",subdirf2(optionfilefiname,"LIJT_"),cpt,k1,nres);
}	}
/* Period (forward stable) prevalence in each health state */	/* Period (forward stable) prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);	fprintf(fichtm,"<br>\n- Convergence to period (stable) prevalence in state %d. Or probability for a person being in state (1 to %d) at different ages, to be alive in state %d some years after. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br>", cpt, nlstate, cpt, subdirf2(optionfilefiname,"P_"),cpt,k1,nres,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));	fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"PIJ_"),subdirf2(optionfilefiname,"PIJ_"));
fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);	fprintf(fichtm,"<img src=\"%s_%d-%d-%d.svg\">" ,subdirf2(optionfilefiname,"P_"),cpt,k1,nres);
}	}
Line 8569 divided by h: <sub>h</sub>P<sub>ij</sub>	Line 8717 divided by h: <sub>h</sub>P<sub>ij</sub>
/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */	/* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
for(cpt=1; cpt<=nlstate;cpt++){	for(cpt=1; cpt<=nlstate;cpt++){
fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \	fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \	from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d). Randomness in cross-sectional prevalence is not taken into \
account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \	account but can visually be appreciated. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);	with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a>", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_"));	fprintf(fichtm," (data from text file <a href=\"%s.txt\">%s.txt</a>)\n<br>",subdirf2(optionfilefiname,"FB_"),subdirf2(optionfilefiname,"FB_"));
fprintf(fichtm," <img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);	fprintf(fichtm," <img src=\"%s_%d-%d-%d.svg\">", subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
Line 13475 int main(int argc, char *argv[])	Line 13623 int main(int argc, char *argv[])
getcwd(pathcd, size);	getcwd(pathcd, size);
#endif	#endif
syscompilerinfo(0);	syscompilerinfo(0);
printf("\nIMaCh prax version %s, %s\n%s",version, copyright, fullversion);	printf("\nIMaCh prax version minfit %s, %s\n%s",version, copyright, fullversion);
if(argc <=1){	if(argc <=1){
printf("\nEnter the parameter file name: ");	printf("\nEnter the parameter file name: ");
if(!fgets(pathr,FILENAMELENGTH,stdin)){	if(!fgets(pathr,FILENAMELENGTH,stdin)){
Line 14468 Interval (in months) between two waves:	Line 14616 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");
#else	#else
printf("Powell\n"); fprintf(ficlog,"Powell\n");	printf("Powell-mort\n"); fprintf(ficlog,"Powell-mort\n");
#endif	#endif
strcpy(filerespow,"POW-MORT_");	strcpy(filerespow,"POW-MORT_");
strcat(filerespow,fileresu);	strcat(filerespow,fileresu);
Line 14571 Interval (in months) between two waves:	Line 14719 Interval (in months) between two waves:

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++)
matcov[i][j]=matcov[j][i];	matcov[i][j]=matcov[j][i];

printf("\nCovariance matrix\n ");	printf("\nCovariance matrix\n ");
fprintf(ficlog,"\nCovariance matrix\n ");	fprintf(ficlog,"\nCovariance matrix\n ");
Line 15025 Please run with mle=-1 to get a correct	Line 15173 Please run with mle=-1 to get a correct
}	}

/* Results */	/* Results */
/* Value of covariate in each resultine will be compututed (if product) and sorted according to model rank */	/* Value of covariate in each resultine will be computed (if product) and sorted according to model rank */
/* It is precov[] because we need the varying age in order to compute the real cov[] of the model equation */	/* It is precov[] because we need the varying age in order to compute the real cov[] of the model equation */
precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1);	precov=matrix(1,MAXRESULTLINESPONE,1,NCOVMAX+1);
endishere=0;	endishere=0;
Line 15428 Please run with mle=-1 to get a correct	Line 15576 Please run with mle=-1 to get a correct
/* */	/* */
if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */	if(i1 != 1 && TKresult[nres]!= k) /* TKresult[nres] is the combination of this nres resultline. All the i1 combinations are not output */
continue;	continue;
printf("\n# model %s \n#****** Result for:", model); /* HERE model is empty */	printf("\n# model=1+age+%s \n#****** Result for:", model); /* HERE model is empty */
fprintf(ficrest,"\n# model %s \n#****** Result for:", model);	fprintf(ficrest,"\n# model=1+age+%s \n#****** Result for:", model);
fprintf(ficlog,"\n# model %s \n#****** Result for:", model);	fprintf(ficlog,"\n# model=1+age+%s \n#****** Result for:", model);
/* It might not be a good idea to mix dummies and quantitative */	/* It might not be a good idea to mix dummies and quantitative */
/* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model \/ /	/* for(j=1;j<=cptcoveff;j++){ /\* j=resultpos. Could be a loop on cptcovs: number of single dummy covariate in the result line as well as in the model \/ /
for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */	for(j=1;j<=cptcovs;j++){ /* j=resultpos. Could be a loop on cptcovs: number of single covariate (dummy or quantitative) in the result line as well as in the model */
Line 15625 Please run with mle=-1 to get a correct	Line 15773 Please run with mle=-1 to get a correct
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);	free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
} /* mle==-3 arrives here for freeing */	} /* mle==-3 arrives here for freeing */
/* endfree:*/	/* endfree:*/
	if(mle!=-3) free_matrix(precov, 1,MAXRESULTLINESPONE,1,NCOVMAX+1); /* Could be elsewhere ?*/
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);

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

Removed from v.1.4
changed lines
	Added in v.1.5