1 #ifndef LOCALREGRESSION_H_
2 #define LOCALREGRESSION_H_
3 
4 #include "MUQ/config.h"
5 
6 #if MUQ_HAS_PARCER
7 #include <parcer/Communicator.h>
8 #endif
9 
10 #include "MUQ/Modeling/ModPiece.h"
11 #include "MUQ/Modeling/Flann/FlannCache.h"
12 
13 #include "MUQ/Approximation/Regression/Regression.h"
14 
15 namespace muq {
16   namespace Approximation {
17     class LocalRegression : public muq::Modeling::ModPiece {
18     public:
19 
20       /**
21 	 @param[in] function The function we wish to approximate with a local polynomial
22 	 @param[in] pt Options for the regression
23        */
24       LocalRegression(std::shared_ptr<muq::Modeling::ModPiece> function, boost::property_tree::ptree& pt);
25 
26 #if MUQ_HAS_PARCER
27       /**
28 	 @param[in] function The function we wish to approximate with a local polynomial
29 	 @param[in] pt Options for the regression
30 	 @param[in] comm The parcer communicator
31        */
32       LocalRegression(std::shared_ptr<muq::Modeling::ModPiece> function, boost::property_tree::ptree& pt, std::shared_ptr<parcer::Communicator> comm);
33 #endif
34 
35 #if MUQ_HAS_PARCER
36       ~LocalRegression();
37 #else
38       ~LocalRegression() = default;
39 #endif
40 
41       /// Add some points to the cache
42       /**
43 	 @param[in] inputs Points to add to the cache
44        */
45       void Add(std::vector<Eigen::VectorXd> const& inputs) const;
46 
47       /// Add a single point to the cache
48       /**
49 	 @param[in] input Point to add to the cache
50 	 \return The function result at the new point
51        */
52       Eigen::VectorXd Add(Eigen::VectorXd const& input) const;
53 
54       /// Get the total size of the cache
55       /**
56 	 \return The cache size
57        */
58       unsigned int CacheSize() const;
59 
60       /// Is a point in the cache?
61       /**
62       @param[in] point We want to know if this point is in the cache
63       \return true: it is in the cache, false: it is not in the cache
64       */
65       bool InCache(Eigen::VectorXd const& point) const;
66 
67       /// A point in the cache
68       /**
69 	 @param[in] index The index of the input
70 	 \return The input point in the cache associated with the index
71        */
72       Eigen::VectorXd CachePoint(unsigned int const index) const;
73 
74       /// Get the poisedness constant
75       /**
76 	 Get the poisedness constant associated with the nearest neighbors of the input point.
77 	 @param[in] input The input point
78 	 \return first: the point where the Lagrange polynomials are maximized, second: the poisedness constant, third: the index of the nearest neighbor closest to the new point
79        */
80       std::tuple<Eigen::VectorXd, double, unsigned int> PoisednessConstant(Eigen::VectorXd const& input) const;
81 
82       /// Get the poisedness constant
83       /**
84 	 Get the poisedness constant associated with the nearest neighbors of the input point.
85 	 @param[in] input The input point
86 	 @param[in] neighbors The nearest neighbors
87 	 \return first: the point where the Lagrange polynomials are maximized, second: the poisedness constant, third: the index of the nearest neighbor closest to the new point
88        */
89       std::tuple<Eigen::VectorXd, double, unsigned int> PoisednessConstant(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd> const& neighbors) const;
90 
91       /// Get the error indicator
92       /**
93 	 Get the error indicator \f$\Lambda \sqrt{k} \Delta^{p+1}\f$
94 	 @param[in] input The input point
95 	 \return first: the error indicator, second: the radius of the ball
96        */
97       std::pair<double, double> ErrorIndicator(Eigen::VectorXd const& input) const;
98 
99       /// Get the error indicator
100       /**
101 	 Get the error indicator \f$\Lambda \sqrt{k} \Delta^{p+1}\f$
102 	 @param[in] input The input point
103 	 @param[in] neighbors The nearest neighbors
104    \return first: the error indicator, second: the radius of the ball
105        */
106       std::pair<double, double> ErrorIndicator(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd> const& neighbors) const;
107 
108       /// Get the number of nearest neighbors
109       /**
110 	 @param[in] input We want the \f$k\f$ nearest neighbors to this point
111 	 @param[out] neighbors The \f$k\f$ nearest neighbors
112        */
113       void NearestNeighbors(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd>& neighbors) const;
114 
115       /// Get the number of nearest neighbors (with result)
116       /**
117 	 @param[in] input We want the \f$k\f$ nearest neighbors to this point
118 	 @param[out] neighbors The \f$k\f$ nearest neighbors
119 	 @param[out] results The corresponding output of the function
120        */
121       void NearestNeighbors(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd>& neighbors, std::vector<Eigen::VectorXd>& result) const;
122 
123       Eigen::VectorXd EvaluateRegressor(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd> const& neighbors, std::vector<Eigen::VectorXd> const& result) const;
124 
125 #if MUQ_HAS_PARCER
126       void Probe() const;
127 #endif
128 
129       /// The number of nearest neighbors to use by the regressor
130       const unsigned int kn;
131 
132       /// Get the centroid of the cache
133       /**
134       \return The cache centroid
135       */
136       Eigen::VectorXd CacheCentroid() const;
137 
138       /// Clear the cache
139       void ClearCache();
140 
141       /// Polynomial order
142       unsigned int Order() const;
143 
144     private:
145 
146       virtual void EvaluateImpl(muq::Modeling::ref_vector<Eigen::VectorXd> const& inputs) override;
147 
148       /// Fit the regression to the nearest neighbors
149       void FitRegression(Eigen::VectorXd const& input) const;
150 
151       /// Set up the regressor
152       void SetUp(std::shared_ptr<muq::Modeling::ModPiece> function, boost::property_tree::ptree& pt);
153 
154       /// A cache containing previous model evaluations
155       std::shared_ptr<muq::Modeling::FlannCache> cache;
156 
157       /// A regressor
158       std::shared_ptr<Regression> reg;
159 
160 #if MUQ_HAS_PARCER
161       std::shared_ptr<parcer::Communicator> comm = nullptr;
162       const int tagSingle = 0;
163 #endif
164     };
165   } // namespace Approximation
166 } // namespace muq
167 
168 #endif
169