svm.h

// Copyright (C) 2007  Davis E. King (davis@dlib.net)
// License: Boost Software License   See LICENSE.txt for the full license.
#ifndef DLIB_SVm_
#define DLIB_SVm_

#include "svm_abstract.h"
#include <cmath>
#include <limits>
#include <sstream>
#include "../matrix.h"
#include "../algs.h"
#include "../serialize.h"
#include "../rand.h"
#include "../std_allocator.h"
#include "function.h"
#include "kernel.h"
#include "../enable_if.h"
#include "../optimization.h"

namespace dlib
{

// ----------------------------------------------------------------------------------------

    class invalid_svm_nu_error : public dlib::error 
    { 
    public: 
        invalid_svm_nu_error(const std::string& msg, double nu_) : dlib::error(msg), nu(nu_) {};
        const double nu;
    };

// ----------------------------------------------------------------------------------------

    template <
        typename T
        >
    typename T::type maximum_nu_impl (
        const T& y
    )
    {
        typedef typename T::type scalar_type;
        // make sure requires clause is not broken
        DLIB_ASSERT(y.nr() > 1 && y.nc() == 1,
            "\ttypedef T::type maximum_nu(y)"
            << "\n\ty should be a column vector with more than one entry"
            << "\n\ty.nr(): " << y.nr() 
            << "\n\ty.nc(): " << y.nc() 
            );

        long pos_count = 0;
        long neg_count = 0;
        for (long r = 0; r < y.nr(); ++r)
        {
            if (y(r) == 1.0)
            {
                ++pos_count;
            }
            else if (y(r) == -1.0)
            {
                ++neg_count;
            }
            else
            {
                // make sure requires clause is not broken
                DLIB_ASSERT(y(r) == -1.0 || y(r) == 1.0,
                       "\ttypedef T::type maximum_nu(y)"
                       << "\n\ty should contain only 1 and 0 entries"
                       << "\n\tr:    " << r 
                       << "\n\ty(r): " << y(r) 
                );
            }
        }
        return static_cast<scalar_type>(2.0*(scalar_type)std::min(pos_count,neg_count)/(scalar_type)y.nr());
    }

    template <
        typename T
        >
    typename T::type maximum_nu (
        const T& y
    )
    {
        return maximum_nu_impl(vector_to_matrix(y));
    }

    template <
        typename T
        >
    typename T::value_type maximum_nu (
        const T& y
    )
    {
        return maximum_nu_impl(vector_to_matrix(y));
    }

// ----------------------------------------------------------------------------------------

    template <
        typename T,
        typename U
        >
    bool is_binary_classification_problem_impl (
        const T& x,
        const U& x_labels
    )
    {
        bool seen_neg_class = false;
        bool seen_pos_class = false;
        if (x.nc() != 1 || x_labels.nc() != 1) return false; 
        if (x.nr() != x_labels.nr()) return false;
        if (x.nr() <= 1) return false;
        for (long r = 0; r < x_labels.nr(); ++r)
        {
            if (x_labels(r) != -1 && x_labels(r) != 1)
                return false;

            if (x_labels(r) == 1)
                seen_pos_class = true;
            if (x_labels(r) == -1)
                seen_neg_class = true;
        }

        return seen_pos_class && seen_neg_class;
    }

    template <
        typename T,
        typename U
        >
    bool is_binary_classification_problem (
        const T& x,
        const U& x_labels
    )
    {
        return is_binary_classification_problem_impl(vector_to_matrix(x), vector_to_matrix(x_labels));
    }

// ----------------------------------------------------------------------------------------

    template <
        typename K,
        typename sample_vector_type,
        typename scalar_vector_type
        >
    class kernel_matrix_cache
    {
    public:
        typedef float scalar_type;
        //typedef typename K::scalar_type scalar_type;
        typedef typename K::sample_type sample_type;
        typedef typename K::mem_manager_type mem_manager_type;

        const sample_vector_type& x;
        const scalar_vector_type& y;
        K kernel_function;

        mutable matrix<scalar_type,0,0,mem_manager_type> cache;
        mutable matrix<scalar_type,0,1,mem_manager_type> diag_cache;
        mutable matrix<long,0,1,mem_manager_type> lookup;
        mutable matrix<long,0,1,mem_manager_type> rlookup;
        mutable long next;

        /*!
        INITIAL VALUE
            - for all valid x:
                - lookup(x) == -1 
                - rlookup(x) == -1 

        CONVENTION
            - if (lookup(c) != -1) then
                - cache(lookup(c),*) == the cached column c of the kernel matrix
                - rlookup(lookup(c)) == c

            - if (rlookup(x) != -1) then
                - lookup(rlookup(x)) == x
                - cache(x,*) == the cached column rlookup(x) of the kernel matrix

            - next == the next row in the cache table to use to cache something 
        !*/

    public:
        kernel_matrix_cache (
            const sample_vector_type& x_,
            const scalar_vector_type& y_,
            K kernel_function_,
            long max_size_megabytes
        ) : x(x_), y(y_), kernel_function(kernel_function_) 
        {
            // figure out how many rows of the kernel matrix we can have
            // with the given amount of memory.
            long max_size = (max_size_megabytes*1024*1024)/(x.nr()*sizeof(scalar_type));
            // don't let it be 0
            if (max_size == 0)
                max_size = 1;
            long size = std::min(max_size,x.nr());

            diag_cache.set_size(x.nr(),1);
            cache.set_size(size,x.nr());
            lookup.set_size(x.nr(),1);
            rlookup.set_size(size,1);
            set_all_elements(lookup,-1);
            set_all_elements(rlookup,-1);
            next = 0;

            for (long i = 0; i < diag_cache.nr(); ++i)
                diag_cache(i) = kernel_function(x(i),x(i));
        }

