Skip to content

Commit

Permalink
refactored LS
Browse files Browse the repository at this point in the history
  • Loading branch information
@teseoch
teseoch committed Oct 4, 2023
1 parent 6600af7 commit e82e514
Show file tree
Hide file tree
Showing 9 changed files with 185 additions and 306 deletions.
121 changes: 21 additions & 100 deletions src/polysolve/nonlinear/line_search/Armijo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,91 +11,52 @@ namespace polysolve::nonlinear::line_search
Armijo::Armijo(spdlog::logger &logger)
: Superclass(logger)
{
this->min_step_size = 1e-7;
this->max_step_size_iter = 20; // std::numeric_limits<int>::max();
min_step_size = 1e-7;
max_step_size_iter = 20; // std::numeric_limits<int>::max();

default_init_step_size = 1.0;
step_ratio = 0.5;
}

double Armijo::line_search(
double Armijo::compute_descent_step_size(
const TVector &x,
const TVector &searchDir,
Problem &objFunc)
Problem &objFunc,
const bool use_grad_norm,
const double old_energy_in,
const double starting_step_size)
{
TVector grad(x.rows());
double f_in, alpha;
double alpha_init = default_alpha_init;

{
POLYSOLVE_SCOPED_TIMER("LS begin");

this->cur_iter = 0;

f_in = objFunc.value(x);
if (std::isnan(f_in))
{
m_logger.error("Original energy in line search is nan!");
return std::nan("");
}

objFunc.gradient(x, grad);

// TODO: removed feature
// alpha_init = std::min(objFunc.heuristic_max_step(searchDir), alpha_init);
}

{
POLYSOLVE_SCOPED_TIMER("LS compute finite energy step size", this->checking_for_nan_inf_time);
alpha_init = this->compute_nan_free_step_size(x, searchDir, objFunc, alpha_init, tau);
if (std::isnan(alpha_init))
return std::nan("");
}

const double nan_free_step_size = alpha_init;

{
POLYSOLVE_SCOPED_TIMER("CCD broad-phase", this->broad_phase_ccd_time);
const TVector x1 = x + alpha_init * searchDir;
objFunc.line_search_begin(x, x1);
}

{
POLYSOLVE_SCOPED_TIMER("CCD narrow-phase", this->ccd_time);
alpha = this->compute_collision_free_step_size(x, searchDir, objFunc, alpha_init);
if (std::isnan(alpha))
return std::nan("");
}

const double collision_free_step_size = alpha;
double f_in = old_energy_in;
double alpha = starting_step_size;

double f;
bool valid;
{
POLYSOLVE_SCOPED_TIMER("energy min in LS", this->classical_line_search_time);
POLYSOLVE_SCOPED_TIMER("energy min in LS", classical_line_search_time);

TVector x1 = x + alpha * searchDir;
{
POLYSOLVE_SCOPED_TIMER("constraint set update in LS", this->constraint_set_update_time);
POLYSOLVE_SCOPED_TIMER("constraint set update in LS", constraint_set_update_time);
objFunc.solution_changed(x1);
}

objFunc.gradient(x, grad);
const bool use_grad_norm = grad.norm() < this->use_grad_norm_tol;
if (use_grad_norm)
f_in = grad.squaredNorm();

f = use_grad_norm ? grad.squaredNorm() : objFunc.value(x1);
const double Cache = c * grad.dot(searchDir);
double f = use_grad_norm ? grad.squaredNorm() : objFunc.value(x1);
const double cache = c * grad.dot(searchDir);
valid = objFunc.is_step_valid(x, x1);

// max_step_size should return a collision free step
// assert(objFunc.is_step_collision_free(x, x1));

while ((std::isinf(f) || std::isnan(f) || f > f_in + alpha * Cache || !valid) && alpha > this->min_step_size && this->cur_iter <= this->max_step_size_iter)
while ((std::isinf(f) || std::isnan(f) || f > f_in + alpha * cache || !valid) && alpha > min_step_size && cur_iter <= max_step_size_iter)
{
alpha *= tau;
alpha *= step_ratio;
x1 = x + alpha * searchDir;

{
POLYSOLVE_SCOPED_TIMER("constraint set update in LS", this->constraint_set_update_time);
POLYSOLVE_SCOPED_TIMER("constraint set update in LS", constraint_set_update_time);
objFunc.solution_changed(x1);
}

Expand All @@ -109,52 +70,12 @@ namespace polysolve::nonlinear::line_search

valid = objFunc.is_step_valid(x, x1);

// max_step_size should return a collision free step
// assert(objFunc.is_step_collision_free(x, x1));

m_logger.trace("ls it: {} f: {} (f_in + alpha * Cache): {} invalid: {} ", this->cur_iter, f, f_in + alpha * Cache, !valid);

this->cur_iter++;
}
}

const double descent_step_size = alpha;
m_logger.trace("ls it: {} f: {} (f_in + alpha * Cache): {} invalid: {} ", cur_iter, f, f_in + alpha * cache, !valid);

if (this->cur_iter >= this->max_step_size_iter || alpha <= this->min_step_size)
{
{
POLYSOLVE_SCOPED_TIMER("LS end");
objFunc.line_search_end();
cur_iter++;
}

m_logger.warn(
"Line search failed to find descent step (f(x)={:g} f(x+αΔx)={:g} α_CCD={:g} α={:g}, ||Δx||={:g} is_step_valid={} iter={:d})",
f_in, f, default_alpha_init, alpha, searchDir.norm(),
valid ? "true" : "false", this->cur_iter);
return std::nan("");
}

// #ifndef NDEBUG
// // -------------
// // CCD safeguard
// // -------------
// {
// POLYSOLVE_SCOPED_TIMER("safeguard in LS");
// alpha = this->compute_debug_collision_free_step_size(x, searchDir, objFunc, alpha, tau);
// }

// const double debug_collision_free_step_size = alpha;
// #endif

{
POLYSOLVE_SCOPED_TIMER("LS end");
objFunc.line_search_end();
}

m_logger.debug(
"Line search finished (nan_free_step_size={} collision_free_step_size={} descent_step_size={} final_step_size={})",
nan_free_step_size, collision_free_step_size, descent_step_size, alpha);

return alpha;
}

