IMP logo
IMP Reference Guide  2.8.0
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-2017 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, Vector<double>& weights,
57  bool use_offset = false, 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, double max_c2 = 4.0,
67  bool use_offset = false) 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,
73  bool use_offset = false) const;
74 
75  private:
76  WeightedFitParameters search_fit_parameters(
77  ProfilesTemp& partial_profiles, double min_c1, double max_c1, double min_c2,
78  double max_c2, double old_chi, Vector<double>& weights, bool use_offset) const;
79 
80  private:
81  IMP_Eigen::MatrixXf W_; // weights matrix
82 
83  // weights matrix multiplied by experimental intensities vector
84  IMP_Eigen::VectorXf Wb_;
85 
86  // intensities
87  IMP_Eigen::MatrixXf A_;
88 };
89 
90 template <typename ScoringFunctionT>
92  const ProfilesTemp& profiles,
93  Vector<double>& weights,
94  bool use_offset, bool nnls) const {
95 
96  // no need to compute weighted profile for ensemble of size 1
97  if (profiles.size() == 1) {
98  weights.resize(1);
99  weights[0] = 1.0;
102  profiles[0], use_offset);
103  }
104 
105  // compute_weights(profiles, weights);
106  int m = profiles.size();
108 
109  // fill in A_
110  WeightedProfileFitter* non_const_this =
111  const_cast<WeightedProfileFitter*>(this);
112 
113  if (A_.cols() != m) non_const_this->A_.resize(n, m);
114  for (int j = 0; j < m; j++) {
115  for (int i = 0; i < n; i++) {
116  (non_const_this->A_)(i, j) = profiles[j]->get_intensity(i);
117  }
118  }
119 
120  IMP_Eigen::VectorXf w;
121  if (!nnls) { // solve least squares
122  IMP_Eigen::JacobiSVD<IMP_Eigen::MatrixXf> svd(W_.asDiagonal() * A_,
123  ComputeThinU | ComputeThinV);
124  w = svd.solve(Wb_);
125  // zero the negatives
126  for (int i = 0; i < w.size(); i++)
127  if (w(i) < 0) w(i) = 0;
128  } else {
129  w = NNLS(W_.asDiagonal() * A_, Wb_);
130  }
131  w /= w.sum();
132 
133  // compute weighted profile
134  IMP_NEW(Profile, weighted_profile,
138  IMP_Eigen::VectorXf wp = A_ * w;
139  weighted_profile->set_qs(profiles[0]->get_qs());
140  weighted_profile->set_intensities(wp);
141  weights.resize(w.size());
142  for (int i = 0; i < w.size(); i++) weights[i] = w[i];
143 
146  weighted_profile, use_offset);
147 }
148 
149 template <typename ScoringFunctionT>
151  ProfilesTemp partial_profiles,
152  double min_c1, double max_c1,
153  double min_c2, double max_c2,
154  bool use_offset) const {
155  Vector<double> weights;
157  search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2, max_c2,
158  std::numeric_limits<double>::max(), weights, use_offset);
159  return fp;
160 }
161 
162 template <typename ScoringFunctionT>
164  ProfilesTemp partial_profiles,
165  const WeightedFitParameters& fp,
166  const std::string fit_file_name,
167  bool use_offset) const {
168  double best_c1 = fp.get_c1();
169  double best_c2 = fp.get_c2();
170 
171  // compute a profile for best c1/c2 combination
172  for (unsigned int i = 0; i < partial_profiles.size(); i++)
173  partial_profiles[i]->sum_partial_profiles(best_c1, best_c2);
174 
175  if (partial_profiles.size() == 1) {
176  // compute scale
177  double offset = 0.0;
178  if (use_offset)
179  offset =
182  partial_profiles[0]);
183 
184  double c =
187  partial_profiles[0], offset);
189  fit_file_name, partial_profiles[0], fp.get_chi(), c, offset);
190  } else {
191 
192  // computed weighted profile
193  IMP_NEW(Profile, weighted_profile,
197 
198  const Vector<double>& weights = fp.get_weights();
199  for (unsigned int i = 0; i < weights.size(); i++)
200  weighted_profile->add(partial_profiles[i], weights[i]);
201 
202  // compute scale
203  double offset = 0.0;
204  if (use_offset)
205  offset =
208  weighted_profile);
209  double c =
212  weighted_profile, offset);
213 
215  fit_file_name, weighted_profile,
216  fp.get_chi(), c, offset);
217  }
218 }
219 
220 template <typename ScoringFunctionT>
222  ProfilesTemp& partial_profiles,
223  double min_c1, double max_c1,
224  double min_c2, double max_c2,
225  double old_chi, Vector<double>& weights,
226  bool use_offset) const {
227  int c1_cells = 10;
228  int c2_cells = 10;
229  if (old_chi < (std::numeric_limits<double>::max() - 1)) { // second iteration
230  c1_cells = 5;
231  c2_cells = 5;
232  }
233 
234  double delta_c1 = (max_c1 - min_c1) / c1_cells;
235  double delta_c2 = (max_c2 - min_c2) / c2_cells;
236 
237  bool last_c1 = false;
238  bool last_c2 = false;
239  if (delta_c1 < 0.0001) {
240  c1_cells = 1;
241  delta_c1 = max_c1 - min_c1;
242  last_c1 = true;
243  }
244  if (delta_c2 < 0.001) {
245  c2_cells = 1;
246  delta_c2 = max_c2 - min_c2;
247  last_c2 = true;
248  }
249  double best_c1(1.0), best_c2(0.0), best_chi(old_chi);
250  bool best_set = false;
251 
252  // c1 iteration
253  double c1(min_c1);
254  for (int i = 0; i <= c1_cells; i++, c1 += delta_c1) {
255  // c2 iteration
256  double c2 = min_c2;
257  for (int j = 0; j <= c2_cells; j++, c2 += delta_c2) {
258  // sum up the profiles for c1/c2 combo
259  for (unsigned int k = 0; k < partial_profiles.size(); k++)
260  partial_profiles[k]->sum_partial_profiles(c1, c2);
261  Vector<double> curr_weights;
262  double curr_chi = compute_score(partial_profiles, curr_weights, use_offset);
263  if (!best_set || curr_chi < best_chi) {
264  best_set = true;
265  best_chi = curr_chi;
266  best_c1 = c1;
267  best_c2 = c2;
268  weights = curr_weights;
269  }
270  }
271  }
272 
273  if (std::fabs(best_chi - old_chi) > 0.005 &&
274  (!(last_c1 && last_c2))) { // refine more
275  min_c1 = std::max(best_c1 - delta_c1, min_c1);
276  max_c1 = std::min(best_c1 + delta_c1, max_c1);
277  min_c2 = std::max(best_c2 - delta_c2, min_c2);
278  max_c2 = std::min(best_c2 + delta_c2, max_c2);
279  return search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2,
280  max_c2, best_chi, weights, use_offset);
281  }
282  return WeightedFitParameters(best_chi, best_c1, best_c2, weights);
283 }
284 
285 IMPSAXS_END_NAMESPACE
286 
287 #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
A more IMP-like version of the std::vector.
Definition: Vector.h:39
#define IMP_NEW(Typename, varname, args)
Declare a ref counted pointer to a new object.
Definition: object_macros.h:64
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
double compute_offset(const Profile *model_profile) const
computes offset for fitting for fitting the modeled profile
Definition: ProfileFitter.h:88
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