Reverse the role of p_e and a in optimization

The p_e can be calculated for a fixed a without logarithmization, so the
minimum should be the true minimum. Also some other minor changes to the
function around were made.

diff --git a/Measurement.h b/Measurement.h
index 85a32ea0feae241ca533d4f0044baf92de3929c0..4b5c596a57495bc2731f7c2a2d5a3b24e81139bb 100644 (file)
--- a/Measurement.h
+++ b/Measurement.h
@@ -17,10 +17,11 @@ class Measurement
                size_t                                  m_size;
 
                float   calc_equilibrium_dist() const;
-               float   calc_decay_speed(float equilibrium_pressure) const;
-               float   get_error(float equilibrium_pressure) const;
-               float   get_error_diff(float equilibrium_pressure, float step = 0.001) const;
-               float   get_equilibrium_pressure(size_t max_iter = 10) const;
+               float   calc_equilibrium_pressure(float equilibrium_pressure) const;
+               float   get_error(float decay_speed) const;
+               float   get_error_diff(float decay_speed, float step) const;
+               float   get_decay_speed(size_t max_iter = 20) const;
+               float   screen_for_minimum(float upper_bound) const;
 
        public:
                Measurement() = delete;

diff --git a/Measurement.ino b/Measurement.ino
index 082a45acea099eba3ee9d86ae36352bec1f6867b..6eadc2d6f42486482767bb29c8cd2714a632396b 100644 (file)
--- a/Measurement.ino
+++ b/Measurement.ino
@@ -61,16 +61,18 @@ void        Measurement::clear()
 // exp(-a * t) >= threshold
 //
 // The model has 2 parameters, and it is difficult to get them analytically
-// using the least squares method. But getting the optimal "a" parameter for a
-// given p_e is quite easy using least squares, and it gives this:
+// using the least squares method. But getting the optimal "p_e" parameter for a
+// given decay_speed "a" is quite easy using least squares, and it gives this:
 //
-// a(p_e) = sum(ln((p_i - p_e) / (p_0 - p_e)) * t_i) / sum((t_i)^2)
+// p_e(a) = sum((p_i - p_0 * exp_i) * (1 - exp_i)) / sum((1 - exp_i)^2)
 //