        inline bool is_cached (
            long r
        ) const
        {
            return (lookup(r) != -1);
        }

        const scalar_type* col(long i) const 
        { 
            if (is_cached(i) == false)
                add_col_to_cache(i);

            // find where this column is in the cache
            long idx = lookup(i);
            if (idx == next)
            {
                // if this column was the next to be replaced
                // then make sure that doesn't happen
                next = (next + 1)%cache.nr();
            }

            return &cache(idx,0); 
        }
        const scalar_type* diag() const { return &diag_cache(0); }

        inline scalar_type operator () (
            long r,
            long c
        ) const
        {
            if (lookup(c) != -1)
            {
                return cache(lookup(c),r);
            }
            else if (r == c)
            {
                return diag_cache(r);
            }
            else if (lookup(r) != -1)
            {
                // the kernel is symmetric so this is legit
                return cache(lookup(r),c);
            }
            else
            {
                add_col_to_cache(c);
                return cache(lookup(c),r);
            }
        }

    private:
        void add_col_to_cache(
            long c
        ) const
        {
            // if the lookup table is pointing to cache(next,*) then clear lookup(next)
            if (rlookup(next) != -1)
                lookup(rlookup(next)) = -1;

            // make the lookup table so that it says c is now cached at the spot indicated by next
            lookup(c) = next;
            rlookup(next) = c;

            // compute this column in the kernel matrix and store it in the cache
            for (long i = 0; i < cache.nc(); ++i)
                cache(next,i) = y(c)*y(i)*kernel_function(x(c),x(i));

            next = (next + 1)%cache.nr();
        }

    };

// ----------------------------------------------------------------------------------------

    template <
        typename dec_funct_type,
        typename in_sample_vector_type,
        typename in_scalar_vector_type
        >
    const matrix<typename dec_funct_type::scalar_type, 1, 2, typename dec_funct_type::mem_manager_type> 
    test_binary_decision_function_impl (
        const dec_funct_type& dec_funct,
        const in_sample_vector_type& x_test,
        const in_scalar_vector_type& y_test
    )
    {
        typedef typename dec_funct_type::scalar_type scalar_type;
        typedef typename dec_funct_type::sample_type sample_type;
        typedef typename dec_funct_type::mem_manager_type mem_manager_type;
        typedef matrix<sample_type,0,1,mem_manager_type> sample_vector_type;
        typedef matrix<scalar_type,0,1,mem_manager_type> scalar_vector_type;

        // make sure requires clause is not broken
        DLIB_ASSERT( is_binary_classification_problem(x_test,y_test) == true,
                    "\tmatrix test_binary_decision_function()"
                    << "\n\t invalid inputs were given to this function"
                    << "\n\t is_binary_classification_problem(x_test,y_test): " 
                    << ((is_binary_classification_problem(x_test,y_test))? "true":"false"));


        // count the number of positive and negative examples
        long num_pos = 0;
        long num_neg = 0;


        long num_pos_correct = 0;
        long num_neg_correct = 0;


        // now test this trained object 
        for (long i = 0; i < x_test.nr(); ++i)
        {
            // if this is a positive example
            if (y_test(i) == +1.0)
            {
                ++num_pos;
                if (dec_funct(x_test(i)) >= 0)
                    ++num_pos_correct;
            }
            else if (y_test(i) == -1.0)
            {
                ++num_neg;
                if (dec_funct(x_test(i)) < 0)
                    ++num_neg_correct;
            }
            else
            {
                throw dlib::error("invalid input labels to the test_binary_decision_function() function");
            }
        }


        matrix<scalar_type, 1, 2, mem_manager_type> res;
        res(0) = (scalar_type)num_pos_correct/(scalar_type)(num_pos); 
        res(1) = (scalar_type)num_neg_correct/(scalar_type)(num_neg); 
        return res;
    }

    template <
        typename dec_funct_type,
        typename in_sample_vector_type,
        typename in_scalar_vector_type
        >
    const matrix<typename dec_funct_type::scalar_type, 1, 2, typename dec_funct_type::mem_manager_type> 
    test_binary_decision_function (
        const dec_funct_type& dec_funct,
        const in_sample_vector_type& x_test,
        const in_scalar_vector_type& y_test
    )
    {
        return test_binary_decision_function_impl(dec_funct,
                                 vector_to_matrix(x_test),
                                 vector_to_matrix(y_test));
    }

// ----------------------------------------------------------------------------------------

