sobol.h

Go to the documentation of this file.

 
#ifndef CMP_SOBOL_H
#define CMP_SOBOL_H
 
#include <cmp_defines.h>
#include <statistics.h>
#include <grid.h>
#include <stdexcept>
#include <utility>
#include <memory>
 
namespace cmp::sobol {
 
struct SobolResults {
    Eigen::VectorXd firstOrder;
 
    Eigen::VectorXd totalOrder;
 
    Eigen::VectorXd secondOrder;
 
    double getFirstOrder(size_t i) const {
        if(i >= firstOrder.size()) throw std::out_of_range("Index out of range");
        return firstOrder(i);
    }
 
    double getTotalOrder(size_t i) const {
        if(i >= totalOrder.size()) throw std::out_of_range("Index out of range");
        return totalOrder(i);
    }
 
    double getSecondOrder(size_t i, size_t j) const {
        if(secondOrder.size() == 0) throw std::logic_error("Second order Sobol indices were not computed.");
        if(i >= firstOrder.size() || j >= firstOrder.size()) throw std::out_of_range("Index out of range");
        if(i == j) throw std::invalid_argument("Indices must be different for second order Sobol index.");
        if(i > j) std::swap(i, j);
        size_t index = i * (2 * firstOrder.size() - i - 1) / 2 + (j - i - 1);
        return secondOrder(index);
    }
};
 
class SobolSaltelli {
  private:
    size_t nObs_ = 0;
 
    size_t nTotalObs_ = 0;
 
    size_t dim_ = 0;
 
    bool secondOrder_ = false;
 
  public:
    SobolSaltelli() = default;
 
    SobolSaltelli(size_t nObs, size_t dim, bool secondOrder = false)
        : nObs_(nObs), dim_(dim), secondOrder_(secondOrder) {
        nTotalObs_ = (2 + dim_ + (secondOrder_ ? dim_ * (dim_ - 1) / 2 : 0)) * nObs_;
    }
 
    size_t nTotalObs() const {
        return nTotalObs_;
    }
 
    size_t dim() const {
        return dim_;
    }
 
    size_t nObs() const {
        return nObs_;
    }
 
    Eigen::MatrixXd evaluationGrid(size_t nObs, std::shared_ptr<cmp::grid::Grid> gridType, bool secondOrder = false);
 
    Eigen::VectorXd sliceSaltelliOutput(const Eigen::VectorXd &Y, const Eigen::VectorXs &idx, size_t nObs, size_t dim, bool secondOrder);
 
    SobolResults compute(const Eigen::VectorXd &Y);
 
    std::pair<SobolResults, SobolResults> computeWithBootstrap(const Eigen::VectorXd &Y, size_t nBootstrap, size_t randomSeed = 42);
 
};
}
#endif // CMP_SOBOL_H