gctl_optimization/lib/optimization/lcg.h

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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 <http://www.gnu.org/licenses/>.
 * 
 * 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<double> &x, array<double> &ax)'. A virtual function "int LCG_Progress(const 
	 * array<double> &m, const double converge, const lcg_para &param, 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<double> zk, gk, dk, Adk;
		array<double> rk, r0_T, pk, vk;
    	array<double> Apx, uk, qk, qk_m, wk;
		array<double> m_new, gk_new;
    	array<double> 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<double> &x, array<double> &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<double> &x, array<double> &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<double> &m, const double converge, 
			const lcg_para &param, 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 &param);

		/**
		 * @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<double> &m, const array<double> &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<double> &m, const array<double> &b, const array<double> &low, const array<double> &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<double> &m, const array<double> &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<double> &m, const array<double> &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<double> &m, const array<double> &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<double> &m, const array<double> &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<double> &m, const array<double> &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<double> &m, const array<double> &B, const array<double> &low, const array<double> &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<double> &m, const array<double> &B, const array<double> &low, const array<double> &hig, std::ostream &ss = std::clog);
	};
}

#endif // _GCTL_LCG_H