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