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 <limits.h>
#include <sys/types.h>
-/* #if defined(__GNUC__) */
-/* #include <sys/utsname.h> /\* Doesn't work on Windows *\/ */
-/* #endif */
+#if defined(__GNUC__) && !defined(_WIN32)
+#include <sys/utsname.h> /* Doesn't work on Windows */
+#endif
#include <sys/stat.h>
#include <errno.h>
/* $State$ */
#include "version.h"
char version[]=__IMACH_VERSION__;
-char copyright[]="April 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024";
+char copyright[]="November 2024,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2020, Nihon University 2021-202, INED 2000-2024";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
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) { /* Added for understanding small changes */
+ /* 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^2 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 */
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 */
/* 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*/
- /*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]]);
- }
+ /*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]]);
}
- ijp++;
+ }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 */
+ /* } */ /* 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[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]]); */
- }
- }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]]); */
- }
+ /* 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;
git://henry.ined.fr / .git / commitdiff