IMP logo
IMP Reference Guide  2.20.2
The Integrative Modeling Platform
bivariate_functions.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/bivariate_functions.h
3  * \brief Classes for general functions
4  *
5  * Copyright 2007-2022 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPISD_BIVARIATE_FUNCTIONS_H
9 #define IMPISD_BIVARIATE_FUNCTIONS_H
10 
11 #include <IMP/isd/isd_config.h>
12 #include <IMP/Particle.h>
13 #include <IMP/isd/Nuisance.h>
14 #include <IMP/isd/Scale.h>
15 #include <IMP/isd/Switching.h>
16 #include <IMP/Object.h>
17 #include <Eigen/Dense>
18 
19 #define IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM 1e-7
20 
21 IMPISD_BEGIN_NAMESPACE
22 
23 //! Base class for functions of two variables
24 class IMPISDEXPORT BivariateFunction : public Object {
25  public:
26  BivariateFunction(std::string str) : Object(str) {}
27 
28  //! evaluate the function at a certain point
29  virtual Floats operator()(const Floats& x1, const Floats& x2) const = 0;
30 
31  //! evaluate the function at a list of points
32  /* returns a NxN matrix of entries where N=xlist.rows()
33  * and entry ij of the matrix is f(xlist(i),xlist(j))
34  */
35  virtual Eigen::MatrixXd operator()(const IMP::FloatsList& xlist)
36  const = 0;
37 
38  //! used for testing only
39  virtual FloatsList operator()(const IMP::FloatsList& xlist,
40  bool stupid) const = 0;
41 
42  //! return true if internal parameters have changed.
43  virtual bool has_changed() const = 0;
44 
45  //! update internal parameters
46  virtual void update() = 0;
47 
48  //! update derivatives of particles
49  virtual void add_to_derivatives(const Floats& x1, const Floats& x2,
50  DerivativeAccumulator& accum) const = 0;
51 
52  //! update derivatives of particles
53  /* add to the given particle the specified derivative
54  * guarantees that the particle_no (starting at 0) matches with
55  * the columns of get_derivative_matrix.
56  */
57  virtual void add_to_particle_derivative(unsigned particle_no, double value,
58  DerivativeAccumulator& accum)
59  const = 0;
60 
61  //! return derivative matrix
62  /* m_ij = d(func(xlist[i],xlist[j]))/dparticle_no
63  * the matrix is NxN where N = xlist.size()
64  */
65  virtual Eigen::MatrixXd get_derivative_matrix(
66  unsigned particle_no, const FloatsList& xlist) const = 0;
67 
68  // for testing purposes
69  virtual FloatsList get_derivative_matrix(unsigned particle_no,
70  const FloatsList& xlist,
71  bool stupid) const = 0;
72 
73  //! return second derivative matrix
74  virtual Eigen::MatrixXd get_second_derivative_matrix(
75  unsigned particle_a, unsigned particle_b,
76  const FloatsList& xlist) const = 0;
77 
78  // for testing purposes
79  virtual FloatsList get_second_derivative_matrix(unsigned particle_a,
80  unsigned particle_b,
81  const FloatsList& xlist,
82  bool stupid) const = 0;
83 
84  //! returns the number of input dimensions
85  virtual unsigned get_ndims_x1() const = 0;
86  virtual unsigned get_ndims_x2() const = 0;
87 
88  //! returns the number of output dimensions
89  virtual unsigned get_ndims_y() const = 0;
90 
91  //! returns the number of particles that this function uses
92  virtual unsigned get_number_of_particles() const = 0;
93 
94  //! returns true if the particle whose index is provided is optimized
95  virtual bool get_particle_is_optimized(unsigned particle_no) const = 0;
96 
97  //! returns the number of particles that are optimized
98  virtual unsigned get_number_of_optimized_particles() const = 0;
99 
100  //! particle manipulation
101  virtual ModelObjectsTemp get_inputs() const = 0;
102 
104 };
105 //
106 //! Covariance function
107 /* \f[w(x,x') = \tau^2 \exp\left(-\frac{|x-x'|^\alpha}{2\lambda^\alpha} +
108  * \delta_{ij} J\f]
109  * \param[in] \f$\tau\f$ ISD Scale
110  * \param[in] \f$\lambda\f$ ISD Scale
111  * \param[in] \f$\alpha\f$ positive double, usually greater than 1.
112  * Default is 2.
113  * \param[in] J is some jitter. Try J=0.01 if you get NANs.
114  * \param[in] cutoff is a positive double indicating when to consider that
115  * values are zero (to avoid computing exponentials). cutoff is relative to the
116  * value when x=x', and affects only the call get_derivative_matrix.
117  */
118 class IMPISDEXPORT Covariance1DFunction : public BivariateFunction {
119  public:
120  Covariance1DFunction(Particle* tau, Particle* ilambda,
121  double alpha = 2.0, double jitter = 0.0,
122  double cutoff = 1e-7)
123  : BivariateFunction("Covariance1DFunction %1%"),
124  alpha_(alpha),
125  tau_(tau),
126  lambda_(ilambda),
127  J_(jitter),
128  cutoff_(cutoff) {
129  IMP_LOG_TERSE("Covariance1DFunction: constructor" << std::endl);
130  lambda_val_ = Scale(ilambda).get_nuisance();
131  tau_val_ = Scale(tau).get_nuisance();
132  do_jitter = (jitter > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
133  alpha_square_ = (std::abs(alpha - 2) < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
134  update();
135  }
136 
137  bool has_changed() const override {
138  double tmpt = Scale(tau_).get_nuisance();
139  double tmpl = Scale(lambda_).get_nuisance();
140  IMP_LOG_VERBOSE("Covariance1DFunction: has_changed(): ");
141  IMP_LOG_VERBOSE(tmpt << " " << tau_val_ << " ");
142  IMP_LOG_VERBOSE(tmpl << " " << lambda_val_ << " ");
143  if ((std::abs(tmpt - tau_val_) > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) ||
144  (std::abs(tmpl - lambda_val_) > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM)) {
145  IMP_LOG_VERBOSE("true" << std::endl);
146  return true;
147  } else {
148  IMP_LOG_VERBOSE("false" << std::endl);
149  return false;
150  }
151  }
152 
153  void update() override {
154  lambda_val_ = Scale(lambda_).get_nuisance();
155  tau_val_ = Scale(tau_).get_nuisance();
156  IMP_LOG_TERSE("Covariance1DFunction: update() tau:= "
157  << tau_val_ << " lambda:=" << lambda_val_ << std::endl);
158  IMP_INTERNAL_CHECK(!IMP::isnan(tau_val_), "tau is nan.");
159  IMP_INTERNAL_CHECK(!IMP::isnan(lambda_val_), "lambda is nan.");
160  IMP_INTERNAL_CHECK(tau_val_ >= 0, "tau is negative.");
161  IMP_INTERNAL_CHECK(lambda_val_ >= 0, "lambda is negative.");
162  IMP_INTERNAL_CHECK(lambda_val_ != 0, "lambda is zero.");
163  }
164 
165  Floats operator()(const Floats& x1, const Floats& x2) const override {
166  IMP_USAGE_CHECK(x1.size() == 1, "expecting a 1-D vector");
167  IMP_USAGE_CHECK(x2.size() == 1, "expecting a 1-D vector");
168  Floats ret(1, get_value(x1[0], x2[0]));
169  return ret;
170  }
171 
172  Eigen::MatrixXd operator()(const IMP::FloatsList& xlist) const override {
173  const unsigned M = xlist.size();
174  Eigen::MatrixXd Mret(M, M);
175  for (unsigned i = 0; i < M; i++) {
176  for (unsigned j = i; j < M; j++) {
177  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
178  IMP_USAGE_CHECK(xlist[j].size() == 1, "expecting a 1-D vector");
179  double x1 = xlist[i][0];
180  double x2 = xlist[j][0];
181  double ret = get_value(x1, x2);
182  Mret(i, j) = ret;
183  if (i != j) Mret(j, i) = ret;
184  }
185  }
186  return Mret;
187  }
188 
189  FloatsList operator()(const IMP::FloatsList& xlist, bool) const override {
190  Eigen::MatrixXd mat((*this)(xlist));
191  FloatsList ret;
192  for (unsigned i = 0; i < xlist.size(); i++)
193  for (unsigned j = 0; j < xlist.size(); j++)
194  ret.push_back(Floats(1, mat(i, j)));
195  return ret;
196  }
197 
198  void add_to_derivatives(const Floats& x1, const Floats& x2,
199  DerivativeAccumulator& accum) const override {
200  // d[w(x1,x2)]/dtau = 2/tau*(w(x1,x2))
201  double val = get_value(x1[0], x2[0]);
202  double tauderiv = 2. / tau_val_ * val;
203  IMP_INTERNAL_CHECK(!IMP::isnan(tauderiv), "tau derivative is nan.");
204  Scale(tau_).add_to_nuisance_derivative(tauderiv, accum);
205  // d[w(x,x')]/dlambda
206  // = w(x,x') ( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
207  double lambdaderiv =
208  val *
209  (alpha_ * std::pow((std::abs(x1[0] - x2[0]) / lambda_val_), alpha_) /
210  (2. * lambda_val_));
211  IMP_INTERNAL_CHECK(!IMP::isnan(lambdaderiv), "lambda derivative is nan.");
212  Scale(lambda_).add_to_nuisance_derivative(lambdaderiv, accum);
213  }
214 
215  void add_to_particle_derivative(unsigned particle_no, double value,
216  DerivativeAccumulator& accum) const override {
217  switch (particle_no) {
218  case 0: // tau
219  IMP_INTERNAL_CHECK(!IMP::isnan(value), "tau derivative is nan.");
220  Scale(tau_).add_to_nuisance_derivative(value, accum);
221  break;
222  case 1: // lambda
223  IMP_INTERNAL_CHECK(!IMP::isnan(value), "lambda derivative is nan.");
224  Scale(lambda_).add_to_nuisance_derivative(value, accum);
225  break;
226  default:
227  IMP_THROW("Invalid particle number", ModelException);
228  }
229  }
230 
231  Eigen::MatrixXd get_derivative_matrix(unsigned particle_no,
232  const FloatsList& xlist) const override {
233  // Strategy: fill in the main diagonal, then fill with zeros
234  // if the value of the function falls below cutoff.
235  // assumes data points are ordered!
236  unsigned N = xlist.size();
237  Eigen::MatrixXd ret(N, N);
238  double diag;
239  switch (particle_no) {
240  case 0: // tau
241  // d[w(x1,x1)]/dtau
242  // = 2/tau*w(x1,x1)
243  diag = get_value(xlist[0][0], xlist[0][0]);
244  diag *= 2. / tau_val_;
245  break;
246  case 1: // lambda
247  // d[w(x,x)]/dlambda
248  //= w(x,x)
249  // *( alpha /(2 lambda^{alpha+1}))
250  diag = 0;
251  break;
252  default:
253  IMP_THROW("Invalid particle number", ModelException);
254  }
256  "derivative matrix is nan on the diagonal.");
257  for (unsigned i = 0; i < N; i++) ret(i, i) = diag;
258  //
259  bool initial_loop = true;
260  double abs_cutoff = cutoff_ * diag;
261  double dmax = -1;
262  for (unsigned i = 0; i < N; i++) {
263  for (unsigned j = i + 1; j < N; j++) {
264  double x1(xlist[i][0]), x2(xlist[j][0]);
265  double val;
266  double dist(std::abs(x1 - x2));
267  // compute all entries as long as the cutoff distance was
268  // not recorded (initial_loop) or as long as the distance is
269  // smaller than the cutoff distance
270  if (initial_loop || dist <= dmax) {
271  switch (particle_no) {
272  case 0: // tau
273  // d[w(x1,x2)]/dtau
274  // = 2/tau*w(x1,x2)
275  val = get_value(xlist[i][0], xlist[j][0]);
276  val *= 2. / tau_val_;
277  break;
278  case 1: // lambda
279  // d[w(x,x')]/dlambda
280  //= w(x,x')
281  // *( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
282  if (dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
283  val = 0;
284  } else {
285  val = get_value(xlist[i][0], xlist[j][0]);
286  val *= alpha_ * std::pow((dist / lambda_val_), alpha_) /
287  (2. * lambda_val_);
288  }
289  break;
290  default:
291  IMP_THROW("Invalid particle number", ModelException);
292  }
293  // the value has been computed and is in val
294  // now check if it is smaller than the cutoff.
295  // If true change the flag and update the distance
296  if (std::abs(val) <= abs_cutoff) {
297  if (initial_loop) {
298  initial_loop = false;
299  dmax = dist;
300  } else if (dist < dmax) {
301  dmax = dist;
302  }
303  }
304  } else { // e.g. initial_loop == false && dist > dmax
305  val = 0;
306  }
308  "derivative matrix is nan at position("
309  << i << "," << j << ").");
310  ret(i, j) = val;
311  ret(j, i) = val;
312  }
313  }
314  return ret;
315  }
316 
317  FloatsList get_derivative_matrix(unsigned particle_no,
318  const FloatsList& xlist, bool) const override {
319  Eigen::MatrixXd mat(get_derivative_matrix(particle_no, xlist));
320  FloatsList ret;
321  for (int i = 0; i < mat.rows(); i++) {
322  Floats line;
323  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
324  ret.push_back(line);
325  }
326  return ret;
327  }
328 
330  unsigned particle_a, unsigned particle_b,
331  const FloatsList& xlist) const override {
332  unsigned N(xlist.size());
333  Eigen::MatrixXd ret(N, N);
334  if (particle_a > 1) IMP_THROW("Invalid particle 1 number", ModelException);
335  if (particle_b > 1) IMP_THROW("Invalid particle 2 number", ModelException);
336  // d2f/dtau2
337  if (particle_a == 0 && particle_b == 0) {
338  for (unsigned i = 0; i < N; i++) {
339  for (unsigned j = i; j < N; j++) {
340  double dist = std::abs(xlist[i][0] - xlist[j][0]);
341  double exponent = std::pow(dist / lambda_val_, alpha_);
342  double expterm = std::exp(-0.5 * exponent);
343  ret(i, j) = 2 * expterm;
344  if (i != j) ret(j, i) = ret(i, j);
345  }
346  }
347  } else if (particle_a == 1 && particle_b == 1) { // d2f/dlambda2
348  for (unsigned i = 0; i < N; i++) {
349  for (unsigned j = i; j < N; j++) {
350  double dist = std::abs(xlist[i][0] - xlist[j][0]);
351  double exponent = std::pow(dist / lambda_val_, alpha_);
352  double expterm = std::exp(-0.5 * exponent);
353  ret(i, j) = tau_val_ * tau_val_ * expterm * exponent /
354  (lambda_val_ * lambda_val_) * alpha_ / 2 *
355  (alpha_ / 2 * exponent - (alpha_ + 1));
356  if (i != j) ret(j, i) = ret(i, j);
357  }
358  }
359  } else { // d2f/d(tau)d(lambda)
360  for (unsigned i = 0; i < N; i++) {
361  for (unsigned j = i; j < N; j++) {
362  double dist = std::abs(xlist[i][0] - xlist[j][0]);
363  double exponent = std::pow(dist / lambda_val_, alpha_);
364  double expterm = std::exp(-0.5 * exponent);
365  ret(i, j) = tau_val_ * alpha_ * expterm / lambda_val_ * exponent;
366  if (i != j) ret(j, i) = ret(i, j);
367  }
368  }
369  }
370  return ret;
371  }
372 
374  unsigned particle_a, unsigned particle_b,
375  const FloatsList& xlist, bool) const override {
376  Eigen::MatrixXd mat(
377  get_second_derivative_matrix(particle_a, particle_b, xlist));
378  FloatsList ret;
379  for (int i = 0; i < mat.rows(); i++) {
380  Floats line;
381  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
382  ret.push_back(line);
383  }
384  return ret;
385  }
386 
387  unsigned get_ndims_x1() const override { return 1; }
388  unsigned get_ndims_x2() const override { return 1; }
389  unsigned get_ndims_y() const override { return 1; }
390 
391  unsigned get_number_of_particles() const override { return 2; }
392 
393  bool get_particle_is_optimized(unsigned particle_no) const override {
394  switch (particle_no) {
395  case 0: // tau
396  return Scale(tau_).get_nuisance_is_optimized();
397  case 1: // lambda
398  return Scale(lambda_).get_nuisance_is_optimized();
399  default:
400  IMP_THROW("Invalid particle number", ModelException);
401  }
402  }
403 
404  unsigned get_number_of_optimized_particles() const override {
405  unsigned count = 0;
406  if (Scale(tau_).get_nuisance_is_optimized()) count++;
407  if (Scale(lambda_).get_nuisance_is_optimized()) count++;
408  return count;
409  }
410 
411  ModelObjectsTemp get_inputs() const override {
412  ModelObjectsTemp ret;
413  ret.push_back(tau_);
414  ret.push_back(lambda_);
415  return ret;
416  }
417 
418  /*IMP_OBJECT_INLINE(Covariance1DFunction,
419  out << "covariance function with alpha = "
420  << alpha_ << std::endl, {});*/
422 
423  private:
424  inline double get_value(double x1, double x2) const {
425  double dist = std::abs(x1 - x2);
426  double ret = dist / lambda_val_;
427  if (alpha_square_) {
428  ret *= ret;
429  } else {
430  ret = std::pow(ret, alpha_);
431  }
432  ret = IMP::square(tau_val_) * std::exp(-0.5 * ret);
433  if (do_jitter && dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
434  ret += IMP::square(tau_val_) * J_;
435  }
437  "function value is nan. tau = "
438  << tau_val_ << " lambda = " << lambda_val_
439  << " q1 = " << x1 << " q2 = " << x2);
440  return ret;
441  }
442 
443  private:
444  double alpha_;
445  Pointer<Particle> tau_, lambda_;
446  double tau_val_, lambda_val_, J_, cutoff_, alpha_square_;
447  bool do_jitter;
448 };
449 
450 IMPISD_END_NAMESPACE
451 
452 #endif /* IMPISD_BIVARIATE_FUNCTIONS_H */
Eigen::MatrixXd operator()(const IMP::FloatsList &xlist) const override
evaluate the function at a list of points
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
A decorator for switching parameters particles.
unsigned get_ndims_x1() const override
returns the number of input dimensions
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const override
update derivatives of particles
A decorator for scale parameters particles.
A decorator for nuisance parameters particles.
Add scale parameter to particle.
Definition: Scale.h:24
#define IMP_LOG_VERBOSE(expr)
Definition: log_macros.h:83
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Set up destructor for a ref counted object.
#define IMP_LOG_TERSE(expr)
Definition: log_macros.h:72
Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const override
return derivative matrix
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Definition: check_macros.h:139
void update() override
update internal parameters
FloatsList operator()(const IMP::FloatsList &xlist, bool) const override
used for testing only
Common base class for heavy weight IMP objects.
Definition: Object.h:111
virtual Eigen::MatrixXd get_second_derivative_matrix(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const =0
return second derivative matrix
virtual void update()=0
update internal parameters
bool has_changed() const override
return true if internal parameters have changed.
Floats operator()(const Floats &x1, const Floats &x2) const override
evaluate the function at a certain point
unsigned get_number_of_particles() const override
returns the number of particles that this function uses
ModelObjectsTemp get_inputs() const override
particle manipulation
Base class for functions of two variables.
bool isnan(const T &a)
Return true if a number is NaN.
Definition: math.h:24
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
Definition: check_macros.h:50
Eigen::MatrixXd get_second_derivative_matrix(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const override
return second derivative matrix
A shared base class to help in debugging and things.
Object(std::string name)
Construct an object with the given name.
Class to handle individual particles of a Model object.
Definition: Particle.h:43
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
unsigned get_ndims_y() const override
returns the number of output dimensions
void add_to_derivatives(const Floats &x1, const Floats &x2, DerivativeAccumulator &accum) const override
update derivatives of particles
An exception which is thrown when the Model has attributes with invalid values.
Definition: exception.h:188
Class for adding derivatives from restraints to the model.
unsigned get_number_of_optimized_particles() const override
returns the number of particles that are optimized
bool get_particle_is_optimized(unsigned particle_no) const override
returns true if the particle whose index is provided is optimized
virtual Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const =0
return derivative matrix