SmolyakEstimator.h

Go to the documentation of this file.

#ifndef SMOLYAKESTIMATOR_H
#define SMOLYAKESTIMATOR_H
 
#include "MUQ/Modeling/ModPiece.h"
#include "MUQ/Modeling/Flann/FlannCache.h"
 
#include "MUQ/Utilities/MultiIndices/MultiIndexSet.h"
 
#include <vector>
#include <boost/property_tree/ptree.hpp>
 
namespace muq {
namespace Approximation {
 
  template<typename EstimateType>
  class SmolyakEstimator {
  public:
 
    SmolyakEstimator(std::shared_ptr<muq::Modeling::ModPiece> const& modelIn);
 
    virtual ~SmolyakEstimator() = default;
 
    virtual EstimateType Compute(std::shared_ptr<muq::Utilities::MultiIndexSet> const& fixedSet,
                                 boost::property_tree::ptree                           options = boost::property_tree::ptree());
 
 
    virtual EstimateType Adapt(boost::property_tree::ptree options);
 
    virtual double Error() const{return globalError;};
 
    virtual unsigned int NumEvals() const{return numEvals;};
 
    virtual std::vector<double> ErrorHistory() const{return errorHistory;};
 
    virtual std::vector<int> EvalHistory() const{return evalHistory;};
 
    virtual std::vector<double> TimeHistory() const{return timeHistory;};
 
    virtual std::vector<std::vector<Eigen::VectorXd>> PointHistory() const;
 
    virtual std::vector<std::vector<std::shared_ptr<muq::Utilities::MultiIndex>>> TermHistory() const{return termHistory;};
 
  protected:
 
    virtual void Reset();
 
    virtual void AddTerms(std::shared_ptr<muq::Utilities::MultiIndexSet> const& fixedSet);
 
    virtual void AddTerms(std::vector<std::shared_ptr<muq::Utilities::MultiIndex>> const& fixedSet);
 
    virtual void UpdateErrors();
 
    virtual bool Refine();
 
    virtual void EvaluatePoints(std::set<unsigned int> const& ptsToEval);
 
    virtual std::vector<Eigen::VectorXd> OneTermPoints(std::shared_ptr<muq::Utilities::MultiIndex> const& multi) = 0;
 
    virtual EstimateType ComputeOneTerm(std::shared_ptr<muq::Utilities::MultiIndex>                const& multi,
                                        std::vector<std::reference_wrapper<const Eigen::VectorXd>> const& modEvals) = 0;
 
    virtual EstimateType ComputeWeightedSum(Eigen::VectorXd const& weights) const;
    virtual EstimateType ComputeWeightedSum() const;
 
 
    virtual EstimateType AddEstimates(double w1, EstimateType const& part1, double w2, EstimateType const& part2) const = 0;
 
    virtual double ComputeMagnitude(EstimateType const& estimate) const = 0;
 
    std::shared_ptr<muq::Modeling::ModPiece> model;
 
    std::shared_ptr<muq::Utilities::MultiIndexSet> termMultis;
 
    std::vector<double> errorHistory;
 
    std::vector<int> evalHistory;
 
    std::vector<double> timeHistory;
 
    std::vector<std::set<unsigned int>> pointHistory;
 
    std::vector<std::vector<std::shared_ptr<muq::Utilities::MultiIndex>>> termHistory;
 
    muq::Modeling::DynamicKDTreeAdaptor<> pointCache;
    std::vector<Eigen::VectorXd> evalCache;
 
    int InCache(Eigen::VectorXd const& input) const;
    Eigen::VectorXd const& GetFromCache(unsigned int index) const{return pointCache.m_data.at(index);};
    int AddToCache(Eigen::VectorXd const& newPt);
    int CacheSize() const{return pointCache.m_data.size();};
 
    const double cacheTol = 10.0*std::numeric_limits<double>::epsilon(); // <- points are considered equal if they are closer than this
 
    double timeTol = std::numeric_limits<double>::infinity();
 
    double errorTol = 0.0;
 
    unsigned int maxNumEvals = std::numeric_limits<unsigned int>::max();
 
    unsigned int numEvals=0;
 
    struct SmolyTerm {
 
      // Value of the tensor product approximation for one term in the Smolyak expansion
      EstimateType val;
 
      /* Weight on this tensor product approximation -- will change as terms are
         added to the Smolyak rule.
      */
      double weight = 0.0;
 
      /* Has val been computed for this term yet? */
      bool isComputed = false;
 
      /* Whether or not this term is part of the "old" set as defined by G&G */
      bool isOld = false;
 
      /* Is this term needed for the estimate or an error indicator? */
      bool isNeeded = false;
 
      // A local error indicator.  See eq (5.1) of Conrad and Marzouk
      double localError = -1.0;
 
      /* A vector containing the indices of points in the evaluation cache
        (see evalCache and pointCache variables) that are needed to compute the
        value for this term.  This is used for lazy evaluation and helps avoid
        reevaluation of the model.
      */
      std::vector<unsigned int> evalInds;
 
      Eigen::VectorXd diffWeights;
 
    };
 
    std::vector<SmolyTerm> terms;
 
    double globalError;
 
    const double nzTol = 10.0*std::numeric_limits<double>::epsilon();
 
  }; // class SmolyakEstimator
 
} // namespace muq
} // namespace Approximation
 
 
#endif