IMP  2.1.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-2013 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): base::Object("ProfileFitter%1%"),
37  exp_profile_(exp_profile) {
38  set_was_used(true);
39  scoring_function_ = new ScoringFunctionT();
40  }
41 
42  ProfileFitter(const Profile* exp_profile,
43  ScoringFunctionT* sf): base::Object("ProfileFitter%1%"),
44  exp_profile_(exp_profile) {
45  set_was_used(true);
46  scoring_function_ = sf;
47  }
48 
49  //! compute fit score
50  Float compute_score(const Profile* model_profile,
51  bool use_offset = false,
52  const std::string fit_file_name = "") const;
53 
54  //! compute fit score in the interval
55  Float compute_score(const Profile* model_profile,
56  Float min_q, Float max_q) const;
57 
58  //! fit experimental profile through optimization of c1 and c2 parameters
59  /**
60  c1 - adjusts the excluded volume, valid range [0.95 - 1.05]
61  c2 - adjusts the density of hydration layer, valid range [-2.0 - 4.0]
62  \param[in] partial_profile partial profiles computed
63  \param[in] min_c1 minimal c1 value
64  \param[in] max_c1 maximal c1 value
65  \param[in] min_c2 minimal c2 value
66  \param[in] max_c2 maximal c2 value
67  \param[in] use_offset use offset in fitting
68  \param[in] fit_file_name write fit in the given filename,
69  nothing is written if empty
70  \return FitParameters (score, c1, c2)
71  */
72  FitParameters fit_profile(Profile* partial_profile,
73  float min_c1=0.95, float max_c1=1.05,
74  float min_c2=-2.0, float max_c2=4.0,
75  bool use_offset = false,
76  const std::string fit_file_name = "") const;
77 
78  //! computes the scaling factor needed for fitting the modeled profile
79  // onto the experimental one
80  /** resampling of the modeled profile is required prior to calling
81  to this function, in order to match the q values of the exp. profile.
82  this is not done by default to minimize the number of calls to resample.
83  */
84  Float compute_scale_factor(const Profile* model_profile,
85  Float offset = 0.0) const {
86  return scoring_function_->compute_scale_factor(exp_profile_, model_profile,
87  offset);
88  }
89 
90  //! computes offset for fitting for fitting the modeled profile
91  // onto the experimental one
92  /** resampling of the modeled profile is required prior to calling
93  to this function, in order to match the q values of the exp. profile.
94  this is not done by default to minimize the number of calls to resample.
95  */
96  Float compute_offset(const Profile* model_profile) const {
97  return scoring_function_->compute_offset(exp_profile_, model_profile);
98  }
99 
100  //! resampling of the modeled profile is required to fit the q
101  // values of the computational profile to the experimental one
102  void resample(const Profile* model_profile, Profile* resampled_profile) const;
103 
104  // writes 3 column fit file, given scale factor c, offset and chi square
105  void write_SAXS_fit_file(const std::string& file_name,
106  const Profile* model_profile,
107  const Float chi_square,
108  const Float c=1, const Float offset=0) const;
109 
110  private:
111  FitParameters search_fit_parameters(Profile* partial_profile,
112  float min_c1, float max_c1,
113  float min_c2, float max_c2,
114  bool use_offset, float old_chi) const;
115 
116 
118  friend class DerivativeCalculator;
119 
120  protected:
121  // experimental saxs profile
123  ScoringFunctionT* scoring_function_;
124 };
125 
126 template<class ScoringFunctionT>
128  Profile* resampled_profile) const
129 {
130  model_profile->resample(exp_profile_, resampled_profile);
131 }
132 
133 template<class ScoringFunctionT>
135  Profile* partial_profile,
136  float min_c1, float max_c1,
137  float min_c2, float max_c2,
138  bool use_offset, float old_chi) const
139 {
140  int c1_cells = 10;
141  int c2_cells = 10;
142  if(old_chi < (std::numeric_limits<float>::max()-1)) { // second iteration
143  c1_cells = 5;
144  c2_cells = 5;
145  }
146 
147  float delta_c1 = (max_c1-min_c1)/c1_cells;
148  float delta_c2 = (max_c2-min_c2)/c2_cells;
149 
150  bool last_c1 = false;
151  bool last_c2 = false;
152  if(delta_c1 < 0.0001) { c1_cells=1; delta_c1 = max_c1-min_c1; last_c1=true; }
153  if(delta_c2 < 0.001) { c2_cells=1; delta_c2 = max_c2-min_c2; last_c2=true; }
154  float best_c1(1.0), best_c2(0.0), best_chi(old_chi);
155  bool best_set = false;
156 
157  float c1(min_c1);
158  for(int i=0; i<=c1_cells; i++, c1+= delta_c1) {
159  float c2 = min_c2;
160  for(int j=0; j<=c2_cells; j++, c2+= delta_c2) {
161  partial_profile->sum_partial_profiles(c1, c2);
162  float curr_chi = compute_score(partial_profile, use_offset);
163  if(!best_set || curr_chi < best_chi) {
164  best_set = true;
165  best_chi = curr_chi;
166  best_c1 = c1;
167  best_c2 = c2;
168  }
169  }
170  }
171 
172  if(std::fabs(best_chi-old_chi) > 0.0001 &&
173  (!(last_c1 && last_c2))) { //refine more
174  min_c1 = std::max(best_c1-delta_c1, min_c1);
175  max_c1 = std::min(best_c1+delta_c1, max_c1);
176  min_c2 = std::max(best_c2-delta_c2, min_c2);
177  max_c2 = std::min(best_c2+delta_c2, max_c2);
178  return search_fit_parameters(partial_profile,
179  min_c1, max_c1, min_c2, max_c2,
180  use_offset, best_chi);
181  }
182  return FitParameters(best_chi, best_c1, best_c2);
183 }
184 
185 template<class ScoringFunctionT>
187  Profile* partial_profile,
188  float min_c1, float max_c1,
189  float min_c2, float max_c2,
190  bool use_offset,
191  const std::string fit_file_name) const {
192 
193  // compute chi value for default c1/c1 (remove?)
194  float default_c1 = 1.0, default_c2 = 0.0;
195  partial_profile->sum_partial_profiles(default_c1, default_c2);
196  float default_chi = compute_score(partial_profile, use_offset);
197 
198  FitParameters fp = search_fit_parameters(partial_profile,
199  min_c1, max_c1, min_c2, max_c2,
200  use_offset,
201  std::numeric_limits<float>::max());
202  float best_c1 = fp.get_c1();
203  float best_c2 = fp.get_c2();
204  fp.set_default_chi(default_chi);
205 
206  // compute a profile for best c1/c2 combination
207  partial_profile->sum_partial_profiles(best_c1, best_c2);
208  compute_score(partial_profile, use_offset, fit_file_name);
209 
210  // std::cout << " Chi = " << best_chi << " c1 = " << best_c1 << " c2 = "
211  // << best_c2 << " default chi = " << default_chi << std::endl;
212  return fp;
213 }
214 
215 template<class ScoringFunctionT>
217  const Profile* model_profile,
218  Float min_q, Float max_q) const
219 {
220  IMP_NEW(Profile, resampled_profile, (exp_profile_->get_min_q(),
221  exp_profile_->get_max_q(),
222  exp_profile_->get_delta_q()));
223  model_profile->resample(exp_profile_, resampled_profile);
224 
225  Float score = scoring_function_->compute_score(exp_profile_,
226  resampled_profile,
227  min_q, max_q);
228  return score;
229 }
230 
231 template<class ScoringFunctionT>
233  const Profile* model_profile,
234  bool use_offset,
235  const std::string fit_file_name) const
236 {
237  IMP_NEW(Profile, resampled_profile, (exp_profile_->get_min_q(),
238  exp_profile_->get_max_q(),
239  exp_profile_->get_delta_q()));
240  model_profile->resample(exp_profile_, resampled_profile);
241 
242  Float score = scoring_function_->compute_score(exp_profile_,
243  resampled_profile,
244  use_offset);
245 
246  if(fit_file_name.length() > 0) {
247  Float offset = 0.0;
248  if(use_offset)
249  offset=scoring_function_->compute_offset(exp_profile_,
250  resampled_profile);
251  Float c= scoring_function_->compute_scale_factor(exp_profile_,
252  resampled_profile,
253  offset);
254  write_SAXS_fit_file(fit_file_name, resampled_profile, score, c, offset);
255  }
256  return score;
257 }
258 
259 template<class ScoringFunctionT>
261  const std::string& file_name,
262  const Profile* model_profile,
263  const Float score,
264  const Float c,
265  const Float offset) const {
266  std::ofstream out_file(file_name.c_str());
267  if (!out_file) {
268  IMP_THROW("Can't open file " << file_name,
269  IOException);
270  }
271 
272  unsigned int profile_size = std::min(model_profile->size(),
273  exp_profile_->size());
274  // header line
275  out_file.precision(15);
276  out_file << "# SAXS profile: number of points = " << profile_size
277  << ", q_min = " << exp_profile_->get_min_q()
278  << ", q_max = " << exp_profile_->get_max_q();
279  out_file << ", delta_q = " << exp_profile_->get_delta_q() << std::endl;
280 
281  out_file.setf(std::ios::showpoint);
282  out_file << "# offset = " << offset << ", scaling c = " << c
283  << ", Chi = " << score << std::endl;
284  out_file << "# q exp_intensity model_intensity"
285  << std::endl;
286 
287  out_file.setf(std::ios::fixed, std::ios::floatfield);
288  // Main data
289  for (unsigned int i = 0; i < profile_size; i++) {
290  out_file.setf(std::ios::left);
291  out_file.width(10);
292  out_file.precision(5);
293  out_file << exp_profile_->get_q(i) << " ";
294 
295  out_file.setf(std::ios::left);
296  out_file.width(15);
297  out_file.precision(8);
298  out_file << exp_profile_->get_intensity(i) << " ";
299 
300  out_file.setf(std::ios::left);
301  out_file.width(15);
302  out_file.precision(8);
303  out_file << model_profile->get_intensity(i)*c - offset << std::endl;
304  }
305  out_file.close();
306 }
307 
308 IMPSAXS_END_NAMESPACE
309 
310 #endif /* IMPSAXS_PROFILE_FITTER_H */
unsigned int size() const
return number of entries in SAXS profile
Definition: Profile.h:173
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:146
#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:84
void sum_partial_profiles(Float c1, Float c2)
computes full profile for given fitting parameters
Copyright 2007-2013 IMP Inventors. All rights reserved.
Common base class for heavy weight IMP objects.
Float compute_offset(const Profile *model_profile) const
computes offset for fitting for fitting the modeled profile
Definition: ProfileFitter.h:96
#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.