EvaluableLinearCombination.hpp

Go to the documentation of this file.

// This file is part of GQCG-GQCP.
//
// Copyright (C) 2017-2020  the GQCG developers
//
// GQCG-GQCP is free software: 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 3 of the License, or
// (at your option) any later version.
//
// GQCG-GQCP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
 
#include "Mathematical/Functions/Function.hpp"
#include "Mathematical/Functions/VectorSpaceArithmetic.hpp"
#include "Utilities/aliases.hpp"
#include "Utilities/complex.hpp"
 
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>
 
#include <vector>
 
 
namespace GQCP {
 
 
template <typename _Coefficient, typename _FunctionType>
class EvaluableLinearCombination:
    public Function<sum_t<_Coefficient, typename _FunctionType::OutputType>, typename _FunctionType::InputType>,
    public VectorSpaceArithmetic<EvaluableLinearCombination<_Coefficient, _FunctionType>, _Coefficient> {
 
public:
    // The type of a coefficient.
    using Coefficient = _Coefficient;
 
    // The type of the function.
    using FunctionType = _FunctionType;
 
    // The return type of the `operator()`.
    using OutputType = sum_t<Coefficient, typename FunctionType::OutputType>;
 
    // The input type of the `operator()`.
    using InputType = typename FunctionType::InputType;
 
    // static_assert(std::is_base_of<Function<OutputType, InputType>, FunctionType>::value, "EvaluableLinearCombination: FunctionType must derive from `Function`.");
 
    // The type of 'this'.
    using Self = EvaluableLinearCombination<Coefficient, FunctionType>;
 
 
protected:
    // The coefficients of the linear combination.
    std::vector<Coefficient> m_coefficients;
 
    // The functions of the linear combination.
    std::vector<FunctionType> m_functions;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    EvaluableLinearCombination(const std::vector<Coefficient>& coefficients, const std::vector<FunctionType>& functions) :
        m_coefficients {coefficients},
        m_functions {functions} {
 
        if (coefficients.size() != functions.size()) {
            throw std::invalid_argument("EvaluableLinearCombination(const std::vector<Coefficient>&, const std::vector<FunctionType>&): The number of coefficients and functions do not match.");
        }
    }
 
 
    EvaluableLinearCombination() :
        EvaluableLinearCombination(std::vector<Coefficient> {}, std::vector<FunctionType> {}) {}
 
 
    EvaluableLinearCombination(const Coefficient coefficient, const FunctionType& function) :
        EvaluableLinearCombination(std::vector<Coefficient> {coefficient}, std::vector<FunctionType> {function}) {}
 
 
    EvaluableLinearCombination(const FunctionType& function) :
        EvaluableLinearCombination(Coefficient {1.0}, function) {}
 
 
    EvaluableLinearCombination(const int zero) :
        EvaluableLinearCombination() {
 
        if (zero != 0) {
            throw std::invalid_argument("EvaluableLinearCombination(const int): Can't convert a non-zero integer to a 'zero' instance.");
        }
    }
 
 
    /*
     *  MARK: Vector space arithmetic
     */
 
    Self& operator+=(const Self& rhs) override {
        this->append(rhs.m_coefficients, rhs.m_functions);
        return *this;
    }
 
 
    Self& operator*=(const Coefficient& a) override {
 
        if (std::abs(a) < 1.0e-16) {  // multiplication by zero returns a 'zero vector'
            this->m_coefficients = std::vector<Coefficient> {};
            this->m_functions = std::vector<FunctionType> {};
            return *this;
        }
 
        std::transform(this->m_coefficients.begin(), this->m_coefficients.end(),
                       this->m_coefficients.begin(), [a](const Coefficient& C) { return C * a; });
 
        return *this;
    }
 
 
    /*
     *  MARK: General information
     */
 
    std::string description() const {
 
        std::string description = "[";
        for (size_t i = 0; i < this->length(); i++) {
 
            // Provide the coefficient and the function's description.
            description += (boost::format("(%|.3f|) %s") % this->coefficient(i) % this->function(i).description()).str();
 
            // Add a '+' when we're not at the end.
            if (i != this->length() - 1) {
                description += " + ";
            } else {
                description += "]";
            }
        }
 
        return description;
    }
 
    size_t length() const { return this->m_coefficients.size(); }
 
 
    /*
     *  MARK: Access
     */
 
    const Coefficient& coefficient(const size_t i) const { return this->m_coefficients[i]; }
 
    const std::vector<Coefficient>& coefficients() const { return this->m_coefficients; }
 
    const FunctionType& function(const size_t i) const { return this->m_functions[i]; }
 
    const std::vector<FunctionType>& functions() const { return this->m_functions; }
 
 
    /*
     *  MARK: Conforming to `Function`
     */
 
    OutputType operator()(const InputType& in) const override {
 
        // Evaluate every function of the linear combination.
        OutputType out {};  // Default initialization of the `OutputType`.
        const auto n = this->length();
        for (size_t i = 0; i < n; i++) {
            out += this->m_coefficients[i] * this->m_functions[i].operator()(in);
        }
 
        return out;
    }
 
 
    /*
     *  MARK: Comparison operators
     */
 
    bool operator==(const Self& other) const {
        return (this->m_coefficients == other.m_coefficients) && (this->m_functions == other.m_functions);
    }
 
 
    /*
     *  MARK: Appending
     */
 
    void append(const Coefficient& coefficient, const FunctionType& function) {
 
        this->m_coefficients.push_back(coefficient);
        this->m_functions.push_back(function);
    }
 
 
    void appendWithThreshold(const Coefficient& coefficient, const FunctionType& function, const double threshold = 1.0e-16) {
 
        // Only enlarge the linear combination for sufficiently large coefficients.
        if (std::abs(coefficient) < threshold) {
            return;
        } else {
            this->append(coefficient, function);
        }
    }
 
 
    void append(const std::vector<Coefficient>& coefficients, const std::vector<FunctionType>& functions) {
 
        if (coefficients.size() != functions.size()) {
            throw std::invalid_argument("EvaluableLinearCombination::append(const std::vector<Coefficient>&, const std::vector<FunctionType>&): The number of coefficients and functions do not match.");
        }
 
        this->m_coefficients.insert(this->m_coefficients.end(), coefficients.begin(), coefficients.end());
        this->m_functions.insert(this->m_functions.end(), functions.begin(), functions.end());
    }
};
 
 
}  // namespace GQCP
 
 
#include <Eigen/Dense>
 
