Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80432945
micpsolver_min.inl
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Aug 31, 17:04
Size
11 KB
Mime Type
text/x-c++
Expires
Mon, Sep 2, 17:04 (2 d)
Engine
blob
Format
Raw Data
Handle
20319220
Attached To
rSPECMICP SpecMiCP / ReactMiCP
micpsolver_min.inl
View Options
/*-------------------------------------------------------
- Module : micpsolver
- File : micpsolver_min.inl
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "micpsolver_min.hpp" // for syntaxic coloration
#include "utils/log.hpp"
#include "ncp_function.hpp"
#include "estimate_cond_number.hpp"
namespace specmicp {
namespace micpsolver {
template <class program_t>
MiCPSolverReturnCode MiCPSolverMin<program_t>::solve(Eigen::VectorXd& x)
{
int cnt = 0;
MiCPSolverReturnCode retcode = MiCPSolverReturnCode::NotConvergedYet;
Eigen::VectorXd update(Eigen::VectorXd::Zero(get_neq()));
setup_residuals(x);
while (retcode == MiCPSolverReturnCode::NotConvergedYet)
{
DEBUG << "Iteration : " << cnt;
DEBUG << "Solution : \n" << x;
get_program()->hook_start_iteration(x, m_residuals.norm());
setup_residuals(x);
SPAM << "Residuals : \n ----- \n" << m_residuals << "\n ----- \n";
retcode = base::check_convergence(cnt, update, x, m_residuals);
get_perf().return_code = retcode;
if (retcode != MiCPSolverReturnCode::NotConvergedYet) break;
++cnt;
get_perf().max_taken = false;
setup_jacobian(x);
search_direction_calculation(x, update);
int termcode = linesearch(update, x);
DEBUG << "Return LineSearch : " << termcode;
base::projection(x);
get_perf().nb_iterations = cnt;
}
return retcode;
}
template <class program_t>
MiCPSolverReturnCode MiCPSolverMin<program_t>::search_direction_calculation(
Eigen::VectorXd& x,
Eigen::VectorXd& update)
{
Eigen::MatrixXd reduced_jacobian;
Eigen::VectorXd reduced_residual;
reduce_system(x, reduced_jacobian, reduced_residual);
DEBUG << reduced_jacobian;
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> solver;
m_gradient_step_taken = false;
solver.compute(reduced_jacobian);
get_perf().nb_factorization += 1;
{ // first condition : is the factorization ok ?
if (solver.info() != Eigen::Success or not solver.isInvertible())
{
DEBUG << "Solver.info : " << solver.info() << " - is invertible : " << solver.isInvertible();
m_gradient_step_taken = true;
goto after_solution; // jump directly to the gradient step
}
}
{ // second condition : is the condition number ok ?
const double cond = estimate_condition_number(solver.matrixR().triangularView<Eigen::Upper>());
DEBUG << "Condition number : " << cond;
if (cond > get_options().condition_limit)
{
m_gradient_step_taken = true;
}
}
after_solution:
// third condition : is the descent condition respected
update = solver.solve(-reduced_residual);
reformulate_result(x, update);
base::reformulate_jacobian_cck(x, m_residuals, m_jacobian);
reformulate_residuals_cck_inplace(x, m_residuals);
m_grad_phicck = m_jacobian.transpose()*m_residuals;
//m_grad_phicck(0) = 0;
if (not m_gradient_step_taken)
{
m_newton_length = base::is_step_too_long(update);
// we compute the descent condition
const double descent_cond = m_grad_phicck.dot(update);
const double norm_update = update.norm();
DEBUG << "grad(phi).dot(update) = " << descent_cond << " compared to : " << (
- get_options().factor_descent_condition * std::pow(norm_update, get_options().power_descent_condition));
if (descent_cond > - get_options().factor_descent_condition * std::pow(norm_update, get_options().power_descent_condition))
{
m_gradient_step_taken = true;
}
}
if (m_gradient_step_taken)
{
INFO << "Full gradient step taken !";
update = - m_grad_phicck;
m_newton_length = base::is_step_too_long(update);
}
return MiCPSolverReturnCode::NotConvergedYet;
}
template <class program_t>
int MiCPSolverMin<program_t>::reduce_system(const Eigen::VectorXd& x,
Eigen::MatrixXd& reduced_jacobian,
Eigen::VectorXd& reduced_residual)
{
reduced_jacobian.resizeLike(m_jacobian); // memory is cheap, we will resize at the end
reduced_jacobian.setZero();
reduced_residual.resizeLike(m_residuals);
// copy identical information
int ideq_reduced = get_neq_free();
reduced_jacobian.block(0, 0, get_neq_free(), get_neq_free()) = m_jacobian.block(0, 0, get_neq_free(), get_neq_free());
// select active degree of freedom
Eigen::VectorXd to_remove(get_neq()-get_neq_free());
for (int dof=get_neq_free(); dof<get_neq(); ++dof)
{
if (x(dof) >= m_residuals(dof))
{
DEBUG << "Mineral to precipitate : " << dof;
reduced_residual(ideq_reduced) = m_residuals(dof);
reduced_jacobian.block(ideq_reduced, 0, 1, get_neq_free()) = m_jacobian.block(dof, 0, 1, get_neq_free());
reduced_jacobian.block(0, ideq_reduced, get_neq_free(), 1) = m_jacobian.block(0, dof, get_neq_free(), 1);
to_remove(dof-get_neq_free()) = 0;
++ideq_reduced;
}
else
{
to_remove(dof-get_neq_free()) = x(dof);
}
}
reduced_residual.block(0, 0, get_neq_free(), 1) -=
m_jacobian.block(0, get_neq_free(), get_neq_free(), get_neq()-get_neq_free())*to_remove;
reduced_jacobian.conservativeResize(ideq_reduced, ideq_reduced);
reduced_residual.conservativeResize(ideq_reduced);
DEBUG << "ideq reduced : " << ideq_reduced;
return ideq_reduced;
}
template <class program_t>
void MiCPSolverMin<program_t>::reformulate_result(const Eigen::VectorXd& x,
Eigen::VectorXd& update)
{
update.conservativeResizeLike(x);
int tot_to_keep = 0;
for (int dof=get_neq_free(); dof<get_neq(); ++dof)
{
if (x(dof) >= m_residuals(dof))
++tot_to_keep;
}
int kept_dof = 1;
for (int dof=get_neq()-1; dof>=get_neq_free(); --dof)
{ // we go backwards to avoid extra copies
if (x(dof) >= m_residuals(dof))
{
update(dof) = update(get_neq_free()+(tot_to_keep-kept_dof));
++kept_dof;
}
else
{
update(dof) = -x(dof);
}
}
}
template <class program_t>
void MiCPSolverMin<program_t>::reformulate_residuals_cck_inplace(const Eigen::VectorXd& x,
Eigen::VectorXd& residuals)
{
for (int i = get_neq_free(); i<get_neq(); ++i)
{
residuals(i) = penalized_fisher_burmeister(x(i), residuals(i), get_options().penalization_factor);
}
}
template <class program_t>
int MiCPSolverMin<program_t>::linesearch(Eigen::VectorXd& p, Eigen::VectorXd& x)
{
// Reference Algo A6.3.1 : Dennis and Schnabel (1983)
DEBUG << "Linesearch";
Eigen::VectorXd xp(get_neq());
Eigen::VectorXd new_res(get_neq());
double fcp;
get_perf().max_taken = false;
int retcode = 2;
const double alpha = get_options().factor_descent_condition;
double newtlen = m_newton_length;
//double newtlen = p.norm();
double init_slope = m_grad_phicck.dot(p);
double rellength = std::abs(p(0));
for (int i=1; i<get_neq(); ++i)
{
rellength = std::max(rellength, std::abs(p(i)));
}
double minlambda = get_options().steptol / rellength;
double lambda = get_program()->max_lambda(x, p);
DEBUG << "Initial lambda : " << lambda;
double lambda_prev = lambda;
// non monotone linesearch
//
// - reference : Munson et al. (2001)
// ------------------------------------
double merit_value = 0.5*m_residuals.squaredNorm();
// // new residual
xp = x + lambda*p;
DEBUG << "update \n" << p <<std::endl;
base::compute_residuals(xp, new_res);
reformulate_residuals_cck_inplace(xp, new_res);
fcp = 0.5*new_res.squaredNorm();
// Skip linesearch if enough progress is done
if (fcp < get_options().coeff_accept_newton_step *merit_value)
{
if (m_max_merits.size() > 0) m_max_merits[m_max_merits.size()-1] = merit_value;
else m_max_merits.push_back(merit_value);
x = xp;
return 0;
}
DEBUG << "Merit value : " << merit_value;
double mmax = merit_value;
if (m_max_merits.size() > 0)
{
mmax = m_max_merits[m_max_merits.size()-1];
}
if (m_max_merits.size() < 4)
{
m_max_merits.push_back(merit_value);
if (merit_value < mmax) merit_value = (3*merit_value + mmax)/4;
}
else if (merit_value < mmax)
{
m_max_merits[3] = merit_value;
merit_value = mmax;
}
if (m_gradient_step_taken)
{
merit_value *= 100;
}
DEBUG << "Merit value used : " << merit_value;
double fc = merit_value;
double fcp_prev;
int cnt = 0;
do
{
DEBUG << "fcp : " << fcp << "\n fc+alin : " << fc+alpha*lambda*init_slope << " # fc : " << fc << std::endl;
if (fcp <= fc + alpha*lambda*init_slope)
{
retcode = 0;
if (lambda ==1 and (newtlen > 0.99 * get_options().maxstep)) {
get_perf().max_taken = true;
}
break;
}
else if (lambda < minlambda)
{
retcode = 1;
break;
}
else
{
double lambdatmp;
if (cnt == 0) { // only a quadratic at the first
lambdatmp = - init_slope / (2*(fcp - fc -init_slope));
}
else
{
const double factor = 1 /(lambda - lambda_prev);
const double x1 = fcp - fc - lambda*init_slope;
const double x2 = fcp_prev - fc - lambda_prev*init_slope;
const double a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev));
const double b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev));
if (a == 0)
{ // cubic interpolation is in fact a quadratic interpolation
lambdatmp = - init_slope/(2*b);
}
else
{
const double disc = b*b-3*a*init_slope;
lambdatmp = (-b+std::sqrt(disc))/(3*a);
}
if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda;
}
lambda_prev = lambda;
fcp_prev = fcp;
if (lambdatmp < 0.1*lambda) {
lambda = 0.1 * lambda;
} else {
lambda = lambdatmp;
}
}
xp = x + lambda*p;
base::compute_residuals(xp, new_res);
reformulate_residuals_cck_inplace(xp, new_res);
fcp = 0.5*new_res.squaredNorm();
++cnt;
} while(retcode == 2 and cnt < 100);
DEBUG << "Lambda : " << lambda;
if (cnt == 100)
{
ERROR << "Too much linesearch iterations ! We stop";
}
x = xp;
p = lambda*p;
return retcode;
}
} // end namespace micpsolver
} // end namespace specmicp
Event Timeline
Log In to Comment