IMP  2.2.1
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-2014 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/kernel/Particle.h>
13 #include <IMP/isd/Nuisance.h>
14 #include <IMP/isd/Scale.h>
15 #include <IMP/isd/Switching.h>
16 #include <IMP/base/Object.h>
17 #include <IMP/algebra/eigen3/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 base::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 IMP_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 IMP_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 IMP_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 kernel::ParticlesTemp get_input_particles() const = 0;
102  virtual ContainersTemp get_input_containers() const = 0;
103 
105 };
106 //
107 //! Covariance function
108 /* \f[w(x,x') = \tau^2 \exp\left(-\frac{|x-x'|^\alpha}{2\lambda^\alpha} +
109  * \delta_{ij} J\f]
110  * \param[in] \f$\tau\f$ ISD Scale
111  * \param[in] \f$\lambda\f$ ISD Scale
112  * \param[in] \f$\alpha\f$ positive double, usually greater than 1.
113  * Default is 2.
114  * \param[in] J is some jitter. Try J=0.01 if you get NANs.
115  * \param[in] cutoff is a positive double indicating when to consider that
116  * values are zero (to avoid computing exponentials). cutoff is relative to the
117  * value when x=x', and affects only the call get_derivative_matrix.
118  */
119 class IMPISDEXPORT Covariance1DFunction : public BivariateFunction {
120  public:
122  double alpha = 2.0, double jitter = 0.0,
123  double cutoff = 1e-7)
124  : BivariateFunction("Covariance1DFunction %1%"),
125  alpha_(alpha),
126  tau_(tau),
127  lambda_(ilambda),
128  J_(jitter),
129  cutoff_(cutoff) {
130  IMP_LOG_TERSE("Covariance1DFunction: constructor" << std::endl);
133  lambda_val_ = Scale(ilambda).get_nuisance();
134  tau_val_ = Scale(tau).get_nuisance();
135  do_jitter = (jitter > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
136  alpha_square_ = (std::abs(alpha - 2) < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
137  update();
138  }
139 
140  bool has_changed() const {
141  double tmpt = Scale(tau_).get_nuisance();
142  double tmpl = Scale(lambda_).get_nuisance();
143  IMP_LOG_VERBOSE("Covariance1DFunction: has_changed(): ");
144  IMP_LOG_VERBOSE(tmpt << " " << tau_val_ << " ");
145  IMP_LOG_VERBOSE(tmpl << " " << lambda_val_ << " ");
146  if ((std::abs(tmpt - tau_val_) > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) ||
147  (std::abs(tmpl - lambda_val_) > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM)) {
148  IMP_LOG_VERBOSE("true" << std::endl);
149  return true;
150  } else {
151  IMP_LOG_VERBOSE("false" << std::endl);
152  return false;
153  }
154  }
155 
156  void update() {
157  lambda_val_ = Scale(lambda_).get_nuisance();
158  tau_val_ = Scale(tau_).get_nuisance();
159  IMP_LOG_TERSE("Covariance1DFunction: update() tau:= "
160  << tau_val_ << " lambda:=" << lambda_val_ << std::endl);
161  IMP_INTERNAL_CHECK(!base::isnan(tau_val_), "tau is nan.");
162  IMP_INTERNAL_CHECK(!base::isnan(lambda_val_), "lambda is nan.");
163  IMP_INTERNAL_CHECK(tau_val_ >= 0, "tau is negative.");
164  IMP_INTERNAL_CHECK(lambda_val_ >= 0, "lambda is negative.");
165  IMP_INTERNAL_CHECK(lambda_val_ != 0, "lambda is zero.");
166  }
167 
168  Floats operator()(const Floats& x1, const Floats& x2) const {
169  IMP_USAGE_CHECK(x1.size() == 1, "expecting a 1-D vector");
170  IMP_USAGE_CHECK(x2.size() == 1, "expecting a 1-D vector");
171  Floats ret(1, get_value(x1[0], x2[0]));
172  return ret;
173  }
174 
175  IMP_Eigen::MatrixXd operator()(const IMP::FloatsList& xlist) const {
176  const unsigned M = xlist.size();
177  IMP_Eigen::MatrixXd Mret(M, M);
178  for (unsigned i = 0; i < M; i++) {
179  for (unsigned j = i; j < M; j++) {
180  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
181  IMP_USAGE_CHECK(xlist[j].size() == 1, "expecting a 1-D vector");
182  double x1 = xlist[i][0];
183  double x2 = xlist[j][0];
184  double ret = get_value(x1, x2);
185  Mret(i, j) = ret;
186  if (i != j) Mret(j, i) = ret;
187  }
188  }
189  return Mret;
190  }
191 
192  FloatsList operator()(const IMP::FloatsList& xlist, bool) const {
193  IMP_Eigen::MatrixXd mat((*this)(xlist));
194  FloatsList ret;
195  for (unsigned i = 0; i < xlist.size(); i++)
196  for (unsigned j = 0; j < xlist.size(); j++)
197  ret.push_back(Floats(1, mat(i, j)));
198  return ret;
199  }
200 
201  void add_to_derivatives(const Floats& x1, const Floats& x2,
202  DerivativeAccumulator& accum) const {
203  // d[w(x1,x2)]/dtau = 2/tau*(w(x1,x2))
204  double val = get_value(x1[0], x2[0]);
205  double tauderiv = 2. / tau_val_ * val;
206  IMP_INTERNAL_CHECK(!base::isnan(tauderiv), "tau derivative is nan.");
207  Scale(tau_).add_to_nuisance_derivative(tauderiv, accum);
208  // d[w(x,x')]/dlambda
209  // = w(x,x') ( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
210  double lambdaderiv =
211  val *
212  (alpha_ * std::pow((std::abs(x1[0] - x2[0]) / lambda_val_), alpha_) /
213  (2. * lambda_val_));
214  IMP_INTERNAL_CHECK(!base::isnan(lambdaderiv), "lambda derivative is nan.");
215  Scale(lambda_).add_to_nuisance_derivative(lambdaderiv, accum);
216  }
217 
218  void add_to_particle_derivative(unsigned particle_no, double value,
219  DerivativeAccumulator& accum) const {
220  switch (particle_no) {
221  case 0: // tau
222  IMP_INTERNAL_CHECK(!base::isnan(value), "tau derivative is nan.");
223  Scale(tau_).add_to_nuisance_derivative(value, accum);
224  break;
225  case 1: // lambda
226  IMP_INTERNAL_CHECK(!base::isnan(value), "lambda derivative is nan.");
227  Scale(lambda_).add_to_nuisance_derivative(value, accum);
228  break;
229  default:
230  IMP_THROW("Invalid particle number", ModelException);
231  }
232  }
233 
234  IMP_Eigen::MatrixXd get_derivative_matrix(unsigned particle_no,
235  const FloatsList& xlist) const {
236  // Strategy: fill in the main diagonal, then fill with zeros
237  // if the value of the function falls below cutoff.
238  // assumes data points are ordered!
239  unsigned N = xlist.size();
240  IMP_Eigen::MatrixXd ret(N, N);
241  double diag;
242  switch (particle_no) {
243  case 0: // tau
244  // d[w(x1,x1)]/dtau
245  // = 2/tau*w(x1,x1)
246  diag = get_value(xlist[0][0], xlist[0][0]);
247  diag *= 2. / tau_val_;
248  break;
249  case 1: // lambda
250  // d[w(x,x)]/dlambda
251  //= w(x,x)
252  // *( alpha /(2 lambda^{alpha+1}))
253  diag = 0;
254  break;
255  default:
256  IMP_THROW("Invalid particle number", ModelException);
257  }
259  "derivative matrix is nan on the diagonal.");
260  for (unsigned i = 0; i < N; i++) ret(i, i) = diag;
261  //
262  bool initial_loop = true;
263  double abs_cutoff = cutoff_ * diag;
264  double dmax = -1;
265  for (unsigned i = 0; i < N; i++) {
266  for (unsigned j = i + 1; j < N; j++) {
267  double x1(xlist[i][0]), x2(xlist[j][0]);
268  double val;
269  double dist(std::abs(x1 - x2));
270  // compute all entries as long as the cutoff distance was
271  // not recorded (initial_loop) or as long as the distance is
272  // smaller than the cutoff distance
273  if (initial_loop || dist <= dmax) {
274  switch (particle_no) {
275  case 0: // tau
276  // d[w(x1,x2)]/dtau
277  // = 2/tau*w(x1,x2)
278  val = get_value(xlist[i][0], xlist[j][0]);
279  val *= 2. / tau_val_;
280  break;
281  case 1: // lambda
282  // d[w(x,x')]/dlambda
283  //= w(x,x')
284  // *( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
285  if (dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
286  val = 0;
287  } else {
288  val = get_value(xlist[i][0], xlist[j][0]);
289  val *= alpha_ * std::pow((dist / lambda_val_), alpha_) /
290  (2. * lambda_val_);
291  }
292  break;
293  default:
294  IMP_THROW("Invalid particle number", ModelException);
295  }
296  // the value has been computed and is in val
297  // now check if it is smaller than the cutoff.
298  // If true change the flag and update the distance
299  if (std::abs(val) <= abs_cutoff) {
300  if (initial_loop) {
301  initial_loop = false;
302  dmax = dist;
303  } else if (dist < dmax) {
304  dmax = dist;
305  }
306  }
307  } else { // e.g. initial_loop == false && dist > dmax
308  val = 0;
309  }
311  "derivative matrix is nan at position("
312  << i << "," << j << ").");
313  ret(i, j) = val;
314  ret(j, i) = val;
315  }
316  }
317  return ret;
318  }
319 
320  FloatsList get_derivative_matrix(unsigned particle_no,
321  const FloatsList& xlist, bool) const {
322  IMP_Eigen::MatrixXd mat(get_derivative_matrix(particle_no, xlist));
323  FloatsList ret;
324  for (int i = 0; i < mat.rows(); i++) {
325  Floats line;
326  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
327  ret.push_back(line);
328  }
329  return ret;
330  }
331 
332  IMP_Eigen::MatrixXd get_second_derivative_matrix(
333  unsigned particle_a, unsigned particle_b, const FloatsList& xlist) const {
334  unsigned N(xlist.size());
335  IMP_Eigen::MatrixXd ret(N, N);
336  if (particle_a > 1) IMP_THROW("Invalid particle 1 number", ModelException);
337  if (particle_b > 1) IMP_THROW("Invalid particle 2 number", ModelException);
338  // d2f/dtau2
339  if (particle_a == 0 && particle_b == 0) {
340  for (unsigned i = 0; i < N; i++) {
341  for (unsigned j = i; j < N; j++) {
342  double dist = std::abs(xlist[i][0] - xlist[j][0]);
343  double exponent = std::pow(dist / lambda_val_, alpha_);
344  double expterm = std::exp(-0.5 * exponent);
345  ret(i, j) = 2 * expterm;
346  if (i != j) ret(j, i) = ret(i, j);
347  }
348  }
349  } else if (particle_a == 1 && particle_b == 1) { // d2f/dlambda2
350  for (unsigned i = 0; i < N; i++) {
351  for (unsigned j = i; j < N; j++) {
352  double dist = std::abs(xlist[i][0] - xlist[j][0]);
353  double exponent = std::pow(dist / lambda_val_, alpha_);
354  double expterm = std::exp(-0.5 * exponent);
355  ret(i, j) = tau_val_ * tau_val_ * expterm * exponent /
356  (lambda_val_ * lambda_val_) * alpha_ / 2 *
357  (alpha_ / 2 * exponent - (alpha_ + 1));
358  if (i != j) ret(j, i) = ret(i, j);
359  }
360  }
361  } else { // d2f/d(tau)d(lambda)
362  for (unsigned i = 0; i < N; i++) {
363  for (unsigned j = i; j < N; j++) {
364  double dist = std::abs(xlist[i][0] - xlist[j][0]);
365  double exponent = std::pow(dist / lambda_val_, alpha_);
366  double expterm = std::exp(-0.5 * exponent);
367  ret(i, j) = tau_val_ * alpha_ * expterm / lambda_val_ * exponent;
368  if (i != j) ret(j, i) = ret(i, j);
369  }
370  }
371  }
372  return ret;
373  }
374 
375  FloatsList get_second_derivative_matrix(unsigned particle_a,
376  unsigned particle_b,
377  const FloatsList& xlist, bool) const {
378  IMP_Eigen::MatrixXd mat(
379  get_second_derivative_matrix(particle_a, particle_b, xlist));
380  FloatsList ret;
381  for (int i = 0; i < mat.rows(); i++) {
382  Floats line;
383  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
384  ret.push_back(line);
385  }
386  return ret;
387  }
388 
389  unsigned get_ndims_x1() const { return 1; }
390  unsigned get_ndims_x2() const { return 1; }
391  unsigned get_ndims_y() const { return 1; }
392 
393  unsigned get_number_of_particles() const { return 2; }
394 
395  bool get_particle_is_optimized(unsigned particle_no) const {
396  switch (particle_no) {
397  case 0: // tau
398  return Scale(tau_).get_nuisance_is_optimized();
399  case 1: // lambda
400  return Scale(lambda_).get_nuisance_is_optimized();
401  default:
402  IMP_THROW("Invalid particle number", ModelException);
403  }
404  }
405 
407  unsigned count = 0;
408  if (Scale(tau_).get_nuisance_is_optimized()) count++;
409  if (Scale(lambda_).get_nuisance_is_optimized()) count++;
410  return count;
411  }
412 
415  ret.push_back(tau_);
416  ret.push_back(lambda_);
417  return ret;
418  }
419 
421  ContainersTemp ret;
422  return ret;
423  }
424 
425  /*IMP_OBJECT_INLINE(Covariance1DFunction,
426  out << "covariance function with alpha = "
427  << alpha_ << std::endl, {});*/
429 
430  private:
431  inline double get_value(double x1, double x2) const {
432  double dist = std::abs(x1 - x2);
433  double ret = dist / lambda_val_;
434  if (alpha_square_) {
435  ret *= ret;
436  } else {
437  ret = std::pow(ret, alpha_);
438  }
439  ret = IMP::square(tau_val_) * std::exp(-0.5 * ret);
440  if (do_jitter && dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
441  ret += IMP::square(tau_val_) * J_;
442  }
444  "function value is nan. tau = "
445  << tau_val_ << " lambda = " << lambda_val_
446  << " q1 = " << x1 << " q2 = " << x2);
447  return ret;
448  }
449 
450  private:
451  double alpha_;
452  base::Pointer<kernel::Particle> tau_, lambda_;
453  double tau_val_, lambda_val_, J_, cutoff_, alpha_square_;
454  bool do_jitter;
455 };
456 
457 IMPISD_END_NAMESPACE
458 
459 #endif /* IMPISD_BIVARIATE_FUNCTIONS_H */
Class for adding derivatives from restraints to the model.
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
A decorator for switching parameters particles.
IMP_Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const
return derivative matrix
void update()
update internal parameters
unsigned get_ndims_y() const
returns the number of output dimensions
A decorator for scale parameters particles.
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Ref counted objects should have private destructors.
ParticlesTemp get_input_particles(const ModelObjectsTemp &mos)
A decorator for nuisance parameters particles.
Add scale parameter to particle.
Definition: Scale.h:24
virtual IMP_Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const =0
return derivative matrix
static Scale decorate_particle(::IMP::kernel::Particle *p)
Definition: Scale.h:29
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
#define IMP_LOG_TERSE(expr)
IMP_Eigen::MatrixXd get_second_derivative_matrix(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const
return second derivative matrix
ContainersTemp get_input_containers(const ModelObjectsTemp &mos)
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const
update derivatives of particles
unsigned get_ndims_x1() const
returns the number of input dimensions
Object(std::string name)
Construct an object with the given name.
virtual void update()=0
update internal parameters
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
IMP::base::Vector< IMP::base::WeakPointer< Container > > ContainersTemp
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
IMP_Eigen::MatrixXd operator()(const IMP::FloatsList &xlist) const
evaluate the function at a list of points
#define IMP_IF_CHECK(level)
Execute the code block if a certain level checks are on.
Base class for functions of two variables.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
bool get_particle_is_optimized(unsigned particle_no) const
returns true if the particle whose index is provided is optimized
Class to handle individual model particles.
Common base class for heavy weight IMP objects.
Definition: base/Object.h:106
Classes to handle individual model particles.
bool isnan(const T &a)
Return true if a number is NaN.
Definition: base/math.h:24
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.
Floats operator()(const Floats &x1, const Floats &x2) const
evaluate the function at a certain point
A shared base class to help in debugging and things.
kernel::ParticlesTemp get_input_particles() const
particle manipulation
virtual IMP_Eigen::MatrixXd get_second_derivative_matrix(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const =0
return second derivative matrix
#define IMP_LOG_VERBOSE(expr)