regressors.hpp

/*
 * superviseddescent: A C++11 implementation of the supervised descent
 *                    optimisation method
 * File: superviseddescent/regressors.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 REGRESSORS_HPP_
#define REGRESSORS_HPP_

#include "cereal/cereal.hpp"
#include "superviseddescent/utils/mat_cerealisation.hpp"

#include "Eigen/Dense"

#include "opencv2/core/core.hpp"

#include <iostream>

namespace superviseddescent {

class Regressor
{
public:
    virtual ~Regressor() {};

    virtual bool learn(cv::Mat data, cv::Mat labels) = 0;

    virtual double test(cv::Mat data, cv::Mat labels) = 0;

    virtual cv::Mat predict(cv::Mat values) = 0;
};

class Regulariser
{
public:
    enum class RegularisationType
    {
        Manual, 
        MatrixNorm, 
    };

    Regulariser(RegularisationType regularisation_type = RegularisationType::Manual, float param = 0.0f, bool regularise_last_row = true) : regularisation_type(regularisation_type), lambda(param), regularise_last_row(regularise_last_row)
    {
    };

    cv::Mat get_matrix(cv::Mat data, int num_training_elements)
    {
        switch (regularisation_type)
        {
        case RegularisationType::Manual:
            // We just take lambda as it was given, no calculation necessary.
            break;
        case RegularisationType::MatrixNorm:
            // The given lambda is the factor we have to multiply the automatic value with:
            lambda = lambda * static_cast<float>(cv::norm(data)) / static_cast<float>(num_training_elements);
            break;
        default:
            break;
        }

        cv::Mat regulariser = cv::Mat::eye(data.rows, data.cols, CV_32FC1) * lambda;

        if (!regularise_last_row) {
            // no lambda for the bias:
            regulariser.at<float>(regulariser.rows - 1, regulariser.cols - 1) = 0.0f;
        }
        return regulariser;
    };

private:
    RegularisationType regularisation_type; 
    float lambda; 
    bool regularise_last_row; 

    friend class cereal::access;
    template<class Archive>
    void serialize(Archive& ar)
    {
        ar(regularisation_type, lambda, regularise_last_row);
    }
};

class PartialPivLUSolver
{
public:
    cv::Mat solve(cv::Mat data, cv::Mat labels, Regulariser regulariser)
    {
        using cv::Mat;
        using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;

        // Map the cv::Mat data and labels to Eigen matrices:
        Eigen::Map<RowMajorMatrixXf> A_Eigen(data.ptr<float>(), data.rows, data.cols);
        Eigen::Map<RowMajorMatrixXf> labels_Eigen(labels.ptr<float>(), labels.rows, labels.cols);

        RowMajorMatrixXf AtA_Eigen = A_Eigen.transpose() * A_Eigen;

        // Note: This is a bit of unnecessary back-and-forth mapping, just for the regularisation:
        Mat AtA_Map(static_cast<int>(AtA_Eigen.rows()), static_cast<int>(AtA_Eigen.cols()), CV_32FC1, AtA_Eigen.data());
        Mat regularisation_matrix = regulariser.get_matrix(AtA_Map, data.rows);
        Eigen::Map<RowMajorMatrixXf> reg_Eigen(regularisation_matrix.ptr<float>(), regularisation_matrix.rows, regularisation_matrix.cols);

        Eigen::DiagonalMatrix<float, Eigen::Dynamic> reg_Eigen_diag(regularisation_matrix.rows);
        Eigen::VectorXf diag_vec(regularisation_matrix.rows);
        for (int i = 0; i < diag_vec.size(); ++i) {
            diag_vec(i) = regularisation_matrix.at<float>(i, i);
        }
        reg_Eigen_diag.diagonal() = diag_vec;
        AtA_Eigen = AtA_Eigen + reg_Eigen_diag.toDenseMatrix();

        // Perform a fast PartialPivLU and use ::solve() (better than inverting):
        Eigen::PartialPivLU<RowMajorMatrixXf> lu_of_AtA(AtA_Eigen);
        RowMajorMatrixXf x_Eigen = lu_of_AtA.solve(A_Eigen.transpose() * labels_Eigen);
        //RowMajorMatrixXf x_Eigen = AtA_Eigen.partialPivLu.solve(A_Eigen.transpose() * labels_Eigen);

        // Map the resulting x back to a cv::Mat by creating a Mat header:
        Mat x(static_cast<int>(x_Eigen.rows()), static_cast<int>(x_Eigen.cols()), CV_32FC1, x_Eigen.data());

        // We have to copy the data because the underlying data is managed by Eigen::Matrix x_Eigen, which will go out of scope after we leave this function:
        return x.clone();
        //return qrOfAtA.isInvertible();
    };
};

class ColPivHouseholderQRSolver
{
public:
    cv::Mat solve(cv::Mat data, cv::Mat labels, Regulariser regulariser)
    {
        using cv::Mat;
        using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
        // Map the cv::Mat data and labels to Eigen matrices:
        Eigen::Map<RowMajorMatrixXf> A_Eigen(data.ptr<float>(), data.rows, data.cols);
        Eigen::Map<RowMajorMatrixXf> labels_Eigen(labels.ptr<float>(), labels.rows, labels.cols);

        RowMajorMatrixXf AtA_Eigen = A_Eigen.transpose() * A_Eigen;

        // Note: This is a bit of unnecessary back-and-forth mapping, just for the regularisation:
        Mat AtA_Map(static_cast<int>(AtA_Eigen.rows()), static_cast<int>(AtA_Eigen.cols()), CV_32FC1, AtA_Eigen.data());
        Mat regularisation_matrix = regulariser.get_matrix(AtA_Map, data.rows);
        Eigen::Map<RowMajorMatrixXf> reg_Eigen(regularisation_matrix.ptr<float>(), regularisation_matrix.rows, regularisation_matrix.cols);

        Eigen::DiagonalMatrix<float, Eigen::Dynamic> reg_Eigen_diag(regularisation_matrix.rows);
        Eigen::VectorXf diag_vec(regularisation_matrix.rows);
        for (int i = 0; i < diag_vec.size(); ++i) {
            diag_vec(i) = regularisation_matrix.at<float>(i, i);
        }
        reg_Eigen_diag.diagonal() = diag_vec;
        AtA_Eigen = AtA_Eigen + reg_Eigen_diag.toDenseMatrix();

        // Perform a ColPivHouseholderQR (faster than FullPivLU) that allows to check for invertibility:
        Eigen::ColPivHouseholderQR<RowMajorMatrixXf> qr_of_AtA(AtA_Eigen);
        auto rankOfAtA = qr_of_AtA.rank();
        if (!qr_of_AtA.isInvertible()) {
            // Eigen may return garbage (their docu is not very specific). Best option is to increase regularisation.
            std::cout << "The regularised AtA is not invertible. We continued learning, but Eigen may return garbage (their docu is not very specific). (The rank is " << std::to_string(rankOfAtA) << ", full rank would be " << std::to_string(AtA_Eigen.rows()) << "). Increase lambda." << std::endl;
        }
        RowMajorMatrixXf AtAInv_Eigen = qr_of_AtA.inverse();

        // x = (AtAReg)^-1 * At * b:
        RowMajorMatrixXf x_Eigen = AtAInv_Eigen * A_Eigen.transpose() * labels_Eigen;

        // Map the resulting x back to a cv::Mat by creating a Mat header:
        Mat x(static_cast<int>(x_Eigen.rows()), static_cast<int>(x_Eigen.cols()), CV_32FC1, x_Eigen.data());

        // We have to copy the data because the underlying data is managed by Eigen::Matrix x_Eigen, which will go out of scope after we leave this function:
        return x.clone();
        //return qrOfAtA.isInvertible();
    };
};

template<class Solver = PartialPivLUSolver>
class LinearRegressor : public Regressor
{

public:
    LinearRegressor(Regulariser regulariser = Regulariser()) : x(), regulariser(regulariser)
    {
    };

    bool learn(cv::Mat data, cv::Mat labels) override
    {
        cv::Mat x = solver.solve(data, labels, regulariser);
        this->x = x;
        return true; // see todo above
    };

    double test(cv::Mat data, cv::Mat labels) override
    {
        cv::Mat predictions;
        for (int i = 0; i < data.rows; ++i) {
            cv::Mat prediction = predict(data.row(i));
            predictions.push_back(prediction);
        }
        return cv::norm(predictions, labels, cv::NORM_L2) / cv::norm(labels, cv::NORM_L2);
    };

    cv::Mat predict(cv::Mat values) override
    {
        cv::Mat prediction = values * x;
        return prediction;
    };
    
    cv::Mat x; 

private:
    Regulariser regulariser; 
    Solver solver; 

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

} /* namespace superviseddescent */
#endif /* REGRESSORS_HPP_ */