IMP  2.3.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-2014 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/kernel/Particle.h>
13 #include <IMP/isd/Nuisance.h>
14 #include <IMP/isd/Scale.h>
15 #include <IMP/isd/Switching.h>
16 #include <IMP/base/Object.h>
17 #include <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 base::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 kernel::ParticlesTemp get_input_particles() const = 0;
102  virtual ContainersTemp get_input_containers() 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 
262  ret.push_back(a_);
263  ret.push_back(b_);
264  return ret;
265  }
266 
268  ContainersTemp ret;
269  return ret;
270  }
271 
272  /*IMP_OBJECT_INLINE(Linear1DFunction, out << "y = " << a_val_
273  << " * x + " << b_val_ << std::endl, {});*/
275 
276  private:
278  double a_val_, b_val_;
279 };
280 
281 //! 1D mean function for SAS data
282 /* Generalized Guinier-Porod model (Hammouda, J. Appl. Cryst., 2010, eq 3 & 4
283  * to which a constant offset is added)
284  * I(q) = A + G/q^s exp(-(q.Rg)^2/(3-s)) for q <= q1
285  * I(q) = A + D/q^d for q > q1
286  * q1 = 1/Rg * ((d-s)(3-s)/2)^(1/2)
287  * D = G exp(-(q1.Rg)^2/(3-s)) q1^(d-s)
288  * only valid for q>0 Rg>0 0<s<d s<3 G>0
289  * s=0 in the globular case.
290  * G, Rg, d, s and A are particles (ISD Scales).
291  */
293  public:
296  kernel::Particle* A)
297  : UnivariateFunction("GeneralizedGuinierPorodFunction %1%"),
298  G_(G),
299  Rg_(Rg),
300  d_(d),
301  s_(s),
302  A_(A) {
303  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: constructor" << std::endl);
304  update();
305  }
306 
307  bool has_changed() const {
308  double tmpG = Scale(G_).get_scale();
309  double tmpRg = Scale(Rg_).get_scale();
310  double tmpd = Scale(d_).get_scale();
311  double tmps = Scale(s_).get_scale();
312  double tmpA = Nuisance(A_).get_nuisance();
313  if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
314  (std::abs(tmpRg - Rg_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
315  (std::abs(tmpd - d_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
316  (std::abs(tmps - s_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
317  (std::abs(tmpA - A_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
318  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: has_changed():");
319  IMP_LOG_TERSE("true" << std::endl);
320  return true;
321  } else {
322  return false;
323  }
324  }
325 
326  void update() {
327  G_val_ = Scale(G_).get_scale();
328  Rg_val_ = Scale(Rg_).get_scale();
329  d_val_ = Scale(d_).get_scale();
330  s_val_ = Scale(s_).get_scale();
331  if (d_val_ == s_val_) {
332  IMP_LOG_TERSE("Warning: d==s !" << std::endl);
333  if (s_val_ > 0.001) {
334  s_val_ -= 0.001;
335  } else {
336  d_val_ += 0.001;
337  }
338  // IMP_THROW("d == s ! ", ModelException);
339  }
340  A_val_ = Nuisance(A_).get_nuisance();
341  q1_param_ = std::sqrt((d_val_ - s_val_) * (3 - s_val_) / 2.);
342  D_param_ = G_val_ * std::exp(-IMP::square(q1_param_) / (3 - s_val_));
343  q1_param_ = q1_param_ / Rg_val_;
344  D_param_ *= std::pow(q1_param_, d_val_ - s_val_);
345  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: update() G:= "
346  << G_val_ << " Rg:=" << Rg_val_ << " d:=" << d_val_
347  << " s:=" << s_val_ << " A:=" << A_val_ << " Q1.Rg ="
348  << q1_param_ * Rg_val_ << " D =" << D_param_ << std::endl);
349  /*std::cout << "cpp"
350  << " 3:Q1 " << q1_param_
351  << " 5:D " << D_param_
352  << " 7:G " << G_val_
353  << " 9:Rg " << Rg_val_
354  << " 11:d " << d_val_
355  << " 13:s " << s_val_
356  << " 15:val " << get_value(0.1)
357  <<std::endl;*/
358  }
359 
360  Floats operator()(const Floats& x) const {
361  IMP_USAGE_CHECK(x.size() == 1, "expecting a 1-D vector");
362  Floats ret(1, get_value(x[0]));
363  return ret;
364  }
365 
366  IMP_Eigen::VectorXd operator()(const FloatsList& xlist) const {
367  unsigned M = xlist.size();
368  IMP_Eigen::VectorXd retlist(M);
369  for (unsigned i = 0; i < M; i++) {
370  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
371  retlist(i) = get_value(xlist[i][0]);
372  }
373  return retlist;
374  }
375 
376  FloatsList operator()(const FloatsList& xlist, bool) const {
377  IMP_Eigen::VectorXd vec((*this)(xlist));
378  FloatsList ret;
379  for (unsigned i = 0; i < xlist.size(); i++)
380  ret.push_back(Floats(1, vec(i)));
381  return ret;
382  }
383 
384  void add_to_derivatives(const Floats& x, DerivativeAccumulator& accum) const {
385  double qval = x[0];
386  double value = get_value(qval) - A_val_;
387  double deriv;
388  // d[f(x)+A]/dG = f(x)/G
389  deriv = value / G_val_;
390  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for G is nan.");
391  Scale(G_).add_to_nuisance_derivative(deriv, accum);
392  if (qval <= q1_param_) {
393  // d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
394  deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
395  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for Rg is nan.");
396  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
397  // d[f(x)]/dd = 0
398  //
399  // d[f(x)]/ds = - f(x) * ( (q Rg / (3-s))^2 + log(q) )
400  deriv = -value *
401  (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
402  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for s is nan.");
403  Scale(s_).add_to_nuisance_derivative(deriv, accum);
404  } else {
405  // d[f(x)]/dRg = f(x) * (s-d)/Rg
406  deriv = value * (s_val_ - d_val_) / Rg_val_;
407  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for Rg is nan.");
408  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
409  // d[f(x)]/dd = f(x) * log(q1/q)
410  deriv = value * std::log(q1_param_ / qval);
411  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for d is nan.");
412  Scale(d_).add_to_nuisance_derivative(deriv, accum);
413  // d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
414  deriv = -value *
415  ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
416  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for d is nan.");
417  Scale(s_).add_to_nuisance_derivative(deriv, accum);
418  }
419  // d[f(x)+A]/dA = 1
420  deriv = 1;
421  Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
422  }
423 
424  void add_to_particle_derivative(unsigned particle_no, double value,
425  DerivativeAccumulator& accum) const {
426  switch (particle_no) {
427  case 0:
428  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for G is nan.");
429  Scale(G_).add_to_scale_derivative(value, accum);
430  break;
431  case 1:
432  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for Rg is nan.");
433  Scale(Rg_).add_to_scale_derivative(value, accum);
434  break;
435  case 2:
436  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for d is nan.");
437  Scale(d_).add_to_scale_derivative(value, accum);
438  break;
439  case 3:
440  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for s is nan.");
441  Scale(s_).add_to_scale_derivative(value, accum);
442  break;
443  case 4:
444  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for A is nan.");
445  Nuisance(A_).add_to_nuisance_derivative(value, accum);
446  break;
447  default:
448  IMP_THROW("Invalid particle number", ModelException);
449  }
450  }
451 
452  IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no,
453  const FloatsList& xlist) const {
454  unsigned N = xlist.size();
455  IMP_Eigen::VectorXd ret(N);
456  switch (particle_no) {
457  case 0: // G
458  // d[f(x)]/dG = f(x)/G
459  ret = ((*this)(xlist) - IMP_Eigen::VectorXd::Constant(N, A_val_)) /
460  G_val_;
461  break;
462  case 1: // Rg
463  for (unsigned i = 0; i < N; i++) {
464  double qval = xlist[i][0];
465  if (qval <= q1_param_) {
466  // d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
467  ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
468  Rg_val_ / (3 - s_val_);
469  } else {
470  // d[f(x)]/dRg = f(x) * (s-d)/Rg
471  ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
472  }
473  }
474  break;
475  case 2: // d
476  for (unsigned i = 0; i < N; i++) {
477  double qval = xlist[i][0];
478  if (qval <= q1_param_) {
479  // d[f(x)]/dd = 0
480  ret(i) = 0;
481  } else {
482  // d[f(x)]/dd = f(x) * log(q1/q)
483  ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
484  }
485  }
486  break;
487  case 3: // s
488  for (unsigned i = 0; i < N; i++) {
489  double qval = xlist[i][0];
490  if (qval <= q1_param_) {
491  // d[f(x)]/ds = - f(x)
492  // * (q^2 Rg^2 + (3-s)^2 log(q)) / (3-s)^2
493  ret(i) =
494  -(get_value(qval) - A_val_) *
495  (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
496  } else {
497  // d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
498  ret(i) =
499  -(get_value(qval) - A_val_) *
500  ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
501  }
502  }
503  break;
504  case 4: // A
505  ret = IMP_Eigen::VectorXd::Constant(N, 1);
506  break;
507  default:
508  IMP_THROW("Invalid particle number", ModelException);
509  }
510  return ret;
511  }
512 
513  FloatsList get_derivative_matrix(const FloatsList& xlist, bool) const {
514  IMP_Eigen::MatrixXd mat(xlist.size(), 5);
515  mat.col(0) = get_derivative_vector(0, xlist);
516  mat.col(1) = get_derivative_vector(1, xlist);
517  mat.col(2) = get_derivative_vector(2, xlist);
518  mat.col(3) = get_derivative_vector(3, xlist);
519  mat.col(4) = get_derivative_vector(4, xlist);
520  FloatsList ret;
521  for (int i = 0; i < mat.rows(); i++) {
522  Floats line;
523  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
524  ret.push_back(line);
525  }
526  return ret;
527  }
528 
529  IMP_Eigen::VectorXd get_second_derivative_vector(
530  unsigned particle_a, unsigned particle_b, const FloatsList& xlist) const {
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 IMP_Eigen::VectorXd::Zero(N);
536  unsigned pa = std::min(particle_a, particle_b);
537  unsigned pb = std::max(particle_a, particle_b);
538  IMP_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() = IMP_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_b,
720  const FloatsList& xlist, bool) const {
721  IMP_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 { return 1; }
733  unsigned get_ndims_y() const { return 1; }
734 
735  unsigned get_number_of_particles() const { return 5; }
736 
737  bool get_particle_is_optimized(unsigned particle_no) const {
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 
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 
775  ContainersTemp ret;
776  return ret;
777  }
778 
779  /*IMP_OBJECT_INLINE(GeneralizedGuinierPorodFunction, out
780  << " G = " << G_val_
781  << " Rg = " << Rg_val_
782  << " d = " << d_val_
783  << " s = " << s_val_
784  << " A = " << A_val_
785  << " Q1.Rg = " << q1_param_*Rg_val_
786  << std::endl, {});*/
788 
789  private:
790  inline double get_value(double q) const {
791  double value;
792  if (q <= q1_param_) {
793  value = A_val_ + G_val_ / std::pow(q, s_val_) *
794  std::exp(-IMP::square(q * Rg_val_) / (3 - s_val_));
795  } else {
796  value = A_val_ + D_param_ / std::pow(q, d_val_);
797  }
798  return value;
799  }
800 
801  base::Pointer<kernel::Particle> G_, Rg_, d_, s_, A_;
802  double G_val_, Rg_val_, d_val_, s_val_, A_val_, q1_param_, D_param_;
803 };
804 
805 IMPISD_END_NAMESPACE
806 
807 #endif /* IMPISD_UNIVARIATE_FUNCTIONS_H */
bool has_changed() const
return true if internal parameters have changed.
Base class for functions of one variable.
Class for adding derivatives from restraints to the model.
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 smart pointer to a reference counted object.
Definition: Pointer.h:87
ParticlesTemp get_input_particles(const ModelObjectsTemp &mos)
Return all the input particles for a given ModelObject.
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
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:141
Linear one-dimensional function.
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
ContainersTemp get_input_containers(const ModelObjectsTemp &mos)
Return all the input particles for a given ModelObject.
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
Object(std::string name)
Construct an object with the given name.
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
IMP::base::Vector< IMP::base::WeakPointer< Container > > ContainersTemp
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
Class to handle individual model particles.
FloatsList operator()(const FloatsList &xlist, bool) const
used for testing only
Common base class for heavy weight IMP objects.
Definition: Object.h:106
Classes to handle individual model particles. (Note that implementation of inline functions in in int...
Floats operator()(const Floats &x) const
evaluate the function at a certain point
kernel::ParticlesTemp get_input_particles() const
particle manipulation
void add_to_derivatives(const Floats &x, DerivativeAccumulator &accum) const
update derivatives of particles
virtual void update()=0
update internal parameters
bool isnan(const T &a)
Return true if a number is NaN.
Definition: math.h:24
#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.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:170
IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
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
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
kernel::ParticlesTemp get_input_particles() const
particle manipulation
void update()
update internal parameters