IMP logo
IMP Reference Guide  2.16.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-2021 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 <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 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 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 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 IMP_OVERRIDE {
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 
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 
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  Eigen::VectorXd operator()(const FloatsList& xlist) const IMP_OVERRIDE {
147  unsigned M = xlist.size();
148  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 IMP_OVERRIDE {
157  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,
165  DerivativeAccumulator& accum) const IMP_OVERRIDE {
166  // d[f(x)]/da = x
167  Nuisance(a_).add_to_nuisance_derivative(x[0], accum);
168  // d[f(x)]/db = 1
169  Nuisance(b_).add_to_nuisance_derivative(1, accum);
170  }
171 
172  void add_to_particle_derivative(unsigned particle_no, double value,
173  DerivativeAccumulator& accum) const IMP_OVERRIDE {
174  switch (particle_no) {
175  case 0:
176  Nuisance(a_).add_to_nuisance_derivative(value, accum);
177  break;
178  case 1:
179  Nuisance(b_).add_to_nuisance_derivative(value, accum);
180  break;
181  default:
182  IMP_THROW("Invalid particle number", ModelException);
183  }
184  }
185 
186  Eigen::VectorXd get_derivative_vector(unsigned particle_no,
187  const FloatsList& xlist) const IMP_OVERRIDE {
188  unsigned N = xlist.size();
189  Eigen::VectorXd ret(N);
190  switch (particle_no) {
191  case 0: // a
192  for (unsigned i = 0; i < N; i++) ret(i) = xlist[i][0];
193  break;
194  case 1: // b
195  ret.setOnes();
196  break;
197  default:
198  IMP_THROW("Invalid particle number", ModelException);
199  }
200  return ret;
201  }
202 
204  const FloatsList& xlist, bool) const IMP_OVERRIDE {
205  Eigen::MatrixXd mat(xlist.size(), 2);
206  mat.col(0) = get_derivative_vector(0, xlist);
207  mat.col(1) = get_derivative_vector(1, xlist);
208  FloatsList ret;
209  for (int i = 0; i < mat.rows(); i++) {
210  Floats line;
211  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
212  ret.push_back(line);
213  }
214  return ret;
215  }
216 
218  unsigned, unsigned, const FloatsList& xlist) const IMP_OVERRIDE {
219  // The Hessian is zero for all particles.
220  unsigned N = xlist.size();
221  Eigen::VectorXd H(Eigen::VectorXd::Zero(N));
222  return H;
223  }
224 
226  unsigned particle_a, unsigned particle_b,
227  const FloatsList& xlist, bool) const IMP_OVERRIDE {
228  Eigen::VectorXd mat(
229  get_second_derivative_vector(particle_a, particle_b, xlist));
230  FloatsList ret;
231  for (int i = 0; i < mat.rows(); i++) {
232  Floats line;
233  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
234  ret.push_back(line);
235  }
236  return ret;
237  }
238 
239  unsigned get_ndims_x() const IMP_OVERRIDE { return 1; }
240  unsigned get_ndims_y() const IMP_OVERRIDE { return 1; }
241 
242  unsigned get_number_of_particles() const IMP_OVERRIDE { return 2; }
243 
244  bool get_particle_is_optimized(unsigned particle_no) const IMP_OVERRIDE {
245  switch (particle_no) {
246  case 0: // a
247  return Nuisance(a_).get_nuisance_is_optimized();
248  case 1: // b
249  return Nuisance(b_).get_nuisance_is_optimized();
250  default:
251  IMP_THROW("Invalid particle number", ModelException);
252  }
253  }
254 
256  unsigned count = 0;
257  if (Nuisance(a_).get_nuisance_is_optimized()) count++;
258  if (Nuisance(b_).get_nuisance_is_optimized()) count++;
259  return count;
260  }
261 
263  ModelObjectsTemp ret;
264  ret.push_back(a_);
265  ret.push_back(b_);
266  return ret;
267  }
268 
269  /*IMP_OBJECT_INLINE(Linear1DFunction, out << "y = " << a_val_
270  << " * x + " << b_val_ << std::endl, {});*/
272 
273  private:
274  Pointer<Particle> a_, b_;
275  double a_val_, b_val_;
276 };
277 
278 //! 1D mean function for SAS data
279 /* Generalized Guinier-Porod model (Hammouda, J. Appl. Cryst., 2010, eq 3 & 4
280  * to which a constant offset is added)
281  * I(q) = A + G/q^s exp(-(q.Rg)^2/(3-s)) for q <= q1
282  * I(q) = A + D/q^d for q > q1
283  * q1 = 1/Rg * ((d-s)(3-s)/2)^(1/2)
284  * D = G exp(-(q1.Rg)^2/(3-s)) q1^(d-s)
285  * only valid for q>0 Rg>0 0<s<d s<3 G>0
286  * s=0 in the globular case.
287  * G, Rg, d, s and A are particles (ISD Scales).
288  */
290  public:
292  Particle* d, Particle* s,
293  Particle* A)
294  : UnivariateFunction("GeneralizedGuinierPorodFunction %1%"),
295  G_(G),
296  Rg_(Rg),
297  d_(d),
298  s_(s),
299  A_(A) {
300  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: constructor" << std::endl);
301  update();
302  }
303 
304  bool has_changed() const IMP_OVERRIDE {
305  double tmpG = Scale(G_).get_scale();
306  double tmpRg = Scale(Rg_).get_scale();
307  double tmpd = Scale(d_).get_scale();
308  double tmps = Scale(s_).get_scale();
309  double tmpA = Nuisance(A_).get_nuisance();
310  if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
311  (std::abs(tmpRg - Rg_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
312  (std::abs(tmpd - d_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
313  (std::abs(tmps - s_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
314  (std::abs(tmpA - A_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
315  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: has_changed():");
316  IMP_LOG_TERSE("true" << std::endl);
317  return true;
318  } else {
319  return false;
320  }
321  }
322 
324  G_val_ = Scale(G_).get_scale();
325  Rg_val_ = Scale(Rg_).get_scale();
326  d_val_ = Scale(d_).get_scale();
327  s_val_ = Scale(s_).get_scale();
328  if (d_val_ == s_val_) {
329  IMP_LOG_TERSE("Warning: d==s !" << std::endl);
330  if (s_val_ > 0.001) {
331  s_val_ -= 0.001;
332  } else {
333  d_val_ += 0.001;
334  }
335  // IMP_THROW("d == s ! ", ModelException);
336  }
337  A_val_ = Nuisance(A_).get_nuisance();
338  q1_param_ = std::sqrt((d_val_ - s_val_) * (3 - s_val_) / 2.);
339  D_param_ = G_val_ * std::exp(-IMP::square(q1_param_) / (3 - s_val_));
340  q1_param_ = q1_param_ / Rg_val_;
341  D_param_ *= std::pow(q1_param_, d_val_ - s_val_);
342  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: update() G:= "
343  << G_val_ << " Rg:=" << Rg_val_ << " d:=" << d_val_
344  << " s:=" << s_val_ << " A:=" << A_val_ << " Q1.Rg ="
345  << q1_param_ * Rg_val_ << " D =" << D_param_ << std::endl);
346  /*std::cout << "cpp"
347  << " 3:Q1 " << q1_param_
348  << " 5:D " << D_param_
349  << " 7:G " << G_val_
350  << " 9:Rg " << Rg_val_
351  << " 11:d " << d_val_
352  << " 13:s " << s_val_
353  << " 15:val " << get_value(0.1)
354  <<std::endl;*/
355  }
356 
358  IMP_USAGE_CHECK(x.size() == 1, "expecting a 1-D vector");
359  Floats ret(1, get_value(x[0]));
360  return ret;
361  }
362 
363  Eigen::VectorXd operator()(const FloatsList& xlist) const IMP_OVERRIDE {
364  unsigned M = xlist.size();
365  Eigen::VectorXd retlist(M);
366  for (unsigned i = 0; i < M; i++) {
367  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
368  retlist(i) = get_value(xlist[i][0]);
369  }
370  return retlist;
371  }
372 
373  FloatsList operator()(const FloatsList& xlist, bool) const IMP_OVERRIDE {
374  Eigen::VectorXd vec((*this)(xlist));
375  FloatsList ret;
376  for (unsigned i = 0; i < xlist.size(); i++)
377  ret.push_back(Floats(1, vec(i)));
378  return ret;
379  }
380 
382  const Floats& x, DerivativeAccumulator& accum) const IMP_OVERRIDE {
383  double qval = x[0];
384  double value = get_value(qval) - A_val_;
385  double deriv;
386  // d[f(x)+A]/dG = f(x)/G
387  deriv = value / G_val_;
388  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for G is nan.");
389  Scale(G_).add_to_nuisance_derivative(deriv, accum);
390  if (qval <= q1_param_) {
391  // d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
392  deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
393  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for Rg is nan.");
394  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
395  // d[f(x)]/dd = 0
396  //
397  // d[f(x)]/ds = - f(x) * ( (q Rg / (3-s))^2 + log(q) )
398  deriv = -value *
399  (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
400  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for s is nan.");
401  Scale(s_).add_to_nuisance_derivative(deriv, accum);
402  } else {
403  // d[f(x)]/dRg = f(x) * (s-d)/Rg
404  deriv = value * (s_val_ - d_val_) / Rg_val_;
405  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for Rg is nan.");
406  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
407  // d[f(x)]/dd = f(x) * log(q1/q)
408  deriv = value * std::log(q1_param_ / qval);
409  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for d is nan.");
410  Scale(d_).add_to_nuisance_derivative(deriv, accum);
411  // d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
412  deriv = -value *
413  ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
414  IMP_INTERNAL_CHECK(!IMP::isnan(deriv), "derivative for d is nan.");
415  Scale(s_).add_to_nuisance_derivative(deriv, accum);
416  }
417  // d[f(x)+A]/dA = 1
418  deriv = 1;
419  Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
420  }
421 
422  void add_to_particle_derivative(unsigned particle_no, double value,
423  DerivativeAccumulator& accum) const IMP_OVERRIDE {
424  switch (particle_no) {
425  case 0:
426  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for G is nan.");
427  Scale(G_).add_to_scale_derivative(value, accum);
428  break;
429  case 1:
430  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for Rg is nan.");
431  Scale(Rg_).add_to_scale_derivative(value, accum);
432  break;
433  case 2:
434  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for d is nan.");
435  Scale(d_).add_to_scale_derivative(value, accum);
436  break;
437  case 3:
438  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for s is nan.");
439  Scale(s_).add_to_scale_derivative(value, accum);
440  break;
441  case 4:
442  IMP_INTERNAL_CHECK(!IMP::isnan(value), "derivative for A is nan.");
443  Nuisance(A_).add_to_nuisance_derivative(value, accum);
444  break;
445  default:
446  IMP_THROW("Invalid particle number", ModelException);
447  }
448  }
449 
450  Eigen::VectorXd get_derivative_vector(unsigned particle_no,
451  const FloatsList& xlist) const IMP_OVERRIDE {
452  unsigned N = xlist.size();
453  Eigen::VectorXd ret(N);
454  switch (particle_no) {
455  case 0: // G
456  // d[f(x)]/dG = f(x)/G
457  ret = ((*this)(xlist) - Eigen::VectorXd::Constant(N, A_val_)) /
458  G_val_;
459  break;
460  case 1: // Rg
461  for (unsigned i = 0; i < N; i++) {
462  double qval = xlist[i][0];
463  if (qval <= q1_param_) {
464  // d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
465  ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
466  Rg_val_ / (3 - s_val_);
467  } else {
468  // d[f(x)]/dRg = f(x) * (s-d)/Rg
469  ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
470  }
471  }
472  break;
473  case 2: // d
474  for (unsigned i = 0; i < N; i++) {
475  double qval = xlist[i][0];
476  if (qval <= q1_param_) {
477  // d[f(x)]/dd = 0
478  ret(i) = 0;
479  } else {
480  // d[f(x)]/dd = f(x) * log(q1/q)
481  ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
482  }
483  }
484  break;
485  case 3: // s
486  for (unsigned i = 0; i < N; i++) {
487  double qval = xlist[i][0];
488  if (qval <= q1_param_) {
489  // d[f(x)]/ds = - f(x)
490  // * (q^2 Rg^2 + (3-s)^2 log(q)) / (3-s)^2
491  ret(i) =
492  -(get_value(qval) - A_val_) *
493  (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
494  } else {
495  // d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
496  ret(i) =
497  -(get_value(qval) - A_val_) *
498  ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
499  }
500  }
501  break;
502  case 4: // A
503  ret = Eigen::VectorXd::Constant(N, 1);
504  break;
505  default:
506  IMP_THROW("Invalid particle number", ModelException);
507  }
508  return ret;
509  }
510 
512  const FloatsList& xlist, bool) const IMP_OVERRIDE {
513  Eigen::MatrixXd mat(xlist.size(), 5);
514  mat.col(0) = get_derivative_vector(0, xlist);
515  mat.col(1) = get_derivative_vector(1, xlist);
516  mat.col(2) = get_derivative_vector(2, xlist);
517  mat.col(3) = get_derivative_vector(3, xlist);
518  mat.col(4) = get_derivative_vector(4, xlist);
519  FloatsList ret;
520  for (int i = 0; i < mat.rows(); i++) {
521  Floats line;
522  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
523  ret.push_back(line);
524  }
525  return ret;
526  }
527 
529  unsigned particle_a, unsigned particle_b,
530  const FloatsList& xlist) const IMP_OVERRIDE {
531  if (particle_a >= 5) IMP_THROW("Invalid particle 1 number", ModelException);
532  if (particle_b >= 5) IMP_THROW("Invalid particle 2 number", ModelException);
533  unsigned N = xlist.size();
534  // hessian involving A is always 0
535  if (particle_a == 4 || particle_b == 4) return Eigen::VectorXd::Zero(N);
536  unsigned pa = std::min(particle_a, particle_b);
537  unsigned pb = std::max(particle_a, particle_b);
538  Eigen::VectorXd ret(N);
539  switch (pa) {
540  case 0: // G
541  switch (pb) {
542  case 0: // G
543  // d2f/dG2 = 0
544  ret.noalias() = Eigen::VectorXd::Zero(N);
545  break;
546 
547  case 1: // Rg
548  // d2f/dGdRg = df/dRg * 1/G
549  ret.noalias() = get_derivative_vector(1, xlist) / G_val_;
550  break;
551 
552  case 2: // d
553  // d2f/dGdd = df/dd * 1/G
554  ret.noalias() = get_derivative_vector(2, xlist) / G_val_;
555  break;
556 
557  case 3: // s
558  // d2f/dGds = df/ds * 1/G
559  ret.noalias() = get_derivative_vector(3, xlist) / G_val_;
560  break;
561 
562  default:
563  IMP_THROW("Invalid particle 2 number", ModelException);
564  break;
565  }
566  break;
567 
568  case 1: // Rg
569  switch (pb) {
570  case 1: // Rg
571  for (unsigned i = 0; i < N; i++) {
572  double qval = xlist[i][0];
573  if (qval <= q1_param_) {
574  // d2[f(x)]/dRg2 =
575  // f(x) * (2 q^2 / (3-s))
576  // * (2 q^2 Rg^2 / (3-s) - 1)
577  ret(i) = (get_value(qval) - A_val_) * 2 * IMP::square(qval) /
578  (3 - s_val_) *
579  (2 * IMP::square(qval * Rg_val_) / (3 - s_val_) - 1);
580  } else {
581  // d2[f(x)]/dRg2=f(x) * (d-s)/Rg^2 * (d-s+1)
582  ret(i) = (get_value(qval) - A_val_) * (d_val_ - s_val_) /
583  IMP::square(Rg_val_) * (d_val_ - s_val_ + 1);
584  }
585  }
586  break;
587 
588  case 2: // d
589  for (unsigned i = 0; i < N; i++) {
590  double qval = xlist[i][0];
591  if (qval <= q1_param_) {
592  // d2[f(x)]/dddRg = 0
593  ret(i) = 0;
594  } else {
595  // d2[f(x)]/dddRg = -f(x)/Rg
596  // - (d-s)/Rg *df(x)/dd
597  double val = (get_value(qval) - A_val_);
598  ret(i) = -val / Rg_val_ - (val * std::log(q1_param_ / qval) *
599  (d_val_ - s_val_) / Rg_val_);
600  }
601  }
602  break;
603 
604  case 3: // s
605  for (unsigned i = 0; i < N; i++) {
606  double qval = xlist[i][0];
607  double val = (get_value(qval) - A_val_);
608  if (qval <= q1_param_) {
609  // d2[f(x)]/dsdRg = -2q^2Rg/(3-s)
610  // * (df(x)/ds + f(x)/(3-s))
611  double deriv =
612  -val * (IMP::square((qval * Rg_val_) / (3 - s_val_)) +
613  std::log(qval));
614  ret(i) = -2 * IMP::square(qval) * Rg_val_ / (3 - s_val_) *
615  (deriv + val / (3 - s_val_));
616  } else {
617  // d2[f(x)]/dsdRg =
618  // 1/Rg * (f(x)- (d-s)*df(x)/ds)
619  double deriv = -val * ((d_val_ - s_val_) / (2 * (3 - s_val_)) +
620  std::log(q1_param_));
621  ret(i) = (val - (d_val_ - s_val_) * deriv) / Rg_val_;
622  }
623  }
624  break;
625 
626  default:
627  IMP_THROW("Invalid particle 2 number", ModelException);
628  break;
629  }
630  break;
631 
632  case 2: // d
633  switch (pb) {
634  case 2: // d
635  for (unsigned i = 0; i < N; i++) {
636  double qval = xlist[i][0];
637  if (qval <= q1_param_) {
638  // d2[f(x)]/dddd = 0
639  ret(i) = 0;
640  } else {
641  // d2[f(x)]/dddd =
642  // f(x)*(log(q1/q)^2 + 1/(2*(d-s)))
643  double val = (get_value(qval) - A_val_);
644  ret(i) = val * (IMP::square(std::log(q1_param_ / qval)) +
645  1 / (2 * (d_val_ - s_val_)));
646  }
647  }
648  break;
649 
650  case 3: // s
651  {
652  double lterm =
653  (d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_);
654  double rterm = 0.5 * (1 / (3 - s_val_) + 1 / (d_val_ - s_val_));
655  for (unsigned i = 0; i < N; i++) {
656  double qval = xlist[i][0];
657  if (qval <= q1_param_) {
658  // d2[f(x)]/ddds = 0
659  ret(i) = 0;
660  } else {
661  // d2[f(x)]/ddds = log(q1/q)*df(x)/ds
662  // - (1/(3-s) + 1/(d-s))*f(x)/2
663  //
664  double val = (get_value(qval) - A_val_);
665  ret(i) = -val * (std::log(q1_param_ / qval) * lterm + rterm);
666  }
667  }
668  } break;
669 
670  default:
671  IMP_THROW("Invalid particle 2 number", ModelException);
672  break;
673  }
674  break;
675 
676  case 3: // s
677  switch (pb) {
678  case 3: // s
679  {
680  double cterm =
681  IMP::square(0.5 * (d_val_ - s_val_) / (3 - s_val_) +
682  std::log(q1_param_)) +
683  0.5 * ((6 - s_val_ - d_val_) / IMP::square(3 - s_val_) +
684  1. / (d_val_ - s_val_));
685  for (unsigned i = 0; i < N; i++) {
686  double qval = xlist[i][0];
687  double val = (get_value(qval) - A_val_);
688  if (qval <= q1_param_) {
689  // d2[f(x)]/dsds = f(x) *
690  // ( [(qRg)^2/(3-s)^2 + log q ]^2
691  // - 2(qRg)^2/(3-s)^3 )
692  double factor = IMP::square((qval * Rg_val_) / (3 - s_val_));
693  ret(i) = val * (IMP::square(factor + std::log(qval)) -
694  2 * factor / (3 - s_val_));
695  } else {
696  // d2[f(x)]/dsds = f(x)
697  // * ( [ (d-s)/(2*(3-s)) + log(q1) ]^2
698  // + 1/2 * [ (6-s-d)/(3-s)^2
699  // + 1/(d-s) ] )
700  ret(i) = val * cterm;
701  }
702  }
703  } break;
704 
705  default:
706  IMP_THROW("Invalid particle 2 number", ModelException);
707  break;
708  }
709  break;
710 
711  default:
712  IMP_THROW("Invalid particle 1 number", ModelException);
713  break;
714  }
715  return ret;
716  }
717 
719  unsigned particle_a, unsigned particle_b,
720  const FloatsList& xlist, bool) const IMP_OVERRIDE {
721  Eigen::VectorXd mat(
722  get_second_derivative_vector(particle_a, particle_b, xlist));
723  FloatsList ret;
724  for (int i = 0; i < mat.rows(); i++) {
725  Floats line;
726  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
727  ret.push_back(line);
728  }
729  return ret;
730  }
731 
732  unsigned get_ndims_x() const IMP_OVERRIDE { return 1; }
733  unsigned get_ndims_y() const IMP_OVERRIDE { return 1; }
734 
735  unsigned get_number_of_particles() const IMP_OVERRIDE { return 5; }
736 
737  bool get_particle_is_optimized(unsigned particle_no) const IMP_OVERRIDE {
738  switch (particle_no) {
739  case 0: // G
740  return Scale(G_).get_scale_is_optimized();
741  case 1: // Rg
742  return Scale(Rg_).get_scale_is_optimized();
743  case 2: // d
744  return Scale(d_).get_scale_is_optimized();
745  case 3: // s
746  return Scale(s_).get_scale_is_optimized();
747  case 4: // A
748  return Nuisance(A_).get_nuisance_is_optimized();
749  default:
750  IMP_THROW("Invalid particle number", ModelException);
751  }
752  }
753 
755  unsigned count = 0;
756  if (Scale(G_).get_scale_is_optimized()) count++;
757  if (Scale(Rg_).get_scale_is_optimized()) count++;
758  if (Scale(d_).get_scale_is_optimized()) count++;
759  if (Scale(s_).get_scale_is_optimized()) count++;
760  if (Nuisance(A_).get_nuisance_is_optimized()) count++;
761  return count;
762  }
763 
765  ModelObjectsTemp ret;
766  ret.push_back(G_);
767  ret.push_back(Rg_);
768  ret.push_back(d_);
769  ret.push_back(s_);
770  ret.push_back(A_);
771  return ret;
772  }
773 
774  /*IMP_OBJECT_INLINE(GeneralizedGuinierPorodFunction, out
775  << " G = " << G_val_
776  << " Rg = " << Rg_val_
777  << " d = " << d_val_
778  << " s = " << s_val_
779  << " A = " << A_val_
780  << " Q1.Rg = " << q1_param_*Rg_val_
781  << std::endl, {});*/
783 
784  private:
785  inline double get_value(double q) const {
786  double value;
787  if (q <= q1_param_) {
788  value = A_val_ + G_val_ / std::pow(q, s_val_) *
789  std::exp(-IMP::square(q * Rg_val_) / (3 - s_val_));
790  } else {
791  value = A_val_ + D_param_ / std::pow(q, d_val_);
792  }
793  return value;
794  }
795 
796  Pointer<Particle> G_, Rg_, d_, s_, A_;
797  double G_val_, Rg_val_, d_val_, s_val_, A_val_, q1_param_, D_param_;
798 };
799 
800 IMPISD_END_NAMESPACE
801 
802 #endif /* IMPISD_UNIVARIATE_FUNCTIONS_H */
bool has_changed() const
return true if internal parameters have changed.
Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
Base class for functions of one variable.
virtual Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const =0
return second derivative vector
Floats operator()(const Floats &x) const
evaluate the function at a certain point
Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const
return second derivative vector
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
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:72
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
Eigen::VectorXd get_second_derivative_vector(unsigned, unsigned, const FloatsList &xlist) const
return second derivative vector
Linear one-dimensional function.
virtual Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const =0
return derivative vector
ModelObjectsTemp get_inputs() const
particle manipulation
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const
for testing purposes
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
Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
bool get_particle_is_optimized(unsigned particle_no) const
returns true if the particle whose index is provided is optimized
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
Eigen::VectorXd operator()(const FloatsList &xlist) const
evaluate the function at a list of points
ModelObjectsTemp get_inputs() const
particle manipulation
unsigned get_ndims_y() const
returns the number of output dimensions
bool has_changed() const
return true if internal parameters have changed.
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
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const
for testing purposes
Eigen::VectorXd operator()(const FloatsList &xlist) const
evaluate the function at a list of points
#define IMP_OVERRIDE
Cause a compile error if this method does not override a parent method.
Class for adding derivatives from restraints to the model.
void update()
update internal parameters