Expand Down
14 changes: 8 additions & 6 deletions src/polysolve/nonlinear/line_search/Armijo.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,16 @@ namespace polysolve::nonlinear::line_search

Armijo(spdlog::logger &logger);

double line_search(
protected:
double compute_descent_step_size(
const TVector &x,
const TVector &searchDir,
Problem &objFunc) override;
const TVector &delta_x,
Problem &objFunc,
const bool use_grad_norm,
const double old_energy_in,
const double starting_step_size) override;

protected:
const double default_alpha_init = 1.0;
private:
const double c = 0.5;
const double tau = 0.5;
};
} // namespace polysolve::nonlinear::line_search
131 changes: 6 additions & 125 deletions src/polysolve/nonlinear/line_search/Backtracking.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,132 +14,29 @@ namespace polysolve::nonlinear::line_search
{
this->min_step_size = 0;
this->max_step_size_iter = 100; // std::numeric_limits<int>::max();
}

double Backtracking::line_search(
const TVector &x,
const TVector &delta_x,
Problem &objFunc)
{
// ----------------
// Begin linesearch
// ----------------

double old_energy, step_size;
{
POLYSOLVE_SCOPED_TIMER("LS begin");

this->cur_iter = 0;

old_energy = objFunc.value(x);
if (std::isnan(old_energy))
{
m_logger.error("Original energy in line search is nan!");
return std::nan("");
}

step_size = 1;
// TODO: removed feature
// objFunc.heuristic_max_step(delta_x);
}

// ----------------------------
// Find finite energy step size
// ----------------------------

{
POLYSOLVE_SCOPED_TIMER("LS compute finite energy step size", this->checking_for_nan_inf_time);
step_size = this->compute_nan_free_step_size(x, delta_x, objFunc, step_size, 0.5);
if (std::isnan(step_size))
return std::nan("");
}

const double nan_free_step_size = step_size;

// -----------------------------
// Find collision-free step size
// -----------------------------

{
POLYSOLVE_SCOPED_TIMER("Line Search Begin - CCD broad-phase", this->broad_phase_ccd_time);
TVector new_x = x + step_size * delta_x;
objFunc.line_search_begin(x, new_x);
}

{
POLYSOLVE_SCOPED_TIMER("CCD narrow-phase", this->ccd_time);
m_logger.trace("Performing narrow-phase CCD");
step_size = this->compute_collision_free_step_size(x, delta_x, objFunc, step_size);
if (std::isnan(step_size))
return std::nan("");
}

const double collision_free_step_size = step_size;

// ----------------------
// Find descent step size
// ----------------------

{
POLYSOLVE_SCOPED_TIMER("energy min in LS", this->classical_line_search_time);
step_size = compute_descent_step_size(x, delta_x, objFunc, old_energy, step_size);
if (std::isnan(step_size))
{
// Superclass::save_sampled_values("failed-line-search-values.csv", x, delta_x, objFunc);
return std::nan("");
}
}

const double descent_step_size = step_size;

// #ifndef NDEBUG
// // -------------
// // CCD safeguard
// // -------------

// {
// POLYSOLVE_SCOPED_TIMER("safeguard in LS");
// step_size = this->compute_debug_collision_free_step_size(x, delta_x, objFunc, step_size, 0.5);
// }

// const double debug_collision_free_step_size = step_size;
// #endif

{
POLYSOLVE_SCOPED_TIMER("LS end");
objFunc.line_search_end();
}

m_logger.debug(
"Line search finished (nan_free_step_size={} collision_free_step_size={} descent_step_size={} final_step_size={})",
nan_free_step_size, collision_free_step_size, descent_step_size, step_size);

return step_size;
default_init_step_size = 1.0;
step_ratio = 0.5;
}