    template <
        typename trainer_type,
        typename in_sample_vector_type,
        typename in_scalar_vector_type
        >
    const matrix<typename trainer_type::scalar_type, 1, 2, typename trainer_type::mem_manager_type> 
    cross_validate_trainer_impl (
        const trainer_type& trainer,
        const in_sample_vector_type& x,
        const in_scalar_vector_type& y,
        const long folds
    )
    {
        typedef typename trainer_type::scalar_type scalar_type;
        typedef typename trainer_type::sample_type sample_type;
        typedef typename trainer_type::mem_manager_type mem_manager_type;
        typedef matrix<sample_type,0,1,mem_manager_type> sample_vector_type;
        typedef matrix<scalar_type,0,1,mem_manager_type> scalar_vector_type;

        // make sure requires clause is not broken
        DLIB_ASSERT(is_binary_classification_problem(x,y) == true &&
                    1 < folds && folds <= x.nr(),
            "\tmatrix cross_validate_trainer()"
            << "\n\t invalid inputs were given to this function"
            << "\n\t x.nr(): " << x.nr() 
            << "\n\t folds:  " << folds 
            << "\n\t is_binary_classification_problem(x,y): " << ((is_binary_classification_problem(x,y))? "true":"false")
            );


        // count the number of positive and negative examples
        long num_pos = 0;
        long num_neg = 0;
        for (long r = 0; r < y.nr(); ++r)
        {
            if (y(r) == +1.0)
                ++num_pos;
            else
                ++num_neg;
        }

        // figure out how many positive and negative examples we will have in each fold
        const long num_pos_test_samples = num_pos/folds; 
        const long num_pos_train_samples = num_pos - num_pos_test_samples; 
        const long num_neg_test_samples = num_neg/folds; 
        const long num_neg_train_samples = num_neg - num_neg_test_samples; 


        typename trainer_type::trained_function_type d;
        sample_vector_type x_test, x_train;
        scalar_vector_type y_test, y_train;
        x_test.set_size (num_pos_test_samples  + num_neg_test_samples);
        y_test.set_size (num_pos_test_samples  + num_neg_test_samples);
        x_train.set_size(num_pos_train_samples + num_neg_train_samples);
        y_train.set_size(num_pos_train_samples + num_neg_train_samples);

        long pos_idx = 0;
        long neg_idx = 0;

        matrix<scalar_type, 1, 2, mem_manager_type> res;
        set_all_elements(res,0);

        for (long i = 0; i < folds; ++i)
        {
            long cur = 0;

            // load up our positive test samples
            while (cur < num_pos_test_samples)
            {
                if (y(pos_idx) == +1.0)
                {
                    x_test(cur) = x(pos_idx);
                    y_test(cur) = +1.0;
                    ++cur;
                }
                pos_idx = (pos_idx+1)%x.nr();
            }

            // load up our negative test samples
            while (cur < x_test.nr())
            {
                if (y(neg_idx) == -1.0)
                {
                    x_test(cur) = x(neg_idx);
                    y_test(cur) = -1.0;
                    ++cur;
                }
                neg_idx = (neg_idx+1)%x.nr();
            }

            // load the training data from the data following whatever we loaded
            // as the testing data
            long train_pos_idx = pos_idx;
            long train_neg_idx = neg_idx;
            cur = 0;

            // load up our positive train samples
            while (cur < num_pos_train_samples)
            {
                if (y(train_pos_idx) == +1.0)
                {
                    x_train(cur) = x(train_pos_idx);
                    y_train(cur) = +1.0;
                    ++cur;
                }
                train_pos_idx = (train_pos_idx+1)%x.nr();
            }

            // load up our negative train samples
            while (cur < x_train.nr())
            {
                if (y(train_neg_idx) == -1.0)
                {
                    x_train(cur) = x(train_neg_idx);
                    y_train(cur) = -1.0;
                    ++cur;
                }
                train_neg_idx = (train_neg_idx+1)%x.nr();
            }

            try
            {
                // do the training and testing
                res += test_binary_decision_function(trainer.train(x_train,y_train),x_test,y_test);
            }
            catch (invalid_svm_nu_error&)
            {
                // Just ignore the error in this case since we are going to
                // interpret an invalid nu value the same as generating a decision
                // function that miss-classifies everything.
            }

        } // for (long i = 0; i < folds; ++i)

        return res/(scalar_type)folds;
    }

    template <
        typename trainer_type,
        typename in_sample_vector_type,
        typename in_scalar_vector_type
        >
    const matrix<typename trainer_type::scalar_type, 1, 2, typename trainer_type::mem_manager_type> 
    cross_validate_trainer (
        const trainer_type& trainer,
        const in_sample_vector_type& x,
        const in_scalar_vector_type& y,
        const long folds
    )
    {
        return cross_validate_trainer_impl(trainer,
                                           vector_to_matrix(x),
                                           vector_to_matrix(y),
                                           folds);
    }

// ----------------------------------------------------------------------------------------

    namespace prob_impl
    {
        template <typename vect_type>
        struct objective
        {
            objective (
                const vect_type& f_,
                const vect_type& t_
            ) : f(f_), t(t_) {}

            double operator() (
                const matrix<double,2,1>& x
            ) const
            {
                const double A = x(0);
                const double B = x(1);

                double res = 0;
                for (unsigned long i = 0; i < f.size(); ++i)
                {
                    const double val = A*f[i]+B;
                    // See the paper "A Note on Platt's Probabilistic Outputs for Support Vector Machines"
                    // for an explanation of why this code looks the way it does (rather than being the 
                    // obvious formula).
                    if (val < 0)
                        res += (t[i] - 1)*val + std::log(1 + std::exp(val));
                    else
                        res += t[i]*val + std::log(1 + std::exp(-val));
                }

                return res;
            }

            const vect_type& f;
            const vect_type& t;
        };

        template <typename vect_type>
        struct der
        {
            der (
                const vect_type& f_,
                const vect_type& t_
            ) : f(f_), t(t_) {}

            matrix<double,2,1> operator() (
                const matrix<double,2,1>& x
            ) const
            {
                const double A = x(0);
                const double B = x(1);

                double derA = 0;
                double derB = 0;

                for (unsigned long i = 0; i < f.size(); ++i)
                {
                    const double val = A*f[i]+B;
                    double p;
                    // compute p = 1/(1+exp(val)) 
                    // but do so in a way that avoids numerical overflow.
                    if (val < 0)
                        p = 1.0/(1 + std::exp(val));
                    else
                        p = std::exp(-val)/(1 + std::exp(-val));

                    derA += f[i]*(t[i] - p);
                    derB +=      (t[i] - p);
                }

                matrix<double,2,1> res;
                res = derA, derB;
                return res;
            }

            const vect_type& f;
            const vect_type& t;
        };

