Simplify step size choice.

1 files changed, 12 insertions, 26 deletions

diff --git a/init.c b/init.c
index ff995b6..c512cca 100644
--- a/init.c
+++ b/init.c
@@ -343,8 +343,7 @@ int td_solve(TDContext *td, double *coeffs_init[3])
         const double delta = 1e-5;
         const double a1 = a0 - delta;
 
-        double cur_amplitude, new_amplitude, inverse_step;
-        int dir;
+        double cur_amplitude, new_amplitude, step;
 
         ret = constraint_eval_alloc(td, a0, &ce);
         if (ret < 0)
@@ -382,16 +381,11 @@ int td_solve(TDContext *td, double *coeffs_init[3])
         }
 
         cur_amplitude = a1;
-        dir = SGN(td->amplitude - cur_amplitude);
-        inverse_step = 1.0 / td->amplitude - 1.0 / cur_amplitude;
-        while (fabs(cur_amplitude - td->amplitude) > DBL_EPSILON) {
-            //double scale_factor;
-            //scale_factor  = cur_amplitude;
-            //cur_amplitude = MIN(td->amplitude, 2 * cur_amplitude);
-            //scale_factor  = cur_amplitude / scale_factor;
+        step = (td->amplitude - cur_amplitude) / 128.0;
 
-            new_amplitude = 1.0 / ((1.0 / cur_amplitude) + inverse_step);
-            if (dir * (td->amplitude - new_amplitude) < 0)
+        while (fabs(cur_amplitude - td->amplitude) > DBL_EPSILON) {
+            new_amplitude = cur_amplitude + step;
+            if (SGN(step) != SGN(td->amplitude - new_amplitude))
                 new_amplitude = td->amplitude;
             tdi_log(&s->logger, 2, "Trying amplitude %g\n", new_amplitude);
 
@@ -404,27 +398,19 @@ int td_solve(TDContext *td, double *coeffs_init[3])
                                     ce, s->coeffs_tmp, 1);
             tdi_constraint_eval_free(&ce);
             if (ret == -EDOM) {
-                inverse_step = 0.5 * inverse_step;
-                //if (fabs(inverse_step) < 1e-2)
-                //    return ret;
+                step *= 0.5;
+                if (fabs(step) < 1e-5)
+                    return ret;
                 continue;
             } else if (ret < 0)
                 return ret;
+
+            if (ret <= s->solver->maxiter / 2)
+                step *= 1.75;
+
             cur_amplitude = new_amplitude;
             memcpy(s->coeffs, s->coeffs_tmp, sizeof(*s->coeffs) * N);
         }
-        //while (1) {
-
-
-        //    if (fabs(cur_amplitude - td->amplitude) < DBL_EPSILON)
-        //        goto finish;
-
-        //    cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], SQR(scale_factor), td->coeffs[0], 1);
-        //    cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], scale_factor, td->coeffs[1], 1);
-        //    cblas_dscal(td->nb_coeffs[0] * td->nb_coeffs[1], scale_factor, td->coeffs[2], 1);
-
-        //}
-
     }
 finish:
     tdi_nlsolve_print_stats(s->solver);