IMP logo
IMP Reference Guide  2.18.0
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 
74  public:
75  //! Constructor (default)
76  RadialDistributionFunction(double bin_size = pr_resolution);
77 
78  //! Constructor from gnom file
79  RadialDistributionFunction(const std::string& file_name);
80 
81  //! scale distribution by a constant
82  void scale(double c);
83 
84  //! add another distribution
85  void add(const RadialDistributionFunction& model_pr);
86 
87  //! print tables
88  void show(std::ostream& out = std::cout) const;
89 
90  //! analogy crystallographic R-factor score
91  double R_factor_score(const RadialDistributionFunction& model_pr,
92  const std::string& file_name = "") const;
93 
94  // //! analogy to chi score \untested{chi_score}
95  // double chi_score(const RadialDistributionFunction& model_pr) const;
96 
97  //! fit the distributions by scaling according to maximum
98  double fit(const RadialDistributionFunction& model_pr,
99  const std::string& file_name = "") const;
100 
101  //! normalize to area = 1.0
102  void normalize();
103 
104  void add_to_distribution(double dist, double value) {
105  unsigned int index = get_index_from_distance(dist);
106  if (index >= size()) {
107  if (capacity() <= index)
108  reserve(2 * index); // to avoid many re-allocations
109  resize(index + 1, 0);
110  max_distance_ = get_distance_from_index(index + 1);
111  }
112  (*this)[index] += value;
113  }
114 
115  private:
116  //! read gnom file
117  void read_pr_file(const std::string& file_name);
118 
119  //! write fit file for the two distributions
120  void write_fit_file(const RadialDistributionFunction& model_pr, double c = 1.0,
121  const std::string& file_name = "") const;
122 };
123 
124 /**
125 Delta Distribution class for calculating the derivatives of SAXS Score
126 this distribution is:
127 sum_i [f_p(0) * f_i(0) * (x_p - x_i)]
128 sum_i [f_p(0) * f_i(0) * (y_p - y_i)]
129 sum_i [f_p(0) * f_i(0) * (z_p - z_i)]
130 */
131 class IMPSAXSEXPORT DeltaDistributionFunction
133  public:
134  //! Constructor
135  DeltaDistributionFunction(const Particles& particles,
136  double max_distance = 0.0,
137  double bin_size = pr_resolution);
138 
139  //! calculates distribution for an atom defined by particle
140  void calculate_derivative_distribution(Particle* particle);
141 
142  //! print tables
143  void show(std::ostream& out = std::cout, std::string prefix = "") const;
144 
145  private:
146  void add_to_distribution(double dist, const algebra::Vector3D& value) {
147  unsigned int index = get_index_from_distance(dist);
148  if (index >= size()) {
149  if (capacity() <= index)
150  reserve(2 * index); // to avoid many re-allocations
151  resize(index + 1, algebra::Vector3D(0.0, 0.0, 0.0));
152  max_distance_ = get_distance_from_index(index + 1);
153  }
154  (*this)[index] += value;
155  }
156 
157  void init() {
158  clear();
159  insert(begin(), get_index_from_distance(max_distance_) + 1,
160  algebra::Vector3D(0.0, 0.0, 0.0));
161  }
162 
163  protected:
164  std::vector<algebra::Vector3D> coordinates_;
165  Vector<double> form_factors_;
166 };
167 
168 IMPSAXS_END_NAMESPACE
169 
170 #endif /* IMPSAXS_DISTRIBUTION_H */
double get_max_distance() const
returns maximal distance value of distribution
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:421
Class to handle individual particles of a Model object.
Definition: Particle.h:41
A class for profile storing and computation.