IMP  2.2.0
The Integrative Modeling Platform
ProfileFitter.h
Go to the documentation of this file.
1 /**
2  * \file IMP/saxs/ProfileFitter.h \brief a class for fitting two profiles
3  *
4  * Copyright 2007-2014 IMP Inventors. All rights reserved.
5  *
6  */
7 
8 #ifndef IMPSAXS_PROFILE_FITTER_H
9 #define IMPSAXS_PROFILE_FITTER_H
10 
11 #include "ChiScore.h"
12 #include "FitParameters.h"
13 #include "Profile.h"
14 #include <IMP/base/Object.h>
15 
16 #include <fstream>
17 
18 IMPSAXS_BEGIN_NAMESPACE
19 
20 /**
21  ProfileFitter is a class for fitting the two profiles with user-defined
22  scoring function that is a template parameter. By default chi score is used.
23  The scoring function template parameter class has to implement three
24  basic functions: compute_score, compute_scale_factor and compute_offset.
25  see ChiScore for example.
26  Currently three scoring functions are implemented:
27  ChiScore, ChiScoreLog and ChiFreeScore.
28  */
29 template <class ScoringFunctionT = ChiScore>
30 class ProfileFitter : public base::Object {
31  public:
32  //! Constructor
33  /**
34  \param[in] exp_profile Experimental profile we want to fit
35  */
36  ProfileFitter(const Profile* exp_profile)
37  : base::Object("ProfileFitter%1%"), exp_profile_(exp_profile) {
38  set_was_used(true);
39  scoring_function_ = new ScoringFunctionT();
40  }
41 
42  ProfileFitter(const Profile* exp_profile, ScoringFunctionT* sf)
43  : base::Object("ProfileFitter%1%"), exp_profile_(exp_profile) {
44  set_was_used(true);
45  scoring_function_ = sf;
46  }
47 
48  //! compute fit score
49  Float compute_score(const Profile* model_profile, bool use_offset = false,
50  const std::string fit_file_name = "") const;
51 
52  //! compute fit score in the interval
53  Float compute_score(const Profile* model_profile, Float min_q,
54  Float max_q) const;
55 
56  //! fit experimental profile through optimization of c1 and c2 parameters
57  /**
58  c1 - adjusts the excluded volume, valid range [0.95 - 1.05]
59  c2 - adjusts the density of hydration layer, valid range [-2.0 - 4.0]
60  \param[in] partial_profile partial profiles computed
61  \param[in] min_c1 minimal c1 value
62  \param[in] max_c1 maximal c1 value
63  \param[in] min_c2 minimal c2 value
64  \param[in] max_c2 maximal c2 value
65  \param[in] use_offset use offset in fitting
66  \param[in] fit_file_name write fit in the given filename,
67  nothing is written if empty
68  \return FitParameters (score, c1, c2)
69  */
70  FitParameters fit_profile(Profile* partial_profile, float min_c1 = 0.95,
71  float max_c1 = 1.05, float min_c2 = -2.0,
72  float max_c2 = 4.0, bool use_offset = false,
73  const std::string fit_file_name = "") const;
74 
75  //! computes the scaling factor needed for fitting the modeled profile
76  // onto the experimental one
77  /** resampling of the modeled profile is required prior to calling
78  to this function, in order to match the q values of the exp. profile.
79  this is not done by default to minimize the number of calls to resample.
80  */
81  Float compute_scale_factor(const Profile* model_profile,
82  Float offset = 0.0) const {
83  return scoring_function_->compute_scale_factor(exp_profile_, model_profile,
84  offset);
85  }
86 
87  //! computes offset for fitting for fitting the modeled profile
88  // onto the experimental one
89  /** resampling of the modeled profile is required prior to calling
90  to this function, in order to match the q values of the exp. profile.
91  this is not done by default to minimize the number of calls to resample.
92  */
93  Float compute_offset(const Profile* model_profile) const {
94  return scoring_function_->compute_offset(exp_profile_, model_profile);
95  }
96 
97  //! resampling of the modeled profile is required to fit the q
98  // values of the computational profile to the experimental one
99  void resample(const Profile* model_profile, Profile* resampled_profile) const;
100 
101  // writes 3 column fit file, given scale factor c, offset and chi square
102  void write_SAXS_fit_file(const std::string& file_name,
103  const Profile* model_profile, const Float chi_square,
104  const Float c = 1, const Float offset = 0) const;
105 
106  private:
107  FitParameters search_fit_parameters(Profile* partial_profile, float min_c1,
108  float max_c1, float min_c2, float max_c2,
109  bool use_offset, float old_chi) const;
110 
112  friend class DerivativeCalculator;
113 
114  protected:
115  // experimental saxs profile
117  ScoringFunctionT* scoring_function_;
118 };
119 
120 template <class ScoringFunctionT>
122  const Profile* model_profile, Profile* resampled_profile) const {
123  model_profile->resample(exp_profile_, resampled_profile);
124 }
125 
126 template <class ScoringFunctionT>
128  Profile* partial_profile, float min_c1, float max_c1, float min_c2,
129  float max_c2, bool use_offset, float old_chi) const {
130  int c1_cells = 10;
131  int c2_cells = 10;
132  if (old_chi < (std::numeric_limits<float>::max() - 1)) { // second iteration
133  c1_cells = 5;
134  c2_cells = 5;
135  }
136 
137  float delta_c1 = (max_c1 - min_c1) / c1_cells;
138  float delta_c2 = (max_c2 - min_c2) / c2_cells;
139 
140  bool last_c1 = false;
141  bool last_c2 = false;
142  if (delta_c1 < 0.0001) {
143  c1_cells = 1;
144  delta_c1 = max_c1 - min_c1;
145  last_c1 = true;
146  }
147  if (delta_c2 < 0.001) {
148  c2_cells = 1;
149  delta_c2 = max_c2 - min_c2;
150  last_c2 = true;
151  }
152  float best_c1(1.0), best_c2(0.0), best_chi(old_chi);
153  bool best_set = false;
154 
155  float c1(min_c1);
156  for (int i = 0; i <= c1_cells; i++, c1 += delta_c1) {
157  float c2 = min_c2;
158  for (int j = 0; j <= c2_cells; j++, c2 += delta_c2) {
159  partial_profile->sum_partial_profiles(c1, c2);
160  float curr_chi = compute_score(partial_profile, use_offset);
161  if (!best_set || curr_chi < best_chi) {
162  best_set = true;
163  best_chi = curr_chi;
164  best_c1 = c1;
165  best_c2 = c2;
166  }
167  }
168  }
169 
170  if (std::fabs(best_chi - old_chi) > 0.0001 &&
171  (!(last_c1 && last_c2))) { // refine more
172  min_c1 = std::max(best_c1 - delta_c1, min_c1);
173  max_c1 = std::min(best_c1 + delta_c1, max_c1);
174  min_c2 = std::max(best_c2 - delta_c2, min_c2);
175  max_c2 = std::min(best_c2 + delta_c2, max_c2);
176  return search_fit_parameters(partial_profile, min_c1, max_c1, min_c2,
177  max_c2, use_offset, best_chi);
178  }
179  return FitParameters(best_chi, best_c1, best_c2);
180 }
181 
182 template <class ScoringFunctionT>
184  Profile* partial_profile, float min_c1, float max_c1, float min_c2,
185  float max_c2, bool use_offset, const std::string fit_file_name) const {
186 
187  // compute chi value for default c1/c1 (remove?)
188  float default_c1 = 1.0, default_c2 = 0.0;
189  partial_profile->sum_partial_profiles(default_c1, default_c2);
190  float default_chi = compute_score(partial_profile, use_offset);
191 
192  FitParameters fp =
193  search_fit_parameters(partial_profile, min_c1, max_c1, min_c2, max_c2,
194  use_offset, std::numeric_limits<float>::max());
195  float best_c1 = fp.get_c1();
196  float best_c2 = fp.get_c2();
197  fp.set_default_chi(default_chi);
198 
199  // compute a profile for best c1/c2 combination
200  partial_profile->sum_partial_profiles(best_c1, best_c2);
201  compute_score(partial_profile, use_offset, fit_file_name);
202 
203  // std::cout << " Chi = " << best_chi << " c1 = " << best_c1 << " c2 = "
204  // << best_c2 << " default chi = " << default_chi << std::endl;
205  return fp;
206 }
207 
208 template <class ScoringFunctionT>
210  const Profile* model_profile, Float min_q, Float max_q) const {
211  IMP_NEW(Profile, resampled_profile,
212  (exp_profile_->get_min_q(), exp_profile_->get_max_q(),
213  exp_profile_->get_delta_q()));
214  model_profile->resample(exp_profile_, resampled_profile);
215 
216  Float score = scoring_function_->compute_score(
217  exp_profile_, resampled_profile, min_q, max_q);
218  return score;
219 }
220 
221 template <class ScoringFunctionT>
223  const Profile* model_profile, bool use_offset,
224  const std::string fit_file_name) const {
225  IMP_NEW(Profile, resampled_profile,
226  (exp_profile_->get_min_q(), exp_profile_->get_max_q(),
227  exp_profile_->get_delta_q()));
228  model_profile->resample(exp_profile_, resampled_profile);
229 
230  Float score = scoring_function_->compute_score(exp_profile_,
231  resampled_profile, use_offset);
232 
233  if (fit_file_name.length() > 0) {
234  Float offset = 0.0;
235  if (use_offset)
236  offset =
237  scoring_function_->compute_offset(exp_profile_, resampled_profile);
238  Float c = scoring_function_->compute_scale_factor(
239  exp_profile_, resampled_profile, offset);
240  write_SAXS_fit_file(fit_file_name, resampled_profile, score, c, offset);
241  }
242  return score;
243 }
244 
245 template <class ScoringFunctionT>
247  const std::string& file_name, const Profile* model_profile,
248  const Float score, const Float c, const Float offset) const {
249  std::ofstream out_file(file_name.c_str());
250  if (!out_file) {
251  IMP_THROW("Can't open file " << file_name, IOException);
252  }
253 
254  unsigned int profile_size =
255  std::min(model_profile->size(), exp_profile_->size());
256  // header line
257  out_file.precision(15);
258  out_file << "# SAXS profile: number of points = " << profile_size
259  << ", q_min = " << exp_profile_->get_min_q()
260  << ", q_max = " << exp_profile_->get_max_q();
261  out_file << ", delta_q = " << exp_profile_->get_delta_q() << std::endl;
262 
263  out_file.setf(std::ios::showpoint);
264  out_file << "# offset = " << offset << ", scaling c = " << c
265  << ", Chi = " << score << std::endl;
266  out_file << "# q exp_intensity model_intensity" << std::endl;
267 
268  out_file.setf(std::ios::fixed, std::ios::floatfield);
269  // Main data
270  for (unsigned int i = 0; i < profile_size; i++) {
271  out_file.setf(std::ios::left);
272  out_file.width(10);
273  out_file.precision(5);
274  out_file << exp_profile_->get_q(i) << " ";
275 
276  out_file.setf(std::ios::left);
277  out_file.width(15);
278  out_file.precision(8);
279  out_file << exp_profile_->get_intensity(i) << " ";
280 
281  out_file.setf(std::ios::left);
282  out_file.width(15);
283  out_file.precision(8);
284  out_file << model_profile->get_intensity(i) * c - offset << std::endl;
285  }
286  out_file.close();
287 }
288 
289 IMPSAXS_END_NAMESPACE
290 
291 #endif /* IMPSAXS_PROFILE_FITTER_H */
unsigned int size() const
return number of entries in SAXS profile
Definition: Profile.h:163
ProfileFitter(const Profile *exp_profile)
Constructor.
Definition: ProfileFitter.h:36
Basic chi score implementation.
void resample(const Profile *exp_profile, Profile *resampled_profile, bool partial_profiles=false) const
return a profile that is sampled on the q values of the exp_profile
A smart pointer to a ref-counted Object that is a class memeber.
Definition: base/Pointer.h:147
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Ref counted objects should have private destructors.
#define IMP_NEW(Typename, varname, args)
Declare a ref counted pointer to a new object.
Float compute_scale_factor(const Profile *model_profile, Float offset=0.0) const
computes the scaling factor needed for fitting the modeled profile
Definition: ProfileFitter.h:81
void sum_partial_profiles(Float c1, Float c2)
computes full profile for given fitting parameters
Copyright 2007-2014 IMP Inventors. All rights reserved.
Common base class for heavy weight IMP objects.
Definition: base/Object.h:106
Float compute_offset(const Profile *model_profile) const
computes offset for fitting for fitting the modeled profile
Definition: ProfileFitter.h:93
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
double Float
Basic floating-point value (could be float, double...)
Definition: base/types.h:20
A shared base class to help in debugging and things.
A class for profile storing and computation.