/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* Geophysical Computational Tools & Library (GCTL)
*
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
*
* GCTL is distributed under a dual licensing scheme. You can redistribute
* it and/or modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation, either version 2
* of the License, or (at your option) any later version. You should have
* received a copy of the GNU Lesser General Public License along with this
* program. If not, see .
*
* If the terms and conditions of the LGPL v.2. would prevent you from using
* the GCTL, please consider the option to obtain a commercial license for a
* fee. These licenses are offered by the GCTL's original author. As a rule,
* licenses are provided "as-is", unlimited in time for a one time fee. Please
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
* to include some description of your company and the realm of its activities.
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#ifndef _GCTL_LCG_H
#define _GCTL_LCG_H
#include "gctl/core.h"
#include "gctl/maths.h"
#include "gctl/algorithm.h"
#include "gctl_optimization_config.h"
#ifdef GCTL_OPTIMIZATION_TOML
#include "toml.hpp"
#endif // GCTL_OPTIMIZATION_TOML
#if defined _WINDOWS || __WIN32__
#include "windows.h"
#endif // _WINDOWS || __WIN32__
namespace gctl
{
/**
* @brief Types of method that could be recognized by the lcg_solver() function.
*/
enum lcg_solver_type
{
/**
* Conjugate gradient method.
*/
LCG_CG,
/**
* Preconditioned conjugate gradient method.
*/
LCG_PCG,
/**
* Conjugate gradient squared method.
*/
LCG_CGS,
/**
* Biconjugate gradient method.
*/
LCG_BICGSTAB,
/**
* Biconjugate gradient method with restart.
*/
LCG_BICGSTAB2,
/**
* Conjugate gradient method with projected gradient for inequality constraints.
* This algorithm comes without non-monotonic linear search for the step length.
*/
LCG_PG,
/**
* Conjugate gradient method with spectral projected gradient for inequality constraints.
* This algorithm comes with non-monotonic linear search for the step length.
*/
LCG_SPG,
};
/**
* @brief return value of the lcg_solver() function
*/
enum lcg_return_code
{
LCG_SUCCESS = 0, ///< The solver function terminated successfully.
LCG_CONVERGENCE = 0, ///< The iteration reached convergence.
LCG_STOP, ///< The iteration is stopped by the monitoring function.
LCG_ALREADY_OPTIMIZIED, ///< The initial solution is already optimized.
// A negative number means a error
LCG_UNKNOWN_ERROR = -1024, ///< Unknown error.
LCG_INVILAD_VARIABLE_SIZE, ///< The variable size is negative
LCG_INVILAD_MAX_ITERATIONS, ///< The maximal iteration times is negative.
LCG_INVILAD_EPSILON, ///< The epsilon is negative.
LCG_INVILAD_RESTART_EPSILON, ///< The restart epsilon is negative.
LCG_REACHED_MAX_ITERATIONS, ///< Iteration reached maximal limit.
LCG_NULL_PRECONDITION_MATRIX, ///< Null precondition matrix.
LCG_NAN_VALUE, ///< Nan value.
LCG_INVALID_POINTER, ///< Invalid pointer.
LCG_INVALID_LAMBDA, ///< Invalid range for lambda.
LCG_INVALID_SIGMA, ///< Invalid range for sigma.
LCG_INVALID_BETA, ///< Invalid range for beta.
LCG_INVALID_MAXIM, ///< Invalid range for maxi_m.
LCG_SIZE_NOT_MATCH, ///< Sizes of m and B do not match
};
/**
* @brief Message type of the LCG algorithms.
*
*/
enum lcg_message_type
{
LCG_THROW, ///< throw error only
LCG_ERROR, ///< display error only
LCG_SOLUTION, ///< display info for evry solution
LCG_ITERATION, ///< display info for every iteration
};
/**
* @brief Parameters of the conjugate gradient methods.
*/
struct lcg_para
{
/**
* Maximal iteration times. The process will continue till the convergence is met
* if this option is set to zero (default).
*/
int max_iterations;
/**
* Epsilon for convergence test.
* This parameter determines the accuracy with which the solution is to be
* found. A minimization terminates when ||g||/max(||g0||, 1.0) <= epsilon or
* sqrt(||g||)/N <= epsilon for the lcg_solver() function, where ||.|| denotes
* the Euclidean (L2) norm. The default value of epsilon is 1e-8.
*/
double epsilon;
/**
* Whether to use absolute mean differences (AMD) between |Ax - B| to evaluate the process.
* The default value is false which means the gradient based evaluating method is used.
* The AMD based method will be used if this variable is set to true. This parameter is only
* applied to the non-constrained methods.
*/
int abs_diff;
/**
* Restart epsilon for the LCG_BICGSTAB2 algorithm. The default value is 1e-6
*/
double restart_epsilon;
/**
* Initial step length for the project gradient method. The default is 1.0
*/
double step;
/**
* multiplier for updating solutions with the spectral projected gradient method. The range of
* this variable is (0, 1). The default is given as 0.95
*/
double sigma;
/**
* descending ratio for conducting the non-monotonic linear search. The range of
* this variable is (0, 1). The default is given as 0.9
*/
double beta;
/**
* The maximal record times of the objective values for the SPG method. The method use the
* objective values from the most recent maxi_m times to preform the non-monotonic linear search.
* The default value is 10.
*/
int maxi_m;
};
/**
* @brief This abstract class implements conjugate gradient algorithms for solving
* a linear system like Ax = B where A is a N*N matrix. As the actual variable
* we need is the product of 'Ax', the kernel matrix 'A' is not declared within
* the class definition. Instead, a pure virtual function is declared as the
* callback interface for calculating the product of 'Ax' as 'void LCG_Ax(const
* array &x, array &ax)'. A virtual function "int LCG_Progress(const
* array &m, const double converge, const lcg_para ¶m, size_t t,
* std::ostream &ss)" could be reloaded for customed convergence tests.
*/
class lcg_solver
{
private:
lcg_para lcg_param_;
size_t lcg_inter_;
lcg_message_type lcg_msg_;
// make them class variables are more suitable for repetitively usages
array zk, gk, dk, Adk;
array rk, r0_T, pk, vk;
array Apx, uk, qk, qk_m, wk;
array m_new, gk_new;
array sk, yk;
/**
* @brief Display info of a given return code. This is a private function
* and can only be called by other class functions.
*
* @param err_code Input retrun code
* @param ss Output stream of runtime info.
*/
void lcg_error_str(lcg_return_code err_code, std::ostream &ss);
public:
lcg_solver(); ///< default constructor
virtual ~lcg_solver(); ///< default de-constructor
/**
* @brief Callback interface for calculating the product of 'A' multipled by an arbitrary vector 'x'.
*
* @param x Multipler
* @param ax Product of Ax
*/
virtual void LCG_Ax(const array &x, array &ax) = 0;
/**
* @brief Callback interface for calculating the product of 'M' multipled by an arbitrary vector 'x'.
* In which 'M' is the inverse of the pre-conditioning matrix. This function must be reloaded for the
* LCG_PCG algorithm.
*
* @param x Multipler
* @param mx Product of Ax
*/
virtual void LCG_Mx(const array &x, array &mx);
/**
* @brief Utility function for monitoring the solving process.
*
* @param m Current solution
* @param converge Current convergence
* @param param Employed parameters
* @param t Current iterative times
* @param ss Output stream of runtime info.
* @return Quit the solver if returned a non-zero value.
*/
virtual int LCG_Progress(const array &m, const double converge,
const lcg_para ¶m, size_t t, std::ostream &ss);
/**
* @brief Set the lcg message object
*
* @param msg Input message type.
*/
void set_lcg_message(lcg_message_type msg);
/**
* @brief Set the lcg report intervals
*
* @param inter Input reprot intervals.
*/
void set_lcg_report_interval(size_t inter);
/**
* @brief Set the lcg para object
*
* @param param Input lcg parameters.
*/
void set_lcg_para(const lcg_para ¶m);
/**
* @brief Return a lcg_para object with default values
*
* @return lcg_para
*/
lcg_para default_lcg_para();
#ifdef GCTL_OPTIMIZATION_TOML
/**
* @brief Set parameters of the conjugate gradient algorithms using a toml file.
* All parameter options must be listed under a top-level table 'lcg'. Available options
* under the 'lcg' table are as declared in the lcg_para structure.
*
* @param filename Input configuration file
*/
void set_lcg_para(std::string filename);
/**
* @brief Set parameters of the conjugate gradient algorithms using a toml::value object.
* All parameter options must be listed under a top-level table 'lcg'. Available options
* under the 'lcg' table are as declared in the lcg_para structure.
*
* @param toml_data Input toml data
*/
void set_lcg_para(const toml::value &toml_data);
#endif // GCTL_OPTIMIZATION_TOML
/**
* @brief The allback interface for all CG algorithms. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param solver_id Selected solver type.
* @param ss Output stream of runtime info.
*/
void LCG_Minimize(array &m, const array &b, lcg_solver_type solver_id = LCG_CG, std::ostream &ss = std::clog);
/**
* @brief The allback interface for all CG algorithms. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param solver_id Selected solver type.
* @param low Lower bound of acceptable solutions.
* @param hig Higher bound of acceptable solutions.
* @param ss Output stream of runtime info.
*/
void LCG_MinimizeConstrained(array &m, const array &b, const array &low, const array &hig, lcg_solver_type solver_id = LCG_PG, std::ostream &ss = std::clog);
/**
* @brief The standalone callback interface for the CG algorithm. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param ss Output stream of runtime info.
*/
void lcg(array &m, const array &B, std::ostream &ss = std::clog);
/**
* @brief The standalone callback interface for the PCG algorithm. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param ss Output stream of runtime info.
*/
void lpcg(array &m, const array &B, std::ostream &ss = std::clog);
/**
* @brief The standalone callback interface for the CGS algorithm. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param ss Output stream of runtime info.
*/
void lcgs(array &m, const array &B, std::ostream &ss = std::clog);
/**
* @brief The standalone callback interface for the BICGSTAB algorithm. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param ss Output stream of runtime info.
*/
void lbicgstab(array &m, const array &B, std::ostream &ss = std::clog);
/**
* @brief The standalone callback interface for the BICGSTAB2 algorithm. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param ss Output stream of runtime info.
*/
void lbicgstab2(array &m, const array &B, std::ostream &ss = std::clog);
/**
* @brief The standalone callback interface for the PG algorithm. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param low Lower bound of acceptable solutions.
* @param hig Higher bound of acceptable solutions.
* @param ss Output stream of runtime info.
*/
void lpg(array &m, const array &B, const array &low, const array &hig, std::ostream &ss = std::clog);
/**
* @brief The standalone callback interface for the SPG algorithm. Set message type to LCG_THROW to suppresses all info outputs.
*
* @param m Initial/Input solution.
* @param B Right hand term of the system system.
* @param low Lower bound of acceptable solutions.
* @param hig Higher bound of acceptable solutions.
* @param ss Output stream of runtime info.
*/
void lspg(array &m, const array &B, const array &low, const array &hig, std::ostream &ss = std::clog);
};
}
#endif // _GCTL_LCG_H