ADMM.hpp

/*
 * This file is part of bogus, a C++ sparse block matrix library.
 *
 * Copyright 2016 Gilles Daviet <gdaviet@gmail.com>
 *
 * 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/.
*/


#ifndef BOGUS_ADMM_HPP
#define BOGUS_ADMM_HPP

#include "ConstrainedSolverBase.hpp"

#include <vector>

namespace bogus
{

namespace admm {
    enum Variant {
        Standard,
        Accelerated
    } ;
}


template < typename BlockMatrixType >
class ADMM : public ConstrainedSolverBase< ADMM<BlockMatrixType >, BlockMatrixType >
{
public:
    typedef ConstrainedSolverBase< ADMM, BlockMatrixType > Base ;

    typedef typename Base::GlobalProblemTraits GlobalProblemTraits ;
    typedef typename GlobalProblemTraits::Scalar Scalar ;

    ADMM( ) : Base() { init() ; }
    explicit ADMM( const BlockObjectBase< BlockMatrixType > & matrix ) : Base()
    { init() ; Base::setMatrix( matrix ) ; }

    template < admm::Variant variant, typename NSLaw, typename ProxOp, typename RhsT, typename ResT >
    Scalar solve(
            const NSLaw &law, const ProxOp& op,
            const RhsT &w, ResT &v, ResT &r ) const ;

    template < typename NSLaw, typename ProxOp, typename RhsT, typename ResT >
    Scalar solve(
            const NSLaw &law, const ProxOp& op,
            const RhsT &w, ResT &v, ResT &r ) const ;

    ADMM& setMatrix( const BlockObjectBase< BlockMatrixType > & matrix )
    {
        m_matrix = &matrix ;
        Base::updateScalings() ;
        return *this ;
    }

    void setStepSize( const Scalar size )
    { m_stepSize = size ; }

    void setDefaultVariant( admm::Variant variant )
    { m_defaultVariant = variant ; }

    Scalar stepSize() const { return m_stepSize ; }

protected:

    typedef typename Base::Index Index ;

    void init()
    {
        m_tol = 1.e-6 ;
        m_maxIters = 300 ;

        m_stepSize = 1.e-2 ;

        m_defaultVariant = admm::Accelerated ;
    }

    using Base::m_matrix ;
    using Base::m_maxIters ;
    using Base::m_tol ;

    Scalar m_stepSize ;

    admm::Variant m_defaultVariant ;

} ;


template< typename ObjectType >
struct QuadraticProxOp
{
    typedef BlockMatrixTraits< ObjectType > BlockTraits;
    typedef typename BlockTraits::Scalar Scalar ;
    typedef ProblemTraits< Scalar > GlobalProblemTraits ;

    typedef BlockObjectBase< ObjectType > LinearOp ;
    typedef typename GlobalProblemTraits::DynVector AffineVec ;


    QuadraticProxOp( const LinearOp& linearExpr, const Scalar lambda,
                     const AffineVec& affinePart )
        : m_coefficient( lambda ), m_linearOp( linearExpr ),
          m_affinePart( affinePart )
    {
    }


    template < typename RhsT, typename ResT >
    void eval( RhsT & rhs, ResT & res ) const ;

    Scalar coefficient() const {
        return m_coefficient ;
    }

private:
    const Scalar     m_coefficient ;
    const LinearOp&  m_linearOp ;
    const AffineVec& m_affinePart ;
};



template < typename BlockMatrixType >
class DualAMA : public ConstrainedSolverBase< DualAMA<BlockMatrixType>, BlockMatrixType >
{
public:
    typedef ConstrainedSolverBase< DualAMA, BlockMatrixType > Base ;

    typedef typename Base::GlobalProblemTraits GlobalProblemTraits ;
    typedef typename GlobalProblemTraits::Scalar Scalar ;

    DualAMA( ) : Base() { init() ; }
    explicit DualAMA( const BlockObjectBase< BlockMatrixType > & matrix ) : Base()
    { init() ; Base::setMatrix( matrix ) ; }

    template < typename NSLaw, typename MatrixT, typename RhsT, typename ResT >
    Scalar solve(
            const NSLaw &law, const BlockObjectBase< MatrixT >& A,
            const RhsT &f, const RhsT &w,
            ResT &v, ResT &r ) const ;

    template < admm::Variant variant, typename NSLaw, typename MatrixT,
               typename RhsT, typename ResT >
    Scalar solve(
            const NSLaw &law, const BlockObjectBase< MatrixT >& A,
            const RhsT &f, const RhsT &w,
            ResT &v, ResT &r ) const ;

    template < admm::Variant variant, typename NSLaw,
               typename AType, typename BType, typename HType, typename PrecondT,
               typename RhsT, typename ORhsT, typename ResT, typename OResT >
    Scalar solveWithLinearConstraints(
            const NSLaw &law,
            const BlockObjectBase< AType >& A,
            const BlockObjectBase< BType >& B,
            const BlockObjectBase< HType >& H,
            const PrecondT& preconditioner,
            const RhsT &f, const ORhsT& k, const RhsT &w,
            ResT &v, OResT &p, ResT &r, Scalar stepRatio = 1 ) const ;

    DualAMA& setMatrix( const BlockObjectBase< BlockMatrixType > & matrix )
    {
        m_matrix = &matrix ;
        Base::updateScalings() ;
        return *this ;
    }

    void setFpStepSize( const Scalar size )
    { m_fpStepSize = size ; }

    void setProjStepSize( const Scalar size )
    { m_projStepSize = size ; }

    void setDefaultVariant( admm::Variant variant )
    { m_defaultVariant = variant ; }

    Scalar fpStepSize() const { return m_fpStepSize ; }
    Scalar projStepSize() const { return m_projStepSize ; }

    // LS

    void setLineSearchIterations( const unsigned lsIterations )
    { m_lsIters = lsIterations ; }

    void setLineSearchOptimisticFactor( const Scalar lsOptimisticFactor )
    { m_lsOptimisticFactor = lsOptimisticFactor ; }

    void setLineSearchPessimisticFactor( const Scalar lsPessimisticFactor )
    { m_lsPessimisticFactor = lsPessimisticFactor ; }

    unsigned lineSearchIterations() const { return m_lsIters ; }
    Scalar lineSearchOptimisticFactor() const { return m_lsOptimisticFactor ; }
    Scalar lineSearchPessimisticFactor() const { return m_lsPessimisticFactor ; }

protected:

    typedef typename Base::Index Index ;

    void init()
    {
        m_tol = 1.e-6 ;
        m_maxIters = 300 ;

        m_fpStepSize = 1.e-1 ;
        m_projStepSize = 1.e-1 ;

        m_lsIters = 6 ;
        m_lsOptimisticFactor = 1.25 ;
        m_lsPessimisticFactor = .5 ;

        m_defaultVariant = admm::Accelerated ;
    }

    using Base::m_matrix ;
    using Base::m_maxIters ;
    using Base::m_tol ;

    Scalar m_fpStepSize ;
    Scalar m_projStepSize ;

    admm::Variant m_defaultVariant ;

    unsigned m_lsIters ;
    Scalar m_lsOptimisticFactor ;
    Scalar m_lsPessimisticFactor ;

} ;

} //namespace bogus


#endif