IMP Reference Guide  develop.548de65454,2020/11/26 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  *
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 {
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() {
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 {
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 {
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 {
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 {
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.");
205  // d[w(x,x')]/dlambda
206  // = w(x,x') ( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
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.");
213  }
214
215  void add_to_particle_derivative(unsigned particle_no, double value,
216  DerivativeAccumulator& accum) const {
217  switch (particle_no) {
218  case 0: // tau
219  IMP_INTERNAL_CHECK(!IMP::isnan(value), "tau derivative is nan.");
221  break;
222  case 1: // lambda
223  IMP_INTERNAL_CHECK(!IMP::isnan(value), "lambda derivative is nan.");
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 {
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 {
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, const FloatsList& xlist) const {
331  unsigned N(xlist.size());
332  Eigen::MatrixXd ret(N, N);
333  if (particle_a > 1) IMP_THROW("Invalid particle 1 number", ModelException);
334  if (particle_b > 1) IMP_THROW("Invalid particle 2 number", ModelException);
335  // d2f/dtau2
336  if (particle_a == 0 && particle_b == 0) {
337  for (unsigned i = 0; i < N; i++) {
338  for (unsigned j = i; j < N; j++) {
339  double dist = std::abs(xlist[i][0] - xlist[j][0]);
340  double exponent = std::pow(dist / lambda_val_, alpha_);
341  double expterm = std::exp(-0.5 * exponent);
342  ret(i, j) = 2 * expterm;
343  if (i != j) ret(j, i) = ret(i, j);
344  }
345  }
346  } else if (particle_a == 1 && particle_b == 1) { // d2f/dlambda2
347  for (unsigned i = 0; i < N; i++) {
348  for (unsigned j = i; j < N; j++) {
349  double dist = std::abs(xlist[i][0] - xlist[j][0]);
350  double exponent = std::pow(dist / lambda_val_, alpha_);
351  double expterm = std::exp(-0.5 * exponent);
352  ret(i, j) = tau_val_ * tau_val_ * expterm * exponent /
353  (lambda_val_ * lambda_val_) * alpha_ / 2 *
354  (alpha_ / 2 * exponent - (alpha_ + 1));
355  if (i != j) ret(j, i) = ret(i, j);
356  }
357  }
358  } else { // d2f/d(tau)d(lambda)
359  for (unsigned i = 0; i < N; i++) {
360  for (unsigned j = i; j < N; j++) {
361  double dist = std::abs(xlist[i][0] - xlist[j][0]);
362  double exponent = std::pow(dist / lambda_val_, alpha_);
363  double expterm = std::exp(-0.5 * exponent);
364  ret(i, j) = tau_val_ * alpha_ * expterm / lambda_val_ * exponent;
365  if (i != j) ret(j, i) = ret(i, j);
366  }
367  }
368  }
369  return ret;
370  }
371
372  FloatsList get_second_derivative_matrix(unsigned particle_a,
373  unsigned particle_b,
374  const FloatsList& xlist, bool) const {
375  Eigen::MatrixXd mat(
376  get_second_derivative_matrix(particle_a, particle_b, xlist));
377  FloatsList ret;
378  for (int i = 0; i < mat.rows(); i++) {
379  Floats line;
380  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
381  ret.push_back(line);
382  }
383  return ret;
384  }
385
386  unsigned get_ndims_x1() const { return 1; }
387  unsigned get_ndims_x2() const { return 1; }
388  unsigned get_ndims_y() const { return 1; }
389
390  unsigned get_number_of_particles() const { return 2; }
391
392  bool get_particle_is_optimized(unsigned particle_no) const {
393  switch (particle_no) {
394  case 0: // tau
395  return Scale(tau_).get_nuisance_is_optimized();
396  case 1: // lambda
397  return Scale(lambda_).get_nuisance_is_optimized();
398  default:
399  IMP_THROW("Invalid particle number", ModelException);
400  }
401  }
402
404  unsigned count = 0;
405  if (Scale(tau_).get_nuisance_is_optimized()) count++;
406  if (Scale(lambda_).get_nuisance_is_optimized()) count++;
407  return count;
408  }
409
411  ModelObjectsTemp ret;
412  ret.push_back(tau_);
413  ret.push_back(lambda_);
414  return ret;
415  }
416
417  /*IMP_OBJECT_INLINE(Covariance1DFunction,
418  out << "covariance function with alpha = "
419  << alpha_ << std::endl, {});*/
421
422  private:
423  inline double get_value(double x1, double x2) const {
424  double dist = std::abs(x1 - x2);
425  double ret = dist / lambda_val_;
426  if (alpha_square_) {
427  ret *= ret;
428  } else {
429  ret = std::pow(ret, alpha_);
430  }
431  ret = IMP::square(tau_val_) * std::exp(-0.5 * ret);
432  if (do_jitter && dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
433  ret += IMP::square(tau_val_) * J_;
434  }
436  "function value is nan. tau = "
437  << tau_val_ << " lambda = " << lambda_val_
438  << " q1 = " << x1 << " q2 = " << x2);
439  return ret;
440  }
441
442  private:
443  double alpha_;
444  Pointer<Particle> tau_, lambda_;
445  double tau_val_, lambda_val_, J_, cutoff_, alpha_square_;
446  bool do_jitter;
447 };
448
449 IMPISD_END_NAMESPACE
450
451 #endif /* IMPISD_BIVARIATE_FUNCTIONS_H */
ModelObjectsTemp get_inputs() const
particle manipulation
bool has_changed() const
return true if internal parameters have changed.
unsigned get_number_of_optimized_particles() const
returns the number of particles that are optimized
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
A decorator for switching parameters particles.
void update()
update internal parameters
unsigned get_ndims_y() const
returns the number of output dimensions
A decorator for scale parameters particles.
A decorator for nuisance parameters particles.
Definition: Scale.h:24
#define IMP_LOG_VERBOSE(expr)
Definition: log_macros.h:83
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Ref counted objects should have private destructors.
#define IMP_LOG_TERSE(expr)
Definition: log_macros.h:72
Eigen::MatrixXd operator()(const IMP::FloatsList &xlist) const
evaluate the function at a list of points
#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
Common base class for heavy weight IMP objects.
Definition: Object.h:106
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const
update derivatives of particles
virtual Eigen::MatrixXd get_second_derivative_matrix(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const =0
return second derivative matrix
unsigned get_ndims_x1() const
returns the number of input dimensions
virtual void update()=0
update internal parameters
Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const
return derivative matrix
void add_to_derivatives(const Floats &x1, const Floats &x2, DerivativeAccumulator &accum) const
update derivatives of particles
FloatsList operator()(const IMP::FloatsList &xlist, bool) const
used for testing only
Base class for functions of two variables.
Eigen::MatrixXd get_second_derivative_matrix(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const
return second derivative matrix
bool get_particle_is_optimized(unsigned particle_no) const
returns true if the particle whose index is provided is optimized
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...
unsigned get_number_of_particles() const
returns the number of particles that this function uses
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
Definition: check_macros.h:50
A shared base class to help in debugging and things.
Object(std::string name)
Construct an object with the given name.
Floats operator()(const Floats &x1, const Floats &x2) const
evaluate the function at a certain point
Class to handle individual particles of a Model object.
Definition: Particle.h:41
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
An exception which is thrown when the Model has attributes with invalid values.
Definition: exception.h:189
Class for adding derivatives from restraints to the model.
virtual Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const =0
return derivative matrix