Homec4science
Diffusion SpecMiCP / ReactMiCP (chloride) Browse

specmicp/src/micpsolver/micpsolver_min.inl
98612e6c3e19chloride

Last Change
Raw File

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
c4science · Help