Added isotonic_regression.

dlib/optimization.h

@@ -16,6 +16,7 @@
 #include "optimization/find_max_factor_graph_nmplp.h"
 #include "optimization/find_max_factor_graph_viterbi.h"
 #include "optimization/find_max_parse_cky.h"
+#include "optimization/isotonic_regression.h"
 
 #endif // DLIB_OPTIMIZATIOn_HEADER
 

dlib/optimization/isotonic_regression.h

+// Copyright (C) 2018  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_ISOTONIC_ReGRESSION_H_
+#define DLIB_ISOTONIC_ReGRESSION_H_
+
+#include "isotonic_regression_abstract.h"
+#include <vector>
+
+namespace dlib
+{
+    class isotonic_regression
+    {
+    public:
+
+        template <
+            typename const_iterator, 
+            typename iterator
+            >
+        void operator() (
+            const_iterator begin,
+            const_iterator end,
+            iterator obegin
+        )
+        {
+            for (auto i = begin; i != end; ++i)
+            {
+                blocks.emplace_back(*i);
+                while (blocks.size() > 1 && prev_block().avg > current_block().avg)
+                {
+                    // merge the last two blocks.
+                    prev_block() = prev_block() + current_block();
+                    blocks.pop_back();
+                }
+            }
+
+            // unpack results to output
+            for (auto& block : blocks)
+            {
+                for (size_t k = 0; k < block.num; ++k)
+                    *obegin++ = block.avg;
+            }
+
+            blocks.clear();
+        }
+
+        void operator() (
+            std::vector<double>& vect
+        ) { (*this)(vect.begin(), vect.end(), vect.begin()); }
+
+    private:
+
+
+        struct block_t
+        {
+            block_t(double val) : num(1), avg(val) {}
+            block_t(size_t n, double val) : num(n), avg(val) {}
+
+            size_t num;
+            double avg;
+
+            inline block_t operator+(const block_t& rhs) const
+            {
+                return block_t(num+rhs.num,
+                    (num*avg + rhs.num*rhs.avg)/(num+rhs.num));
+            }
+        };
+
+        inline block_t& prev_block() { return blocks[blocks.size()-2]; }
+        inline block_t& current_block() { return blocks.back(); }
+
+        std::vector<block_t> blocks;
+    };
+}
+
+#endif // DLIB_ISOTONIC_ReGRESSION_H_
+
+

dlib/optimization/isotonic_regression_abstract.h

+// Copyright (C) 2018  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_ISOTONIC_ReGRESSION_ABSTRACT_H_
+#ifdef DLIB_ISOTONIC_ReGRESSION_ABSTRACT_H_
+
+#include <vector>
+
+namespace dlib
+{
+    class isotonic_regression
+    {
+        /*!
+            WHAT THIS OBJECT REPRESENTS
+                This object is a tool for performing 1-D isotonic regression. That is, it
+                finds the least squares fit of a non-parametric curve to some user supplied
+                data, subject to the constraint that the fitted curve is non-decreasing.
+
+                This is done using the fast O(n) pool adjacent violators algorithm.
+        !*/
+
+    public:
+
+        template <
+            typename const_iterator, 
+            typename iterator
+            >
+        void operator() (
+            const_iterator begin,
+            const_iterator end,
+            iterator obegin
+        );
+        /*!
+            requires
+                - [begin,end) is an iterator range of float or doubles.
+                - obegin points to an iterator range at least std::distance(begin,end) in
+                  size which is capable of storing double or float values.
+            ensures
+                - Given the range of real values stored in [begin,end), this method performs isotonic regression
+                  on this data and writes the results to obegin.  To be specific:
+                    - let IN refer to the input values stored in the iterator range [begin,end).
+                    - let OUT refer to the output values stored in the iterator range [obegin, obegin+std::distance(begin,end)).
+                    - This function populates OUT with values such that the sum_i of
+                      (IN[i]-OUT[i])^2 is minimized, subject to the constraint that 
+                      OUT[i] <= OUT[i+1], i.e. that OUT is monotonic.
+                - It is OK for [begin,end) to overlap with the range pointed to by obegin.
+                  That is, this function can run in-place.
+        !*/
+
+        void operator() (
+            std::vector<double>& vect
+        ) { (*this)(vect.begin(), vect.end(), vect.begin()); }
+        /*!
+            ensures
+                - performs in-place isotonic regression.  Therefore, #vect will contain the
+                  isotonic regression of vect.
+                - #vect.size() == vect.size()
+        !*/
+
+    };
+}
+
+#endif // DLIB_ISOTONIC_ReGRESSION_ABSTRACT_H_
+
+
+

dlib/test/CMakeLists.txt

@@ -73,6 +73,7 @@ set (tests
    image.cpp
    iosockstream.cpp
    is_same_object.cpp
+   isotonic_regression.cpp
    kcentroid.cpp
    kernel_matrix.cpp
    kmeans.cpp

dlib/test/isotonic_regression.cpp

+// Copyright (C) 2018  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+
+
+#include <dlib/optimization.h>
+#include <dlib/global_optimization.h>
+#include <sstream>
+#include <string>
+#include <cstdlib>
+#include <ctime>
+#include <vector>
+
+#include "tester.h"
+
+
+namespace  
+{
+
+    using namespace test;
+    using namespace dlib;
+    using namespace std;
+
+    logger dlog("test.isotonic_regression");
+
+// ----------------------------------------------------------------------------------------
+
+    class optimization_tester : public tester
+    {
+    public:
+        optimization_tester (
+        ) :
+            tester ("test_isotonic_regression",
+                    "Runs tests on the isotonic_regression object.")
+        {}
+
+        void perform_test (
+        )
+        {
+            dlib::rand rnd;
+
+            for (int round = 0; round < 100; ++round)
+            {
+                print_spinner();
+                std::vector<double> vect;
+                for (int i = 0; i < 5; ++i)
+                    vect.push_back(put_in_range(-1,1,rnd.get_random_gaussian()));
+
+
+                auto f = [&](const matrix<double,0,1>& x)
+                {
+                    double dist = 0;
+                    double sum = 0;
+                    for (long i = 0; i < x.size(); ++i)
+                    {
+                        sum += x(i);
+                        dist += (sum-vect[i])*(sum-vect[i]);
+                    }
+                    return dist;
+                };
+
+                auto objval = [vect](const matrix<double,0,1>& x)
+                {
+                    return sum(squared(mat(vect)-x));
+                };
+
+                auto is_monotonic = [](const matrix<double,0,1>& x)
+                {
+                    for (long i = 1; i < x.size(); ++i)
+                    {
+                        if (x(i-1) > x(i))
+                            return false;
+                    }
+                    return true;
+                };
+
+                matrix<double,0,1> lower(5), upper(5);
+                lower = 0;
+                lower(0) = -4;
+                upper = 4;
+                // find the solution with find_min_global() and then check that it matches
+                auto result = find_min_global(f, lower, upper, max_function_calls(40));
+
+                for (long i = 1; i < result.x.size(); ++i)
+                    result.x(i) += result.x(i-1);
+
+                isotonic_regression mr;
+                mr(vect);
+
+                dlog << LINFO << "err: "<<  objval(mat(vect)) - objval(result.x);
+
+                DLIB_CASSERT(is_monotonic(mat(vect)));
+                DLIB_CASSERT(is_monotonic(result.x));
+                // isotonic_regression should be at least as good as find_min_global().
+                DLIB_CASSERT(objval(mat(vect)) - objval(result.x) < 1e-13);
+            }
+
+        }
+    } a;
+
+}
+
+
+