8 #ifndef IMPISD_BIVARIATE_FUNCTIONS_H
9 #define IMPISD_BIVARIATE_FUNCTIONS_H
11 #include <IMP/isd/isd_config.h>
17 #include <Eigen/Dense>
19 #define IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM 1e-7
21 IMPISD_BEGIN_NAMESPACE
40 bool stupid)
const = 0;
43 virtual bool has_changed()
const = 0;
46 virtual void update() = 0;
49 virtual void add_to_derivatives(
const Floats& x1,
const Floats& x2,
57 virtual void add_to_particle_derivative(
unsigned particle_no,
double value,
65 virtual Eigen::MatrixXd get_derivative_matrix(
66 unsigned particle_no,
const FloatsList& xlist)
const = 0;
69 virtual FloatsList get_derivative_matrix(
unsigned particle_no,
71 bool stupid)
const = 0;
74 virtual Eigen::MatrixXd get_second_derivative_matrix(
75 unsigned particle_a,
unsigned particle_b,
79 virtual FloatsList get_second_derivative_matrix(
unsigned particle_a,
82 bool stupid)
const = 0;
85 virtual unsigned get_ndims_x1()
const = 0;
86 virtual unsigned get_ndims_x2()
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;
121 double alpha = 2.0,
double jitter = 0.0,
122 double cutoff = 1e-7)
129 IMP_LOG_TERSE(
"Covariance1DFunction: constructor" << std::endl);
130 lambda_val_ =
Scale(ilambda).get_nuisance();
131 tau_val_ =
Scale(tau).get_nuisance();
132 do_jitter = (jitter > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
133 alpha_square_ = (std::abs(alpha - 2) < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM);
138 double tmpt =
Scale(tau_).get_nuisance();
139 double tmpl =
Scale(lambda_).get_nuisance();
143 if ((std::abs(tmpt - tau_val_) > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) ||
144 (std::abs(tmpl - lambda_val_) > IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM)) {
154 lambda_val_ =
Scale(lambda_).get_nuisance();
155 tau_val_ =
Scale(tau_).get_nuisance();
157 << tau_val_ <<
" lambda:=" << lambda_val_ << std::endl);
168 Floats ret(1, get_value(x1[0], x2[0]));
173 const unsigned M = xlist.size();
174 Eigen::MatrixXd Mret(M, M);
175 for (
unsigned i = 0; i < M; i++) {
176 for (
unsigned j = i; j < M; j++) {
179 double x1 = xlist[i][0];
180 double x2 = xlist[j][0];
181 double ret = get_value(x1, x2);
183 if (i != j) Mret(j, i) = ret;
190 Eigen::MatrixXd mat((*
this)(xlist));
192 for (
unsigned i = 0; i < xlist.size(); i++)
193 for (
unsigned j = 0; j < xlist.size(); j++)
194 ret.push_back(
Floats(1, mat(i, j)));
201 double val = get_value(x1[0], x2[0]);
202 double tauderiv = 2. / tau_val_ * val;
204 Scale(tau_).add_to_nuisance_derivative(tauderiv, accum);
209 (alpha_ * std::pow((std::abs(x1[0] - x2[0]) / lambda_val_), alpha_) /
212 Scale(lambda_).add_to_nuisance_derivative(lambdaderiv, accum);
217 switch (particle_no) {
220 Scale(tau_).add_to_nuisance_derivative(value, accum);
224 Scale(lambda_).add_to_nuisance_derivative(value, accum);
236 unsigned N = xlist.size();
237 Eigen::MatrixXd ret(N, N);
239 switch (particle_no) {
243 diag = get_value(xlist[0][0], xlist[0][0]);
244 diag *= 2. / tau_val_;
256 "derivative matrix is nan on the diagonal.");
257 for (
unsigned i = 0; i < N; i++) ret(i, i) = diag;
259 bool initial_loop =
true;
260 double abs_cutoff = cutoff_ * diag;
262 for (
unsigned i = 0; i < N; i++) {
263 for (
unsigned j = i + 1; j < N; j++) {
264 double x1(xlist[i][0]), x2(xlist[j][0]);
266 double dist(std::abs(x1 - x2));
270 if (initial_loop || dist <= dmax) {
271 switch (particle_no) {
275 val = get_value(xlist[i][0], xlist[j][0]);
276 val *= 2. / tau_val_;
282 if (dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
285 val = get_value(xlist[i][0], xlist[j][0]);
286 val *= alpha_ * std::pow((dist / lambda_val_), alpha_) /
296 if (std::abs(val) <= abs_cutoff) {
298 initial_loop =
false;
300 }
else if (dist < dmax) {
308 "derivative matrix is nan at position("
309 << i <<
"," << j <<
").");
318 const FloatsList& xlist,
bool)
const override {
321 for (
int i = 0; i < mat.rows(); i++) {
323 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
330 unsigned particle_a,
unsigned particle_b,
332 unsigned N(xlist.size());
333 Eigen::MatrixXd ret(N, N);
337 if (particle_a == 0 && particle_b == 0) {
338 for (
unsigned i = 0; i < N; i++) {
339 for (
unsigned j = i; j < N; j++) {
340 double dist = std::abs(xlist[i][0] - xlist[j][0]);
341 double exponent = std::pow(dist / lambda_val_, alpha_);
342 double expterm = std::exp(-0.5 * exponent);
343 ret(i, j) = 2 * expterm;
344 if (i != j) ret(j, i) = ret(i, j);
347 }
else if (particle_a == 1 && particle_b == 1) {
348 for (
unsigned i = 0; i < N; i++) {
349 for (
unsigned j = i; j < N; j++) {
350 double dist = std::abs(xlist[i][0] - xlist[j][0]);
351 double exponent = std::pow(dist / lambda_val_, alpha_);
352 double expterm = std::exp(-0.5 * exponent);
353 ret(i, j) = tau_val_ * tau_val_ * expterm * exponent /
354 (lambda_val_ * lambda_val_) * alpha_ / 2 *
355 (alpha_ / 2 * exponent - (alpha_ + 1));
356 if (i != j) ret(j, i) = ret(i, j);
360 for (
unsigned i = 0; i < N; i++) {
361 for (
unsigned j = i; j < N; j++) {
362 double dist = std::abs(xlist[i][0] - xlist[j][0]);
363 double exponent = std::pow(dist / lambda_val_, alpha_);
364 double expterm = std::exp(-0.5 * exponent);
365 ret(i, j) = tau_val_ * alpha_ * expterm / lambda_val_ * exponent;
366 if (i != j) ret(j, i) = ret(i, j);
374 unsigned particle_a,
unsigned particle_b,
375 const FloatsList& xlist,
bool)
const override {
379 for (
int i = 0; i < mat.rows(); i++) {
381 for (
int j = 0; j < mat.cols(); j++) line.push_back(mat(i, j));
388 unsigned get_ndims_x2()
const override {
return 1; }
394 switch (particle_no) {
396 return Scale(tau_).get_nuisance_is_optimized();
398 return Scale(lambda_).get_nuisance_is_optimized();
406 if (
Scale(tau_).get_nuisance_is_optimized()) count++;
407 if (
Scale(lambda_).get_nuisance_is_optimized()) count++;
414 ret.push_back(lambda_);
424 inline double get_value(
double x1,
double x2)
const {
425 double dist = std::abs(x1 - x2);
426 double ret = dist / lambda_val_;
430 ret = std::pow(ret, alpha_);
432 ret = IMP::square(tau_val_) * std::exp(-0.5 * ret);
433 if (do_jitter && dist < IMP_ISD_BIVARIATE_FUNCTIONS_MINIMUM) {
434 ret += IMP::square(tau_val_) * J_;
437 "function value is nan. tau = "
438 << tau_val_ <<
" lambda = " << lambda_val_
439 <<
" q1 = " << x1 <<
" q2 = " << x2);
445 Pointer<Particle> tau_, lambda_;
446 double tau_val_, lambda_val_, J_, cutoff_, alpha_square_;
Eigen::MatrixXd operator()(const IMP::FloatsList &xlist) const override
evaluate the function at a list of points
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
A decorator for switching parameters particles.
unsigned get_ndims_x1() const override
returns the number of input dimensions
void add_to_particle_derivative(unsigned particle_no, double value, DerivativeAccumulator &accum) const override
update derivatives of particles
A decorator for scale parameters particles.
A decorator for nuisance parameters particles.
Add scale parameter to particle.
#define IMP_LOG_VERBOSE(expr)
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Set up destructor for a ref counted object.
#define IMP_LOG_TERSE(expr)
Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const override
return derivative matrix
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
void update() override
update internal parameters
FloatsList operator()(const IMP::FloatsList &xlist, bool) const override
used for testing only
Common base class for heavy weight IMP objects.
virtual Eigen::MatrixXd get_second_derivative_matrix(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const =0
return second derivative matrix
virtual void update()=0
update internal parameters
bool has_changed() const override
return true if internal parameters have changed.
Floats operator()(const Floats &x1, const Floats &x2) const override
evaluate the function at a certain point
unsigned get_number_of_particles() const override
returns the number of particles that this function uses
ModelObjectsTemp get_inputs() const override
particle manipulation
Base class for functions of two variables.
bool isnan(const T &a)
Return true if a number is NaN.
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
Eigen::MatrixXd get_second_derivative_matrix(unsigned particle_a, unsigned particle_b, const FloatsList &xlist) const override
return second derivative matrix
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.
unsigned get_ndims_y() const override
returns the number of output dimensions
void add_to_derivatives(const Floats &x1, const Floats &x2, DerivativeAccumulator &accum) const override
update derivatives of particles
An exception which is thrown when the Model has attributes with invalid values.
Class for adding derivatives from restraints to the model.
unsigned get_number_of_optimized_particles() const override
returns the number of particles that are optimized
bool get_particle_is_optimized(unsigned particle_no) const override
returns true if the particle whose index is provided is optimized
virtual Eigen::MatrixXd get_derivative_matrix(unsigned particle_no, const FloatsList &xlist) const =0
return derivative matrix