-// Using this we can then numerically optimize p_e. Interval halving is used
-// because it is foolproof, the precision needed is low and the evaluation is
-// not computationally demanding.
+// where exp_i = exp(-a * t_i).
+//
+// Using this we can then numerically optimize a. Interval halving is used
+// because it is foolproof and the evaluation is not computationally demanding.
 void   Measurement::wait_for_equilibrium(float threshold = 0.01)
 {
+       clear();
        while (m_size < 20 || threshold < calc_equilibrium_dist())
        {
                get_measurement();
@@ -80,67 +82,107 @@ void       Measurement::wait_for_equilibrium(float threshold = 0.01)
 
 float  Measurement::calc_equilibrium_dist() const
 {
-       float   decay_speed{calc_decay_speed(get_equilibrium_pressure())};
-
-       return (exp(-decay_speed * m_time[m_size - 1]));
+       return (exp(-get_decay_speed() * m_time[m_size - 1]));
 }
 
-float  Measurement::calc_decay_speed(float equilibrium_pressure) const
+float  Measurement::calc_equilibrium_pressure(float decay_speed) const
 {
-       float   time_squares_sum{0};
-       float   time_logp_sum{0};
+       float   numerator_sum{0};
+       float   denominator_sum{0};
+       float   exponential_term{};
 
        for (size_t i{0}; i < m_size; ++i)
        {
-               time_squares_sum += static_cast<float>(m_time[i]) * m_time[i];
-               time_logp_sum += m_time[i] * log((m_pressure[i] - equilibrium_pressure)
-                               / (m_pressure[0] - equilibrium_pressure));
+               exponential_term = exp(-decay_speed * m_time[i]);
+               numerator_sum += (m_pressure[i] - m_pressure[0] * exponential_term)
+                       * (1 - exponential_term);
+               denominator_sum += (1 - exponential_term) * (1 - exponential_term);
        }
-       return (- time_logp_sum / time_squares_sum);
+       return (numerator_sum / denominator_sum);
 }
 
-float  Measurement::get_error(float eq_pressure) const
+float  Measurement::get_error(float decay_speed) const
 {
-       float   decay_speed{calc_decay_speed(eq_pressure)};
+       float   eq_pressure{calc_equilibrium_pressure(decay_speed)};
        float   residual{};
        float   error{0};
 
-       for (size_t i{0}; i < m_size; ++i)
+       if (decay_speed > 0)
        {
-               residual = m_pressure[i] - (m_pressure[0] - eq_pressure) * exp(- decay_speed * m_time[i]) - eq_pressure;
-               error += residual * residual;
+               for (size_t i{0}; i < m_size; ++i)
+               {
+                       residual = m_pressure[i] - (m_pressure[0] - eq_pressure) * exp(- decay_speed * m_time[i]) - eq_pressure;
+                       error += residual * residual;
+               }
+       }
+       else
+       {
+               for (size_t i{0}; i < m_size; ++i)
+               {
+                       residual = m_pressure[i] - eq_pressure;
+                       error += residual * residual;
+               }
        }
        return (error);
 }
 
-float  Measurement::get_error_diff(float eq_pressure, float step) const
+float  Measurement::get_error_diff(float decay_speed, float step) const
 {
        float   error_minus{};
        float   error_plus{};
 
-       error_plus = get_error(eq_pressure + step);
-       if (eq_pressure < step)
+       error_plus = get_error(decay_speed + step);
+       if (decay_speed < step)
        {
                error_minus = get_error(0);
-               return ((error_plus - error_minus) / (step + eq_pressure));
+               return ((error_plus - error_minus) / (step + decay_speed));
        }
-       error_minus = get_error(eq_pressure - step);
+       error_minus = get_error(decay_speed - step);
        return ((error_plus - error_minus) / (2 * step));
 }
 
-float  Measurement::get_equilibrium_pressure(size_t max_iter) const
+float  Measurement::screen_for_minimum(float upper_bound) const
 {
-       float   min_arg{0};
-       float   min_val{get_error_diff(min_arg)};
-       float   max_arg{1.1 * m_sensor.m_max_pressure};
-       float   max_val{get_error_diff(max_arg)};
+       float   test_arg{upper_bound};
+       float   last_arg{};
+       float   test_value{};
+       float   last_value{};
+
+       test_value = get_error(test_arg);
+       do
+       {
+               last_value = test_value;
+               last_arg = test_arg;
+               test_arg /= 10;
+               test_value = get_error(test_arg);
+       }
+       while (test_value <= last_value);
+       return (last_arg);
+}
+
+float  Measurement::get_decay_speed(size_t max_iter) const
+{
+       float   min_arg{};
+       float   min_val{};
+       float   max_arg{};
+       float   max_val{};
        float   test{};
        float   test_val{};
 
+       min_arg = screen_for_minimum(1);
+       max_arg = min_arg;
+       do
+       {
+               min_arg /= 2;
+               max_arg *= 2;
+               min_val = get_error_diff(min_arg, min_arg / 1000);
+               max_val = get_error_diff(max_arg, max_arg / 1000);
+       }
+       while ((min_val > 0 && max_val > 0) || (min_val < 0 && max_val < 0));
        for (; max_iter > 0; --max_iter)
        {
                test = (max_arg + min_arg) / 2;
-               test_val = get_error_diff(test);
+               test_val = get_error_diff(test, test / 1000);
                if (test_val == 0)
                        return (test);
                else if ((test_val < 0 && min_val < 0) || (test_val > 0 && min_val > 0))