MUQ  0.4.3
LocalRegression.cpp
Go to the documentation of this file.
2 
3 #if MUQ_HAS_PARCER
4 #include <parcer/Eigen.h>
5 #endif
6 
7 namespace pt = boost::property_tree;
8 using namespace muq::Modeling;
9 using namespace muq::Approximation;
10 
11 LocalRegression::LocalRegression(std::shared_ptr<ModPiece> function, pt::ptree& pt) : ModPiece(function->inputSizes, function->outputSizes), kn(pt.get<unsigned int>("NumNeighbors")) {
12  SetUp(function, pt);
13 }
14 
15 #if MUQ_HAS_PARCER
16 LocalRegression::LocalRegression(std::shared_ptr<muq::Modeling::ModPiece> function, boost::property_tree::ptree& pt, std::shared_ptr<parcer::Communicator> comm) : ModPiece(function->inputSizes, function->outputSizes), kn(pt.get<unsigned int>("NumNeighbors")), comm(comm) /*tagSingle(comm->GetSize()+1), tagMulti(comm->GetSize()+2)*/ {
17  SetUp(function, pt);
18 }
19 
21  if( comm ) {
22  // needs to be destroyed on all processors
23  comm->Barrier();
24 
25  // clear any messages
26  Probe();
27  }
28 }
29 #endif
30 
31 void LocalRegression::SetUp(std::shared_ptr<muq::Modeling::ModPiece> function, boost::property_tree::ptree& pt) {
32  // can only have one input and output
33  assert(inputSizes.size()==1);
34  assert(outputSizes.size()==1);
35 
36  // create a cache of model evaluations
37  cache = std::make_shared<FlannCache>(function);
38 
39  // create a regression object
40  pt.put<std::string>("PolynomialBasis", pt.get<std::string>("PolynomialBasis", "Legendre")); // set default to Legendre
41  pt.put<unsigned int>("Order", pt.get<unsigned int>("Order", 2)); // set default order to 2
42  pt.put<double>("MaxPoisednessRadius", pt.get<double>("MaxPoisednessRadius", 1.0));
43  pt.put<unsigned int>("InputSize", function->inputSizes(0));
44  reg = std::make_shared<Regression>(pt);
45 }
46 
48  assert(cache);
49  cache = std::make_shared<FlannCache>(cache->Function());
50 }
51 
52 void LocalRegression::FitRegression(Eigen::VectorXd const& input) const {
53  // find the nearest neighbors
54  std::vector<Eigen::VectorXd> neighbors;
55  std::vector<Eigen::VectorXd> result;
56  cache->NearestNeighbors(input, kn, neighbors, result);
57 
58  // fit the regression
59  reg->Fit(neighbors, result, input);
60 }
61 
63 #if MUQ_HAS_PARCER
64  Probe();
65 #endif
66 
67  // fit the regressor
68  FitRegression(inputs[0]);
69 
70  // evaluate the regressor
71  outputs.resize(1);
72  outputs[0] = (Eigen::VectorXd)boost::any_cast<Eigen::MatrixXd const&>(reg->Evaluate(inputs[0].get()) [0]).col(0);
73 }
74 
75 unsigned int LocalRegression::CacheSize() const {
76  assert(cache);
77  return cache->Size();
78 }
79 
80 Eigen::VectorXd LocalRegression::CachePoint(unsigned int const index) const {
81  assert(cache);
82  return cache->at(index);
83 }
84 
85 bool LocalRegression::InCache(Eigen::VectorXd const& point) const {
86  assert(cache);
87  return cache->InCache(point)>=0;
88 }
89 
90 Eigen::VectorXd LocalRegression::Add(Eigen::VectorXd const& input) const {
91  assert(cache);
92  const Eigen::VectorXd result = cache->Add(input);
93 
94 #if MUQ_HAS_PARCER
95  if( comm ) {
96  for( unsigned int i=0; i<comm->GetSize(); ++i ) {
97  if( i==comm->GetRank() ) { continue; }
98 
99  comm->Send(std::pair<Eigen::VectorXd, Eigen::VectorXd>(input, result), i, tagSingle);
100  }
101  Probe();
102  }
103 #endif
104 
105  return result;
106 }
107 
108 void LocalRegression::Add(std::vector<Eigen::VectorXd> const& inputs) const {
109  for( auto it : inputs ) { Add(it); }
110 }
111 
112 std::tuple<Eigen::VectorXd, double, unsigned int> LocalRegression::PoisednessConstant(Eigen::VectorXd const& input) const {
113  // find the nearest neighbors
114  std::vector<Eigen::VectorXd> neighbors;
115  cache->NearestNeighbors(input, kn, neighbors);
116 
117  return PoisednessConstant(input, neighbors);
118 }
119 
120 std::tuple<Eigen::VectorXd, double, unsigned int> LocalRegression::PoisednessConstant(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd> const& neighbors) const {
121  assert(neighbors.size()>=kn);
122  assert(reg);
123  std::pair<Eigen::VectorXd, double> lambda = reg->PoisednessConstant(neighbors, input, kn);
124 
125  double dist = 0.0;
126  unsigned int index = 0;
127  for( unsigned int i=0; i<kn; ++i ) {
128  const double d = (input-neighbors[i]).norm();
129  if( d<dist ) { dist = d; index = i; }
130  }
131 
132  return std::tuple<Eigen::VectorXd, double, unsigned int>(lambda.first, lambda.second, index);
133 }
134 
135 std::pair<double, double> LocalRegression::ErrorIndicator(Eigen::VectorXd const& input) const {
136  // find the nearest neighbors
137  std::vector<Eigen::VectorXd> neighbors;
138  cache->NearestNeighbors(input, kn, neighbors);
139 
140  return ErrorIndicator(input, neighbors);
141 }
142 
143 unsigned int LocalRegression::Order() const { return reg->order; }
144 
145 std::pair<double, double> LocalRegression::ErrorIndicator(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd> const& neighbors) const {
146  assert(neighbors.size()==kn);
147 
148  // create a local factorial function (caution: may be problematic if n is sufficiently large)
149  //std::function<int(int)> factorial = [&factorial](int const n) { return ((n==2 || n==1)? n : n*factorial(n-1)); };
150 
151  // compute the radius
152  double delta = 0.0;
153  for( auto n : neighbors) { delta = std::max(delta, (n-input).norm()); }
154  /*Eigen::ArrayXd radius = Eigen::ArrayXd::Zero(input.size());
155  for( auto n : neighbors) {
156  //std::cout << n.transpose() << std::endl;
157  radius = radius.max((n-input).array().abs());
158  }*/
159 
160  return std::pair<double, double>(((double)reg->order+1.0)*std::log(delta), delta);
161 }
162 
163 void LocalRegression::NearestNeighbors(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd>& neighbors) const {
164  assert(cache);
165  cache->NearestNeighbors(input, kn, neighbors);
166 }
167 
168 void LocalRegression::NearestNeighbors(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd>& neighbors, std::vector<Eigen::VectorXd>& result) const {
169  assert(cache);
170  cache->NearestNeighbors(input, kn, neighbors, result);
171 }
172 
173 Eigen::VectorXd LocalRegression::EvaluateRegressor(Eigen::VectorXd const& input, std::vector<Eigen::VectorXd> const& neighbors, std::vector<Eigen::VectorXd> const& result) const {
174 #if MUQ_HAS_PARCER
175  Probe();
176 #endif
177 
178  // fit the regression
179  reg->Fit(neighbors, result, input);
180 
181  // evaluate the regressor
182  return (Eigen::VectorXd)boost::any_cast<Eigen::MatrixXd const&>(reg->Evaluate(input) [0]).col(0);
183 }
184 
185 #if MUQ_HAS_PARCER
187  if( !comm ) { return; }
188  //std::cout << "START PROBE!!! ";// /*<< comm->GetRank() */<< std::endl;
189 
190  for( unsigned int i=0; i<comm->GetSize(); ++i ) {
191  // std::cout << "CHECKING " << i << " of " << comm->GetSize() << " rank: " << comm->GetRank();// << std::endl;
192  //std::cout << "START RECIEVING " << i+1 << " of " << 2 << std::endl;
193  if( i==comm->GetRank() ) { continue; }
194  //std::cout << "START RECIEVING " << i << " of " << comm->GetSize() << " rank: " << comm->GetRank();// << std::endl;
195 
196  //{ // get single adds
197  //parcer::RecvRequest recvReq;
198  //bool enter = comm->Iprobe(i, tagSingle, recvReq);
199  while( true ) {
200  parcer::RecvRequest recvReq;
201  if( !comm->Iprobe(i, tagSingle, recvReq) ) { break; }
202  //std::cout << "PROBED ";
203  // get the point
204  comm->Irecv(i, tagSingle, recvReq);
205  //std::cout << "RECIEVED, proc: " << std::endl << std::flush;
206  //std::cout << "RECIEVED, proc: " << comm->GetRank() << std::endl << std::flush;
207  const std::pair<Eigen::VectorXd, Eigen::VectorXd> point = recvReq.GetObject<std::pair<Eigen::VectorXd, Eigen::VectorXd> >();
208  //std::cout << std::flush;
209  //auto point = recvReq.GetObject<std::pair<Eigen::VectorXd, Eigen::VectorXd> >();
210  //std::cout << "GOT OBJECT!" << " rank: " << comm->GetRank() << " ";
211  //std::cout << "GOT OBJECT!" << std::endl;
212 
213  // add the point
214  cache->Add(point.first, point.second);
215 
216  //std::cout << "ADDED!" << " rank: " << comm->GetRank() << " ";
217 
218  recvReq.Clear();
219 
220  //std::cout << "CLEARED!" << " rank: " << comm->GetRank() << " ";
221 
222  //enter = comm->Iprobe(i, tagSingle, recvReq);
223 
224  // std::cout << "PROBED!" << " rank: " << comm->GetRank();
225  }
226  //}
227 
228  //std::cout << "DONE WITH THIS BIT" << " rank: " << comm->GetRank() << " ";
229  //std::cout << "DONE RECIEVING " << i << " of " << comm->GetSize() << " ";// << std::endl;
230  }
231 
232 // std::cout << "DONE WITH PROBE!!!" << std::endl;
233 }
234 #endif
235 
236 Eigen::VectorXd LocalRegression::CacheCentroid() const {
237  assert(cache);
238  return cache->Centroid();
239 }
std::pair< double, double > ErrorIndicator(Eigen::VectorXd const &input) const
Get the error indicator.
Eigen::VectorXd CacheCentroid() const
Get the centroid of the cache.
virtual void EvaluateImpl(muq::Modeling::ref_vector< Eigen::VectorXd > const &inputs) override
void SetUp(std::shared_ptr< muq::Modeling::ModPiece > function, boost::property_tree::ptree &pt)
Set up the regressor.
unsigned int Order() const
Polynomial order.
Eigen::VectorXd EvaluateRegressor(Eigen::VectorXd const &input, std::vector< Eigen::VectorXd > const &neighbors, std::vector< Eigen::VectorXd > const &result) const
Eigen::VectorXd CachePoint(unsigned int const index) const
A point in the cache.
unsigned int CacheSize() const
Get the total size of the cache.
std::shared_ptr< muq::Modeling::FlannCache > cache
A cache containing previous model evaluations.
std::tuple< Eigen::VectorXd, double, unsigned int > PoisednessConstant(Eigen::VectorXd const &input) const
Get the poisedness constant.
std::shared_ptr< parcer::Communicator > comm
void Add(std::vector< Eigen::VectorXd > const &inputs) const
Add some points to the cache.
std::shared_ptr< Regression > reg
A regressor.
const unsigned int kn
The number of nearest neighbors to use by the regressor.
bool InCache(Eigen::VectorXd const &point) const
Is a point in the cache?
void FitRegression(Eigen::VectorXd const &input) const
Fit the regression to the nearest neighbors.
void NearestNeighbors(Eigen::VectorXd const &input, std::vector< Eigen::VectorXd > &neighbors) const
Get the number of nearest neighbors.
Provides an abstract interface for defining vector-valued model components.
Definition: ModPiece.h:148
const Eigen::VectorXi inputSizes
Definition: ModPiece.h:469
std::vector< Eigen::VectorXd > outputs
Definition: ModPiece.h:503
const Eigen::VectorXi outputSizes
Definition: ModPiece.h:472
std::vector< std::reference_wrapper< const T > > ref_vector
A vector of references to something ...
Definition: WorkPiece.h:37
auto get(const nlohmann::detail::iteration_proxy_value< IteratorType > &i) -> decltype(i.key())
Definition: json.h:3956