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