        template <typename vect_type>
        struct hessian 
        {
            hessian (
                const vect_type& f_,
                const vect_type& t_
            ) : f(f_), t(t_) {}

            matrix<double,2,2> operator() (
                const matrix<double,2,1>& x
            ) const
            {
                const double A = x(0);
                const double B = x(1);

                matrix<double,2,2> h;
                h = 0;

                for (unsigned long i = 0; i < f.size(); ++i)
                {
                    const double val = A*f[i]+B;
                    // compute pp = 1/(1+exp(val)) and
                    // compute pn = 1 - pp
                    // but do so in a way that avoids numerical overflow and catastrophic cancellation.
                    double pp, pn;
                    if (val < 0)
                    {
                        const double temp = std::exp(val);
                        pp = 1.0/(1 + temp);
                        pn = temp*pp; 
                    }
                    else
                    {
                        const double temp = std::exp(-val);
                        pn = 1.0/(1 + temp);
                        pp = temp*pn; 
                    }

                    h(0,0) += f[i]*f[i]*pp*pn;
                    const double temp2 = f[i]*pp*pn;
                    h(0,1) += temp2;
                    h(1,0) += temp2;
                    h(1,1) += pp*pn;
                }

                return h;
            }

            const vect_type& f;
            const vect_type& t;
        };
    }

// ----------------------------------------------------------------------------------------

    template <
        typename trainer_type,
        typename in_sample_vector_type,
        typename in_scalar_vector_type
        >
    const probabilistic_decision_function<typename trainer_type::kernel_type> train_probabilistic_decision_function_impl (
        const trainer_type& trainer,
        const in_sample_vector_type& x,
        const in_scalar_vector_type& y,
        const long folds
    )
    {
        typedef typename trainer_type::sample_type sample_type;
        typedef typename trainer_type::scalar_type scalar_type;
        typedef typename trainer_type::mem_manager_type mem_manager_type;
        typedef typename trainer_type::kernel_type K;


        /*
            This function fits a sigmoid function to the output of the 
            svm trained by svm_nu_trainer.  The technique used is the one
            described in the papers:
                
                Probabilistic Outputs for Support Vector Machines and
                Comparisons to Regularized Likelihood Methods by 
                John C. Platt.  March 26, 1999

                A Note on Platt's Probabilistic Outputs for Support Vector Machines
                by Hsuan-Tien Lin, Chih-Jen Lin, and Ruby C. Weng
        */

        // make sure requires clause is not broken
        DLIB_ASSERT(is_binary_classification_problem(x,y) == true &&
                    1 < folds && folds <= x.nr(),
            "\tprobabilistic_decision_function train_probabilistic_decision_function()"
            << "\n\t invalid inputs were given to this function"
            << "\n\t x.nr(): " << x.nr() 
            << "\n\t y.nr(): " << y.nr() 
            << "\n\t x.nc(): " << x.nc() 
            << "\n\t y.nc(): " << y.nc() 
            << "\n\t folds:  " << folds 
            << "\n\t is_binary_classification_problem(x,y): " << ((is_binary_classification_problem(x,y))? "true":"false")
            );

        // count the number of positive and negative examples
        const long num_pos = sum(y > 0);
        const long num_neg = sum(y < 0);

        // figure out how many positive and negative examples we will have in each fold
        const long num_pos_test_samples = num_pos/folds; 
        const long num_pos_train_samples = num_pos - num_pos_test_samples; 
        const long num_neg_test_samples = num_neg/folds; 
        const long num_neg_train_samples = num_neg - num_neg_test_samples; 

        decision_function<K> d;
        typename decision_function<K>::sample_vector_type x_test, x_train;
        typename decision_function<K>::scalar_vector_type y_test, y_train;
        x_test.set_size (num_pos_test_samples  + num_neg_test_samples);
        y_test.set_size (num_pos_test_samples  + num_neg_test_samples);
        x_train.set_size(num_pos_train_samples + num_neg_train_samples);
        y_train.set_size(num_pos_train_samples + num_neg_train_samples);

        typedef std_allocator<scalar_type, mem_manager_type> alloc_scalar_type_vector;
        typedef std::vector<scalar_type, alloc_scalar_type_vector > dvector;

        dvector out;
        dvector target;

        long pos_idx = 0;
        long neg_idx = 0;

        const scalar_type prior0 = num_pos_test_samples*folds; 
        const scalar_type prior1 = num_neg_test_samples*folds; 
        const scalar_type hi_target = (prior1+1)/(prior1+2);
        const scalar_type lo_target = 1.0/(prior0+2);

        for (long i = 0; i < folds; ++i)
        {
            long cur = 0;

            // load up our positive test samples
            while (cur < num_pos_test_samples)
            {
                if (y(pos_idx) == +1.0)
                {
                    x_test(cur) = x(pos_idx);
                    y_test(cur) = +1.0;
                    ++cur;
                }
                pos_idx = (pos_idx+1)%x.nr();
            }

            // load up our negative test samples
            while (cur < x_test.nr())
            {
                if (y(neg_idx) == -1.0)
                {
                    x_test(cur) = x(neg_idx);
                    y_test(cur) = -1.0;
                    ++cur;
                }
                neg_idx = (neg_idx+1)%x.nr();
            }

            // load the training data from the data following whatever we loaded
            // as the testing data
            long train_pos_idx = pos_idx;
            long train_neg_idx = neg_idx;
            cur = 0;

            // load up our positive train samples
            while (cur < num_pos_train_samples)
            {
                if (y(train_pos_idx) == +1.0)
                {
                    x_train(cur) = x(train_pos_idx);
                    y_train(cur) = +1.0;
                    ++cur;
                }
                train_pos_idx = (train_pos_idx+1)%x.nr();
            }

            // load up our negative train samples
            while (cur < x_train.nr())
            {
                if (y(train_neg_idx) == -1.0)
                {
                    x_train(cur) = x(train_neg_idx);
                    y_train(cur) = -1.0;
                    ++cur;
                }
                train_neg_idx = (train_neg_idx+1)%x.nr();
            }

            // do the training
            d = trainer.train (x_train,y_train);

            // now test this fold 
            for (long i = 0; i < x_test.nr(); ++i)
            {
                out.push_back(d(x_test(i)));
                // if this was a positive example
                if (y_test(i) == +1.0)
                {
                    target.push_back(hi_target);
                }
                else if (y_test(i) == -1.0)
                {
                    target.push_back(lo_target);
                }
                else
                {
                    throw dlib::error("invalid input labels to the train_probabilistic_decision_function() function");
                }
            }

        } // for (long i = 0; i < folds; ++i)

        // Now find the maximum likelihood parameters of the sigmoid.  

        prob_impl::objective<dvector> obj(out, target);
        prob_impl::der<dvector> obj_der(out, target);
        prob_impl::hessian<dvector> obj_hessian(out, target);

        matrix<double,2,1> val;
        val = 0;
        find_min(newton_search_strategy(obj_hessian),
                 gradient_norm_stop_strategy(),
                 obj,
                 obj_der,
                 val,
                 0);

        const double A = val(0);
        const double B = val(1);

        return probabilistic_decision_function<K>( A, B, trainer.train(x,y) );
    }

