superviseddescent.hpp

/*
 * superviseddescent: A C++11 implementation of the supervised descent
 *                    optimisation method
 * File: superviseddescent/superviseddescent.hpp
 *
 * Copyright 2014, 2015 Patrik Huber
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once

#ifndef SUPERVISEDDESCENT_HPP_
#define SUPERVISEDDESCENT_HPP_

#include "superviseddescent/utils/ThreadPool.h"

#include "cereal/cereal.hpp"
#include "cereal/types/vector.hpp"

#include "opencv2/core/core.hpp"

namespace superviseddescent {

inline void no_eval(const cv::Mat& current_predictions)
{
};

class NoNormalisation
{
public:
    inline cv::Mat operator()(cv::Mat params) {
        return cv::Mat::ones(1, params.cols, params.type());
    };
};

template<class RegressorType, class NormalisationStrategy = NoNormalisation>
class SupervisedDescentOptimiser
{
public:
    SupervisedDescentOptimiser() = default;

    SupervisedDescentOptimiser(std::vector<RegressorType> regressors, NormalisationStrategy normalisation = NoNormalisation()) : regressors(std::move(regressors)), normalisation_strategy(std::move(normalisation))
    {
    };

    template<class ProjectionFunction>
    void train(cv::Mat parameters, cv::Mat initialisations, cv::Mat templates, ProjectionFunction projection)
    {
        return train(parameters, initialisations, templates, projection, no_eval);
    };

    template<class ProjectionFunction, class OnTrainingEpochCallback>
    void train(cv::Mat parameters, cv::Mat initialisations, cv::Mat templates, ProjectionFunction projection, OnTrainingEpochCallback on_training_epoch_callback)
    {
        using cv::Mat;
        Mat current_x = initialisations;
        for (size_t regressor_level = 0; regressor_level < regressors.size(); ++regressor_level) {
            // 1) Project current parameters x to feature space:
            // Enqueue all tasks in a thread pool:
            auto concurent_threads_supported = std::thread::hardware_concurrency();
            if (concurent_threads_supported == 0) {
                concurent_threads_supported = 4;
            }
            utils::ThreadPool thread_pool(concurent_threads_supported);
            std::vector<std::future<typename std::result_of<ProjectionFunction(Mat, size_t, int)>::type>> results; // will be float or Mat. I might remove float for the sake of code clarity, as it's only useful for very simple examples.
            results.reserve(current_x.rows);
            for (int sample_index = 0; sample_index < current_x.rows; ++sample_index) {
                results.emplace_back(
                    thread_pool.enqueue(projection, current_x.row(sample_index), regressor_level, sample_index)
                );
            }
            // Gather the results from all threads and store the features:
            Mat features;
            for (auto&& result : results) {
                features.push_back(result.get());
            }
            // Set the observed values, depending on if a template y is used:
            Mat observed_values;
            if (templates.empty()) { // unknown template training case
                observed_values = features;
            }
            else { // known template
                observed_values = features - templates;
            }
            //Mat b = currentX - parameters; // currentX - x;
            Mat b; // Todo: reserve() for speedup. Also below with x_k.
            // Apply the normalisation strategy to each sample in b:
            for (int sample_index = 0; sample_index < current_x.rows; ++sample_index) {
                cv::Mat update_step = current_x.row(sample_index) - parameters.row(sample_index);
                update_step = update_step.mul(normalisation_strategy(current_x.row(sample_index)));
                b.push_back(update_step);
            }
            // 2) Learn using that data:
            regressors[regressor_level].learn(observed_values, b);
            // 3) Apply the learned regressor and use the predictions to learn the next regressor in next loop iteration:
            Mat x_k; // x_k = currentX - R * (h(currentX) - y):
            for (int sample_index = 0; sample_index < current_x.rows; ++sample_index) {
                // No need to re-extract the features, we already did so in step 1)
                cv::Mat update_step = regressors[regressor_level].predict(observed_values.row(sample_index));
                update_step = update_step.mul(1 / normalisation_strategy(current_x.row(sample_index))); // Need to multiply the regressor prediction with the IED of the current prediction
                x_k.push_back(Mat(current_x.row(sample_index) - update_step));
            }
            current_x = x_k;
            on_training_epoch_callback(current_x);
        }
    };

    template<class ProjectionFunction>
    cv::Mat test(cv::Mat initialisations, cv::Mat templates, ProjectionFunction projection)
    {
        return test(initialisations, templates, projection, no_eval);
    };

    template<class ProjectionFunction, class OnRegressorIterationCallback>
    cv::Mat test(cv::Mat initialisations, cv::Mat templates, ProjectionFunction projection, OnRegressorIterationCallback on_regressor_iteration_callback)
    {
        using cv::Mat;
        Mat current_x = initialisations;
        for (size_t regressor_level = 0; regressor_level < regressors.size(); ++regressor_level) {
            // Enqueue all tasks in a thread pool:
            auto concurent_threads_supported = std::thread::hardware_concurrency();
            if (concurent_threads_supported == 0) {
                concurent_threads_supported = 4;
            }
            utils::ThreadPool thread_pool(concurent_threads_supported);
            std::vector<std::future<typename std::result_of<ProjectionFunction(Mat, size_t, int)>::type>> results; // will be float or Mat. I might remove float for the sake of code clarity, as it's only useful for very simple examples.
            results.reserve(current_x.rows);
            for (int sample_index = 0; sample_index < current_x.rows; ++sample_index) {
                results.emplace_back(
                    thread_pool.enqueue(projection, current_x.row(sample_index), regressor_level, sample_index)
                    );
            }
            // Gather the results from all threads and store the features:
            Mat features;
            for (auto&& result : results) {
                features.push_back(result.get());
            }

            Mat observed_values;
            if (templates.empty()) { // unknown template training case
                observed_values = features;
            }
            else { // known template
                observed_values = features - templates;
            }
            Mat x_k;
            // Calculate x_k = currentX - R * (h(currentX) - y):
            for (int sample_index = 0; sample_index < current_x.rows; ++sample_index) {
                cv::Mat update_step = regressors[regressor_level].predict(observed_values.row(sample_index));
                update_step = update_step.mul(1 / normalisation_strategy(current_x.row(sample_index))); // Need to multiply the regressor prediction with the IED of the current prediction
                x_k.push_back(Mat(current_x.row(sample_index) - update_step));
                //x_k.push_back(Mat(currentX.row(sampleIndex) - regressors[regressorLevel].predict(observedValues.row(sampleIndex)))); // we need Mat() because the subtraction yields a (non-persistent) MatExpr
            }
            current_x = x_k;
            on_regressor_iteration_callback(current_x);
        }
        return current_x; // Return the final predictions
    };

    template<class ProjectionFunction>
    cv::Mat predict(cv::Mat initialisations, cv::Mat templates, ProjectionFunction projection)
    {
        using cv::Mat;
        Mat current_x = initialisations;
        for (size_t r = 0; r < regressors.size(); ++r) {
            // calculate x_k = currentX - R * (h(currentX) - y):
            Mat observed_values;
            if (templates.empty()) { // unknown template training case
                observed_values = projection(current_x, r);
            }
            else { // known template
                observed_values = projection(current_x, r) - templates;
            }
            cv::Mat update_step = regressors[r].predict(observed_values);
            update_step = update_step.mul(1 / normalisation_strategy(current_x)); // Need to multiply the regressor prediction with the IED of the current prediction
            Mat x_k = current_x - update_step;
            //Mat x_k = currentX - regressors[r].predict(observedValues);
            current_x = x_k;
        }
        return current_x;
    };

private:
    std::vector<RegressorType> regressors; 
    NormalisationStrategy normalisation_strategy; 

    friend class cereal::access;
    template<class Archive>
    void serialize(Archive& ar)
    {
        ar(regressors, normalisation_strategy);
    };
};

} /* namespace superviseddescent */
#endif /* SUPERVISEDDESCENT_HPP_ */