12 #ifndef IMPSAXS_WEIGHTED_PROFILE_FITTER_H
13 #define IMPSAXS_WEIGHTED_PROFILE_FITTER_H
19 #include <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 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;
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);
122 Eigen::JacobiSVD<Eigen::MatrixXf> svd(W_.asDiagonal() * A_,
124 | Eigen::ComputeThinV);
127 for (
int i = 0; i < w.size(); i++)
128 if (w(i) < 0) w(i) = 0;
130 w =
NNLS(W_.asDiagonal() * A_, Wb_);
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];
147 weighted_profile, use_offset);
150 template <
typename ScoringFunctionT>
153 double min_c1,
double max_c1,
154 double min_c2,
double max_c2,
155 bool use_offset)
const {
158 search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2, max_c2,
159 std::numeric_limits<double>::max(), weights, use_offset);
163 template <
typename ScoringFunctionT>
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();
173 for (
unsigned int i = 0; i < partial_profiles.size(); i++)
174 partial_profiles[i]->sum_partial_profiles(best_c1, best_c2);
176 if (partial_profiles.size() == 1) {
183 partial_profiles[0]);
188 partial_profiles[0], offset);
190 fit_file_name, partial_profiles[0], fp.get_chi_square(), c, offset);
200 for (
unsigned int i = 0; i < weights.size(); i++)
201 weighted_profile->add(partial_profiles[i], weights[i]);
213 weighted_profile, offset);
216 fit_file_name, weighted_profile,
217 fp.get_chi_square(), c, offset);
221 template <
typename ScoringFunctionT>
224 double min_c1,
double max_c1,
225 double min_c2,
double max_c2,
227 bool use_offset)
const {
230 if (old_chi < (std::numeric_limits<double>::max() - 1)) {
235 double delta_c1 = (max_c1 - min_c1) / c1_cells;
236 double delta_c2 = (max_c2 - min_c2) / c2_cells;
238 bool last_c1 =
false;
239 bool last_c2 =
false;
240 if (delta_c1 < 0.0001) {
242 delta_c1 = max_c1 - min_c1;
245 if (delta_c2 < 0.001) {
247 delta_c2 = max_c2 - min_c2;
250 double best_c1(1.0), best_c2(0.0), best_chi(old_chi);
251 bool best_set =
false;
255 for (
int i = 0; i <= c1_cells; i++, c1 += delta_c1) {
258 for (
int j = 0; j <= c2_cells; j++, c2 += delta_c2) {
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) {
269 weights = curr_weights;
274 if (std::fabs(best_chi - old_chi) > 0.005 &&
275 (!(last_c1 && last_c2))) {
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);
283 return WeightedFitParameters(best_chi, best_c1, best_c2, weights);
286 IMPSAXS_END_NAMESPACE
unsigned int size() const
return number of entries in SAXS profile
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
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.
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