    template <
        typename trainer_type,
        typename in_sample_vector_type,
        typename in_scalar_vector_type
        >
    const probabilistic_decision_function<typename trainer_type::kernel_type> train_probabilistic_decision_function (
        const trainer_type& trainer,
        const in_sample_vector_type& x,
        const in_scalar_vector_type& y,
        const long folds
    )
    {
        return train_probabilistic_decision_function_impl(trainer,
                                                          vector_to_matrix(x),
                                                          vector_to_matrix(y),
                                                          folds);
    }

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    template <
        typename T,
        typename U,
        typename rand_type 
        >
    typename enable_if<is_matrix<T>,void>::type randomize_samples (
        T& t,
        U& u,
        rand_type& r
    )
    {
        // make sure requires clause is not broken
        DLIB_ASSERT(is_vector(t) && is_vector(u) && u.size() == t.size(),
            "\t randomize_samples(t,u)"
            << "\n\t invalid inputs were given to this function"
            << "\n\t t.size(): " << t.size()
            << "\n\t u.size(): " << u.size()
            << "\n\t is_vector(t): " << (is_vector(t)? "true" : "false")
            << "\n\t is_vector(u): " << (is_vector(u)? "true" : "false")
            );

        long n = t.size()-1;
        while (n > 0)
        {
            // put a random integer into idx
            unsigned long idx = r.get_random_32bit_number();

            // make idx be less than n
            idx %= n;

            // swap our randomly selected index into the n position
            exchange(t(idx), t(n));
            exchange(u(idx), u(n));

            --n;
        }
    }

// ----------------------------------------------------------------------------------------

    template <
        typename T,
        typename U,
        typename rand_type
        >
    typename disable_if<is_matrix<T>,void>::type randomize_samples (
        T& t,
        U& u,
        rand_type& r
    )
    {
        // make sure requires clause is not broken
        DLIB_ASSERT(u.size() == t.size(),
            "\t randomize_samples(t,u)"
            << "\n\t invalid inputs were given to this function"
            << "\n\t t.size(): " << t.size()
            << "\n\t u.size(): " << u.size()
            );

        long n = t.size()-1;
        while (n > 0)
        {
            // put a random integer into idx
            unsigned long idx = r.get_random_32bit_number();

            // make idx be less than n
            idx %= n;

            // swap our randomly selected index into the n position
            exchange(t[idx], t[n]);
            exchange(u[idx], u[n]);

            --n;
        }
    }

// ----------------------------------------------------------------------------------------

    template <
        typename T,
        typename U
        >
    typename disable_if<is_rand<U>,void>::type randomize_samples (
        T& t,
        U& u
    )
    {
        rand::kernel_1a r;
        randomize_samples(t,u,r);
    }

// ----------------------------------------------------------------------------------------

    template <
        typename T,
        typename rand_type
        >
    typename enable_if_c<is_matrix<T>::value && is_rand<rand_type>::value,void>::type randomize_samples (
        T& t,
        rand_type& r
    )
    {
        // make sure requires clause is not broken
        DLIB_ASSERT(is_vector(t),
            "\t randomize_samples(t)"
            << "\n\t invalid inputs were given to this function"
            << "\n\t is_vector(t): " << (is_vector(t)? "true" : "false")
            );

        long n = t.size()-1;
        while (n > 0)
        {
            // put a random integer into idx
            unsigned long idx = r.get_random_32bit_number();

            // make idx be less than n
            idx %= n;

            // swap our randomly selected index into the n position
            exchange(t(idx), t(n));

            --n;
        }
    }

// ----------------------------------------------------------------------------------------

    template <
        typename T,
        typename rand_type
        >
    typename disable_if_c<(is_matrix<T>::value==true)||(is_rand<rand_type>::value==false),void>::type randomize_samples (
        T& t,
        rand_type& r
    )
    {
        long n = t.size()-1;
        while (n > 0)
        {
            // put a random integer into idx
            unsigned long idx = r.get_random_32bit_number();

            // make idx be less than n
            idx %= n;

            // swap our randomly selected index into the n position
            exchange(t[idx], t[n]);

            --n;
        }
    }

// ----------------------------------------------------------------------------------------

