Commit 2a4beb5f authored by Max Lyon's avatar Max Lyon
Browse files

Merge branch 'marinom/REFORM-922-merge-new-sanitization' into 'master'

WIP REFORM-922 merge new sanitization

See merge request !55
parents bdd20c88 e12b50c0
Pipeline #15360 passed with stages
in 4 minutes and 22 seconds
......@@ -134,10 +134,11 @@ protected:
void pre_factorize(NProblemInterface* _quadratic_problem, const SMatrixD& _A, const VectorD& _b)
{
const int n = _quadratic_problem->n_unknowns();
const int m = _A.rows();
const int m = static_cast<int>(_A.rows());
const int nf = n+m;
// get hessian of quadratic problem
SMatrixD H(n,n);
_quadratic_problem->eval_hessian(x_ls_.data(), H);
......@@ -146,21 +147,30 @@ protected:
std::vector< Triplet > trips;
trips.reserve(H.nonZeros() + 2*_A.nonZeros());
// add elements of H
for (int k=0; k<H.outerSize(); ++k)
for (SMatrixD::InnerIterator it(H,k); it; ++it)
trips.push_back(Triplet(it.row(),it.col(),it.value()));
const auto add_triplet = [&trips](
const size_t _i, const size_t _j, const double _v)
{
trips.emplace_back(static_cast<int>(_i), static_cast<int>(_j), _v);
};
// add elements of H
for (int k = 0; k < H.outerSize(); ++k)
{
for (SMatrixD::InnerIterator it(H, k); it; ++it)
add_triplet(it.row(), it.col(), it.value());
}
// add elements of _A
for (int k=0; k<_A.outerSize(); ++k)
for (SMatrixD::InnerIterator it(_A,k); it; ++it)
for (int k = 0; k < _A.outerSize(); ++k)
{
for (SMatrixD::InnerIterator it(_A, k); it; ++it)
{
// insert _A block below
trips.push_back(Triplet(it.row()+n,it.col(),it.value()));
add_triplet(it.row() + n, it.col(), it.value());
// insert _A^T block right
trips.push_back(Triplet(it.col(),it.row()+n,it.value()));
add_triplet(it.col(), it.row() + n, it.value());
}
}
// create KKT matrix
SMatrixD KKT(nf,nf);
......@@ -175,8 +185,8 @@ protected:
// DEB_line(2, "-> re-try with regularized constraints...");
std::cerr << "Eigen::SparseLU reported problem while factoring KKT system: " << lu_solver_.lastErrorMessage() << std::endl;
for( int i=0; i<m; ++i)
trips.push_back(Triplet(n+i,n+i,1e-9));
for (int i = 0; i < m; ++i)
add_triplet(n + i, n + i, 1e-9);
// create KKT matrix
KKT.setFromTriplets(trips.begin(), trips.end());
......
......@@ -37,7 +37,7 @@ LinearConstraint::LinearConstraint(const SVectorNC& _coeffs, const double _b, co
int LinearConstraint::n_unknowns()
{
return coeffs_.innerSize();
return static_cast<int>(coeffs_.innerSize());
}
void LinearConstraint::resize(const std::size_t _n)
......@@ -103,7 +103,8 @@ void LinearConstraint::eval_gradient( const double* _x, SVectorNC& _g )
void LinearConstraint::eval_hessian ( const double* _x, SMatrixNC& _h )
{
_h.clear();
_h.resize(coeffs_.innerSize(), coeffs_.innerSize());
_h.resize(static_cast<unsigned int>(coeffs_.innerSize()),
static_cast<unsigned int>(coeffs_.innerSize()));
}
......
......@@ -46,7 +46,7 @@ solve(NProblemGmmInterface* _problem)
// check for convergence
if( gmm::vect_norm2(g) < eps_)
{
DEB_line(verbosity_, "Newton Solver converged after " << i << " iterations");
DEB_line(2, "Newton Solver converged after " << i << " iterations");
_problem->store_result(P(x));
return true;
}
......@@ -80,7 +80,7 @@ solve(NProblemGmmInterface* _problem)
f = f_new;
improvement = true;
DEB_line(verbosity_, "energy improved to " << f);
DEB_line(6, "energy improved to " << f);
}
}
......@@ -97,7 +97,7 @@ solve(NProblemGmmInterface* _problem)
else
{
_problem->store_result(P(x));
DEB_line(verbosity_, "Newton solver reached max regularization but did not "
DEB_warning(2, "Newton solver reached max regularization but did not "
"converge");
converged_ = false;
return false;
......@@ -105,7 +105,7 @@ solve(NProblemGmmInterface* _problem)
}
}
_problem->store_result(P(x));
DEB_line(verbosity_, "Newton Solver did not converge!!! after iterations.");
DEB_warning(2, "Newton Solver did not converge!!! after iterations.");
converged_ = false;
return false;
......@@ -123,7 +123,7 @@ solve(NProblemGmmInterface* _problem)
int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
const VectorD& _b)
{
DEB_time_func(verbosity_);
DEB_time_func_def;
converged_ = false;
const double KKT_res_eps = 1e-6;
......@@ -136,7 +136,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
// number of constraints
size_t m = _b.size();
DEB_line(verbosity_, "optimize via Newton with " << n << " unknowns and " << m <<
DEB_line(2, "optimize via Newton with " << n << " unknowns and " << m <<
" linear constraints");
// initialize vectors of unknowns
......@@ -191,12 +191,14 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
{
DEB_warning(2, "Numerical issues in KKT system");
DEB_warning_if(!fact_ok, 2, "Factorization not ok");
DEB_warning_if(kkt_res2 > KKT_res_eps, 2, "KKT Residuum too high: " << kkt_res2);
DEB_warning_if(constraint_res2 > max_allowed_constraint_violation2, 2, "Constraint violation too high: " << constraint_res2);
DEB_line_if(
kkt_res2 > KKT_res_eps, 3, "KKT Residuum too high: " << kkt_res2);
DEB_line_if(constraint_res2 > max_allowed_constraint_violation2, 3,
"Constraint violation too high: " << constraint_res2);
// alternate hessian and constraints regularization
if(reg_iters % 2 == 0 || regularize_constraints >= regularize_constraints_limit)
{
DEB_line(verbosity_, "residual ^ 2 " << kkt_res2 << "->regularize hessian");
DEB_line(2, "residual ^ 2 " << kkt_res2 << "->regularize hessian");
if(regularize_hessian == 0.0)
regularize_hessian = 1e-6;
else
......@@ -204,7 +206,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
}
else
{
DEB_line(verbosity_, "residual^2 " << kkt_res2 << " -> regularize constraints");
DEB_line(2, "residual^2 " << kkt_res2 << " -> regularize constraints");
if(regularize_constraints == 0.0)
regularize_constraints = 1e-8;
else
......@@ -249,17 +251,16 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
regularize_constraints_limit = regularize_constraints;
}
DEB_line(verbosity_, "iter: " << iter
<< ", f(x) = " << fx << ", t = " << t << " (tmax=" << t_max << ")"
<< (t < t_max ? " _clamped_" : "")
<< ", eps = [Newton decrement] = " << newton_decrement
<< ", constraint violation prior = " << rhs.tail(m).norm()
<< ", after = " << (_b - _A*x).norm()
<< ", KKT residual^2 = " << kkt_res2);
DEB_line(4,
"iter: " << iter << ", f(x) = " << fx << ", t = " << t
<< " (tmax=" << t_max << ")" << (t < t_max ? " _clamped_" : "")
<< ", eps = [Newton decrement] = " << newton_decrement
<< ", constraint violation prior = " << rhs.tail(m).norm()
<< ", after = " << (_b - _A * x).norm()
<< ", KKT residual^2 = " << kkt_res2);
// converged?
if(newton_decrement < eps_ || std::abs(t) < eps_ls_ )
if (newton_decrement < eps_ || std::abs(t) < eps_ls_)
{
converged_ = true;
break;
......
......@@ -57,8 +57,12 @@ public:
typedef Eigen::Triplet<double> Triplet;
/// Default constructor
NewtonSolver(const double _eps = 1e-6, const double _eps_line_search = 1e-9, const int _max_iters = 200, const double _alpha_ls=0.2, const double _beta_ls = 0.6)
: eps_(_eps), eps_ls_(_eps_line_search), max_iters_(_max_iters), alpha_ls_(_alpha_ls), beta_ls_(_beta_ls), verbosity_(2), solver_type_(LS_EigenLU), constant_hessian_structure_(false)
NewtonSolver(const double _eps = 1e-6, const double _eps_line_search = 1e-9,
const int _max_iters = 200, const double _alpha_ls = 0.2,
const double _beta_ls = 0.6)
: eps_(_eps), eps_ls_(_eps_line_search), max_iters_(_max_iters),
alpha_ls_(_alpha_ls), beta_ls_(_beta_ls), solver_type_(LS_EigenLU),
constant_hessian_structure_(false)
{
//#if COMISO_SUITESPARSE_AVAILABLE
// solver_type_ = LS_Umfpack;
......@@ -84,12 +88,8 @@ public:
solver_type_ = _st;
}
// set verbosity level of the solver. Lower numbers are more verbose
void set_verbosity(int _verbosity) { verbosity_ = _verbosity; }
bool converged() { return converged_; }
protected:
bool factorize(NProblemInterface* _problem, const SMatrixD& _A,
......@@ -138,7 +138,6 @@ private:
int max_iters_;
double alpha_ls_;
double beta_ls_;
int verbosity_;
VectorD x_ls_;
......
......@@ -4,6 +4,7 @@
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <CoMISo/Utils/CoMISoError.hh>
#include <Base/Debug/DebOut.hh>
#include <vector>
namespace COMISO {
......@@ -159,8 +160,8 @@ double ExactConstraintSatisfaction::get_delta()
void ExactConstraintSatisfaction::IREF_Gaussian(SP_Matrix_R& A, Eigen::VectorXi& b, const Eigen::VectorXd& x){
number_pivots_ = 0;
int rows = A.rows(); //number of rows
int cols = A.cols(); //number of columns
int rows = static_cast<int>(A.rows()); //number of rows
int cols = static_cast<int>(A.cols()); //number of columns
int col_index = 0; //save the next column where we search for a pivot
for (int k = 0; k < rows; k++)
......@@ -208,7 +209,7 @@ void ExactConstraintSatisfaction::IREF_Gaussian(SP_Matrix_R& A, Eigen::VectorXi&
eliminate_row(A, b, i, k, col_p);
}
}
A.prune(0.0 , 0);
A.prune(0, 0);
}
void ExactConstraintSatisfaction::IRREF_Jordan(SP_Matrix_R& A, Eigen::VectorXi& b)
......@@ -234,14 +235,14 @@ void ExactConstraintSatisfaction::IRREF_Jordan(SP_Matrix_R& A, Eigen::VectorXi&
break;
if (it.value() == 0)
continue;
eliminate_row(A, b, it.row(), k, pivot_col);
eliminate_row(A, b, static_cast<int>(it.row()), k, pivot_col);
}
}
// A.makeCompressed();
// A.finalize();
A.prune(0.0, 0);
A.prune(0, 0);
}
//-------------------------------------Evaluation--------------------------------------------------------------------
......@@ -349,7 +350,7 @@ void ExactConstraintSatisfaction::evaluate(const ExactConstraintSatisfaction::SP
SP_Matrix_C A = _A; //change the matrix type to allow easier iteration
int cols = A.cols();
int cols = static_cast<int>(A.cols());
for(int k = cols -1; k >= 0; k--)
{
auto pivot_row = get_pivot_row_new(A, _A, k);
......@@ -377,7 +378,7 @@ void ExactConstraintSatisfaction::evaluate(const ExactConstraintSatisfaction::SP
int ExactConstraintSatisfaction::get_pivot_col_student(SP_Matrix_R& _A, int k, int& col_index)
{
int pivot_row = -1;
int cols = _A.cols();
int cols = static_cast<int>(_A.cols());
for(; col_index < cols; col_index++)
{
......@@ -402,8 +403,8 @@ int ExactConstraintSatisfaction::get_pivot_col_student(SP_Matrix_R& _A, int k, i
int ExactConstraintSatisfaction::get_pivot_col_new(ExactConstraintSatisfaction::SP_Matrix_R& _A, int k, int& col_index)
{
int cols = _A.cols();
int rows = _A.rows();
int cols = static_cast<int>(_A.cols());
int rows = static_cast<int>(_A.rows());
for(; col_index < cols; col_index++)
{
for (int i = k; i < rows; ++i)
......@@ -490,7 +491,7 @@ std::vector<std::pair<int, double> > ExactConstraintSatisfaction::get_dot_produc
{
DEB_enter_func;
std::vector<std::pair<int, double>> S;
int cols = A.cols();
int cols = static_cast<int>(A.cols());
for(int i = k+1; i < cols; i++)
{ //construct the list S to do the dot Product
std::pair<int, double> pair;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment