12 #ifndef IMPSAXS_WEIGHTED_PROFILE_FITTER_H
13 #define IMPSAXS_WEIGHTED_PROFILE_FITTER_H
19 #include <IMP/algebra/eigen3/Eigen/Dense>
21 IMPSAXS_BEGIN_NAMESPACE
27 template <
typename ScoringFunctionT = ChiScore>
37 W_(exp_profile->size(), 1),
38 Wb_(exp_profile->size()),
39 A_(exp_profile->size(), 2) {
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));
46 Wb_ = W_.asDiagonal() * Wb_;
57 bool use_offset =
false,
bool NNLS =
true)
const;
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;
72 const std::string fit_file_name,
73 bool use_offset =
false)
const;
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;
81 IMP_Eigen::MatrixXf W_;
84 IMP_Eigen::VectorXf Wb_;
87 IMP_Eigen::MatrixXf A_;
90 template <
typename ScoringFunctionT>
94 bool use_offset,
bool nnls)
const {
97 if (profiles.size() == 1) {
102 profiles[0], use_offset);
106 int m = profiles.size();
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);
120 IMP_Eigen::VectorXf w;
122 IMP_Eigen::JacobiSVD<IMP_Eigen::MatrixXf> svd(W_.asDiagonal() * A_,
123 ComputeThinU | ComputeThinV);
126 for (
int i = 0; i < w.size(); i++)
127 if (w(i) < 0) w(i) = 0;
129 w =
NNLS(W_.asDiagonal() * A_, Wb_);
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];
146 weighted_profile, use_offset);
149 template <
typename ScoringFunctionT>
152 double min_c1,
double max_c1,
153 double min_c2,
double max_c2,
154 bool use_offset)
const {
157 search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2, max_c2,
158 std::numeric_limits<double>::max(), weights, use_offset);
162 template <
typename ScoringFunctionT>
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();
172 for (
unsigned int i = 0; i < partial_profiles.size(); i++)
173 partial_profiles[i]->sum_partial_profiles(best_c1, best_c2);
175 if (partial_profiles.size() == 1) {
182 partial_profiles[0]);
187 partial_profiles[0], offset);
189 fit_file_name, partial_profiles[0], fp.get_chi(), c, offset);
199 for (
unsigned int i = 0; i < weights.size(); i++)
200 weighted_profile->add(partial_profiles[i], weights[i]);
212 weighted_profile, offset);
215 fit_file_name, weighted_profile,
216 fp.get_chi(), c, offset);
220 template <
typename ScoringFunctionT>
223 double min_c1,
double max_c1,
224 double min_c2,
double max_c2,
226 bool use_offset)
const {
229 if (old_chi < (std::numeric_limits<double>::max() - 1)) {
234 double delta_c1 = (max_c1 - min_c1) / c1_cells;
235 double delta_c2 = (max_c2 - min_c2) / c2_cells;
237 bool last_c1 =
false;
238 bool last_c2 =
false;
239 if (delta_c1 < 0.0001) {
241 delta_c1 = max_c1 - min_c1;
244 if (delta_c2 < 0.001) {
246 delta_c2 = max_c2 - min_c2;
249 double best_c1(1.0), best_c2(0.0), best_chi(old_chi);
250 bool best_set =
false;
254 for (
int i = 0; i <= c1_cells; i++, c1 += delta_c1) {
257 for (
int j = 0; j <= c2_cells; j++, c2 += delta_c2) {
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) {
268 weights = curr_weights;
273 if (std::fabs(best_chi - old_chi) > 0.005 &&
274 (!(last_c1 && last_c2))) {
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);
282 return WeightedFitParameters(best_chi, best_c1, best_c2, weights);
285 IMPSAXS_END_NAMESPACE
unsigned int size() const
return number of entries in SAXS profile
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.
#define IMP_NEW(Typename, varname, args)
Declare a ref counted pointer to a new object.
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
double compute_offset(const Profile *model_profile) const
computes offset for fitting for fitting the modeled profile
Parameters of a weighted fit, from WeightedProfileFitter.
Fit two profiles with user-defined scoring function as a template parameter.
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