IMP  2.3.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);
131  lambda_val_ = Scale(ilambda).get_nuisance();
132  tau_val_ = Scale(tau).get_nuisance();
133  do_jitter = (jitter > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
134  alpha_square_ = (std::abs(alpha - 2) < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
135  update();
136  }
137 
138  bool has_changed() const {
139  double tmpt = Scale(tau_).get_nuisance();
140  double tmpl = Scale(lambda_).get_nuisance();
141  IMP_LOG_VERBOSE("Covariance1DFunction: has_changed(): ");
142  IMP_LOG_VERBOSE(tmpt << " " << tau_val_ << " ");
143  IMP_LOG_VERBOSE(tmpl << " " << lambda_val_ << " ");
144  if ((std::abs(tmpt - tau_val_) > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) ||
145  (std::abs(tmpl - lambda_val_) > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM)) {
146  IMP_LOG_VERBOSE("true" << std::endl);
147  return true;
148  } else {
149  IMP_LOG_VERBOSE("false" << std::endl);
150  return false;
151  }
152  }
153 
154  void update() {
155  lambda_val_ = Scale(lambda_).get_nuisance();
156  tau_val_ = Scale(tau_).get_nuisance();
157  IMP_LOG_TERSE("Covariance1DFunction: update() tau:= "
158  << tau_val_ << " lambda:=" << lambda_val_ << std::endl);
159  IMP_INTERNAL_CHECK(!base::isnan(tau_val_), "tau is nan.");
160  IMP_INTERNAL_CHECK(!base::isnan(lambda_val_), "lambda is nan.");
161  IMP_INTERNAL_CHECK(tau_val_ >= 0, "tau is negative.");
162  IMP_INTERNAL_CHECK(lambda_val_ >= 0, "lambda is negative.");
163  IMP_INTERNAL_CHECK(lambda_val_ != 0, "lambda is zero.");
164  }
165 
166  Floats operator()(const Floats& x1, const Floats& x2) const {
167  IMP_USAGE_CHECK(x1.size() == 1, "expecting a 1-D vector");
168  IMP_USAGE_CHECK(x2.size() == 1, "expecting a 1-D vector");
169  Floats ret(1, get_value(x1[0], x2[0]));
170  return ret;
171  }
172 
173  IMP_Eigen::MatrixXd operator()(const IMP::FloatsList& xlist) const {
174  const unsigned M = xlist.size();
175  IMP_Eigen::MatrixXd Mret(M, M);
176  for (unsigned i = 0; i < M; i++) {
177  for (unsigned j = i; j < M; j++) {
178  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
179  IMP_USAGE_CHECK(xlist[j].size() == 1, "expecting a 1-D vector");
180  double x1 = xlist[i][0];
181  double x2 = xlist[j][0];
182  double ret = get_value(x1, x2);
183  Mret(i, j) = ret;
184  if (i != j) Mret(j, i) = ret;
185  }
186  }
187  return Mret;
188  }
189 
190  FloatsList operator()(const IMP::FloatsList& xlist, bool) const {
191  IMP_Eigen::MatrixXd mat((*this)(xlist));
192  FloatsList ret;
193  for (unsigned i = 0; i < xlist.size(); i++)
194  for (unsigned j = 0; j < xlist.size(); j++)
195  ret.push_back(Floats(1, mat(i, j)));
196  return ret;
197  }
198 
199  void add_to_derivatives(const Floats& x1, const Floats& x2,
200  DerivativeAccumulator& accum) const {
201  // d[w(x1,x2)]/dtau = 2/tau*(w(x1,x2))
202  double val = get_value(x1[0], x2[0]);
203  double tauderiv = 2. / tau_val_ * val;
204  IMP_INTERNAL_CHECK(!base::isnan(tauderiv), "tau derivative is nan.");
205  Scale(tau_).add_to_nuisance_derivative(tauderiv, accum);
206  // d[w(x,x')]/dlambda
207  // = w(x,x') ( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
208  double lambdaderiv =
209  val *
210  (alpha_ * std::pow((std::abs(x1[0] - x2[0]) / lambda_val_), alpha_) /
211  (2. * lambda_val_));
212  IMP_INTERNAL_CHECK(!base::isnan(lambdaderiv), "lambda derivative is nan.");
213  Scale(lambda_).add_to_nuisance_derivative(lambdaderiv, accum);
214  }
215 
216  void add_to_particle_derivative(unsigned particle_no, double value,
217  DerivativeAccumulator& accum) const {
218  switch (particle_no) {
219  case 0: // tau
220  IMP_INTERNAL_CHECK(!base::isnan(value), "tau derivative is nan.");
221  Scale(tau_).add_to_nuisance_derivative(value, accum);
222  break;
223  case 1: // lambda
224  IMP_INTERNAL_CHECK(!base::isnan(value), "lambda derivative is nan.");
225  Scale(lambda_).add_to_nuisance_derivative(value, accum);
226  break;
227  default:
228  IMP_THROW("Invalid particle number", ModelException);
229  }
230  }
231 
232  IMP_Eigen::MatrixXd get_derivative_matrix(unsigned particle_no,
233  const FloatsList& xlist) const {
234  // Strategy: fill in the main diagonal, then fill with zeros
235  // if the value of the function falls below cutoff.
236  // assumes data points are ordered!
237  unsigned N = xlist.size();
238  IMP_Eigen::MatrixXd ret(N, N);
239  double diag;
240  switch (particle_no) {
241  case 0: // tau
242  // d[w(x1,x1)]/dtau
243  // = 2/tau*w(x1,x1)
244  diag = get_value(xlist[0][0], xlist[0][0]);
245  diag *= 2. / tau_val_;
246  break;
247  case 1: // lambda
248  // d[w(x,x)]/dlambda
249  //= w(x,x)
250  // *( alpha /(2 lambda^{alpha+1}))
251  diag = 0;
252  break;
253  default:
254  IMP_THROW("Invalid particle number", ModelException);
255  }
257  "derivative matrix is nan on the diagonal.");
258  for (unsigned i = 0; i < N; i++) ret(i, i) = diag;
259  //
260  bool initial_loop = true;
261  double abs_cutoff = cutoff_ * diag;
262  double dmax = -1;
263  for (unsigned i = 0; i < N; i++) {
264  for (unsigned j = i + 1; j < N; j++) {
265  double x1(xlist[i][0]), x2(xlist[j][0]);
266  double val;
267  double dist(std::abs(x1 - x2));
268  // compute all entries as long as the cutoff distance was
269  // not recorded (initial_loop) or as long as the distance is
270  // smaller than the cutoff distance
271  if (initial_loop || dist <= dmax) {
272  switch (particle_no) {
273  case 0: // tau
274  // d[w(x1,x2)]/dtau
275  // = 2/tau*w(x1,x2)
276  val = get_value(xlist[i][0], xlist[j][0]);
277  val *= 2. / tau_val_;
278  break;
279  case 1: // lambda
280  // d[w(x,x')]/dlambda
281  //= w(x,x')
282  // *( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
283  if (dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
284  val = 0;
285  } else {
286  val = get_value(xlist[i][0], xlist[j][0]);
287  val *= alpha_ * std::pow((dist / lambda_val_), alpha_) /
288  (2. * lambda_val_);
289  }
290  break;
291  default:
292  IMP_THROW("Invalid particle number", ModelException);
293  }
294  // the value has been computed and is in val
295  // now check if it is smaller than the cutoff.
296  // If true change the flag and update the distance
297  if (std::abs(val) <= abs_cutoff) {
298  if (initial_loop) {
299  initial_loop = false;
300  dmax = dist;
301  } else if (dist < dmax) {
302  dmax = dist;
303  }
304  }
305  } else { // e.g. initial_loop == false && dist > dmax
306  val = 0;
307  }
309  "derivative matrix is nan at position("
310  << i << "," << j << ").");
311  ret(i, j) = val;
312  ret(j, i) = val;
313  }
314  }
315  return ret;
316  }
317 
318  FloatsList get_derivative_matrix(unsigned particle_no,
319  const FloatsList& xlist, bool) const {
320  IMP_Eigen::MatrixXd mat(get_derivative_matrix(particle_no, xlist));
321  FloatsList ret;
322  for (int i = 0; i < mat.rows(); i++) {
323  Floats line;
324  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
325  ret.push_back(line);
326  }
327  return ret;
328  }
329 
330  IMP_Eigen::MatrixXd get_second_derivative_matrix(
331  unsigned particle_a, unsigned particle_b, const FloatsList& xlist) const {
332  unsigned N(xlist.size());
333  IMP_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 
373  FloatsList get_second_derivative_matrix(unsigned particle_a,
374  unsigned particle_b,
375  const FloatsList& xlist, bool) const {
376  IMP_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 { return 1; }
388  unsigned get_ndims_x2() const { return 1; }
389  unsigned get_ndims_y() const { return 1; }
390 
391  unsigned get_number_of_particles() const { return 2; }
392 
393  bool get_particle_is_optimized(unsigned particle_no) const {
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 
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 
413  ret.push_back(tau_);
414  ret.push_back(lambda_);
415  return ret;
416  }
417 
419  ContainersTemp ret;
420  return ret;
421  }
422 
423  /*IMP_OBJECT_INLINE(Covariance1DFunction,
424  out << "covariance function with alpha = "
425  << alpha_ << std::endl, {});*/
427 
428  private:
429  inline double get_value(double x1, double x2) const {
430  double dist = std::abs(x1 - x2);
431  double ret = dist / lambda_val_;
432  if (alpha_square_) {
433  ret *= ret;
434  } else {
435  ret = std::pow(ret, alpha_);
436  }
437  ret = IMP::square(tau_val_) * std::exp(-0.5 * ret);
438  if (do_jitter && dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
439  ret += IMP::square(tau_val_) * J_;
440  }
442  "function value is nan. tau = "
443  << tau_val_ << " lambda = " << lambda_val_
444  << " q1 = " << x1 << " q2 = " << x2);
445  return ret;
446  }
447 
448  private:
449  double alpha_;
450  base::Pointer<kernel::Particle> tau_, lambda_;
451  double tau_val_, lambda_val_, J_, cutoff_, alpha_square_;
452  bool do_jitter;
453 };
454 
455 IMPISD_END_NAMESPACE
456 
457 #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
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
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.
ParticlesTemp get_input_particles(const ModelObjectsTemp &mos)
Return all the input particles for a given ModelObject.
A decorator for nuisance parameters particles.
Add scale parameter to particle.
Definition: Scale.h:24
#define IMP_LOG_VERBOSE(expr)
Definition: log_macros.h:94
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Ref counted objects should have private destructors.
virtual IMP_Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const =0
return derivative matrix
#define IMP_LOG_TERSE(expr)
Definition: log_macros.h:83
#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:141
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)
Return all the input particles for a given ModelObject.
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
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
Base class for functions of two variables.
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: Object.h:106
Classes to handle individual model particles. (Note that implementation of inline functions in in int...
bool isnan(const T &a)
Return true if a number is NaN.
Definition: 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.
Definition: check_macros.h:50
A shared base class to help in debugging and things.
Floats operator()(const Floats &x1, const Floats &x2) const
evaluate the function at a certain point
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:170
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