linear_forms_block.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//linear_forms_block.hpp is free software: you can redistribute it and/or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//linear_forms_block.hpp 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 General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with linear_forms_block.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
// Copyright(C) Bertini2 Development Team
//
// See <http://www.gnu.org/licenses/> for a copy of the license,
// as well as COPYING.  Bertini2 is provided with permitted
// additional terms in the b2/licenses/ directory.

 
#pragma once
 
#include <stdexcept>
#include <string>
#include <tuple>
 
#include "bertini2/num_traits.hpp"
#include "bertini2/eigen_extensions.hpp"
#include "bertini2/system/blocks/describe.hpp"
 
namespace bertini {
namespace blocks {

class LinearFormsBlock
{
public:
    LinearFormsBlock() : num_vars_(0), precision_(DefaultPrecision()) {}

    LinearFormsBlock(size_t num_vars, Mat<complex_mp> coefficients)
        : num_vars_(num_vars), coefficients_highest_precision_(std::move(coefficients)),
          precision_(DefaultPrecision())
    {
        assert(static_cast<size_t>(coefficients_highest_precision_.cols()) == num_vars_ + 1 &&
               "a linear-forms coefficient matrix must have num_vars+1 columns");
        BuildWorking();
    }

    LinearFormsBlock(size_t num_vars, Mat<complex_mp> coefficients, bool homogenized)
        : num_vars_(num_vars), coefficients_highest_precision_(std::move(coefficients)),
          precision_(DefaultPrecision())
    {
        homogeneous_ = homogenized;
        if (static_cast<size_t>(coefficients_highest_precision_.cols()) != num_vars_ + (homogenized ? 0 : 1))
            throw std::invalid_argument("LinearFormsBlock: the coefficient matrix has "
                + std::to_string(coefficients_highest_precision_.cols()) + " columns, but a block on "
                + std::to_string(num_vars_) + " variables that " + (homogenized ? "is" : "is not")
                + " homogenized needs " + std::to_string(num_vars_ + (homogenized ? 0 : 1)));
        BuildWorking();
    }

    size_t NumFunctions() const { return static_cast<size_t>(coefficients_highest_precision_.rows()); }

    std::vector<int> Degrees() const { return std::vector<int>(NumFunctions(), 1); }
    std::vector<int> Degrees(VariableGroup const&) const { return Degrees(); }

    bool IsPolynomial(VariableGroup const&) const { return true; }
    bool IsHomogeneous(VariableGroup const&) const { return homogeneous_; }

    void Homogenize(VariableGroup const& /*group*/, std::shared_ptr<node::Variable> const& /*hom_var*/)
    {
        if (homogeneous_)
            throw std::runtime_error("LinearFormsBlock::Homogenize: block is already homogenized "
                "(multiple affine variable groups are not yet supported for linear-forms blocks)");
        const auto& M = coefficients_highest_precision_;
        Mat<complex_mp> Mh(M.rows(), M.cols());                                   // same #cols: n+1
        Mh.col(0) = M.col(static_cast<Eigen::Index>(num_vars_));                    // constant -> h column (front)
        Mh.rightCols(static_cast<Eigen::Index>(num_vars_)) =
            M.leftCols(static_cast<Eigen::Index>(num_vars_));                       // original variable columns
        coefficients_highest_precision_ = Mh;
        num_vars_ += 1;          // h is now a real variable; all n+1 columns are variable columns
        homogeneous_ = true;
        BuildWorking();
    }

    size_t NumVariables() const { return num_vars_; }

    Mat<complex_mp> const& Coefficients() const { return coefficients_highest_precision_; }
    bool IsHomogenized() const { return homogeneous_; }

    void Describe(std::ostream& out, size_t& row, VariableGroup const& vars, bool verbose) const
    {
        auto const& M = coefficients_highest_precision_;
        const Eigen::Index n = M.rows();
        if (n == 0)
            return;
        const Eigen::Index cap   = static_cast<Eigen::Index>(describe_detail::kTerseRowCap);
        const Eigen::Index shown = (verbose || n <= cap) ? n : cap;
 
        // placeholder line per form: f_k = c.[x, y, 1]
        for (Eigen::Index r = 0; r < shown; ++r)
        {
            out << "  f_" << row++ << " = c.";
            describe_detail::PrintAugmentedVars(out, vars, num_vars_, homogeneous_);
            out << "\n";
        }
 
        // the actual coefficients below, as a named-expression-style legend (row r <-> f_k above)
        const int sig = describe_detail::CoeffSig(verbose);
        out << "    c =\n";
        for (Eigen::Index r = 0; r < shown; ++r)
        {
            out << "      [ ";
            for (Eigen::Index c = 0; c < M.cols(); ++c)
            {
                if (c) out << ", ";
                describe_detail::PrintCoeff(out, M(r, c), sig);
            }
            out << " ]\n";
        }
 
        if (shown < n)
        {
            row += static_cast<size_t>(n - shown);     // keep the global row index correct
            out << "    ... (" << (n - shown) << " more form" << (n - shown == 1 ? "" : "s")
                << "; describe(verbose=True) for all)\n";
        }
    }

