IMP logo
IMP Reference Guide  2.8.0
The Integrative Modeling Platform
univariate_functions.h
Go to the documentation of this file.
1 /**
2  * \file IMP/isd/univariate_functions.h
3  * \brief Classes for general functions
4  *
5  * Copyright 2007-2017 IMP Inventors. All rights reserved.
6  */
7 
8 #ifndef IMPISD_UNIVARIATE_FUNCTIONS_H
9 #define IMPISD_UNIVARIATE_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 <IMP/algebra/eigen3/Eigen/Dense>
18 
19 #define IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM 1e-7
20 
21 IMPISD_BEGIN_NAMESPACE
22 
23 //! Base class for functions of one variable
24 class IMPISDEXPORT UnivariateFunction : public Object {
25  public:
26  UnivariateFunction(std::string str) : Object(str) {}
27 
28  //! evaluate the function at a certain point
29  virtual Floats operator()(const Floats& x) const = 0;
30 
31  //! evaluate the function at a list of points
32  virtual IMP_Eigen::VectorXd operator()(const IMP::FloatsList& xlist)
33  const = 0;
34 
35  //! used for testing only
36  virtual FloatsList operator()(const IMP::FloatsList& xlist,
37  bool stupid) const = 0;
38 
39  //! return true if internal parameters have changed.
40  virtual bool has_changed() const = 0;
41 
42  //! update internal parameters
43  virtual void update() = 0;
44 
45  //! update derivatives of particles
46  /* add to each particle the derivative of the function
47  * times the weight of the DA.
48  */
49  virtual void add_to_derivatives(const Floats& x,
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 vector
62  /* This function returns a column from the matrix m:
63  * m_ij = d(func(xlist[i]))/dparticle_j
64  * The matrix has N rows and M columns
65  * where N = xlist.size() and M is the number of particles
66  * associated to this function.
67  */
68  virtual IMP_Eigen::VectorXd get_derivative_vector(
69  unsigned particle_no, const FloatsList& xlist) const = 0;
70 
71  //! for testing purposes
72  virtual FloatsList get_derivative_matrix(const FloatsList& xlist,
73  bool stupid) const = 0;
74 
75  //! return second derivative vector
76  virtual IMP_Eigen::VectorXd get_second_derivative_vector(
77  unsigned particle_a, unsigned particle_b,
78  const FloatsList& xlist) const = 0;
79 
80  //! for testing purposes
81  virtual FloatsList get_second_derivative_vector(unsigned particle_a,
82  unsigned particle_b,
83  const FloatsList& xlist,
84  bool stupid) const = 0;
85 
86  //! returns the number of input dimensions
87  virtual unsigned get_ndims_x() const = 0;
88 
89  //! returns the number of output dimensions
90  virtual unsigned get_ndims_y() const = 0;
91 
92  //! returns the number of particles that this function uses
93  virtual unsigned get_number_of_particles() const = 0;
94 
95  //! returns true if the particle whose index is provided is optimized
96  virtual bool get_particle_is_optimized(unsigned particle_no) const = 0;
97 
98  //! returns the number of particles that are optimized
99  virtual unsigned get_number_of_optimized_particles() const = 0;
100 
101  //! particle manipulation
102  virtual ModelObjectsTemp get_inputs() const = 0;
103 
105 };
106 //
107 //! Linear one-dimensional function
108 /* f(x) = a*x + b, where a,b are ISD nuisances, and f(x) and x are doubles.
109  */
110 class IMPISDEXPORT Linear1DFunction : public UnivariateFunction {
111  public:
113  : UnivariateFunction("Linear1DFunction %1%"), a_(a), b_(b) {
114  IMP_LOG_TERSE("Linear1DFunction: constructor" << std::endl);
115  a_val_ = Nuisance(a).get_nuisance();
116  b_val_ = Nuisance(b).get_nuisance();
117  update();
118  }
119 
120  bool has_changed() const {
121  double tmpa = Nuisance(a_).get_nuisance();
122  double tmpb = Nuisance(b_).get_nuisance();
123  if ((std::abs(tmpa - a_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
124  (std::abs(tmpb - b_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
125  IMP_LOG_TERSE("Linear1DFunction: has_changed():");
126  IMP_LOG_TERSE("true" << std::endl);
127  return true;
128  } else {
129  return false;
130  }
131  }
132 
133  void update() {
134  a_val_ = Nuisance(a_).get_nuisance();
135  b_val_ = Nuisance(b_).get_nuisance();
136  IMP_LOG_TERSE("Linear1DFunction: update() a:= " << a_val_ << " b:="
137  << b_val_ << std::endl);
138  }
139 
140  Floats operator()(const Floats& x) const {
141  IMP_USAGE_CHECK(x.size() == 1, "expecting a 1-D vector");
142  Floats ret(1, a_val_ * x[0] + b_val_);
143  return ret;
144  }
145 
146  IMP_Eigen::VectorXd operator()(const FloatsList& xlist) const {
147  unsigned M = xlist.size();
148  IMP_Eigen::VectorXd retlist(M);
149  for (unsigned i = 0; i < M; i++) {
150  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
151  retlist(i) = a_val_ * xlist[i][0] + b_val_;
152  }
153  return retlist;
154  }
155 
156  FloatsList operator()(const FloatsList& xlist, bool) const {
157  IMP_Eigen::VectorXd vec((*this)(xlist));
158  FloatsList ret;
159  for (unsigned i = 0; i < xlist.size(); i++)
160  ret.push_back(Floats(1, vec(i)));
161  return ret;
162  }
163 
164  void add_to_derivatives(const Floats& x, DerivativeAccumulator& accum) const {
165  // d[f(x)]/da = x
166  Nuisance(a_).add_to_nuisance_derivative(x[0], accum);
167  // d[f(x)]/db = 1
168  Nuisance(b_).add_to_nuisance_derivative(1, accum);
169  }
170 
171  void add_to_particle_derivative(unsigned particle_no, double value,
172  DerivativeAccumulator& accum) const {
173  switch (particle_no) {
174  case 0:
175  Nuisance(a_).add_to_nuisance_derivative(value, accum);
176  break;
177  case 1:
178  Nuisance(b_).add_to_nuisance_derivative(value, accum);
179  break;
180  default:
181  IMP_THROW("Invalid particle number", ModelException);
182  }
183  }
184 
185  IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no,
186  const FloatsList& xlist) const {
187  unsigned N = xlist.size();
188  IMP_Eigen::VectorXd ret(N);
189  switch (particle_no) {
190  case 0: // a
191  for (unsigned i = 0; i < N; i++) ret(i) = xlist[i][0];
192  break;
193  case 1: // b
194  ret.setOnes();
195  break;
196  default:
197  IMP_THROW("Invalid particle number", ModelException);
198  }
199  return ret;
200  }
201 
202  FloatsList get_derivative_matrix(const FloatsList& xlist, bool) const {
203  IMP_Eigen::MatrixXd mat(xlist.size(), 2);
204  mat.col(0) = get_derivative_vector(0, xlist);
205  mat.col(1) = get_derivative_vector(1, xlist);
206  FloatsList ret;
207  for (int i = 0; i < mat.rows(); i++) {
208  Floats line;
209  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
210  ret.push_back(line);
211  }
212  return ret;
213  }
214 
215  IMP_Eigen::VectorXd get_second_derivative_vector(
216  unsigned, unsigned, const FloatsList& xlist) const {
217  // The Hessian is zero for all particles.
218  unsigned N = xlist.size();
219  IMP_Eigen::VectorXd H(IMP_Eigen::VectorXd::Zero(N));
220  return H;
221  }
222 
224  unsigned particle_b,
225  const FloatsList& xlist, bool) const {
226  IMP_Eigen::VectorXd mat(
227  get_second_derivative_vector(particle_a, particle_b, xlist));
228  FloatsList ret;
229  for (int i = 0; i < mat.rows(); i++) {
230  Floats line;
231  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
232  ret.push_back(line);
233  }
234  return ret;
235  }
236 
237  unsigned get_ndims_x() const { return 1; }
238  unsigned get_ndims_y() const { return 1; }
239 
240  unsigned get_number_of_particles() const { return 2; }
241 
242  bool get_particle_is_optimized(unsigned particle_no) const {
243  switch (particle_no) {
244  case 0: // a
245  return Nuisance(a_).get_nuisance_is_optimized();
246  case 1: // b
247  return Nuisance(b_).get_nuisance_is_optimized();
248  default:
249  IMP_THROW("Invalid particle number", ModelException);
250  }
251  }
252 
254  unsigned count = 0;
255  if (Nuisance(a_).get_nuisance_is_optimized()) count++;
256  if (Nuisance(b_).get_nuisance_is_optimized()) count++;
257  return count;
258  }
259 
261  ModelObjectsTemp ret;
262  ret.push_back(a_);
263  ret.push_back(b_);
264  return ret;
265  }
266 
267  /*IMP_OBJECT_INLINE(Linear1DFunction, out << "y = " << a_val_
268  << " * x + " << b_val_ << std::endl, {});*/
270 
271  private:
272  Pointer<Particle> a_, b_;
273  double a_val_, b_val_;
274 };
275 
276 //! 1D mean function for SAS data
277 /* Generalized Guinier-Porod model (Hammouda, J. Appl. Cryst., 2010, eq 3 & 4
278  * to which a constant offset is added)
279  * I(q) = A + G/q^s exp(-(q.Rg)^2/(3-s)) for q <= q1
280  * I(q) = A + D/q^d for q > q1
281  * q1 = 1/Rg * ((d-s)(3-s)/2)^(1/2)
282  * D = G exp(-(q1.Rg)^2/(3-s)) q1^(d-s)
283  * only valid for q>0 Rg>0 0<s<d s<3 G>0
284  * s=0 in the globular case.
285  * G, Rg, d, s and A are particles (ISD Scales).
286  */
288  public:
290  Particle* d, Particle* s,
291  Particle* A)
292  : UnivariateFunction("GeneralizedGuinierPorodFunction %1%"),
293  G_(G),
294  Rg_(Rg),
295  d_(d),
296  s_(s),
297  A_(A) {
298  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: constructor" << std::endl);
299  update();
300  }
301 
302  bool has_changed() const {
303  double tmpG = Scale(G_).get_scale();
304  double tmpRg = Scale(Rg_).get_scale();
305  double tmpd = Scale(d_).get_scale();
306  double tmps = Scale(s_).get_scale();
307  double tmpA = Nuisance(A_).get_nuisance();
308  if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
309  (std::abs(tmpRg - Rg_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
310  (std::abs(tmpd - d_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
311  (std::abs(tmps - s_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
312  (std::abs(tmpA - A_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
313  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: has_changed():");
314  IMP_LOG_TERSE("true" << std::endl);
315  return true;
316  } else {
317  return false;
318  }
319  }
320 
321  void update() {
322  G_val_ = Scale(G_).get_scale();
323  Rg_val_ = Scale(Rg_).get_scale();
324  d_val_ = Scale(d_).get_scale();
325  s_val_ = Scale(s_).get_scale();
326  if (d_val_ == s_val_) {
327  IMP_LOG_TERSE("Warning: d==s !" << std::endl);
328  if (s_val_ > 0.001) {
329  s_val_ -= 0.001;
330  } else {
331  d_val_ += 0.001;
332  }
333  // IMP_THROW("d == s ! ", ModelException);
334  }
335  A_val_ = Nuisance(A_).get_nuisance();
336  q1_param_ = std::sqrt((d_val_ - s_val_) * (3 - s_val_) / 2.);
337  D_param_ = G_val_ * std::exp(-IMP::square(q1_param_) / (3 - s_val_));
338  q1_param_ = q1_param_ / Rg_val_;
339  D_param_ *= std::pow(q1_param_, d_val_ - s_val_);
340  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: update() G:= "
341  << G_val_ << " Rg:=" << Rg_val_ << " d:=" << d_val_
342  << " s:=" << s_val_ << " A:=" << A_val_ << " Q1.Rg ="
343  << q1_param_ * Rg_val_ << " D =" << D_param_ << std::endl);
344  /*std::cout << "cpp"
345  << " 3:Q1 " << q1_param_
346  << " 5:D " << D_param_
347  << " 7:G " << G_val_
348  << " 9:Rg " << Rg_val_
349  << " 11:d " << d_val_
350  << " 13:s " << s_val_
351  << " 15:val " << get_value(0.1)
352  <<std::endl;*/
353  }
354 
355  Floats operator()(const Floats& x) const {
356  IMP_USAGE_CHECK(x.size() == 1, "expecting a 1-D vector");
357  Floats ret(1, get_value(x[0]));
358  return ret;
359  }
360 
361  IMP_Eigen::VectorXd operator()(const FloatsList& xlist) const {
362  unsigned M = xlist.size();
363  IMP_Eigen::VectorXd retlist(M);
364  for (unsigned i = 0; i < M; i++) {
365  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
366  retlist(i) = get_value(xlist[i][0]);
367  }
368  return retlist;
369  }
370 
371  FloatsList operator()(const FloatsList& xlist, bool) const {
372  IMP_Eigen::VectorXd vec((*this)(xlist));
373  FloatsList ret;
374  for (unsigned i = 0; i < xlist.size(); i++)
375  ret.push_back(Floats(1, vec(i)));
376  return ret;
377  }
378 
379  void add_to_derivatives(const Floats& x, DerivativeAccumulator& accum) const {
380  double qval = x[0];
381  double value = get_value(qval) - A_val_;
382  double deriv;
383  // d[f(x)+A]/dG = f(x)/G
384  deriv = value / G_val_;
385  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for G is nan.");
386  Scale(G_).add_to_nuisance_derivative(deriv, accum);
387  if (qval <= q1_param_) {
388  // d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
389  deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
390  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for Rg is nan.");
391  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
392  // d[f(x)]/dd = 0
393  //
394  // d[f(x)]/ds = - f(x) * ( (q Rg / (3-s))^2 + log(q) )
395  deriv = -value *
396  (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
397  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for s is nan.");
398  Scale(s_).add_to_nuisance_derivative(deriv, accum);
399  } else {
400  // d[f(x)]/dRg = f(x) * (s-d)/Rg
401  deriv = value * (s_val_ - d_val_) / Rg_val_;
402  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for Rg is nan.");
403  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
404  // d[f(x)]/dd = f(x) * log(q1/q)
405  deriv = value * std::log(q1_param_ / qval);
406  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for d is nan.");
407  Scale(d_).add_to_nuisance_derivative(deriv, accum);
408  // d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
409  deriv = -value *
410  ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
411  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for d is nan.");
412  Scale(s_).add_to_nuisance_derivative(deriv, accum);
413  }
414  // d[f(x)+A]/dA = 1
415  deriv = 1;
416  Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
417  }
418 
419  void add_to_particle_derivative(unsigned particle_no, double value,
420  DerivativeAccumulator& accum) const {
421  switch (particle_no) {
422  case 0:
423  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for G is nan.");
424  Scale(G_).add_to_scale_derivative(value, accum);
425  break;
426  case 1:
427  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for Rg is nan.");
428  Scale(Rg_).add_to_scale_derivative(value, accum);
429  break;
430  case 2:
431  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for d is nan.");
432  Scale(d_).add_to_scale_derivative(value, accum);
433  break;
434  case 3:
435  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for s is nan.");
436  Scale(s_).add_to_scale_derivative(value, accum);
437  break;
438  case 4:
439  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for A is nan.");
440  Nuisance(A_).add_to_nuisance_derivative(value, accum);
441  break;
442  default:
443  IMP_THROW("Invalid particle number", ModelException);
444  }
445  }
446 
447  IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no,
448  const FloatsList& xlist) const {
449  unsigned N = xlist.size();
450  IMP_Eigen::VectorXd ret(N);
451  switch (particle_no) {
452  case 0: // G
453  // d[f(x)]/dG = f(x)/G
454  ret = ((*this)(xlist) - IMP_Eigen::VectorXd::Constant(N, A_val_)) /
455  G_val_;
456  break;
457  case 1: // Rg
458  for (unsigned i = 0; i < N; i++) {
459  double qval = xlist[i][0];
460  if (qval <= q1_param_) {
461  // d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
462  ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
463  Rg_val_ / (3 - s_val_);
464  } else {
465  // d[f(x)]/dRg = f(x) * (s-d)/Rg
466  ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
467  }
468  }
469  break;
470  case 2: // d
471  for (unsigned i = 0; i < N; i++) {
472  double qval = xlist[i][0];
473  if (qval <= q1_param_) {
474  // d[f(x)]/dd = 0
475  ret(i) = 0;
476  } else {
477  // d[f(x)]/dd = f(x) * log(q1/q)
478  ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
479  }
480  }
481  break;
482  case 3: // s
483  for (unsigned i = 0; i < N; i++) {
484  double qval = xlist[i][0];
485  if (qval <= q1_param_) {
486  // d[f(x)]/ds = - f(x)
487  // * (q^2 Rg^2 + (3-s)^2 log(q)) / (3-s)^2
488  ret(i) =
489  -(get_value(qval) - A_val_) *
490  (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
491  } else {
492  // d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
493  ret(i) =
494  -(get_value(qval) - A_val_) *
495  ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
496  }
497  }
498  break;
499  case 4: // A
500  ret = IMP_Eigen::VectorXd::Constant(N, 1);
501  break;
502  default:
503  IMP_THROW("Invalid particle number", ModelException);
504  }
505  return ret;
506  }
507 
508  FloatsList get_derivative_matrix(const FloatsList& xlist, bool) const {
509  IMP_Eigen::MatrixXd mat(xlist.size(), 5);
510  mat.col(0) = get_derivative_vector(0, xlist);
511  mat.col(1) = get_derivative_vector(1, xlist);
512  mat.col(2) = get_derivative_vector(2, xlist);
513  mat.col(3) = get_derivative_vector(3, xlist);
514  mat.col(4) = get_derivative_vector(4, xlist);
515  FloatsList ret;
516  for (int i = 0; i < mat.rows(); i++) {
517  Floats line;
518  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
519  ret.push_back(line);
520  }
521  return ret;
522  }
523 
524  IMP_Eigen::VectorXd get_second_derivative_vector(
525  unsigned particle_a, unsigned particle_b, const FloatsList& xlist) const {
526  if (particle_a >= 5) IMP_THROW("Invalid particle 1 number", ModelException);
527  if (particle_b >= 5) IMP_THROW("Invalid particle 2 number", ModelException);
528  unsigned N = xlist.size();
529  // hessian involving A is always 0
530  if (particle_a == 4 || particle_b == 4) return IMP_Eigen::VectorXd::Zero(N);
531  unsigned pa = std::min(particle_a, particle_b);
532  unsigned pb = std::max(particle_a, particle_b);
533  IMP_Eigen::VectorXd ret(N);
534  switch (pa) {
535  case 0: // G
536  switch (pb) {
537  case 0: // G
538  // d2f/dG2 = 0
539  ret.noalias() = IMP_Eigen::VectorXd::Zero(N);
540  break;
541 
542  case 1: // Rg
543  // d2f/dGdRg = df/dRg * 1/G
544  ret.noalias() = get_derivative_vector(1, xlist) / G_val_;
545  break;
546 
547  case 2: // d
548  // d2f/dGdd = df/dd * 1/G
549  ret.noalias() = get_derivative_vector(2, xlist) / G_val_;
550  break;
551 
552  case 3: // s
553  // d2f/dGds = df/ds * 1/G
554  ret.noalias() = get_derivative_vector(3, xlist) / G_val_;
555  break;
556 
557  default:
558  IMP_THROW("Invalid particle 2 number", ModelException);
559  break;
560  }
561  break;
562 
563  case 1: // Rg
564  switch (pb) {
565  case 1: // Rg
566  for (unsigned i = 0; i < N; i++) {
567  double qval = xlist[i][0];
568  if (qval <= q1_param_) {
569  // d2[f(x)]/dRg2 =
570  // f(x) * (2 q^2 / (3-s))
571  // * (2 q^2 Rg^2 / (3-s) - 1)
572  ret(i) = (get_value(qval) - A_val_) * 2 * IMP::square(qval) /
573  (3 - s_val_) *
574  (2 * IMP::square(qval * Rg_val_) / (3 - s_val_) - 1);
575  } else {
576  // d2[f(x)]/dRg2=f(x) * (d-s)/Rg^2 * (d-s+1)
577  ret(i) = (get_value(qval) - A_val_) * (d_val_ - s_val_) /
578  IMP::square(Rg_val_) * (d_val_ - s_val_ + 1);
579  }
580  }
581  break;
582 
583  case 2: // d
584  for (unsigned i = 0; i < N; i++) {
585  double qval = xlist[i][0];
586  if (qval <= q1_param_) {
587  // d2[f(x)]/dddRg = 0
588  ret(i) = 0;
589  } else {
590  // d2[f(x)]/dddRg = -f(x)/Rg
591  // - (d-s)/Rg *df(x)/dd
592  double val = (get_value(qval) - A_val_);
593  ret(i) = -val / Rg_val_ - (val * std::log(q1_param_ / qval) *
594  (d_val_ - s_val_) / Rg_val_);
595  }
596  }
597  break;
598 
599  case 3: // s
600  for (unsigned i = 0; i < N; i++) {
601  double qval = xlist[i][0];
602  double val = (get_value(qval) - A_val_);
603  if (qval <= q1_param_) {
604  // d2[f(x)]/dsdRg = -2q^2Rg/(3-s)
605  // * (df(x)/ds + f(x)/(3-s))
606  double deriv =
607  -val * (IMP::square((qval * Rg_val_) / (3 - s_val_)) +
608  std::log(qval));
609  ret(i) = -2 * IMP::square(qval) * Rg_val_ / (3 - s_val_) *
610  (deriv + val / (3 - s_val_));
611  } else {
612  // d2[f(x)]/dsdRg =
613  // 1/Rg * (f(x)- (d-s)*df(x)/ds)
614  double deriv = -val * ((d_val_ - s_val_) / (2 * (3 - s_val_)) +
615  std::log(q1_param_));
616  ret(i) = (val - (d_val_ - s_val_) * deriv) / Rg_val_;
617  }
618  }
619  break;
620 
621  default:
622  IMP_THROW("Invalid particle 2 number", ModelException);
623  break;
624  }
625  break;
626 
627  case 2: // d
628  switch (pb) {
629  case 2: // d
630  for (unsigned i = 0; i < N; i++) {
631  double qval = xlist[i][0];
632  if (qval <= q1_param_) {
633  // d2[f(x)]/dddd = 0
634  ret(i) = 0;
635  } else {
636  // d2[f(x)]/dddd =
637  // f(x)*(log(q1/q)^2 + 1/(2*(d-s)))
638  double val = (get_value(qval) - A_val_);
639  ret(i) = val * (IMP::square(std::log(q1_param_ / qval)) +
640  1 / (2 * (d_val_ - s_val_)));
641  }
642  }
643  break;
644 
645  case 3: // s
646  {
647  double lterm =
648  (d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_);
649  double rterm = 0.5 * (1 / (3 - s_val_) + 1 / (d_val_ - s_val_));
650  for (unsigned i = 0; i < N; i++) {
651  double qval = xlist[i][0];
652  if (qval <= q1_param_) {
653  // d2[f(x)]/ddds = 0
654  ret(i) = 0;
655  } else {
656  // d2[f(x)]/ddds = log(q1/q)*df(x)/ds
657  // - (1/(3-s) + 1/(d-s))*f(x)/2
658  //
659  double val = (get_value(qval) - A_val_);
660  ret(i) = -val * (std::log(q1_param_ / qval) * lterm + rterm);
661  }
662  }
663  } break;
664 
665  default:
666  IMP_THROW("Invalid particle 2 number", ModelException);
667  break;
668  }
669  break;
670 
671  case 3: // s
672  switch (pb) {
673  case 3: // s
674  {
675  double cterm =
676  IMP::square(0.5 * (d_val_ - s_val_) / (3 - s_val_) +
677  std::log(q1_param_)) +
678  0.5 * ((6 - s_val_ - d_val_) / IMP::square(3 - s_val_) +
679  1. / (d_val_ - s_val_));
680  for (unsigned i = 0; i < N; i++) {
681  double qval = xlist[i][0];
682  double val = (get_value(qval) - A_val_);
683  if (qval <= q1_param_) {
684  // d2[f(x)]/dsds = f(x) *
685  // ( [(qRg)^2/(3-s)^2 + log q ]^2
686  // - 2(qRg)^2/(3-s)^3 )
687  double factor = IMP::square((qval * Rg_val_) / (3 - s_val_));
688  ret(i) = val * (IMP::square(factor + std::log(qval)) -
689  2 * factor / (3 - s_val_));
690  } else {
691  // d2[f(x)]/dsds = f(x)
692  // * ( [ (d-s)/(2*(3-s)) + log(q1) ]^2
693  // + 1/2 * [ (6-s-d)/(3-s)^2
694  // + 1/(d-s) ] )
695  ret(i) = val * cterm;
696  }
697  }
698  } break;
699 
700  default:
701  IMP_THROW("Invalid particle 2 number", ModelException);
702  break;
703  }
704  break;
705 
706  default:
707  IMP_THROW("Invalid particle 1 number", ModelException);
708  break;
709  }
710  return ret;
711  }
712 
714  unsigned particle_b,
715  const FloatsList& xlist, bool) const {
716  IMP_Eigen::VectorXd mat(
717  get_second_derivative_vector(particle_a, particle_b, xlist));
718  FloatsList ret;
719  for (int i = 0; i < mat.rows(); i++) {
720  Floats line;
721  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
722  ret.push_back(line);
723  }
724  return ret;
725  }
726 
727  unsigned get_ndims_x() const { return 1; }
728  unsigned get_ndims_y() const { return 1; }
729 
730  unsigned get_number_of_particles() const { return 5; }
731 
732  bool get_particle_is_optimized(unsigned particle_no) const {
733  switch (particle_no) {
734  case 0: // G
735  return Scale(G_).get_scale_is_optimized();
736  case 1: // Rg
737  return Scale(Rg_).get_scale_is_optimized();
738  case 2: // d
739  return Scale(d_).get_scale_is_optimized();
740  case 3: // s
741  return Scale(s_).get_scale_is_optimized();
742  case 4: // A
743  return Nuisance(A_).get_nuisance_is_optimized();
744  default:
745  IMP_THROW("Invalid particle number", ModelException);
746  }
747  }
748 
750  unsigned count = 0;
751  if (Scale(G_).get_scale_is_optimized()) count++;
752  if (Scale(Rg_).get_scale_is_optimized()) count++;
753  if (Scale(d_).get_scale_is_optimized()) count++;
754  if (Scale(s_).get_scale_is_optimized()) count++;
755  if (Nuisance(A_).get_nuisance_is_optimized()) count++;
756  return count;
757  }
758 
760  ModelObjectsTemp ret;
761  ret.push_back(G_);
762  ret.push_back(Rg_);
763  ret.push_back(d_);
764  ret.push_back(s_);
765  ret.push_back(A_);
766  return ret;
767  }
768 
769  /*IMP_OBJECT_INLINE(GeneralizedGuinierPorodFunction, out
770  << " G = " << G_val_
771  << " Rg = " << Rg_val_
772  << " d = " << d_val_
773  << " s = " << s_val_
774  << " A = " << A_val_
775  << " Q1.Rg = " << q1_param_*Rg_val_
776  << std::endl, {});*/
778 
779  private:
780  inline double get_value(double q) const {
781  double value;
782  if (q <= q1_param_) {
783  value = A_val_ + G_val_ / std::pow(q, s_val_) *
784  std::exp(-IMP::square(q * Rg_val_) / (3 - s_val_));
785  } else {
786  value = A_val_ + D_param_ / std::pow(q, d_val_);
787  }
788  return value;
789  }
790 
791  Pointer<Particle> G_, Rg_, d_, s_, A_;
792  double G_val_, Rg_val_, d_val_, s_val_, A_val_, q1_param_, D_param_;
793 };
794 
795 IMPISD_END_NAMESPACE
796 
797 #endif /* IMPISD_UNIVARIATE_FUNCTIONS_H */
bool has_changed() const
return true if internal parameters have changed.
Base class for functions of one variable.
Floats operator()(const Floats &x) const
evaluate the function at a certain point
FloatsList get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist, bool) const
for testing purposes
FloatsList operator()(const FloatsList &xlist, bool) const
used for testing only
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
unsigned get_ndims_y() const
returns the number of output dimensions
IMP_Eigen::VectorXd operator()(const FloatsList &xlist) const
evaluate the function at a list of points
unsigned get_ndims_x() const
returns the number of input dimensions
A decorator for switching parameters particles.
A decorator for scale parameters particles.
A decorator for nuisance parameters particles.
void add_to_derivatives(const Floats &x, DerivativeAccumulator &accum) const
update derivatives of particles
Add scale parameter to particle.
Definition: Scale.h:24
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Ref counted objects should have private destructors.
bool get_particle_is_optimized(unsigned particle_no) const
returns true if the particle whose index is provided is optimized
unsigned get_number_of_particles() const
returns the number of particles that this function uses
#define IMP_LOG_TERSE(expr)
Definition: log_macros.h:83
A smart pointer to a reference counted object.
Definition: Pointer.h:87
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const
update derivatives of particles
#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
Linear one-dimensional function.
ModelObjectsTemp get_inputs() const
particle manipulation
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const
for testing purposes
IMP_Eigen::VectorXd get_second_derivative_vector(unsigned, unsigned, const FloatsList &xlist) const
return second derivative vector
Common base class for heavy weight IMP objects.
Definition: Object.h:106
unsigned get_number_of_particles() const
returns the number of particles that this function uses
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const
update derivatives of particles
FloatsList get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist, bool) const
for testing purposes
Add nuisance parameter to particle.
Definition: Nuisance.h:25
bool get_particle_is_optimized(unsigned particle_no) const
returns true if the particle whose index is provided is optimized
virtual IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const =0
return derivative vector
void update()
update internal parameters
unsigned get_number_of_optimized_particles() const
returns the number of particles that are optimized
unsigned get_number_of_optimized_particles() const
returns the number of particles that are optimized
FloatsList operator()(const FloatsList &xlist, bool) const
used for testing only
bool isnan(const T &a)
Return true if a number is NaN.
Definition: math.h:24
Floats operator()(const Floats &x) const
evaluate the function at a certain point
void add_to_derivatives(const Floats &x, DerivativeAccumulator &accum) const
update derivatives of particles
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
virtual void update()=0
update internal parameters
#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.
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
IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
ModelObjectsTemp get_inputs() const
particle manipulation
IMP_Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const
return second derivative vector
unsigned get_ndims_y() const
returns the number of output dimensions
bool has_changed() const
return true if internal parameters have changed.
IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
An exception which is thrown when the Model has attributes with invalid values.
Definition: exception.h:189
unsigned get_ndims_x() const
returns the number of input dimensions
virtual IMP_Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const =0
return second derivative vector
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const
for testing purposes
IMP_Eigen::VectorXd operator()(const FloatsList &xlist) const
evaluate the function at a list of points
Class for adding derivatives from restraints to the model.
void update()
update internal parameters