/*
 *  Make GQCP::EvaluableLinearCombination<FunctionType> an Eigen scalar type.
 */
 
namespace Eigen {
 
 
template <typename Coefficient, typename FunctionType>
struct NumTraits<GQCP::EvaluableLinearCombination<Coefficient, FunctionType>>:
    public NumTraits<double> {  // permits to get the epsilon, dummy_precision, lowest, highest functions
 
    using Real = GQCP::EvaluableLinearCombination<Coefficient, FunctionType>;
    using NonInteger = GQCP::EvaluableLinearCombination<Coefficient, FunctionType>;
    using Nested = GQCP::EvaluableLinearCombination<Coefficient, FunctionType>;
 
    enum {
        IsComplex = 0,
        IsInteger = 0,
        IsSigned = 1,
        RequireInitialization = 1,
        ReadCost = 1,
        AddCost = 5,
        MulCost = 1000  // just put something big
    };
};
 
 
// Enable the scalar product of a EvaluableLinearCombination<Coefficient, FunctionType> with its own Coefficient (both sides).
 
template <typename Coefficient, typename FunctionType>
struct ScalarBinaryOpTraits<GQCP::EvaluableLinearCombination<Coefficient, FunctionType>, Coefficient,
                            Eigen::internal::scalar_product_op<GQCP::EvaluableLinearCombination<Coefficient, FunctionType>, Coefficient>> {
 
    using ReturnType = GQCP::EvaluableLinearCombination<Coefficient, FunctionType>;
};
 
template <typename Coefficient, typename FunctionType>
struct ScalarBinaryOpTraits<Coefficient, GQCP::EvaluableLinearCombination<Coefficient, FunctionType>,
                            Eigen::internal::scalar_product_op<Coefficient, GQCP::EvaluableLinearCombination<Coefficient, FunctionType>>> {
 
    using ReturnType = GQCP::EvaluableLinearCombination<Coefficient, FunctionType>;
};
 
 
}  // namespace Eigen