prec.h

Go to the documentation of this file.

 
 
#ifndef CHEBYSHEV_PREC_H
#define CHEBYSHEV_PREC_H
 
#include <string>
#include <vector>
#include <map>
#include <iostream>
#include <memory>
#include <mutex>
#include <thread>
 
#include "./prec/prec_structures.h"
#include "./core/output.h"
#include "./core/random.h"
#include "./prec/estimator.h"
#include "./prec/fail.h"
 
 
namespace chebyshev {
 
namespace prec {
 
 
    struct prec_settings {
        
        std::string moduleName = "unknown";
 
        unsigned int defaultIterations = CHEBYSHEV_PREC_ITER;
 
        FailFunction defaultFailFunction = fail::fail_on_max_err();
 
        long double defaultTolerance = CHEBYSHEV_PREC_TOLERANCE;
 
        std::vector<std::string> outputFiles {};
 
        std::vector<std::string> estimateColumns = {
            "name", "meanErr", "rmsErr", "maxErr", "failed"
        };
 
        std::vector<std::string> estimateOutputFiles {};
 
        std::vector<std::string> equationColumns = {
            "name", "difference", "tolerance", "failed"
        };
 
        std::vector<std::string> equationOutputFiles {};
 
        std::map<std::string, bool> pickedTests {};
 
        bool multithreading = true;
 
    };
 
 
    class prec_context {
    private:
 
        std::map<std::string, std::vector<estimate_result>> estimateResults {};
 
        std::mutex estimateMutex;
 
        std::vector<std::thread> estimateThreads;
 
        std::map<std::string, std::vector<equation_result>> equationResults {};
 
        bool wasTerminated {false};
 
    public:
 
        prec_settings settings;
 
        std::shared_ptr<output::output_context> output;
 
