generic_model.hpp

#ifndef CPPAD_CG_GENERIC_MODEL_INCLUDED
#define CPPAD_CG_GENERIC_MODEL_INCLUDED
/* --------------------------------------------------------------------------
 *  CppADCodeGen: C++ Algorithmic Differentiation with Source Code Generation:
 *    Copyright (C) 2012 Ciengis
 *    Copyright (C) 2020 Joao Leal
 *
 *  CppADCodeGen is distributed under multiple licenses:
 *
 *   - Eclipse Public License Version 1.0 (EPL1), and
 *   - GNU General Public License Version 3 (GPL3).
 *
 *  EPL1 terms and conditions can be found in the file "epl-v10.txt", while
 *  terms and conditions for the GPL3 can be found in the file "gpl3.txt".
 * ----------------------------------------------------------------------------
 * Author: Joao Leal
 */

namespace CppAD {
namespace cg {

template<class Base>
class GenericModel {
protected:
    CGAtomicGenericModel<Base>* _atomic;
    // whether or not to evaluate forward mode of atomics during a reverse sweep
    bool _evalAtomicForwardOne4CppAD;
public:

    GenericModel() :
        _atomic(nullptr),
        _evalAtomicForwardOne4CppAD(true) {
    }

    inline GenericModel(GenericModel&& other) noexcept :
            _atomic(other._atomic),
            _evalAtomicForwardOne4CppAD(other._evalAtomicForwardOne4CppAD) {
        other._atomic = nullptr;
    }

    inline virtual ~GenericModel() {
        delete _atomic;
    }

    virtual const std::string& getName() const = 0;


    virtual bool isJacobianSparsityAvailable() = 0;

    // Jacobian sparsity
    virtual std::vector<std::set<size_t> > JacobianSparsitySet() = 0;
    virtual std::vector<bool> JacobianSparsityBool() = 0;
    virtual void JacobianSparsity(std::vector<size_t>& equations,
                                  std::vector<size_t>& variables) = 0;

    virtual bool isHessianSparsityAvailable() = 0;

    virtual std::vector<std::set<size_t> > HessianSparsitySet() = 0;
    virtual std::vector<bool> HessianSparsityBool() = 0;
    virtual void HessianSparsity(std::vector<size_t>& rows,
                                 std::vector<size_t>& cols) = 0;

    virtual bool isEquationHessianSparsityAvailable() = 0;

    virtual std::vector<std::set<size_t> > HessianSparsitySet(size_t i) = 0;

    virtual std::vector<bool> HessianSparsityBool(size_t i) = 0;

    virtual void HessianSparsity(size_t i,
                                 std::vector<size_t>& rows,
                                 std::vector<size_t>& cols) = 0;

    virtual size_t Domain() const = 0;

    virtual size_t Range() const = 0;

    virtual const std::vector<std::string>& getAtomicFunctionNames() = 0;

    virtual bool addAtomicFunction(atomic_base<Base>& atomic) = 0;

    virtual bool addExternalModel(GenericModel<Base>& atomic) = 0;

    inline void setAtomicEvalForwardOne4CppAD(bool evalForwardOne4CppAD) {
        _evalAtomicForwardOne4CppAD = evalForwardOne4CppAD;
    }

    inline bool isAtomicEvalForwardOne4CppAD() const {
        return _evalAtomicForwardOne4CppAD;
    }

    /***********************************************************************
     *                        Forward zero
     **********************************************************************/

    virtual bool isForwardZeroAvailable() = 0;

    template<typename VectorBase>
    inline VectorBase ForwardZero(const VectorBase& x) {
        VectorBase dep(Range());
        this->ForwardZero(ArrayView<const Base>(&x[0], x.size()),
                          ArrayView<Base>(&dep[0], dep.size()));
        return dep;
    }

    virtual void ForwardZero(const CppAD::vector<bool>& vx,
                             CppAD::vector<bool>& vy,
                             ArrayView<const Base> tx,
                             ArrayView<Base> ty) = 0;

    template<typename VectorBase>
    inline void ForwardZero(const VectorBase& x,
                            VectorBase& dep) {
        dep.resize(Range());
        this->ForwardZero(ArrayView<const Base>(&x[0], x.size()),
                          ArrayView<Base>(&dep[0], dep.size()));
    }

    virtual void ForwardZero(ArrayView<const Base> x,
                             ArrayView<Base> dep) = 0;

    virtual void ForwardZero(const std::vector<const Base*> &x,
                             ArrayView<Base> dep) = 0;

    /***********************************************************************
     *                        Dense Jacobian
     **********************************************************************/

    virtual bool isJacobianAvailable() = 0;

