8 #ifndef IMPSAXS_PROFILE_FITTER_H
9 #define IMPSAXS_PROFILE_FITTER_H
18 IMPSAXS_BEGIN_NAMESPACE
28 template <
typename ScoringFunctionT = ChiScore>
36 :
Object(
"ProfileFitter%1%"), exp_profile_(exp_profile) {
38 scoring_function_ =
new ScoringFunctionT();
42 :
Object(
"ProfileFitter%1%"), exp_profile_(exp_profile) {
44 scoring_function_ = sf;
48 double compute_score(
const Profile* model_profile,
bool use_offset =
false,
49 const std::string fit_file_name =
"")
const;
65 FitParameters fit_profile(Profile* partial_profile,
double min_c1 = 0.95,
66 double max_c1 = 1.05,
double min_c2 = -2.0,
67 double max_c2 = 4.0,
bool use_offset =
false,
68 const std::string fit_file_name =
"")
const;
77 double offset = 0.0)
const {
78 return scoring_function_->compute_scale_factor(exp_profile_, model_profile,
89 return scoring_function_->compute_offset(exp_profile_, model_profile);
94 void resample(
const Profile* model_profile,
Profile* resampled_profile)
const;
97 void write_SAXS_fit_file(
const std::string& file_name,
98 const Profile* model_profile,
const double chi_square,
99 const double c = 1,
const double offset = 0)
const;
103 double max_c1,
double min_c2,
double max_c2,
104 bool use_offset,
double old_chi)
const;
115 template <
typename ScoringFunctionT>
118 model_profile->
resample(exp_profile_, resampled_profile);
121 template <
typename ScoringFunctionT>
123 Profile* partial_profile,
double min_c1,
double max_c1,
double min_c2,
124 double max_c2,
bool use_offset,
double old_chi)
const {
127 if (old_chi < (std::numeric_limits<double>::max() - 1)) {
132 double delta_c1 = (max_c1 - min_c1) / c1_cells;
133 double delta_c2 = (max_c2 - min_c2) / c2_cells;
135 bool last_c1 =
false;
136 bool last_c2 =
false;
137 if (delta_c1 < 0.0001) {
139 delta_c1 = max_c1 - min_c1;
142 if (delta_c2 < 0.001) {
144 delta_c2 = max_c2 - min_c2;
147 double best_c1(1.0), best_c2(0.0), best_chi(old_chi);
148 bool best_set =
false;
151 for (
int i = 0; i <= c1_cells; i++, c1 += delta_c1) {
153 for (
int j = 0; j <= c2_cells; j++, c2 += delta_c2) {
155 double curr_chi = compute_score(partial_profile, use_offset);
156 if (!best_set || curr_chi < best_chi) {
165 if (std::fabs(best_chi - old_chi) > 0.0001 &&
166 (!(last_c1 && last_c2))) {
167 min_c1 = std::max(best_c1 - delta_c1, min_c1);
168 max_c1 = std::min(best_c1 + delta_c1, max_c1);
169 min_c2 = std::max(best_c2 - delta_c2, min_c2);
170 max_c2 = std::min(best_c2 + delta_c2, max_c2);
171 return search_fit_parameters(partial_profile, min_c1, max_c1, min_c2,
172 max_c2, use_offset, best_chi);
174 return FitParameters(best_chi, best_c1, best_c2);
177 template <
typename ScoringFunctionT>
179 Profile* partial_profile,
double min_c1,
double max_c1,
double min_c2,
180 double max_c2,
bool use_offset,
const std::string fit_file_name)
const {
183 double default_c1 = 1.0, default_c2 = 0.0;
185 double default_chi = compute_score(partial_profile, use_offset);
188 search_fit_parameters(partial_profile, min_c1, max_c1, min_c2, max_c2,
189 use_offset, std::numeric_limits<double>::max());
190 double best_c1 = fp.get_c1();
191 double best_c2 = fp.get_c2();
192 fp.set_default_chi_square(default_chi);
196 compute_score(partial_profile, use_offset, fit_file_name);
201 template <
typename ScoringFunctionT>
203 const Profile* model_profile,
bool use_offset,
204 const std::string fit_file_name)
const {
206 (exp_profile_->get_min_q(), exp_profile_->get_max_q(),
207 exp_profile_->get_delta_q()));
208 model_profile->
resample(exp_profile_, resampled_profile);
210 double score = scoring_function_->compute_score(exp_profile_,
213 if (fit_file_name.length() > 0) {
217 scoring_function_->compute_offset(exp_profile_, resampled_profile);
218 double c = scoring_function_->compute_scale_factor(exp_profile_,
221 write_SAXS_fit_file(fit_file_name, resampled_profile, score, c, offset);
226 template <
typename ScoringFunctionT>
228 const std::string& file_name,
const Profile* model_profile,
229 const double score,
const double c,
const double offset)
const {
230 std::ofstream out_file(file_name.c_str());
235 unsigned int profile_size =
236 std::min(model_profile->
size(), exp_profile_->size());
238 out_file.precision(15);
239 out_file <<
"# SAXS profile: number of points = " << profile_size
240 <<
", q_min = " << exp_profile_->get_min_q()
241 <<
", q_max = " << exp_profile_->get_max_q();
242 out_file <<
", delta_q = " << exp_profile_->get_delta_q() << std::endl;
244 out_file.setf(std::ios::showpoint);
245 out_file <<
"# offset = " << offset <<
", scaling c = " << c
246 <<
", Chi^2 = " << score << std::endl;
247 out_file <<
"# q exp_intensity model_intensity error" << std::endl;
249 out_file.setf(std::ios::fixed, std::ios::floatfield);
251 for (
unsigned int i = 0; i < profile_size; i++) {
252 out_file.setf(std::ios::left);
254 out_file.precision(8);
255 out_file << exp_profile_->get_q(i) <<
" ";
257 out_file.setf(std::ios::left);
259 out_file.precision(8);
260 out_file << exp_profile_->get_intensity(i) <<
" ";
262 out_file.setf(std::ios::left);
264 out_file.precision(8);
265 out_file << model_profile->get_intensity(i) * c - offset <<
" ";
267 out_file.setf(std::ios::left);
269 out_file.precision(8);
270 out_file << exp_profile_->get_error(i) << std::endl;
276 std::string file_name2 = file_name.substr(0, file_name.length()-4);
277 file_name2 +=
".fit";
278 std::ofstream out_file2(file_name2.c_str());
280 IMP_THROW(
"Can't open file " << file_name2, IOException);
284 out_file2.precision(15);
285 out_file2 <<
"# SAXS profile: number of points = " << profile_size
286 <<
", q_min = " << exp_profile_->get_min_q()
287 <<
", q_max = " << exp_profile_->get_max_q();
288 out_file2 <<
", delta_q = " << exp_profile_->get_delta_q() << std::endl;
290 out_file2.setf(std::ios::showpoint);
291 out_file2 <<
"# offset = " << offset <<
", scaling c = " << c
292 <<
", Chi^2 = " << score << std::endl;
293 out_file2 <<
"# q exp_intensity error model_intensity" << std::endl;
295 out_file2.setf(std::ios::fixed, std::ios::floatfield);
297 for (
unsigned int i = 0; i < profile_size; i++) {
298 out_file2.setf(std::ios::left);
300 out_file2.precision(8);
301 out_file2 << exp_profile_->get_q(i) <<
" ";
303 out_file2.setf(std::ios::left);
305 out_file2.precision(8);
306 out_file2 << exp_profile_->get_intensity(i) <<
" ";
308 out_file2.setf(std::ios::left);
310 out_file2.precision(8);
311 out_file2 << exp_profile_->get_error(i) <<
" ";
313 out_file2.setf(std::ios::left);
315 out_file2.precision(8);
316 out_file2 << model_profile->get_intensity(i) * c - offset << std::endl;
323 IMPSAXS_END_NAMESPACE
unsigned int size() const
return number of entries in SAXS profile
ProfileFitter(const Profile *exp_profile)
Constructor.
Basic chi score implementation.
Parameters of a fit, from ProfileFitter.
#define IMP_REF_COUNTED_DESTRUCTOR(Name)
Ref counted objects should have private destructors.
#define IMP_NEW(Typename, varname, args)
Declare a ref counted pointer to a new object.
An input/output exception.
Common base class for heavy weight IMP objects.
A smart pointer to a ref-counted Object that is a class member.
Copyright 2007-2019 IMP Inventors. All rights reserved.
double compute_offset(const Profile *model_profile) const
computes offset for fitting for fitting the modeled profile
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
A shared base class to help in debugging and things.
void resample(const Profile *exp_profile, Profile *resampled_profile) const
return a profile that is sampled on the q values of the exp_profile
Fit two profiles with user-defined scoring function as a template parameter.
void sum_partial_profiles(double c1, double c2, bool check_cashed=true)
computes full profile for given fitting parameters
double compute_scale_factor(const Profile *model_profile, double offset=0.0) const
computes the scaling factor needed for fitting the modeled profile
A class for profile storing and computation.