IMP logo
IMP Reference Guide  2.6.1
The Integrative Modeling Platform
WeightedProfileFitter.h
Go to the documentation of this file.
1 /**
2  * \file IMP/saxs/WeightedProfileFitter.h
3  * \brief Fitting of multiple profiles to the experimental one.
4  * The weights of the profiles are computed analytically using
5  * non-negative least squares fitting (NNLS)
6  *
7  * \authors Dina Schneidman
8  * Copyright 2007-2016 IMP Inventors. All rights reserved.
9  *
10  */
11 
12 #ifndef IMPSAXS_WEIGHTED_PROFILE_FITTER_H
13 #define IMPSAXS_WEIGHTED_PROFILE_FITTER_H
14 
15 #include "ProfileFitter.h"
16 #include "ChiScore.h"
17 #include "WeightedFitParameters.h"
18 #include "nnls.h"
19 #include <IMP/algebra/eigen3/Eigen/Dense>
20 
21 IMPSAXS_BEGIN_NAMESPACE
22 
23 //! Fitting of multiple profiles to the experimental one.
24 /** The weights of the profiles are computed analytically using
25  non-negative least squares fitting (NNLS).
26 */
27 template <typename ScoringFunctionT = ChiScore>
28 class WeightedProfileFitter : public ProfileFitter<ScoringFunctionT> {
29 
30  public:
31  //! Constructor
32  /**
33  \param[in] exp_profile Experimental profile we want to fit
34  */
35  WeightedProfileFitter(const Profile* exp_profile) :
36  ProfileFitter<ScoringFunctionT>(exp_profile),
37  W_(exp_profile->size(), 1),
38  Wb_(exp_profile->size()),
39  A_(exp_profile->size(), 2) {
40 
41  IMP_Eigen::VectorXf b(exp_profile->size());
42  for (unsigned int i = 0; i < exp_profile->size(); i++) {
43  Wb_(i) = exp_profile->get_intensity(i);
44  W_(i) = 1.0 / (exp_profile->get_error(i));
45  }
46  Wb_ = W_.asDiagonal() * Wb_;
47  }
48 
49  //! compute a weighted score that minimizes chi
50  /**
51  it is assumed that the q values of the profiles are the same as
52  the q values of the experimental profile. Use Profile::resample to resample
53  if(NNLS = true, solve non-negative least squares, otherwise solve just
54  least squares, that may return negative weights to be discarded later
55  */
56  double compute_score(const ProfilesTemp& profiles,
57  Vector<double>& weights, bool NNLS = true) const;
58 
59  //! fit profiles by optimization of c1/c2 and weights
60  /**
61  it is assumed that the q values of the profiles are the same as
62  the q values of the experimental profile. Use Profile::resample to resample
63  */
64  WeightedFitParameters fit_profile(ProfilesTemp partial_profiles,
65  double min_c1 = 0.95, double max_c1 = 1.05,
66  double min_c2 = -2.0,
67  double max_c2 = 4.0) const;
68 
69  //! write a fit file
70  void write_fit_file(ProfilesTemp partial_profiles,
71  const WeightedFitParameters& fp,
72  const std::string fit_file_name) const;
73 
74  private:
75  WeightedFitParameters search_fit_parameters(
76  ProfilesTemp& partial_profiles, double min_c1, double max_c1, double min_c2,
77  double max_c2, double old_chi, Vector<double>& weights) const;
78 
79  private:
80  IMP_Eigen::MatrixXf W_; // weights matrix
81 
82  // weights matrix multiplied by experimental intensities vector
83  IMP_Eigen::VectorXf Wb_;
84 
85  // intensities
86  IMP_Eigen::MatrixXf A_;
87 };
88 
89 template <typename ScoringFunctionT>
91  const ProfilesTemp& profiles,
92  Vector<double>& weights,
93  bool nnls) const {
94 
95  // no need to compute weighted profile for ensemble of size 1
96  if (profiles.size() == 1) {
97  weights.resize(1);
98  weights[0] = 1.0;
101  profiles[0]);
102  }
103 
104  // compute_weights(profiles, weights);
105  int m = profiles.size();
107 
108  // fill in A_
109  WeightedProfileFitter* non_const_this =
110  const_cast<WeightedProfileFitter*>(this);
111 
112  if (A_.cols() != m) non_const_this->A_.resize(n, m);
113  for (int j = 0; j < m; j++) {
114  for (int i = 0; i < n; i++) {
115  (non_const_this->A_)(i, j) = profiles[j]->get_intensity(i);
116  }
117  }
118 
119  IMP_Eigen::VectorXf w;
120  if (!nnls) { // solve least squares
121  IMP_Eigen::JacobiSVD<IMP_Eigen::MatrixXf> svd(W_.asDiagonal() * A_,
122  ComputeThinU | ComputeThinV);
123  w = svd.solve(Wb_);
124  // zero the negatives
125  for (int i = 0; i < w.size(); i++)
126  if (w(i) < 0) w(i) = 0;
127  } else {
128  w = NNLS(W_.asDiagonal() * A_, Wb_);
129  }
130  w /= w.sum();
131 
132  // compute weighted profile
133  IMP_NEW(Profile, weighted_profile,
137  IMP_Eigen::VectorXf wp = A_ * w;
138  weighted_profile->set_qs(profiles[0]->get_qs());
139  weighted_profile->set_intensities(wp);
140  weights.resize(w.size());
141  for (int i = 0; i < w.size(); i++) weights[i] = w[i];
142 
145  weighted_profile);
146 }
147 
148 template <typename ScoringFunctionT>
150  ProfilesTemp partial_profiles,
151  double min_c1, double max_c1,
152  double min_c2, double max_c2) const {
153  Vector<double> weights;
155  search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2, max_c2,
156  std::numeric_limits<double>::max(), weights);
157  return fp;
158 }
159 
160 template <typename ScoringFunctionT>
162  ProfilesTemp partial_profiles,
163  const WeightedFitParameters& fp,
164  const std::string fit_file_name) const {
165  double best_c1 = fp.get_c1();
166  double best_c2 = fp.get_c2();
167 
168  // compute a profile for best c1/c2 combination
169  for (unsigned int i = 0; i < partial_profiles.size(); i++)
170  partial_profiles[i]->sum_partial_profiles(best_c1, best_c2);
171 
172  if (partial_profiles.size() == 1) {
173  // compute scale
174  double c =
177  partial_profiles[0]);
179  fit_file_name, partial_profiles[0], fp.get_chi(), c);
180  } else {
181 
182  // computed weighted profile
183  IMP_NEW(Profile, weighted_profile,
187 
188  const Vector<double>& weights = fp.get_weights();
189  for (unsigned int i = 0; i < weights.size(); i++)
190  weighted_profile->add(partial_profiles[i], weights[i]);
191 
192  // compute scale
193  double c =
196  weighted_profile);
197 
199  fit_file_name, weighted_profile, fp.get_chi(), c);
200  }
201 }
202 
203 template <typename ScoringFunctionT>
205  ProfilesTemp& partial_profiles,
206  double min_c1, double max_c1,
207  double min_c2, double max_c2,
208  double old_chi, Vector<double>& weights) const {
209  int c1_cells = 10;
210  int c2_cells = 10;
211  if (old_chi < (std::numeric_limits<double>::max() - 1)) { // second iteration
212  c1_cells = 5;
213  c2_cells = 5;
214  }
215 
216  double delta_c1 = (max_c1 - min_c1) / c1_cells;
217  double delta_c2 = (max_c2 - min_c2) / c2_cells;
218 
219  bool last_c1 = false;
220  bool last_c2 = false;
221  if (delta_c1 < 0.0001) {
222  c1_cells = 1;
223  delta_c1 = max_c1 - min_c1;
224  last_c1 = true;
225  }
226  if (delta_c2 < 0.001) {
227  c2_cells = 1;
228  delta_c2 = max_c2 - min_c2;
229  last_c2 = true;
230  }
231  double best_c1(1.0), best_c2(0.0), best_chi(old_chi);
232  bool best_set = false;
233 
234  // c1 iteration
235  double c1(min_c1);
236  for (int i = 0; i <= c1_cells; i++, c1 += delta_c1) {
237  // c2 iteration
238  double c2 = min_c2;
239  for (int j = 0; j <= c2_cells; j++, c2 += delta_c2) {
240  // sum up the profiles for c1/c2 combo
241  for (unsigned int k = 0; k < partial_profiles.size(); k++)
242  partial_profiles[k]->sum_partial_profiles(c1, c2);
243  Vector<double> curr_weights;
244  double curr_chi = compute_score(partial_profiles, curr_weights);
245  if (!best_set || curr_chi < best_chi) {
246  best_set = true;
247  best_chi = curr_chi;
248  best_c1 = c1;
249  best_c2 = c2;
250  weights = curr_weights;
251  }
252  }
253  }
254 
255  if (std::fabs(best_chi - old_chi) > 0.005 &&
256  (!(last_c1 && last_c2))) { // refine more
257  min_c1 = std::max(best_c1 - delta_c1, min_c1);
258  max_c1 = std::min(best_c1 + delta_c1, max_c1);
259  min_c2 = std::max(best_c2 - delta_c2, min_c2);
260  max_c2 = std::min(best_c2 + delta_c2, max_c2);
261  return search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2,
262  max_c2, best_chi, weights);
263  }
264  return WeightedFitParameters(best_chi, best_c1, best_c2, weights);
265 }
266 
267 IMPSAXS_END_NAMESPACE
268 
269 #endif /* IMPSAXS_WEIGHTED_PROFILE_FITTER_H */
unsigned int size() const
return number of entries in SAXS profile
Definition: Profile.h:171
Basic chi score implementation.
WeightedProfileFitter(const Profile *exp_profile)
Constructor.
a class for fitting two profiles
#define IMP_NEW(Typename, varname, args)
Declare a ref counted pointer to a new object.
Definition: object_macros.h:62
Fitting of multiple profiles to the experimental one.
IMP_Eigen::VectorXf NNLS(const IMP_Eigen::MatrixXf &A, const IMP_Eigen::VectorXf &b)
non-negative least square fitting for profile weight solving problem
Definition: nnls.h:21
Parameters of a weighted fit, from WeightedProfileFitter.
Fit two profiles with user-defined scoring function as a template parameter.
Definition: ProfileFitter.h:29
double compute_score(const Profile *model_profile, bool use_offset=false, const std::string fit_file_name="") const
compute fit score
double compute_scale_factor(const Profile *model_profile, double offset=0.0) const
computes the scaling factor needed for fitting the modeled profile
Definition: ProfileFitter.h:76