Lemma/Modules/Optimization/include/cg.h

/* This file is part of Lemma, a geophysical modelling and inversion API */

/* This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

/**
  @file
  @author   Trevor Irons
  @version  $Id: cg.h 87 2013-09-05 22:44:05Z tirons $
 **/

#ifndef  CG_INC
#define  CG_INC


#include <iostream>
#include <fstream>
#include "lemma.h"

namespace Lemma {

    /** Port of netlib.org http://www.netlib.org/templates/cpp//cg.h
     *  Solves the symmetric postive definite system Ax=b.
     *  No preconditioner is used.
     *  Iterative template routine -- CG
     *
     *  CG solves the symmetric positive definite linear
     *  system Ax=b using the Conjugate Gradient method.
     *
     *  CG follows the algorithm described on p. 15 in the
     *  SIAM Templates book.
     *
     *  The return value indicates convergence within max_iter (input)
     *  iterations (0), or no convergence within max_iter iterations (1).
     *
     *  Upon successful return, output arguments have the following values:
     *
     *         x  --  approximate solution to Ax = b
     *  max_iter  --  the number of iterations performed before the
     *                tolerance was reached
     *       tol  --  the residual after the final iteration
     *  @param[in] A is a Real matrix A to be solved.
     *  @param[in] x0 is the starting model.
     *  @param[in] b is the right hand side.
     */
    template < typename Scalar >
    VectorXr CG(const MatrixXr &A, const VectorXr &x0, const VectorXr &b,
        int &max_iter, Scalar &tol) {

        VectorXr p, q;
        Scalar alpha, beta, rho, rho_1(0);
        VectorXr x = x0;
        Scalar normb = b.norm( );
        VectorXr r = b - A*x;

        if (normb == 0.0) {
            normb = 1;
        }

        Scalar resid = r.norm() / normb;
        if (resid <= tol) {
            tol = resid;
            max_iter = 0;
            return x;
        }

        for (int i = 1; i <= max_iter; i++) {

            rho = r.transpose()*r;

            if (i == 1) {
                p = r;
            }
            else {
                beta = rho / rho_1;
                p = r + beta * p;
            }

            q = A*p;
            alpha = rho / p.dot(q);

            x += alpha * p;
            r -= alpha * q;

            if ((resid = r.norm() / normb) <= tol) {
                tol = resid;
                max_iter = i;
                return x;
            }
            rho_1 = rho;
        }
        tol = resid;
        std::cerr << "CG FAILED TO REACH CONVERGENCE\n";
        return x;
    }

