SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
Basis.cpp
1 #include <stdexcept>
2 #include "Validator/ObjectRequirement.h"
3 #include "Drivers/DriverException.h"
4 #include "Basis.h"
5 #include "schema.h"
6 
7 namespace SSAGES
8 {
9  BasisFunction* BasisFunction::Build(const Json::Value& json, const std::string& path, unsigned int nbins)
10  {
11  auto type = json.get("type","none").asString();
12  if(type == "Legendre")
13  return Legendre::Build(json, path, nbins);
14  else if (type == "Chebyshev")
15  return Chebyshev::Build(json, path, nbins);
16  else if (type == "Fourier")
17  return Fourier::Build(json, path, nbins);
18  else
19  throw std::invalid_argument("Invalid basis set type \"" + type + "\".");
20  }
21 
22  Chebyshev* Chebyshev::Build(const Json::Value& json,
23  const std::string& path,
24  unsigned int nbin)
25  {
26  Json::ObjectRequirement validator;
27  Json::Value schema;
28  Json::CharReaderBuilder rbuilder;
29  Json::CharReader* reader = rbuilder.newCharReader();
30 
31  reader->parse(JsonSchema::ChebyshevBasis.c_str(),
32  JsonSchema::ChebyshevBasis.c_str() + JsonSchema::ChebyshevBasis.size(),
33  &schema, NULL);
34  validator.Parse(schema, path);
35 
36  //Validate Inputs
37  validator.Validate(json, path);
38  if(validator.HasErrors())
39  throw BuildException(validator.GetErrors());
40 
41  return new Chebyshev(
42  json["polynomial_order"].asInt(),
43  json["lower_bound"].asDouble(),
44  json["upper_bound"].asDouble(),
45  nbin
46  );
47  }
48 
49  Fourier* Fourier::Build(const Json::Value& json,
50  const std::string& path,
51  unsigned int nbin)
52  {
53  Json::ObjectRequirement validator;
54  Json::Value schema;
55  Json::CharReaderBuilder rbuilder;
56  Json::CharReader* reader = rbuilder.newCharReader();
57 
58  reader->parse(JsonSchema::FourierBasis.c_str(),
59  JsonSchema::FourierBasis.c_str() + JsonSchema::FourierBasis.size(),
60  &schema, NULL);
61  validator.Parse(schema, path);
62 
63  //Validate Inputs
64  validator.Validate(json, path);
65  if(validator.HasErrors())
66  throw BuildException(validator.GetErrors());
67 
68  return new Fourier(
69  json["polynomial_order"].asInt(),
70  json["lower_bound"].asDouble(),
71  json["upper_bound"].asDouble(),
72  nbin
73  );
74  }
75  Legendre* Legendre::Build(const Json::Value& json,
76  const std::string& path,
77  unsigned int nbin)
78  {
79  Json::ObjectRequirement validator;
80  Json::Value schema;
81  Json::CharReaderBuilder rbuilder;
82  Json::CharReader* reader = rbuilder.newCharReader();
83 
84  reader->parse(JsonSchema::LegendreBasis.c_str(),
85  JsonSchema::LegendreBasis.c_str() + JsonSchema::LegendreBasis.size(),
86  &schema, NULL);
87  validator.Parse(schema, path);
88 
89  //Validate Inputs
90  validator.Validate(json, path);
91  if(validator.HasErrors())
92  throw BuildException(validator.GetErrors());
93 
94  return new Legendre(
95  json["polynomial_order"].asInt(),
96  nbin
97  );
98  }
99 
101  {
102  double coeff_size = 1;
103  //Get the size of the number of coefficients + 1
104  for(size_t i = 0; i < functions_.size(); ++i)
105  {
106  coeff_size *= functions_[i]->GetOrder()+1;
107  }
108 
109  std::vector<int> jdx(functions_.size(), 0);
110  Map temp_map(jdx,0.0);
111  //Initialize the mapping for the coeff function
112  for(size_t i = 0; i < coeff_size; ++i)
113  {
114  for(size_t j = 0; j < jdx.size(); ++j)
115  {
116  if(jdx[j] > 0 && jdx[j] % (functions_[j]->GetOrder()+1) == 0)
117  {
118  if(j != functions_.size() - 1)
119  {
120  jdx[j+1]++;
121  jdx[j] = 0;
122  }
123  }
124  temp_map.map[j] = jdx[j];
125  temp_map.value = 0;
126  }
127  coeff_.push_back(temp_map);
128  jdx[0]++;
129  }
130  }
131 
133  {
134  for (size_t i=0; i<functions_.size(); i++)
135  {
136  unsigned int poly = functions_[i]->GetOrder()+1;
137  unsigned int nbin = functions_[i]->GetBins();
138  std::vector<double> dervs(nbin*poly,0);
139  std::vector<double> vals(nbin*poly,0);
140 
141  for (size_t j=0; j<poly; j++) {
142  for(size_t k=0; k<nbin; k++) {
143  double x = ((functions_[i]->GetRange())*k
144  - functions_[i]->GetLower())/nbin + functions_[i]->GetLower();
145  vals[k+j*nbin] = functions_[i]->Evaluate(x,j);
146  dervs[k+j*nbin] = functions_[i]->EvalGrad(x,j);
147  }
148  }
149  BasisLUT TempLUT(vals,dervs);
150  lookup_.push_back(TempLUT);
151  }
152  }
153 
154  void BasisEvaluator::UpdateBias(Grid<double> *bias, Grid<std::vector<double>> *grad)
155  {
156  double basis;
157  double temp;
158  size_t j = 0;
159  int nbins;
160 
161  for(Grid<double>::iterator it = bias->begin(); it != bias->end(); ++it, ++j)
162  {
163  std::vector<double> tmp_grad (functions_.size(),0);
164  double tmp_bias = 0;
165  for (size_t i = 1; i < coeff_.size(); ++i)
166  {
167  basis = 1.0;
168  for (size_t l = 0; l < functions_.size(); l++)
169  {
170  temp = 1.0;
171  //For the gradients we only evaluate the gradient for diagonal terms
172  //Off diagonals are the basis set value
173  for (size_t k = 0; k < functions_.size(); k++)
174  {
175  nbins = bias->GetNumPoints(k);
176  temp *= l == k ? lookup_[k].derivs[it.index(k) + coeff_[i].map[k]*nbins] * functions_[l]->GetRange() / (bias->GetUpper(l) - bias->GetLower(l))
177  : lookup_[k].values[it.index(k) + coeff_[i].map[k]*nbins];
178  }
179  tmp_grad[l] -= coeff_[i].value * temp;
180  nbins = bias->GetNumPoints(l);
181  basis *= lookup_[l].values[it.index(l) + coeff_[i].map[l]*nbins];
182  }
183  //Update the bias values
184  tmp_bias += coeff_[i].value * basis;
185  //Store the gradient values
186  }
187  grad->at(it.coordinates()) = tmp_grad;
188  *it = tmp_bias;
189  }
190  }
191 
192  //Calculates the inner product of the basis set and the biased histogram
193  //This function then returns the coefficients from this calculation
194  double BasisEvaluator::UpdateCoeff(const std::vector<double> &array, Grid<unsigned int> *hist)
195  {
196  double coeffTemp = 0;
197  double sum = 0;
198 
199  for(auto& coeff : coeff_)
200  {
201  // The method uses a standard integration
202  size_t j = 0;
203  coeffTemp = coeff.value;
204  coeff.value = 0.0;
205  for(Grid<unsigned int>::iterator it2 = hist->begin(); it2 != hist->end(); ++it2, ++j)
206  {
207 
208  double weight = std::pow(2.0,functions_.size());
209 
210  // This adds in a trap-rule type weighting which lowers error significantly at the boundaries
211  for(size_t k = 0; k < functions_.size(); ++k)
212  {
213  if( it2.index(k) == 0 ||
214  it2.index(k) == hist->GetNumPoints(k)-1)
215  weight /= 2.0;
216  }
217  /*The numerical integration of the biased histogram across the entirety of CV space
218  *All calculations include the normalization as well
219  */
220  double basis = 1.0;
221 
222  for(size_t l = 0; l < functions_.size(); l++)
223  {
224  int nbins = hist->GetNumPoints(l);
225  basis *= lookup_[l].values[it2.index(l) + coeff.map[l]*nbins]*functions_[l]->GetRange()/nbins;
226  //Normalize the values by the associated value
227  basis *= functions_[l]->GetNorm(coeff.map[l])*functions_[l]->Weight(it2.coordinate(l));
228  }
229  coeff.value += basis * array[j] * weight/std::pow(2.0,functions_.size());
230  }
231  coeffTemp -= coeff.value;
232  //std::cout<<coeffTemp<<std::endl;
233  sum += fabs(coeffTemp);
234  }
235  return sum;
236  }
237 };
const std::vector< double > GetLower() const
Return the lower edges of the Grid.
Definition: GridBase.h:226
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
static Fourier * Build(const Json::Value &json, const std::string &path, unsigned int nbins)
Build the Fourier polynomial.
Definition: Basis.cpp:49
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
void CoeffInit()
Initializes coefficient lookup vector.
Definition: Basis.cpp:100
void BasisInit()
Creates lookup table for basis functions.
Definition: Basis.cpp:132
Map for histogram and coefficients.
Definition: Basis.h:41
iterator begin()
Return iterator at first grid point.
Definition: Grid.h:527
Exception to be thrown when building the Driver fails.
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
Abstract class for all BasisFunction inheritance.
Definition: Basis.h:61
const std::vector< double > GetUpper() const
Return the upper edges of the Grid.
Definition: GridBase.h:257
Look-up table for basis functions.
Definition: Basis.h:16
const std::vector< int > GetNumPoints() const
Get the number of points for all dimensions.
Definition: GridBase.h:200
Requirements on an object.
static Legendre * Build(const Json::Value &json, const std::string &path, unsigned int nbins)
Build the Legendre polynomial.
Definition: Basis.cpp:75
static Chebyshev * Build(const Json::Value &json, const std::string &path, unsigned int nbins)
Build the Chebyshev polynomial.
Definition: Basis.cpp:22
double UpdateCoeff(const std::vector< double > &array, Grid< unsigned int > *hist)
Definition: Basis.cpp:194
Defines the class of Legendre polynomials.
Definition: Basis.h:241
void UpdateBias(Grid< double > *bias, Grid< std::vector< double >> *grad)
Outputs the basis projection at a specific coordinate.
Definition: Basis.cpp:154
unsigned int GetOrder()
Gets the order of the current polynomial.
Definition: Basis.h:96
static BasisFunction * Build(const Json::Value &json, const std::string &path, unsigned int nbins)
Build BasisFunction from JSON value.
Definition: Basis.cpp:9
Defines the class of Chebyshev polynomials.
Definition: Basis.h:183
Defines the class of Fourier polynomials.
Definition: Basis.h:284
double GetLower()
Gets the lower bound on the CV.
Definition: Basis.h:105
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.
iterator end()
Return iterator after last valid grid point.
Definition: Grid.h:540