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;
114 a_val_ =
Nuisance(a).get_nuisance();
115 b_val_ =
Nuisance(b).get_nuisance();
120 double tmpa =
Nuisance(a_).get_nuisance();
121 double tmpb =
Nuisance(b_).get_nuisance();
122 if ((std::abs(tmpa - a_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
123 (std::abs(tmpb - b_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
133 a_val_ =
Nuisance(a_).get_nuisance();
134 b_val_ =
Nuisance(b_).get_nuisance();
135 IMP_LOG_TERSE(
"Linear1DFunction: update() a:= " << a_val_ <<
" b:="
136 << b_val_ << std::endl);
141 Floats ret(1, a_val_ * x[0] + b_val_);
146 unsigned M = xlist.size();
147 IMP_Eigen::VectorXd retlist(M);
148 for (
unsigned i = 0; i < M; i++) {
150 retlist(i) = a_val_ * xlist[i][0] + b_val_;
156 IMP_Eigen::VectorXd vec((*
this)(xlist));
158 for (
unsigned i = 0; i < xlist.size(); i++)
159 ret.push_back(
Floats(1, vec(i)));
165 Nuisance(a_).add_to_nuisance_derivative(x[0], accum);
167 Nuisance(b_).add_to_nuisance_derivative(1, accum);
172 switch (particle_no) {
174 Nuisance(a_).add_to_nuisance_derivative(value, accum);
177 Nuisance(b_).add_to_nuisance_derivative(value, accum);
186 unsigned N = xlist.size();
187 IMP_Eigen::VectorXd ret(N);
188 switch (particle_no) {
190 for (
unsigned i = 0; i < N; i++) ret(i) = xlist[i][0];
202 IMP_Eigen::MatrixXd mat(xlist.size(), 2);
206 for (
int i = 0; i < mat.rows(); i++) {
208 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
215 unsigned,
unsigned,
const FloatsList& xlist)
const {
217 unsigned N = xlist.size();
218 IMP_Eigen::VectorXd H(IMP_Eigen::VectorXd::Zero(N));
225 IMP_Eigen::VectorXd mat(
228 for (
int i = 0; i < mat.rows(); i++) {
230 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
242 switch (particle_no) {
244 return Nuisance(a_).get_nuisance_is_optimized();
246 return Nuisance(b_).get_nuisance_is_optimized();
254 if (
Nuisance(a_).get_nuisance_is_optimized()) count++;
255 if (
Nuisance(b_).get_nuisance_is_optimized()) count++;
272 double a_val_, b_val_;
297 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: constructor" << std::endl);
302 double tmpG =
Scale(G_).get_scale();
303 double tmpRg =
Scale(Rg_).get_scale();
304 double tmpd =
Scale(d_).get_scale();
305 double tmps =
Scale(s_).get_scale();
306 double tmpA =
Nuisance(A_).get_nuisance();
307 if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
308 (std::abs(tmpRg - Rg_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
309 (std::abs(tmpd - d_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
310 (std::abs(tmps - s_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
311 (std::abs(tmpA - A_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
312 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: has_changed():");
321 G_val_ =
Scale(G_).get_scale();
322 Rg_val_ =
Scale(Rg_).get_scale();
323 d_val_ =
Scale(d_).get_scale();
324 s_val_ =
Scale(s_).get_scale();
325 if (d_val_ == s_val_) {
327 if (s_val_ > 0.001) {
334 A_val_ =
Nuisance(A_).get_nuisance();
335 q1_param_ = std::sqrt((d_val_ - s_val_) * (3 - s_val_) / 2.);
336 D_param_ = G_val_ * std::exp(-IMP::square(q1_param_) / (3 - s_val_));
337 q1_param_ = q1_param_ / Rg_val_;
338 D_param_ *= std::pow(q1_param_, d_val_ - s_val_);
339 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: update() G:= "
340 << G_val_ <<
" Rg:=" << Rg_val_ <<
" d:=" << d_val_
341 <<
" s:=" << s_val_ <<
" A:=" << A_val_ <<
" Q1.Rg ="
342 << q1_param_ * Rg_val_ <<
" D =" << D_param_ << std::endl);
356 Floats ret(1, get_value(x[0]));
361 unsigned M = xlist.size();
362 IMP_Eigen::VectorXd retlist(M);
363 for (
unsigned i = 0; i < M; i++) {
365 retlist(i) = get_value(xlist[i][0]);
371 IMP_Eigen::VectorXd vec((*
this)(xlist));
373 for (
unsigned i = 0; i < xlist.size(); i++)
374 ret.push_back(
Floats(1, vec(i)));
380 double value = get_value(qval) - A_val_;
383 deriv = value / G_val_;
385 Scale(G_).add_to_nuisance_derivative(deriv, accum);
386 if (qval <= q1_param_) {
388 deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
390 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
395 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
397 Scale(s_).add_to_nuisance_derivative(deriv, accum);
400 deriv = value * (s_val_ - d_val_) / Rg_val_;
402 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
404 deriv = value * std::log(q1_param_ / qval);
406 Scale(d_).add_to_nuisance_derivative(deriv, accum);
409 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
411 Scale(s_).add_to_nuisance_derivative(deriv, accum);
415 Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
420 switch (particle_no) {
423 Scale(G_).add_to_scale_derivative(value, accum);
427 Scale(Rg_).add_to_scale_derivative(value, accum);
431 Scale(d_).add_to_scale_derivative(value, accum);
435 Scale(s_).add_to_scale_derivative(value, accum);
439 Nuisance(A_).add_to_nuisance_derivative(value, accum);
448 unsigned N = xlist.size();
449 IMP_Eigen::VectorXd ret(N);
450 switch (particle_no) {
453 ret = ((*this)(xlist) - IMP_Eigen::VectorXd::Constant(N, A_val_)) /
457 for (
unsigned i = 0; i < N; i++) {
458 double qval = xlist[i][0];
459 if (qval <= q1_param_) {
461 ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
462 Rg_val_ / (3 - s_val_);
465 ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
470 for (
unsigned i = 0; i < N; i++) {
471 double qval = xlist[i][0];
472 if (qval <= q1_param_) {
477 ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
482 for (
unsigned i = 0; i < N; i++) {
483 double qval = xlist[i][0];
484 if (qval <= q1_param_) {
488 -(get_value(qval) - A_val_) *
489 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
493 -(get_value(qval) - A_val_) *
494 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
499 ret = IMP_Eigen::VectorXd::Constant(N, 1);
508 IMP_Eigen::MatrixXd mat(xlist.size(), 5);
515 for (
int i = 0; i < mat.rows(); i++) {
517 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
524 unsigned particle_a,
unsigned particle_b,
const FloatsList& xlist)
const {
527 unsigned N = xlist.size();
529 if (particle_a == 4 || particle_b == 4)
return IMP_Eigen::VectorXd::Zero(N);
530 unsigned pa = std::min(particle_a, particle_b);
531 unsigned pb = std::max(particle_a, particle_b);
532 IMP_Eigen::VectorXd ret(N);
538 ret.noalias() = IMP_Eigen::VectorXd::Zero(N);
565 for (
unsigned i = 0; i < N; i++) {
566 double qval = xlist[i][0];
567 if (qval <= q1_param_) {
571 ret(i) = (get_value(qval) - A_val_) * 2 * IMP::square(qval) /
573 (2 * IMP::square(qval * Rg_val_) / (3 - s_val_) - 1);
576 ret(i) = (get_value(qval) - A_val_) * (d_val_ - s_val_) /
577 IMP::square(Rg_val_) * (d_val_ - s_val_ + 1);
583 for (
unsigned i = 0; i < N; i++) {
584 double qval = xlist[i][0];
585 if (qval <= q1_param_) {
591 double val = (get_value(qval) - A_val_);
592 ret(i) = -val / Rg_val_ - (val * std::log(q1_param_ / qval) *
593 (d_val_ - s_val_) / Rg_val_);
599 for (
unsigned i = 0; i < N; i++) {
600 double qval = xlist[i][0];
601 double val = (get_value(qval) - A_val_);
602 if (qval <= q1_param_) {
606 -val * (IMP::square((qval * Rg_val_) / (3 - s_val_)) +
608 ret(i) = -2 * IMP::square(qval) * Rg_val_ / (3 - s_val_) *
609 (deriv + val / (3 - s_val_));
613 double deriv = -val * ((d_val_ - s_val_) / (2 * (3 - s_val_)) +
614 std::log(q1_param_));
615 ret(i) = (val - (d_val_ - s_val_) * deriv) / Rg_val_;
629 for (
unsigned i = 0; i < N; i++) {
630 double qval = xlist[i][0];
631 if (qval <= q1_param_) {
637 double val = (get_value(qval) - A_val_);
638 ret(i) = val * (IMP::square(std::log(q1_param_ / qval)) +
639 1 / (2 * (d_val_ - s_val_)));
647 (d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_);
648 double rterm = 0.5 * (1 / (3 - s_val_) + 1 / (d_val_ - s_val_));
649 for (
unsigned i = 0; i < N; i++) {
650 double qval = xlist[i][0];
651 if (qval <= q1_param_) {
658 double val = (get_value(qval) - A_val_);
659 ret(i) = -val * (std::log(q1_param_ / qval) * lterm + rterm);
675 IMP::square(0.5 * (d_val_ - s_val_) / (3 - s_val_) +
676 std::log(q1_param_)) +
677 0.5 * ((6 - s_val_ - d_val_) / IMP::square(3 - s_val_) +
678 1. / (d_val_ - s_val_));
679 for (
unsigned i = 0; i < N; i++) {
680 double qval = xlist[i][0];
681 double val = (get_value(qval) - A_val_);
682 if (qval <= q1_param_) {
686 double factor = IMP::square((qval * Rg_val_) / (3 - s_val_));
687 ret(i) = val * (IMP::square(factor + std::log(qval)) -
688 2 * factor / (3 - s_val_));
694 ret(i) = val * cterm;
715 IMP_Eigen::VectorXd mat(
718 for (
int i = 0; i < mat.rows(); i++) {
720 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
732 switch (particle_no) {
734 return Scale(G_).get_scale_is_optimized();
736 return Scale(Rg_).get_scale_is_optimized();
738 return Scale(d_).get_scale_is_optimized();
740 return Scale(s_).get_scale_is_optimized();
742 return Nuisance(A_).get_nuisance_is_optimized();
750 if (
Scale(G_).get_scale_is_optimized()) count++;
751 if (
Scale(Rg_).get_scale_is_optimized()) count++;
752 if (
Scale(d_).get_scale_is_optimized()) count++;
753 if (
Scale(s_).get_scale_is_optimized()) count++;
754 if (
Nuisance(A_).get_nuisance_is_optimized()) count++;
779 inline double get_value(
double q)
const {
781 if (q <= q1_param_) {
782 value = A_val_ + G_val_ / std::pow(q, s_val_) *
783 std::exp(-IMP::square(q * Rg_val_) / (3 - s_val_));
785 value = A_val_ + D_param_ / std::pow(q, d_val_);
790 Pointer<Particle> G_, Rg_, d_, s_, A_;
791 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.
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 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)
A smart pointer to a reference counted object.
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.
ModelObjectsTemp get_inputs() const
particle manipulation
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
Common base class for heavy weight IMP objects.
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
FloatsList get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist, bool) const
for testing purposes
Add nuisance parameter to particle.
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
FloatsList operator()(const FloatsList &xlist, bool) const
used for testing only
bool isnan(const T &a)
Return true if a number is NaN.
Floats operator()(const Floats &x) const
evaluate the function at a certain point
void add_to_derivatives(const Floats &x, DerivativeAccumulator &accum) const
update derivatives of particles
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
virtual void update()=0
update internal parameters
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
A shared base class to help in debugging and things.
Object(std::string name)
Construct an object with the given name.
Class to handle individual model particles.
#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
ModelObjectsTemp get_inputs() const
particle manipulation
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
An exception which is thrown when the Model has attributes with invalid values.
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
Class for adding derivatives from restraints to the model.
void update()
update internal parameters