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