8 #ifndef IMPISD_UNIVARIATE_FUNCTIONS_H
9 #define IMPISD_UNIVARIATE_FUNCTIONS_H
11 #include <IMP/isd/isd_config.h>
17 #include <IMP/algebra/eigen3/Eigen/Dense>
19 #define IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM 1e-7
21 IMPISD_BEGIN_NAMESPACE
37 bool stupid)
const = 0;
40 virtual bool has_changed()
const = 0;
43 virtual void update() = 0;
49 virtual void add_to_derivatives(
const Floats& x,
57 virtual void add_to_particle_derivative(
unsigned particle_no,
double value,
67 virtual IMP_Eigen::VectorXd get_derivative_vector(
68 unsigned particle_no,
const FloatsList& xlist)
const = 0;
72 bool stupid)
const = 0;
75 virtual IMP_Eigen::VectorXd get_second_derivative_vector(
76 unsigned particle_a,
unsigned particle_b,
80 virtual FloatsList get_second_derivative_vector(
unsigned particle_a,
83 bool stupid)
const = 0;
86 virtual unsigned get_ndims_x()
const = 0;
89 virtual unsigned get_ndims_y()
const = 0;
92 virtual unsigned get_number_of_particles()
const = 0;
95 virtual bool get_particle_is_optimized(
unsigned particle_no)
const = 0;
98 virtual unsigned get_number_of_optimized_particles()
const = 0;
115 a_val_ =
Nuisance(a).get_nuisance();
116 b_val_ =
Nuisance(b).get_nuisance();
121 double tmpa =
Nuisance(a_).get_nuisance();
122 double tmpb =
Nuisance(b_).get_nuisance();
123 if ((std::abs(tmpa - a_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
124 (std::abs(tmpb - b_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
134 a_val_ =
Nuisance(a_).get_nuisance();
135 b_val_ =
Nuisance(b_).get_nuisance();
136 IMP_LOG_TERSE(
"Linear1DFunction: update() a:= " << a_val_ <<
" b:="
137 << b_val_ << std::endl);
142 Floats ret(1, a_val_ * x[0] + b_val_);
147 unsigned M = xlist.size();
148 IMP_Eigen::VectorXd retlist(M);
149 for (
unsigned i = 0; i < M; i++) {
151 retlist(i) = a_val_ * xlist[i][0] + b_val_;
157 IMP_Eigen::VectorXd vec((*
this)(xlist));
159 for (
unsigned i = 0; i < xlist.size(); i++)
160 ret.push_back(
Floats(1, vec(i)));
166 Nuisance(a_).add_to_nuisance_derivative(x[0], accum);
168 Nuisance(b_).add_to_nuisance_derivative(1, accum);
173 switch (particle_no) {
175 Nuisance(a_).add_to_nuisance_derivative(value, accum);
178 Nuisance(b_).add_to_nuisance_derivative(value, accum);
181 IMP_THROW(
"Invalid particle number", ModelException);
187 unsigned N = xlist.size();
188 IMP_Eigen::VectorXd ret(N);
189 switch (particle_no) {
191 for (
unsigned i = 0; i < N; i++) ret(i) = xlist[i][0];
197 IMP_THROW(
"Invalid particle number", ModelException);
203 IMP_Eigen::MatrixXd mat(xlist.size(), 2);
207 for (
int i = 0; i < mat.rows(); i++) {
209 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
216 unsigned,
unsigned,
const FloatsList& xlist)
const {
218 unsigned N = xlist.size();
219 IMP_Eigen::VectorXd H(IMP_Eigen::VectorXd::Zero(N));
226 IMP_Eigen::VectorXd mat(
229 for (
int i = 0; i < mat.rows(); i++) {
231 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
243 switch (particle_no) {
245 return Nuisance(a_).get_nuisance_is_optimized();
247 return Nuisance(b_).get_nuisance_is_optimized();
249 IMP_THROW(
"Invalid particle number", ModelException);
255 if (
Nuisance(a_).get_nuisance_is_optimized()) count++;
256 if (
Nuisance(b_).get_nuisance_is_optimized()) count++;
278 double a_val_, b_val_;
303 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: constructor" << std::endl);
308 double tmpG =
Scale(G_).get_scale();
309 double tmpRg =
Scale(Rg_).get_scale();
310 double tmpd =
Scale(d_).get_scale();
311 double tmps =
Scale(s_).get_scale();
312 double tmpA =
Nuisance(A_).get_nuisance();
313 if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
314 (std::abs(tmpRg - Rg_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
315 (std::abs(tmpd - d_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
316 (std::abs(tmps - s_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
317 (std::abs(tmpA - A_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
318 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: has_changed():");
327 G_val_ =
Scale(G_).get_scale();
328 Rg_val_ =
Scale(Rg_).get_scale();
329 d_val_ =
Scale(d_).get_scale();
330 s_val_ =
Scale(s_).get_scale();
331 if (d_val_ == s_val_) {
333 if (s_val_ > 0.001) {
340 A_val_ =
Nuisance(A_).get_nuisance();
341 q1_param_ = std::sqrt((d_val_ - s_val_) * (3 - s_val_) / 2.);
342 D_param_ = G_val_ * std::exp(-IMP::square(q1_param_) / (3 - s_val_));
343 q1_param_ = q1_param_ / Rg_val_;
344 D_param_ *= std::pow(q1_param_, d_val_ - s_val_);
345 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: update() G:= "
346 << G_val_ <<
" Rg:=" << Rg_val_ <<
" d:=" << d_val_
347 <<
" s:=" << s_val_ <<
" A:=" << A_val_ <<
" Q1.Rg ="
348 << q1_param_ * Rg_val_ <<
" D =" << D_param_ << std::endl);
362 Floats ret(1, get_value(x[0]));
367 unsigned M = xlist.size();
368 IMP_Eigen::VectorXd retlist(M);
369 for (
unsigned i = 0; i < M; i++) {
371 retlist(i) = get_value(xlist[i][0]);
377 IMP_Eigen::VectorXd vec((*
this)(xlist));
379 for (
unsigned i = 0; i < xlist.size(); i++)
380 ret.push_back(
Floats(1, vec(i)));
386 double value = get_value(qval) - A_val_;
389 deriv = value / G_val_;
391 Scale(G_).add_to_nuisance_derivative(deriv, accum);
392 if (qval <= q1_param_) {
394 deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
396 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
401 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
403 Scale(s_).add_to_nuisance_derivative(deriv, accum);
406 deriv = value * (s_val_ - d_val_) / Rg_val_;
408 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
410 deriv = value * std::log(q1_param_ / qval);
412 Scale(d_).add_to_nuisance_derivative(deriv, accum);
415 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
417 Scale(s_).add_to_nuisance_derivative(deriv, accum);
421 Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
426 switch (particle_no) {
429 Scale(G_).add_to_scale_derivative(value, accum);
433 Scale(Rg_).add_to_scale_derivative(value, accum);
437 Scale(d_).add_to_scale_derivative(value, accum);
441 Scale(s_).add_to_scale_derivative(value, accum);
445 Nuisance(A_).add_to_nuisance_derivative(value, accum);
448 IMP_THROW(
"Invalid particle number", ModelException);
454 unsigned N = xlist.size();
455 IMP_Eigen::VectorXd ret(N);
456 switch (particle_no) {
459 ret = ((*this)(xlist) - IMP_Eigen::VectorXd::Constant(N, A_val_)) /
463 for (
unsigned i = 0; i < N; i++) {
464 double qval = xlist[i][0];
465 if (qval <= q1_param_) {
467 ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
468 Rg_val_ / (3 - s_val_);
471 ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
476 for (
unsigned i = 0; i < N; i++) {
477 double qval = xlist[i][0];
478 if (qval <= q1_param_) {
483 ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
488 for (
unsigned i = 0; i < N; i++) {
489 double qval = xlist[i][0];
490 if (qval <= q1_param_) {
494 -(get_value(qval) - A_val_) *
495 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
499 -(get_value(qval) - A_val_) *
500 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
505 ret = IMP_Eigen::VectorXd::Constant(N, 1);
508 IMP_THROW(
"Invalid particle number", ModelException);
514 IMP_Eigen::MatrixXd mat(xlist.size(), 5);
521 for (
int i = 0; i < mat.rows(); i++) {
523 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
530 unsigned particle_a,
unsigned particle_b,
const FloatsList& xlist)
const {
531 if (particle_a >= 5)
IMP_THROW(
"Invalid particle 1 number", ModelException);
532 if (particle_b >= 5)
IMP_THROW(
"Invalid particle 2 number", ModelException);
533 unsigned N = xlist.size();
535 if (particle_a == 4 || particle_b == 4)
return IMP_Eigen::VectorXd::Zero(N);
536 unsigned pa = std::min(particle_a, particle_b);
537 unsigned pb = std::max(particle_a, particle_b);
538 IMP_Eigen::VectorXd ret(N);
544 ret.noalias() = IMP_Eigen::VectorXd::Zero(N);
563 IMP_THROW(
"Invalid particle 2 number", ModelException);
571 for (
unsigned i = 0; i < N; i++) {
572 double qval = xlist[i][0];
573 if (qval <= q1_param_) {
577 ret(i) = (get_value(qval) - A_val_) * 2 * IMP::square(qval) /
579 (2 * IMP::square(qval * Rg_val_) / (3 - s_val_) - 1);
582 ret(i) = (get_value(qval) - A_val_) * (d_val_ - s_val_) /
583 IMP::square(Rg_val_) * (d_val_ - s_val_ + 1);
589 for (
unsigned i = 0; i < N; i++) {
590 double qval = xlist[i][0];
591 if (qval <= q1_param_) {
597 double val = (get_value(qval) - A_val_);
598 ret(i) = -val / Rg_val_ - (val * std::log(q1_param_ / qval) *
599 (d_val_ - s_val_) / Rg_val_);
605 for (
unsigned i = 0; i < N; i++) {
606 double qval = xlist[i][0];
607 double val = (get_value(qval) - A_val_);
608 if (qval <= q1_param_) {
612 -val * (IMP::square((qval * Rg_val_) / (3 - s_val_)) +
614 ret(i) = -2 * IMP::square(qval) * Rg_val_ / (3 - s_val_) *
615 (deriv + val / (3 - s_val_));
619 double deriv = -val * ((d_val_ - s_val_) / (2 * (3 - s_val_)) +
620 std::log(q1_param_));
621 ret(i) = (val - (d_val_ - s_val_) * deriv) / Rg_val_;
627 IMP_THROW(
"Invalid particle 2 number", ModelException);
635 for (
unsigned i = 0; i < N; i++) {
636 double qval = xlist[i][0];
637 if (qval <= q1_param_) {
643 double val = (get_value(qval) - A_val_);
644 ret(i) = val * (IMP::square(std::log(q1_param_ / qval)) +
645 1 / (2 * (d_val_ - s_val_)));
653 (d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_);
654 double rterm = 0.5 * (1 / (3 - s_val_) + 1 / (d_val_ - s_val_));
655 for (
unsigned i = 0; i < N; i++) {
656 double qval = xlist[i][0];
657 if (qval <= q1_param_) {
664 double val = (get_value(qval) - A_val_);
665 ret(i) = -val * (std::log(q1_param_ / qval) * lterm + rterm);
671 IMP_THROW(
"Invalid particle 2 number", ModelException);
681 IMP::square(0.5 * (d_val_ - s_val_) / (3 - s_val_) +
682 std::log(q1_param_)) +
683 0.5 * ((6 - s_val_ - d_val_) / IMP::square(3 - s_val_) +
684 1. / (d_val_ - s_val_));
685 for (
unsigned i = 0; i < N; i++) {
686 double qval = xlist[i][0];
687 double val = (get_value(qval) - A_val_);
688 if (qval <= q1_param_) {
692 double factor = IMP::square((qval * Rg_val_) / (3 - s_val_));
693 ret(i) = val * (IMP::square(factor + std::log(qval)) -
694 2 * factor / (3 - s_val_));
700 ret(i) = val * cterm;
706 IMP_THROW(
"Invalid particle 2 number", ModelException);
712 IMP_THROW(
"Invalid particle 1 number", ModelException);
721 IMP_Eigen::VectorXd mat(
724 for (
int i = 0; i < mat.rows(); i++) {
726 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
738 switch (particle_no) {
740 return Scale(G_).get_scale_is_optimized();
742 return Scale(Rg_).get_scale_is_optimized();
744 return Scale(d_).get_scale_is_optimized();
746 return Scale(s_).get_scale_is_optimized();
748 return Nuisance(A_).get_nuisance_is_optimized();
750 IMP_THROW(
"Invalid particle number", ModelException);
756 if (
Scale(G_).get_scale_is_optimized()) count++;
757 if (
Scale(Rg_).get_scale_is_optimized()) count++;
758 if (
Scale(d_).get_scale_is_optimized()) count++;
759 if (
Scale(s_).get_scale_is_optimized()) count++;
760 if (
Nuisance(A_).get_nuisance_is_optimized()) count++;
790 inline double get_value(
double q)
const {
792 if (q <= q1_param_) {
793 value = A_val_ + G_val_ / std::pow(q, s_val_) *
794 std::exp(-IMP::square(q * Rg_val_) / (3 - s_val_));
796 value = A_val_ + D_param_ / std::pow(q, d_val_);
801 base::Pointer<kernel::Particle> G_, Rg_, d_, s_, A_;
802 double G_val_, Rg_val_, d_val_, s_val_, A_val_, q1_param_, D_param_;
bool has_changed() const
return true if internal parameters have changed.
Base class for functions of one variable.
Class for adding derivatives from restraints to the model.
Floats operator()(const Floats &x) const
evaluate the function at a certain point
FloatsList get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist, bool) const
for testing purposes
FloatsList operator()(const FloatsList &xlist, bool) const
used for testing only
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
unsigned get_ndims_y() const
returns the number of output dimensions
IMP_Eigen::VectorXd operator()(const FloatsList &xlist) const
evaluate the function at a list of points
unsigned get_ndims_x() const
returns the number of input dimensions
A decorator for switching parameters particles.
A decorator for scale parameters particles.
A smart pointer to a reference counted object.
ParticlesTemp get_input_particles(const ModelObjectsTemp &mos)
Return all the input particles for a given ModelObject.
A decorator for nuisance parameters particles.
void add_to_derivatives(const Floats &x, DerivativeAccumulator &accum) const
update derivatives of particles
Add scale parameter to particle.
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Ref counted objects should have private destructors.
bool get_particle_is_optimized(unsigned particle_no) const
returns true if the particle whose index is provided is optimized
1D mean function for SAS data
unsigned get_number_of_particles() const
returns the number of particles that this function uses
#define IMP_LOG_TERSE(expr)
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const
update derivatives of particles
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Linear one-dimensional function.
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const
for testing purposes
IMP_Eigen::VectorXd get_second_derivative_vector(unsigned, unsigned, const FloatsList &xlist) const
return second derivative vector
ContainersTemp get_input_containers(const ModelObjectsTemp &mos)
Return all the input particles for a given ModelObject.
unsigned get_number_of_particles() const
returns the number of particles that this function uses
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const
update derivatives of particles
Object(std::string name)
Construct an object with the given name.
FloatsList get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist, bool) const
for testing purposes
Add nuisance parameter to particle.
IMP::base::Vector< IMP::base::WeakPointer< Container > > ContainersTemp
bool get_particle_is_optimized(unsigned particle_no) const
returns true if the particle whose index is provided is optimized
virtual IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const =0
return derivative vector
void update()
update internal parameters
unsigned get_number_of_optimized_particles() const
returns the number of particles that are optimized
unsigned get_number_of_optimized_particles() const
returns the number of particles that are optimized
Class to handle individual model particles.
FloatsList operator()(const FloatsList &xlist, bool) const
used for testing only
Common base class for heavy weight IMP objects.
Classes to handle individual model particles. (Note that implementation of inline functions in in int...
Floats operator()(const Floats &x) const
evaluate the function at a certain point
kernel::ParticlesTemp get_input_particles() const
particle manipulation
void add_to_derivatives(const Floats &x, DerivativeAccumulator &accum) const
update derivatives of particles
virtual void update()=0
update internal parameters
bool isnan(const T &a)
Return true if a number is NaN.
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
A shared base class to help in debugging and things.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
IMP_Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const
return second derivative vector
unsigned get_ndims_y() const
returns the number of output dimensions
bool has_changed() const
return true if internal parameters have changed.
IMP_Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
unsigned get_ndims_x() const
returns the number of input dimensions
virtual IMP_Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const =0
return second derivative vector
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const
for testing purposes
IMP_Eigen::VectorXd operator()(const FloatsList &xlist) const
evaluate the function at a list of points
kernel::ParticlesTemp get_input_particles() const
particle manipulation
void update()
update internal parameters