estimator.h

Go to the documentation of this file.

 
 
#ifndef CHEBYSHEV_ESTIMATOR
#define CHEBYSHEV_ESTIMATOR
 
 
#include <functional>
#include <cmath>
 
#include "../core/common.h"
#include "../core/random.h"
#include "./prec_structures.h"
 
 
namespace chebyshev {
namespace prec {
 
 
    namespace estimator {
 
 
        template<typename FloatType = real_t>
        inline auto quadrature1D() {
 
            return [](
                EndoFunction<FloatType> funcApprox,
                EndoFunction<FloatType> funcExpected,
                estimate_options<FloatType, FloatType> options) {
 
                if(options.domain.size() != 1) {
                    throw std::runtime_error(
                        "estimator::quadrature1D only works on mono-dimensional domains"
                    );
                }
 
                interval domain = options.domain[0];
 
                FloatType sum = 0;
                FloatType sumSqr = 0;
                FloatType sumAbs = 0;
                FloatType max = 0;
 
                const FloatType length = domain.length();
                const FloatType dx = length / options.iterations;
                FloatType x;
                FloatType coeff;
 
                FloatType diff = std::abs(funcApprox(domain.a) - funcExpected(domain.a));
 
                sum += diff;
                sumSqr += diff * diff;
                sumAbs += std::abs(funcExpected(domain.a));
                max = diff;
 
                for (unsigned int i = 1; i < options.iterations; ++i) {
 
                    x = domain.a + i * dx;
                    diff = std::abs(funcApprox(x) - funcExpected(x));
 
                    if(diff > max)
                        max = diff;
 
                    if(i % 2 == 0)
                        coeff = 2;
                    else
                        coeff = 4;
 
                    sum += coeff * diff;
                    sumSqr += coeff * diff * diff;
                    sumAbs += coeff * funcExpected(x);
                }
 
                diff = std::abs(funcApprox(domain.b) - funcExpected(domain.b));
 
                sum += diff;
                sumSqr += diff * diff;
                sumAbs += std::abs(funcExpected(domain.b));
                
                if(diff > max)
                    max = diff;
 
                estimate_result res {};
                res.absErr = sum;
                res.maxErr = max;
                res.meanErr = (sum * dx / 3.0) / length;
                res.rmsErr = std::sqrt((sumSqr * dx / 3.0) / length);
                res.relErr = std::abs((sum * dx / 3.0) / (sumAbs * dx / 3.0));
                
                return res;
            };
        }
 
 
        template <
            typename IntType = int,
            typename ReturnType = IntType
        >
        inline auto discrete1D() {
 
            // Return a one-dimensional discrete estimator as a lambda function
            return [](
                std::function<IntType(ReturnType)> funcApprox,
                std::function<IntType(ReturnType)> funcExpected,
                estimate_options<IntType, ReturnType> options) {
 
                if(options.domain.size() != 1) {
                    throw std::runtime_error(
                        "estimator::discrete1D only works on mono-dimensional domains"
                    );
                }
 
                IntType extreme1 = IntType(std::ceil(options.domain[0].a));
                IntType extreme2 = IntType(std::floor(options.domain[0].b));
 
                IntType lower = extreme1 < extreme2 ? extreme1 : extreme2;
                IntType upper = extreme1 > extreme2 ? extreme1 : extreme2;
 
                prec_t maxErr = 0;
                prec_t sumDiff = 0;
                prec_t sumSqr = 0;
                prec_t sumAbs = 0;
                uint64_t totalPoints = 0;
 
                for (IntType n = lower; n <= upper; ++n) {
 
                    const ReturnType resExpected = funcExpected(n);
                    const ReturnType resApprox = funcApprox(n);
 
                    const prec_t diff = (prec_t) resExpected > resApprox ?
                        (resExpected - resApprox) : (resApprox - resExpected);
 
                    maxErr = std::max(maxErr, diff);
                    sumDiff += diff;
                    sumSqr += diff * diff;
                    sumAbs += std::abs((prec_t) funcExpected(n));
                    totalPoints++;
                }
 
                estimate_result res {};
                res.absErr = sumAbs;
                res.maxErr = maxErr;
                res.meanErr = totalPoints > 0 ? (sumDiff / totalPoints) : 0;
                res.rmsErr = totalPoints > 0 ? (std::sqrt(sumSqr) / totalPoints) : 0;
                res.relErr = sumDiff / sumAbs;
                return res;
            };
        }
 
 
        template<typename FloatType = real_t>
        inline auto montecarlo1D(std::shared_ptr<random::random_context> rnd_ctx) {
 
            // Return a one-dimensional Monte Carlo estimator as a lambda function
            return [rnd_ctx](
                EndoFunction<FloatType> funcApprox,
                EndoFunction<FloatType> funcExpected,
                estimate_options<FloatType, FloatType> options) {
 
                random::random_source rnd = rnd_ctx->get_rnd();
 
                if(options.domain.size() != 1) {
                    throw std::runtime_error(
                        "estimator::montecarlo1D only works on mono-dimensional domains"
                    );
                }
 
                FloatType sum = 0;
                FloatType sumSqr = 0;
                FloatType sumAbs = 0;
                FloatType max = 0;
                const FloatType length = options.domain[0].length();
 
                for (int i = 0; i < options.iterations; ++i) {
                    
                    FloatType x = rnd.uniform(options.domain[0].a, options.domain[0].b);
                    const FloatType diff = std::abs(funcApprox(x) - funcExpected(x));
 
                    max = std::max(max, diff);
                    sum += diff;
                    sumSqr += diff * diff;
                    sumAbs += std::abs(funcExpected(x));
                }
 
                estimate_result res {};
                res.maxErr = max;
                res.meanErr = sum / options.iterations;
                res.absErr = sum * (length / options.iterations);
                res.rmsErr = std::sqrt(sumSqr * (length / options.iterations));
                res.relErr = sum / sumAbs;
 
                return res;
            };
        }
 
 
        template <
            typename FloatType = real_t,
            typename Vector = std::vector<FloatType>
        >
        inline auto montecarlo(
            std::shared_ptr<random::random_context> rnd_ctx, unsigned int dimensions) {
 
            // Return an n-dimensional Monte Carlo estimator as a lambda function
            return [rnd_ctx, dimensions](
                std::function<FloatType(Vector)> funcApprox,
                std::function<FloatType(Vector)> funcExpected,
                estimate_options<FloatType, std::vector<FloatType>> options) {
 
                random::random_source rnd = rnd_ctx->get_rnd();
 
                if(options.domain.size() != dimensions) {
                    throw std::runtime_error(
                        "The estimation domain's dimension does not match "
                        "the instantiated number of dimensions "
                        "in estimator::montecarlo"
                    );
                }
 
                FloatType sum = 0;
                FloatType sumSqr = 0;
                FloatType sumAbs = 0;
                FloatType max = 0;
 
                // Compute the measure of a multi-interval
                FloatType volume = 1;
                for (interval k : options.domain)
                    volume *= k.length();
 
                Vector x (dimensions);
 
                for (unsigned int i = 0; i < options.iterations; ++i) {
                    
                    rnd.uniform(x, options.domain);
 
                    const FloatType diff = std::abs(funcApprox(x) - funcExpected(x));
 
                    if(max < diff)
                        max = diff;
 
                    sum += diff;
                    sumSqr += diff * diff;
                    sumAbs += std::abs(funcExpected(x));
                }
 
                estimate_result res {};
                res.maxErr = max;
                res.meanErr = sum / options.iterations;
                res.absErr = sum * (volume / options.iterations);
                res.rmsErr = std::sqrt(sumSqr * (volume / options.iterations));
                res.relErr = sum / sumAbs;
 
                return res;
            };
        }
    }
}}
 
 
#endif