8 #ifndef IMPISD_UNIVARIATE_FUNCTIONS_H
9 #define IMPISD_UNIVARIATE_FUNCTIONS_H
11 #include <IMP/isd/isd_config.h>
17 #include <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,
68 virtual Eigen::VectorXd get_derivative_vector(
69 unsigned particle_no,
const FloatsList& xlist)
const = 0;
73 bool stupid)
const = 0;
76 virtual Eigen::VectorXd get_second_derivative_vector(
77 unsigned particle_a,
unsigned particle_b,
81 virtual FloatsList get_second_derivative_vector(
unsigned particle_a,
84 bool stupid)
const = 0;
87 virtual unsigned get_ndims_x()
const = 0;
90 virtual unsigned get_ndims_y()
const = 0;
93 virtual unsigned get_number_of_particles()
const = 0;
96 virtual bool get_particle_is_optimized(
unsigned particle_no)
const = 0;
99 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 Eigen::VectorXd retlist(M);
149 for (
unsigned i = 0; i < M; i++) {
151 retlist(i) = a_val_ * xlist[i][0] + b_val_;
157 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);
187 unsigned N = xlist.size();
188 Eigen::VectorXd ret(N);
189 switch (particle_no) {
191 for (
unsigned i = 0; i < N; i++) ret(i) = xlist[i][0];
203 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 Eigen::VectorXd H(Eigen::VectorXd::Zero(N));
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();
255 if (
Nuisance(a_).get_nuisance_is_optimized()) count++;
256 if (
Nuisance(b_).get_nuisance_is_optimized()) count++;
273 double a_val_, b_val_;
298 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: constructor" << std::endl);
303 double tmpG =
Scale(G_).get_scale();
304 double tmpRg =
Scale(Rg_).get_scale();
305 double tmpd =
Scale(d_).get_scale();
306 double tmps =
Scale(s_).get_scale();
307 double tmpA =
Nuisance(A_).get_nuisance();
308 if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
309 (std::abs(tmpRg - Rg_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
310 (std::abs(tmpd - d_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
311 (std::abs(tmps - s_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
312 (std::abs(tmpA - A_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
313 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: has_changed():");
322 G_val_ =
Scale(G_).get_scale();
323 Rg_val_ =
Scale(Rg_).get_scale();
324 d_val_ =
Scale(d_).get_scale();
325 s_val_ =
Scale(s_).get_scale();
326 if (d_val_ == s_val_) {
328 if (s_val_ > 0.001) {
335 A_val_ =
Nuisance(A_).get_nuisance();
336 q1_param_ = std::sqrt((d_val_ - s_val_) * (3 - s_val_) / 2.);
337 D_param_ = G_val_ * std::exp(-IMP::square(q1_param_) / (3 - s_val_));
338 q1_param_ = q1_param_ / Rg_val_;
339 D_param_ *= std::pow(q1_param_, d_val_ - s_val_);
340 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: update() G:= "
341 << G_val_ <<
" Rg:=" << Rg_val_ <<
" d:=" << d_val_
342 <<
" s:=" << s_val_ <<
" A:=" << A_val_ <<
" Q1.Rg ="
343 << q1_param_ * Rg_val_ <<
" D =" << D_param_ << std::endl);
357 Floats ret(1, get_value(x[0]));
362 unsigned M = xlist.size();
363 Eigen::VectorXd retlist(M);
364 for (
unsigned i = 0; i < M; i++) {
366 retlist(i) = get_value(xlist[i][0]);
372 Eigen::VectorXd vec((*
this)(xlist));
374 for (
unsigned i = 0; i < xlist.size(); i++)
375 ret.push_back(
Floats(1, vec(i)));
381 double value = get_value(qval) - A_val_;
384 deriv = value / G_val_;
386 Scale(G_).add_to_nuisance_derivative(deriv, accum);
387 if (qval <= q1_param_) {
389 deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
391 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
396 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
398 Scale(s_).add_to_nuisance_derivative(deriv, accum);
401 deriv = value * (s_val_ - d_val_) / Rg_val_;
403 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
405 deriv = value * std::log(q1_param_ / qval);
407 Scale(d_).add_to_nuisance_derivative(deriv, accum);
410 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
412 Scale(s_).add_to_nuisance_derivative(deriv, accum);
416 Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
421 switch (particle_no) {
424 Scale(G_).add_to_scale_derivative(value, accum);
428 Scale(Rg_).add_to_scale_derivative(value, accum);
432 Scale(d_).add_to_scale_derivative(value, accum);
436 Scale(s_).add_to_scale_derivative(value, accum);
440 Nuisance(A_).add_to_nuisance_derivative(value, accum);
449 unsigned N = xlist.size();
450 Eigen::VectorXd ret(N);
451 switch (particle_no) {
454 ret = ((*this)(xlist) - Eigen::VectorXd::Constant(N, A_val_)) /
458 for (
unsigned i = 0; i < N; i++) {
459 double qval = xlist[i][0];
460 if (qval <= q1_param_) {
462 ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
463 Rg_val_ / (3 - s_val_);
466 ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
471 for (
unsigned i = 0; i < N; i++) {
472 double qval = xlist[i][0];
473 if (qval <= q1_param_) {
478 ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
483 for (
unsigned i = 0; i < N; i++) {
484 double qval = xlist[i][0];
485 if (qval <= q1_param_) {
489 -(get_value(qval) - A_val_) *
490 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
494 -(get_value(qval) - A_val_) *
495 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
500 ret = Eigen::VectorXd::Constant(N, 1);
509 Eigen::MatrixXd mat(xlist.size(), 5);
516 for (
int i = 0; i < mat.rows(); i++) {
518 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
525 unsigned particle_a,
unsigned particle_b,
const FloatsList& xlist)
const {
528 unsigned N = xlist.size();
530 if (particle_a == 4 || particle_b == 4)
return Eigen::VectorXd::Zero(N);
531 unsigned pa = std::min(particle_a, particle_b);
532 unsigned pb = std::max(particle_a, particle_b);
533 Eigen::VectorXd ret(N);
539 ret.noalias() = Eigen::VectorXd::Zero(N);
566 for (
unsigned i = 0; i < N; i++) {
567 double qval = xlist[i][0];
568 if (qval <= q1_param_) {
572 ret(i) = (get_value(qval) - A_val_) * 2 * IMP::square(qval) /
574 (2 * IMP::square(qval * Rg_val_) / (3 - s_val_) - 1);
577 ret(i) = (get_value(qval) - A_val_) * (d_val_ - s_val_) /
578 IMP::square(Rg_val_) * (d_val_ - s_val_ + 1);
584 for (
unsigned i = 0; i < N; i++) {
585 double qval = xlist[i][0];
586 if (qval <= q1_param_) {
592 double val = (get_value(qval) - A_val_);
593 ret(i) = -val / Rg_val_ - (val * std::log(q1_param_ / qval) *
594 (d_val_ - s_val_) / Rg_val_);
600 for (
unsigned i = 0; i < N; i++) {
601 double qval = xlist[i][0];
602 double val = (get_value(qval) - A_val_);
603 if (qval <= q1_param_) {
607 -val * (IMP::square((qval * Rg_val_) / (3 - s_val_)) +
609 ret(i) = -2 * IMP::square(qval) * Rg_val_ / (3 - s_val_) *
610 (deriv + val / (3 - s_val_));
614 double deriv = -val * ((d_val_ - s_val_) / (2 * (3 - s_val_)) +
615 std::log(q1_param_));
616 ret(i) = (val - (d_val_ - s_val_) * deriv) / Rg_val_;
630 for (
unsigned i = 0; i < N; i++) {
631 double qval = xlist[i][0];
632 if (qval <= q1_param_) {
638 double val = (get_value(qval) - A_val_);
639 ret(i) = val * (IMP::square(std::log(q1_param_ / qval)) +
640 1 / (2 * (d_val_ - s_val_)));
648 (d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_);
649 double rterm = 0.5 * (1 / (3 - s_val_) + 1 / (d_val_ - s_val_));
650 for (
unsigned i = 0; i < N; i++) {
651 double qval = xlist[i][0];
652 if (qval <= q1_param_) {
659 double val = (get_value(qval) - A_val_);
660 ret(i) = -val * (std::log(q1_param_ / qval) * lterm + rterm);
676 IMP::square(0.5 * (d_val_ - s_val_) / (3 - s_val_) +
677 std::log(q1_param_)) +
678 0.5 * ((6 - s_val_ - d_val_) / IMP::square(3 - s_val_) +
679 1. / (d_val_ - s_val_));
680 for (
unsigned i = 0; i < N; i++) {
681 double qval = xlist[i][0];
682 double val = (get_value(qval) - A_val_);
683 if (qval <= q1_param_) {
687 double factor = IMP::square((qval * Rg_val_) / (3 - s_val_));
688 ret(i) = val * (IMP::square(factor + std::log(qval)) -
689 2 * factor / (3 - s_val_));
695 ret(i) = val * cterm;
719 for (
int i = 0; i < mat.rows(); i++) {
721 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
733 switch (particle_no) {
735 return Scale(G_).get_scale_is_optimized();
737 return Scale(Rg_).get_scale_is_optimized();
739 return Scale(d_).get_scale_is_optimized();
741 return Scale(s_).get_scale_is_optimized();
743 return Nuisance(A_).get_nuisance_is_optimized();
751 if (
Scale(G_).get_scale_is_optimized()) count++;
752 if (
Scale(Rg_).get_scale_is_optimized()) count++;
753 if (
Scale(d_).get_scale_is_optimized()) count++;
754 if (
Scale(s_).get_scale_is_optimized()) count++;
755 if (
Nuisance(A_).get_nuisance_is_optimized()) count++;
780 inline double get_value(
double q)
const {
782 if (q <= q1_param_) {
783 value = A_val_ + G_val_ / std::pow(q, s_val_) *
784 std::exp(-IMP::square(q * Rg_val_) / (3 - s_val_));
786 value = A_val_ + D_param_ / std::pow(q, d_val_);
791 Pointer<Particle> G_, Rg_, d_, s_, A_;
792 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.
Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
Base class for functions of one variable.
virtual Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const =0
return second derivative vector
Floats operator()(const Floats &x) const
evaluate the function at a certain point
Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const
return second derivative vector
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
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.
Eigen::VectorXd get_second_derivative_vector(unsigned, unsigned, const FloatsList &xlist) const
return second derivative vector
Linear one-dimensional function.
virtual Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const =0
return derivative vector
ModelObjectsTemp get_inputs() const
particle manipulation
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const
for testing purposes
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.
Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const
return derivative vector
bool get_particle_is_optimized(unsigned particle_no) const
returns true if the particle whose index is provided is optimized
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 particles of a Model object.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Eigen::VectorXd operator()(const FloatsList &xlist) const
evaluate the function at a list of points
ModelObjectsTemp get_inputs() const
particle manipulation
unsigned get_ndims_y() const
returns the number of output dimensions
bool has_changed() const
return true if internal parameters have changed.
An exception which is thrown when the Model has attributes with invalid values.
unsigned get_ndims_x() const
returns the number of input dimensions
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const
for testing purposes
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