    template <
        typename T
        >
    void randomize_samples (
        T& t
    )
    {
        rand::kernel_1a r;
        randomize_samples(t,r);
    }

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    template <
        typename K 
        >
    class svm_nu_trainer
    {
    public:
        typedef K kernel_type;
        typedef typename kernel_type::scalar_type scalar_type;
        typedef typename kernel_type::sample_type sample_type;
        typedef typename kernel_type::mem_manager_type mem_manager_type;
        typedef decision_function<kernel_type> trained_function_type;

        svm_nu_trainer (
        ) :
            nu(0.1),
            cache_size(200),
            eps(0.001)
        {
        }

        svm_nu_trainer (
            const kernel_type& kernel_, 
            const scalar_type& nu_
        ) :
            kernel_function(kernel_),
            nu(nu_),
            cache_size(200),
            eps(0.001)
        {
            // make sure requires clause is not broken
            DLIB_ASSERT(0 < nu && nu <= 1,
                "\tsvm_nu_trainer::svm_nu_trainer(kernel,nu)"
                << "\n\t invalid inputs were given to this function"
                << "\n\t nu: " << nu 
                );
        }

        void set_cache_size (
            long cache_size_
        )
        {
            // make sure requires clause is not broken
            DLIB_ASSERT(cache_size_ > 0,
                "\tvoid svm_nu_trainer::set_cache_size(cache_size_)"
                << "\n\t invalid inputs were given to this function"
                << "\n\t cache_size: " << cache_size_ 
                );
            cache_size = cache_size_;
        }

        long get_cache_size (
        ) const
        {
            return cache_size;
        }

        void set_epsilon (
            scalar_type eps_
        )
        {
            // make sure requires clause is not broken
            DLIB_ASSERT(eps_ > 0,
                "\tvoid svm_nu_trainer::set_epsilon(eps_)"
                << "\n\t invalid inputs were given to this function"
                << "\n\t eps: " << eps_ 
                );
            eps = eps_;
        }

        const scalar_type get_epsilon (
        ) const
        { 
            return eps;
        }

        void set_kernel (
            const kernel_type& k
        )
        {
            kernel_function = k;
        }

        const kernel_type& get_kernel (
        ) const
        {
            return kernel_function;
        }

        void set_nu (
            scalar_type nu_
        )
        {
            // make sure requires clause is not broken
            DLIB_ASSERT(0 < nu_ && nu_ <= 1,
                "\tvoid svm_nu_trainer::set_nu(nu_)"
                << "\n\t invalid inputs were given to this function"
                << "\n\t nu: " << nu_ 
                );
            nu = nu_;
        }

        const scalar_type get_nu (
        ) const
        {
            return nu;
        }

        template <
            typename in_sample_vector_type,
            typename in_scalar_vector_type
            >
        const decision_function<kernel_type> train (
            const in_sample_vector_type& x,
            const in_scalar_vector_type& y
        ) const
        {
            return do_train(vector_to_matrix(x), vector_to_matrix(y));
        }

        void swap (
            svm_nu_trainer& item
        )
        {
            exchange(kernel_function, item.kernel_function);
            exchange(nu,              item.nu);
            exchange(cache_size,      item.cache_size);
            exchange(eps,             item.eps);
        }

    private:

    // ------------------------------------------------------------------------------------

        template <
            typename in_sample_vector_type,
            typename in_scalar_vector_type
            >
        const decision_function<kernel_type> do_train (
            const in_sample_vector_type& x,
            const in_scalar_vector_type& y
        ) const
        {
            typedef typename K::scalar_type scalar_type;
            typedef typename decision_function<K>::sample_vector_type sample_vector_type;
            typedef typename decision_function<K>::scalar_vector_type scalar_vector_type;

            // make sure requires clause is not broken
            DLIB_ASSERT(is_binary_classification_problem(x,y) == true,
                "\tdecision_function svm_nu_trainer::train(x,y)"
                << "\n\t invalid inputs were given to this function"
                << "\n\t x.nr(): " << x.nr() 
                << "\n\t y.nr(): " << y.nr() 
                << "\n\t x.nc(): " << x.nc() 
                << "\n\t y.nc(): " << y.nc() 
                << "\n\t is_binary_classification_problem(x,y): " << ((is_binary_classification_problem(x,y))? "true":"false")
                );


            const scalar_type tau = 1e-12;
            scalar_vector_type df; // delta f(alpha)
            scalar_vector_type alpha;

            kernel_matrix_cache<K, in_sample_vector_type, in_scalar_vector_type> Q(x,y,kernel_function,cache_size);
            typedef typename kernel_matrix_cache<K, in_sample_vector_type, in_scalar_vector_type>::scalar_type cache_type;

            alpha.set_size(x.nr());
            df.set_size(x.nr());

            // now initialize alpha
            set_initial_alpha(y, nu, alpha);


            set_all_elements(df, 0);
            // initialize df.  Compute df = Q*alpha
            for (long r = 0; r < df.nr(); ++r)
            {
                if (alpha(r) != 0)
                {
                    const cache_type* Q_r = Q.col(r);
                    for (long c = 0; c < alpha.nr(); ++c)
                    {
                        df(c) += alpha(r)*Q_r[c];
                    }
                }
            }

            // now perform the actual optimization of alpha
            long i, j;
            while (find_working_group(y,alpha,Q,df,tau,eps,i,j))
            {
                const scalar_type old_alpha_i = alpha(i);
                const scalar_type old_alpha_j = alpha(j);

                optimize_working_pair(y,alpha,Q,df,tau,i,j);

                // update the df vector now that we have modified alpha(i) and alpha(j)
                scalar_type delta_alpha_i = alpha(i) - old_alpha_i;
                scalar_type delta_alpha_j = alpha(j) - old_alpha_j;

                const cache_type* Q_i = Q.col(i);
                const cache_type* Q_j = Q.col(j);

                for(long k = 0; k < df.nr(); ++k)
                    df(k) += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;

            }

            scalar_type rho, b;
            calculate_rho_and_b(y,alpha,df,rho,b);
            alpha = pointwise_multiply(alpha,y)/rho;

            // count the number of support vectors
            long sv_count = 0;
            for (long i = 0; i < alpha.nr(); ++i)
            {
                if (alpha(i) != 0)
                    ++sv_count;
            }

            scalar_vector_type sv_alpha;
            sample_vector_type support_vectors;

            // size these column vectors so that they have an entry for each support vector
            sv_alpha.set_size(sv_count);
            support_vectors.set_size(sv_count);

            // load the support vectors and their alpha values into these new column matrices
            long idx = 0;
            for (long i = 0; i < alpha.nr(); ++i)
            {
                if (alpha(i) != 0)
                {
                    sv_alpha(idx) = alpha(i);
                    support_vectors(idx) = x(i);
                    ++idx;
                }
            }

            // now return the decision function
            return decision_function<K> (sv_alpha, b, kernel_function, support_vectors);
        }

