From 1d6c98fb35d9343a3788f656edffe91f33928907 Mon Sep 17 00:00:00 2001 From: nbrouard Date: Thu, 5 Jun 2025 09:58:42 +0200 Subject: [PATCH] * src/imach.c: Trying to debug age*age but very difficult because praxis jumps anywhere and lli is corrupted. Therefore trying to test isinf or isnan which is time consuming. However it looks wrong with Zachs example. --- src/Makefile | 13 +-- src/imach.c | 236 +++++++++++++++++++++++++++++++++++++-------------- 2 files changed, 180 insertions(+), 69 deletions(-) diff --git a/src/Makefile b/src/Makefile index 95f74bf..102385a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -2,7 +2,7 @@ #VERSION=0.99s7 #VERSION=$(shell echo `grep IMACH_VERSION__ version.h | echo 'titi'`) VERSION=$(shell echo `grep IMACH_VERSION__ version.h | awk 'BEGIN { FS = "[ \t\n\"]+" } { print $$3 }'`) -VERSION=0.99s12 +#VERSION=0.99s12 OSTYPE = $(shell echo $$OSTYPE) COPYRIGHT=Copyright (C) 2002-2024 INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121) - Intel Software 2016-19 @@ -23,9 +23,9 @@ IMACHSETUPVERSION=icl # make _macosx=1 imachopt # -#CC= gcc -v +CC= gcc -v #CC=$(GCC) -CC=clang +#CC=clang GCCOPT=$(CC) CCOPT=$(GCCOPT) #GCC= gcc @@ -85,12 +85,15 @@ endif ifdef _linux -CC=clang +CC=gcc +#CC=clang #CFLAGS= -g -DUNIX -DDEBUGHESS #CFLAGS= -g -DDEBUG -DFIXMNBRAK #CFLAGS= -g -DDEBUG -CFLAGS= -g +#CFLAGS= -g -DLINMINORIGINAL -DDEBUGLINMIN -DDEBUG -DDEBUGGOMP +CFLAGS= -g -DDEBUGPRAX -DDEBUGPRAXSS LFLAGS= -g -lm +STRIP= strip INLOPT= -I/usr/local/include LNLOPT= -lm -L/usr/local/lib -lnlopt CFLAGSOPT= -O3 -g diff --git a/src/imach.c b/src/imach.c index ef8910f..b8a4a87 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1584,6 +1584,7 @@ static double maxarg1,maxarg2; #define rint(a) floor(a+0.5) /* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */ #define mytinydouble 1.0e-16 +#define myhugeexpdouble 100. /* exp(100)=2.68*e43 /* #define DEQUAL(a,b) (fabs((a)-(b)) fb and fb < fc. *fa=(*func)(*ax); /* xta[j]=pcom[j]+(*ax)*xicom[j]; fa=f(xta[j])*/ *fb=(*func)(*bx); /* xtb[j]=pcom[j]+(*bx)*xicom[j]; fb=f(xtb[j]) */ +#ifdef DEBUGLINMIN + printf(" DEBUGLINMIN mnbrak *fa=%12.7f *fb=%12.7f, *ax=%12.7f *bx=%12.7f\n", *fa, *fb, *ax, *bx); +#endif /* while(*fb != *fb){ /\* *ax should be ok, reducing distance to *ax *\/ */ /* printf("Warning mnbrak *fb = %lf, *bx=%lf *ax=%lf *fa==%lf iter=%d\n",*fb, *bx, *ax, *fa, iterscale++); */ @@ -3081,6 +3085,7 @@ static double flin(double l, int j) /* double l; */ { int i; + double f2deb; /* #ifndef MSDOS */ /* double tflin[N]; */ /* #endif */ @@ -3121,6 +3126,11 @@ static double flin(double l, int j) return (*fun)((tflin-1), n); #else /* return (*fun)(tflin, n);*/ + f2deb= (*fun)(tflin); + if(isnan(f2deb)){ + printf("Isnan flin, f2deb=%14.8e \n",f2deb); + } + printf("flin, f2deb=%14.8e \n",f2deb); return (*fun)(tflin); #endif } @@ -3157,19 +3167,44 @@ void minny(int j, int nits, double *d2, double *x1, double f1, int fk) matprint(" min vectors:",v,n,n); #endif /* find step size */ +#ifdef DEBUGPRAX + printf(" Find step size\n"); +#endif s = 0.; /* for (i=0; i s) t2 = s; - if (t2 < small_windows) t2 = small_windows; - if (t2 > 0.01*h) t2 = 0.01 * h; + if (dz && t2 > s){ + t2 = s; +#ifdef DEBUGPRAXSS + printf(" additional flin dz&&t2>s s=%14.7f t2=%14.7f \n",s,t2); +#endif + } + if (t2 < small_windows){ + t2 = small_windows; +#ifdef DEBUGPRAXSS + printf(" additional flin t2 < small_windows s=%14.7f t2=%14.7f \n",s,t2); +#endif + } + if (t2 > 0.01*h){ + t2 = 0.01 * h; + } if (fk && f1 <= fm) { +#ifdef DEBUGPRAXSS + printf(" additional flin fk && f1 <= fm X1=%14.7f t2=%14.7f f1=%14.7f fm=%14.7f fk=%d\n",*x1,t2,f1,fm,fk); +#endif xm = *x1; fm = f1; } @@ -3202,6 +3237,11 @@ L0: /*L0 loop or next */ #endif f2 = flin(x2, j); #ifdef DEBUGPRAX + if(isnan(f2)){ + printf("Isnan flin, f2=%14.8e \n",f2); + printf(" second additional second flin x2=%14.8e x1=%14.8e f0=%14.8e f1=%18.12e dirj=%d\n",x2,*x1,f0,f1,j); + f2 = flin(x2, j); + } printf(" additional second flin x2=%16.10e x1=%16.10e f1=%18.12e f0=%18.10e f2=%18.10e fm=%18.10e\n",x2, *x1, f1,f0,f2,fm); #endif if (f2 <= fm) { @@ -3213,13 +3253,13 @@ L0: /*L0 loop or next */ #ifdef DEBUGPRAX double d11,d12; d11=(f1-f0)/(*x1);d12=(f2-f0)/x2; - printf(" d11=%18.12e d12=%18.12e d11-d12=%18.12e x1-x2=%18.12e (d11-d12)/(x2-(*x1))=%18.12e\n", d11 ,d12, d11-d12, x2-(*x1), (d11-d12)/(x2-(*x1))); + printf(" d11=%18.12e d12=%18.12e d11-d12=%18.12e x1-x2=%18.12e (d11-d12)/((*x1)-x2)=%18.12e\n", d11 ,d12, d11-d12, x2-(*x1), (d11-d12)/((*x1)-x2)); printf(" original computing f1=%18.12e *d2=%16.10e f0=%18.12e f1-f0=%16.10e f2-f0=%16.10e\n",f1,*d2,f0,f1-f0, f2-f0); - double ff1=7.783920622852e+04; - double f1mf0=9.0344736236e-05; - *d2 = (f1mf0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); + /* double ff1=7.783920622852e+04; */ + /* double f1mf0=9.0344736236e-05; */ + /* *d2 = (f1mf0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); */ /* *d2 = (ff1-f0)/ (*x1)/((*x1)-x2) - (f2-f0)/x2/((*x1)-x2); */ - printf(" simpliff computing *d2=%16.10e f1mf0=%18.12e,f1=f0+f1mf0=%18.12e\n",*d2,f1mf0,f0+f1mf0); + /* printf(" simpliff computing *d2=%16.10e f1mf0=%18.12e,f1=f0+f1mf0=%18.12e\n",*d2,f1mf0,f0+f1mf0); */ *d2 = ((f1-f0)/ (*x1) - (f2-f0)/x2)/((*x1)-x2); printf(" overlifi computing *d2=%16.10e\n",*d2); #endif @@ -3252,7 +3292,7 @@ L1: /* L1 or try loop */ #endif f2 = flin(x2, j); /* x[i]+x2*v[i][j] */ #ifdef DEBUGPRAX - printf(" after flin f0=%14.8e f1=%14.8e f2=%14.8e fm=%14.8e\n",f0,f1,f2, fm); + printf(" after flin f0=%14.8e f1=%14.8e f2=%14.8e fm=%14.8e x2=%14.8e *x1=%14.8e\n",f0,f1,f2, fm, x2, *x1); #endif if ((k < nits) && (f2 > f0)) { #ifdef DEBUGPRAX @@ -4398,6 +4438,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret /* Convergence test will use last linmin estimation (fret) and compare to former iteration (fp) */ /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit */ /* New value of last point Pn is not computed, P(n-1) */ +#ifdef FLATSUP for(j=1;j<=n;j++) { if(flatdir[j] >0){ printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]); @@ -4406,6 +4447,7 @@ void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret /* printf("\n"); */ /* fprintf(ficlog,"\n"); */ } +#endif /* FLATSUP */ /* if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /\* Did we reach enough precision? *\/ */ if (2.0*fabs(fp-(*fret)) <= ftol) { /* Did we reach enough precision? */ /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */ @@ -5071,19 +5113,50 @@ double **pmij(double **ps, double *cov, int ncovmodel, double *x, int nlstate ) s1=0; for(j=1; j= myhugeexpdouble) /* should return close to HUGE_VAL */ + s1+=exp(myhugeexpdouble); + else + s1+=exp(ps[i][j]); /* In fact sums pij/pii */ +#endif } for(j=i+1; j<=nlstate+ndeath; j++){ /* printf("debug2 %d %d ps=%lf exp(ps)=%lf \n",i,j,ps[i][j],exp(ps[i][j])); */ +#ifdef INFINITYORIGINAL s1+=exp(ps[i][j]); /* In fact sums pij/pii */ +#else + if(ps[i][j] >= myhugeexpdouble) /* should return close to HUGE_VAL */ + s1+=exp(myhugeexpdouble); + else + s1+=exp(ps[i][j]); /* In fact sums pij/pii */ +#endif } /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */ - ps[i][i]=1./(s1+1.); + ps[i][i]=1./(s1+1.); /* in fact pii */ /* Computing other pijs */ - for(j=1; j= myhugeexpdouble) /* should return close to HUGE_VAL */ + ps[i][j]= exp(myhugeexpdouble)*ps[i][i];/* in fact pij, Bug valgrind */ + else + ps[i][j]= exp(ps[i][j])*ps[i][i];/* in fact pij, Bug valgrind */ +#endif + } for(j=i+1; j<=nlstate+ndeath; j++) - ps[i][j]= exp(ps[i][j])*ps[i][i]; + { +#ifdef INFINITYORIGINAL + ps[i][j]= exp(ps[i][j])*ps[i][i];/* in fact pij, Bug valgrind */ +#else + if(ps[i][j] >= myhugeexpdouble) /* should return close to HUGE_VAL */ + ps[i][j]= exp(myhugeexpdouble)*ps[i][i];/* in fact pij, Bug valgrind */ + else + ps[i][j]= exp(ps[i][j])*ps[i][i];/* in fact pij, Bug valgrind */ +#endif + } /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */ } /* end i */ @@ -5832,15 +5905,15 @@ double func( double *x) it should be very low but not zero otherwise the log go to infinity. */ -/* #ifdef INFINITYORIGINAL */ -/* lli=log(out[s1][s2] - savm[s1][s2]); */ -/* #else */ -/* if ((out[s1][s2] - savm[s1][s2]) < mytinydouble) */ -/* lli=log(mytinydouble); */ -/* else */ -/* lli=log(out[s1][s2] - savm[s1][s2]); */ -/* #endif */ +#ifdef INFINITYORIGINAL lli=log(out[s1][s2] - savm[s1][s2]); +#else + if ((out[s1][s2] - savm[s1][s2]) < mytinydouble) + lli=log(mytinydouble); + else + lli=log(out[s1][s2] - savm[s1][s2]); +#endif + /* lli=log(out[s1][s2] - savm[s1][s2]); */ } else if ( s2==-1 ) { /* alive */ for (j=1,survp=0. ; j<=nlstate; j++) @@ -5859,7 +5932,19 @@ double func( double *x) /* lli= log(survp); */ /* } */ else{ +#ifdef INFINITYORIGINAL lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */ +#else + if((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2] (double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */ } /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/ @@ -5867,6 +5952,12 @@ double func( double *x) /* printf("num[i], i=%d, bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */ ipmx +=1; sw += weight[i]; +#ifdef DEBUGPRAX + if (isinf(lli)){ + printf("isnan (lli)\n"); + printf("isnana lli=%14.8e i=%d mw[mi][i]=%d s1=%d s2=%d out[s1][s2]=%14.8e savm[s1][s2]=%14.8e\n",lli,i,mw[mi][i], s1, s2,out[s1][s2],savm[s1][s2] ); + } +#endif ll[s[mw[mi][i]][i]] += 2*weight[i]*lli; /* if (lli < log(mytinydouble)){ */ /* printf("Close to inf lli = %.10lf < %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */ @@ -5976,7 +6067,15 @@ double func( double *x) s1=s[mw[mi][i]][i]; s2=s[mw[mi+1][i]][i]; if( s2 > nlstate){ +#ifdef INFINITYORIGINAL lli=log(out[s1][s2] - savm[s1][s2]); +#else + if ((out[s1][s2] - savm[s1][s2]) < mytinydouble) + lli=log(mytinydouble); + else + lli=log(out[s1][s2] - savm[s1][s2]); +#endif + /* lli=log(out[s1][s2] - savm[s1][s2]); */ } else if ( s2==-1 ) { /* alive */ for (j=1,survp=0. ; j<=nlstate; j++) survp += out[s1][j]; @@ -6345,8 +6444,15 @@ double funcone( double *x) /* bias is positive if real duration * is higher than the multiple of stepm and negative otherwise. */ - if( s2 > nlstate && (mle <5) ){ /* Jackson */ + if( s2 > nlstate && (mle <5) ){ /* Jackson *//* if p13 is very close to 1, out - savm is close to 0 or negative */ +#ifdef INFINITYORIGINAL lli=log(out[s1][s2] - savm[s1][s2]); +#else + if ((out[s1][s2] - savm[s1][s2]) < mytinydouble) + lli=log(mytinydouble); + else + lli=log(out[s1][s2] - savm[s1][s2]); +#endif } else if ( s2==-1 ) { /* alive */ for (j=1,survp=0. ; j<=nlstate; j++) survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j]; @@ -7149,9 +7255,10 @@ void date2dmy(double date,double *day, double *month, double *year){ void freqsummary(char fileres[], double p[], double pstart[], int iagemin, int iagemax, int **s, double **agev, int nlstate, int imx, \ int *Tvaraff, int *invalidvarcomb, int **nbcode, int *ncodemax,double **mint,double **anint, char strstart[], \ int firstpass, int lastpass, int stepm, int weightopt, char model[]) -{ /* Some frequencies as well as proposing some starting values */ +{ + /* Some frequencies as well as proposing some starting values */ /* Frequencies of any combination of dummy covariate used in the model equation */ - int i, m, jk, j1, bool, z1,j, nj, nl, k, iv, jj=0, s1=1, s2=1; + int i, m, jk, j1, ibool, z1,j, nj, nl, k, iv, jj=0, s1=1, s2=1; int iind=0, iage=0; int mi; /* Effective wave */ int first; @@ -7297,7 +7404,7 @@ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age /* For that combination of covariates j1 (V4=1 V3=0 for example), we count and print the frequencies in one pass */ for (iind=1; iind<=imx; iind++) { /* For each individual iind */ - bool=1; + ibool=1; if(j !=0){ if(anyvaryingduminmodel==0){ /* If All fixed covariates */ if (cptcoveff >0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */ @@ -7312,9 +7419,9 @@ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age printf("Error Tvaraff[z1]=%d<1 or >=%d, cptcoveff=%d model=1+age+%s\n",Tvaraff[z1],NCOVMAX, cptcoveff, model); if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]){ /* for combination j1 of covariates */ /* Tests if the value of the covariate z1 for this individual iind responded to combination j1 (V4=1 V3=0) */ - bool=0; /* bool should be equal to 1 to be selected, one covariate value failed */ - /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", */ - /* bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),*/ + ibool=0; /* ibool should be equal to 1 to be selected, one covariate value failed */ + /* printf("ibool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", */ + /* ibool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),*/ /* j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/ /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/ } /* Onlyf fixed */ @@ -7322,7 +7429,7 @@ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age } /* cptcoveff > 0 */ } /* end any */ }/* end j==0 */ - if (bool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ + if (ibool==1){ /* We selected an individual iind satisfying combination j1 (V4=1 V3=0) or all fixed covariates */ /* for(m=firstpass; m<=lastpass; m++){ */ for(mi=firstpass; miDatafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age if (cotvar[m][iv][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) /* iv=1 to ntv, right modality. If covariate's value is -1, we don't select. It differs from the constant and age model which counts them. */ - bool=0; /* not selected */ + ibool=0; /* not selected */ }else if( Fixed[Tmodelind[z1]]== 0) { /* fixed */ /* i1=Tvaraff[z1]; */ /* i2=TnsdVar[i1]; */ @@ -7343,14 +7450,14 @@ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age /* i4=covar[i1][iind]; */ /* if(i4 != i3){ */ if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,TnsdVar[Tvaraff[z1]])]) { /* Bug valgrind */ - bool=0; + ibool=0; } } } }/* Some are varying covariates, we tried to speed up if all fixed covariates in the model, avoiding waves loop */ } /* end j==0 */ - /* bool =0 we keep that guy which corresponds to the combination of dummy values */ - if(bool==1){ /*Selected */ + /* ibool =0 we keep that guy which corresponds to the combination of dummy values */ + if(ibool==1){ /*Selected */ /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind] and mw[mi+1][iind]. dh depends on stepm. */ agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/ @@ -7388,15 +7495,15 @@ Title=%s
Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */ } }else{ - bool=1; - }/* end bool 2 */ + ibool=1; + }/* end ibool 2 */ } /* end m */ /* for (z1=1; z1<= nqfveff; z1++) { /\* Quantitative variables, calculating mean *\/ */ /* idq[z1]=idq[z1]+weight[iind]; */ /* meanq[z1]+=covar[ncovcol+z1][iind]*weight[iind]; /\* Computes mean of quantitative with selected filter *\/ */ /* stdq[z1]+=covar[ncovcol+z1][iind]*covar[ncovcol+z1][iind]*weight[iind]*weight[iind]; /\* *weight[iind];*\/ /\* Computes mean of quantitative with selected filter *\/ */ /* } */ - } /* end bool */ + } /* end ibool */ } /* end iind = 1 to imx */ /* prop[s][age] is fed for any initial and valid live state as well as freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */ @@ -7829,7 +7936,7 @@ void prevalence(double ***probs, double agemin, double agemax, int **s, double * We still use firstpass and lastpass as another selection. */ - int i, m, jk, j1, bool, z1,j, iv; + int i, m, jk, j1, ibool, z1,j, iv; int mi; /* Effective wave */ int iage; double agebegin; /*, ageend;*/ @@ -7860,7 +7967,7 @@ void prevalence(double ***probs, double agemin, double agemax, int **s, double * fprintf(ficlog,"Prevalence combination of varying and fixed dummies %d\n",j1); int ierroragemin=0;int ierroragemax=0; for (i=1; i<=imx; i++) { /* Each individual */ - bool=1; + ibool=1; /* for(m=firstpass; m<=lastpass; m++){/\* Other selection (we can limit to certain interviews*\/ */ for(mi=firstpass; mi=firstpass && m <=lastpass){ @@ -7899,7 +8006,7 @@ void prevalence(double ***probs, double agemin, double agemax, int **s, double * } /* end valid statuses */ } /* end selection of dates */ } /* end selection of waves */ - } /* end bool */ + } /* end ibool */ } /* end wave */ } /* end individual */ if(ierroragemin>=1 || ierroragemax>=1){ @@ -12430,11 +12537,10 @@ double gompertz(double x[]) } L=0.0; /* agegomp=AGEGOMP; */ - /* for (i=0; i<=imx; i++) - if (wav[i]>0) printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d %d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]);*/ - + /* for (i=1; i<=imx; i++) */ + /* printf("i=%d ageex=%lf agecens=%lf agedc=%lf cens=%d wav=%d\n" ,i,ageexmed[i],agecens[i],agedc[i],cens[i],wav[i]); */ for (i=1;i<=imx ; i++) { - /* mu(a)=mu(agecomp)*exp(teta*(age-agegomp)) + /* mu(a)=mu(agegomp)*exp(teta*(age-agegomp)) mu(a)=x[1]*exp(x[2]*(age-agegomp)); x[1] and x[2] are per year. * L= Product mu(agedeces)exp(-\int_ageexam^agedc mu(u) du ) for a death between agedc (in month) * and agedc +1 month, cens[i]=0: log(x[1]/YEARM) @@ -12445,19 +12551,21 @@ double gompertz(double x[]) if (cens[i] == 1){ A=-x[1]/(x[2])*(exp(x[2]*(agecens[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))); } else if (cens[i] == 0){ - A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) - +log(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); + A=-x[1]/(x[2])*(exp(x[2]*(agedc[i]-agegomp))-exp(x[2]*(ageexmed[i]-agegomp))) + +log(fabs(x[1])/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); /* +log(x[1]/YEARM) +x[2]*(agedc[i]-agegomp)+log(YEARM); */ /* To be seen */ } else printf("Gompertz cens[%d] neither 1 nor 0\n",i); /*if (wav[i] > 1 && agecens[i] > 15) {*/ /* ??? */ L=L+A*weight[i]; - /* printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]);*/ +#ifdef DEBUGGOMP + printf("\ni=%d A=%f L=%lf x[1]=%lf x[2]=%lf ageex=%lf agecens=%lf cens=%d agedc=%lf weight=%lf\n",i,A,L,x[1],x[2],ageexmed[i]*12,agecens[i]*12,cens[i],agedc[i]*12,weight[i]); +#endif } } - - /*printf("x1=%2.9f x2=%2.9f x3=%2.9f L=%f\n",x[1],x[2],x[3],L);*/ - +#ifdef DEBUGGOMP + printf("x1=%2.9lf x2=%2.9lf x3=%2.9lf L=%lf\n",x[1],x[2],x[3],L); +#endif return -2*L*num/sump; } @@ -16037,16 +16145,16 @@ Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf
\n",\ flatdir=ivector(1,npar); for (j=1;j<=npar;j++) flatdir[j]=0; #endif /*LINMINORIGINAL */ - /* powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); */ - /* double h0=0.25; */ - macheps=pow(16.0,-13.0); - printf("Praxis Gegenfurtner mle=%d\n",mle); - fprintf(ficlog, "Praxis Gegenfurtner mle=%d\n", mle);fflush(ficlog); - /* ffmin = praxis(ftol,macheps, h0, npar, prin, p, gompertz); */ - /* For the Gompertz we use only two parameters */ - int _npar=2; - ffmin = praxis(ftol,macheps, h0, _npar, 4, p, gompertz); - printf("End Praxis\n"); + powell(p,ximort,NDIM,ftol,&iter,&fret,gompertz); + /* /\* double h0=0.25; *\/ */ + /* macheps=pow(16.0,-13.0); */ + /* printf("Praxis Gegenfurtner mle=%d\n",mle); */ + /* fprintf(ficlog, "Praxis Gegenfurtner mle=%d\n", mle);fflush(ficlog); */ + /* /\* ffmin = praxis(ftol,macheps, h0, npar, prin, p, gompertz); *\/ */ + /* /\* For the Gompertz we use only two parameters *\/ */ + /* int _npar=2; */ + /* ffmin = praxis(ftol,macheps, h0, _npar, 4, p, gompertz); */ + /* printf("End Praxis\n"); */ fclose(ficrespow); #ifdef LINMINORIGINAL #else -- 2.43.0