IMP logo
IMP Reference Guide  2.14.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-2020 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 
22 IMPSAXS_BEGIN_NAMESPACE
23 
24 class RadialDistributionFunction;
25 
26 /**
27  Basic profile class, can be initialized from the input file
28  (experimental or theoretical) or computed from a set of Model
29  Particles (theoretical)
30 */
31 class IMPSAXSEXPORT Profile : public Object {
32  public:
33  // Constructors
34 
35  //! init from file
36  /**
37  \param[in] file_name profile file name
38  \param[in] fit_file if true, intensities are read from column 3
39  \param[in] max_q read till maximal q value = max_q, or all if max_q<=0
40  \param[in] units gets 1, 2, or 3 for unknown q units, 1/A, or 1/nm
41  */
42  Profile(const std::string& file_name, bool fit_file = false, double max_q = 0.0, int units = 1);
43 
44  //! init for theoretical profile
45  Profile(double qmin = 0.0, double qmax = 0.5, double delta = 0.005);
46 
47  // Various ways to compute a profile
48 
49  //! computes theoretical profile
50  void calculate_profile(const Particles& particles,
51  FormFactorType ff_type = HEAVY_ATOMS,
52  bool reciprocal = false) {
53  if (!reciprocal)
54  calculate_profile_real(particles, ff_type);
55  else
56  calculate_profile_reciprocal(particles, ff_type);
57  }
58 
59  //! compute profile for fitting with hydration layer and excluded volume
60  /**
61  A partial profile is a pre-computed profile, where intensities are
62  split into 6 parts that can be summed up into a regular profile
63  given a pair of c1/c2 values by sum_partial_profiles function.
64  see FoXS paper for details.
65  */
66  void calculate_profile_partial(const Particles& particles,
67  const Vector<double>& surface = Vector<double>(),
68  FormFactorType ff_type = HEAVY_ATOMS);
69 
70  //! compute profile for fitting with hydration layer and excluded volume
71  void calculate_profile_partial(const Particles& particles1,
72  const Particles& particles2,
73  const Vector<double>& surface1 = Vector<double>(),
74  const Vector<double>& surface2 = Vector<double>(),
75  FormFactorType ff_type = HEAVY_ATOMS);
76 
77  void calculate_profile_reciprocal_partial(const Particles& particles,
78  const Vector<double>& surface = Vector<double>(),
79  FormFactorType ff_type = HEAVY_ATOMS);
80 
81 
82  //! computes theoretical profile contribution from inter-molecular
83  //! interactions between the particles
84  void calculate_profile(const Particles& particles1,
85  const Particles& particles2,
86  FormFactorType ff_type = HEAVY_ATOMS) {
87  calculate_profile_real(particles1, particles2, ff_type);
88  }
89 
90  //! calculate Intensity at zero (= squared number of electrons)
91  double calculate_I0(const Particles& particles,
92  FormFactorType ff_type = HEAVY_ATOMS);
93 
94  //! calculate profile for any type of Particles that have coordinates
95  void calculate_profile_constant_form_factor(const Particles& particles,
96  double form_factor = 1.0);
97 
98 
99  // computes theoretical profile faster for cyclically symmetric particles
100  // assumes that the units particles are ordered one after another in the
101  // input particles vector (n - symmetry order)
102  void calculate_profile_symmetric(const Particles& particles,
103  unsigned int n,
104  FormFactorType ff_type = HEAVY_ATOMS);
105 
106  //! convert to real space P(r) function P(r) = 1/2PI^2 Sum(I(q)*qr*sin(qr))
107  void profile_2_distribution(RadialDistributionFunction& rd,
108  double max_distance) const;
109 
110  //! convert to reciprocal space I(q) = Sum(P(r)*sin(qr)/qr)
111  void distribution_2_profile(const RadialDistributionFunction& r_dist);
112 
113  //! return a profile that is sampled on the q values of the exp_profile
114  void resample(const Profile* exp_profile, Profile* resampled_profile) const;
115 
116  //! downsample the profile to a given number of points
117  void downsample(Profile* downsampled_profile,
118  unsigned int point_number) const;
119 
120  //! compute radius of gyration with Guinier approximation
121  /** ln[I(q)]=ln[I(0)] - (q^2*rg^2)/3
122  \param[in] end_q_rg determines the range of profile used for approximation:
123  i.e. q*rg < end_q_rg. Use 1.3 for globular proteins, 0.8 for elongated
124  */
125  double radius_of_gyration(double end_q_rg = 1.3) const;
126 
127  //! calculate mean intensity
128  double mean_intensity() const;
129 
130 
131  //! reads SAXS profile from file
132  /**
133  \param[in] file_name profile file name
134  \param[in] fit_file if true, intensities are read from column 3
135  \param[in] max_q read till maximal q value = max_q, or all if max_q<=0
136  \param[in] units gets 1, 2, or 3 for unknown q units, 1/A, or 1/nm
137  */
138  void read_SAXS_file(const std::string& file_name, bool fit_file = false, double max_q = 0.0, int units = 1);
139 
140  //! print to file
141  /** \param[in] file_name output file name
142  \param[in] max_q output till maximal q value = max_q, or all if max_q<=0
143  */
144  void write_SAXS_file(const std::string& file_name, double max_q = 0.0) const;
145 
146  //! read a partial profile from file (7 columns)
147  void read_partial_profiles(const std::string& file_name);
148 
149  //! write a partial profile to file
150  void write_partial_profiles(const std::string& file_name) const;
151 
152  // Access functions
153 
154  //! return sampling resolution
155  double get_delta_q() const { return delta_q_; }
156 
157  //! return minimal sampling point
158  double get_min_q() const { return min_q_; }
159 
160  //! return maximal sampling point
161  double get_max_q() const { return max_q_; }
162 
163  double get_intensity(unsigned int i) const { return intensity_(i); }
164  double get_q(unsigned int i) const { return q_(i); }
165  double get_error(unsigned int i) const { return error_(i); }
166  double get_weight(unsigned int i) const {
167  IMP_UNUSED(i);
168  return 1.0;
169  }
170 
171  const Eigen::VectorXf& get_qs() const { return q_; }
172  const Eigen::VectorXf& get_intensities() const { return intensity_; }
173  const Eigen::VectorXf& get_errors() const { return error_; }
174 
175  double get_average_radius() const { return average_radius_; }
176 
177  //! return number of entries in SAXS profile
178  unsigned int size() const { return q_.size(); }
179 
180  //! checks the sampling of experimental profile
181  bool is_uniform_sampling() const;
182 
183  bool is_partial_profile() const { return (partial_profiles_.size()>0); }
184 
185  std::string get_name() const { return name_; }
186 
187  unsigned int get_id() const { return id_; }
188 
189  // Modifiers
190  void set_qs(const Eigen::VectorXf& q) { q_ = q; }
191  void set_intensities(const Eigen::VectorXf& i) { intensity_ = i; }
192  void set_errors(const Eigen::VectorXf& e) { error_ = e; }
193 
194  void set_intensity(unsigned int i, double iq) { intensity_(i) = iq; }
195 
196  //! required for reciprocal space calculation
197  void set_ff_table(FormFactorTable* ff_table) { ff_table_ = ff_table; }
198 
199  void set_average_radius(double r) { average_radius_ = r; }
200 
201  void set_average_volume(double v) { average_volume_ = v; }
202 
203  void set_name(std::string name) { name_ = name; }
204 
205  void set_id(unsigned int id) { id_ = id; }
206 
207  void set_beam_profile(std::string beam_profile_file) {
208  beam_profile_ = new Profile(beam_profile_file);
209  }
210 
211  //! add simulated error
212  void add_errors();
213 
214  //! add simulated noise
215  void add_noise(double percentage = 0.03);
216 
217  //! computes full profile for given fitting parameters
218  void sum_partial_profiles(double c1, double c2, bool check_cashed = true);
219 
220  //! add another profile - useful for rigid bodies
221  void add(const Profile* other_profile, double weight = 1.0);
222 
223  //! add partial profiles
224  void add_partial_profiles(const Profile* other_profile, double weight = 1.0);
225 
226  //! add other profiles - useful for weighted ensembles
227  void add(const Vector<Profile*>& profiles,
228  const Vector<double>& weights = Vector<double>());
229 
230  //! add other partial profiles
231  void add_partial_profiles(const Vector<Profile*>& profiles,
232  const Vector<double>& weights = Vector<double>());
233 
234  //! background adjustment option
235  void background_adjust(double start_q);
236 
237  //! scale
238  void scale(double c);
239 
240  //! offset profile by c, I(q) = I(q) - c
241  void offset(double c);
242 
243  // copy error bars from the matching experimental profile
244  void copy_errors(const Profile* exp_profile);
245 
246  // parameter for E^2(q), used in faster calculation
247  static const double modulation_function_parameter_;
248 
250 
251  protected:
252  void init(unsigned int size = 0, unsigned int partial_profiles_size = 0);
253 
254  private:
255  void calculate_profile_reciprocal(const Particles& particles,
256  FormFactorType ff_type = HEAVY_ATOMS);
257 
258  void calculate_profile_reciprocal(const Particles& particles1,
259  const Particles& particles2,
260  FormFactorType ff_type = HEAVY_ATOMS);
261 
262  void calculate_profile_real(const Particles& particles,
263  FormFactorType ff_type = HEAVY_ATOMS);
264 
265  void calculate_profile_real(const Particles& particles1,
266  const Particles& particles2,
267  FormFactorType ff_type = HEAVY_ATOMS);
268 
269  void squared_distribution_2_profile(const RadialDistributionFunction& r_dist);
270 
271  void squared_distributions_2_partial_profiles(
272  const Vector<RadialDistributionFunction>& r_dist);
273 
274  double radius_of_gyration_fixed_q(double end_q) const;
275 
276  double find_max_q(const std::string& file_name) const;
277 
278  protected:
279  Eigen::VectorXf q_; // q sampling points
280  Eigen::VectorXf intensity_;
281  Eigen::VectorXf error_; // error bar of each point
282 
283  double min_q_, max_q_; // minimal and maximal q values in the profile
284  double delta_q_; // profile sampling resolution
285  FormFactorTable* ff_table_; // pointer to form factors table
286 
287  // stores the intensity split into 6 for c1/c2 enumeration
288  std::vector<Eigen::VectorXf> partial_profiles_;
289  double c1_, c2_;
290 
291  bool experimental_; // experimental profile read from file
292  double average_radius_; // average radius of the particles
293  double average_volume_; // average volume
294 
295  // mapping from q values to vector index for fast profile resampling
296  std::map<double, unsigned int> q_mapping_;
297 
298  std::string name_; // file name
299  unsigned int id_; // identifier
300 
301  Profile* beam_profile_;
302 };
303 
305 
306 IMPSAXS_END_NAMESPACE
307 
308 #endif /* IMPSAXS_PROFILE_H */
unsigned int size() const
return number of entries in SAXS profile
Definition: Profile.h:178
void calculate_profile(const Particles &particles, FormFactorType ff_type=HEAVY_ATOMS, bool reciprocal=false)
computes theoretical profile
Definition: Profile.h:50
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
Common base class for heavy weight IMP objects.
Definition: Object.h:106
Various general useful macros for IMP.
#define IMP_UNUSED(variable)
A class for computation of atomic and residue level form factors for SAXS calculations.
double get_delta_q() const
return sampling resolution
Definition: Profile.h:155
#define IMP_OBJECTS(Name, PluralName)
Define the types for storing lists of object pointers.
Definition: object_macros.h:44
double get_min_q() const
return minimal sampling point
Definition: Profile.h:158
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:197
void calculate_profile(const Particles &particles1, const Particles &particles2, FormFactorType ff_type=HEAVY_ATOMS)
Definition: Profile.h:84
computes distribution functions
double get_max_q() const
return maximal sampling point
Definition: Profile.h:161
double radius_of_gyration(const Particles &particles)
compute radius_of_gyration
Definition: saxs/utility.h:74