IMP logo
IMP Reference Guide  2.19.0
The Integrative Modeling Platform
Profile.h
Go to the documentation of this file.
1 /**
2  * \file IMP/saxs/Profile.h
3  * \brief A class for profile storing and computation
4  *
5  * Copyright 2007-2023 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPSAXS_PROFILE_H
10 #define IMPSAXS_PROFILE_H
11 
12 #include <IMP/saxs/saxs_config.h>
13 #include <IMP/Object.h>
14 #include <IMP/warning_macros.h>
15 
16 #include "FormFactorTable.h"
17 #include "Distribution.h"
18 
19 #include <iostream>
20 #include <vector>
21 #include <cereal/access.hpp>
22 #include <cereal/types/base_class.hpp>
23 #include <cereal/types/vector.hpp>
24 
25 IMPSAXS_BEGIN_NAMESPACE
26 
27 class RadialDistributionFunction;
28 
29 /**
30  Basic profile class, can be initialized from the input file
31  (experimental or theoretical) or computed from a set of Model
32  Particles (theoretical)
33 */
34 class IMPSAXSEXPORT Profile : public Object {
35  public:
36  // Constructors
37 
38  //! init from file
39  /**
40  \param[in] file_name profile file name
41  \param[in] fit_file if true, intensities are read from column 3
42  \param[in] max_q read till maximal q value = max_q, or all if max_q<=0
43  \param[in] units gets 1, 2, or 3 for unknown q units, 1/A, or 1/nm
44  */
45  Profile(const std::string& file_name, bool fit_file = false, double max_q = 0.0, int units = 1);
46 
47  //! init for theoretical profile
48  Profile(double qmin = 0.0, double qmax = 0.5, double delta = 0.005);
49 
50  // Various ways to compute a profile
51 
52  //! computes theoretical profile
53  void calculate_profile(const Particles& particles,
54  FormFactorType ff_type = HEAVY_ATOMS,
55  bool reciprocal = false) {
56  if (!reciprocal)
57  calculate_profile_real(particles, ff_type);
58  else
59  calculate_profile_reciprocal(particles, ff_type);
60  }
61 
62  //! compute profile for fitting with hydration layer and excluded volume
63  /**
64  A partial profile is a pre-computed profile, where intensities are
65  split into 6 parts that can be summed up into a regular profile
66  given a pair of c1/c2 values by sum_partial_profiles function.
67  see FoXS paper for details.
68  */
69  void calculate_profile_partial(const Particles& particles,
70  const Vector<double>& surface = Vector<double>(),
71  FormFactorType ff_type = HEAVY_ATOMS);
72 
73  //! compute profile for fitting with hydration layer and excluded volume
74  void calculate_profile_partial(const Particles& particles1,
75  const Particles& particles2,
76  const Vector<double>& surface1 = Vector<double>(),
77  const Vector<double>& surface2 = Vector<double>(),
78  FormFactorType ff_type = HEAVY_ATOMS);
79 
80  void calculate_profile_reciprocal_partial(const Particles& particles,
81  const Vector<double>& surface = Vector<double>(),
82  FormFactorType ff_type = HEAVY_ATOMS);
83 
84 
85  //! computes theoretical profile contribution from inter-molecular
86  //! interactions between the particles
87  void calculate_profile(const Particles& particles1,
88  const Particles& particles2,
89  FormFactorType ff_type = HEAVY_ATOMS) {
90  calculate_profile_real(particles1, particles2, ff_type);
91  }
92 
93  //! calculate Intensity at zero (= squared number of electrons)
94  double calculate_I0(const Particles& particles,
95  FormFactorType ff_type = HEAVY_ATOMS);
96 
97  //! calculate profile for any type of Particles that have coordinates
98  void calculate_profile_constant_form_factor(const Particles& particles,
99  double form_factor = 1.0);
100 
101 
102  // computes theoretical profile faster for cyclically symmetric particles
103  // assumes that the units particles are ordered one after another in the
104  // input particles vector (n - symmetry order)
105  void calculate_profile_symmetric(const Particles& particles,
106  unsigned int n,
107  FormFactorType ff_type = HEAVY_ATOMS);
108 
109  //! convert to real space P(r) function P(r) = 1/2PI^2 Sum(I(q)*qr*sin(qr))
110  void profile_2_distribution(RadialDistributionFunction& rd,
111  double max_distance) const;
112 
113  //! convert to reciprocal space I(q) = Sum(P(r)*sin(qr)/qr)
114  void distribution_2_profile(const RadialDistributionFunction& r_dist);
115 
116  //! return a profile that is sampled on the q values of the exp_profile
117  void resample(const Profile* exp_profile, Profile* resampled_profile) const;
118 
119  //! downsample the profile to a given number of points
120  void downsample(Profile* downsampled_profile,
121  unsigned int point_number) const;
122 
123  //! compute radius of gyration with Guinier approximation
124  /** ln[I(q)]=ln[I(0)] - (q^2*rg^2)/3
125  \param[in] end_q_rg determines the range of profile used for approximation:
126  i.e. q*rg < end_q_rg. Use 1.3 for globular proteins, 0.8 for elongated
127  */
128  double radius_of_gyration(double end_q_rg = 1.3) const;
129 
130  //! calculate mean intensity
131  double mean_intensity() const;
132 
133 
134  //! reads SAXS profile from file
135  /**
136  \param[in] file_name profile file name
137  \param[in] fit_file if true, intensities are read from column 3
138  \param[in] max_q read till maximal q value = max_q, or all if max_q<=0
139  \param[in] units gets 1, 2, or 3 for unknown q units, 1/A, or 1/nm
140  */
141  void read_SAXS_file(const std::string& file_name, bool fit_file = false, double max_q = 0.0, int units = 1);
142 
143  //! print to file
144  /** \param[in] file_name output file name
145  \param[in] max_q output till maximal q value = max_q, or all if max_q<=0
146  */
147  void write_SAXS_file(const std::string& file_name, double max_q = 0.0) const;
148 
149  //! read a partial profile from file (7 columns)
150  void read_partial_profiles(const std::string& file_name);
151 
152  //! write a partial profile to file
153  void write_partial_profiles(const std::string& file_name) const;
154 
155  // Access functions
156 
157  //! return sampling resolution
158  double get_delta_q() const { return delta_q_; }
159 
160  //! return minimal sampling point
161  double get_min_q() const { return min_q_; }
162 
163  //! return maximal sampling point
164  double get_max_q() const { return max_q_; }
165 
166  double get_intensity(unsigned int i) const { return intensity_(i); }
167  double get_q(unsigned int i) const { return q_(i); }
168  double get_error(unsigned int i) const { return error_(i); }
169  double get_weight(unsigned int i) const {
170  IMP_UNUSED(i);
171  return 1.0;
172  }
173 
174  const Eigen::VectorXf& get_qs() const { return q_; }
175  const Eigen::VectorXf& get_intensities() const { return intensity_; }
176  const Eigen::VectorXf& get_errors() const { return error_; }
177 
178  double get_average_radius() const { return average_radius_; }
179 
180  //! return number of entries in SAXS profile
181  unsigned int size() const { return q_.size(); }
182 
183  //! checks the sampling of experimental profile
184  bool is_uniform_sampling() const;
185 
186  bool is_partial_profile() const { return (partial_profiles_.size()>0); }
187 
188  std::string get_name() const { return name_; }
189 
190  unsigned int get_id() const { return id_; }
191 
192  // Modifiers
193  void set_qs(const Eigen::VectorXf& q) { q_ = q; }
194  void set_intensities(const Eigen::VectorXf& i) { intensity_ = i; }
195  void set_errors(const Eigen::VectorXf& e) { error_ = e; }
196 
197  void set_intensity(unsigned int i, double iq) { intensity_(i) = iq; }
198 
199  //! required for reciprocal space calculation
200  void set_ff_table(FormFactorTable* ff_table) { ff_table_ = ff_table; }
201 
202  void set_average_radius(double r) { average_radius_ = r; }
203 
204  void set_average_volume(double v) { average_volume_ = v; }
205 
206  void set_name(std::string name) { name_ = name; }
207 
208  void set_id(unsigned int id) { id_ = id; }
209 
210  void set_beam_profile(std::string beam_profile_file) {
211  beam_profile_ = new Profile(beam_profile_file);
212  }
213 
214  //! add simulated error
215  void add_errors();
216 
217  //! add simulated noise
218  void add_noise(double percentage = 0.03);
219 
220  //! computes full profile for given fitting parameters
221  void sum_partial_profiles(double c1, double c2, bool check_cashed = true);
222 
223  //! add another profile - useful for rigid bodies
224  void add(const Profile* other_profile, double weight = 1.0);
225 
226  //! add partial profiles
227  void add_partial_profiles(const Profile* other_profile, double weight = 1.0);
228 
229  //! add other profiles - useful for weighted ensembles
230  void add(const Vector<Profile*>& profiles,
231  const Vector<double>& weights = Vector<double>());
232 
233  //! add other partial profiles
234  void add_partial_profiles(const Vector<Profile*>& profiles,
235  const Vector<double>& weights = Vector<double>());
236 
237  //! background adjustment option
238  void background_adjust(double start_q);
239 
240  //! scale
241  void scale(double c);
242 
243  //! offset profile by c, I(q) = I(q) - c
244  void offset(double c);
245 
246  // copy error bars from the matching experimental profile
247  void copy_errors(const Profile* exp_profile);
248 
249  // parameter for E^2(q), used in faster calculation
250  static const double modulation_function_parameter_;
251 
253 
254  protected:
255  void init(unsigned int size = 0, unsigned int partial_profiles_size = 0);
256 
257  private:
258  void calculate_profile_reciprocal(const Particles& particles,
259  FormFactorType ff_type = HEAVY_ATOMS);
260 
261  void calculate_profile_reciprocal(const Particles& particles1,
262  const Particles& particles2,
263  FormFactorType ff_type = HEAVY_ATOMS);
264 
265  void calculate_profile_real(const Particles& particles,
266  FormFactorType ff_type = HEAVY_ATOMS);
267 
268  void calculate_profile_real(const Particles& particles1,
269  const Particles& particles2,
270  FormFactorType ff_type = HEAVY_ATOMS);
271 
272  void squared_distribution_2_profile(const RadialDistributionFunction& r_dist);
273 
274  void squared_distributions_2_partial_profiles(
275  const Vector<RadialDistributionFunction>& r_dist);
276 
277  double radius_of_gyration_fixed_q(double end_q) const;
278 
279  double find_max_q(const std::string& file_name) const;
280 
281  protected:
282  Eigen::VectorXf q_; // q sampling points
283  Eigen::VectorXf intensity_;
284  Eigen::VectorXf error_; // error bar of each point
285 
286  double min_q_, max_q_; // minimal and maximal q values in the profile
287  double delta_q_; // profile sampling resolution
288  FormFactorTable* ff_table_; // pointer to form factors table
289 
290  // stores the intensity split into 6 for c1/c2 enumeration
291  std::vector<Eigen::VectorXf> partial_profiles_;
292  double c1_, c2_;
293 
294  bool experimental_; // experimental profile read from file
295  double average_radius_; // average radius of the particles
296  double average_volume_; // average volume
297 
298  // mapping from q values to vector index for fast profile resampling
299  mutable std::map<double, unsigned int> q_mapping_;
300 
301  std::string name_; // file name
302  unsigned int id_; // identifier
303 
304  Pointer<Profile> beam_profile_;
305  private:
306  friend class cereal::access;
307  template<class Archive> void serialize(Archive &ar) {
308  ar(cereal::base_class<Object>(this), q_, intensity_, error_,
309  min_q_, max_q_, delta_q_, partial_profiles_, c1_, c2_,
310  experimental_, average_radius_, average_volume_,
311  name_, id_, beam_profile_);
312  if (std::is_base_of<cereal::detail::InputArchiveBase, Archive>::value) {
313  // q_mapping_ is regenerated when needed in resample()
314  q_mapping_.clear();
315  bool default_ff_table;
316  ar(default_ff_table);
317  if (default_ff_table) {
318  ff_table_ = get_default_form_factor_table();
319  } else {
320  ff_table_ = nullptr;
321  }
322  } else {
323  if (ff_table_ == nullptr) {
324  ar(false);
325  } else if (ff_table_ == get_default_form_factor_table()) {
326  ar(true);
327  } else {
328  IMP_THROW("Serialization of profiles using non-default form "
329  "factors is not supported", IMP::ValueException);
330  }
331  }
332  }
333  IMP_OBJECT_SERIALIZE_DECL(Profile);
334 };
335 
337 
338 IMPSAXS_END_NAMESPACE
339 
340 namespace cereal {
341  template<class Archive, typename _Scalar, int _Rows, int _Cols,
342  int _Options, int _MaxRows, int _MaxCols>
343  inline void serialize(
344  Archive &ar, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows,
345  _MaxCols> &matrix) {
346  int rows, cols;
347  if (std::is_base_of<cereal::detail::OutputArchiveBase, Archive>::value) {
348  rows = matrix.rows();
349  cols = matrix.cols();
350  }
351  ar(rows, cols);
352 
353  if (std::is_base_of<cereal::detail::InputArchiveBase, Archive>::value) {
354  if (rows != matrix.rows() || cols != matrix.cols()) {
355  matrix.resize(rows, cols);
356  }
357  }
358  auto mat_data = cereal::binary_data(matrix.data(),
359  rows * cols * sizeof(_Scalar));
360  if (matrix.size() != 0) {
361  ar(mat_data);
362  }
363  }
364 }
365 
366 #endif /* IMPSAXS_PROFILE_H */
unsigned int size() const
return number of entries in SAXS profile
Definition: Profile.h:181
void calculate_profile(const Particles &particles, FormFactorType ff_type=HEAVY_ATOMS, bool reciprocal=false)
computes theoretical profile
Definition: Profile.h:53
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
FormFactorTable * get_default_form_factor_table()
A smart pointer to a reference counted object.
Definition: Pointer.h:87
Common base class for heavy weight IMP objects.
Definition: Object.h:111
Macros to control compiler warnings.
#define IMP_UNUSED(variable)
A class for computation of atomic and residue level form factors for SAXS calculations.
#define IMP_OBJECT_SERIALIZE_DECL(Name)
Declare methods needed for serialization of Object pointers.
Definition: object_macros.h:95
double get_delta_q() const
return sampling resolution
Definition: Profile.h:158
#define IMP_OBJECTS(Name, PluralName)
Define the types for storing lists of object pointers.
Definition: object_macros.h:44
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
Definition: check_macros.h:50
double get_min_q() const
return minimal sampling point
Definition: Profile.h:161
A shared base class to help in debugging and things.
FormFactorType
type of the form factors for profile calculations
void set_ff_table(FormFactorTable *ff_table)
required for reciprocal space calculation
Definition: Profile.h:200
void calculate_profile(const Particles &particles1, const Particles &particles2, FormFactorType ff_type=HEAVY_ATOMS)
Definition: Profile.h:87
computes distribution functions
double get_max_q() const
return maximal sampling point
Definition: Profile.h:164
An exception for an invalid value being passed to IMP.
Definition: exception.h:136
double radius_of_gyration(const Particles &particles)
compute radius_of_gyration
Definition: saxs/utility.h:74