diff --git a/src/fitting/models/gaussian_model.cpp b/src/fitting/models/gaussian_model.cpp index f1b59449..dbd9920a 100644 --- a/src/fitting/models/gaussian_model.cpp +++ b/src/fitting/models/gaussian_model.cpp @@ -109,7 +109,7 @@ Fit_Parameters Gaussian_Model::_generate_default_fit_parameters( fit_params.add_parameter(Fit_Param(STR_COMPTON_FWHM_CORR, (T_real)1.0, (T_real)4.0, (T_real)1.0, (T_real)0.1, E_Bound_Type::LIMITED_LO_HI)); fit_params.add_parameter(Fit_Param(STR_COMPTON_AMPLITUDE, MIN_COUNTS_LIMIT_LOG, MAX_COUNTS_LIMIT_LOG, (T_real)5.0, STEP_COUNTS_LIMIT_LOG, E_Bound_Type::LIMITED_LO_HI)); fit_params.add_parameter(Fit_Param(STR_COMPTON_F_STEP, (T_real)0.0, (T_real)1.0, (T_real)0.0, (T_real)0.1, E_Bound_Type::FIXED)); - fit_params.add_parameter(Fit_Param(STR_COMPTON_F_TAIL, (T_real)0.0, (T_real)1.0, (T_real)0.1, (T_real)0.1, E_Bound_Type::LIMITED_LO)); + fit_params.add_parameter(Fit_Param(STR_COMPTON_F_TAIL, (T_real)0.0, (T_real)1.0, (T_real)0.1, (T_real)0.1, E_Bound_Type::LIMITED_LO_HI)); fit_params.add_parameter(Fit_Param(STR_COMPTON_GAMMA, (T_real)0.1, (T_real)10., (T_real)1.0, (T_real)0.1, E_Bound_Type::FIXED)); fit_params.add_parameter(Fit_Param(STR_COMPTON_HI_F_TAIL, (T_real)0.0, (T_real)1.0, (T_real)0.013, (T_real)0.0000001, E_Bound_Type::LIMITED_LO_HI)); fit_params.add_parameter(Fit_Param(STR_COMPTON_HI_GAMMA, (T_real)0.1, (T_real)3., (T_real)1.0, (T_real)0.01, E_Bound_Type::FIXED)); @@ -199,7 +199,7 @@ void Gaussian_Model::set_fit_params_preset(Fit_Params_Preset preset) _fit_parameters[STR_COMPTON_FWHM_CORR].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_AMPLITUDE].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_F_STEP].bound_type = E_Bound_Type::FIXED; - _fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO; + _fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_GAMMA].bound_type = E_Bound_Type::FIXED; _fit_parameters[STR_COMPTON_HI_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_HI_GAMMA].bound_type = E_Bound_Type::FIXED; @@ -237,12 +237,12 @@ void Gaussian_Model::set_fit_params_preset(Fit_Params_Preset preset) _fit_parameters[STR_COMPTON_FWHM_CORR].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_AMPLITUDE].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_F_STEP].bound_type = E_Bound_Type::FIXED; - _fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO; + _fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_GAMMA].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_HI_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_HI_GAMMA].bound_type = E_Bound_Type::LIMITED_LO_HI; - _fit_parameters[STR_SNIP_WIDTH].bound_type = E_Bound_Type::FIT; + _fit_parameters[STR_SNIP_WIDTH].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_F_STEP_OFFSET].bound_type = E_Bound_Type::FIXED; _fit_parameters[STR_F_STEP_LINEAR].bound_type = E_Bound_Type::FIXED; @@ -313,7 +313,7 @@ void Gaussian_Model::set_fit_params_preset(Fit_Params_Preset preset) _fit_parameters[STR_COMPTON_FWHM_CORR].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_AMPLITUDE].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_F_STEP].bound_type = E_Bound_Type::FIXED; - _fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO; + _fit_parameters[STR_COMPTON_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_GAMMA].bound_type = E_Bound_Type::FIXED; _fit_parameters[STR_COMPTON_HI_F_TAIL].bound_type = E_Bound_Type::LIMITED_LO_HI; _fit_parameters[STR_COMPTON_HI_GAMMA].bound_type = E_Bound_Type::FIXED; @@ -586,7 +586,7 @@ const Spectra Gaussian_Model::model_spectrum_element(const Fit_P for (size_t idx = 0; idx < energy_ratios.size(); idx++) { const Element_Energy_Ratio& er_struct = energy_ratios.at(idx); - T_real sigma = std::sqrt(std::pow((fitp->at(STR_FWHM_OFFSET).value / (T_real)2.3548), (T_real)2.0) + (er_struct.energy) * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value); + T_real sigma = std::sqrt( std::pow((fitp->at(STR_FWHM_OFFSET).value / (T_real)2.3548), (T_real)2.0) + (er_struct.energy * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value ) ); T_real f_step = std::abs( er_struct.mu_fraction * ( fitp->at(STR_F_STEP_OFFSET).value + (fitp->at(STR_F_STEP_LINEAR).value * er_struct.energy))); T_real f_tail = std::abs( fitp->at(STR_F_TAIL_OFFSET).value + (fitp->at(STR_F_TAIL_LINEAR).value * er_struct.mu_fraction)); T_real kb_f_tail = std::abs( fitp->at(STR_KB_F_TAIL_OFFSET).value + (fitp->at(STR_KB_F_TAIL_LINEAR).value * er_struct.mu_fraction)); @@ -725,7 +725,19 @@ const ArrayTr Gaussian_Model::peak(T_real gain, T_real sigma, co template const ArrayTr Gaussian_Model::step(T_real gain, T_real sigma, const ArrayTr& delta_energy, T_real peak_E) const { - ArrayTrcounts = delta_energy.unaryExpr([gain, sigma, peak_E](T_real v) { return gain / (T_real)2.0 / peak_E * std::erfc(v / ((T_real)(M_SQRT2)*sigma)); }); + /* + ArrayTrcounts = delta_energy.unaryExpr([gain, sigma, peak_E](T_real v) + { + T_real val = ((T_real)(M_SQRT2) * sigma); + val = std::erfc(v / val); + val = peak_E * val; + val = (T_real)2.0 / val; + return gain / val; + }); + return counts; + */ + return delta_energy.unaryExpr([gain, sigma, peak_E](T_real v) { return gain / (T_real)2.0 / peak_E * std::erfc(v / ((T_real)(M_SQRT2)*sigma)); }); + /* ArrayTrcounts(delta_energy.size()); for (unsigned int i = 0; i < delta_energy.size(); i++) @@ -733,7 +745,6 @@ const ArrayTr Gaussian_Model::step(T_real gain, T_real sigma, co counts[i] = gain / (T_real)2.0 / peak_E * std::erfc((T_real)delta_energy[i] / ((T_real)(M_SQRT2) * sigma)); } */ - return counts; } @@ -764,7 +775,7 @@ const ArrayTr Gaussian_Model::elastic_peak(const Fit_Parameters< { ArrayTr counts(ev.size()); counts.setZero(); - T_real sigma = std::sqrt( std::pow( (fitp->at(STR_FWHM_OFFSET).value / (T_real)2.3548), (T_real)2.0 ) + fitp->at(STR_COHERENT_SCT_ENERGY).value * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value ); + T_real sigma = std::sqrt( std::pow( (fitp->at(STR_FWHM_OFFSET).value / (T_real)2.3548), (T_real)2.0 ) + (fitp->at(STR_COHERENT_SCT_ENERGY).value * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value ) ); if(false == std::isfinite(sigma)) { logE << "sigma = " << sigma << "\n"; @@ -798,7 +809,7 @@ const ArrayTr Gaussian_Model::compton_peak(const Fit_Parameters< T_real compton_E = fitp->at(STR_COHERENT_SCT_ENERGY).value/((T_real)1.0 +(fitp->at(STR_COHERENT_SCT_ENERGY).value / (T_real)511.0 ) * ((T_real)1.0 -std::cos( fitp->at(STR_COMPTON_ANGLE).value * (T_real)2.0 * (T_real)(M_PI) / (T_real)360.0 ))); - T_real sigma = std::sqrt( std::pow( (fitp->at(STR_FWHM_OFFSET).value/(T_real)2.3548), (T_real)62.0) + compton_E * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value ); + T_real sigma = std::sqrt( std::pow( (fitp->at(STR_FWHM_OFFSET).value/(T_real)2.3548), (T_real)2.0) + (compton_E * (T_real)2.96 * fitp->at(STR_FWHM_FANOPRIME).value ) ); if(false == std::isfinite(sigma)) { logE << "sigma = " << sigma << "\n"; diff --git a/src/fitting/optimizers/mpfit_optimizer.cpp b/src/fitting/optimizers/mpfit_optimizer.cpp index 1c29ce54..3542996c 100644 --- a/src/fitting/optimizers/mpfit_optimizer.cpp +++ b/src/fitting/optimizers/mpfit_optimizer.cpp @@ -321,7 +321,7 @@ void MPFit_Optimizer::set_options(std::unordered_map -void MPFit_Optimizer::_fill_limits(Fit_Parameters *fit_params , std::vector > &par) +bool MPFit_Optimizer::_fill_limits(Fit_Parameters *fit_params , std::vector > &par) { for (auto itr = fit_params->begin(); itr != fit_params->end(); itr++) { @@ -346,13 +346,11 @@ void MPFit_Optimizer::_fill_limits(Fit_Parameters *fit_params , || fit.bound_type == E_Bound_Type::LIMITED_LO || fit.bound_type == E_Bound_Type::LIMITED_LO_HI) { - if (fit.max_val == fit.min_val) - { - fit.max_val += (T_real)1.0; - fit.min_val -= (T_real)1.0; - (*fit_params)[itr->first].max_val += (T_real)1.0; - (*fit_params)[itr->first].min_val -= (T_real)1.0; - } + if(fit.min_val == fit.max_val) + { + logE<<"Fit parameter "<::_fill_limits(Fit_Parameters *fit_params , par[fit.opt_array_index].deriv_abstol = (T_real)0.00001; // Absolute tolerance for derivative debug printout } } + return true; } //----------------------------------------------------------------------------- @@ -483,7 +482,7 @@ OPTIMIZER_OUTCOME MPFit_Optimizer::minimize(Fit_Parameters*fit_p std::vector resid(energy_range.count()); std::vector covar(fitp_arr.size() * fitp_arr.size()); - size_t total_itr = num_itr * (fitp_arr.size() + 1); + size_t total_itr = num_itr; fill_user_data(ud, fit_params, spectra, elements_to_fit, model, energy_range, status_callback, total_itr, use_weights); /* @@ -521,7 +520,10 @@ OPTIMIZER_OUTCOME MPFit_Optimizer::minimize(Fit_Parameters*fit_p _options.maxfev = _options.maxiter * ((int)fitp_arr.size() + 1); - _fill_limits(fit_params, par); + if(false == _fill_limits(fit_params, par)) + { + return OPTIMIZER_OUTCOME::FAILED; + } mp_result result; result.xerror = &perror[0]; @@ -669,7 +671,10 @@ OPTIMIZER_OUTCOME MPFit_Optimizer::minimize_func(Fit_Parameters std::vector > par; par.resize(fitp_arr.size()); - _fill_limits(fit_params, par); + if(false == _fill_limits(fit_params, par)) + { + return OPTIMIZER_OUTCOME::FAILED; + } mp_result result; result.xerror = &perror[0]; @@ -821,7 +826,10 @@ OPTIMIZER_OUTCOME MPFit_Optimizer::minimize_quantification(Fit_Parameter std::vector > par; par.resize(fitp_arr.size()); - _fill_limits(fit_params, par); + if(false == _fill_limits(fit_params, par)) + { + return OPTIMIZER_OUTCOME::FAILED; + } this->_last_outcome = mpfit(quantification_residuals_mpfit, (int)ud.quant_map.size(), (int)fitp_arr.size(), &fitp_arr[0], &par[0], &_options, (void *) &ud, &result); logI << "\nOutcome: " << optimizer_outcome_to_str(this->_outcome_map[this->_last_outcome]) << "\nNum iter: " << result.niter << "\n Norm of the residue vector: " << *result.resid << "\n"; diff --git a/src/fitting/optimizers/mpfit_optimizer.h b/src/fitting/optimizers/mpfit_optimizer.h index a3fb3859..034a2f2c 100644 --- a/src/fitting/optimizers/mpfit_optimizer.h +++ b/src/fitting/optimizers/mpfit_optimizer.h @@ -95,9 +95,9 @@ class DLL_EXPORT MPFit_Optimizer: public Optimizer private: - void _fill_limits(Fit_Parameters *fit_params, std::vector > &par); + bool _fill_limits(Fit_Parameters *fit_params, std::vector > &par); - struct mp_config _options; + struct mp_config _options; }; diff --git a/src/io/file/aps/aps_fit_params_import.h b/src/io/file/aps/aps_fit_params_import.h index 09a15421..1591d2bd 100644 --- a/src/io/file/aps/aps_fit_params_import.h +++ b/src/io/file/aps/aps_fit_params_import.h @@ -1076,7 +1076,7 @@ DLL_EXPORT bool save_parameters_override(std::string path, Params_Override -DLL_EXPORT bool save_fit_parameters_override(std::string path, Fit_Parameters fit_params, std::string result) +DLL_EXPORT bool save_fit_parameters_override(const std::string& path, const Fit_Parameters& fit_params, const std::string& result) { std::ofstream out_stream(path);