8 #ifndef IMPSAXS_PROFILE_FITTER_H
9 #define IMPSAXS_PROFILE_FITTER_H
17 IMPSAXS_BEGIN_NAMESPACE
28 template<
class ScoringFunctionT = ChiScore>
29 class ProfileFitter:
public base::RefCounted {
35 ProfileFitter(
const Profile& exp_profile): exp_profile_(exp_profile) {
36 scoring_function_ =
new ScoringFunctionT();
40 Float compute_score(
const Profile& model_profile,
41 bool use_offset =
false,
42 const std::string fit_file_name =
"")
const;
45 Float compute_score(
const Profile& model_profile,
56 FitParameters fit_profile(Profile& partial_profile,
57 float min_c1=0.95,
float max_c1=1.05,
58 float min_c2=-2.0,
float max_c2=4.0,
59 bool use_offset =
false,
60 const std::string fit_file_name =
"")
const;
68 Float compute_scale_factor(
const Profile& model_profile,
69 Float offset = 0.0)
const {
70 return scoring_function_->compute_scale_factor(exp_profile_, model_profile,
80 Float compute_offset(
const Profile& model_profile)
const {
81 return scoring_function_->compute_offset(exp_profile_, model_profile);
86 void resample(
const Profile& model_profile, Profile& resampled_profile)
const;
90 FitParameters search_fit_parameters(Profile& partial_profile,
91 float min_c1,
float max_c1,
92 float min_c2,
float max_c2,
93 bool use_offset,
float old_chi)
const;
96 void write_SAXS_fit_file(
const std::string& file_name,
97 const Profile& model_profile,
98 const Float chi_square,
102 friend class DerivativeCalculator;
105 const Profile exp_profile_;
106 ScoringFunctionT* scoring_function_;
109 template<
class ScoringFunctionT>
110 void ProfileFitter<ScoringFunctionT>::resample(
const Profile& model_profile,
111 Profile& resampled_profile)
const
114 std::map<float, unsigned int> q_mapping;
115 for (
unsigned int k=0; k<model_profile.size(); k++) {
116 q_mapping[model_profile.get_q(k)] = k;
119 for (
unsigned int k=0; k<exp_profile_.size(); k++) {
120 Float q = exp_profile_.get_q(k);
121 if(q>model_profile.get_max_q())
break;
122 std::map<float, unsigned int>::iterator it = q_mapping.lower_bound(q);
123 if(it == q_mapping.end())
break;
124 unsigned int i = it->second;
126 resampled_profile.add_entry(q, model_profile.get_intensity(i));
128 Float delta_q = model_profile.get_q(i)-model_profile.get_q(i-1);
129 if(delta_q <= 1.0e-16) {
130 resampled_profile.add_entry(q, model_profile.get_intensity(i));
132 Float alpha = (q - model_profile.get_q(i-1)) / delta_q;
133 if(alpha > 1.0) alpha = 1.0;
134 Float intensity = model_profile.get_intensity(i-1)
135 + (alpha)*(model_profile.get_intensity(i)
136 - model_profile.get_intensity(i-1));
137 resampled_profile.add_entry(q, intensity);
143 template<
class ScoringFunctionT>
144 FitParameters ProfileFitter<ScoringFunctionT>::search_fit_parameters(
145 Profile& partial_profile,
146 float min_c1,
float max_c1,
147 float min_c2,
float max_c2,
148 bool use_offset,
float old_chi)
const
152 if(old_chi < (std::numeric_limits<float>::max()-1)) {
157 float delta_c1 = (max_c1-min_c1)/c1_cells;
158 float delta_c2 = (max_c2-min_c2)/c2_cells;
160 bool last_c1 =
false;
161 bool last_c2 =
false;
162 if(delta_c1 < 0.0001) { c1_cells=1; delta_c1 = max_c1-min_c1; last_c1=
true; }
163 if(delta_c2 < 0.001) { c2_cells=1; delta_c2 = max_c2-min_c2; last_c2=
true; }
164 float best_c1(1.0), best_c2(0.0), best_chi(old_chi);
165 bool best_set =
false;
168 for(
int i=0; i<=c1_cells; i++, c1+= delta_c1) {
170 for(
int j=0; j<=c2_cells; j++, c2+= delta_c2) {
171 partial_profile.sum_partial_profiles(c1, c2, partial_profile);
172 float curr_chi = compute_score(partial_profile, use_offset);
173 if(!best_set || curr_chi < best_chi) {
182 if(std::fabs(best_chi-old_chi) > 0.0001 &&
183 (!(last_c1 && last_c2))) {
184 min_c1 = std::max(best_c1-delta_c1, min_c1);
185 max_c1 = std::min(best_c1+delta_c1, max_c1);
186 min_c2 = std::max(best_c2-delta_c2, min_c2);
187 max_c2 = std::min(best_c2+delta_c2, max_c2);
188 return search_fit_parameters(partial_profile,
189 min_c1, max_c1, min_c2, max_c2,
190 use_offset, best_chi);
192 return FitParameters(best_chi, best_c1, best_c2);
195 template<
class ScoringFunctionT>
196 FitParameters ProfileFitter<ScoringFunctionT>::fit_profile(
197 Profile& partial_profile,
198 float min_c1,
float max_c1,
199 float min_c2,
float max_c2,
201 const std::string fit_file_name)
const {
204 float default_c1 = 1.0, default_c2 = 0.0;
205 partial_profile.sum_partial_profiles(default_c1, default_c2, partial_profile);
206 float default_chi = compute_score(partial_profile, use_offset);
208 FitParameters fp = search_fit_parameters(partial_profile,
209 min_c1, max_c1, min_c2, max_c2,
211 std::numeric_limits<float>::max());
212 float best_c1 = fp.get_c1();
213 float best_c2 = fp.get_c2();
214 fp.set_default_chi(default_chi);
217 partial_profile.sum_partial_profiles(best_c1, best_c2, partial_profile);
218 compute_score(partial_profile, use_offset, fit_file_name);
225 template<
class ScoringFunctionT>
226 Float ProfileFitter<ScoringFunctionT>::compute_score(
227 const Profile& model_profile,
230 Profile resampled_profile(exp_profile_.get_min_q(),
231 exp_profile_.get_max_q(),
232 exp_profile_.get_delta_q());
233 resample(model_profile, resampled_profile);
235 Float score = scoring_function_->compute_score(exp_profile_,
241 template<
class ScoringFunctionT>
242 Float ProfileFitter<ScoringFunctionT>::compute_score(
243 const Profile& model_profile,
245 const std::string fit_file_name)
const
247 Profile resampled_profile(exp_profile_.get_min_q(),
248 exp_profile_.get_max_q(),
249 exp_profile_.get_delta_q());
250 resample(model_profile, resampled_profile);
252 Float score = scoring_function_->compute_score(exp_profile_,
253 resampled_profile, use_offset);
255 if(fit_file_name.length() > 0) {
258 offset=scoring_function_->compute_offset(exp_profile_, resampled_profile);
259 Float c= scoring_function_->compute_scale_factor(exp_profile_,
260 resampled_profile, offset);
261 write_SAXS_fit_file(fit_file_name, resampled_profile, score, c, offset);
266 template<
class ScoringFunctionT>
267 void ProfileFitter<ScoringFunctionT>::write_SAXS_fit_file(
268 const std::string& file_name,
269 const Profile& model_profile,
272 const Float offset)
const {
273 std::ofstream out_file(file_name.c_str());
275 IMP_THROW(
"Can't open file " << file_name,
279 unsigned int profile_size = std::min(model_profile.size(),
280 exp_profile_.size());
282 out_file.precision(15);
283 out_file <<
"# SAXS profile: number of points = " << profile_size
284 <<
", q_min = " << exp_profile_.get_min_q()
285 <<
", q_max = " << exp_profile_.get_max_q();
286 out_file <<
", delta_q = " << exp_profile_.get_delta_q() << std::endl;
288 out_file.setf(std::ios::showpoint);
289 out_file <<
"# offset = " << offset <<
", scaling c = " << c
290 <<
", Chi = " << score << std::endl;
291 out_file <<
"# q exp_intensity model_intensity"
294 out_file.setf(std::ios::fixed, std::ios::floatfield);
296 for (
unsigned int i = 0; i < profile_size; i++) {
297 out_file.setf(std::ios::left);
299 out_file.precision(5);
300 out_file << exp_profile_.get_q(i) <<
" ";
302 out_file.setf(std::ios::left);
304 out_file.precision(8);
305 out_file << exp_profile_.get_intensity(i) <<
" ";
307 out_file.setf(std::ios::left);
309 out_file.precision(8);
310 out_file << model_profile.get_intensity(i)*c - offset << std::endl;
315 IMPSAXS_END_NAMESPACE