IMP  2.0.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-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 "FormFactorTable.h"
14 #include "Distribution.h"
15 
16 #include <iostream>
17 #include <vector>
18 
19 IMPSAXS_BEGIN_NAMESPACE
20 
21 class RadialDistributionFunction;
22 
23 /**
24  Basic profile class, can be initialized from the input file
25  (experimental or theoretical) or computed from a set of Model
26  Particles (theoretical)
27 */
28 class IMPSAXSEXPORT Profile {
29 public:
30  //! init from file
31  Profile(const String& file_name);
32 
33  //! init for theoretical profile
34  Profile(Float qmin = 0.0, Float qmax = 0.5, Float delta = 0.005);
35 
36 private:
37  class IntensityEntry {
38  public:
39  IntensityEntry() : q_(0.0), intensity_(0.0), error_(1.0), weight_(1.0) {}
40  IntensityEntry(Float q) : q_(q),intensity_(0.0),error_(1.0),weight_(1.0) {}
41  IntensityEntry(Float q, Float intensity, Float error)
42  : q_(q), intensity_(intensity), error_(error), weight_(1.0) {}
43 
44  Float q_;
45  Float intensity_;
46  Float error_;
47  Float weight_;
48  };
49 
50  friend std::ostream& operator<<(std::ostream& q, const IntensityEntry& e);
51 
52  friend std::istream& operator>>(std::istream& q, IntensityEntry& e);
53 
54 public:
55  //! computes theoretical profile
56  void calculate_profile(const Particles& particles,
57  FormFactorType ff_type = HEAVY_ATOMS,
58  bool reciprocal=false,
59  bool variance=false,
60  double variance_tau=0.1) {
61  IMP_USAGE_CHECK(!(reciprocal && variance),
62  "variance not implemented in reciprocal calculation");
63  if(!reciprocal) calculate_profile_real(particles, ff_type,
64  variance, variance_tau);
65  else calculate_profile_reciprocal(particles, ff_type);
66  }
67 
68  //! compute profile for fitting with hydration layer and excluded volume
69  void calculate_profile_partial(const Particles& particles,
70  const Floats& surface = Floats(),
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 Floats& surface1 = Floats(),
77  const Floats& surface2 = Floats(),
78  FormFactorType ff_type = HEAVY_ATOMS);
79 
80  void calculate_profile_reciprocal_partial(const Particles& particles,
81  const Floats& surface = Floats(),
82  FormFactorType ff_type = HEAVY_ATOMS);
83 
84  //! computes theoretical profile contribution from iter-molecular
85  //! interactions between the particles
86  void calculate_profile(const Particles& particles1,
87  const Particles& particles2,
88  FormFactorType ff_type = HEAVY_ATOMS,
89  bool variance=false,
90  double variance_tau=0.1) {
91  calculate_profile_real(particles1, particles2, ff_type,
92  variance, variance_tau);
93  }
94 
95  //! calculate Intensity at zero (= squared number of electrons)
96  Float calculate_I0(const Particles& particles,
97  FormFactorType ff_type = HEAVY_ATOMS);
98 
99  //! calculate profile for any type of Particles that have coordinates
100  void calculate_profile_constant_form_factor(const Particles& particles,
101  Float form_factor = 1.0);
102 
103  // computes theoretical profile faster for cyclically symmetric particles
104  // assumes that the units particles are ordered one after another in the
105  // input particles vector (n - symmetry order)
106  void calculate_profile_symmetric(const Particles& particles, 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  Float 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  //! add another profile - useful for rigid bodies
117  void add(const Profile& other_profile, Float weight = 1.0);
118 
119  //! add partial profiles
120  void add_partial_profiles(const Profile& other_profile, Float weight = 1.0);
121 
122  //! background adjustment option
123  void background_adjust(double start_q);
124 
125  //! scale
126  void scale(Float c);
127 
128  //! offset profile by c, I(q) = I(q) - c
129  void offset(Float c);
130 
131  //! reads SAXS profile from file
132  void read_SAXS_file(const String& file_name);
133 
134  //! print to file
135  /** \param[in] file_name output file name
136  \param[in] max_q output till maximal q value = max_q, or all if max_q<=0
137  */
138  void write_SAXS_file(const String& file_name, Float max_q=0.0) const;
139 
140  void write_partial_profiles(const String& file_name) const;
141 
142  // compute radius of gyration with Guinier approximation
143  // ln[I(q)]=ln[I(0)] - (q^2*rg^2)/3
144  // end_q_rg determines the range of profile used for approximation:
145  // i.e. q*rg < end_q_rg. Use 1.3 for globular proteins, 0.8 for elongated
146  double radius_of_gyration(double end_q_rg = 1.3) const;
147 
148  //! return sampling resolution
149  Float get_delta_q() const { return delta_q_; }
150 
151  //! return minimal sampling point
152  Float get_min_q() const { return min_q_; }
153 
154  //! return maximal sampling point
155  Float get_max_q() const { return max_q_; }
156 
157  //! return number of entries in SAXS profile
158  unsigned int size() const { return profile_.size(); }
159 
160  // Profile access functions
161  Float get_intensity(unsigned int i) const { return profile_[i].intensity_; }
162  Float get_q(unsigned int i) const { return profile_[i].q_; }
163  Float get_error(unsigned int i) const { return profile_[i].error_; }
164  Float get_weight(unsigned int i) const { return profile_[i].weight_; }
165  Float get_variance(unsigned int i, unsigned int j) const
166  { unsigned a=std::min(i,j); unsigned b=std::max(i,j);
167  return variances_[a][b-a]; }
168 
169  Float get_average_radius() const { return average_radius_; }
170 
171  void set_intensity(unsigned int i, Float iq) { profile_[i].intensity_ = iq; }
172 
173  //! required for reciprocal space calculation
174  void set_ff_table(FormFactorTable* ff_table) { ff_table_ = ff_table; }
175 
176  void set_average_radius(Float r) { average_radius_ = r; }
177 
178  void set_average_volume(Float v) { average_volume_ = v; }
179 
180  //! add intensity entry to profile
181  void add_entry(Float q, Float intensity, Float error=1.0) {
182  profile_.push_back(IntensityEntry(q, intensity, error));
183  }
184 
185  //! checks the sampling of experimental profile
186  bool is_uniform_sampling() const;
187 
188  //! add simulated error
189  void add_errors();
190 
191  //! add simulated noise
192  void add_noise(Float percentage = 0.03);
193 
194  //! computes full profile for given fitting parameters
195  void sum_partial_profiles(Float c1, Float c2, Profile& out_profile);
196 
197  // parameter for E^2(q), used in faster calculation
198  static const Float modulation_function_parameter_;
199 
200  private:
201  void init(bool variance = false);
202 
203  void calculate_profile_reciprocal(const Particles& particles,
204  FormFactorType ff_type = HEAVY_ATOMS);
205 
206  void calculate_profile_reciprocal(const Particles& particles1,
207  const Particles& particles2,
208  FormFactorType ff_type = HEAVY_ATOMS);
209 
210  void calculate_profile_real(const Particles& particles,
211  FormFactorType ff_type = HEAVY_ATOMS,
212  bool variance = false,
213  double variance_tau = 0.1);
214 
215  void calculate_profile_real(const Particles& particles1,
216  const Particles& particles2,
217  FormFactorType ff_type = HEAVY_ATOMS,
218  bool variance = false,
219  double variance_tau = 0.1);
220 
221  void squared_distribution_2_profile(const RadialDistributionFunction& r_dist,
222  const RadialDistributionFunction& r_dist2,
223  bool variance=false, double variance_tau=0.1);
224 
225  void squared_distributions_2_partial_profiles(
226  const std::vector<RadialDistributionFunction>& r_dist);
227 
228  double radius_of_gyration_fixed_q(double end_q) const;
229 
230  protected:
231  std::vector<IntensityEntry> profile_; // the profile
232  std::vector<std::vector<double> > variances_; //profile variances
233  Float min_q_, max_q_; // minimal and maximal s values in the profile
234  Float delta_q_; // profile sampling resolution
235  FormFactorTable* ff_table_; // pointer to form factors table
236  std::vector<Profile> partial_profiles_;
237  bool experimental_; // experimental profile read from file
238  Float average_radius_; // average radius of the particles
239  Float average_volume_; // average volume
240 };
241 
242 IMPSAXS_END_NAMESPACE
243 
244 #endif /* IMPSAXS_PROFILE_H */