Commit 934e72a7 authored by Lefebvre, Jordan P's avatar Lefebvre, Jordan P

Applying clang-format formatting to radixmin.

parent 0397fd8f
Pipeline #48314 failed with stages
in 1 minute
......@@ -47,6 +47,7 @@ ELSE()
# Process subpackages
#
TRIBITS_PROCESS_SUBPACKAGES()
TRIBITS_PACKAGE_DEF()
TRIBITS_PACKAGE_POSTPROCESS()
ENDIF()
......
This diff is collapsed.
#include <memory> // std::shared_ptr
#include <algorithm> // std::find
#include "radixmin/problem.hh" // ProblemInterface
#include "radixmin/solver.hh" // Solver
#include "radixbug/bug.hh"
#include "radixmin/problem.hh" // ProblemInterface
#include "radixmin/solver.hh" // Solver
#include <algorithm> // std::find
#include <memory> // std::shared_ptr
#pragma once
namespace radixmin
{
namespace radixmin {
/** An adapter, e.g., active set. Member functions default to just passing
* calls from m_solver to m_problem and vice versa. */
class Adapter : public ProblemInterface, public Solver
{
class Adapter : public ProblemInterface, public Solver {
protected:
// inherits m_problem from Solver
/// need a handle on the down-chain solver to call its step(), etc
SPSolver m_solver;
// inherits m_problem from Solver
/// need a handle on the down-chain solver to call its step(), etc
SPSolver m_solver;
public:
Adapter(SPProblemInterface problem);
void check_consistency() const override;
void set_solver(SPSolver solver) { m_solver = solver; }
// implement the ProblemInterface interface
FGH obj(const Array& x, int flags=F+G+H) override;
LinearSystem icon() const override { return m_problem->icon(); }
LinearSystem econ() const override { return m_problem->econ(); }
Array x() const override { return m_problem->x(); }
void x(const Array& newval) override { m_problem->x(newval); }
double gradient_tolerance() const final override {
return m_problem->gradient_tolerance(); }
// implement the Solver interface
void refresh() override;
bool converging_rapidly() const override;
bool converged() const override;
void step() override;
Adapter(SPProblemInterface problem);
void check_consistency() const override;
void set_solver(SPSolver solver) { m_solver = solver; }
// implement the ProblemInterface interface
FGH obj(const Array &x, int flags = F + G + H) override;
LinearSystem icon() const override { return m_problem->icon(); }
LinearSystem econ() const override { return m_problem->econ(); }
Array x() const override { return m_problem->x(); }
void x(const Array &newval) override { m_problem->x(newval); }
double gradient_tolerance() const final override {
return m_problem->gradient_tolerance();
}
// implement the Solver interface
void refresh() override;
bool converging_rapidly() const override;
bool converged() const override;
void step() override;
};
typedef std::shared_ptr<Adapter> SPAdapter;
......@@ -48,31 +47,30 @@ typedef std::shared_ptr<Adapter> SPAdapter;
/** The SVDAdapter uses singular value decomposition to handle equality
* constraints. The new problem has fewer variables than the original
* problem. */
class SVDAdapter : public Adapter
{
class SVDAdapter : public Adapter {
protected:
/// svd of m_problem's econ
SVD m_svd;
/// svd of m_problem's econ
SVD m_svd;
public:
SVDAdapter(SPProblemInterface problem);
void check_consistency() const override;
// implement the ProblemInterface interface
/// reverse-transform x when calling obj
FGH obj(const Array& y_noninflu, int flags=F+G+H) override;
/// the modified problem has no icon
LinearSystem icon() const override;
/// the modified problem has no econ
LinearSystem econ() const override;
/// x is transformed when getting x
Array x() const override;
/// x is reverse-transformed when setting x
void x(const Array& newval) override;
// implement the Solver interface
void refresh() override;
SVDAdapter(SPProblemInterface problem);
void check_consistency() const override;
// implement the ProblemInterface interface
/// reverse-transform x when calling obj
FGH obj(const Array &y_noninflu, int flags = F + G + H) override;
/// the modified problem has no icon
LinearSystem icon() const override;
/// the modified problem has no econ
LinearSystem econ() const override;
/// x is transformed when getting x
Array x() const override;
/// x is reverse-transformed when setting x
void x(const Array &newval) override;
// implement the Solver interface
void refresh() override;
};
typedef std::shared_ptr<SVDAdapter> SPSVDAdapter;
......@@ -81,52 +79,50 @@ typedef std::shared_ptr<SVDAdapter> SPSVDAdapter;
/** An active set is a list of which icons are temporarily made into icons.
* Used by ActiveSetAdapter. */
class ActiveSet
{
class ActiveSet {
protected:
size_t TRUE_ECON_SENTINEL = size_t(-1);
class EconNotActive : public std::runtime_error
{
using std::runtime_error::runtime_error;
};
/// icons from problem
LinearSystem m_icon;
/// econs from problem plus activated icons
LinearSystem m_econ;
/// m_index_map[k_econ] = k_icon (or TRUE_ECON)
VectorSizeT m_index_map;
size_t TRUE_ECON_SENTINEL = size_t(-1);
class EconNotActive : public std::runtime_error {
using std::runtime_error::runtime_error;
};
/// icons from problem
LinearSystem m_icon;
/// econs from problem plus activated icons
LinearSystem m_econ;
/// m_index_map[k_econ] = k_icon (or TRUE_ECON)
VectorSizeT m_index_map;
public:
ActiveSet(SPProblemInterface problem);
void check_consistency() const;
// get functions
const LinearSystem& icon() const { return m_icon; }
const LinearSystem& econ() const { return m_econ; }
const VectorSizeT& index_map() const { return m_index_map; }
// get status
/** Look up the k_econ corresponding to an active k_icon, or throw
* EconNotActive */
size_t econ_index(size_t k_icon);
/// Whether an icon is active
bool contains(size_t k_icon) const;
/// True for active icons, false for inactive icons
VectorBool active_flags() const;
/// False for active icons, true for inactive icons
VectorBool inactive_flags() const;
// set status
/// Add an icon to the active set
size_t activate(size_t k_icon);
/// Remove an icon from the active set
size_t deactivate(size_t k_icon);
// get icon filtered by status
LinearSystem active_icon() const;
LinearSystem inactive_icon() const;
ActiveSet(SPProblemInterface problem);
void check_consistency() const;
// get functions
const LinearSystem &icon() const { return m_icon; }
const LinearSystem &econ() const { return m_econ; }
const VectorSizeT &index_map() const { return m_index_map; }
// get status
/** Look up the k_econ corresponding to an active k_icon, or throw
* EconNotActive */
size_t econ_index(size_t k_icon);
/// Whether an icon is active
bool contains(size_t k_icon) const;
/// True for active icons, false for inactive icons
VectorBool active_flags() const;
/// False for active icons, true for inactive icons
VectorBool inactive_flags() const;
// set status
/// Add an icon to the active set
size_t activate(size_t k_icon);
/// Remove an icon from the active set
size_t deactivate(size_t k_icon);
// get icon filtered by status
LinearSystem active_icon() const;
LinearSystem inactive_icon() const;
};
//~ #####################################################################
......@@ -136,37 +132,35 @@ public:
* "active". Active constraints are turned into econs, which are enforced
* using a different adapter, e.g., SVDAdapter. An icon is only declared
* inactive if the gradient pushes x away from the boundary. */
class ActiveSetAdapter : public Adapter
{
class ActiveSetAdapter : public Adapter {
protected:
/// which icons are being enforced as econs
ActiveSet m_active_set;
/// which icons are being enforced as econs
ActiveSet m_active_set;
public:
ActiveSetAdapter(SPProblemInterface problem);
void check_consistency() const override;
/** Adjust x so that it obeys the icon. Update the active_set as well.
* Returns true if any changes were made. */
bool project_inbounds(ActiveSet& active_set, Array& x);
// implement the ProblemInterface interface
bool variable_econ() const override { return true; }
/// obj forces x to obey the icon
FGH obj(const Array& x, int flags=F+G+H) override;
/// the modified problem has no icon
LinearSystem icon() const override;
/// the econs vary depending on which icons are active
LinearSystem econ() const override;
/// when setting x, update m_active_set and force x inbounds
void x(const Array& newval) override;
// implement the Solver interface
void refresh() override;
ActiveSetAdapter(SPProblemInterface problem);
void check_consistency() const override;
/** Adjust x so that it obeys the icon. Update the active_set as well.
* Returns true if any changes were made. */
bool project_inbounds(ActiveSet &active_set, Array &x);
// implement the ProblemInterface interface
bool variable_econ() const override { return true; }
/// obj forces x to obey the icon
FGH obj(const Array &x, int flags = F + G + H) override;
/// the modified problem has no icon
LinearSystem icon() const override;
/// the econs vary depending on which icons are active
LinearSystem econ() const override;
/// when setting x, update m_active_set and force x inbounds
void x(const Array &newval) override;
// implement the Solver interface
void refresh() override;
};
typedef std::shared_ptr<SVDAdapter> SPSVDAdapter;
} // end namespace
} // namespace radixmin
......@@ -2,35 +2,31 @@
#pragma once
namespace radixmin
{
namespace radixmin {
/** A Chain connects a Problem, Adapters, and a Solver and provides the
* solve() function (which repeatedly calls step()). */
template<
class solver_t,
class econ_adapter_t=SVDAdapter,
class icon_adapter_t=ActiveSetAdapter>
class SolverChain
{
template <class solver_t, class econ_adapter_t = SVDAdapter,
class icon_adapter_t = ActiveSetAdapter>
class SolverChain {
protected:
typedef std::vector<SPSolver> VectorSolver;
typedef std::vector<SPSolver> VectorSolver;
SPProblem m_problem;
VectorSolver m_solvers;
SPProblem m_problem;
VectorSolver m_solvers;
public:
SolverChain(SPProblem problem);
SolverChain(SPProblem problem);
void check_consistency() const;
void check_consistency() const;
// get functions
const SPProblem problem() const { return m_problem; }
const VectorSolver solvers() const { return m_solvers; }
// get functions
const SPProblem problem() const { return m_problem; }
const VectorSolver solvers() const { return m_solvers; }
void solve(size_t max_iter=200);
void solve(size_t max_iter = 200);
};
} // end namespace
} // namespace radixmin
#include "chain.i.hh" // the templated implementations
#include "chain.i.hh" // the templated implementations
......@@ -2,82 +2,69 @@
#pragma once
namespace radixmin
{
namespace radixmin {
// solver type, econ_adapter_type, etc
template<class s_t, class e_t, class i_t>
template <class s_t, class e_t, class i_t>
SolverChain<s_t, e_t, i_t>::SolverChain(SPProblem problem)
: m_problem(problem)
{
// build a chain by connecting adapters and a solver to the end
SPProblemInterface last_problem = m_problem;
// each adapter will need a handle on its solver
SPAdapter last_adapter;
: m_problem(problem) {
// build a chain by connecting adapters and a solver to the end
SPProblemInterface last_problem = m_problem;
// each adapter will need a handle on its solver
SPAdapter last_adapter;
// add an icon adapter
if (last_problem->icon().nrows() != 0)
{
auto adapter = std::make_shared<i_t>(last_problem);
if (last_adapter)
{
last_adapter->set_solver(adapter);
}
m_solvers.push_back(adapter);
last_problem = adapter;
last_adapter = adapter;
}
// add an econ adapter
if (last_problem->econ().nrows() != 0 || last_problem->variable_econ())
{
auto adapter = std::make_shared<e_t>(last_problem);
if (last_adapter)
{
last_adapter->set_solver(adapter);
}
m_solvers.push_back(adapter);
last_problem = adapter;
last_adapter = adapter;
// add an icon adapter
if (last_problem->icon().nrows() != 0) {
auto adapter = std::make_shared<i_t>(last_problem);
if (last_adapter) {
last_adapter->set_solver(adapter);
}
// add a solver
auto solver = std::make_shared<s_t>(last_problem);
if (last_adapter)
{
last_adapter->set_solver(solver);
m_solvers.push_back(adapter);
last_problem = adapter;
last_adapter = adapter;
}
// add an econ adapter
if (last_problem->econ().nrows() != 0 || last_problem->variable_econ()) {
auto adapter = std::make_shared<e_t>(last_problem);
if (last_adapter) {
last_adapter->set_solver(adapter);
}
m_solvers.push_back(solver);
m_solvers.push_back(adapter);
last_problem = adapter;
last_adapter = adapter;
}
// add a solver
auto solver = std::make_shared<s_t>(last_problem);
if (last_adapter) {
last_adapter->set_solver(solver);
}
m_solvers.push_back(solver);
}
template<class s_t, class e_t, class i_t>
void SolverChain<s_t, e_t, i_t>::check_consistency() const
{
radix_check(m_problem);
m_problem->check_consistency();
template <class s_t, class e_t, class i_t>
void SolverChain<s_t, e_t, i_t>::check_consistency() const {
radix_check(m_problem);
m_problem->check_consistency();
radix_block(
for (auto solver : m_solvers)
{
radix_check(solver);
solver->check_consistency();
})
radix_block(for (auto solver
: m_solvers) {
radix_check(solver);
solver->check_consistency();
})
}
template<class s_t, class e_t, class i_t>
void SolverChain<s_t, e_t, i_t>::solve(size_t max_iter)
{
this->check_consistency();
SPProblem problem = m_problem;
SPSolver solver = m_solvers[0];
template <class s_t, class e_t, class i_t>
void SolverChain<s_t, e_t, i_t>::solve(size_t max_iter) {
this->check_consistency();
SPProblem problem = m_problem;
SPSolver solver = m_solvers[0];
for (size_t i = 0; i < max_iter; ++i)
{
solver->step();
if (solver->converged())
{
break;
}
for (size_t i = 0; i < max_iter; ++i) {
solver->step();
if (solver->converged()) {
break;
}
}
}
} // end namespace
} // namespace radixmin
#include "radixmin/problems/leastsq.hh"
#include <math.h> // std::exp
#include "radixbug/bug.hh" // radix_check, etc
#include "radixbug/bug.hh" // radix_check, etc
#include <math.h> // std::exp
namespace radixmin {
namespace radixmin
{
FGH LeastSquaresProblem::obj_(const Array &x, int flags) const {
this->check_consistency();
const size_t I = m_data.size();
const size_t J = x.size();
FGH LeastSquaresProblem::obj_(
const Array& x, int flags) const
{
this->check_consistency();
const size_t I = m_data.size();
const size_t J = x.size();
// chi2 G requires model F+G; chi2 H requires model F+G+H
int model_flags = (flags & (F + G + H) ? F : 0) + (flags & (G + H) ? G : 0) +
(flags & (H) ? H : 0);
// chi2 G requires model F+G; chi2 H requires model F+G+H
int model_flags =
(flags & (F+G+H) ? F : 0) +
(flags & ( G+H) ? G : 0) +
(flags & ( H) ? H : 0);
// calculate scaled residuals and their derivatives w/rt x
VectorFGH scaled_residuals = this->model(x, model_flags);
radix_check(scaled_residuals.size() == I);
for (size_t i = 0; i < I; ++i)
{
if (model_flags & F)
{
scaled_residuals[i].f -= m_data[i];
scaled_residuals[i].f /= m_data_stdev[i];
}
if (model_flags & G)
{
radix_check(scaled_residuals[i].g.size() == J);
scaled_residuals[i].g /= m_data_stdev[i];
}
if (model_flags & H && scaled_residuals[i].h.nrows() > 0)
{
radix_check(scaled_residuals[i].h.nrows() == J);
radix_check(scaled_residuals[i].h.ncols() == J);
scaled_residuals[i].h /= m_data_stdev[i];
}
// calculate scaled residuals and their derivatives w/rt x
VectorFGH scaled_residuals = this->model(x, model_flags);
radix_check(scaled_residuals.size() == I);
for (size_t i = 0; i < I; ++i) {
if (model_flags & F) {
scaled_residuals[i].f -= m_data[i];
scaled_residuals[i].f /= m_data_stdev[i];
}
if (model_flags & G) {
radix_check(scaled_residuals[i].g.size() == J);
scaled_residuals[i].g /= m_data_stdev[i];
}
if (model_flags & H && scaled_residuals[i].h.nrows() > 0) {
radix_check(scaled_residuals[i].h.nrows() == J);
radix_check(scaled_residuals[i].h.ncols() == J);
scaled_residuals[i].h /= m_data_stdev[i];
}
}
// calculate chi2 and its derivatives w/rt x
FGH chi2;
if (flags & F)
{
chi2.f = 0.;
for (size_t i = 0; i < I; ++i)
{
chi2.f += scaled_residuals[i].f * scaled_residuals[i].f;
}
// calculate chi2 and its derivatives w/rt x
FGH chi2;
if (flags & F) {
chi2.f = 0.;
for (size_t i = 0; i < I; ++i) {
chi2.f += scaled_residuals[i].f * scaled_residuals[i].f;
}
if (flags & G)
{
chi2.g = Array(J, 0.);
for (size_t i = 0; i < I; ++i)
{
chi2.g += scaled_residuals[i].g * scaled_residuals[i].f * 2.;
}
}
if (flags & G) {
chi2.g = Array(J, 0.);
for (size_t i = 0; i < I; ++i) {
chi2.g += scaled_residuals[i].g * scaled_residuals[i].f * 2.;
}
if (flags & H)
{
chi2.h = Matrix(J, J, 0.);
for (size_t i = 0; i < I; ++i)
{
chi2.h += scaled_residuals[i].g.outer(scaled_residuals[i].g) * 2.;
if (scaled_residuals[i].h.nrows() > 0)
{
chi2.h += scaled_residuals[i].h * scaled_residuals[i].f * 2.;
}
}
}
if (flags & H) {
chi2.h = Matrix(J, J, 0.);
for (size_t i = 0; i < I; ++i) {
chi2.h += scaled_residuals[i].g.outer(scaled_residuals[i].g) * 2.;
if (scaled_residuals[i].h.nrows() > 0) {
chi2.h += scaled_residuals[i].h * scaled_residuals[i].f * 2.;
}
}
return chi2;
}
return chi2;
}
VectorFGH QuadraticLSProblem::model(const Array& x, int flags) const
{
const size_t I = 3;
const size_t J = 3;
VectorFGH QuadraticLSProblem::model(const Array &x, int flags) const {
const size_t I = 3;
const size_t J = 3;
std::vector<Matrix> coeffs(I, Matrix(J, J));
coeffs[0][0][0] = 1.;
coeffs[1][1][1] = 2.;
coeffs[2][2][2] = 3.;
std::vector<Matrix> coeffs(I, Matrix(J, J));
coeffs[0][0][0] = 1.;
coeffs[1][1][1] = 2.;
coeffs[2][2][2] = 3.;
VectorFGH fgh(I);
for (size_t i = 0; i < I; ++i)
{
if (flags & F)
{
fgh[i].f = coeffs[i].dot(x).dot(x);
}
if (flags & G)
{
fgh[i].g = coeffs[i].dot(x) * 2.;
}
if (flags & H)
{
fgh[i].h = coeffs[i] * 2.;
}
VectorFGH fgh(I);
for (size_t i = 0; i < I; ++i) {
if (flags & F) {
fgh[i].f = coeffs[i].dot(x).dot(x);
}
if (flags & G) {
fgh[i].g = coeffs[i].dot(x) * 2.;
}
if (flags & H) {
fgh[i].h = coeffs[i] * 2.;
}
}
return fgh;
return fgh;
}
VectorFGH ExpLSCProblem::model(const Array& x, int flags) const
{
const size_t I = m_data.size();
const size_t J = x.size();
radix_check(I == J);
VectorFGH ExpLSCProblem::model(const Array &x, int flags) const {
const size_t I = m_data.size();
const size_t J = x.size();
radix_check(I == J);
VectorFGH fgh(I);
for (size_t i = 0; i < I; ++i)
{
double f = std::exp(x[i]);
VectorFGH fgh(I);
for (size_t i = 0; i < I; ++i) {
double f = std::exp(x[i]);
if (flags & F)
{
fgh[i].f = f;
}
if (flags & G)
{
fgh[i].g = Array(J, 0.);
fgh[i].g[i] = f;
}
if (flags & H)
{
fgh[i].h = Matrix(J, J, 0.);
fgh[i].h[i][i] = f;
}
if (flags & F) {
fgh[i].f = f;
}
if (flags & G) {
fgh[i].g = Array(J, 0.);
fgh[i].g[i] = f;
}
if (flags & H) {
fgh[i].h = Matrix(J, J, 0.);
fgh[i].h[i][i] = f;
}
}
return fgh;
return fgh;
}
} // end namespace
} // namespace radixmin
#include "radixmin/problems/misc.hh"
namespace radixmin {
namespace radixmin
{
FGH RosenbrockProblem::obj_(
const Array& x, int flags) const
{
this->check_consistency();
// define a single term in the summation
auto calc_term = [flags](double x0, double x1)
{
FGH fgh;