CMonteCarlo.h

Go to the documentation of this file.

/* +------------------------------------------------------------------------+
   |                     Mobile Robot Programming Toolkit (MRPT)            |
   |                          http://www.mrpt.org/                          |
   |                                                                        |
   | Copyright (c) 2005-2018, Individual contributors, see AUTHORS file     |
   | See: http://www.mrpt.org/Authors - All rights reserved.                |
   | Released under BSD License. See details in http://www.mrpt.org/License |
   +------------------------------------------------------------------------+ */
#ifndef _MRPT_MONTE_CARLO_H_
#define _MRPT_MONTE_CARLO_H_

#include <map>
#include <vector>
#include <numeric>
#include <algorithm>
#include <stdexcept>
#include <functional>
#include <mrpt/random.h>
#include <mrpt/system/CTicTac.h>

namespace mrpt::math
{
/** Montecarlo simulation for experiments in 1D.
     Template arguments are:
        - T: base type, i.e., if an experiment needs to generate random points,
 then T may be a TPoint3D, and so on.
        - NUM: the numeric type used to represent the error. Usually, double.
        - OTHER: an intermediate type, used especially when testing inverse
 functions. Leave as int or double if you don't use it.

    HOW TO USE THIS CLASS:
        - Create an instance of the class.
        - Refill the "valueGenerator" member with an appropriate function.
        - If your experiment calculates the error directly from the base value,
 then refill the "errorFun1" member.
        - Otherwise, if your experiment involves the calculation of some value
 whom with the experimental function is compared, refill "intermediateFun" and
 "errorFun2".
        - Refill only on of the alternatives.
 * \ingroup mrpt_math_grp
*/
template <typename T, typename NUM, typename OTHER>
class CMonteCarlo
{
   private:
    mrpt::random::CRandomGenerator gen;
    class CStatisticalAnalyzer
    {
       private:
        Eigen::Matrix<NUM, Eigen::Dynamic, 1> data;

       public:
        template <typename VEC>
        inline CStatisticalAnalyzer(const VEC& v1) : data(v1.begin(), v1.end())
        {
        }
        template <typename VEC>
        inline void setData(const VEC& v1)
        {
            data.assign(v1.begin(), v1.end());
        }
        template <typename VEC>
        inline void getData(VEC& v1) const
        {
            v1.assign(data.begin(), data.end());
        }
        template <typename VEC1, typename VEC2>
        inline void getDistribution(
            VEC1& vx, VEC2& vy, const NUM width = 1.0) const
        {
            std::vector<double> vvx, vvy;
            getDistribution(vvx, vvy, width);
            vx.assign(vvx.begin(), vvx.end());
            vy.assign(vvy.begin(), vvy.end());
        }
        // Function overload, not specialization (GCC complains otherwise):
        inline void getDistribution(
            std::vector<double>& vx, std::vector<double>& vy,
            const NUM width = 1.0) const
        {
            CHistogram hist(
                CHistogram::createWithFixedWidth(
                    0, *max_element(data.begin(), data.end()), width));
            hist.add(data);
            hist.getHistogram(vx, vy);
        }
    };

   public:
    // TODO: use templates for function types.
    // Random generator.
    T (*valueGenerator)(mrpt::random::CRandomGenerator&);
    // Computes automatically the error (without an intermediate type)
    NUM (*errorFun1)(const T&);

    OTHER (*intermediateFun)(const T&);
    NUM (*errorFun2)(const T&, const OTHER&);
    inline CMonteCarlo()
        : gen(),
          valueGenerator(nullptr),
          errorFun1(nullptr),
          intermediateFun(nullptr),
          errorFun2(nullptr)
    {
    }
    NUM doExperiment(size_t N, double& time, bool showInWindow = false)
    {
        if (!valueGenerator)
            throw std::logic_error("Value generator function is not set.");
        std::vector<T> baseData(N);
        std::vector<NUM> errorData(N);
        mrpt::system::CTicTac meter;
        for (size_t i = 0; i < N; ++i) baseData[i] = valueGenerator(gen);
        if (errorFun1)
        {
            meter.Tic();
            std::transform(
                baseData.begin(), baseData.end(), errorData.begin(), errorFun1);
            time = meter.Tac();
        }
        else
        {
            if (!intermediateFun || !errorFun2)
                throw std::logic_error(
                    "Experiment-related functions are not set.");
            std::vector<OTHER> intermediate(N);
            transform(
                baseData.begin(), baseData.end(), intermediate.begin(),
                intermediateFun);
            meter.Tic();
            for (size_t i = 0; i < N; ++i)
                errorData[i] = errorFun2(baseData[i], intermediate[i]);
            time = meter.Tac();
        }
        NUM res = accumulate(errorData.begin(), errorData.end(), NUM(0)) /
                  errorData.size();
#if 0
        if (showInWindow)
        {
            CStatisticalAnalyzer st(errorData);
            mrpt::gui::CDisplayWindowPlots wnd("Error results from Monte Carlo simulation");
            std::vector<NUM> errorX,errorY;
            st.getDistribution(errorX,errorY,0.1);
            wnd.plot(errorX,errorY,"b-","Plot1");
            NUM maxVal=*std::max_element(errorY.begin(),errorY.end());
            const std::vector<NUM> dx{res, res}, dy{.0, maxVal};
            wnd.plot(dx,dy,"r-","Plot2");
            while (wnd.isOpen()) {};
        }
#endif
        return res;
    }
};
}
#endif