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_;
65 double min_c1 = 0.95,
double max_c1 = 1.05,
67 double max_c2 = 4.0)
const;
72 const std::string fit_file_name)
const;
76 ProfilesTemp& partial_profiles,
double min_c1,
double max_c1,
double min_c2,
80 IMP_Eigen::MatrixXf W_;
83 IMP_Eigen::VectorXf Wb_;
86 IMP_Eigen::MatrixXf A_;
89 template <
typename ScoringFunctionT>
96 if (profiles.size() == 1) {
105 int m = profiles.size();
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);
119 IMP_Eigen::VectorXf w;
121 IMP_Eigen::JacobiSVD<IMP_Eigen::MatrixXf> svd(W_.asDiagonal() * A_,
122 ComputeThinU | ComputeThinV);
125 for (
int i = 0; i < w.size(); i++)
126 if (w(i) < 0) w(i) = 0;
128 w =
NNLS(W_.asDiagonal() * A_, Wb_);
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];
148 template <
typename ScoringFunctionT>
151 double min_c1,
double max_c1,
152 double min_c2,
double max_c2)
const {
155 search_fit_parameters(partial_profiles, min_c1, max_c1, min_c2, max_c2,
156 std::numeric_limits<double>::max(), weights);
160 template <
typename ScoringFunctionT>
164 const std::string fit_file_name)
const {
165 double best_c1 = fp.get_c1();
166 double best_c2 = fp.get_c2();
169 for (
unsigned int i = 0; i < partial_profiles.size(); i++)
170 partial_profiles[i]->sum_partial_profiles(best_c1, best_c2);
172 if (partial_profiles.size() == 1) {
177 partial_profiles[0]);
179 fit_file_name, partial_profiles[0], fp.get_chi(), c);
189 for (
unsigned int i = 0; i < weights.size(); i++)
190 weighted_profile->add(partial_profiles[i], weights[i]);
199 fit_file_name, weighted_profile, fp.get_chi(), c);
203 template <
typename ScoringFunctionT>
206 double min_c1,
double max_c1,
207 double min_c2,
double max_c2,
211 if (old_chi < (std::numeric_limits<double>::max() - 1)) {
216 double delta_c1 = (max_c1 - min_c1) / c1_cells;
217 double delta_c2 = (max_c2 - min_c2) / c2_cells;
219 bool last_c1 =
false;
220 bool last_c2 =
false;
221 if (delta_c1 < 0.0001) {
223 delta_c1 = max_c1 - min_c1;
226 if (delta_c2 < 0.001) {
228 delta_c2 = max_c2 - min_c2;
231 double best_c1(1.0), best_c2(0.0), best_chi(old_chi);
232 bool best_set =
false;
236 for (
int i = 0; i <= c1_cells; i++, c1 += delta_c1) {
239 for (
int j = 0; j <= c2_cells; j++, c2 += delta_c2) {
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) {
250 weights = curr_weights;
255 if (std::fabs(best_chi - old_chi) > 0.005 &&
256 (!(last_c1 && last_c2))) {
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);
264 return WeightedFitParameters(best_chi, best_c1, best_c2, weights);
267 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
#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
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