    template<typename VectorBase>
    inline VectorBase Jacobian(const VectorBase& x) {
        VectorBase jac(Range() * Domain());
        Jacobian(ArrayView<const Base>(&x[0], x.size()),
                 ArrayView<Base>(&jac[0], jac.size()));
        return jac;
    }

    template<typename VectorBase>
    inline void Jacobian(const VectorBase& x,
                         VectorBase& jac) {
        jac.resize(Range() * Domain());
        Jacobian(ArrayView<const Base>(&x[0], x.size()),
                 ArrayView<Base>(&jac[0], jac.size()));
    }

    virtual void Jacobian(ArrayView<const Base> x,
                          ArrayView<Base> jac) = 0;

    /***********************************************************************
     *                        Dense Hessian
     **********************************************************************/

    virtual bool isHessianAvailable() = 0;


    template<typename VectorBase>
    inline VectorBase Hessian(const VectorBase& x,
                              const VectorBase& w) {
        VectorBase hess(Domain() * Domain());
        this->Hessian(ArrayView<const Base>(&x[0], x.size()),
                      ArrayView<const Base>(&w[0], w.size()),
                      ArrayView<Base>(&hess[0], hess.size()));
        return hess;
    }

    template<typename VectorBase>
    inline void Hessian(const VectorBase& x,
                        const VectorBase& w,
                        VectorBase& hess) {
        hess.resize(Domain() * Domain());
        this->Hessian(ArrayView<const Base>(&x[0], x.size()),
                      ArrayView<const Base>(&w[0], w.size()),
                      ArrayView<Base>(&hess[0], hess.size()));
    }

    template<typename VectorBase>
    inline VectorBase Hessian(const VectorBase& x,
                              size_t i) {
        CPPADCG_ASSERT_KNOWN(i < Range(), "Invalid equation index")

        VectorBase w(Range());
        w[i] = 1.0;
        VectorBase hess(Domain() * Domain());
        this->Hessian(ArrayView<const Base>(&x[0], x.size()),
                      ArrayView<const Base>(&w[0], w.size()),
                      ArrayView<Base>(&hess[0], hess.size()));
        return hess;
    }

    virtual void Hessian(ArrayView<const Base> x,
                         ArrayView<const Base> w,
                         ArrayView<Base> hess) = 0;

    /***********************************************************************
     *                        Forward one
     **********************************************************************/

    virtual bool isForwardOneAvailable() = 0;

    template<typename VectorBase>
    inline VectorBase ForwardOne(const VectorBase& tx) {
        size_t m = Range();
        const size_t k = 1;
        VectorBase ty((k + 1) * m);

        this->ForwardOne(ArrayView<const Base>(&tx[0], tx.size()),
                         ArrayView<Base>(&ty[0], ty.size()));

        VectorBase dy(m);
        for (size_t i = 0; i < m; i++) {
            dy[i] = ty[i * (k + 1) + k];
        }

        return dy;
    }

    virtual void ForwardOne(ArrayView<const Base> tx,
                            ArrayView<Base> ty) = 0;

    virtual bool isSparseForwardOneAvailable() = 0;

    virtual void ForwardOne(ArrayView<const Base> x,
                            size_t tx1Nnz, const size_t idx[], const Base tx1[],
                            ArrayView<Base> ty1) = 0;

    /***********************************************************************
     *                        Reverse one
     **********************************************************************/

    virtual bool isReverseOneAvailable() = 0;

    template<typename VectorBase>
    inline VectorBase ReverseOne(const VectorBase& tx,
                                 const VectorBase& ty,
                                 const VectorBase& py) {
        const size_t k = 0;
        VectorBase px((k + 1) * Domain());
        this->ReverseOne(tx, ty, px, py);
        return px;
    }

    template<typename VectorBase>
    inline void ReverseOne(const VectorBase& tx,
                           const VectorBase& ty,
                           VectorBase& px,
                           const VectorBase& py) {
        this->ReverseOne(ArrayView<const Base>(&tx[0], tx.size()),
                         ArrayView<const Base>(&ty[0], ty.size()),
                         ArrayView<Base>(&px[0], px.size()),
                         ArrayView<const Base>(&py[0], py.size()));
    }

    virtual bool isSparseReverseOneAvailable() = 0;

    virtual void ReverseOne(ArrayView<const Base> tx,
                            ArrayView<const Base> ty,
                            ArrayView<Base> px,
                            ArrayView<const Base> py) = 0;

    virtual void ReverseOne(ArrayView<const Base> x,
                            ArrayView<Base> px,
                            size_t pyNnz, const size_t idx[], const Base py[]) = 0;

    /***********************************************************************
     *                        Reverse two
     **********************************************************************/