double Backtracking::compute_descent_step_size(
const TVector &x,
const TVector &delta_x,
Problem &objFunc,
const double old_energy_in,
const bool use_grad_norm,
const double old_energy,
const double starting_step_size)
{
double step_size = starting_step_size;

TVector grad(x.rows());
objFunc.gradient(x, grad);
const bool use_grad_norm = grad.norm() < this->use_grad_norm_tol;
if (grad.norm() < 1e-30)
return 1;

const double old_energy = use_grad_norm ? grad.squaredNorm() : old_energy_in;

// Find step that reduces the energy
double cur_energy = std::nan("");
bool is_step_valid = false;
while (step_size > this->min_step_size && this->cur_iter < this->max_step_size_iter)
while (step_size > min_step_size && cur_iter < max_step_size_iter)
{
this->iterations++;
iterations++;

TVector new_x = x + step_size * delta_x;

Expand Down Expand Up @@ -169,7 +66,6 @@ namespace polysolve::nonlinear::line_search

m_logger.trace("ls it: {} delta: {} invalid: {} ", this->cur_iter, (cur_energy - old_energy), !is_step_valid);

// if (!std::isfinite(cur_energy) || (cur_energy >= old_energy && fabs(cur_energy - old_energy) > 1e-12) || !is_step_valid)
if (!std::isfinite(cur_energy) || cur_energy >= old_energy || !is_step_valid)
{
step_size /= 2.0;
Expand All @@ -183,21 +79,6 @@ namespace polysolve::nonlinear::line_search
this->cur_iter++;
}

if (this->cur_iter >= this->max_step_size_iter || step_size <= this->min_step_size)
{
m_logger.warn(
"Line search failed to find descent step (f(x)={:g} f(x+αΔx)={:g} α_CCD={:g} α={:g}, ||Δx||={:g} is_step_valid={} use_grad_norm={} iter={:d})",
old_energy, cur_energy, starting_step_size, step_size, delta_x.norm(),
is_step_valid, use_grad_norm, this->cur_iter);
objFunc.solution_changed(x);
#ifndef NDEBUG
// tolerance for rounding error due to multithreading
assert(abs(old_energy_in - objFunc.value(x)) < 1e-15);
#endif
objFunc.line_search_end();
return std::nan("");
}

return step_size;
}
} // namespace polysolve::nonlinear::line_search
8 changes: 2 additions & 6 deletions src/polysolve/nonlinear/line_search/Backtracking.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,17 +13,13 @@ namespace polysolve::nonlinear::line_search

Backtracking(spdlog::logger &logger);

double line_search(
const TVector &x,
const TVector &delta_x,
Problem &objFunc) override;

protected:
double compute_descent_step_size(
const TVector &x,
const TVector &delta_x,
Problem &objFunc,
const bool use_grad_norm,
const double old_energy_in,
const double starting_step_size = 1);
const double starting_step_size) override;
};
} // namespace polysolve::nonlinear::line_search
13 changes: 9 additions & 4 deletions src/polysolve/nonlinear/line_search/CppOptArmijo.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,12 +18,17 @@ namespace polysolve::nonlinear::line_search
{
}

double line_search(
protected:
double compute_descent_step_size(
const TVector &x,
const TVector &searchDir,
Problem &objFunc) override
const TVector &delta_x,
Problem &objFunc,
const bool,
const double,
const double starting_step_size) override
{
return cppoptlib::Armijo<Problem, 1>::linesearch(x, searchDir, objFunc);
const double tmp = cppoptlib::Armijo<Problem, 1>::linesearch(x, delta_x, objFunc);
return std::min(starting_step_size, tmp);
}
};
} // namespace polysolve::nonlinear::line_search
Loading

0 comments on commit e82e514

Please sign in to comment.