Skip to content

Commit

Permalink
Merge pull request #174 from aglowacki/master
Browse files Browse the repository at this point in the history 
Changed residual calc
  • Loading branch information
@aglowacki
aglowacki authored Dec 5, 2023
2 parents 473405b + ebb6d02 commit 6ef1958
Show file tree
Hide file tree
Showing 4 changed files with 150 additions and 77 deletions.
18 changes: 12 additions & 6 deletions src/fitting/optimizers/lmfit_optimizer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ void residuals_lmfit( const T_real *par, int m_dat, const void *data, T_real *fv
// Calculate residuals
for (int i = 0; i < m_dat; i++ )
{
fvec[i] = pow((ud->spectra[i] - ud->spectra_model[i]),2) * ud->weights[i];
fvec[i] = abs(ud->spectra[i] - ud->spectra_model[i]) * ud->weights[i];
if (std::isfinite(fvec[i]) == false)
{
//logE << "\n\n\n";
Expand Down Expand Up @@ -124,7 +124,7 @@ void general_residuals_lmfit( const T_real *par, int m_dat, const void *data, T_
// Calculate residuals
for (int i = 0; i < m_dat; i++ )
{
fvec[i] = pow(( ud->spectra[i] - ud->spectra_model[i] ),2) * ud->weights[i];
fvec[i] = abs( ud->spectra[i] - ud->spectra_model[i] ) * ud->weights[i];
if (std::isfinite(fvec[i]) == false)
{
fvec[i] = ud->spectra[i];
Expand Down Expand Up @@ -163,7 +163,7 @@ void quantification_residuals_lmfit( const T_real *par, int m_dat, const void *d
}
else
{
fvec[idx] = itr.second.e_cal_ratio - result_map[itr.first];
fvec[idx] = abs(itr.second.e_cal_ratio - result_map[itr.first]);
}
idx++;
}
Expand Down Expand Up @@ -229,7 +229,8 @@ std::unordered_map<std::string, T_real> LMFit_Optimizer<T_real>::get_options()
{STR_OPT_EPSILON, _options.epsilon},
{STR_OPT_STEP, _options.stepbound},
{STR_OPT_SCALE_DIAG, _options.scale_diag},
{STR_OPT_MAXITER, _options.patience}
{STR_OPT_MAXITER, _options.patience},
{STR_OPT_VERBOSE_LEVEL, _options.verbosity}
};
return opts;
}
Expand Down Expand Up @@ -267,6 +268,10 @@ void LMFit_Optimizer<T_real>::set_options(std::unordered_map<std::string, T_real
{
_options.patience = (int)opt.at(STR_OPT_MAXITER);
}
if (opt.count(STR_OPT_VERBOSE_LEVEL) > 0)
{
_options.verbosity = (int)opt.at(STR_OPT_VERBOSE_LEVEL);
}
}

// ----------------------------------------------------------------------------
Expand All @@ -289,6 +294,7 @@ OPTIMIZER_OUTCOME LMFit_Optimizer<T_real>::minimize(Fit_Parameters<T_real> *fit_
}
std::vector<T_real> perror(fitp_arr.size());

//_options.verbosity = 2;
size_t total_itr = _options.patience * (fitp_arr.size() + 1);
fill_user_data(ud, fit_params, spectra, elements_to_fit, model, energy_range, status_callback, total_itr, use_weights);

Expand Down Expand Up @@ -398,7 +404,7 @@ OPTIMIZER_OUTCOME LMFit_Optimizer<T_real>::minimize_func(Fit_Parameters<T_real>*
std::vector<T_real> perror(fitp_arr.size());

lm_status_struct<T_real> status;

_options.verbosity = 0;
lmmin((int)fitp_arr.size(), &fitp_arr[0], (int)energy_range.count(), (const void*) &ud, general_residuals_lmfit, &_options, &status );
this->_last_outcome = status.outcome;
fit_params->from_array(fitp_arr);
Expand Down Expand Up @@ -444,7 +450,7 @@ OPTIMIZER_OUTCOME LMFit_Optimizer<T_real>::minimize_quantification(Fit_Parameter
}
ud.quantification_model = quantification_model;
ud.fit_parameters = fit_params;

//_options.verbosity = 2;
std::vector<T_real> fitp_arr = fit_params->to_array();
if (fitp_arr.size() == 0)
{
Expand Down
37 changes: 28 additions & 9 deletions src/fitting/optimizers/mpfit_optimizer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,8 @@ int residuals_mpfit(int m, int params_size, T_real *params, T_real *dy, T_real *
//Calculate residuals
for (int i=0; i<m; i++)
{
dy[i] = pow((ud->spectra[i] - ud->spectra_model[i]),2) * ud->weights[i];
//dy[i] = pow((ud->spectra[i] - ud->spectra_model[i]),2) * ud->weights[i];
dy[i] = abs(ud->spectra[i] - ud->spectra_model[i]) * ud->weights[i];
if (std::isfinite(dy[i]) == false)
{
//logE << "\n\n\n";
Expand Down Expand Up @@ -132,7 +133,7 @@ int gen_residuals_mpfit(int m, int params_size, T_real *params, T_real *dy, T_re
// Calculate residuals
for (int i=0; i<m; i++)
{
dy[i] = pow(( ud->spectra[i] - ud->spectra_model[i] ),2) * ud->weights[i];
dy[i] = abs( ud->spectra[i] - ud->spectra_model[i] ) * ud->weights[i];
if (std::isfinite(dy[i]) == false)
{
dy[i] = ud->spectra[i];
Expand Down Expand Up @@ -173,7 +174,7 @@ int quantification_residuals_mpfit(int m, int params_size, T_real *params, T_rea
}
else
{
dy[idx] = itr.second.e_cal_ratio - result_map[itr.first];
dy[idx] = abs(itr.second.e_cal_ratio - result_map[itr.first]);
}
idx++;
}
Expand Down Expand Up @@ -305,14 +306,14 @@ void MPFit_Optimizer<T_real>::_fill_limits(Fit_Parameters<T_real> *fit_params ,
if (fit.value > fit.max_val)
{
logW << itr->first << " value (" << fit.value << ") > max_val(" << fit.max_val << ") : setting value = max_val\n";
fit.value = fit.max_val;
(*fit_params)[itr->first].value = fit.max_val;
fit.value = fit.max_val - fit.step_size;
(*fit_params)[itr->first].value = fit.value;

}
if (fit.value < fit.min_val)
{
logW << itr->first << " value (" << fit.value << ") < min_val(" << fit.min_val << ") : setting value = min_val\n";
fit.value = fit.min_val;
fit.value = fit.min_val + fit.step_size;
(*fit_params)[itr->first].value = fit.min_val;
}
if (fit.bound_type == E_Bound_Type::LIMITED_HI
Expand Down Expand Up @@ -410,6 +411,24 @@ std::string MPFit_Optimizer<T_real>::detailed_outcome(int info)
return "Xtol is too small. no further improvement in the approximate solution x is possible.";
case 8:
return "Gtol is too small. fvec is orthogonal to the columns of the jacobian to machine precision.";
case MP_ERR_NAN:
return "User function produced non-finite values.";
case MP_ERR_FUNC:
return "No user function was supplied.";
case MP_ERR_NPOINTS:
return "No user data points were supplied.";
case MP_ERR_NFREE:
return "No free parameters.";
case MP_ERR_MEMORY:
return "Memory allocation error.";
case MP_ERR_INITBOUNDS:
return "Initial values inconsistent w constraints.";
case MP_ERR_BOUNDS:
return "Initial constraints inconsistent.";
case MP_ERR_PARAM:
return "General input parameter error.";
case MP_ERR_DOF:
return "Not enough degrees of freedom.";
default:
return "Unknown info status";
}
Expand Down Expand Up @@ -485,7 +504,7 @@ OPTIMIZER_OUTCOME MPFit_Optimizer<T_real>::minimize(Fit_Parameters<T_real>*fit_p

this->_last_outcome = mpfit(residuals_mpfit<T_real>, (int)energy_range.count(), (int)fitp_arr.size(), &fitp_arr[0], &par[0], &_options, (void *) &ud, &result);

logI<< detailed_outcome(this->_last_outcome);
logI<< detailed_outcome(this->_last_outcome) << "\n";

fit_params->from_array(fitp_arr);

Expand Down Expand Up @@ -728,7 +747,7 @@ OPTIMIZER_OUTCOME MPFit_Optimizer<T_real>::minimize_quantification(Fit_Parameter
ud.fit_parameters = fit_params;

std::vector<T_real> fitp_arr = fit_params->to_array();
if (fitp_arr.size() == 0)
if (fitp_arr.size() == 0 || ud.quant_map.size() == 0)
{
return OPTIMIZER_OUTCOME::STOPPED;
}
Expand Down Expand Up @@ -781,7 +800,7 @@ OPTIMIZER_OUTCOME MPFit_Optimizer<T_real>::minimize_quantification(Fit_Parameter
this->_last_outcome = mpfit(quantification_residuals_mpfit<T_real>, (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";

logI << detailed_outcome(this->_last_outcome);
logI << detailed_outcome(this->_last_outcome) << "\n";

fit_params->from_array(fitp_arr);

Expand Down
3 changes: 2 additions & 1 deletion src/fitting/optimizers/optimizer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ using namespace fitting::models;
//LM
#define STR_OPT_SCALE_DIAG "scale_diag"
#define STR_OPT_MAXITER "maxiter"
#define STR_OPT_VERBOSE_LEVEL "verbose_level"
//MP
#define STR_OPT_COVTOL "covtol"

Expand Down Expand Up @@ -185,7 +186,7 @@ void fill_user_data(User_Data<T_real> &ud,
weights = convolve1d(weights, 5);
weights = Eigen::abs(weights);
weights /= weights.maxCoeff();
weights = weights.unaryExpr([](T_real v) { return std::isfinite(v) ? v : (T_real)0.0; });
weights = weights.unaryExpr([](T_real v) { return std::isfinite(v) ? v : (T_real)1.0; });
ud.weights = weights.segment(energy_range.min, energy_range.count());
/*
ArrayTr<T_real> weights = (*spectra);
Expand Down
Loading

0 comments on commit 6ef1958

Please sign in to comment.