    virtual bool isReverseTwoAvailable() = 0;

    template<typename VectorBase>
    inline VectorBase ReverseTwo(const VectorBase& tx,
                                 const VectorBase& ty,
                                 const VectorBase& py) {
        const size_t k = 1;
        VectorBase px((k + 1) * Domain());
        this->ReverseTwo(tx, ty, px, py);
        return px;
    }

    template<typename VectorBase>
    inline void ReverseTwo(const VectorBase& tx,
                           const VectorBase& ty,
                           VectorBase& px,
                           const VectorBase& py) {
        this->ReverseTwo(ArrayView<const Base>(&tx[0], tx.size()),
                         ArrayView<const Base>(&ty[0], ty.size()),
                         ArrayView<Base>(&px[0], px.size()),
                         ArrayView<const Base>(&py[0], py.size()));
    }

    virtual bool isSparseReverseTwoAvailable() = 0;

    virtual void ReverseTwo(ArrayView<const Base> tx,
                            ArrayView<const Base> ty,
                            ArrayView<Base> px,
                            ArrayView<const Base> py) = 0;
    virtual void ReverseTwo(ArrayView<const Base> x,
                            size_t tx1Nnz, const size_t idx[], const Base tx1[],
                            ArrayView<Base> px2,
                            ArrayView<const Base> py2) = 0;

    /***********************************************************************
     *                        Sparse Jacobians
     **********************************************************************/

    virtual bool isSparseJacobianAvailable() = 0;

    template<typename VectorBase>
    inline VectorBase SparseJacobian(const VectorBase& x) {
        VectorBase jac(Range() * Domain());
        SparseJacobian(ArrayView<const Base>(&x[0], x.size()),
                       ArrayView<Base>(&jac[0], jac.size()));
        return jac;
    }

    template<typename VectorBase>
    inline void SparseJacobian(const VectorBase& x,
                               VectorBase& jac) {
        jac.resize(Range() * Domain());
        SparseJacobian(ArrayView<const Base>(&x[0], x.size()),
                       ArrayView<Base>(&jac[0], jac.size()));
    }

    virtual void SparseJacobian(ArrayView<const Base> x,
                                ArrayView<Base> jac) = 0;

    virtual void SparseJacobian(const std::vector<Base> &x,
                                std::vector<Base>& jac,
                                std::vector<size_t>& row,
                                std::vector<size_t>& col) = 0;

    virtual void SparseJacobian(ArrayView<const Base> x,
                                ArrayView<Base> jac,
                                size_t const** row,
                                size_t const** col) = 0;

    virtual void SparseJacobian(const std::vector<const Base*>& x,
                                ArrayView<Base> jac,
                                size_t const** row,
                                size_t const** col) = 0;

    /***********************************************************************
     *                        Sparse Hessians
     **********************************************************************/

    virtual bool isSparseHessianAvailable() = 0;

    template<typename VectorBase>
    inline VectorBase SparseHessian(const VectorBase& x,
                                    const VectorBase& w) {
        VectorBase hess(Domain() * Domain());
        SparseHessian(ArrayView<const Base>(&x[0], x.size()),
                      ArrayView<const Base>(&w[0], w.size()),
                      ArrayView<Base>(&hess[0], hess.size()));
        return hess;
    }

    template<typename VectorBase>
    inline void SparseHessian(const VectorBase& x,
                              const VectorBase& w,
                              VectorBase& hess) {
        hess.resize(Domain() * Domain());
        SparseHessian(ArrayView<const Base>(&x[0], x.size()),
                      ArrayView<const Base>(&w[0], w.size()),
                      ArrayView<Base>(&hess[0], hess.size()));
    }

    virtual void SparseHessian(ArrayView<const Base> x,
                               ArrayView<const Base> w,
                               ArrayView<Base> hess) = 0;

    virtual void SparseHessian(const std::vector<Base>& x,
                               const std::vector<Base>& w,
                               std::vector<Base>& hess,
                               std::vector<size_t>& row,
                               std::vector<size_t>& col) = 0;

    virtual void SparseHessian(ArrayView<const Base> x,
                               ArrayView<const Base> w,
                               ArrayView<Base> hess,
                               size_t const** row,
                               size_t const** col) = 0;

    virtual void SparseHessian(const std::vector<const Base*>& x,
                               ArrayView<const Base> w,
                               ArrayView<Base> hess,
                               size_t const** row,
                               size_t const** col) = 0;

    virtual CGAtomicGenericModel<Base>& asAtomic() {
        if (_atomic == nullptr) {
            _atomic = new CGAtomicGenericModel<Base>(*this);
        }
        return *_atomic;
    }
};

} // END cg namespace
} // END CppAD namespace

#endif