    // Specialised routine that appliex mod(Ax) so that b and x are in real, but
    // A is complex.
    template < typename Scalar >
    VectorXr CG_ModulusAx(const MatrixXcr &A, const VectorXr &x0, const VectorXr &b,
        int &max_iter, Scalar &tol) {

        VectorXr p, q;
        Scalar beta, rho, rho_1(0);
        Scalar alpha;
        VectorXr x = x0;
        Scalar normb = b.norm( );
        VectorXr r = b.array() - (A*x).array().abs();

        if (normb == 0.0) {
            normb = 1;
        }

        Scalar resid = r.norm() / normb;
        if (resid <= tol) {
            tol = resid;
            max_iter = 0;
            return x;
        }

        for (int i = 1; i <= max_iter; i++) {

            rho = r.dot(r); //conjugate().transpose()*r;

            if (i == 1) {
                p = r;
            }
            else {
                beta = rho / rho_1;
                p = r + beta * p;
            }

            q = (A*p).array().abs();
            alpha = rho / p.dot(q);

            x += alpha * p;
            r -= alpha * q;

            if ((resid = r.norm() / normb) <= tol) {
                tol = resid;
                max_iter = i;
                return x;
            }
            rho_1 = rho;
        }

        tol = resid;
        std::cerr << "CG FAILED TO REACH CONVERGENCE\n";
        return x;
    }

// Preconditioned version of above
//*****************************************************************
// Iterative template routine -- CG
//
// CG solves the symmetric positive definite linear
// system Ax=b using the Conjugate Gradient method.
//
// CG follows the algorithm described on p. 15 in the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
//        x  --  approximate solution to Ax = b
// max_iter  --  the number of iterations performed before the
//               tolerance was reached
//      tol  --  the residual after the final iteration
//
//*****************************************************************

//template < class Matrix, class Vector, class Preconditioner, class Real >
//int
//CG(const Matrix &A, Vector &x, const Vector &b,
//   const Preconditioner &M, int &max_iter, Real &tol)
//{
    #include <limits>
    template <typename Preconditioner>
    VectorXr CG(const MatrixXr &A, const VectorXr &x0, const VectorXr &b,
        const Preconditioner &M, int &max_iter, Real &tol) {
        //const Eigen::SparseMatrix<Real> &M, int &max_iter, Real &tol) {

        VectorXr p, z, q;
        VectorXr x = x0;
        Real alpha(0), beta(0), rho(0), rho_1(0);
        Real normb;
        VectorXr r;
        #ifdef LEMMAUSEOMP
        #pragma omp parallel sections
        {
        #endif
            #ifdef LEMMAUSEOMP
            #pragma omp section
            #endif
            {
                normb = b.norm();
            }
            #ifdef LEMMAUSEOMP
            #pragma omp section
            #endif
            {
                r = b - A*x;
            }
        #ifdef LEMMAUSEOMP
        }
        #endif
        if (normb <= std::numeric_limits<double>::epsilon() ) {
            normb = 1;
        }
        Real resid = r.norm() / normb;
        if (resid <= tol) {
            tol = resid;
            max_iter = 0;
            return x;
        }
        // todo do 0th loop manually, gets rid of if statement
        for (int i = 1; i <= max_iter; i++) {
            z = M.solve(r);
            //z = M*r;
            rho = r.transpose()*z;
            if (i == 1) {
                p = z;
            } else {
                beta = rho / rho_1;
                p = z + beta * p;
            }
            q = A*p;
            alpha = rho / p.dot(q);
            #ifdef LEMMAUSEOMP
            #pragma omp parallel sections
            {
            #endif
                #ifdef LEMMAUSEOMP
                #pragma omp section
                #endif
                {
                    x += alpha * p;
                }
                #ifdef LEMMAUSEOMP
                #pragma omp section
                #endif
                {
                    r -= alpha * q;
                }
            #ifdef LEMMAUSEOMP
            }
            #endif
            if ((resid = r.norm() / normb) <= tol) {
               tol = resid;
               max_iter = i;
               return x;
            }
            rho_1 = rho;
        }
        tol = resid;
        std::cerr << "Preconditioned CG failed to converge\n";
        return x;
    }

    template < typename Scalar >
    VectorXr CGJ(const MatrixXr &A, const VectorXr &x0, const VectorXr &b,
        const Eigen::SparseMatrix<Scalar> &M, int &max_iter, Scalar &tol) {

        VectorXr p, z, q;
        VectorXr x = x0;
        Scalar alpha(0), beta(0), rho(0), rho_1(0);
        Scalar normb;
        VectorXr r;
        #ifdef LEMMAUSEOMP
        #pragma omp parallel sections
        {
        #endif
            #ifdef LEMMAUSEOMP
            #pragma omp section
            #endif
            {
                normb = b.norm();
            }
            #ifdef LEMMAUSEOMP
            #pragma omp section
            #endif
            {
                r = b - A*x;
            }
        #ifdef LEMMAUSEOMP
        }
        #endif
        if (normb <= std::numeric_limits<double>::epsilon() ) {
            normb = 1;
        }
        Scalar resid = r.norm() / normb;
        if (resid <= tol) {
            tol = resid;
            max_iter = 0;
            return x;
        }
        // todo do 0th loop manually, gets rid of if statement
        for (int i = 1; i <= max_iter; i++) {
            //z = M.solve(r);
            z = M*r;
            rho = r.transpose()*z;
            if (i == 1) {
                p = z;
            } else {
                beta = rho / rho_1;
                p = z + beta * p;
            }
            q = A*p;
            alpha = rho / p.dot(q);
            #ifdef LEMMAUSEOMP
            #pragma omp parallel sections
            {
            #endif
                #ifdef LEMMAUSEOMP
                #pragma omp section
                #endif
                {
                    x += alpha * p;
                }
                #ifdef LEMMAUSEOMP
                #pragma omp section
                #endif
                {
                    r -= alpha * q;
                }
            #ifdef LEMMAUSEOMP
            }
            #endif
            if ((resid = r.norm() / normb) <= tol) {
               tol = resid;
               max_iter = i;
               return x;
            }
            rho_1 = rho;
        }
        tol = resid;
        std::cerr << "Preconditioned CG failed to converge\n";
        return x;
    }