        std::shared_ptr<random::random_context> random;
 
 
        inline void setup(
            const std::string& moduleName,
            int argc = 0, const char** argv = nullptr) {
 
            // Initialize other modules
            settings = prec_settings();
            output = std::make_shared<output::output_context>();
            random = std::make_shared<random::random_context>();
 
            // Initialize list of picked tests
            if(argc && argv)
                for (int i = 1; i < argc; ++i)
                    settings.pickedTests[argv[i]] = true;
 
            std::cout << "Starting precision testing of the ";
            std::cout << moduleName << " module ..." << std::endl;
 
            settings.moduleName = moduleName;
            wasTerminated = false;
        }
 
 
        inline void terminate(bool exit = false) {
 
            // Ensure all test cases have been executed
            this->wait_results();
            // std::lock_guard<std::mutex> lock(estimateMutex);
 
            unsigned int totalTests = 0;
            unsigned int failedTests = 0;
 
            for (const auto& pair : estimateResults) {
                for (const auto& testCase : pair.second) {
                    totalTests++;
                    failedTests += testCase.failed ? 1 : 0;
                }
            }
 
            for (const auto& pair : equationResults) {
                for (const auto& testCase : pair.second) {
                    totalTests++;
                    failedTests += testCase.failed ? 1 : 0;
                }
            }
 
 
            // Ensure that an output file is specified
            if(  output->settings.outputToFile &&
                !output->settings.outputFiles.size() &&
                !settings.estimateOutputFiles.size() &&
                !settings.equationOutputFiles.size() &&
                !settings.outputFiles.size()) {
 
                settings.outputFiles = { settings.moduleName + "_results" };
            }
 
 
            // Print test results
            std::vector<std::string> outputFiles;
 
            // Print estimate results
            outputFiles  = settings.outputFiles;
            outputFiles.insert(
                outputFiles.end(),
                settings.estimateOutputFiles.begin(),
                settings.estimateOutputFiles.end()
            );
 
            output->print_results(
                estimateResults, settings.estimateColumns, outputFiles
            );
 
 
            // Print equation results
            outputFiles  = settings.outputFiles;
            outputFiles.insert(
                outputFiles.end(),
                settings.equationOutputFiles.begin(),
                settings.equationOutputFiles.end()
            );
 
            output->print_results(
                equationResults, settings.equationColumns, outputFiles
            );
 
 
            // Print overall test results
            std::cout << "Finished testing " << settings.moduleName << '\n';
            std::cout << totalTests << " total tests, ";
            std::cout << failedTests << " failed";
 
            // Print proportion of failed test, avoiding division by zero
            if (totalTests > 0) {
 
                const double percent = (failedTests / (double) totalTests) * 100;
                std::cout << " (" << std::setprecision(3) << percent << "%)" << std::endl;
                
            } else {
                std::cout << "\nNo tests were run!" << std::endl;
            }
 
            if(exit) {
                output->terminate();
                std::exit(failedTests);
            }
 
            wasTerminated = true;
        }
 
 
        prec_context(const std::string& moduleName, int argc, const char** argv) {
            setup(moduleName, argc, argv);
        }
 
 
        ~prec_context() {
            if (!wasTerminated)
                terminate();
        }
 
 
        prec_context(const prec_context& other) {
 
            std::lock_guard<std::mutex> lock(estimateMutex);
            estimateResults = other.estimateResults;
            equationResults = other.equationResults;
            settings = other.settings;
            output = other.output;
            random = other.random;
        }
 
 
        inline prec_context& operator=(const prec_context& other) {
 
            std::lock_guard<std::mutex> lock(estimateMutex);
            estimateResults = other.estimateResults;
            equationResults = other.equationResults;
            settings = other.settings;
            output = other.output;
            random = other.random;
 
            return *this;
        }
 
 
        template <
            typename R, typename ...Args,
            typename Function1 = std::function<R(Args...)>,
            typename Function2 = Function1
        >
        inline void estimate(
            const std::string& name,
            Function1 funcApprox,
            Function2 funcExpected,
            estimate_options<R, Args...> opt) {
 
            // Skip the test case if any tests have been picked
            // and this one was not picked.
            if(settings.pickedTests.size())
                if(settings.pickedTests.find(name) == settings.pickedTests.end())
                    return;
 
            auto task = [this, name, funcApprox, funcExpected, opt]() {
 
                // Use the estimator to estimate error integrals.
                auto res = opt.estimator(funcApprox, funcExpected, opt);
 
                res.name = name;
                res.domain = opt.domain;
                res.tolerance = opt.tolerance;
                res.quiet = opt.quiet;
                res.iterations = opt.iterations;
 
                // Use the fail function to determine whether the test failed.
                res.failed = opt.fail(res);
 
                std::lock_guard<std::mutex> lock(estimateMutex);
                estimateResults[name].push_back(res);
            };
 
            settings.multithreading ?
                estimateThreads.emplace_back(task) : task();
        }
 
 
        template <
            typename R, typename ...Args,
            typename Function1 = std::function<R(Args...)>,
            typename Function2 = Function1
        >
        inline void estimate(
            const std::string& name,
            Function1 funcApprox,
            Function2 funcExpected,
            std::vector<interval> domain,
            long double tolerance, unsigned int iterations,
            FailFunction fail,
            Estimator<R, Args...> estimator,
            bool quiet = false) {
 
            estimate_options<R, Args...> opt {};
            opt.domain = domain;
            opt.tolerance = tolerance;
            opt.iterations = iterations;
            opt.fail = fail;
            opt.estimator = estimator;
            opt.quiet = quiet;
 
            estimate(name, funcApprox, funcExpected, opt);
        }
 
 
        inline void estimate(
            const std::string& name,
            EndoFunction<double> funcApprox,
            EndoFunction<double> funcExpected,
            interval domain,
            long double tolerance = get_nan(),
            unsigned int iterations = 0,
            FailFunction fail = fail::fail_on_max_err(),
            Estimator<double, double> estimator = estimator::quadrature1D<double>(),
            bool quiet = false) {
 
            if (tolerance != tolerance)
                tolerance = settings.defaultTolerance;
 
            if (iterations == 0)
                iterations = settings.defaultIterations;
 
            estimate_options<double, double> opt {};
            opt.domain = { domain };
            opt.tolerance = tolerance;
            opt.iterations = iterations;
            opt.fail = fail;
            opt.estimator = estimator;
            opt.quiet = quiet;
 
            estimate(name, funcApprox, funcExpected, opt);
        }
 
 
        template <
            typename Type, typename Identity = EndoFunction<Type>
        >
        inline void identity(
            const std::string& name,
            Identity id,
            const estimate_options<Type, Type>& opt) {
 
            // Apply the identity function
            EndoFunction<Type> funcApprox = [id](Type x) -> Type {
                return id(x);
            };
 
            // And compare it to the identity
            EndoFunction<Type> funcExpected = [](Type x) -> Type {
                return x;
            };
 
            estimate(name, funcApprox, funcExpected, opt);
        }
 
 
        template <
            typename Type, typename Involution = EndoFunction<Type>
        >
        inline void involution(
            const std::string& name,
            Involution invol,
            const estimate_options<Type, Type>& opt) {
 
            // Apply the involution two times
            EndoFunction<Type> funcApprox = [invol](Type x) -> Type {
                return invol(invol(x));
            };
 
            // And compare it to the identity
            EndoFunction<Type> funcExpected = [](Type x) -> Type {
                return x;
            };
 
            estimate(name, funcApprox, funcExpected, opt);
        }
 
 
        template <
            typename Type, typename Involution = EndoFunction<Type>
        >
        inline void idempotence(
            const std::string& name,
            Involution idem,
            const estimate_options<Type, Type>& opt) {
 
            // Apply the idem two times
            EndoFunction<Type> funcApprox = [idem](Type x) -> Type {
                return idem(idem(x));
            };
 
            // And compare it to the identity
            EndoFunction<Type> funcExpected = [idem](Type x) -> Type {
                return idem(x);
            };
 
            estimate(name, funcApprox, funcExpected, opt);
        }
 
 
        template <
            typename InputType, typename OutputType = InputType,
            typename Homogeneous = std::function<OutputType(InputType)>
        >
        inline void homogeneous(
            const std::string& name,
            Homogeneous hom,
            const estimate_options<InputType, OutputType>& opt,
            OutputType zero_element = OutputType(0.0)) {
 
            // Apply the homogeneous function
            std::function<OutputType(InputType)> funcApprox =
                [=](InputType x) -> OutputType {
                    return hom(x);
                };
 
            // And compare it to the zero element
            std::function<OutputType(InputType)> funcExpected =
                [=](InputType x) -> OutputType {
                    return zero_element;
                };
 
            estimate(name, funcApprox, funcExpected, opt);
        }
 
 
        template<typename Type = double>
        inline void equals(
            const std::string& name,
            const Type& evaluated, const Type& expected,
            equation_options<Type> opt = equation_options<Type>()) {
 
            // Skip the test case if any tests have been picked
            // and this one was not picked.
            if(settings.pickedTests.size())
                if(settings.pickedTests.find(name) == settings.pickedTests.end())
                    return;
 
            equation_result res {};
 
            long double diff = opt.distance(evaluated, expected);
 
            // Mark the test as failed if the
            // distance between the two values
            // is bigger than the tolerance.
            res.failed = (diff > opt.tolerance);
 
            res.name = name;
            res.difference = diff;
            res.tolerance = opt.tolerance;
            res.quiet = opt.quiet;
 
            // Register the result of the equation by name
            equationResults[name].push_back(res);
        }
 
 
        template<typename Type = double>
        inline void equals(
            const std::string& name,
            const Type& evaluated, const Type& expected,
            long double tolerance,
            DistanceFunction<Type> distance,
            bool quiet = false) {
 
            equation_options<Type> opt {};
            opt.tolerance = tolerance;
            opt.distance = distance;
            opt.quiet = quiet;
 
            equals(name, evaluated, expected, opt);
        }
 
 
        inline void equals(
            const std::string& name,
            long double evaluated, long double expected,
            long double tolerance = get_nan(),
            bool quiet = false) {
 
            if (tolerance != tolerance)
                tolerance = settings.defaultTolerance;
 
            // Skip the test case if any tests have been picked
            // and this one was not picked.
            if(settings.pickedTests.size())
                if(settings.pickedTests.find(name) == settings.pickedTests.end())
                    return;
 
            equation_result res {};
 
            long double diff = distance::abs_distance(evaluated, expected);
 
            // Mark the test as failed if the
            // distance between the two values
            // is bigger than the tolerance.
            res.failed = (diff > tolerance);
 
            res.name = name;
            res.difference = diff;
            res.tolerance = tolerance;
            res.quiet = quiet;
 
            res.evaluated = evaluated;
            res.expected = expected;
 
            // Register the result of the equation by name
            equationResults[name].push_back(res);
        }
 
 
        template<typename Type>
        inline void equals(
            const std::string& name,
            std::vector<std::array<Type, 2>> values,
            long double tolerance = get_nan(),
            bool quiet = false) {
 
            if (tolerance != tolerance)
                tolerance = settings.defaultTolerance;
 
            // Skip the test case if any tests have been picked
            // and this one was not picked.
            if(settings.pickedTests.size())
                if(settings.pickedTests.find(name) == settings.pickedTests.end())
                    return;
 
            for (const auto& v : values)
                equals(name, v[0], v[1], tolerance, quiet);
        }
 
 
        inline void wait_results() {
 
            for (auto& t : estimateThreads)
                if (t.joinable())
                    t.join();
 
            estimateThreads.clear();
        }
 
 
        inline std::vector<estimate_result> get_estimate(const std::string& name) {
 
            this->wait_results();
 
            return estimateResults[name];
        }
 
 
        inline estimate_result get_estimate(const std::string& name, unsigned int i) {
 
            this->wait_results();
 
            return estimateResults[name].at(i);
        }
 
 
        inline std::vector<equation_result> get_equation(const std::string& name) {
            return equationResults[name];
        }
 
 
        inline equation_result get_equation(const std::string& name, unsigned int i) {
            return equationResults[name].at(i);
        }
    
    };
 
 
    prec_context make_context(const std::string& moduleName,
            int argc = 0, const char** argv = nullptr) {
        
        return prec_context(moduleName, argc, argv);
    }
 
 
}}
 
#endif