IMP logo
IMP Reference Guide  2.16.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-2021 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 <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  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  Eigen::MatrixXf W_; // weights matrix
82 
83  // weights matrix multiplied by experimental intensities vector
84  Eigen::VectorXf Wb_;
85 
86  // intensities
87  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  Eigen::VectorXf w;
121  if (!nnls) { // solve least squares
122  Eigen::JacobiSVD<Eigen::MatrixXf> svd(W_.asDiagonal() * A_,
123  Eigen::ComputeThinU
124  | Eigen::ComputeThinV);
125  w = svd.solve(Wb_);
126  // zero the negatives
127  for (int i = 0; i < w.size(); i++)
128  if (w(i) < 0) w(i) = 0;
129  } else {
130  w = NNLS(W_.asDiagonal() * A_, Wb_);
131  }
132  w /= w.sum();
133 
134  // compute weighted profile
135  IMP_NEW(Profile, weighted_profile,
139  Eigen::VectorXf wp = A_ * w;
140  weighted_profile->set_qs(profiles[0]->get_qs());
141  weighted_profile->set_intensities(wp);
142  weights.resize(w.size());
143  for (int i = 0; i < w.size(); i++) weights[i] = w[i];
144 
147  weighted_profile, use_offset);
148 }
149 
150 template <typename ScoringFunctionT>
152  ProfilesTemp partial_profiles,
153  double min_c1, double max_c1,
154  double min_c2, double max_c2,
155  bool use_offset) const {
156  Vector<double> weights;
158  search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2, max_c2,
159  std::numeric_limits<double>::max(), weights, use_offset);
160  return fp;
161 }
162 
163 template <typename ScoringFunctionT>
165  ProfilesTemp partial_profiles,
166  const WeightedFitParameters& fp,
167  const std::string fit_file_name,
168  bool use_offset) const {
169  double best_c1 = fp.get_c1();
170  double best_c2 = fp.get_c2();
171 
172  // compute a profile for best c1/c2 combination
173  for (unsigned int i = 0; i < partial_profiles.size(); i++)
174  partial_profiles[i]->sum_partial_profiles(best_c1, best_c2);
175 
176  if (partial_profiles.size() == 1) {
177  // compute scale
178  double offset = 0.0;
179  if (use_offset)
180  offset =
183  partial_profiles[0]);
184 
185  double c =
188  partial_profiles[0], offset);
190  fit_file_name, partial_profiles[0], fp.get_chi_square(), c, offset);
191  } else {
192 
193  // computed weighted profile
194  IMP_NEW(Profile, weighted_profile,
198 
199  const Vector<double>& weights = fp.get_weights();
200  for (unsigned int i = 0; i < weights.size(); i++)
201  weighted_profile->add(partial_profiles[i], weights[i]);
202 
203  // compute scale
204  double offset = 0.0;
205  if (use_offset)
206  offset =
209  weighted_profile);
210  double c =
213  weighted_profile, offset);
214 
216  fit_file_name, weighted_profile,
217  fp.get_chi_square(), c, offset);
218  }
219 }
220 
221 template <typename ScoringFunctionT>
223  ProfilesTemp& partial_profiles,
224  double min_c1, double max_c1,
225  double min_c2, double max_c2,
226  double old_chi, Vector<double>& weights,
227  bool use_offset) const {
228  int c1_cells = 10;
229  int c2_cells = 10;
230  if (old_chi < (std::numeric_limits<double>::max() - 1)) { // second iteration
231  c1_cells = 5;
232  c2_cells = 5;
233  }
234 
235  double delta_c1 = (max_c1 - min_c1) / c1_cells;
236  double delta_c2 = (max_c2 - min_c2) / c2_cells;
237 
238  bool last_c1 = false;
239  bool last_c2 = false;
240  if (delta_c1 < 0.0001) {
241  c1_cells = 1;
242  delta_c1 = max_c1 - min_c1;
243  last_c1 = true;
244  }
245  if (delta_c2 < 0.001) {
246  c2_cells = 1;
247  delta_c2 = max_c2 - min_c2;
248  last_c2 = true;
249  }
250  double best_c1(1.0), best_c2(0.0), best_chi(old_chi);
251  bool best_set = false;
252 
253  // c1 iteration
254  double c1(min_c1);
255  for (int i = 0; i <= c1_cells; i++, c1 += delta_c1) {
256  // c2 iteration
257  double c2 = min_c2;
258  for (int j = 0; j <= c2_cells; j++, c2 += delta_c2) {
259  // sum up the profiles for c1/c2 combo
260  for (unsigned int k = 0; k < partial_profiles.size(); k++)
261  partial_profiles[k]->sum_partial_profiles(c1, c2);
262  Vector<double> curr_weights;
263  double curr_chi = compute_score(partial_profiles, curr_weights, use_offset);
264  if (!best_set || curr_chi < best_chi) {
265  best_set = true;
266  best_chi = curr_chi;
267  best_c1 = c1;
268  best_c2 = c2;
269  weights = curr_weights;
270  }
271  }
272  }
273 
274  if (std::fabs(best_chi - old_chi) > 0.005 &&
275  (!(last_c1 && last_c2))) { // refine more
276  min_c1 = std::max(best_c1 - delta_c1, min_c1);
277  max_c1 = std::min(best_c1 + delta_c1, max_c1);
278  min_c2 = std::max(best_c2 - delta_c2, min_c2);
279  max_c2 = std::min(best_c2 + delta_c2, max_c2);
280  return search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2,
281  max_c2, best_chi, weights, use_offset);
282  }
283  return WeightedFitParameters(best_chi, best_c1, best_c2, weights);
284 }
285 
286 IMPSAXS_END_NAMESPACE
287 
288 #endif /* IMPSAXS_WEIGHTED_PROFILE_FITTER_H */
unsigned int size() const
return number of entries in SAXS profile
Definition: Profile.h:178
Basic chi score implementation.
WeightedProfileFitter(const Profile *exp_profile)
Constructor.
Eigen::VectorXf NNLS(const Eigen::MatrixXf &A, const Eigen::VectorXf &b)
non-negative least square fitting for profile weight solving problem
Definition: nnls.h:19
a class for fitting two profiles
A more IMP-like version of the std::vector.
Definition: Vector.h:40
#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.
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