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