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)));
167 Nuisance(a_).add_to_nuisance_derivative(x[0], accum);
169 Nuisance(b_).add_to_nuisance_derivative(1, accum);
174 switch (particle_no) {
176 Nuisance(a_).add_to_nuisance_derivative(value, accum);
179 Nuisance(b_).add_to_nuisance_derivative(value, accum);
188 unsigned N = xlist.size();
189 Eigen::VectorXd ret(N);
190 switch (particle_no) {
192 for (
unsigned i = 0; i < N; i++) ret(i) = xlist[i][0];
204 const FloatsList& xlist,
bool)
const override {
205 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 override {
220 unsigned N = xlist.size();
221 Eigen::VectorXd H(Eigen::VectorXd::Zero(N));
226 unsigned particle_a,
unsigned particle_b,
227 const FloatsList& xlist,
bool)
const override {
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();
257 if (
Nuisance(a_).get_nuisance_is_optimized()) count++;
258 if (
Nuisance(b_).get_nuisance_is_optimized()) count++;
275 double a_val_, b_val_;
300 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: constructor" << std::endl);
305 double tmpG =
Scale(G_).get_scale();
306 double tmpRg =
Scale(Rg_).get_scale();
307 double tmpd =
Scale(d_).get_scale();
308 double tmps =
Scale(s_).get_scale();
309 double tmpA =
Nuisance(A_).get_nuisance();
310 if ((std::abs(tmpG - G_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
311 (std::abs(tmpRg - Rg_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
312 (std::abs(tmpd - d_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
313 (std::abs(tmps - s_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM) ||
314 (std::abs(tmpA - A_val_) > IMP_ISD_UNIVARIATE_FUNCTIONS_MINIMUM)) {
315 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: has_changed():");
324 G_val_ =
Scale(G_).get_scale();
325 Rg_val_ =
Scale(Rg_).get_scale();
326 d_val_ =
Scale(d_).get_scale();
327 s_val_ =
Scale(s_).get_scale();
328 if (d_val_ == s_val_) {
330 if (s_val_ > 0.001) {
337 A_val_ =
Nuisance(A_).get_nuisance();
338 q1_param_ = std::sqrt((d_val_ - s_val_) * (3 - s_val_) / 2.);
339 D_param_ = G_val_ * std::exp(-IMP::square(q1_param_) / (3 - s_val_));
340 q1_param_ = q1_param_ / Rg_val_;
341 D_param_ *= std::pow(q1_param_, d_val_ - s_val_);
342 IMP_LOG_TERSE(
"GeneralizedGuinierPorodFunction: update() G:= "
343 << G_val_ <<
" Rg:=" << Rg_val_ <<
" d:=" << d_val_
344 <<
" s:=" << s_val_ <<
" A:=" << A_val_ <<
" Q1.Rg ="
345 << q1_param_ * Rg_val_ <<
" D =" << D_param_ << std::endl);
359 Floats ret(1, get_value(x[0]));
364 unsigned M = xlist.size();
365 Eigen::VectorXd retlist(M);
366 for (
unsigned i = 0; i < M; i++) {
368 retlist(i) = get_value(xlist[i][0]);
374 Eigen::VectorXd vec((*
this)(xlist));
376 for (
unsigned i = 0; i < xlist.size(); i++)
377 ret.push_back(
Floats(1, vec(i)));
384 double value = get_value(qval) - A_val_;
387 deriv = value / G_val_;
389 Scale(G_).add_to_nuisance_derivative(deriv, accum);
390 if (qval <= q1_param_) {
392 deriv = -value * 2 * IMP::square(qval) * Rg_val_ / (3 - s_val_);
394 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
399 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
401 Scale(s_).add_to_nuisance_derivative(deriv, accum);
404 deriv = value * (s_val_ - d_val_) / Rg_val_;
406 Scale(Rg_).add_to_nuisance_derivative(deriv, accum);
408 deriv = value * std::log(q1_param_ / qval);
410 Scale(d_).add_to_nuisance_derivative(deriv, accum);
413 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
415 Scale(s_).add_to_nuisance_derivative(deriv, accum);
419 Nuisance(A_).add_to_nuisance_derivative(deriv, accum);
424 switch (particle_no) {
427 Scale(G_).add_to_scale_derivative(value, accum);
431 Scale(Rg_).add_to_scale_derivative(value, accum);
435 Scale(d_).add_to_scale_derivative(value, accum);
439 Scale(s_).add_to_scale_derivative(value, accum);
443 Nuisance(A_).add_to_nuisance_derivative(value, accum);
452 unsigned N = xlist.size();
453 Eigen::VectorXd ret(N);
454 switch (particle_no) {
457 ret = ((*this)(xlist) - Eigen::VectorXd::Constant(N, A_val_)) /
461 for (
unsigned i = 0; i < N; i++) {
462 double qval = xlist[i][0];
463 if (qval <= q1_param_) {
465 ret(i) = -(get_value(qval) - A_val_) * 2 * IMP::square(qval) *
466 Rg_val_ / (3 - s_val_);
469 ret(i) = (get_value(qval) - A_val_) * (s_val_ - d_val_) / Rg_val_;
474 for (
unsigned i = 0; i < N; i++) {
475 double qval = xlist[i][0];
476 if (qval <= q1_param_) {
481 ret(i) = (get_value(qval) - A_val_) * std::log(q1_param_ / qval);
486 for (
unsigned i = 0; i < N; i++) {
487 double qval = xlist[i][0];
488 if (qval <= q1_param_) {
492 -(get_value(qval) - A_val_) *
493 (IMP::square((qval * Rg_val_) / (3 - s_val_)) + std::log(qval));
497 -(get_value(qval) - A_val_) *
498 ((d_val_ - s_val_) / (2 * (3 - s_val_)) + std::log(q1_param_));
503 ret = Eigen::VectorXd::Constant(N, 1);
512 const FloatsList& xlist,
bool)
const override {
513 Eigen::MatrixXd mat(xlist.size(), 5);
520 for (
int i = 0; i < mat.rows(); i++) {
522 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
529 unsigned particle_a,
unsigned particle_b,
533 unsigned N = xlist.size();
535 if (particle_a == 4 || particle_b == 4)
return Eigen::VectorXd::Zero(N);
536 unsigned pa = std::min(particle_a, particle_b);
537 unsigned pb = std::max(particle_a, particle_b);
538 Eigen::VectorXd ret(N);
544 ret.noalias() = Eigen::VectorXd::Zero(N);
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_;
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);
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;
719 unsigned particle_a,
unsigned particle_b,
720 const FloatsList& xlist,
bool)
const override {
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();
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++;
785 inline double get_value(
double q)
const {
787 if (q <= q1_param_) {
788 value = A_val_ + G_val_ / std::pow(q, s_val_) *
789 std::exp(-IMP::square(q * Rg_val_) / (3 - s_val_));
791 value = A_val_ + D_param_ / std::pow(q, d_val_);
796 Pointer<Particle> G_, Rg_, d_, s_, A_;
797 double G_val_, Rg_val_, d_val_, s_val_, A_val_, q1_param_, D_param_;
Eigen::VectorXd get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const override
return second 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
unsigned get_ndims_x() const override
returns the number of input dimensions
bool has_changed() const override
return true if internal parameters have changed.
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const override
update derivatives of particles
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
void update() override
update internal parameters
A decorator for switching parameters particles.
A decorator for scale parameters particles.
Floats operator()(const Floats &x) const override
evaluate the function at a certain point
void add_to_derivatives(const Floats &x, DerivativeAccumulator &accum) const override
update derivatives of particles
A decorator for nuisance parameters particles.
void update() override
update internal parameters
Add scale parameter to particle.
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Set up destructor for a ref counted object.
1D mean function for SAS data
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const override
for testing purposes
unsigned get_number_of_optimized_particles() const override
returns the number of particles that are optimized
#define IMP_LOG_TERSE(expr)
Eigen::VectorXd operator()(const FloatsList &xlist) const override
evaluate the function at a list of points
A smart pointer to a reference counted object.
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
unsigned get_ndims_y() const override
returns the number of output dimensions
Linear one-dimensional function.
virtual Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const =0
return derivative vector
Common base class for heavy weight IMP objects.
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const override
update derivatives of particles
Add nuisance parameter to particle.
unsigned get_number_of_particles() const override
returns the number of particles that this function uses
unsigned get_number_of_particles() const override
returns the number of particles that this function uses
bool get_particle_is_optimized(unsigned particle_no) const override
returns true if the particle whose index is provided is optimized
Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const override
return derivative vector
FloatsList operator()(const FloatsList &xlist, bool) const override
used for testing only
bool isnan(const T &a)
Return true if a number is NaN.
Floats operator()(const Floats &x) const override
evaluate the function at a certain point
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
virtual void update()=0
update internal parameters
unsigned get_number_of_optimized_particles() const override
returns the number of particles that are optimized
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
A shared base class to help in debugging and things.
void add_to_derivatives(const Floats &x, DerivativeAccumulator &accum) const override
update derivatives of particles
Eigen::VectorXd get_second_derivative_vector(unsigned, unsigned, const FloatsList &xlist) const override
return second derivative vector
Object(std::string name)
Construct an object with the given name.
Class to handle individual particles of a Model object.
FloatsList get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist, bool) const override
for testing purposes
FloatsList get_derivative_matrix(const FloatsList &xlist, bool) const override
for testing purposes
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Eigen::VectorXd get_derivative_vector(unsigned particle_no, const FloatsList &xlist) const override
return derivative vector
Eigen::VectorXd operator()(const FloatsList &xlist) const override
evaluate the function at a list of points
ModelObjectsTemp get_inputs() const override
particle manipulation
bool get_particle_is_optimized(unsigned particle_no) const override
returns true if the particle whose index is provided is optimized
An exception which is thrown when the Model has attributes with invalid values.
FloatsList operator()(const FloatsList &xlist, bool) const override
used for testing only
bool has_changed() const override
return true if internal parameters have changed.
FloatsList get_second_derivative_vector(unsigned particle_a, unsigned particle_b, const FloatsList &xlist, bool) const override
for testing purposes
ModelObjectsTemp get_inputs() const override
particle manipulation
unsigned get_ndims_y() const override
returns the number of output dimensions
Class for adding derivatives from restraints to the model.
unsigned get_ndims_x() const override
returns the number of input dimensions