The same imach parameter file can be used but the option for mle should be -3.
- Agnès, who wrote this part of the code, tried to keep most of the
+ Agnès, who wrote this part of the code, tried to keep most of the
former routines in order to include the new code within the former code.
The output is very simple: only an estimate of the intercept and of
- Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
- Institut national d'études démographiques, Paris.
+ Authors: Nicolas Brouard (brouard@ined.fr) and Agnès Lièvre (lievre@ined.fr).
+ Institut national d'études démographiques, Paris.
This software have been partly granted by Euro-REVES, a concerted action
from the European Union.
It is copyrighted identically to a GNU software product, ie programme and
#include <io.h>
#include <windows.h>
#include <tchar.h>
+#include <direct.h>
#else
#include <unistd.h>
#endif
#include <limits.h>
#include <sys/types.h>
-#if defined(__GNUC__)
+#if defined(__GNUC__) && !defined(_WIN32)
#include <sys/utsname.h> /* Doesn't work on Windows */
#endif
/* From a file name with (full) path (either Unix or Windows) we extract the directory (dirc)
the name of the file (name), its extension only (ext) and its first part of the name (finame)
*/
- char *ss; /* pointer */
- int l1=0, l2=0; /* length counters */
+ char *ss, *sso; /* pointer */
+ int l1=0, l2=0, lss=0, lsso=0; /* length counters */
l1 = strlen(path ); /* length of path */
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
- ss= strrchr( path, DIRSEPARATOR ); /* find last / */
- if ( ss == NULL ) { /* no directory, so determine current directory */
+ if((ss= strrchr( path, DIRSEPARATOR )) != NULL){/* find last of other / */
+ lss=strlen(ss);
+ printf(" strlen lss=%d path=%s\n", lss, path);
+ }
+ while((sso= strrchr( path, ODIRSEPARATOR )) != NULL){ /* find last of other / */
+ lsso=strlen(sso);
+ printf(" strlen lsso=%d OLDpath=%s\n", lsso, path);
+ path[l1-lsso]=DIRSEPARATOR;
+ printf(" strlen lsso=%d NEWpath=%s\n", lsso, path);
+ ss=sso;
+ printf(" NEWss=%s\n", ss);
+ }
+ fflush(stdout);
+ if ( ss == NULL ){ /* no directory, so determine current directory */
strcpy( name, path ); /* we got the fullname name because no directory */
/*if(strrchr(path, ODIRSEPARATOR )==NULL)
printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
/* get current working directory */
/* extern char* getcwd ( char *buf , int len);*/
#ifdef WIN32
- if (_getcwd( dirc, FILENAME_MAX ) == NULL ) {
+ if (_getcwd( dirc, FILENAME_MAX ) == NULL )
#else
- if (getcwd(dirc, FILENAME_MAX) == NULL) {
+ if (getcwd(dirc, FILENAME_MAX) == NULL)
#endif
return( GLOCK_ERROR_GETCWD );
- }
+
/* got dirc from getcwd*/
printf(" DIRC = %s \n",dirc);
} else { /* strip directory from path */
{
/* Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is
* between ax and cx, and f(bx) is less than both f(ax) and f(cx) ), this routine isolates
- * the minimum to a fractional precision of about tol using Brent’s method. The abscissa of
+ * the minimum to a fractional precision of about tol using Brent’s method. The abscissa of
* the minimum is returned as xmin, and the minimum function value is returned as brent , the
* returned function value.
*/
if (prin) print2();
mloop:
- biter++; /* Added to count the loops */
+ biter++; /* Added to count the loops */
/* sf = d[0]; */
/* s = d[0] = 0.0; */
- printf("\n Big iteration %d \n",biter);
- fprintf(ficlog,"\n Big iteration %d \n",biter);
- sf = d[1];
+ printf("\n Big iteration %d \n",biter);
+ fprintf(ficlog,"\n Big iteration %d \n",biter);
+ sf = d[1];
s = d[1] = 0.0;
/* minimize along first direction V(*,1) */
-#ifdef DEBUGPRAX
+#ifdef DEBUGPRAXS
printf(" Minimize along the first direction V(*,1). illc=%d\n",illc);
/* fprintf(ficlog," Minimize along the first direction V(*,1).\n"); */
#endif
The inner loop starts here.
*/
#ifdef DEBUGPRAX
- printf(" The inner loop here from k=%d to n=%d.\n",k,n);
- /* fprintf(ficlog," The inner loop here from k=%d to n=%d.\n",k,n); */
+ printf(" The inner loop here from k=%d to n=%d.\n",k,n);
+ /* fprintf(ficlog," The inner loop here from k=%d to n=%d.\n",k,n); */
#endif
- /* for (i=0; i<n; i++) */
- for (i=1; i<=n; i++)
- y[i] = x[i];
- sf = fx;
+ /* for (i=0; i<n; i++) */
+ for (i=1; i<=n; i++)
+ y[i] = x[i];
+ sf = fx;
#ifdef DEBUGPRAX
- printf(" illc=%d and kt=%d and ktm=%d\n", illc, kt, ktm);
+ printf(" illc=%d and kt=%d and ktm=%d\n", illc, kt, ktm);
+#endif
+ illc = illc || (kt > 0);
+ nextL1:
+ kl = k;
+ df = 0.0;
+ if (illc) { /* random step to get off resolution valley */
+#ifdef DEBUGPRAXS
+ printf(" A random step follows, to avoid resolution valleys.\n");
#endif
- illc = illc || (kt > 0);
-next:
- kl = k;
- df = 0.0;
- if (illc) { /* random step to get off resolution valley */
#ifdef DEBUGPRAX
- printf(" A random step follows, to avoid resolution valleys.\n");
- matprint(" before rand, vectors:",v,n,n);
+ matprint(" before rand, vectors:",v,n,n);
#endif
- for (i=1; i<=n; i++) {
+ for (i=1; i<=n; i++) {
#ifdef NOBRENTRAND
- r = drandom();
+ r = drandom();
#else
- seed=i;
- /* seed=i+1; */
+ seed=i;
+ /* seed=i+1; */
#ifdef DEBUGRAND
- printf(" Random seed=%d, brent i=%d",seed,i); /* YYYY i=5 j=1 vji= -0.0001170073 */
+ printf(" Random seed=%d, brent i=%d",seed,i); /* YYYY i=5 j=1 vji= -0.0001170073 */
#endif
- r = randbrent ( &seed );
+ r = randbrent ( &seed );
#endif
#ifdef DEBUGRAND
- printf(" Random r=%.7g \n",r);
+ printf(" Random r=%.7g \n",r);
#endif
- z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (r - 0.5);
- /* z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (drandom() - 0.5); */
-
- s = z[i];
- for (j=1; j <= n; j++)
- x[j] += s * v[j][i];
- }
+ z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (r - 0.5);
+ /* z[i] = (0.1 * ldt + t2 * pow(10.0,(double)kt)) * (drandom() - 0.5); */
+
+ s = z[i];
+ for (j=1; j <= n; j++)
+ x[j] += s * v[j][i];
+ }
#ifdef DEBUGRAND
- matprint(" after rand, vectors:",v,n,n);
+ matprint(" after rand, vectors:",v,n,n);
#endif
#ifdef NR_SHIFT
- fx = (*fun)((x-1), n);
+ fx = (*fun)((x-1), n);
#else
- fx = (*fun)(x);
+ fx = (*fun)(x);
#endif
- /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/
- nf++;
- }
- /* minimize along non-conjugate directions */
-#ifdef DEBUGPRAX
- printf(" Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n);
- /* fprintf(ficlog," Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n); */
+ /* fx = (*func) ( (x-1) ); *//* This for func which is computed from x[1] and not from x[0] xm1=(x-1)*/
+ nf++;
+ }
+ /* minimize along non-conjugate directions */
+#ifdef DEBUGPRAXS
+ printf(" Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n);
+ /* fprintf(ficlog," Minimize along the 'non-conjugate' directions (dots printed) V(*,%d),...,V(*,%d).\n",k,n); */
#endif
- /* for (k2=k; k2<n; k2++) { /\* Be careful here k2 <=n ? *\/ */
- for (k2=k; k2<=n; k2++) { /* Be careful here k2 <=n ? */
- sl = fx;
- s = 0.0;
-#ifdef DEBUGPRAX
- printf(" Minimize along the 'NON-CONJUGATE' true direction k2=%14d fx=%14.7f\n",k2, fx);
- matprint(" before min vectors:",v,n,n);
+ /* for (k2=k; k2<n; k2++) { /\* Be careful here k2 <=n ? *\/ */
+ for (k2=k; k2<=n; k2++) { /* Be careful here k2 <=n ? */
+ sl = fx;
+ s = 0.0;
+#ifdef DEBUGPRAXS
+ printf(" Minimize along the 'NON-CONJUGATE' true direction k2=%14d fx=%20.15f\n",k2, fx);
#endif
- /* min(k2, 2, &d[k2], &s, fx, 0); */
- /* jsearch=k2-1; */
- /* min(jsearch, 2, &d[jsearch], &s, fx, 0); */
- minny(k2, 2, &d[k2], &s, fx, 0);
#ifdef DEBUGPRAX
- printf(" . D(%d)=%14.7f d[k2]=%14.7f z[k2]=%14.7f illc=%14d fx=%14.7f\n",k2,d[k2],d[k2],z[k2],illc,fx);
+ matprint(" before min vectors:",v,n,n);
#endif
- if (illc) {
- /* double szk = s + z[k2]; */
- /* s = d[k2] * szk*szk; */
- double szk = s + z[k2];
- s = d[k2] * szk*szk;
- }
- else
- s = sl - fx;
- /* if (df < s) { */
- if (df <= s) {
- df = s;
- kl = k2;
+ /* min(k2, 2, &d[k2], &s, fx, 0); */
+ /* jsearch=k2-1; */
+ /* min(jsearch, 2, &d[jsearch], &s, fx, 0); */
+ minny(k2, 2, &d[k2], &s, fx, 0);
#ifdef DEBUGPRAX
- printf(" df=%.7g and choose kl=%d \n",df,kl); /* UUUU */
+ printf(" . D(%d)=%14.7f d[k2]=%14.7f z[k2]=%14.7f illc=%14d fx=%14.7f\n",k2,d[k2],d[k2],z[k2],illc,fx);
#endif
- }
- } /* end loop k2 */
- /*
- If there was not much improvement on the first try, set
- ILLC = true and start the inner loop again.
- */
+ if (illc) {
+ /* double szk = s + z[k2]; */
+ /* s = d[k2] * szk*szk; */
+ double szk = s + z[k2];
+ s = d[k2] * szk*szk;
+ }
+ else
+ s = sl - fx;
+ /* if (df < s) { */
+ if (df <= s) {
+ df = s;
+ kl = k2;
#ifdef DEBUGPRAX
- printf(" If there was not much improvement on the first try, set ILLC = true and start the inner loop again. illc=%d\n",illc);
- /* fprintf(ficlog," If there was not much improvement on the first try, set ILLC = true and start the inner loop again.\n"); */
+ printf(" df=%.7g and choose kl=%d \n",df,kl); /* UUUU */
#endif
- if (!illc && (df < fabs(100.0 * (macheps) * fx))) {
-#ifdef DEBUGPRAX
- printf("\n NO SUCCESS because DF is small, starts inner loop with same K(=%d), fabs( 100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e > df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc);
+ }
+ } /* end loop k2 */
+ /*
+ If there was not much improvement on the first try, set
+ ILLC = true and start the inner loop again.
+ */
+#ifdef DEBUGPRAXS
+ printf(" If there was not much improvement on the first try, set ILLC = true and start the inner loop again. illc=%d\n",illc);
+ /* fprintf(ficlog," If there was not much improvement on the first try, set ILLC = true and start the inner loop again.\n"); */
#endif
- illc = 1;
- goto next;
- }
-#ifdef DEBUGPRAX
- printf("\n SUCCESS, BREAKS inner loop K(=%d) because DF is big, fabs( 100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e <= df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc);
+ if (!illc && (df < fabs(100.0 * (macheps) * fx))) {
+#ifdef DEBUGPRAXS
+ printf("\n NO SUCCESS IILC = FALSE SO TRY ONCE WITH ILLC = TRUE because DF is small, starts inner loop with same K(=%d), fabs( 100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e > df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc);
#endif
-
- /* if ((k == 1) && (prin > 1)){ /\* be careful k=2 *\/ */
- if ((k == 2) && (prin > 1)){ /* be careful k=2 */
-#ifdef DEBUGPRAX
- printf(" NEW D The second difference array d:\n" );
- /* fprintf(ficlog, " NEW D The second difference array d:\n" ); */
+ illc = 1;
+ goto nextL1;
+ }
+#ifdef DEBUGPRAXS
+ printf("\n SUCCESS, BREAKS inner loop K(=%d) because DF is big, fabs( 100.0 * machep(=%.10e) * fx(=%.9e) )=%.9e <= df(=%.9e) break illc=%d\n", k, macheps, fx, fabs ( 100.0 * macheps * fx ), df, illc);
#endif
- vecprint(" NEW D The second difference array d:",d,n);
- }
- /* minimize along conjugate directions */
- /*
- Minimize along the "conjugate" directions V(*,1),...,V(*,K-1).
- */
-#ifdef DEBUGPRAX
- printf("Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1);
- /* fprintf(ficlog,"Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1); */
+
+ /* if ((k == 1) && (prin > 1)){ /\* be careful k=2 *\/ */
+ if ((k == 2) && (prin > 1)){ /* be careful k=2 */
+#ifdef DEBUGPRAXS
+ printf(" NEW D The second difference array d:\n" );
+ /* fprintf(ficlog, " NEW D The second difference array d:\n" ); */
#endif
- /* for (k2=0; k2<=k-1; k2++) { */
- for (k2=1; k2<=k-1; k2++) {
- s = 0.0;
- /* min(k2-1, 2, &d[k2-1], &s, fx, 0); */
- minny(k2, 2, &d[k2], &s, fx, 0);
- }
- f1 = fx;
- fx = sf;
- lds = 0.0;
- /* for (i=0; i<n; i++) { */
- for (i=1; i<=n; i++) {
- sl = x[i];
- x[i] = y[i];
- y[i] = sl - y[i];
- sl = y[i];
- lds = lds + sl*sl;
- }
- lds = sqrt(lds);
-#ifdef DEBUGPRAX
- printf("Minimization done 'conjugate', shifted all points, computed lds=%.8f\n",lds);
+ vecprint(" NEW D The second difference array d:",d,n);
+ }
+ /* minimize along conjugate directions */
+ /*
+ Minimize along the "conjugate" directions V(*,1),...,V(*,K-1).
+ */
+#ifdef DEBUGPRAXS
+ printf("Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1);
+ /* fprintf(ficlog,"Minimize along the 'conjugate' directions V(*,1),...,V(*,K-1=%d).\n",k-1); */
+#endif
+ /* for (k2=0; k2<=k-1; k2++) { */
+ for (k2=1; k2<=k-1; k2++) {
+ s = 0.0;
+ /* min(k2-1, 2, &d[k2-1], &s, fx, 0); */
+ minny(k2, 2, &d[k2], &s, fx, 0);
+ }
+ f1 = fx;
+ fx = sf;
+ lds = 0.0;
+ /* for (i=0; i<n; i++) { */
+ for (i=1; i<=n; i++) {
+ sl = x[i];
+ x[i] = y[i];
+ y[i] = sl - y[i];
+ sl = y[i];
+ lds = lds + sl*sl;
+ }
+ lds = sqrt(lds);
+#ifdef DEBUGPRAXS
+ printf("Minimization done 'conjugate', shifted all points, computed lds=%.8f\n",lds);
#endif
- /*
- Discard direction V(*,kl).
-
- If no random step was taken, V(*,KL) is the "non-conjugate"
- direction along which the greatest improvement was made.
- */
- if (lds > small_windows) {
-#ifdef DEBUGPRAX
+ /*
+ Discard direction V(*,kl).
+
+ If no random step was taken, V(*,KL) is the "non-conjugate"
+ direction along which the greatest improvement was made.
+ */
+ if (lds > small_windows) {
+#ifdef DEBUGPRAXS
printf("lds big enough to throw direction V(*,kl=%d). If no random step was taken, V(*,KL) is the 'non-conjugate' direction along which the greatest improvement was made.\n",kl);
- matprint(" before shift new conjugate vectors:",v,n,n);
#endif
- for (i=kl-1; i>=k; i--) {
- /* for (j=0; j < n; j++) */
- for (j=1; j <= n; j++)
- /* v[j][i+1] = v[j][i]; */ /* This is v[j][i+1]=v[j][i] i=kl-1 to k */
- v[j][i+1] = v[j][i]; /* This is v[j][i+1]=v[j][i] i=kl-1 to k */
- /* v[j][i+1] = v[j][i]; */
- /* d[i+1] = d[i];*/ /* last is d[k+1]= d[k] */
- d[i+1] = d[i]; /* last is d[k]= d[k-1] */
- }
#ifdef DEBUGPRAX
- matprint(" after shift new conjugate vectors:",v,n,n);
-#endif /* d[k] = 0.0; */
- d[k] = 0.0;
- for (i=1; i <= n; i++)
- v[i][k] = y[i] / lds;
- /* v[i][k] = y[i] / lds; */
+ matprint(" before shift new conjugate vectors:",v,n,n);
+#endif
+ for (i=kl-1; i>=k; i--) {
+ /* for (j=0; j < n; j++) */
+ for (j=1; j <= n; j++)
+ /* v[j][i+1] = v[j][i]; */ /* This is v[j][i+1]=v[j][i] i=kl-1 to k */
+ v[j][i+1] = v[j][i]; /* This is v[j][i+1]=v[j][i] i=kl-1 to k */
+ /* v[j][i+1] = v[j][i]; */
+ /* d[i+1] = d[i];*/ /* last is d[k+1]= d[k] */
+ d[i+1] = d[i]; /* last is d[k]= d[k-1] */
+ }
#ifdef DEBUGPRAX
- printf("Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector: (new x) - (old x). d2=%14.7g lds=%.10f\n",k,d[k],lds);
- /* fprintf(ficlog,"Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector: (new x) - (old x).\n",k); */
- matprint(" before min new conjugate vectors:",v,n,n);
+ matprint(" after shift new conjugate vectors:",v,n,n);
+#endif /* d[k] = 0.0; */
+ d[k] = 0.0;
+ for (i=1; i <= n; i++)
+ v[i][k] = y[i] / lds;
+ /* v[i][k] = y[i] / lds; */
+#ifdef DEBUGPRAXS
+ printf("Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector: (new x) - (old x). d2=%14.7g lds=%.10f\n",k,d[k],lds);
+ /* fprintf(ficlog,"Minimize along the new 'conjugate' direction V(*,k=%d), which is the normalized vector: (new x) - (old x).\n",k); */
#endif
- /* min(k-1, 4, &d[k-1], &lds, f1, 1); */
- minny(k, 4, &d[k], &lds, f1, 1);
#ifdef DEBUGPRAX
- printf(" after min d(k)=%d %.7g lds=%14f\n",k,d[k],lds);
- matprint(" after min vectors:",v,n,n);
+ matprint(" before min new conjugate vectors:",v,n,n);
+#endif
+ /* min(k-1, 4, &d[k-1], &lds, f1, 1); */
+ minny(k, 4, &d[k], &lds, f1, 1);
+#ifdef DEBUGPRAXS
+ printf(" after min d(k)=%d %.7g lds=%14f\n",k,d[k],lds);
#endif
- if (lds <= 0.0) {
- lds = -lds;
#ifdef DEBUGPRAX
- printf(" lds changed sign lds=%.14f k=%d\n",lds,k);
+ matprint(" after min vectors:",v,n,n);
+#endif
+ if (lds <= 0.0) {
+ lds = -lds;
+#ifdef DEBUGPRAXS
+ printf(" lds changed sign lds=%.14f k=%d\n",lds,k);
#endif
- /* for (i=0; i<n; i++) */
- /* v[i][k] = -v[i][k]; */
- for (i=1; i<=n; i++)
- v[i][k] = -v[i][k];
- }
+ /* for (i=0; i<n; i++) */
+ /* v[i][k] = -v[i][k]; */
+ for (i=1; i<=n; i++)
+ v[i][k] = -v[i][k];
}
- ldt = ldfac * ldt;
- if (ldt < lds)
- ldt = lds;
- if (prin > 0){
-#ifdef DEBUGPRAX
- printf(" k=%d",k);
- /* fprintf(ficlog," k=%d",k); */
+ } /* End if lds > small_windows */
+#ifdef DEBUGPRAXS
+ printf(" lds=%20.15f <= small_windows=%20.15f, ldt=%20.15f\n",lds, small_windows,ldt);
#endif
- print2();/* n, x, prin, fx, nf, nl ); */
- }
- t2 = 0.0;
- /* for (i=0; i<n; i++) */
- for (i=1; i<=n; i++)
- t2 += x[i]*x[i];
- t2 = m2 * sqrt(t2) + t;
- /*
- See whether the length of the step taken since starting the
- inner loop exceeds half the tolerance.
- */
-#ifdef DEBUGPRAX
- printf("See if step length exceeds half the tolerance.\n"); /* ZZZZZ */
- /* fprintf(ficlog,"See if step length exceeds half the tolerance.\n"); */
+ ldt = ldfac * ldt;
+#ifdef DEBUGPRAXS
+ printf(" lds=%20.15f, New ldt(=lds*ldfac)=%20.15f, ldfac=%20.15f\n",lds, ldt, ldfac);
#endif
- if (ldt > (0.5 * t2))
- kt = 0;
- else
- kt++;
-#ifdef DEBUGPRAX
- printf("if kt=%d >? ktm=%d gotoL2 loop\n",kt,ktm);
+ if (ldt < lds)
+ ldt = lds;
+ if (prin > 0){
+#ifdef DEBUGPRAXS
+ printf(" k=%d and print2 after, ldt may have changed if lower than lds=%20.15f\n",k,ldt);
+ /* fprintf(ficlog," k=%d",k); */
#endif
- if (kt > ktm){
- if ( 0 < prin ){
- /* printf("\nr8vec_print\n X:\n"); */
- /* fprintf(ficlog,"\nr8vec_print\n X:\n"); */
- vecprint ("END X:", x, n );
- }
- goto fret;
+ print2();/* n, x, prin, fx, nf, nl ); */
+#ifdef DEBUGPRAXS
+ printf(" k=%d After print2\n",k);
+#endif
+ }
+ t2 = 0.0;
+ /* for (i=0; i<n; i++) */
+ for (i=1; i<=n; i++)
+ t2 += x[i]*x[i];
+ t2 = m2 * sqrt(t2) + t;
+ /*
+ See whether the length of the step taken since starting the
+ inner loop exceeds half the tolerance.
+ */
+#ifdef DEBUGPRAXS
+ printf("See if step length ldt=%10e exceeds half the tolerance t2=%10e.\n",ldt,t2); /* ZZZZZ */
+ /* fprintf(ficlog,"See if step length exceeds half the tolerance.\n"); */
+#endif
+ if (ldt > (0.5 * t2)){
+ kt = 0;
+#ifdef DEBUGPRAXS
+ printf(" therefore kt is back to zero and back to L2, kt=%d ktm=%d\n",kt,ktm);
+#endif
+ }else{
+ kt++;
+#ifdef DEBUGPRAXS
+ printf(" therefore kt has been increased =%d ktm=%d\n",kt,ktm);
+#endif
+ }
+#ifdef DEBUGPRAXS
+ printf("Test if kt=%d >? ktm=%d gotoL2 loop else END X\n",kt,ktm);
+#endif
+ if (kt > ktm){
+ if ( 0 < prin ){
+ /* printf("\nr8vec_print\n X:\n"); */
+ /* fprintf(ficlog,"\nr8vec_print\n X:\n"); */
+ vecprint ("END X:", x, n );
}
+ goto L2;
+ }
+#ifdef DEBUGPRAXS
+ printf(" L2 loop vectors: k=%d and n=%d\n",k, n);
+#endif
#ifdef DEBUGPRAX
- matprint(" end of L2 loop vectors:",v,n,n);
+ matprint(" L2 loop vectors:",v,n,n);
#endif
-
+
}
- /* printf("The inner loop ends here.\n"); */
- /* fprintf(ficlog,"The inner loop ends here.\n"); */
+#ifdef DEBUGPRAXS
+ printf("The inner loop ends here.\n");
+ fprintf(ficlog,"The inner loop ends here.\n");
+#endif
/*
The inner loop ends here.
Try quadratic extrapolation in case we are in a curved valley.
*/
-#ifdef DEBUGPRAX
+#ifdef DEBUGPRAXS
printf("Try QUAD ratic extrapolation in case we are in a curved valley.\n");
#endif
/* try quadratic extrapolation in case */
s = vlarge;
/* for (i=0; i<n; i++) { */
for (i=1; i<=n; i++) {
- sl = 0.0;
- /* for (j=0; j < n; j++) */
- for (j=1; j <= n; j++)
- sl += v[i][j]*v[i][j];
- z[i] = sqrt(sl);
- if (z[i] < m4)
- z[i] = m4;
- if (s > z[i])
- s = z[i];
+ sl = 0.0;
+ /* for (j=0; j < n; j++) */
+ for (j=1; j <= n; j++)
+ sl += v[i][j]*v[i][j];
+ z[i] = sqrt(sl);
+ if (z[i] < m4)
+ z[i] = m4;
+ if (s > z[i])
+ s = z[i];
}
/* for (i=0; i<n; i++) { */
for (i=1; i<=n; i++) {
- sl = s / z[i];
- z[i] = 1.0 / sl;
- if (z[i] > scbd) {
- sl = 1.0 / scbd;
- z[i] = scbd;
- }
+ sl = s / z[i];
+ z[i] = 1.0 / sl;
+ if (z[i] > scbd) {
+ sl = 1.0 / scbd;
+ z[i] = scbd;
+ }
}
}
for (i=1; i<=n; i++)
- /* for (j=0; j<=i-1; j++) { */
- /* for (j=1; j<=i; j++) { */
- for (j=1; j<=i-1; j++) {
- s = v[i][j];
- v[i][j] = v[j][i];
- v[j][i] = s;
- }
+ /* for (j=0; j<=i-1; j++) { */
+ /* for (j=1; j<=i; j++) { */
+ for (j=1; j<=i-1; j++) {
+ s = v[i][j];
+ v[i][j] = v[j][i];
+ v[j][i] = s;
+ }
#ifdef DEBUGPRAX
- printf(" Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n");
+ printf(" Calculate a new set of orthogonal directions before repeating the main loop.\n Transpose V for MINFIT:...\n");
#endif
- /*
- 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.
- */
- #ifdef DEBUGPRAX
- printf(" MINFIT 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");
+ /*
+ 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.
+ */
+#ifdef DEBUGPRAXS
+ printf(" MINFIT 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");
#endif
minfit(n, macheps, vsmall, v, d);
Unscale the axes.
*/
if (scbd > 1.0) {
-#ifdef DEBUGPRAX
- printf(" Unscale the axes.\n");
+#ifdef DEBUGPRAXS
+ printf(" Unscale the axes.\n");
#endif
/* for (i=0; i<n; i++) { */
for (i=1; i<=n; i++) {
- s = z[i];
- /* for (j=0; j<n; j++) */
- for (j=1; j<=n; j++)
- v[i][j] *= s;
+ s = z[i];
+ /* for (j=0; j<n; j++) */
+ for (j=1; j<=n; j++)
+ v[i][j] *= s;
}
/* for (i=0; i<n; i++) { */
for (i=1; i<=n; i++) {
- s = 0.0;
- /* for (j=0; j<n; j++) */
- for (j=1; j<=n; j++)
- s += v[j][i]*v[j][i];
- s = sqrt(s);
- d[i] *= s;
- s = 1.0 / s;
- /* for (j=0; j<n; j++) */
- for (j=1; j<=n; j++)
- v[j][i] *= s;
+ s = 0.0;
+ /* for (j=0; j<n; j++) */
+ for (j=1; j<=n; j++)
+ s += v[j][i]*v[j][i];
+ s = sqrt(s);
+ d[i] *= s;
+ s = 1.0 / s;
+ /* for (j=0; j<n; j++) */
+ for (j=1; j<=n; j++)
+ v[j][i] *= s;
}
}
/* for (i=0; i<n; i++) { */
/* illc = (m2 * d[0]) > dmin; */
illc = (m2 * d[1]) > dmin;
#ifdef DEBUGPRAX
- printf(" The ratio of the smallest to largest eigenvalue determines whether\n the system is ill conditioned=%d . dmin=%.10lf < m2=%.10lf * d[1]=%.10lf \n",illc, dmin,m2, d[1]);
+ printf(" The ratio of the smallest to largest eigenvalue determines whether\n the system is ill conditioned=%d . dmin=%.10lf < m2=%.10lf * d[1]=%.10lf \n",illc, dmin,m2, d[1]);
#endif
if ((prin > 2) && (scbd > 1.0))
- vecprint("\n The scale factors:",z,n);
+ vecprint("\n The scale factors:",z,n);
if (prin > 2)
- vecprint(" Principal values (EIGEN VALUES OF A) of the quadratic form:",d,n);
+ vecprint(" Principal values (EIGEN VALUES OF A) of the quadratic form:",d,n);
if (prin > 2)
matprint(" The principal axes (EIGEN VECTORS OF A:",v,n, n);
-
+
if ((maxfun > 0) && (nf > maxfun)) {
- if (prin)
- printf("\n... maximum number of function calls reached ...\n");
- goto fret;
+ if (prin)
+ printf("\n... maximum number of function calls reached, exit ...\n");
+ goto L2;
}
#ifdef DEBUGPRAX
printf("Goto main loop\n");
#endif
goto mloop; /* back to main loop */
-fret:
+L2:
if (prin > 0) {
vecprint("\n X:", x, n);
/* printf("\n... ChiSq reduced to %20.10e ...\n", fx); */
/* 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 */
/* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del or directest <0 */
- /* also lamda^2=(f1-f2)^2/mu² is a parasite solution of powell */
+ /* also lamda^2=(f1-f2)^2/mu² is a parasite solution of powell */
/* For powell, inclusion of this average direction is only if t(del)<0 or del inbetween mu^2 and lambda^2 */
/* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
/* Even if f3 <f1, directest can be negative and t >0 */
- /* mu² and del² are equal when f3=f1 */
- /* 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 : lambda² < mu^2 < del then t is negative and directest >0 */
- /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */
+ /* mu² and del² are equal when f3=f1 */
+ /* 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 : lambda² < mu^2 < del then t is negative and directest >0 */
+ /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0 */
#ifdef NRCORIGINAL
t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/
#else
* FixedV[itv] 1 1 1 0 1 0 1 0 1 0 0
*? FixedV[itv] 1 1 1 0 1 0 1 0 1 0 1 0 1 0
* Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4
- * Tvar[TvarFind[itv]] [0]=? ?ncovv 1 Ã ncovvt]
+ * Tvar[TvarFind[itv]] [0]=? ?ncovv 1 Ã\83 ncovvt]
* Not a fixed cotvar[mw][itv][i] 6 7 6 2 7, 2, 6, 3, 7, 3, 6, 4, 7, 4}
* fixed covar[itv] [6] [7] [6][2]
*/
fprintf(ficlog, "Praxis Gegenfurtner\n");fflush(ficlog);
/* praxis ( ftol, h0, npar, prin, p1, func ); */
/* fmin = praxis(1.e-5,macheps, h, n, prin, x, func); */
+printf(" ftol=%lg, macheps=%lg, h0=%lg, npar=%d, prin=%d \n",ftol, macheps, h0, npar, prin);
+fprintf(ficlog," ftol=%lg, macheps=%lg, h0=%lg, npar=%d, prin=%d \n",ftol, macheps, h0, npar, prin);
ffmin = praxis(ftol,macheps, h0, npar, prin, p, func);
printf("End Praxis\n");
#endif /* FLATSUP */
fprintf(ficgp,"\nset out;unset log\n");
/* fprintf(ficgp,"\nset out \"%s.svg\"; replot; set out; # bug gnuplot",subdirf2(optionfilefiname,"ILK_")); */
+ fprintf(ficlog," withing printinggnuplot\n");fflush(ficlog);
strcpy(dirfileres,optionfilefiname);
/* } /\* k1 *\/ */
} /* cpt */
-
+ fprintf(ficlog," withing printinggnuplot 2 eme\n");fflush(ficlog);
+
/*2 eme*/
/* for (k1=1; k1<= m ; k1 ++){ */
for(nres=1; nres <= nresult; nres++){ /* For each resultline */
} /* end nres */
/* } /\* k1 end 2 eme*\/ */
+ fprintf(ficlog," withing printinggnuplot 3 eme\n");fflush(ficlog);
/*3eme*/
/* for (k1=1; k1<= m ; k1 ++){ */
}
} /* end nres */
/* } /\* end kl 3eme *\/ */
+ fprintf(ficlog," withing printinggnuplot 4 eme\n");fflush(ficlog);
/* 4eme */
/* Survival functions (period) from state i in state j by initial state i */
} /* end cpt state*/
} /* end nres */
/* } /\* end covariate k1 *\/ */
+ fprintf(ficlog," withing printinggnuplot 5 eme\n");fflush(ficlog);
/* 5eme */
/* Survival functions (period) from state i in state j by final state j */
/* } /\* end covariate *\/ */
} /* end nres */
+ fprintf(ficlog," withing printinggnuplot 6 eme\n");fflush(ficlog);
/* 6eme */
/* CV preval stable (period) for each covariate */
/* for (k1=1; k1<= m ; k1 ++) /\* For each covariate combination if any *\/ */
} /* end cpt state*/
} /* end covariate */
+ fprintf(ficlog," withing printinggnuplot 7 eme\n");fflush(ficlog);
/* 7eme */
if(prevbcast == 1){
fprintf(ficgp,", '' ");
/* l=(nlstate+ndeath)*(i-1)+1; */
l=(nlstate+ndeath)*(cpt-1)+1; /* fixed for i; cpt=1 1, cpt=2 1+ nlstate+ndeath, 1+2*(nlstate+ndeath) */
- /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a veÌ\81rifier *\/ */
- /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a veÌ\81rifier *\/ */
+ /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a veÃ\8cÂ\81rifier *\/ */
+ /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a veÃ\8cÂ\81rifier *\/ */
fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+i-1); /* To be verified */
/* for (j=2; j<= nlstate ; j ++) */
/* fprintf(ficgp,"+$%d",k+l+j-1); */
} /* end cpt state*/
} /* end covariate */
} /* End if prevbcast */
+ fprintf(ficlog," withing printinggnuplot 8 eme\n");fflush(ficlog);
/* 8eme */
if(prevfcast==1){
} /* end covariate */
} /* End if prevbcast */
+ fprintf(ficlog," withing printinggnuplot92 eme\n");fflush(ficlog);
/* 9eme writing MLE parameters */
fprintf(ficgp,"\n##############\n#9eme MLE estimated parameters\n#############\n");
}
}
fprintf(ficgp,"##############\n#\n");
+ fprintf(ficlog," withing printinggnuplot 10 eme\n");fflush(ficlog);
/*goto avoid;*/
/* 10eme Graphics of probabilities or incidences using written MLE parameters */
if(cptcovdageprod >0){
/* if(j==Tprod[ijp]) { */ /* not necessary */
/* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
- if(ijp <=cptcovprod) { /* Product Vn*Vm and age*VN*Vm*/
- if(DummyV[Tvardk[ijp][1]]==0){/* Vn is dummy */
- if(DummyV[Tvardk[ijp][2]]==0){/* Vn and Vm are dummy */
- /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
- fprintf(ficgp,"+p%d*%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]);
- }else{ /* Vn is dummy and Vm is quanti */
- /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
- fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]);
- }
- }else{ /* age* Vn*Vm Vn is quanti HERE */
- if(DummyV[Tvard[ijp][2]]==0){
- fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[ijp][2]],Tqinvresult[nres][Tvardk[ijp][1]]);
- }else{ /* Both quanti */
- fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]);
- }
+ if(ijp <=cptcovprod) { /* Product Vn*Vm and age*VN*Vm*/
+ /*for (i = 1; i <= cptcovt; i++)printf("DEB i=%d,Tvardk[i][1]=%d, Tvardk[i][2]=%d\n", i, Tvardk[i][1], Tvardk[i][2]);*/
+ /*for (i = 1; i <= cptcovprod; i++)printf("DEB i=%d,Tvard[i][1]=%d, Tvard[i][2]=%d\n", i, Tvard[i][1], Tvard[i][2]);*/
+ /* printf("DEBUG ijp=%d Tvard[ijp][1]=%d \n", ijp, Tvard[ijp][1]);*/
+ if(DummyV[Tvardk[j][1]]==0){/* Vn is dummy */
+ if(DummyV[Tvardk[j][2]]==0){/* Vn and Vm are dummy */
+ /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
+ fprintf(ficgp,"+p%d*%d*%d*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[j][1]],Tinvresult[nres][Tvardk[j][2]]);
+ }else{ /* Vn is dummy and Vm is quanti */
+ /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+ fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[j][1]],Tqinvresult[nres][Tvardk[j][2]]);
+ }
+ }else{ /* age* Vn*Vm Vn is quanti HERE */
+ if(DummyV[Tvard[ijp][2]]==0){
+ fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvardk[j][2]],Tqinvresult[nres][Tvardk[j][1]]);
+ }else{ /* Both quanti */
+ fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvardk[j][1]],Tqinvresult[nres][Tvardk[j][2]]);
}
- ijp++;
}
- /* } */ /* end Tprod */
+ ijp++;
+ }
+ /* } */ /* end Tprod */
}
break;
case 0:
case 3:
if(cptcovdageprod >0){
/* if(j==Tprod[ijp]) { /\* *\/ */
- /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
- if(ijp <=cptcovprod) { /* Product */
- if(DummyV[Tvardk[ijp][1]]==0){/* Vn is dummy */
- if(DummyV[Tvardk[ijp][2]]==0){/* Vn and Vm are dummy */
- /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
- fprintf(ficgp,"+p%d*%d*%d*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][1]],Tinvresult[nres][Tvardk[ijp][2]]);
- /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); */
- }else{ /* Vn is dummy and Vm is quanti */
- /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
- fprintf(ficgp,"+p%d*%d*%f*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]);
- /* fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */
- }
- }else{ /* Vn*Vm Vn is quanti */
- if(DummyV[Tvardk[ijp][2]]==0){
- fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[ijp][2]],Tqinvresult[nres][Tvardk[ijp][1]]);
- /* fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); */
- }else{ /* Both quanti */
- fprintf(ficgp,"+p%d*%f*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvardk[ijp][1]],Tqinvresult[nres][Tvardk[ijp][2]]);
- /* fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */
- }
+ /* printf("Tprod[%d]=%d, j=%d\n", ij, Tprod[ijp], j); */
+ if(ijp <=cptcovprod) { /* Product */
+ if(DummyV[Tvardk[j][1]]==0){/* Vn is dummy */
+ if(DummyV[Tvardk[j][2]]==0){/* Vn and Vm are dummy */
+ /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],nbcode[Tvard[ijp][2]][codtabm(k1,j)]); */
+ fprintf(ficgp,"+p%d*%d*%d*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[j][1]],Tinvresult[nres][Tvardk[j][2]]);
+ /* fprintf(ficgp,"+p%d*%d*%d",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tinvresult[nres][Tvard[ijp][2]]); */
+ }else{ /* Vn is dummy and Vm is quanti */
+ /* fprintf(ficgp,"+p%d*%d*%f",i+j+2+nagesqr-1,nbcode[Tvard[ijp][1]][codtabm(k1,j)],Tqinvresult[nres][Tvard[ijp][2]]); */
+ fprintf(ficgp,"+p%d*%d*%f*x",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[j][1]],Tqinvresult[nres][Tvardk[j][2]]);
+ /* fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */
}
- ijp++;
+ }else{ /* Vn*Vm Vn is quanti */
+ if(DummyV[Tvardk[ijp][2]]==0){
+ fprintf(ficgp,"+p%d*%d*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tinvresult[nres][Tvardk[j][2]],Tqinvresult[nres][Tvardk[j][1]]);
+ /* fprintf(ficgp,"+p%d*%d*%f*x",i+j+2+nagesqr-1,Tinvresult[nres][Tvard[ijp][2]],Tqinvresult[nres][Tvard[ijp][1]]); */
+ }else{ /* Both quanti */
+ fprintf(ficgp,"+p%d*%f*%f",k3+(cpt-1)*ncovmodel+1+j+nagesqr,Tqinvresult[nres][Tvardk[j][1]],Tqinvresult[nres][Tvardk[j][2]]);
+ /* fprintf(ficgp,"+p%d*%f*%f*x",i+j+2+nagesqr-1,Tqinvresult[nres][Tvard[ijp][1]],Tqinvresult[nres][Tvard[ijp][2]]); */
+ }
}
+ ijp++;
+ }
/* } /\* end Tprod *\/ */
} /* end if */
break;
} /* end resultline */
} /* end ng */
/* avoid: */
- fflush(ficgp);
+ fflush(ficgp);
+ fprintf(ficlog," end printinggnuplot\n");fflush(ficlog);
+
} /* end gnuplot */
return 1;
}
for(k=1; k<=j;k++){ /* Loop on any covariate of the RESULT LINE */
+ cutl(stra, strb, resultsav, ' '); /* keeps in strb after the first ' ' (stra is the rest of the resultline to be analyzed in the next loop *//* resultsav= "V4=1 V5=25.1 V3=0" stra= "V5=25.1 V3=0" strb= "V4=1" */
if(nbocc(resultsav,'=') >1){
- cutl(stra,strb,resultsav,' '); /* keeps in strb after the first ' ' (stra is the rest of the resultline to be analyzed in the next loop *//* resultsav= "V4=1 V5=25.1 V3=0" stra= "V5=25.1 V3=0" strb= "V4=1" */
/* If resultsav= "V4= 1 V5=25.1 V3=0" with a blank then strb="V4=" and stra="1 V5=25.1 V3=0" */
cutl(strc,strd,strb,'='); /* strb:"V4=1" strc="1" strd="V4" */
/* If a blank, then strc="V4=" and strd='\0' */
#if defined __INTEL_COMPILER
#if defined(__GNUC__)
- struct utsname sysInfo; /* For Intel on Linux and OS/X */
+ struct utsname sysInfo; /* For Intel on Linux and OS/X but not not on linux mingw-w64*/
#endif
#elif defined(__GNUC__)
#ifndef __APPLE__
-#include <gnu/libc-version.h> /* Only on gnu */
+#ifndef __MINGW32__
+#include <gnu/libc-version.h> /* Only on gnu */
+/* #include <features.h> */ /* Only on gnu */
+ #endif
#endif
+#ifndef __MINGW32__
struct utsname sysInfo;
int cross = CROSS;
if (cross){
if(logged) fprintf(ficlog, "Cross-");
}
#endif
+#endif
printf("Compiled with:");if(logged)fprintf(ficlog,"Compiled with:");
#if defined(__clang__)
#endif
#if defined(__GNUC__) || defined(__GNUG__)
printf(" GNU GCC/G++");if(logged)fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */
+#if defined(__MINGW32__)
+ printf(" with MINGW");if(logged)fprintf(ficlog," with MINGW");/* With MINGW. --------------------------------------------- */
+#endif
#endif
#if defined(__HP_cc) || defined(__HP_aCC)
printf(" Hewlett-Packard C/aC++");if(logged)fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */
printf(" using GNU C version %d.\n", __GNUC_VERSION__);
if(logged) fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__);
+#ifndef __MINGW32__
if (uname(&sysInfo) != -1) {
printf("Running on: %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
if(logged) fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
}
else
perror("uname() error");
+#endif
//#ifndef __INTEL_COMPILER
-#if !defined (__INTEL_COMPILER) && !defined(__APPLE__)
+#if !defined (__INTEL_COMPILER) && !defined(__APPLE__)&& !defined(__MINGW32__) /* MING hasn't yet libc version */
printf("GNU libc version: %s\n", gnu_get_libc_version());
if(logged) fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version());
#endif
<title>IMaCh %s</title></head>\n\
<body><font size=\"7\"><a href=http:/euroreves.ined.fr/imach>IMaCh for Interpolated Markov Chain</a> </font><br>\n\
<font size=\"3\">Sponsored by Copyright (C) 2002-2015 <a href=http://www.ined.fr>INED</a>\
--EUROREVES-Institut de longévité-2013-2022-Japan Society for the Promotion of Sciences 日本学術振興会 \
+-EUROREVES-Institut de longeÌ\81viteÌ\81-2013-2022-Japan Society for the Promotion of Sciences 日本å¦è¡“振興会 \
(<a href=https://www.jsps.go.jp/english/e-grants/>Grant-in-Aid for Scientific Research 25293121</a>) - \
<a href=https://software.intel.com/en-us>Intel Software 2015-2018</a></font><br> \n", optionfilehtm);
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d, ftolpl=%e\n",ageminpar,agemaxpar,bage,fage, estepm, ftolpl);
-
+ fflush(ficlog);
/* Other stuffs, more or less useful */
while(fgets(line, MAXLINE, ficpar)) {
/* If line starts with a # it is a comment */
fprintf(ficres,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
fprintf(ficlog,"begin-prev-date=%.lf/%.lf/%.lf end-prev-date=%.lf/%.lf/%.lf mov_average=%d\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2,mobilav);
}
-
+ fflush(ficlog);
while(fgets(line, MAXLINE, ficpar)) {
/* If line starts with a # it is a comment */
if (line[0] == '#') {
break;
}
-
+ fflush(ficlog);
dateprev1=anprev1+(mprev1-1)/12.+(jprev1-1)/365.;
dateprev2=anprev2+(mprev2-1)/12.+(jprev2-1)/365.;
fprintf(ficparo,"pop_based=%d\n",popbased);
fprintf(ficres,"pop_based=%d\n",popbased);
}
-
+ fflush(ficlog);
/* Results */
/* 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 */
fprintf(ficlog,"ERROR: Current version of IMaCh limits the number of resultlines to %d, you used %d\nYou can use the 'r' parameter file '%s' which uses option mle=0 to get other results. ",MAXRESULTLINESPONE-1,nresult,rfileres);
goto end;
}
-
+ fflush(ficlog);
if(!decoderesult(resultline, nresult)){ /* Fills TKresult[nresult] combination and Tresult[nresult][k4+1] combination values */
fprintf(ficparo,"result: %s\n",resultline);
fprintf(ficres,"result: %s\n",resultline);
decoderesult(".",nresult ); /* No covariate */
} /* End switch parameterline */
}while(endishere==0); /* End do */
-
+ fflush(ficlog);
/* freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvaraff,nbcode, ncodemax,mint,anint); */
/* ,dateprev1,dateprev2,jprev1, mprev1,anprev1,jprev2, mprev2,anprev2); */
datebackf=dateintmean-yrbproj;
date2dmy(datebackf,&jbackf, &mbackf, &anbackf);
}
-
+ fprintf(ficlog," before printinggnuplot\n");fflush(ficlog);
printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, prevbcast, pathc,p, (int)anprojd-bage, (int)anbackd-fage);/* HERE valgrind Tvard*/
}
printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,prevbcast, estepm, \
jprev1,mprev1,anprev1,dateprev1, dateprojd, datebackd,jprev2,mprev2,anprev2,dateprev2,dateprojf, datebackf);
-
+ fprintf(ficlog," after printinghtml\n");fflush(ficlog);
/*------------ free_vector -------------*/
/* chdir(path); */
/* Other results (useful)*/
+ fprintf(ficlog," before prevalence_limit\n");fflush(ficlog);
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */
#ifdef WIN32
if (_chdir(pathcd) != 0)
printf("Can't move to directory %s!\n",path);
+ /*if(_getcwd(pathcd,MAXLINE) > 0)*/
if(_getcwd(pathcd,MAXLINE) > 0)
#else
if(chdir(pathcd) != 0)
/*strcat(plotcmd,CHARSEPARATOR);*/
sprintf(plotcmd,"gnuplot");
#ifdef _WIN32
- sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);
+ sprintf(plotcmd,"\"gnuplot.exe\"");
+ /*sprintf(plotcmd,"\"%sgnuplot.exe\"",pathimach);*/ /* If gnuplot is in the path */
+ printf(" Win32 plotcmd=%s\n",plotcmd);
#endif
if(!stat(plotcmd,&info)){
printf("Error or gnuplot program not found: '%s'\n",plotcmd);fflush(stdout);
strcpy(plotcmd,GNUPLOTPROGRAM);
if(!stat(plotcmd,&info)){
printf("Error gnuplot program not found: '%s'\n",plotcmd);fflush(stdout);
- }else
+ }else{
strcpy(pplotcmd,plotcmd);
+ }
+ printf(" _unix pplotcmd=%s\n",pplotcmd);
#endif
}else
strcpy(pplotcmd,plotcmd);
if((outcmd=system(plotcmd)) != 0){
printf("Error in 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");
+#ifdef _WIN32
+ sprintf(plotcmd,"%sgnuplot.exe %s", pathimach, optionfilegnuplot);
+#else
sprintf(plotcmd,"%sgnuplot %s", pathimach, optionfilegnuplot);
+#endif
if((outcmd=system(plotcmd)) != 0){
printf("\n Still a problem with gnuplot command %s, err=%d\n", plotcmd, outcmd);
strcpy(plotcmd,pplotcmd);
git://henry.ined.fr / .git / commitdiff