CMonteCarlo.h

Go to the documentation of this file.

/* +------------------------------------------------------------------------+
   |                     Mobile Robot Programming Toolkit (MRPT)            |
   |                          http://www.mrpt.org/                          |
   |                                                                        |
   | Copyright (c) 2005-2017, 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/utils/CTicTac.h>

namespace mrpt
{
namespace 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_base_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::utils::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;
        }
};
}
}  // End of namespaces
#endif