Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/micpsolver/micpsolver_min.inl
64b3616eea48reactivetransport

Last Change
Raw File

micpsolver_min.inl
View Options

/*-------------------------------------------------------
- Module : micpsolver
- File : micpsolver_min.inl
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#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);
get_perf().current_residual = m_residuals.norm();
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);
sanitize(x);
int termcode = linesearch(update, x);
get_perf().current_update = update.norm();
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_second_cond; // 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_second_cond:
// third condition : is the descent condition respected
update = solver.solve(-reduced_residual);
m_grad_phicck = reduced_jacobian.transpose()*reduced_residual;
const double descent_cond = m_grad_phicck.dot(update);
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>::sanitize(Eigen::VectorXd& x)
{
if (x(0) <=0) x(0) = 1;
for (int dof=get_neq_free(); dof<get_neq(); ++dof)
{
if (x(dof) <=0) {x(dof) =0; DEBUG<< "sanitize dof : " << 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<p.rows(); ++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
//reformulate_result(x, p);
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 << "cnt : " << cnt << " - lambda : " << lambda;
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)
{
lambda = get_program()->max_lambda(x, p);
xp = x + lambda*p;
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
DEBUG << "not disc : " << - init_slope/(2*b);
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;
}
DEBUG << "lambdatmp : " << lambdatmp;
lambda_prev = lambda;
fcp_prev = fcp;
if (not std::isfinite(lambdatmp))
{
lambda = get_program()->max_lambda(x, p);
xp = x + lambda*p;
retcode = 1;
break;
} else if ((lambdatmp < 0.1*lambda)) {
lambda = 0.1 * lambda;
} else {
lambda = lambdatmp;
}
DEBUG << "lambda end : " << lambda;
}
xp = x + lambda*p;
//sanitize(xp);
DEBUG << "xp : " << std::endl << xp;
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