IMP  2.0.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-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/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 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(Particle * a, 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  }
124 
125  bool has_changed() const {
126  double tmpa = Nuisance(a_).get_nuisance();
127  double tmpb = Nuisance(b_).get_nuisance();
128  if ((std::abs(tmpa - a_val_) >
129  IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)
130  || (std::abs(tmpb - b_val_) >
131  IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM))
132  {
133  IMP_LOG_TERSE( "Linear1DFunction: has_changed():");
134  IMP_LOG_TERSE( "true" << std::endl);
135  return true;
136  } else {
137  return false;
138  }
139  }
140 
141  void update() {
142  a_val_ = Nuisance(a_).get_nuisance();
143  b_val_ = Nuisance(b_).get_nuisance();
144  IMP_LOG_TERSE( "Linear1DFunction: update() a:= "
145  << a_val_ << " b:=" << b_val_ << std::endl);
146  }
147 
148  Floats operator()(const Floats& x) const {
149  IMP_USAGE_CHECK(x.size() == 1, "expecting a 1-D vector");
150  Floats ret(1,a_val_*x[0]+b_val_);
151  return ret;
152  }
153 
154  Eigen::VectorXd operator()(const FloatsList& xlist) const
155  {
156  unsigned M = xlist.size();
157  Eigen::VectorXd retlist(M);
158  for (unsigned i = 0; i < M; i++)
159  {
160  IMP_USAGE_CHECK(xlist[i].size() == 1,
161  "expecting a 1-D vector");
162  retlist(i) = a_val_*xlist[i][0]+b_val_;
163  }
164  return retlist;
165  }
166 
167  FloatsList operator()(const FloatsList& xlist, bool) const
168  {
169  Eigen::VectorXd vec((*this)(xlist));
170  FloatsList ret;
171  for (unsigned i=0; i<xlist.size(); i++)
172  ret.push_back(Floats(1,vec(i)));
173  return ret;
174  }
175 
176  void add_to_derivatives(const Floats& x,
177  DerivativeAccumulator &accum) const
178  {
179  //d[f(x)]/da = x
180  Nuisance(a_).add_to_nuisance_derivative(x[0], accum);
181  //d[f(x)]/db = 1
182  Nuisance(b_).add_to_nuisance_derivative(1, accum);
183  }
184 
185  void add_to_particle_derivative(unsigned particle_no,
186  double value, DerivativeAccumulator &accum) const
187  {
188  switch (particle_no)
189  {
190  case 0:
191  Nuisance(a_).add_to_nuisance_derivative(value, accum);
192  break;
193  case 1:
194  Nuisance(b_).add_to_nuisance_derivative(value, accum);
195  break;
196  default:
197  IMP_THROW("Invalid particle number", ModelException);
198  }
199  }
200 
201  Eigen::VectorXd get_derivative_vector(
202  unsigned particle_no, const FloatsList& xlist) const
203  {
204  unsigned N=xlist.size();
205  Eigen::VectorXd ret(N);
206  switch(particle_no)
207  {
208  case 0: // a
209  for (unsigned i=0; i<N; i++)
210  ret(i) = xlist[i][0];
211  break;
212  case 1: // b
213  ret.setOnes();
214  break;
215  default:
216  IMP_THROW("Invalid particle number", ModelException);
217  }
218  return ret;
219  }
220 
221  FloatsList get_derivative_matrix(
222  const FloatsList& xlist, bool) const
223  {
224  Eigen::MatrixXd mat(xlist.size(),2);
225  mat.col(0) = get_derivative_vector(0, xlist);
226  mat.col(1) = get_derivative_vector(1, xlist);
227  FloatsList ret;
228  for (int i=0; i<mat.rows(); i++)
229  {
230  Floats line;
231  for (int j=0; j<mat.cols(); j++)
232  line.push_back(mat(i,j));
233  ret.push_back(line);
234  }
235  return ret;
236  }
237 
238  Eigen::VectorXd get_second_derivative_vector(
239  unsigned, unsigned, const FloatsList& xlist) const
240  {
241  // The Hessian is zero for all particles.
242  unsigned N=xlist.size();
243  Eigen::VectorXd H(Eigen::VectorXd::Zero(N));
244  return H;
245  }
246 
247  FloatsList get_second_derivative_vector(
248  unsigned particle_a, unsigned particle_b,
249  const FloatsList& xlist, bool) const
250  {
251  Eigen::VectorXd mat( get_second_derivative_vector(
252  particle_a, particle_b, xlist));
253  FloatsList ret;
254  for (int i=0; i<mat.rows(); i++)
255  {
256  Floats line;
257  for (int j=0; j<mat.cols(); j++)
258  line.push_back(mat(i,j));
259  ret.push_back(line);
260  }
261  return ret;
262  }
263 
264  unsigned get_ndims_x() const {return 1;}
265  unsigned get_ndims_y() const {return 1;}
266 
267  unsigned get_number_of_particles() const { return 2; }
268 
269  bool get_particle_is_optimized(unsigned particle_no) const
270  {
271  switch(particle_no)
272  {
273  case 0: //a
274  return Nuisance(a_).get_nuisance_is_optimized();
275  case 1: //b
276  return Nuisance(b_).get_nuisance_is_optimized();
277  default:
278  IMP_THROW("Invalid particle number", ModelException);
279  }
280  }
281 
282  unsigned get_number_of_optimized_particles() const
283  {
284  unsigned count = 0;
285  if (Nuisance(a_).get_nuisance_is_optimized()) count++;
286  if (Nuisance(b_).get_nuisance_is_optimized()) count++;
287  return count;
288  }
289 
290  ParticlesTemp get_input_particles() const
291  {
292  ParticlesTemp ret;
293  ret.push_back(a_);
294  ret.push_back(b_);
295  return ret;
296  }
297 
298  ContainersTemp get_input_containers() const
299  {
300  ContainersTemp ret;
301  return ret;
302  }
303 
304  IMP_OBJECT_INLINE(Linear1DFunction, out << "y = " << a_val_
305  << " * x + " << b_val_ << std::endl, {});
306 
307  private:
308  Pointer<Particle> a_,b_;
309  double a_val_, b_val_;
310 
311 };
312 
313 //! 1D mean function for SAS data
314 /* Generalized Guinier-Porod model (Hammouda, J. Appl. Cryst., 2010, eq 3 & 4
315  * to which a constant offset is added)
316  * I(q) = A + G/q^s exp(-(q.Rg)^2/(3-s)) for q <= q1
317  * I(q) = A + D/q^d for q > q1
318  * q1 = 1/Rg * ((d-s)(3-s)/2)^(1/2)
319  * D = G exp(-(q1.Rg)^2/(3-s)) q1^(d-s)
320  * only valid for q>0 Rg>0 0<s<d s<3 G>0
321  * s=0 in the globular case.
322  * G, Rg, d, s and A are particles (ISD Scales).
323  */
324 class IMPISDEXPORT GeneralizedGuinierPorodFunction : public UnivariateFunction
325 {
326  public:
327  GeneralizedGuinierPorodFunction(Particle * G, Particle * Rg,
328  Particle * d, Particle * s, Particle * A)
329  : UnivariateFunction("GeneralizedGuinierPorodFunction %1%"),
330  G_(G), Rg_(Rg), d_(d), s_(s), A_(A)
331  {
332  IMP_LOG_TERSE( "GeneralizedGuinierPorodFunction: constructor"
333  << std::endl);
334  IMP_IF_CHECK(USAGE_AND_INTERNAL) { Scale::decorate_particle(G); }
335  IMP_IF_CHECK(USAGE_AND_INTERNAL) { Scale::decorate_particle(Rg); }
336  IMP_IF_CHECK(USAGE_AND_INTERNAL) { Scale::decorate_particle(d); }
337  IMP_IF_CHECK(USAGE_AND_INTERNAL) { Scale::decorate_particle(s); }
338  IMP_IF_CHECK(USAGE_AND_INTERNAL) { Nuisance::decorate_particle(A); }
339  update();
340  }
341 
342  bool has_changed() const {
343  double tmpG = Scale(G_).get_scale();
344  double tmpRg = Scale(Rg_).get_scale();
345  double tmpd = Scale(d_).get_scale();
346  double tmps = Scale(s_).get_scale();
347  double tmpA = Nuisance(A_).get_nuisance();
348  if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)
349  || (std::abs(tmpRg - Rg_val_) >
350  IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)
351  || (std::abs(tmpd - d_val_) >
352  IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)
353  || (std::abs(tmps - s_val_) >
354  IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)
355  || (std::abs(tmpA - A_val_) >
356  IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM))
357  {
358  IMP_LOG_TERSE(
359  "GeneralizedGuinierPorodFunction: has_changed():");
360  IMP_LOG_TERSE( "true" << std::endl);
361  return true;
362  } else {
363  return false;
364  }
365  }
366 
367  void update() {
368  G_val_ = Scale(G_).get_scale();
369  Rg_val_ = Scale(Rg_).get_scale();
370  d_val_ = Scale(d_).get_scale();
371  s_val_ = Scale(s_).get_scale();
372  A_val_ = Nuisance(A_).get_nuisance();
373  q1_param_ = std::sqrt((d_val_-s_val_)*(3-s_val_)/2.);
374  D_param_ = G_val_ *std::exp(-IMP::square(q1_param_)/(3-s_val_));
375  q1_param_ = q1_param_ / Rg_val_;
376  D_param_ *= std::pow(q1_param_,d_val_-s_val_);
377  IMP_LOG_TERSE( "GeneralizedGuinierPorodFunction: update() G:= "
378  << G_val_
379  << " Rg:=" << Rg_val_
380  << " d:=" << d_val_
381  << " s:=" << s_val_
382  << " A:=" << A_val_
383  << " Q1.Rg =" << q1_param_*Rg_val_
384  << " D =" << D_param_
385  << std::endl);
386  /*std::cout << "cpp"
387  << " 3:Q1 " << q1_param_
388  << " 5:D " << D_param_
389  << " 7:G " << G_val_
390  << " 9:Rg " << Rg_val_
391  << " 11:d " << d_val_
392  << " 13:s " << s_val_
393  << " 15:val " << get_value(0.1)
394  <<std::endl;*/
395  }
396 
397  Floats operator()(const Floats& x) const {
398  IMP_USAGE_CHECK(x.size() == 1, "expecting a 1-D vector");
399  Floats ret(1,get_value(x[0]));
400  return ret;
401  }
402 
403  Eigen::VectorXd operator()(const FloatsList& xlist) const
404  {
405  unsigned M=xlist.size();
406  Eigen::VectorXd retlist(M);
407  for (unsigned i = 0; i < M; i++)
408  {
409  IMP_USAGE_CHECK(xlist[i].size() == 1,
410  "expecting a 1-D vector");
411  retlist(i) = get_value(xlist[i][0]);
412  }
413  return retlist;
414  }
415 
416  FloatsList operator()(const FloatsList& xlist, bool) const
417  {
418  Eigen::VectorXd vec((*this)(xlist));
419  FloatsList ret;
420  for (unsigned i=0; i<xlist.size(); i++)
421  ret.push_back(Floats(1,vec(i)));
422  return ret;
423  }
424 
425  void add_to_derivatives(const Floats& x,
426  DerivativeAccumulator &accum) const
427  {
428  double qval = x[0];
429  double value = get_value(qval) - A_val_;
430  double deriv;
431  //d[f(x)+A]/dG = f(x)/G
432  deriv = value/G_val_;
433  IMP_INTERNAL_CHECK(!base::isnan(deriv),
434  "derivative for G is nan.");
435  Scale(G_).add_to_nuisance_derivative(deriv, accum);
436  if (qval <= q1_param_)
437  {
438  //d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
439  deriv = - value * 2*IMP::square(qval)*Rg_val_/(3-s_val_);
440  IMP_INTERNAL_CHECK(!base::isnan(deriv),
441  "derivative for Rg is nan.");
442  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
443  //d[f(x)]/dd = 0
444  //
445  //d[f(x)]/ds = - f(x) * ( (q Rg / (3-s))^2 + log(q) )
446  deriv = - value
447  * (IMP::square((qval*Rg_val_)/(3-s_val_)) + std::log(qval));
448  IMP_INTERNAL_CHECK(!base::isnan(deriv),
449  "derivative for s is nan.");
450  Scale(s_).add_to_nuisance_derivative(deriv, accum);
451  } else {
452  //d[f(x)]/dRg = f(x) * (s-d)/Rg
453  deriv = value * (s_val_-d_val_)/Rg_val_;
454  IMP_INTERNAL_CHECK(!base::isnan(deriv),
455  "derivative for Rg is nan.");
456  Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
457  //d[f(x)]/dd = f(x) * log(q1/q)
458  deriv = value * std::log(q1_param_/qval);
459  IMP_INTERNAL_CHECK(!base::isnan(deriv),
460  "derivative for d is nan.");
461  Scale(d_).add_to_nuisance_derivative(deriv, accum);
462  //d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
463  deriv = - value * ( (d_val_-s_val_)/(2*(3-s_val_))
464  + std::log(q1_param_) );
465  IMP_INTERNAL_CHECK(!base::isnan(deriv),
466  "derivative for d is nan.");
467  Scale(s_).add_to_nuisance_derivative(deriv, accum);
468  }
469  //d[f(x)+A]/dA = 1
470  deriv = 1;
471  Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
472  }
473 
474  void add_to_particle_derivative(unsigned particle_no,
475  double value, DerivativeAccumulator &accum) const
476  {
477  switch (particle_no)
478  {
479  case 0:
480  IMP_INTERNAL_CHECK(!base::isnan(value),
481  "derivative for G is nan.");
482  Scale(G_).add_to_scale_derivative(value, accum);
483  break;
484  case 1:
485  IMP_INTERNAL_CHECK(!base::isnan(value),
486  "derivative for Rg is nan.");
487  Scale(Rg_).add_to_scale_derivative(value, accum);
488  break;
489  case 2:
490  IMP_INTERNAL_CHECK(!base::isnan(value),
491  "derivative for d is nan.");
492  Scale(d_).add_to_scale_derivative(value, accum);
493  break;
494  case 3:
495  IMP_INTERNAL_CHECK(!base::isnan(value),
496  "derivative for s is nan.");
497  Scale(s_).add_to_scale_derivative(value, accum);
498  break;
499  case 4:
500  IMP_INTERNAL_CHECK(!base::isnan(value),
501  "derivative for A is nan.");
502  Nuisance(A_).add_to_nuisance_derivative(value, accum);
503  break;
504  default:
505  IMP_THROW("Invalid particle number", ModelException);
506  }
507  }
508 
509  Eigen::VectorXd get_derivative_vector(
510  unsigned particle_no, const FloatsList& xlist) const
511  {
512  unsigned N=xlist.size();
513  Eigen::VectorXd ret(N);
514  switch(particle_no)
515  {
516  case 0: // G
517  //d[f(x)]/dG = f(x)/G
518  ret = ((*this)(xlist)
519  -Eigen::VectorXd::Constant(N,A_val_))/G_val_;
520  break;
521  case 1: // Rg
522  for (unsigned i=0; i<N; i++)
523  {
524  double qval = xlist[i][0];
525  if (qval <= q1_param_)
526  {
527  //d[f(x)]/dRg = - f(x) * 2 q^2 Rg / (3-s)
528  ret(i) = - (get_value(qval) - A_val_)
529  * 2*IMP::square(qval)*Rg_val_/(3-s_val_);
530  } else {
531  //d[f(x)]/dRg = f(x) * (s-d)/Rg
532  ret(i) = (get_value(qval)-A_val_)
533  * (s_val_-d_val_)/Rg_val_;
534  }
535  }
536  break;
537  case 2: // d
538  for (unsigned i=0; i<N; i++)
539  {
540  double qval = xlist[i][0];
541  if (qval <= q1_param_)
542  {
543  //d[f(x)]/dd = 0
544  ret(i) = 0;
545  } else {
546  //d[f(x)]/dd = f(x) * log(q1/q)
547  ret(i) = (get_value(qval)-A_val_)
548  * std::log(q1_param_/qval);
549  }
550  }
551  break;
552  case 3: // s
553  for (unsigned i=0; i<N; i++)
554  {
555  double qval = xlist[i][0];
556  if (qval <= q1_param_)
557  {
558  //d[f(x)]/ds = - f(x)
559  // * (q^2 Rg^2 + (3-s)^2 log(q)) / (3-s)^2
560  ret(i) = - (get_value(qval) - A_val_)
561  * (IMP::square((qval*Rg_val_)
562  /(3-s_val_))
563  + std::log(qval));
564  } else {
565  //d[f(x)]/ds = - f(x) * ( (d-s)/(2(3-s)) + log(q1) )
566  ret(i) = - (get_value(qval) - A_val_)
567  * ( (d_val_-s_val_)/(2*(3-s_val_))
568  + std::log(q1_param_) );
569  }
570  }
571  break;
572  case 4: // A
573  ret = Eigen::VectorXd::Constant(N,1);
574  break;
575  default:
576  IMP_THROW("Invalid particle number", ModelException);
577  }
578  return ret;
579  }
580 
581  FloatsList get_derivative_matrix(
582  const FloatsList& xlist, bool) const
583  {
584  Eigen::MatrixXd mat(xlist.size(),5);
585  mat.col(0) = get_derivative_vector(0, xlist);
586  mat.col(1) = get_derivative_vector(1, xlist);
587  mat.col(2) = get_derivative_vector(2, xlist);
588  mat.col(3) = get_derivative_vector(3, xlist);
589  mat.col(4) = get_derivative_vector(4, xlist);
590  FloatsList ret;
591  for (int i=0; i<mat.rows(); i++)
592  {
593  Floats line;
594  for (int j=0; j<mat.cols(); j++)
595  line.push_back(mat(i,j));
596  ret.push_back(line);
597  }
598  return ret;
599  }
600 
601  Eigen::VectorXd get_second_derivative_vector(
602  unsigned particle_a, unsigned particle_b,
603  const FloatsList& xlist) const
604  {
605  if (particle_a >=5)
606  IMP_THROW("Invalid particle 1 number", ModelException);
607  if (particle_b >= 5)
608  IMP_THROW("Invalid particle 2 number", ModelException);
609  unsigned N=xlist.size();
610  //hessian involving A is always 0
611  if (particle_a == 4 || particle_b == 4)
612  return Eigen::VectorXd::Zero(N);
613  unsigned pa = std::min(particle_a,particle_b);
614  unsigned pb = std::max(particle_a,particle_b);
615  Eigen::VectorXd ret(N);
616  switch(pa)
617  {
618  case 0: //G
619  switch(pb)
620  {
621  case 0: //G
622  //d2f/dG2 = 0
623  ret.noalias() = Eigen::VectorXd::Zero(N);
624  break;
625 
626  case 1: //Rg
627  //d2f/dGdRg = df/dRg * 1/G
628  ret.noalias() =
629  get_derivative_vector(1, xlist)/G_val_;
630  break;
631 
632  case 2: //d
633  //d2f/dGdd = df/dd * 1/G
634  ret.noalias() =
635  get_derivative_vector(2, xlist)/G_val_;
636  break;
637 
638  case 3: //s
639  //d2f/dGds = df/ds * 1/G
640  ret.noalias() =
641  get_derivative_vector(3, xlist)/G_val_;
642  break;
643 
644  default:
645  IMP_THROW("Invalid particle 2 number",
646  ModelException);
647  break;
648  }
649  break;
650 
651  case 1: // Rg
652  switch(pb)
653  {
654  case 1: //Rg
655  for (unsigned i=0; i<N; i++)
656  {
657  double qval = xlist[i][0];
658  if (qval <= q1_param_)
659  {
660  //d2[f(x)]/dRg2 =
661  //f(x) * (2 q^2 / (3-s))
662  // * (2 q^2 Rg^2 / (3-s) - 1)
663  ret(i) = (get_value(qval) - A_val_)
664  * 2*IMP::square(qval)/(3-s_val_)
665  * (2*IMP::square(qval*Rg_val_)/(3-s_val_) -1);
666  } else {
667  //d2[f(x)]/dRg2=f(x) * (d-s)/Rg^2 * (d-s+1)
668  ret(i) = (get_value(qval) - A_val_)
669  * (d_val_-s_val_)/IMP::square(Rg_val_)
670  * (d_val_-s_val_+1);
671  }
672  }
673  break;
674 
675  case 2: //d
676  for (unsigned i=0; i<N; i++)
677  {
678  double qval = xlist[i][0];
679  if (qval <= q1_param_)
680  {
681  //d2[f(x)]/dddRg = 0
682  ret(i) = 0;
683  } else {
684  //d2[f(x)]/dddRg = -f(x)/Rg
685  // - (d-s)/Rg *df(x)/dd
686  double val=(get_value(qval)-A_val_);
687  ret(i) = -val/Rg_val_
688  - (val*std::log(q1_param_/qval)
689  * (d_val_-s_val_)/Rg_val_);
690  }
691  }
692  break;
693 
694  case 3: //s
695  for (unsigned i=0; i<N; i++)
696  {
697  double qval = xlist[i][0];
698  double val=(get_value(qval)-A_val_);
699  if (qval <= q1_param_)
700  {
701  //d2[f(x)]/dsdRg = -2q^2Rg/(3-s)
702  // * (df(x)/ds + f(x)/(3-s))
703  double deriv = - val
704  * (IMP::square((qval*Rg_val_)
705  /(3-s_val_))
706  + std::log(qval));
707  ret(i) =
708  -2*IMP::square(qval)*Rg_val_/(3-s_val_)
709  *(deriv + val/(3-s_val_));
710  } else {
711  //d2[f(x)]/dsdRg =
712  // 1/Rg * (f(x)- (d-s)*df(x)/ds)
713  double deriv = - val
714  * ( (d_val_-s_val_)/(2*(3-s_val_))
715  + std::log(q1_param_) );
716  ret(i) = (val - (d_val_-s_val_)*deriv)
717  /Rg_val_;
718  }
719  }
720  break;
721 
722  default:
723  IMP_THROW("Invalid particle 2 number",
724  ModelException);
725  break;
726  }
727  break;
728 
729  case 2: //d
730  switch(pb)
731  {
732  case 2: //d
733  for (unsigned i=0; i<N; i++)
734  {
735  double qval = xlist[i][0];
736  if (qval <= q1_param_)
737  {
738  //d2[f(x)]/dddd = 0
739  ret(i) = 0;
740  } else {
741  //d2[f(x)]/dddd =
742  // f(x)*(log(q1/q)^2 + 1/(2*(d-s)))
743  double val=(get_value(qval)-A_val_);
744  ret(i) = val * (
745  IMP::square(std::log(q1_param_/qval))
746  +1/(2*(d_val_-s_val_)));
747  }
748  }
749  break;
750 
751  case 3: //s
752  {
753  double lterm=(d_val_-s_val_)/(2*(3-s_val_))
754  + std::log(q1_param_);
755  double rterm=0.5*(1/(3-s_val_) + 1/(d_val_-s_val_));
756  for (unsigned i=0; i<N; i++)
757  {
758  double qval = xlist[i][0];
759  if (qval <= q1_param_)
760  {
761  //d2[f(x)]/ddds = 0
762  ret(i) = 0;
763  } else {
764  //d2[f(x)]/ddds = log(q1/q)*df(x)/ds
765  // - (1/(3-s) + 1/(d-s))*f(x)/2
766  //
767  double val=(get_value(qval)-A_val_);
768  ret(i) = - val * (
769  std::log(q1_param_/qval)*lterm + rterm);
770  }
771  }
772  }
773  break;
774 
775  default:
776  IMP_THROW("Invalid particle 2 number",
777  ModelException);
778  break;
779  }
780  break;
781 
782  case 3: //s
783  switch(pb)
784  {
785  case 3: //s
786  {
787  double cterm =
788  IMP::square( 0.5*(d_val_-s_val_)/(3-s_val_)
789  + std::log(q1_param_) )
790  + 0.5*((6 - s_val_ - d_val_)/IMP::square(3-s_val_)
791  + 1./(d_val_-s_val_));
792  for (unsigned i=0; i<N; i++)
793  {
794  double qval = xlist[i][0];
795  double val=(get_value(qval)-A_val_);
796  if (qval <= q1_param_)
797  {
798  //d2[f(x)]/dsds = f(x) *
799  // ( [(qRg)^2/(3-s)^2 + log q ]^2
800  // - 2(qRg)^2/(3-s)^3 )
801  double factor =
802  IMP::square((qval*Rg_val_)/(3-s_val_));
803  ret(i) = val *
804  (IMP::square(factor + std::log(qval))
805  - 2*factor/(3-s_val_));
806  } else {
807  //d2[f(x)]/dsds = f(x)
808  // * ( [ (d-s)/(2*(3-s)) + log(q1) ]^2
809  // + 1/2 * [ (6-s-d)/(3-s)^2
810  // + 1/(d-s) ] )
811  ret(i) = val * cterm;
812  }
813  }
814  }
815  break;
816 
817  default:
818  IMP_THROW("Invalid particle 2 number",
819  ModelException);
820  break;
821  }
822  break;
823 
824  default:
825  IMP_THROW("Invalid particle 1 number",
826  ModelException);
827  break;
828  }
829  return ret;
830  }
831 
832  FloatsList get_second_derivative_vector(
833  unsigned particle_a, unsigned particle_b,
834  const FloatsList& xlist, bool) const
835  {
836  Eigen::VectorXd mat( get_second_derivative_vector(
837  particle_a, particle_b, xlist));
838  FloatsList ret;
839  for (int i=0; i<mat.rows(); i++)
840  {
841  Floats line;
842  for (int j=0; j<mat.cols(); j++)
843  line.push_back(mat(i,j));
844  ret.push_back(line);
845  }
846  return ret;
847  }
848 
849  unsigned get_ndims_x() const {return 1;}
850  unsigned get_ndims_y() const {return 1;}
851 
852  unsigned get_number_of_particles() const { return 5; }
853 
854  bool get_particle_is_optimized(unsigned particle_no) const
855  {
856  switch(particle_no)
857  {
858  case 0: //G
859  return Scale(G_).get_scale_is_optimized();
860  case 1: //Rg
861  return Scale(Rg_).get_scale_is_optimized();
862  case 2: //d
863  return Scale(d_).get_scale_is_optimized();
864  case 3: //s
865  return Scale(s_).get_scale_is_optimized();
866  case 4: //A
867  return Nuisance(A_).get_nuisance_is_optimized();
868  default:
869  IMP_THROW("Invalid particle number", ModelException);
870  }
871  }
872 
873  unsigned get_number_of_optimized_particles() const
874  {
875  unsigned count = 0;
876  if (Scale(G_).get_scale_is_optimized()) count++;
877  if (Scale(Rg_).get_scale_is_optimized()) count++;
878  if (Scale(d_).get_scale_is_optimized()) count++;
879  if (Scale(s_).get_scale_is_optimized()) count++;
880  if (Nuisance(A_).get_nuisance_is_optimized()) count++;
881  return count;
882  }
883 
884  ParticlesTemp get_input_particles() const
885  {
886  ParticlesTemp ret;
887  ret.push_back(G_);
888  ret.push_back(Rg_);
889  ret.push_back(d_);
890  ret.push_back(s_);
891  ret.push_back(A_);
892  return ret;
893  }
894 
895  ContainersTemp get_input_containers() const
896  {
897  ContainersTemp ret;
898  return ret;
899  }
900 
901  IMP_OBJECT_INLINE(GeneralizedGuinierPorodFunction, out
902  << " G = " << G_val_
903  << " Rg = " << Rg_val_
904  << " d = " << d_val_
905  << " s = " << s_val_
906  << " A = " << A_val_
907  << " Q1.Rg = " << q1_param_*Rg_val_
908  << std::endl, {});
909 
910  private:
911 
912  inline double get_value(double q) const
913  {
914  double value;
915  if (q <= q1_param_)
916  {
917  value = A_val_ + G_val_/std::pow(q,s_val_)
918  * std::exp(- IMP::square(q*Rg_val_)/(3-s_val_));
919  } else {
920  value = A_val_ + D_param_/std::pow(q,d_val_);
921  }
922  return value;
923  }
924 
925  Pointer<Particle> G_,Rg_,d_,s_,A_;
926  double G_val_,Rg_val_,d_val_,s_val_,A_val_,q1_param_,D_param_;
927 
928 };
929 
930 IMPISD_END_NAMESPACE
931 
932 #endif /* IMPISD_UNIVARIATE_FUNCTIONS_H */