IMP  2.0.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/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 ParticlesTemp get_input_particles() const = 0;
109  virtual ContainersTemp get_input_containers() const = 0;
110 
111  IMP_REF_COUNTED_DESTRUCTOR(BivariateFunction);
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:
129  Covariance1DFunction(Particle* tau, Particle* ilambda,
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);
136  IMP_IF_CHECK(USAGE_AND_INTERNAL) { Scale::decorate_particle(tau); }
137  IMP_IF_CHECK(USAGE_AND_INTERNAL) { Scale::decorate_particle(ilambda);}
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  }
144 
145  bool has_changed() const {
146  double tmpt = Scale(tau_).get_nuisance();
147  double tmpl = Scale(lambda_).get_nuisance();
148  if ((std::abs(tmpt - tau_val_) >
149  IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM)
150  || (std::abs(tmpl - lambda_val_) >
151  IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM))
152  {
153  IMP_LOG_TERSE( "Covariance1DFunction: has_changed():");
154  IMP_LOG_TERSE( "true" << std::endl);
155  return true;
156  } else {
157  return false;
158  }
159  }
160 
161  void update() {
162  lambda_val_= Scale(lambda_).get_nuisance();
163  tau_val_= Scale(tau_).get_nuisance();
164  IMP_LOG_TERSE( "Covariance1DFunction: update() tau:= "
165  << tau_val_ << " lambda:=" << lambda_val_
166  << std::endl);
167  IMP_INTERNAL_CHECK(!base::isnan(tau_val_),
168  "tau is nan.");
169  IMP_INTERNAL_CHECK(!base::isnan(lambda_val_),
170  "lambda is nan.");
171  IMP_INTERNAL_CHECK(tau_val_ >= 0, "tau is negative.");
172  IMP_INTERNAL_CHECK(lambda_val_ >= 0, "lambda is negative.");
173  IMP_INTERNAL_CHECK(lambda_val_ != 0, "lambda is zero.");
174  }
175 
176  Floats operator()(const Floats& x1,
177  const Floats& x2) const
178  {
179  IMP_USAGE_CHECK(x1.size() == 1, "expecting a 1-D vector");
180  IMP_USAGE_CHECK(x2.size() == 1, "expecting a 1-D vector");
181  Floats ret(1, get_value(x1[0],x2[0]));
182  return ret;
183  }
184 
185  Eigen::MatrixXd operator()(const IMP::FloatsList& xlist) const
186  {
187  const unsigned M=xlist.size();
188  Eigen::MatrixXd Mret(M,M);
189  for (unsigned i=0; i<M; i++)
190  {
191  for (unsigned j=i; j<M; j++)
192  {
193  IMP_USAGE_CHECK(xlist[i].size() == 1,
194  "expecting a 1-D vector");
195  IMP_USAGE_CHECK(xlist[j].size() == 1,
196  "expecting a 1-D vector");
197  double x1 = xlist[i][0];
198  double x2 = xlist[j][0];
199  double ret = get_value(x1,x2);
200  Mret(i,j) = ret;
201  if (i != j) Mret(j,i) = ret;
202  }
203  }
204  return Mret;
205  }
206 
207  FloatsList operator()(const IMP::FloatsList& xlist, bool) const
208  {
209  Eigen::MatrixXd mat((*this)(xlist));
210  FloatsList ret;
211  for (unsigned i=0; i<xlist.size(); i++)
212  for (unsigned j=0; j<xlist.size(); j++)
213  ret.push_back(Floats(1,mat(i,j)));
214  return ret;
215  }
216 
217  void add_to_derivatives(const Floats& x1,
218  const Floats& x2,
219  DerivativeAccumulator &accum) const
220  {
221  //d[w(x1,x2)]/dtau = 2/tau*(w(x1,x2))
222  double val = get_value(x1[0],x2[0]);
223  double tauderiv = 2./tau_val_ * val;
224  IMP_INTERNAL_CHECK(!base::isnan(tauderiv),
225  "tau derivative is nan.");
226  Scale(tau_).add_to_nuisance_derivative(tauderiv, accum);
227  //d[w(x,x')]/dlambda
228  // = w(x,x') ( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
229  double lambdaderiv =
230  val * (alpha_ *
231  std::pow((std::abs(x1[0]-x2[0])/lambda_val_),alpha_)
232  /(2.*lambda_val_));
233  IMP_INTERNAL_CHECK(!base::isnan(lambdaderiv),
234  "lambda derivative is nan.");
235  Scale(lambda_).add_to_nuisance_derivative(lambdaderiv, accum);
236  }
237 
238  void add_to_particle_derivative(unsigned particle_no,
239  double value, DerivativeAccumulator &accum) const
240  {
241  switch (particle_no)
242  {
243  case 0: //tau
244  IMP_INTERNAL_CHECK(!base::isnan(value),
245  "tau derivative is nan.");
246  Scale(tau_).add_to_nuisance_derivative(value, accum);
247  break;
248  case 1: //lambda
249  IMP_INTERNAL_CHECK(!base::isnan(value),
250  "lambda derivative is nan.");
251  Scale(lambda_).add_to_nuisance_derivative(value, accum);
252  break;
253  default:
254  IMP_THROW("Invalid particle number", ModelException);
255  }
256  }
257 
258  Eigen::MatrixXd get_derivative_matrix(
259  unsigned particle_no,
260  const FloatsList& xlist) const
261  {
262  // Strategy: fill in the main diagonal, then fill with zeros
263  // if the value of the function falls below cutoff.
264  // assumes data points are ordered!
265  unsigned N=xlist.size();
266  Eigen::MatrixXd ret(N,N);
267  double diag;
268  switch (particle_no)
269  {
270  case 0: //tau
271  //d[w(x1,x1)]/dtau
272  // = 2/tau*w(x1,x1)
273  diag = get_value(xlist[0][0],xlist[0][0]);
274  diag *= 2./tau_val_;
275  break;
276  case 1: //lambda
277  //d[w(x,x)]/dlambda
278  //= w(x,x)
279  // *( alpha /(2 lambda^{alpha+1}))
280  diag = 0;
281  break;
282  default:
283  IMP_THROW("Invalid particle number",
284  ModelException);
285  }
286  IMP_INTERNAL_CHECK(!base::isnan(diag),
287  "derivative matrix is nan on the diagonal.");
288  for (unsigned i=0; i<N; i++) ret(i,i) = diag;
289  //
290  bool initial_loop=true;
291  double abs_cutoff = cutoff_*diag;
292  double dmax=-1;
293  for (unsigned i=0; i<N; i++)
294  {
295  for (unsigned j=i+1; j<N; j++)
296  {
297  double x1(xlist[i][0]), x2(xlist[j][0]);
298  double val;
299  double dist(std::abs(x1-x2));
300  //compute all entries as long as the cutoff distance was
301  //not recorded (initial_loop) or as long as the distance is
302  //smaller than the cutoff distance
303  if (initial_loop || dist <= dmax)
304  {
305  switch (particle_no)
306  {
307  case 0: //tau
308  //d[w(x1,x2)]/dtau
309  // = 2/tau*w(x1,x2)
310  val = get_value(xlist[i][0],xlist[j][0]);
311  val *= 2./tau_val_;
312  break;
313  case 1: //lambda
314  //d[w(x,x')]/dlambda
315  //= w(x,x')
316  // *( alpha |x'-x|^alpha/(2 lambda^{alpha+1}))
317  if (dist<IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
318  val = 0;
319  } else {
320  val = get_value(xlist[i][0],xlist[j][0]);
321  val *= alpha_ *
322  std::pow(
323  (dist/lambda_val_), alpha_)
324  /(2.*lambda_val_);
325  }
326  break;
327  default:
328  IMP_THROW("Invalid particle number",
329  ModelException);
330  }
331  // the value has been computed and is in val
332  //now check if it is smaller than the cutoff.
333  //If true change the flag and update the distance
334  if (std::abs(val) <= abs_cutoff)
335  {
336  if (initial_loop)
337  {
338  initial_loop = false;
339  dmax = dist;
340  } else if (dist < dmax)
341  {
342  dmax = dist;
343  }
344  }
345  } else { // e.g. initial_loop == false && dist > dmax
346  val = 0;
347  }
348  IMP_INTERNAL_CHECK(!base::isnan(val),
349  "derivative matrix is nan at position("
350  << i << "," << j << ").");
351  ret(i,j) = val;
352  ret(j,i) = val;
353  }
354  }
355  return ret;
356  }
357 
358  FloatsList get_derivative_matrix(
359  unsigned particle_no,
360  const FloatsList& xlist, bool) const
361  {
362  Eigen::MatrixXd mat(get_derivative_matrix(particle_no, xlist));
363  FloatsList ret;
364  for (int i=0; i<mat.rows(); i++)
365  {
366  Floats line;
367  for (int j=0; j<mat.cols(); j++)
368  line.push_back(mat(i,j));
369  ret.push_back(line);
370  }
371  return ret;
372  }
373 
374  Eigen::MatrixXd get_second_derivative_matrix(
375  unsigned particle_a, unsigned particle_b,
376  const FloatsList& xlist) const
377  {
378  unsigned N(xlist.size());
379  Eigen::MatrixXd ret(N,N);
380  if (particle_a > 1)
381  IMP_THROW("Invalid particle 1 number", ModelException);
382  if (particle_b > 1)
383  IMP_THROW("Invalid particle 2 number", ModelException);
384  //d2f/dtau2
385  if (particle_a == 0 && particle_b == 0)
386  {
387  for (unsigned i=0; i<N; i++)
388  {
389  for (unsigned j=i; j<N; j++)
390  {
391  double dist = std::abs(xlist[i][0]-xlist[j][0]);
392  double exponent = std::pow( dist/lambda_val_ , alpha_);
393  double expterm = std::exp(-0.5*exponent);
394  ret(i,j) = 2*expterm;
395  if (i!=j) ret(j,i) = ret(i,j);
396  }
397  }
398  } else if (particle_a == 1 && particle_b == 1) { //d2f/dlambda2
399  for (unsigned i=0; i<N; i++)
400  {
401  for (unsigned j=i; j<N; j++)
402  {
403  double dist = std::abs(xlist[i][0]-xlist[j][0]);
404  double exponent = std::pow( dist/lambda_val_ , alpha_);
405  double expterm = std::exp(-0.5*exponent);
406  ret(i,j) = tau_val_*tau_val_*expterm
407  *exponent/(lambda_val_*lambda_val_)*alpha_/2
408  * (alpha_/2 * exponent - (alpha_+1));
409  if (i!=j) ret(j,i) = ret(i,j);
410  }
411  }
412  } else { // d2f/d(tau)d(lambda)
413  for (unsigned i=0; i<N; i++)
414  {
415  for (unsigned j=i; j<N; j++)
416  {
417  double dist = std::abs(xlist[i][0]-xlist[j][0]);
418  double exponent = std::pow( dist/lambda_val_ , alpha_);
419  double expterm = std::exp(-0.5*exponent);
420  ret(i,j) = tau_val_ * alpha_ * expterm / lambda_val_
421  * exponent;
422  if (i!=j) ret(j,i) = ret(i,j);
423  }
424  }
425  }
426  return ret;
427  }
428 
429  FloatsList get_second_derivative_matrix(
430  unsigned particle_a, unsigned particle_b,
431  const FloatsList& xlist, bool) const
432  {
433  Eigen::MatrixXd mat( get_second_derivative_matrix(
434  particle_a, particle_b, xlist));
435  FloatsList ret;
436  for (int i=0; i<mat.rows(); i++)
437  {
438  Floats line;
439  for (int j=0; j<mat.cols(); j++)
440  line.push_back(mat(i,j));
441  ret.push_back(line);
442  }
443  return ret;
444  }
445 
446  unsigned get_ndims_x1() const {return 1;}
447  unsigned get_ndims_x2() const {return 1;}
448  unsigned get_ndims_y() const {return 1;}
449 
450  unsigned get_number_of_particles() const { return 2; }
451 
452  bool get_particle_is_optimized(unsigned particle_no) const
453  {
454  switch(particle_no)
455  {
456  case 0: //tau
457  return Scale(tau_).get_nuisance_is_optimized();
458  case 1: //lambda
459  return Scale(lambda_).get_nuisance_is_optimized();
460  default:
461  IMP_THROW("Invalid particle number", ModelException);
462  }
463  }
464 
465  unsigned get_number_of_optimized_particles() const
466  {
467  unsigned count=0;
468  if (Scale(tau_).get_nuisance_is_optimized()) count++;
469  if (Scale(lambda_).get_nuisance_is_optimized()) count++;
470  return count;
471  }
472 
473  ParticlesTemp get_input_particles() const
474  {
475  ParticlesTemp ret;
476  ret.push_back(tau_);
477  ret.push_back(lambda_);
478  return ret;
479  }
480 
481  ContainersTemp get_input_containers() const
482  {
483  ContainersTemp ret;
484  return ret;
485  }
486 
487 
488  IMP_OBJECT_INLINE(Covariance1DFunction,
489  out << "covariance function with alpha = "
490  << alpha_ << std::endl, {});
491 
492 
493  private:
494 
495  inline double get_value(double x1, double x2) const
496  {
497  double dist = std::abs(x1-x2);
498  double ret = dist /lambda_val_ ;
499  if (alpha_square_)
500  {
501  ret *= ret;
502  } else {
503  ret = std::pow(ret, alpha_);
504  }
505  ret = IMP::square(tau_val_) *std::exp(-0.5*ret);
506  if (do_jitter && dist<IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM)
507  {
508  ret += J_;
509  }
510  IMP_INTERNAL_CHECK(!base::isnan(ret),
511  "function value is nan. tau = "
512  << tau_val_ << " lambda = " << lambda_val_
513  << " q1 = " << x1 << " q2 = " << x2);
514  return ret;
515  }
516 
517  private:
518  double alpha_;
519  Pointer<Particle> tau_,lambda_;
520  double tau_val_,lambda_val_,J_,cutoff_,alpha_square_;
521  bool do_jitter;
522 
523 };
524 
525 IMPISD_END_NAMESPACE
526 
527 #endif /* IMPISD_BIVARIATE_FUNCTIONS_H */