    // ------------------------------------------------------------------------------------

        template <
            typename scalar_type,
            typename scalar_vector_type,
            typename scalar_vector_type2
            >
        inline void set_initial_alpha (
            const scalar_vector_type& y,
            const scalar_type nu,
            scalar_vector_type2& alpha
        ) const
        {
            set_all_elements(alpha,0);
            const scalar_type l = y.nr();
            scalar_type temp = nu*l/2;
            long num = (long)std::floor(temp);
            long num_total = (long)std::ceil(temp);

            bool has_slack = false;
            int count = 0;
            for (int i = 0; i < alpha.nr(); ++i)
            {
                if (y(i) == 1)
                {
                    if (count < num)
                    {
                        ++count;
                        alpha(i) = 1;
                    }
                    else 
                    {
                        has_slack = true;
                        if (num_total > num)
                        {
                            ++count;
                            alpha(i) = temp - std::floor(temp);
                        }
                        break;
                    }
                }
            }

            if (count != num_total || has_slack == false)
            {
                std::ostringstream sout;
                sout << "Invalid nu of " << nu << ".  It is required that: 0 < nu < " << 2*(scalar_type)count/y.nr();
                throw invalid_svm_nu_error(sout.str(),nu);
            }

            has_slack = false;
            count = 0;
            for (int i = 0; i < alpha.nr(); ++i)
            {
                if (y(i) == -1)
                {
                    if (count < num)
                    {
                        ++count;
                        alpha(i) = 1;
                    }
                    else 
                    {
                        has_slack = true;
                        if (num_total > num)
                        {
                            ++count;
                            alpha(i) = temp - std::floor(temp);
                        }
                        break;
                    }
                }
            }

            if (count != num_total || has_slack == false)
            {
                std::ostringstream sout;
                sout << "Invalid nu of " << nu << ".  It is required that: 0 < nu < " << 2*(scalar_type)count/y.nr();
                throw invalid_svm_nu_error(sout.str(),nu);
            }
        }

    // ------------------------------------------------------------------------------------

        template <
            typename sample_vector_type,
            typename scalar_vector_type,
            typename scalar_vector_type2,
            typename scalar_type
            >
        inline bool find_working_group (
            const scalar_vector_type2& y,
            const scalar_vector_type& alpha,
            const kernel_matrix_cache<K,sample_vector_type, scalar_vector_type2>& Q,
            const scalar_vector_type& df,
            const scalar_type tau,
            const scalar_type eps,
            long& i_out,
            long& j_out
        ) const
        {
            using namespace std;
            long ip = -1;
            long jp = -1;
            long in = -1;
            long jn = -1;


            typedef typename kernel_matrix_cache<K, sample_vector_type, scalar_vector_type2>::scalar_type cache_type;

            scalar_type ip_val = -numeric_limits<scalar_type>::infinity();
            scalar_type jp_val = numeric_limits<scalar_type>::infinity();
            scalar_type in_val = -numeric_limits<scalar_type>::infinity();
            scalar_type jn_val = numeric_limits<scalar_type>::infinity();

            // loop over the alphas and find the maximum ip and in indices.
            for (long i = 0; i < alpha.nr(); ++i)
            {
                if (y(i) == 1)
                {
                    if (alpha(i) < 1.0)
                    {
                        if (-df(i) > ip_val)
                        {
                            ip_val = -df(i);
                            ip = i;
                        }
                    }
                }
                else
                {
                    if (alpha(i) > 0.0)
                    {
                        if (df(i) > in_val)
                        {
                            in_val = df(i);
                            in = i;
                        }
                    }
                }
            }

            scalar_type Mp = numeric_limits<scalar_type>::infinity();
            scalar_type Mn = numeric_limits<scalar_type>::infinity();
            scalar_type bp = -numeric_limits<scalar_type>::infinity();
            scalar_type bn = -numeric_limits<scalar_type>::infinity();

            // As a speed hack, pull out pointers to the columns of the
            // kernel matrix we will be using below rather than accessing
            // them through the Q(r,c) syntax.
            const cache_type* Q_ip = 0;
            const cache_type* Q_in = 0;
            const cache_type* Q_diag = Q.diag();
            if (ip != -1)
                Q_ip = Q.col(ip);
            if (in != -1)
                Q_in = Q.col(in);


            // now we need to find the minimum jp and jn indices
            for (long j = 0; j < alpha.nr(); ++j)
            {
                if (y(j) == 1)
                {
                    if (alpha(j) > 0.0)
                    {
                        scalar_type b = ip_val + df(j);
                        if (-df(j) < Mp)
                            Mp = -df(j);

                        if (b > 0)
                        {
                            bp = b;
                            scalar_type a = Q_ip[ip] + Q_diag[j] - 2*Q_ip[j]; 
                            if (a <= 0)
                                a = tau;
                            scalar_type temp = -b*b/a;
                            if (temp < jp_val)
                            {
                                jp_val = temp;
                                jp = j;
                            }
                        }
                    }
                }
                else
                {
                    if (alpha(j) < 1.0)
                    {
                        scalar_type b = in_val - df(j);
                        if (df(j) < Mn)
                            Mn = df(j);

                        if (b > 0)
                        {
                            bn = b;
                            scalar_type a = Q_in[in] + Q_diag[j] - 2*Q_in[j]; 
                            if (a <= 0)
                                a = tau;
                            scalar_type temp = -b*b/a;
                            if (temp < jn_val)
                            {
                                jn_val = temp;
                                jn = j;
                            }
                        }
                    }
                }
            }

            // if we are at the optimal point then return false so the caller knows
            // to stop optimizing
            if (std::max(ip_val - Mp, in_val - Mn) < eps)
                return false;

            if (jp_val < jn_val)
            {
                i_out = ip;
                j_out = jp;
            }
            else
            {
                i_out = in;
                j_out = jn;
            }

            if (j_out >= 0 && i_out >= 0)
                return true;
            else
                return false;
        }

