IMP  2.2.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-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);
117  a_val_ = Nuisance(a).get_nuisance();
118  b_val_ = Nuisance(b).get_nuisance();
119  update();
120  }
121 
122  bool has_changed() const {
123  double tmpa = Nuisance(a_).get_nuisance();
124  double tmpb = Nuisance(b_).get_nuisance();
125  if ((std::abs(tmpa - a_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
126  (std::abs(tmpb - b_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
127  IMP_LOG_TERSE("Linear1DFunction: has_changed():");
128  IMP_LOG_TERSE("true" << std::endl);
129  return true;
130  } else {
131  return false;
132  }
133  }
134 
135  void update() {
136  a_val_ = Nuisance(a_).get_nuisance();
137  b_val_ = Nuisance(b_).get_nuisance();
138  IMP_LOG_TERSE("Linear1DFunction: update() a:= " << a_val_ << " b:="
139  << b_val_ << std::endl);
140  }
141 
142  Floats operator()(const Floats& x) const {
143  IMP_USAGE_CHECK(x.size() == 1, "expecting a 1-D vector");
144  Floats ret(1, a_val_ * x[0] + b_val_);
145  return ret;
146  }
147 
148  IMP_Eigen::VectorXd operator()(const FloatsList& xlist) const {
149  unsigned M = xlist.size();
150  IMP_Eigen::VectorXd retlist(M);
151  for (unsigned i = 0; i < M; i++) {
152  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
153  retlist(i) = a_val_ * xlist[i][0] + b_val_;
154  }
155  return retlist;
156  }
157 
158  FloatsList operator()(const FloatsList& xlist, bool) const {
159  IMP_Eigen::VectorXd vec((*this)(xlist));
160  FloatsList ret;
161  for (unsigned i = 0; i < xlist.size(); i++)
162  ret.push_back(Floats(1, vec(i)));
163  return ret;
164  }
165 
166  void add_to_derivatives(const Floats& x, DerivativeAccumulator& accum) const {
167  // d[f(x)]/da = x
168  Nuisance(a_).add_to_nuisance_derivative(x[0], accum);
169  // d[f(x)]/db = 1
170  Nuisance(b_).add_to_nuisance_derivative(1, accum);
171  }
172 
173  void add_to_particle_derivative(unsigned particle_no, double value,
174  DerivativeAccumulator& accum) const {
175  switch (particle_no) {
176  case 0:
177  Nuisance(a_).add_to_nuisance_derivative(value, accum);
178  break;
179  case 1:
180  Nuisance(b_).add_to_nuisance_derivative(value, accum);
181  break;
182  default:
183  IMP_THROW("Invalid particle number", ModelException);
184  }
185  }
186 
187  IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no,
188  const FloatsList& xlist) const {
189  unsigned N = xlist.size();
190  IMP_Eigen::VectorXd ret(N);
191  switch (particle_no) {
192  case 0: // a
193  for (unsigned i = 0; i < N; i++) ret(i) = xlist[i][0];
194  break;
195  case 1: // b
196  ret.setOnes();
197  break;
198  default:
199  IMP_THROW("Invalid particle number", ModelException);
200  }
201  return ret;
202  }
203 
204  FloatsList get_derivative_matrix(const FloatsList& xlist, bool) const {
205  IMP_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 
217  IMP_Eigen::VectorXd get_second_derivative_vector(
218  unsigned, unsigned, const FloatsList& xlist) const {
219  // The Hessian is zero for all particles.
220  unsigned N = xlist.size();
221  IMP_Eigen::VectorXd H(IMP_Eigen::VectorXd::Zero(N));
222  return H;
223  }
224 
226  unsigned particle_b,
227  const FloatsList& xlist, bool) const {
228  IMP_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 { return 1; }
240  unsigned get_ndims_y() const { return 1; }
241 
242  unsigned get_number_of_particles() const { return 2; }
243 
244  bool get_particle_is_optimized(unsigned particle_no) const {
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 
264  ret.push_back(a_);
265  ret.push_back(b_);
266  return ret;
267  }
268 
270  ContainersTemp ret;
271  return ret;
272  }
273 
274  /*IMP_OBJECT_INLINE(Linear1DFunction, out << "y = " << a_val_
275  << " * x + " << b_val_ << std::endl, {});*/
277 
278  private:
280  double a_val_, b_val_;
281 };
282 
283 //! 1D mean function for SAS data
284 /* Generalized Guinier-Porod model (Hammouda, J. Appl. Cryst., 2010, eq 3 & 4
285  * to which a constant offset is added)
286  * I(q) = A + G/q^s exp(-(q.Rg)^2/(3-s)) for q <= q1
287  * I(q) = A + D/q^d for q > q1
288  * q1 = 1/Rg * ((d-s)(3-s)/2)^(1/2)
289  * D = G exp(-(q1.Rg)^2/(3-s)) q1^(d-s)
290  * only valid for q>0 Rg>0 0<s<d s<3 G>0
291  * s=0 in the globular case.
292  * G, Rg, d, s and A are particles (ISD Scales).
293  */
295  public:
298  kernel::Particle* A)
299  : UnivariateFunction("GeneralizedGuinierPorodFunction %1%"),
300  G_(G),
301  Rg_(Rg),
302  d_(d),
303  s_(s),
304  A_(A) {
305  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: constructor" << std::endl);
311  update();
312  }
313 
314  bool has_changed() const {
315  double tmpG = Scale(G_).get_scale();
316  double tmpRg = Scale(Rg_).get_scale();
317  double tmpd = Scale(d_).get_scale();
318  double tmps = Scale(s_).get_scale();
319  double tmpA = Nuisance(A_).get_nuisance();
320  if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
321  (std::abs(tmpRg - Rg_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
322  (std::abs(tmpd - d_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
323  (std::abs(tmps - s_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
324  (std::abs(tmpA - A_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
325  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: has_changed():");
326  IMP_LOG_TERSE("true" << std::endl);
327  return true;
328  } else {
329  return false;
330  }
331  }
332 
333  void update() {
334  G_val_ = Scale(G_).get_scale();
335  Rg_val_ = Scale(Rg_).get_scale();
336  d_val_ = Scale(d_).get_scale();
337  s_val_ = Scale(s_).get_scale();
338  if (d_val_ == s_val_) {
339  IMP_LOG_TERSE("Warning: d==s !" << std::endl);
340  if (s_val_ > 0.001) {
341  s_val_ -= 0.001;
342  } else {
343  d_val_ += 0.001;
344  }
345  // IMP_THROW("d == s ! ", ModelException);
346  }
347  A_val_ = Nuisance(A_).get_nuisance();
348  q1_param_ = std::sqrt((d_val_ - s_val_) * (3 - s_val_) / 2.);
349  D_param_ = G_val_ * std::exp(-IMP::square(q1_param_) / (3 - s_val_));
350  q1_param_ = q1_param_ / Rg_val_;
351  D_param_ *= std::pow(q1_param_, d_val_ - s_val_);
352  IMP_LOG_TERSE("GeneralizedGuinierPorodFunction: update() G:= "
353  << G_val_ << " Rg:=" << Rg_val_ << " d:=" << d_val_
354  << " s:=" << s_val_ << " A:=" << A_val_ << " Q1.Rg ="
355  << q1_param_ * Rg_val_ << " D =" << D_param_ << std::endl);
356  /*std::cout << "cpp"
357  << " 3:Q1 " << q1_param_
358  << " 5:D " << D_param_
359  << " 7:G " << G_val_
360  << " 9:Rg " << Rg_val_
361  << " 11:d " << d_val_
362  << " 13:s " << s_val_
363  << " 15:val " << get_value(0.1)
364  <<std::endl;*/
365  }
366 
367  Floats operator()(const Floats& x) const {
368  IMP_USAGE_CHECK(x.size() == 1, "expecting a 1-D vector");
369  Floats ret(1, get_value(x[0]));
370  return ret;
371  }
372 
373  IMP_Eigen::VectorXd operator()(const FloatsList& xlist) const {
374  unsigned M = xlist.size();
375  IMP_Eigen::VectorXd retlist(M);
376  for (unsigned i = 0; i < M; i++) {
377  IMP_USAGE_CHECK(xlist[i].size() == 1, "expecting a 1-D vector");
378  retlist(i) = get_value(xlist[i][0]);
379  }
380  return retlist;
381  }
382 
383  FloatsList operator()(const FloatsList& xlist, bool) const {
384  IMP_Eigen::VectorXd vec((*this)(xlist));
385  FloatsList ret;
386  for (unsigned i = 0; i < xlist.size(); i++)
387  ret.push_back(Floats(1, vec(i)));
388  return ret;
389  }
390 
391  void add_to_derivatives(const Floats& x, DerivativeAccumulator& accum) const {
392  double qval = x[0];
393  double value = get_value(qval) - A_val_;
394  double deriv;
395  // d[f(x)+A]/dG = f(x)/G
396  deriv = value / G_val_;
397  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for G is nan.");
398  Scale(G_).add_to_nuisance_derivative(deriv, accum);
399  if (qval <= q1_param_) {
400  // d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
401  deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
402  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for Rg is nan.");
403  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
404  // d[f(x)]/dd = 0
405  //
406  // d[f(x)]/ds = - f(x) * ( (q Rg / (3-s))^2 + log(q) )
407  deriv = -value *
408  (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
409  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for s is nan.");
410  Scale(s_).add_to_nuisance_derivative(deriv, accum);
411  } else {
412  // d[f(x)]/dRg = f(x) * (s-d)/Rg
413  deriv = value * (s_val_ - d_val_) / Rg_val_;
414  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for Rg is nan.");
415  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
416  // d[f(x)]/dd = f(x) * log(q1/q)
417  deriv = value * std::log(q1_param_ / qval);
418  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for d is nan.");
419  Scale(d_).add_to_nuisance_derivative(deriv, accum);
420  // d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
421  deriv = -value *
422  ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
423  IMP_INTERNAL_CHECK(!base::isnan(deriv), "derivative for d is nan.");
424  Scale(s_).add_to_nuisance_derivative(deriv, accum);
425  }
426  // d[f(x)+A]/dA = 1
427  deriv = 1;
428  Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
429  }
430 
431  void add_to_particle_derivative(unsigned particle_no, double value,
432  DerivativeAccumulator& accum) const {
433  switch (particle_no) {
434  case 0:
435  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for G is nan.");
436  Scale(G_).add_to_scale_derivative(value, accum);
437  break;
438  case 1:
439  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for Rg is nan.");
440  Scale(Rg_).add_to_scale_derivative(value, accum);
441  break;
442  case 2:
443  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for d is nan.");
444  Scale(d_).add_to_scale_derivative(value, accum);
445  break;
446  case 3:
447  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for s is nan.");
448  Scale(s_).add_to_scale_derivative(value, accum);
449  break;
450  case 4:
451  IMP_INTERNAL_CHECK(!base::isnan(value), "derivative for A is nan.");
452  Nuisance(A_).add_to_nuisance_derivative(value, accum);
453  break;
454  default:
455  IMP_THROW("Invalid particle number", ModelException);
456  }
457  }
458 
459  IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no,
460  const FloatsList& xlist) const {
461  unsigned N = xlist.size();
462  IMP_Eigen::VectorXd ret(N);
463  switch (particle_no) {
464  case 0: // G
465  // d[f(x)]/dG = f(x)/G
466  ret = ((*this)(xlist) - IMP_Eigen::VectorXd::Constant(N, A_val_)) /
467  G_val_;
468  break;
469  case 1: // Rg
470  for (unsigned i = 0; i < N; i++) {
471  double qval = xlist[i][0];
472  if (qval <= q1_param_) {
473  // d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
474  ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
475  Rg_val_ / (3 - s_val_);
476  } else {
477  // d[f(x)]/dRg = f(x) * (s-d)/Rg
478  ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
479  }
480  }
481  break;
482  case 2: // d
483  for (unsigned i = 0; i < N; i++) {
484  double qval = xlist[i][0];
485  if (qval <= q1_param_) {
486  // d[f(x)]/dd = 0
487  ret(i) = 0;
488  } else {
489  // d[f(x)]/dd = f(x) * log(q1/q)
490  ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
491  }
492  }
493  break;
494  case 3: // s
495  for (unsigned i = 0; i < N; i++) {
496  double qval = xlist[i][0];
497  if (qval <= q1_param_) {
498  // d[f(x)]/ds = - f(x)
499  // * (q^2 Rg^2 + (3-s)^2 log(q)) / (3-s)^2
500  ret(i) =
501  -(get_value(qval) - A_val_) *
502  (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
503  } else {
504  // d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
505  ret(i) =
506  -(get_value(qval) - A_val_) *
507  ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
508  }
509  }
510  break;
511  case 4: // A
512  ret = IMP_Eigen::VectorXd::Constant(N, 1);
513  break;
514  default:
515  IMP_THROW("Invalid particle number", ModelException);
516  }
517  return ret;
518  }
519 
520  FloatsList get_derivative_matrix(const FloatsList& xlist, bool) const {
521  IMP_Eigen::MatrixXd mat(xlist.size(), 5);
522  mat.col(0) = get_derivative_vector(0, xlist);
523  mat.col(1) = get_derivative_vector(1, xlist);
524  mat.col(2) = get_derivative_vector(2, xlist);
525  mat.col(3) = get_derivative_vector(3, xlist);
526  mat.col(4) = get_derivative_vector(4, xlist);
527  FloatsList ret;
528  for (int i = 0; i < mat.rows(); i++) {
529  Floats line;
530  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
531  ret.push_back(line);
532  }
533  return ret;
534  }
535 
536  IMP_Eigen::VectorXd get_second_derivative_vector(
537  unsigned particle_a, unsigned particle_b, const FloatsList& xlist) const {
538  if (particle_a >= 5) IMP_THROW("Invalid particle 1 number", ModelException);
539  if (particle_b >= 5) IMP_THROW("Invalid particle 2 number", ModelException);
540  unsigned N = xlist.size();
541  // hessian involving A is always 0
542  if (particle_a == 4 || particle_b == 4) return IMP_Eigen::VectorXd::Zero(N);
543  unsigned pa = std::min(particle_a, particle_b);
544  unsigned pb = std::max(particle_a, particle_b);
545  IMP_Eigen::VectorXd ret(N);
546  switch (pa) {
547  case 0: // G
548  switch (pb) {
549  case 0: // G
550  // d2f/dG2 = 0
551  ret.noalias() = IMP_Eigen::VectorXd::Zero(N);
552  break;
553 
554  case 1: // Rg
555  // d2f/dGdRg = df/dRg * 1/G
556  ret.noalias() = get_derivative_vector(1, xlist) / G_val_;
557  break;
558 
559  case 2: // d
560  // d2f/dGdd = df/dd * 1/G
561  ret.noalias() = get_derivative_vector(2, xlist) / G_val_;
562  break;
563 
564  case 3: // s
565  // d2f/dGds = df/ds * 1/G
566  ret.noalias() = get_derivative_vector(3, xlist) / G_val_;
567  break;
568 
569  default:
570  IMP_THROW("Invalid particle 2 number", ModelException);
571  break;
572  }
573  break;
574 
575  case 1: // Rg
576  switch (pb) {
577  case 1: // Rg
578  for (unsigned i = 0; i < N; i++) {
579  double qval = xlist[i][0];
580  if (qval <= q1_param_) {
581  // d2[f(x)]/dRg2 =
582  // f(x) * (2 q^2 / (3-s))
583  // * (2 q^2 Rg^2 / (3-s) - 1)
584  ret(i) = (get_value(qval) - A_val_) * 2 * IMP::square(qval) /
585  (3 - s_val_) *
586  (2 * IMP::square(qval * Rg_val_) / (3 - s_val_) - 1);
587  } else {
588  // d2[f(x)]/dRg2=f(x) * (d-s)/Rg^2 * (d-s+1)
589  ret(i) = (get_value(qval) - A_val_) * (d_val_ - s_val_) /
590  IMP::square(Rg_val_) * (d_val_ - s_val_ + 1);
591  }
592  }
593  break;
594 
595  case 2: // d
596  for (unsigned i = 0; i < N; i++) {
597  double qval = xlist[i][0];
598  if (qval <= q1_param_) {
599  // d2[f(x)]/dddRg = 0
600  ret(i) = 0;
601  } else {
602  // d2[f(x)]/dddRg = -f(x)/Rg
603  // - (d-s)/Rg *df(x)/dd
604  double val = (get_value(qval) - A_val_);
605  ret(i) = -val / Rg_val_ - (val * std::log(q1_param_ / qval) *
606  (d_val_ - s_val_) / Rg_val_);
607  }
608  }
609  break;
610 
611  case 3: // s
612  for (unsigned i = 0; i < N; i++) {
613  double qval = xlist[i][0];
614  double val = (get_value(qval) - A_val_);
615  if (qval <= q1_param_) {
616  // d2[f(x)]/dsdRg = -2q^2Rg/(3-s)
617  // * (df(x)/ds + f(x)/(3-s))
618  double deriv =
619  -val * (IMP::square((qval * Rg_val_) / (3 - s_val_)) +
620  std::log(qval));
621  ret(i) = -2 * IMP::square(qval) * Rg_val_ / (3 - s_val_) *
622  (deriv + val / (3 - s_val_));
623  } else {
624  // d2[f(x)]/dsdRg =
625  // 1/Rg * (f(x)- (d-s)*df(x)/ds)
626  double deriv = -val * ((d_val_ - s_val_) / (2 * (3 - s_val_)) +
627  std::log(q1_param_));
628  ret(i) = (val - (d_val_ - s_val_) * deriv) / Rg_val_;
629  }
630  }
631  break;
632 
633  default:
634  IMP_THROW("Invalid particle 2 number", ModelException);
635  break;
636  }
637  break;
638 
639  case 2: // d
640  switch (pb) {
641  case 2: // d
642  for (unsigned i = 0; i < N; i++) {
643  double qval = xlist[i][0];
644  if (qval <= q1_param_) {
645  // d2[f(x)]/dddd = 0
646  ret(i) = 0;
647  } else {
648  // d2[f(x)]/dddd =
649  // f(x)*(log(q1/q)^2 + 1/(2*(d-s)))
650  double val = (get_value(qval) - A_val_);
651  ret(i) = val * (IMP::square(std::log(q1_param_ / qval)) +
652  1 / (2 * (d_val_ - s_val_)));
653  }
654  }
655  break;
656 
657  case 3: // s
658  {
659  double lterm =
660  (d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_);
661  double rterm = 0.5 * (1 / (3 - s_val_) + 1 / (d_val_ - s_val_));
662  for (unsigned i = 0; i < N; i++) {
663  double qval = xlist[i][0];
664  if (qval <= q1_param_) {
665  // d2[f(x)]/ddds = 0
666  ret(i) = 0;
667  } else {
668  // d2[f(x)]/ddds = log(q1/q)*df(x)/ds
669  // - (1/(3-s) + 1/(d-s))*f(x)/2
670  //
671  double val = (get_value(qval) - A_val_);
672  ret(i) = -val * (std::log(q1_param_ / qval) * lterm + rterm);
673  }
674  }
675  } break;
676 
677  default:
678  IMP_THROW("Invalid particle 2 number", ModelException);
679  break;
680  }
681  break;
682 
683  case 3: // s
684  switch (pb) {
685  case 3: // s
686  {
687  double cterm =
688  IMP::square(0.5 * (d_val_ - s_val_) / (3 - s_val_) +
689  std::log(q1_param_)) +
690  0.5 * ((6 - s_val_ - d_val_) / IMP::square(3 - s_val_) +
691  1. / (d_val_ - s_val_));
692  for (unsigned i = 0; i < N; i++) {
693  double qval = xlist[i][0];
694  double val = (get_value(qval) - A_val_);
695  if (qval <= q1_param_) {
696  // d2[f(x)]/dsds = f(x) *
697  // ( [(qRg)^2/(3-s)^2 + log q ]^2
698  // - 2(qRg)^2/(3-s)^3 )
699  double factor = IMP::square((qval * Rg_val_) / (3 - s_val_));
700  ret(i) = val * (IMP::square(factor + std::log(qval)) -
701  2 * factor / (3 - s_val_));
702  } else {
703  // d2[f(x)]/dsds = f(x)
704  // * ( [ (d-s)/(2*(3-s)) + log(q1) ]^2
705  // + 1/2 * [ (6-s-d)/(3-s)^2
706  // + 1/(d-s) ] )
707  ret(i) = val * cterm;
708  }
709  }
710  } break;
711 
712  default:
713  IMP_THROW("Invalid particle 2 number", ModelException);
714  break;
715  }
716  break;
717 
718  default:
719  IMP_THROW("Invalid particle 1 number", ModelException);
720  break;
721  }
722  return ret;
723  }
724 
726  unsigned particle_b,
727  const FloatsList& xlist, bool) const {
728  IMP_Eigen::VectorXd mat(
729  get_second_derivative_vector(particle_a, particle_b, xlist));
730  FloatsList ret;
731  for (int i = 0; i < mat.rows(); i++) {
732  Floats line;
733  for (int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
734  ret.push_back(line);
735  }
736  return ret;
737  }
738 
739  unsigned get_ndims_x() const { return 1; }
740  unsigned get_ndims_y() const { return 1; }
741 
742  unsigned get_number_of_particles() const { return 5; }
743 
744  bool get_particle_is_optimized(unsigned particle_no) const {
745  switch (particle_no) {
746  case 0: // G
747  return Scale(G_).get_scale_is_optimized();
748  case 1: // Rg
749  return Scale(Rg_).get_scale_is_optimized();
750  case 2: // d
751  return Scale(d_).get_scale_is_optimized();
752  case 3: // s
753  return Scale(s_).get_scale_is_optimized();
754  case 4: // A
755  return Nuisance(A_).get_nuisance_is_optimized();
756  default:
757  IMP_THROW("Invalid particle number", ModelException);
758  }
759  }
760 
762  unsigned count = 0;
763  if (Scale(G_).get_scale_is_optimized()) count++;
764  if (Scale(Rg_).get_scale_is_optimized()) count++;
765  if (Scale(d_).get_scale_is_optimized()) count++;
766  if (Scale(s_).get_scale_is_optimized()) count++;
767  if (Nuisance(A_).get_nuisance_is_optimized()) count++;
768  return count;
769  }
770 
773  ret.push_back(G_);
774  ret.push_back(Rg_);
775  ret.push_back(d_);
776  ret.push_back(s_);
777  ret.push_back(A_);
778  return ret;
779  }
780 
782  ContainersTemp ret;
783  return ret;
784  }
785 
786  /*IMP_OBJECT_INLINE(GeneralizedGuinierPorodFunction, out
787  << " G = " << G_val_
788  << " Rg = " << Rg_val_
789  << " d = " << d_val_
790  << " s = " << s_val_
791  << " A = " << A_val_
792  << " Q1.Rg = " << q1_param_*Rg_val_
793  << std::endl, {});*/
795 
796  private:
797  inline double get_value(double q) const {
798  double value;
799  if (q <= q1_param_) {
800  value = A_val_ + G_val_ / std::pow(q, s_val_) *
801  std::exp(-IMP::square(q * Rg_val_) / (3 - s_val_));
802  } else {
803  value = A_val_ + D_param_ / std::pow(q, d_val_);
804  }
805  return value;
806  }
807 
808  base::Pointer<kernel::Particle> G_, Rg_, d_, s_, A_;
809  double G_val_, Rg_val_, d_val_, s_val_, A_val_, q1_param_, D_param_;
810 };
811 
812 IMPISD_END_NAMESPACE
813 
814 #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
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.
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Ref counted objects should have private destructors.
A smart pointer to a reference counted object.
Definition: base/Pointer.h:87
ParticlesTemp get_input_particles(const ModelObjectsTemp &mos)
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
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
static Scale decorate_particle(::IMP::kernel::Particle *p)
Definition: Scale.h:29
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const
update derivatives of particles
#define IMP_LOG_TERSE(expr)
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)
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
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
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
static Nuisance decorate_particle(::IMP::kernel::Particle *p)
Definition: Nuisance.h:30
virtual IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const =0
return derivative vector
void update()
update internal parameters
#define IMP_IF_CHECK(level)
Execute the code block if a certain level checks are on.
unsigned get_number_of_optimized_particles() const
returns the number of particles that are optimized
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
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: base/Object.h:106
Classes to handle individual model particles.
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: base/math.h:24
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
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
A shared base class to help in debugging and things.
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