    bool DependsOnPathVariable() const { return false; }

    bool HasConstantJacobian() const { return true; }

    void Differentiate() const {}

    unsigned Precision() const { return precision_; }

    void Precision(unsigned new_precision) const
    {
        // short-circuit when already materialized here.  Each holder of mp values keeps its
        // own "materialized at" tag; System deliberately keeps none and simply fans out on
        // every evaluation, which is cheap precisely because of this early return (ADR-0057).
        if (precision_==new_precision)
            return;
        if (new_precision > DoublePrecision())
        {
            auto& wm = std::get<Mat<complex_mp>>(coefficients_working_);
            for (Eigen::Index r = 0; r < wm.rows(); ++r)
                for (Eigen::Index c = 0; c < wm.cols(); ++c)
                {
                    wm(r, c).precision(new_precision);
                    if (new_precision > precision_)
                        wm(r, c) = coefficients_highest_precision_(r, c);
                }
        }
        precision_ = new_precision;
    }

    template <typename T>
    void EvalInPlace(Eigen::Ref<Vec<T>> result, Vec<T> const& vars, T const& /*path_value*/) const
    {
        SyncPrecision(vars);
        // affine: f(x) = W * [x ; 1] (the trailing 1 carries each row's constant in the last
        // column).  homogeneous (post-Homogenize): every column is a variable column, so it is
        // just W * vars (the old constant is now the homogenizing variable's coefficient).
        if (homogeneous_)
            result.noalias() = Working<T>() * vars;
        else
            result.noalias() = Working<T>() * Augment<T>(vars);
    }

    template <typename T>
    void JacobianInPlace(Eigen::Ref<Mat<T>> J, Vec<T> const& vars, T const& /*path_value*/) const
    {
        SyncPrecision(vars);
        // homogeneous: every column is d f / d x.  affine: drop the trailing constant column.
        if (homogeneous_)
            J = Working<T>();
        else
            J = Working<T>().leftCols(static_cast<Eigen::Index>(num_vars_));
    }

    template <typename T>
    void TimeDerivInPlace(Eigen::Ref<Vec<T>> result, Vec<T> const& /*vars*/, T const& /*path_value*/) const
    {
        result.setZero();
    }
 
private:
    template <typename T>
    void SyncPrecision(Vec<T> const& vars) const
    {
        if constexpr (!std::is_same<T, complex_dbl>::value)
        {
            if (vars.size() > 0)
            {
                const unsigned p = bertini::Precision(vars(0));
                if (p != precision_)
                    Precision(p);
            }
        }
    }
 
    template <typename T>
    const Mat<T>& Working() const
    {
        return std::get<Mat<T>>(coefficients_working_);
    }
 
    template <typename T>
    Vec<T> Augment(Vec<T> const& vars) const
    {
        Vec<T> aug(static_cast<Eigen::Index>(num_vars_ + 1));
        aug.head(static_cast<Eigen::Index>(num_vars_)) = vars;
        T one(1);
        if constexpr (!std::is_same<T, complex_dbl>::value)
            one.precision(precision_);
        aug(static_cast<Eigen::Index>(num_vars_)) = one;
        return aug;
    }
 
    void BuildWorking() const
    {
        const auto& M = coefficients_highest_precision_;
        auto& wd = std::get<Mat<complex_dbl>>(coefficients_working_);
        auto& wm = std::get<Mat<complex_mp>>(coefficients_working_);
        wd.resize(M.rows(), M.cols());
        wm.resize(M.rows(), M.cols());
        for (Eigen::Index r = 0; r < M.rows(); ++r)
            for (Eigen::Index c = 0; c < M.cols(); ++c)
            {
                wd(r, c) = complex_dbl(M(r, c));
                wm(r, c) = M(r, c);
            }
    }
 
    size_t num_vars_;
    bool homogeneous_ = false; 
    Mat<complex_mp> coefficients_highest_precision_; 
    mutable std::tuple<Mat<complex_dbl>, Mat<complex_mp>> coefficients_working_;
    mutable unsigned precision_;
 
    friend class boost::serialization::access;
 
    template <typename Archive>
    void serialize(Archive& ar, const unsigned /*version*/)
    {
        ar & num_vars_;
        ar & homogeneous_;
        ar & precision_;
        ar & coefficients_highest_precision_;
        ar & std::get<0>(coefficients_working_);
        ar & std::get<1>(coefficients_working_);
    }
};
 
} // namespace blocks
} // namespace bertini