    // ------------------------------------------------------------------------------------

        template <
            typename scalar_vector_type,
            typename scalar_vector_type2,
            typename scalar_type
            >
        void calculate_rho_and_b(
            const scalar_vector_type2& y,
            const scalar_vector_type& alpha,
            const scalar_vector_type& df,
            scalar_type& rho, 
            scalar_type& b
        ) const
        {
            using namespace std;
            long num_p_free = 0;
            long num_n_free = 0;
            scalar_type sum_p_free = 0;
            scalar_type sum_n_free = 0;

            scalar_type upper_bound_p = -numeric_limits<scalar_type>::infinity();
            scalar_type upper_bound_n = -numeric_limits<scalar_type>::infinity();
            scalar_type lower_bound_p = numeric_limits<scalar_type>::infinity();
            scalar_type lower_bound_n = numeric_limits<scalar_type>::infinity();

            for(long i = 0; i < alpha.nr(); ++i)
            {
                if(y(i) == 1)
                {
                    if(alpha(i) == 1)
                    {
                        if (df(i) > upper_bound_p)
                            upper_bound_p = df(i);
                    }
                    else if(alpha(i) == 0)
                    {
                        if (df(i) < lower_bound_p)
                            lower_bound_p = df(i);
                    }
                    else
                    {
                        ++num_p_free;
                        sum_p_free += df(i);
                    }
                }
                else
                {
                    if(alpha(i) == 1)
                    {
                        if (df(i) > upper_bound_n)
                            upper_bound_n = df(i);
                    }
                    else if(alpha(i) == 0)
                    {
                        if (df(i) < lower_bound_n)
                            lower_bound_n = df(i);
                    }
                    else
                    {
                        ++num_n_free;
                        sum_n_free += df(i);
                    }
                }
            }

            scalar_type r1,r2;
            if(num_p_free > 0)
                r1 = sum_p_free/num_p_free;
            else
                r1 = (upper_bound_p+lower_bound_p)/2;

            if(num_n_free > 0)
                r2 = sum_n_free/num_n_free;
            else
                r2 = (upper_bound_n+lower_bound_n)/2;

            rho = (r1+r2)/2;
            b = (r1-r2)/2/rho;
        }

    // ------------------------------------------------------------------------------------

        template <
            typename sample_vector_type,
            typename scalar_vector_type,
            typename scalar_vector_type2,
            typename scalar_type
            >
        inline void optimize_working_pair (
            const scalar_vector_type2& ,
            scalar_vector_type& alpha,
            const kernel_matrix_cache<K, sample_vector_type, scalar_vector_type2>& Q,
            const scalar_vector_type& df,
            const scalar_type tau,
            const long i,
            const long j
        ) const
        {
            scalar_type quad_coef = Q(i,i)+Q(j,j)-2*Q(j,i);
            if (quad_coef <= 0)
                quad_coef = tau;
            scalar_type delta = (df(i)-df(j))/quad_coef;
            scalar_type sum = alpha(i) + alpha(j);
            alpha(i) -= delta;
            alpha(j) += delta;

            if(sum > 1)
            {
                if(alpha(i) > 1)
                {
                    alpha(i) = 1;
                    alpha(j) = sum - 1;
                }
                else if(alpha(j) > 1)
                {
                    alpha(j) = 1;
                    alpha(i) = sum - 1;
                }
            }
            else
            {
                if(alpha(j) < 0)
                {
                    alpha(j) = 0;
                    alpha(i) = sum;
                }
                else if(alpha(i) < 0)
                {
                    alpha(i) = 0;
                    alpha(j) = sum;
                }
            }
        }

    // ------------------------------------------------------------------------------------

        kernel_type kernel_function;
        scalar_type nu;
        long cache_size;
        scalar_type eps;
    }; // end of class svm_nu_trainer

// ----------------------------------------------------------------------------------------

    template <typename K>
    void swap (
        svm_nu_trainer<K>& a,
        svm_nu_trainer<K>& b
    ) { a.swap(b); }

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
}

#endif // DLIB_SVm_