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;
117 a_val_ =
Nuisance(a).get_nuisance();
118 b_val_ =
Nuisance(b).get_nuisance();
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)) {
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);
144 Floats ret(1, a_val_ * x[0] + b_val_);
149 unsigned M = xlist.size();
150 IMP_Eigen::VectorXd retlist(M);
151 for (
unsigned i = 0; i < M; i++) {
153 retlist(i) = a_val_ * xlist[i][0] + b_val_;
159 IMP_Eigen::VectorXd vec((*
this)(xlist));
161 for (
unsigned i = 0; i < xlist.size(); i++)
162 ret.push_back(
Floats(1, vec(i)));
168 Nuisance(a_).add_to_nuisance_derivative(x[0], accum);
170 Nuisance(b_).add_to_nuisance_derivative(1, accum);
175 switch (particle_no) {
177 Nuisance(a_).add_to_nuisance_derivative(value, accum);
180 Nuisance(b_).add_to_nuisance_derivative(value, accum);
183 IMP_THROW(
"Invalid particle number", ModelException);
189 unsigned N = xlist.size();
190 IMP_Eigen::VectorXd ret(N);
191 switch (particle_no) {
193 for (
unsigned i = 0; i < N; i++) ret(i) = xlist[i][0];
199 IMP_THROW(
"Invalid particle number", ModelException);
205 IMP_Eigen::MatrixXd mat(xlist.size(), 2);
209 for (
int i = 0; i < mat.rows(); i++) {
211 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
218 unsigned,
unsigned,
const FloatsList& xlist)
const {
220 unsigned N = xlist.size();
221 IMP_Eigen::VectorXd H(IMP_Eigen::VectorXd::Zero(N));
228 IMP_Eigen::VectorXd mat(
231 for (
int i = 0; i < mat.rows(); i++) {
233 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
245 switch (particle_no) {
247 return Nuisance(a_).get_nuisance_is_optimized();
249 return Nuisance(b_).get_nuisance_is_optimized();
251 IMP_THROW(
"Invalid particle number", ModelException);
257 if (
Nuisance(a_).get_nuisance_is_optimized()) count++;
258 if (
Nuisance(b_).get_nuisance_is_optimized()) count++;
280 double a_val_, b_val_;
305 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: constructor" << std::endl);
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():");
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_) {
340 if (s_val_ > 0.001) {
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);
369 Floats ret(1, get_value(x[0]));
374 unsigned M = xlist.size();
375 IMP_Eigen::VectorXd retlist(M);
376 for (
unsigned i = 0; i < M; i++) {
378 retlist(i) = get_value(xlist[i][0]);
384 IMP_Eigen::VectorXd vec((*
this)(xlist));
386 for (
unsigned i = 0; i < xlist.size(); i++)
387 ret.push_back(
Floats(1, vec(i)));
393 double value = get_value(qval) - A_val_;
396 deriv = value / G_val_;
398 Scale(G_).add_to_nuisance_derivative(deriv, accum);
399 if (qval <= q1_param_) {
401 deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
403 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
408 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
410 Scale(s_).add_to_nuisance_derivative(deriv, accum);
413 deriv = value * (s_val_ - d_val_) / Rg_val_;
415 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
417 deriv = value * std::log(q1_param_ / qval);
419 Scale(d_).add_to_nuisance_derivative(deriv, accum);
422 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
424 Scale(s_).add_to_nuisance_derivative(deriv, accum);
428 Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
433 switch (particle_no) {
436 Scale(G_).add_to_scale_derivative(value, accum);
440 Scale(Rg_).add_to_scale_derivative(value, accum);
444 Scale(d_).add_to_scale_derivative(value, accum);
448 Scale(s_).add_to_scale_derivative(value, accum);
452 Nuisance(A_).add_to_nuisance_derivative(value, accum);
455 IMP_THROW(
"Invalid particle number", ModelException);
461 unsigned N = xlist.size();
462 IMP_Eigen::VectorXd ret(N);
463 switch (particle_no) {
466 ret = ((*this)(xlist) - IMP_Eigen::VectorXd::Constant(N, A_val_)) /
470 for (
unsigned i = 0; i < N; i++) {
471 double qval = xlist[i][0];
472 if (qval <= q1_param_) {
474 ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
475 Rg_val_ / (3 - s_val_);
478 ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
483 for (
unsigned i = 0; i < N; i++) {
484 double qval = xlist[i][0];
485 if (qval <= q1_param_) {
490 ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
495 for (
unsigned i = 0; i < N; i++) {
496 double qval = xlist[i][0];
497 if (qval <= q1_param_) {
501 -(get_value(qval) - A_val_) *
502 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
506 -(get_value(qval) - A_val_) *
507 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
512 ret = IMP_Eigen::VectorXd::Constant(N, 1);
515 IMP_THROW(
"Invalid particle number", ModelException);
521 IMP_Eigen::MatrixXd mat(xlist.size(), 5);
528 for (
int i = 0; i < mat.rows(); i++) {
530 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
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();
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);
551 ret.noalias() = IMP_Eigen::VectorXd::Zero(N);
570 IMP_THROW(
"Invalid particle 2 number", ModelException);
578 for (
unsigned i = 0; i < N; i++) {
579 double qval = xlist[i][0];
580 if (qval <= q1_param_) {
584 ret(i) = (get_value(qval) - A_val_) * 2 * IMP::square(qval) /
586 (2 * IMP::square(qval * Rg_val_) / (3 - s_val_) - 1);
589 ret(i) = (get_value(qval) - A_val_) * (d_val_ - s_val_) /
590 IMP::square(Rg_val_) * (d_val_ - s_val_ + 1);
596 for (
unsigned i = 0; i < N; i++) {
597 double qval = xlist[i][0];
598 if (qval <= q1_param_) {
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_);
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_) {
619 -val * (IMP::square((qval * Rg_val_) / (3 - s_val_)) +
621 ret(i) = -2 * IMP::square(qval) * Rg_val_ / (3 - s_val_) *
622 (deriv + val / (3 - s_val_));
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_;
634 IMP_THROW(
"Invalid particle 2 number", ModelException);
642 for (
unsigned i = 0; i < N; i++) {
643 double qval = xlist[i][0];
644 if (qval <= q1_param_) {
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_)));
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_) {
671 double val = (get_value(qval) - A_val_);
672 ret(i) = -val * (std::log(q1_param_ / qval) * lterm + rterm);
678 IMP_THROW(
"Invalid particle 2 number", ModelException);
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_) {
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_));
707 ret(i) = val * cterm;
713 IMP_THROW(
"Invalid particle 2 number", ModelException);
719 IMP_THROW(
"Invalid particle 1 number", ModelException);
728 IMP_Eigen::VectorXd mat(
731 for (
int i = 0; i < mat.rows(); i++) {
733 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
745 switch (particle_no) {
747 return Scale(G_).get_scale_is_optimized();
749 return Scale(Rg_).get_scale_is_optimized();
751 return Scale(d_).get_scale_is_optimized();
753 return Scale(s_).get_scale_is_optimized();
755 return Nuisance(A_).get_nuisance_is_optimized();
757 IMP_THROW(
"Invalid particle number", ModelException);
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++;
797 inline double get_value(
double q)
const {
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_));
803 value = A_val_ + D_param_ / std::pow(q, d_val_);
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_;
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.
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.
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
static Scale decorate_particle(::IMP::kernel::Particle *p)
#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.
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)
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.
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.
#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