IMP  2.1.0
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-2013 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 <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 {
26  public:
27 
28  BivariateFunction(std::string str) : Object(str) {}
29 
30  //! evaluate the function at a certain point
31  virtual Floats operator()
32  (const Floats& x1,
33  const Floats& x2) const = 0;
34 
35  //! evaluate the function at a list of points
36  /* returns a NxN matrix of entries where N=xlist.rows()
37  * and entry ij of the matrix is f(xlist(i),xlist(j))
38  */
39  virtual Eigen::MatrixXd operator()
40  (const IMP::FloatsList& xlist) const = 0;
41 
42  //! used for testing only
43  virtual FloatsList operator() (const IMP::FloatsList& xlist,
44  bool stupid) const = 0;
45 
46  //! return true if internal parameters have changed.
47  virtual bool has_changed() const = 0;
48 
49  //! update internal parameters
50  virtual void update() = 0;
51 
52  //! update derivatives of particles
53  virtual void add_to_derivatives(
54  const Floats& x1,
55  const Floats& x2,
56  DerivativeAccumulator &accum) const = 0;
57 
58  //! update derivatives of particles
59  /* add to the given particle the specified derivative
60  * guarantees that the particle_no (starting at 0) matches with
61  * the columns of get_derivative_matrix.
62  */
63  virtual void add_to_particle_derivative(unsigned particle_no,
64  double value, DerivativeAccumulator &accum) const = 0;
65 
66  //! return derivative matrix
67  /* m_ij = d(func(xlist[i],xlist[j]))/dparticle_no
68  * the matrix is NxN where N = xlist.size()
69  */
70  virtual Eigen::MatrixXd get_derivative_matrix(
71  unsigned particle_no,
72  const FloatsList& xlist) const = 0;
73 
74  //for testing purposes
75  virtual FloatsList get_derivative_matrix(
76  unsigned particle_no,
77  const FloatsList& xlist,
78  bool stupid) const = 0;
79 
80  //! return second derivative matrix
81  virtual Eigen::MatrixXd get_second_derivative_matrix(
82  unsigned particle_a, unsigned particle_b,
83  const FloatsList& xlist) const = 0;
84 
85  //for testing purposes
86  virtual FloatsList get_second_derivative_matrix(
87  unsigned particle_a, unsigned particle_b,
88  const FloatsList& xlist,
89  bool stupid) const = 0;
90 
91  //! returns the number of input dimensions
92  virtual unsigned get_ndims_x1() const = 0;
93  virtual unsigned get_ndims_x2() const = 0;
94 
95  //! returns the number of output dimensions
96  virtual unsigned get_ndims_y() const = 0;
97 
98  //! returns the number of particles that this function uses
99  virtual unsigned get_number_of_particles() const = 0;
100 
101  //! returns true if the particle whose index is provided is optimized
102  virtual bool get_particle_is_optimized(unsigned particle_no) const = 0;
103 
104  //! returns the number of particles that are optimized
105  virtual unsigned get_number_of_optimized_particles() const = 0;
106 
107  //! particle manipulation
108  virtual kernel::ParticlesTemp get_input_particles() const = 0;
109  virtual ContainersTemp get_input_containers() const = 0;
110 
112 };
113 //
114 //! Covariance function
115 /* \f[w(x,x') = \tau^2 \exp\left(-\frac{|x-x'|^\alpha}{2\lambda^\alpha} +
116  * \delta_{ij} J\f]
117  * \param[in] \f$\tau\f$ ISD Scale
118  * \param[in] \f$\lambda\f$ ISD Scale
119  * \param[in] \f$\alpha\f$ positive double, usually greater than 1.
120  * Default is 2.
121  * \param[in] J is some jitter. Try J=0.01 if you get NANs.
122  * \param[in] cutoff is a positive double indicating when to consider that
123  * values are zero (to avoid computing exponentials). cutoff is relative to the
124  * value when x=x', and affects only the call get_derivative_matrix.
125  */
126 class IMPISDEXPORT Covariance1DFunction : public BivariateFunction
127 {
128  public:
130  double alpha=2.0, double jitter =0.0, double cutoff=1e-7) :
131  BivariateFunction("Covariance1DFunction %1%"), alpha_(alpha),
132  tau_(tau), lambda_(ilambda), J_(jitter),
133  cutoff_(cutoff)
134  {
135  IMP_LOG_TERSE( "Covariance1DFunction: constructor" << std::endl);
138  lambda_val_= Scale(ilambda).get_nuisance();
139  tau_val_= Scale(tau).get_nuisance();
140  do_jitter = (jitter>IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
141  alpha_square_ = (std::abs(alpha-2) <
142  IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
143  update();
144  }
145 
146  bool has_changed() const {
147  double tmpt = Scale(tau_).get_nuisance();
148  double tmpl = Scale(lambda_).get_nuisance();
149  IMP_LOG_VERBOSE( "Covariance1DFunction: has_changed(): ");
150  IMP_LOG_VERBOSE( tmpt << " " << tau_val_ << " " );
151  IMP_LOG_VERBOSE( tmpl << " " << lambda_val_ << " " );
152  if ((std::abs(tmpt - tau_val_) >
153  IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM)
154  || (std::abs(tmpl - lambda_val_) >
155  IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM))
156  {
157  IMP_LOG_VERBOSE( "true" << std::endl);
158  return true;
159  } else {
160  IMP_LOG_VERBOSE( "false" << std::endl);
161  return false;
162  }
163  }
164 
165  void update() {
166  lambda_val_= Scale(lambda_).get_nuisance();
167  tau_val_= Scale(tau_).get_nuisance();
168  IMP_LOG_TERSE( "Covariance1DFunction: update() tau:= "
169  << tau_val_ << " lambda:=" << lambda_val_
170  << std::endl);
171  IMP_INTERNAL_CHECK(!base::isnan(tau_val_),
172  "tau is nan.");
173  IMP_INTERNAL_CHECK(!base::isnan(lambda_val_),
174  "lambda is nan.");
175  IMP_INTERNAL_CHECK(tau_val_ >= 0, "tau is negative.");
176  IMP_INTERNAL_CHECK(lambda_val_ >= 0, "lambda is negative.");
177  IMP_INTERNAL_CHECK(lambda_val_ != 0, "lambda is zero.");
178  }
179 
181  const Floats& x2) const
182  {
183  IMP_USAGE_CHECK(x1.size() == 1, "expecting a 1-D vector");
184  IMP_USAGE_CHECK(x2.size() == 1, "expecting a 1-D vector");
185  Floats ret(1, get_value(x1[0],x2[0]));
186  return ret;
187  }
188 
189  Eigen::MatrixXd operator()(const IMP::FloatsList& xlist) const
190  {
191  const unsigned M=xlist.size();
192  Eigen::MatrixXd Mret(M,M);
193  for (unsigned i=0; i<M; i++)
194  {
195  for (unsigned j=i; j<M; j++)
196  {
197  IMP_USAGE_CHECK(xlist[i].size() == 1,
198  "expecting a 1-D vector");
199  IMP_USAGE_CHECK(xlist[j].size() == 1,
200  "expecting a 1-D vector");
201  double x1 = xlist[i][0];
202  double x2 = xlist[j][0];
203  double ret = get_value(x1,x2);
204  Mret(i,j) = ret;
205  if (i != j) Mret(j,i) = ret;
206  }
207  }
208  return Mret;
209  }
210 
211  FloatsList operator()(const IMP::FloatsList& xlist, bool) const
212  {
213  Eigen::MatrixXd mat((*this)(xlist));
214  FloatsList ret;
215  for (unsigned i=0; i<xlist.size(); i++)
216  for (unsigned j=0; j<xlist.size(); j++)
217  ret.push_back(Floats(1,mat(i,j)));
218  return ret;
219  }
220 
221  void add_to_derivatives(const Floats& x1,
222  const Floats& x2,
223  DerivativeAccumulator &accum) const
224  {
225  //d[w(x1,x2)]/dtau = 2/tau*(w(x1,x2))
226  double val = get_value(x1[0],x2[0]);
227  double tauderiv = 2./tau_val_ * val;
228  IMP_INTERNAL_CHECK(!base::isnan(tauderiv),
229  "tau derivative is nan.");
230  Scale(tau_).add_to_nuisance_derivative(tauderiv, accum);
231  //d[w(x,x')]/dlambda
232  // = w(x,x') ( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
233  double lambdaderiv =
234  val * (alpha_ *
235  std::pow((std::abs(x1[0]-x2[0])/lambda_val_),alpha_)
236  /(2.*lambda_val_));
237  IMP_INTERNAL_CHECK(!base::isnan(lambdaderiv),
238  "lambda derivative is nan.");
239  Scale(lambda_).add_to_nuisance_derivative(lambdaderiv, accum);
240  }
241 
242  void add_to_particle_derivative(unsigned particle_no,
243  double value, DerivativeAccumulator &accum) const
244  {
245  switch (particle_no)
246  {
247  case 0: //tau
249  "tau derivative is nan.");
250  Scale(tau_).add_to_nuisance_derivative(value, accum);
251  break;
252  case 1: //lambda
254  "lambda derivative is nan.");
255  Scale(lambda_).add_to_nuisance_derivative(value, accum);
256  break;
257  default:
258  IMP_THROW("Invalid particle number", ModelException);
259  }
260  }
261 
262  Eigen::MatrixXd get_derivative_matrix(
263  unsigned particle_no,
264  const FloatsList& xlist) const
265  {
266  // Strategy: fill in the main diagonal, then fill with zeros
267  // if the value of the function falls below cutoff.
268  // assumes data points are ordered!
269  unsigned N=xlist.size();
270  Eigen::MatrixXd ret(N,N);
271  double diag;
272  switch (particle_no)
273  {
274  case 0: //tau
275  //d[w(x1,x1)]/dtau
276  // = 2/tau*w(x1,x1)
277  diag = get_value(xlist[0][0],xlist[0][0]);
278  diag *= 2./tau_val_;
279  break;
280  case 1: //lambda
281  //d[w(x,x)]/dlambda
282  //= w(x,x)
283  // *( alpha /(2 lambda^{alpha+1}))
284  diag = 0;
285  break;
286  default:
287  IMP_THROW("Invalid particle number",
288  ModelException);
289  }
291  "derivative matrix is nan on the diagonal.");
292  for (unsigned i=0; i<N; i++) ret(i,i) = diag;
293  //
294  bool initial_loop=true;
295  double abs_cutoff = cutoff_*diag;
296  double dmax=-1;
297  for (unsigned i=0; i<N; i++)
298  {
299  for (unsigned j=i+1; j<N; j++)
300  {
301  double x1(xlist[i][0]), x2(xlist[j][0]);
302  double val;
303  double dist(std::abs(x1-x2));
304  //compute all entries as long as the cutoff distance was
305  //not recorded (initial_loop) or as long as the distance is
306  //smaller than the cutoff distance
307  if (initial_loop || dist <= dmax)
308  {
309  switch (particle_no)
310  {
311  case 0: //tau
312  //d[w(x1,x2)]/dtau
313  // = 2/tau*w(x1,x2)
314  val = get_value(xlist[i][0],xlist[j][0]);
315  val *= 2./tau_val_;
316  break;
317  case 1: //lambda
318  //d[w(x,x')]/dlambda
319  //= w(x,x')
320  // *( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
321  if (dist<IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
322  val = 0;
323  } else {
324  val = get_value(xlist[i][0],xlist[j][0]);
325  val *= alpha_ *
326  std::pow(
327  (dist/lambda_val_), alpha_)
328  /(2.*lambda_val_);
329  }
330  break;
331  default:
332  IMP_THROW("Invalid particle number",
333  ModelException);
334  }
335  // the value has been computed and is in val
336  //now check if it is smaller than the cutoff.
337  //If true change the flag and update the distance
338  if (std::abs(val) <= abs_cutoff)
339  {
340  if (initial_loop)
341  {
342  initial_loop = false;
343  dmax = dist;
344  } else if (dist < dmax)
345  {
346  dmax = dist;
347  }
348  }
349  } else { // e.g. initial_loop == false && dist > dmax
350  val = 0;
351  }
353  "derivative matrix is nan at position("
354  << i << "," << j << ").");
355  ret(i,j) = val;
356  ret(j,i) = val;
357  }
358  }
359  return ret;
360  }
361 
363  unsigned particle_no,
364  const FloatsList& xlist, bool) const
365  {
366  Eigen::MatrixXd mat(get_derivative_matrix(particle_no, xlist));
367  FloatsList ret;
368  for (int i=0; i<mat.rows(); i++)
369  {
370  Floats line;
371  for (int j=0; j<mat.cols(); j++)
372  line.push_back(mat(i,j));
373  ret.push_back(line);
374  }
375  return ret;
376  }
377 
379  unsigned particle_a, unsigned particle_b,
380  const FloatsList& xlist) const
381  {
382  unsigned N(xlist.size());
383  Eigen::MatrixXd ret(N,N);
384  if (particle_a > 1)
385  IMP_THROW("Invalid particle 1 number", ModelException);
386  if (particle_b > 1)
387  IMP_THROW("Invalid particle 2 number", ModelException);
388  //d2f/dtau2
389  if (particle_a == 0 && particle_b == 0)
390  {
391  for (unsigned i=0; i<N; i++)
392  {
393  for (unsigned j=i; j<N; j++)
394  {
395  double dist = std::abs(xlist[i][0]-xlist[j][0]);
396  double exponent = std::pow( dist/lambda_val_ , alpha_);
397  double expterm = std::exp(-0.5*exponent);
398  ret(i,j) = 2*expterm;
399  if (i!=j) ret(j,i) = ret(i,j);
400  }
401  }
402  } else if (particle_a == 1 && particle_b == 1) { //d2f/dlambda2
403  for (unsigned i=0; i<N; i++)
404  {
405  for (unsigned j=i; j<N; j++)
406  {
407  double dist = std::abs(xlist[i][0]-xlist[j][0]);
408  double exponent = std::pow( dist/lambda_val_ , alpha_);
409  double expterm = std::exp(-0.5*exponent);
410  ret(i,j) = tau_val_*tau_val_*expterm
411  *exponent/(lambda_val_*lambda_val_)*alpha_/2
412  * (alpha_/2 * exponent - (alpha_+1));
413  if (i!=j) ret(j,i) = ret(i,j);
414  }
415  }
416  } else { // d2f/d(tau)d(lambda)
417  for (unsigned i=0; i<N; i++)
418  {
419  for (unsigned j=i; j<N; j++)
420  {
421  double dist = std::abs(xlist[i][0]-xlist[j][0]);
422  double exponent = std::pow( dist/lambda_val_ , alpha_);
423  double expterm = std::exp(-0.5*exponent);
424  ret(i,j) = tau_val_ * alpha_ * expterm / lambda_val_
425  * exponent;
426  if (i!=j) ret(j,i) = ret(i,j);
427  }
428  }
429  }
430  return ret;
431  }
432 
434  unsigned particle_a, unsigned particle_b,
435  const FloatsList& xlist, bool) const
436  {
437  Eigen::MatrixXd mat( get_second_derivative_matrix(
438  particle_a, particle_b, xlist));
439  FloatsList ret;
440  for (int i=0; i<mat.rows(); i++)
441  {
442  Floats line;
443  for (int j=0; j<mat.cols(); j++)
444  line.push_back(mat(i,j));
445  ret.push_back(line);
446  }
447  return ret;
448  }
449 
450  unsigned get_ndims_x1() const {return 1;}
451  unsigned get_ndims_x2() const {return 1;}
452  unsigned get_ndims_y() const {return 1;}
453 
454  unsigned get_number_of_particles() const { return 2; }
455 
456  bool get_particle_is_optimized(unsigned particle_no) const
457  {
458  switch(particle_no)
459  {
460  case 0: //tau
461  return Scale(tau_).get_nuisance_is_optimized();
462  case 1: //lambda
463  return Scale(lambda_).get_nuisance_is_optimized();
464  default:
465  IMP_THROW("Invalid particle number", ModelException);
466  }
467  }
468 
470  {
471  unsigned count=0;
472  if (Scale(tau_).get_nuisance_is_optimized()) count++;
473  if (Scale(lambda_).get_nuisance_is_optimized()) count++;
474  return count;
475  }
476 
478  {
480  ret.push_back(tau_);
481  ret.push_back(lambda_);
482  return ret;
483  }
484 
486  {
487  ContainersTemp ret;
488  return ret;
489  }
490 
491 
492  /*IMP_OBJECT_INLINE(Covariance1DFunction,
493  out << "covariance function with alpha = "
494  << alpha_ << std::endl, {});*/
496 
497 
498  private:
499 
500  inline double get_value(double x1, double x2) const
501  {
502  double dist = std::abs(x1-x2);
503  double ret = dist /lambda_val_ ;
504  if (alpha_square_)
505  {
506  ret *= ret;
507  } else {
508  ret = std::pow(ret, alpha_);
509  }
510  ret = IMP::square(tau_val_) *std::exp(-0.5*ret);
511  if (do_jitter && dist<IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM)
512  {
513  ret +=IMP::square(tau_val_) * J_;
514  }
516  "function value is nan. tau = "
517  << tau_val_ << " lambda = " << lambda_val_
518  << " q1 = " << x1 << " q2 = " << x2);
519  return ret;
520  }
521 
522  private:
523  double alpha_;
524  base::Pointer<kernel::Particle> tau_,lambda_;
525  double tau_val_,lambda_val_,J_,cutoff_,alpha_square_;
526  bool do_jitter;
527 
528 };
529 
530 IMPISD_END_NAMESPACE
531 
532 #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.
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
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.
Eigen::MatrixXd operator()(const IMP::FloatsList &xlist) const
evaluate the function at a list of points
#define IMP_LOG_TERSE(expr)
ContainersTemp get_input_containers(const ModelObjectsTemp &mos)
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
#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
#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.
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
Class to handle individual model particles.
Common base class for heavy weight IMP objects.
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
#define IMP_LOG_VERBOSE(expr)
virtual Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const =0
return derivative matrix