IMP logo
IMP Reference Guide  develop.5651aa123e,2024/07/23
The Integrative Modeling Platform
saxs/distribution.h
Go to the documentation of this file.
1 /**
2  * \file IMP/saxs/Distribution.h \brief computes distribution functions
3  *
4  * Distribution - base distance distribution class
5  * RadialDistributionFunction required for calculation of SAXS profile
6  * DeltaDistributionFunction requires for chi-square derivatives
7  *
8  * Copyright 2007-2022 IMP Inventors. All rights reserved.
9  *
10  */
11 
12 #ifndef IMPSAXS_DISTRIBUTION_H
13 #define IMPSAXS_DISTRIBUTION_H
14 
15 #include <IMP/saxs/saxs_config.h>
16 #include "Profile.h"
17 #include <IMP/Particle.h>
18 
19 #include <iostream>
20 
21 IMPSAXS_BEGIN_NAMESPACE
22 
23 namespace { // anonymous
24 static const double pr_resolution = 0.5;
25 }
26 
27 /**
28 base class for distribution classes
29 */
30 template <class ValueT>
31 class Distribution : public std::vector<ValueT> {
32  public:
33  //! Constructor
34  Distribution(double bin_size = pr_resolution) { init(bin_size); }
35 
36  //! returns maximal distance value of distribution
37  double get_max_distance() const { return max_distance_; }
38 
39  //! returns bin size
40  double get_bin_size() const { return bin_size_; }
41 
42  unsigned int get_index_from_distance(double dist) const {
43  return algebra::get_rounded(dist * one_over_bin_size_);
44  }
45  double get_distance_from_index(unsigned int index) const {
46  return index * bin_size_;
47  }
48 
49  protected:
50  void init(double bin_size) {
51  // clear();
52  bin_size_ = bin_size;
53  one_over_bin_size_ = 1.0 / bin_size_; // for faster calculation
54  max_distance_ = 50.0; // start with ~50A (by default)
55  std::vector<ValueT>::reserve(get_index_from_distance(max_distance_) + 1);
56  }
57 
58  protected:
59  double bin_size_, one_over_bin_size_; // resolution of discretization
60  double max_distance_; // parameter for maximum r value for p(r) function
61 };
62 
63 #ifdef SWIG
64 %template(FloatDistribution) Distribution<double>;
65 %template(VectorDistribution) Distribution<algebra::Vector3D>;
66 #endif
67 
68 /**
69  Radial Distribution class for calculating SAXS Profile
70  this is distance distribution multiplied by form factors of atoms
71 */
72 class IMPSAXSEXPORT RadialDistributionFunction : public Distribution<double> {
73  mutable std::vector<double> sqrt_distances_;
74 
75  public:
76  //! Constructor (default)
77  RadialDistributionFunction(double bin_size = pr_resolution);
78 
79  //! Constructor from gnom file
80  RadialDistributionFunction(const std::string& file_name);
81 
82  //! Get square root of distance for each distribution point
83  /** This is cached, and points to internal storage, so should not be used
84  after the distribution is modified or destroyed. */
85  const std::vector<double> &get_square_root_distances() const {
86  size_t sz = size();
87  if (sz > sqrt_distances_.size()) {
88  sqrt_distances_.reserve(sz);
89  for (unsigned int r = sqrt_distances_.size(); r < sz; ++r) {
90  sqrt_distances_.push_back(sqrt(get_distance_from_index(r)));
91  }
92  }
93  return sqrt_distances_;
94  }
95 
96  //! scale distribution by a constant
97  void scale(double c);
98 
99  //! add another distribution
100  void add(const RadialDistributionFunction& model_pr);
101 
102  //! print tables
103  void show(std::ostream& out = std::cout) const;
104 
105  //! analogy crystallographic R-factor score
106  double R_factor_score(const RadialDistributionFunction& model_pr,
107  const std::string& file_name = "") const;
108 
109  // //! analogy to chi score \untested{chi_score}
110  // double chi_score(const RadialDistributionFunction& model_pr) const;
111 
112  //! fit the distributions by scaling according to maximum
113  double fit(const RadialDistributionFunction& model_pr,
114  const std::string& file_name = "") const;
115 
116  //! normalize to area = 1.0
117  void normalize();
118 
119  void add_to_distribution(double dist, double value) {
120  unsigned int index = get_index_from_distance(dist);
121  if (index >= size()) {
122  if (capacity() <= index)
123  reserve(2 * index); // to avoid many re-allocations
124  resize(index + 1, 0);
125  max_distance_ = get_distance_from_index(index + 1);
126  }
127  (*this)[index] += value;
128  }
129 
130  private:
131  //! read gnom file
132  void read_pr_file(const std::string& file_name);
133 
134  //! write fit file for the two distributions
135  void write_fit_file(const RadialDistributionFunction& model_pr, double c = 1.0,
136  const std::string& file_name = "") const;
137 };
138 
139 /**
140 Delta Distribution class for calculating the derivatives of SAXS Score
141 this distribution is:
142 sum_i [f_p(0) * f_i(0) * (x_p - x_i)]
143 sum_i [f_p(0) * f_i(0) * (y_p - y_i)]
144 sum_i [f_p(0) * f_i(0) * (z_p - z_i)]
145 */
146 class IMPSAXSEXPORT DeltaDistributionFunction
148  public:
149  //! Constructor
150  DeltaDistributionFunction(const Particles& particles,
151  double max_distance = 0.0,
152  double bin_size = pr_resolution);
153 
154  //! calculates distribution for an atom defined by particle
155  void calculate_derivative_distribution(Particle* particle);
156 
157  //! print tables
158  void show(std::ostream& out = std::cout, std::string prefix = "") const;
159 
160  private:
161  void add_to_distribution(double dist, const algebra::Vector3D& value) {
162  unsigned int index = get_index_from_distance(dist);
163  if (index >= size()) {
164  if (capacity() <= index)
165  reserve(2 * index); // to avoid many re-allocations
166  resize(index + 1, algebra::Vector3D(0.0, 0.0, 0.0));
167  max_distance_ = get_distance_from_index(index + 1);
168  }
169  (*this)[index] += value;
170  }
171 
172  void init() {
173  clear();
174  insert(begin(), get_index_from_distance(max_distance_) + 1,
175  algebra::Vector3D(0.0, 0.0, 0.0));
176  }
177 
178  protected:
179  std::vector<algebra::Vector3D> coordinates_;
180  Vector<double> form_factors_;
181 };
182 
183 IMPSAXS_END_NAMESPACE
184 
185 #endif /* IMPSAXS_DISTRIBUTION_H */
double get_max_distance() const
returns maximal distance value of distribution
const std::vector< double > & get_square_root_distances() const
Get square root of distance for each distribution point.
int get_rounded(const T &x)
Rounds a number to next integer.
double get_bin_size() const
returns bin size
Distribution(double bin_size=pr_resolution)
Constructor.
Classes to handle individual model particles. (Note that implementation of inline functions is in int...
void show(Hierarchy h, std::ostream &out=std::cout)
Print out a molecular hierarchy.
VectorD< 3 > Vector3D
Definition: VectorD.h:408
Class to handle individual particles of a Model object.
Definition: Particle.h:43
A class for profile storing and computation.