diff --git a/src/fitting/optimizers/lmfit_optimizer.cpp b/src/fitting/optimizers/lmfit_optimizer.cpp index c5495869..92010c21 100644 --- a/src/fitting/optimizers/lmfit_optimizer.cpp +++ b/src/fitting/optimizers/lmfit_optimizer.cpp @@ -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"; @@ -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]; @@ -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++; } @@ -229,7 +229,8 @@ std::unordered_map LMFit_Optimizer::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; } @@ -267,6 +268,10 @@ void LMFit_Optimizer::set_options(std::unordered_map 0) + { + _options.verbosity = (int)opt.at(STR_OPT_VERBOSE_LEVEL); + } } // ---------------------------------------------------------------------------- @@ -289,6 +294,7 @@ OPTIMIZER_OUTCOME LMFit_Optimizer::minimize(Fit_Parameters *fit_ } std::vector 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); @@ -398,7 +404,7 @@ OPTIMIZER_OUTCOME LMFit_Optimizer::minimize_func(Fit_Parameters* std::vector perror(fitp_arr.size()); lm_status_struct 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); @@ -444,7 +450,7 @@ OPTIMIZER_OUTCOME LMFit_Optimizer::minimize_quantification(Fit_Parameter } ud.quantification_model = quantification_model; ud.fit_parameters = fit_params; - + //_options.verbosity = 2; std::vector fitp_arr = fit_params->to_array(); if (fitp_arr.size() == 0) { diff --git a/src/fitting/optimizers/mpfit_optimizer.cpp b/src/fitting/optimizers/mpfit_optimizer.cpp index fa8ce710..044f891d 100644 --- a/src/fitting/optimizers/mpfit_optimizer.cpp +++ b/src/fitting/optimizers/mpfit_optimizer.cpp @@ -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; ispectra[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"; @@ -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; ispectra[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]; @@ -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++; } @@ -305,14 +306,14 @@ void MPFit_Optimizer::_fill_limits(Fit_Parameters *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 @@ -410,6 +411,24 @@ std::string MPFit_Optimizer::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"; } @@ -485,7 +504,7 @@ OPTIMIZER_OUTCOME MPFit_Optimizer::minimize(Fit_Parameters*fit_p this->_last_outcome = mpfit(residuals_mpfit, (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); @@ -728,7 +747,7 @@ OPTIMIZER_OUTCOME MPFit_Optimizer::minimize_quantification(Fit_Parameter ud.fit_parameters = fit_params; std::vector 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; } @@ -781,7 +800,7 @@ OPTIMIZER_OUTCOME MPFit_Optimizer::minimize_quantification(Fit_Parameter 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"; - logI << detailed_outcome(this->_last_outcome); + logI << detailed_outcome(this->_last_outcome) << "\n"; fit_params->from_array(fitp_arr); diff --git a/src/fitting/optimizers/optimizer.h b/src/fitting/optimizers/optimizer.h index 08836e5a..33b7d532 100644 --- a/src/fitting/optimizers/optimizer.h +++ b/src/fitting/optimizers/optimizer.h @@ -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" @@ -185,7 +186,7 @@ void fill_user_data(User_Data &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 weights = (*spectra); diff --git a/src/support/lmfit_6.1/lmmin.hpp b/src/support/lmfit_6.1/lmmin.hpp index 7249249f..f6c810d7 100644 --- a/src/support/lmfit_6.1/lmmin.hpp +++ b/src/support/lmfit_6.1/lmmin.hpp @@ -67,8 +67,10 @@ void lm_print_pars(int nout, const _T* par, FILE* fout) { int i; for (i = 0; i < nout; ++i) - fprintf(fout, " %16.9g", par[i]); - fprintf(fout, "\n"); + { + logit_s << par[i] << " "; + } + logit_s<<"\n"; } /******************************************************************************/ @@ -228,7 +230,9 @@ void lm_qrsolv(const int n, _T* r, const int ldr, const int* Pivot, for (j = 0; j < n; j++) { for (i = j; i < n; i++) - r[j*ldr+i] = r[i*ldr+j]; + { + r[j * ldr + i] = r[i * ldr + j]; + } x[j] = r[j*ldr+j]; W[j] = qtb[j]; } @@ -242,7 +246,9 @@ void lm_qrsolv(const int n, _T* r, const int ldr, const int* Pivot, if (diag[Pivot[j]] != 0) { for (k = j; k < n; k++) + { Sdiag[k] = 0; + } Sdiag[j] = diag[Pivot[j]]; /*** The transformations to eliminate the row of D modify only @@ -300,17 +306,21 @@ void lm_qrsolv(const int n, _T* r, const int ldr, const int* Pivot, W[j] = 0; } - for (j = nsing-1; j >= 0; j--) { + for (j = nsing - 1; j >= 0; j--) { sum = 0; - for (i = j+1; i < nsing; i++) - sum += r[j*ldr+i] * W[i]; + for (i = j + 1; i < nsing; i++) + { + sum += r[j * ldr + i] * W[i]; + } W[j] = (W[j] - sum) / Sdiag[j]; } /*** Permute the components of z back to components of x. ***/ for (j = 0; j < n; j++) + { x[Pivot[j]] = W[j]; + } } /*** lm_qrsolv. ***/ @@ -426,13 +436,16 @@ void lm_lmpar(const int n, _T* r, const int ldr, const int* Pivot, } for (j = 0; j < n; j++) + { x[Pivot[j]] = aux[j]; - + } /*** Initialize the iteration counter, evaluate the function at the origin, and test for acceptance of the Gauss-Newton direction. ***/ for (j = 0; j < n; j++) + { xdi[j] = diag[j] * x[j]; + } dxnorm = lm_enorm(n, xdi); fp = dxnorm - delta; if (fp <= p1 * delta) { @@ -450,12 +463,15 @@ void lm_lmpar(const int n, _T* r, const int ldr, const int* Pivot, parl = 0; if (nsing >= n) { for (j = 0; j < n; j++) + { aux[j] = diag[Pivot[j]] * xdi[Pivot[j]] / dxnorm; - + } for (j = 0; j < n; j++) { sum = 0; for (i = 0; i < j; i++) - sum += r[j*ldr+i] * aux[i]; + { + sum += r[j * ldr + i] * aux[i]; + } aux[j] = (aux[j] - sum) / r[j+ldr*j]; } temp = lm_enorm(n, aux); @@ -467,7 +483,9 @@ void lm_lmpar(const int n, _T* r, const int ldr, const int* Pivot, for (j = 0; j < n; j++) { sum = 0; for (i = 0; i <= j; i++) - sum += r[j*ldr+i] * qtb[i]; + { + sum += r[j * ldr + i] * qtb[i]; + } aux[j] = sum / diag[Pivot[j]]; } gnorm = lm_enorm(n, aux); @@ -492,13 +510,16 @@ void lm_lmpar(const int n, _T* r, const int ldr, const int* Pivot, *par = MAX(LM_DWARF, (_T)0.001 * paru); temp = sqrt(*par); for (j = 0; j < n; j++) + { aux[j] = temp * diag[j]; - + } lm_qrsolv(n, r, ldr, Pivot, aux, qtb, x, Sdiag, xdi); /* return values are r, x, Sdiag */ for (j = 0; j < n; j++) + { xdi[j] = diag[j] * x[j]; /* used as output */ + } dxnorm = lm_enorm(n, xdi); fp_old = fp; fp = dxnorm - delta; @@ -518,12 +539,15 @@ void lm_lmpar(const int n, _T* r, const int ldr, const int* Pivot, /** Compute the Newton correction. **/ for (j = 0; j < n; j++) + { aux[j] = diag[Pivot[j]] * xdi[Pivot[j]] / dxnorm; - + } for (j = 0; j < n; j++) { aux[j] = aux[j] / Sdiag[j]; - for (i = j+1; i < n; i++) - aux[i] -= r[j*ldr+i] * aux[j]; + for (i = j + 1; i < n; i++) + { + aux[i] -= r[j * ldr + i] * aux[j]; + } } temp = lm_enorm(n, aux); parc = fp / delta / temp / temp; @@ -616,10 +640,13 @@ void lm_qrfac(const int m, const int n, _T* A, int* Pivot, _T* Rdiag, /** Bring the column of largest norm into the pivot position. **/ kmax = j; - for (k = j+1; k < n; k++) + for (k = j + 1; k < n; k++) + { if (Rdiag[k] > Rdiag[kmax]) + { kmax = k; - + } + } if (kmax != j) { /* Swap columns j and kmax. */ k = Pivot[j]; @@ -648,7 +675,9 @@ void lm_qrfac(const int m, const int n, _T* A, int* Pivot, _T* Rdiag, if (A[j*m+j] < 0) ajnorm = -ajnorm; for (i = j; i < m; i++) - A[j*m+i] /= ajnorm; + { + A[j * m + i] /= ajnorm; + } A[j*m+j] += 1; /** Apply the Householder transformation U_w := 1 - 2*w_j.w_j/|w_j|^2 @@ -657,15 +686,17 @@ void lm_qrfac(const int m, const int n, _T* A, int* Pivot, _T* Rdiag, /* Compute scalar product w_j * a_j. */ sum = 0; for (i = j; i < m; i++) - sum += A[j*m+i] * A[k*m+i]; - + { + sum += A[j * m + i] * A[k * m + i]; + } /* Normalization is simplified by the coincidence |w_j|^2=2w_jj. */ temp = sum / A[j*m+j]; /* Carry out transform U_w_j * a_k. */ for (i = j; i < m; i++) - A[k*m+i] -= temp * A[j*m+i]; - + { + A[k * m + i] -= temp * A[j * m + i]; + } /* No idea what happens here. */ if (Rdiag[k] != 0) { temp = A[m*k+j] / Rdiag[k]; @@ -822,39 +853,32 @@ void lmmin(const int n, _T* x, const int m, const void* data, /*** Check input parameters for errors. ***/ if (n <= 0) { - fprintf(stderr, "lmmin: invalid number of parameters %i\n", n); + logE<<"lmmin: invalid number of parameters "<outcome = 10; return; } if (m < n) { - fprintf(stderr, "lmmin: number of data points (%i) " - "smaller than number of parameters (%i)\n", - m, n); + logE << "lmmin: number of data points (" << m << ") smaller than number of parameters (" << n << ")\n"; S->outcome = 10; return; } if (C->ftol < 0 || C->xtol < 0 || C->gtol < 0) { - fprintf(stderr, - "lmmin: negative tolerance (at least one of %g %g %g)\n", - C->ftol, C->xtol, C->gtol); + logE << "lmmin: negative tolerance (at least one of " << C->ftol << " " << C->xtol << " " << C->gtol << ")\n"; S->outcome = 10; return; } if (maxfev <= 0) { - fprintf(stderr, "lmmin: nonpositive function evaluations limit %i\n", - maxfev); + logE << "lmmin: nonpositive function evaluations limit " << maxfev << "\n"; S->outcome = 10; return; } if (C->stepbound <= 0) { - fprintf(stderr, "lmmin: nonpositive stepbound %g\n", C->stepbound); + logE << "lmmin: nonpositive stepbound "<< C->stepbound<<"\n"; S->outcome = 10; return; } if (C->scale_diag != 0 && C->scale_diag != 1) { - fprintf(stderr, "lmmin: logical variable scale_diag=%i, " - "should be 0 or 1\n", - C->scale_diag); + logE << "lmmin: logical variable scale_diag=" << C->scale_diag << ", should be 0 or 1\n"; S->outcome = 10; return; } @@ -892,25 +916,32 @@ void lmmin(const int n, _T* x, const int m, const void* data, /* Initialize diag. */ if (!C->scale_diag) + { for (j = 0; j < n; j++) + { diag[j] = 1; - + } + } /*** Evaluate function at starting point and calculate norm. ***/ if (C->verbosity) { - fprintf(msgfile, "lmmin start "); + logI<<"lmmin start \n"; lm_print_pars(nout, x, msgfile); } (*evaluate)(x, m, data, fvec, &(S->userbreak)); if (C->verbosity > 4) + { for (i = 0; i < m; ++i) - fprintf(msgfile, " fvec[%4i] = %18.8g\n", i, fvec[i]); + { + logI << " fvec[" << i << "] = " << fvec[i] << "\n"; + } + } S->nfev = 1; if (S->userbreak) goto terminate; fnorm = lm_enorm(m, fvec); if (C->verbosity) - fprintf(msgfile, " fnorm = %18.8g\n", fnorm); + logI << " fnorm = "<< fnorm<<"\n"; if (!std::isfinite(fnorm)) { S->outcome = 12; /* nan */ @@ -934,7 +965,9 @@ void lmmin(const int n, _T* x, const int m, const void* data, if (S->userbreak) goto terminate; for (i = 0; i < m; i++) - fjac[j*m+i] = (wf[i] - fvec[i]) / step; + { + fjac[j * m + i] = (wf[i] - fvec[i]) / step; + } x[j] = temp; /* restore */ } if (C->verbosity >= 10) { @@ -943,7 +976,9 @@ void lmmin(const int n, _T* x, const int m, const void* data, for (i = 0; i < m; i++) { printf(" "); for (j = 0; j < n; j++) - printf("%.5e ", fjac[j*m+i]); + { + printf("%.5e ", fjac[j * m + i]); + } printf("\n"); } } @@ -982,10 +1017,14 @@ void lmmin(const int n, _T* x, const int m, const void* data, if (temp3 != 0) { sum = 0; for (i = j; i < m; i++) - sum += fjac[j*m+i] * wf[i]; + { + sum += fjac[j * m + i] * wf[i]; + } temp = -sum / temp3; for (i = j; i < m; i++) - wf[i] += fjac[j*m+i] * temp; + { + wf[i] += fjac[j * m + i] * temp; + } } fjac[j*m+j] = wa1[j]; qtf[j] = wf[j]; @@ -998,7 +1037,9 @@ void lmmin(const int n, _T* x, const int m, const void* data, continue; sum = 0; for (i = 0; i <= j; i++) - sum += fjac[j*m+i] * qtf[i]; + { + sum += fjac[j * m + i] * qtf[i]; + } gnorm = MAX(gnorm, std::fabs(sum / wa2[Pivot[j]] / fnorm)); } @@ -1012,24 +1053,28 @@ void lmmin(const int n, _T* x, const int m, const void* data, if (C->scale_diag) { /* diag := norms of the columns of the initial Jacobian */ for (j = 0; j < n; j++) + { diag[j] = wa2[j] ? wa2[j] : 1; + } /* xnorm := || D x || */ for (j = 0; j < n; j++) + { wa3[j] = diag[j] * x[j]; + } xnorm = lm_enorm(n, wa3); if (C->verbosity >= 2) { - fprintf(msgfile, "lmmin diag "); + logI << "lmmin diag "; lm_print_pars(nout, x, msgfile); // xnorm - fprintf(msgfile, " xnorm = %18.8g\n", xnorm); + logI<<" xnorm = "<verbosity >= 3) { - fprintf(msgfile, " o i lmpar prered" - " ratio dirder delta" - " pnorm fnorm"); + logI << " o i lmpar prered ratio dirder delta pnorm fnorm\n"; for (i = 0; i < nout; ++i) - fprintf(msgfile, " p%i", i); - fprintf(msgfile, "\n"); + { + logit_s << " p" << i; + } + logit_s<<"\n"; } } else { xnorm = lm_enorm(n, x); @@ -1105,16 +1150,15 @@ void lmmin(const int n, _T* x, const int m, const void* data, ratio = prered ? actred / prered : 0; if (C->verbosity == 2) { - fprintf(msgfile, "lmmin (%i:%i) ", outer, inner); + logI << "lmmin ("<verbosity >= 3) { - printf("%3i %2i %9.2g %9.2g %14.6g" - " %9.2g %10.3e %10.3e %21.15e", - outer, inner, lmpar, prered, ratio, - dirder, delta, pnorm, fnorm1); + logit_s << outer << " " << inner << " " << lmpar << " " << prered << " " << ratio << " " << dirder << " " << delta << " " << pnorm << " " << fnorm1; for (i = 0; i < nout; ++i) - fprintf(msgfile, " %16.9g", wa2[i]); - fprintf(msgfile, "\n"); + { + logit_s << " " << wa2[i]; + } + logit_s<< "\n"; } /* Update the step bound. */ @@ -1147,10 +1191,14 @@ void lmmin(const int n, _T* x, const int m, const void* data, } } else { for (j = 0; j < n; j++) + { x[j] = wa2[j]; + } } for (i = 0; i < m; i++) + { fvec[i] = wf[i]; + } xnorm = lm_enorm(n, wa2); if (!std::isfinite(xnorm)) { S->outcome = 12; /* nan */ @@ -1200,12 +1248,11 @@ void lmmin(const int n, _T* x, const int m, const void* data, terminate: S->fnorm = lm_enorm(m, fvec); if (C->verbosity >= 2) - printf("lmmin outcome (%i) xnorm %g ftol %g xtol %g\n", S->outcome, - xnorm, C->ftol, C->xtol); + logI << "lmmin outcome (" << S->outcome << ") xnorm " << xnorm << " ftol " << C->ftol << " xtol " << C->xtol << "\n"; if (C->verbosity & 1) { - fprintf(msgfile, "lmmin final "); + logI << "lmmin final "; lm_print_pars(nout, x, msgfile); // S->fnorm, - fprintf(msgfile, " fnorm = %18.8g\n", S->fnorm); + logI << " fnorm = " << S->fnorm << "\n"; } if (S->userbreak) /* user-requested break */ S->outcome = 11;