    ///////////////////////////////////////////////////////
    // Log Barrier
    /** Solves one iteration, using implicit multiplication of At*b, and ATA.
     *  @param[in] G is the sensitivity matrix to be inverted.
     *  @param[in] WdTWd is the data weighting matrix, usually sparse.
     *  @param[in] WmTWm is the model weighting matrix, usually sparse.
     *  @param[in] X1 is a vector of the inverse of the current model X. Used in
     *  log barrier. Should be computed as X1 = VectorXr( 1. / (x.array() - minVal ))
     *  @param[in] Y1 is a vector of the inverse of the current model X. Used in
     *  log barrier. Should be computed as X1 = VectorXr( 1. / (maxVal - x.array() ))
     *  @param[in] X2 is the analaogous inverse of X^2 diagonal matrix, stored
     *  as a Vector.
     *  @param[in] Y2 is the analaogous inverse of (maxVal-X)^2 diagonal matrix, stored
     *  as a Vector.
     *  @param[in] D2 is the datamisfit vector, formally (d_predicted - d_observed).
     *  @param[in,out] Mk is the input / output solution vector. The input value is
     *  an itital guess, and output is the solution.
     *  @param[in] Mr is the reference model.
     *  @param[in] BETA is the regularisation (Tikhonov parameter) to apply
     *  @param[in] P is a preconditioner of G. The product Pb ~ x
     *  @param[in,out] max_iter is the number of iterations.
     *  @param[in,out] tol is the desired tolerance on input, and achieved on
     *  output.
     */
    template < typename  Scalar >
    int implicit_log_CG(const MatrixXr& G, const MatrixXr& WdTWd, const VectorXr& WmTWm,
				const VectorXr& X1, const VectorXr& X2,
				const VectorXr& Y1, const VectorXr& Y2,
				const VectorXr& D2, VectorXr& Mk, const VectorXr& Mr,
				const Scalar& BETA, const MatrixXr& P,
                int& max_iter, Scalar& tol) {

        // TODO add WdTWD to this!

        Scalar resid;
        VectorXr p, z, q;
        Scalar alpha(0), beta(0), rho(0), rho_1(0);

        // Calculate 'B'
        VectorXr delM = Mk - Mr;
        VectorXr B = (-G.transpose()*D2).array() - BETA*((WmTWm.asDiagonal()*delM).array())
                     + X1.array() + Y1.array();

        Scalar normb = B.norm();

        // Implicit calc of A*x
        VectorXr Y = BETA*(WmTWm.asDiagonal()*Mk).array() +
                (X2.asDiagonal()*Mk).array() + (Y2.asDiagonal()*Mk).array();
        VectorXr Z = G*Mk;
        VectorXr U = G.transpose()*Z;

        VectorXr r = B - (Y + U);

        if (normb == 0.0) normb = 1;

        if ((resid = r.norm() / normb) <= tol) {
            tol = resid;
            max_iter = 0;
            return 0;
        }

        for (int i = 1; i <= max_iter; i++) {
            //z = M.solve(r); // we can solve directly z = P*r
            z = P*r;
            rho = r.dot(z);
            if (i == 1)  p = z;
            else {
                beta = rho / rho_1;
                p = beta * p;
 	            p = p+z;
            }
            Y = BETA*(WmTWm.array()*p.array()) + X2.array()*p.array() + Y2.array()*p.array();
            Z = G*p;
            U = G.transpose()*Z;
            q = Y+U;

            alpha = rho / p.dot(q);

            Mk = Mk + alpha * p;
            r  = r - (alpha * q);

            if ((resid = r.norm() / normb) <= tol) {
                tol = resid;
                max_iter = i;
                return 0;
            }
            rho_1 = rho;
        }
        tol = resid;
        return 1;
    }
} // end of namespace Lemma

#endif   // ----- #ifndef CG_INC  -----