slice.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//slice.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.
//
//slice.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 slice.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.

 
 
#ifndef BERTINI_SLICE_HPP
#define BERTINI_SLICE_HPP
 
#include <vector>
 
#include "bertini2/function_tree.hpp"
#include "bertini2/function_tree/reintern.hpp"
#include "bertini2/num_traits.hpp"
#include "bertini2/eigen_extensions.hpp"
#include "bertini2/system/blocks/linear_forms_block.hpp"
 
 
namespace bertini {
 
    // Slice::AddTo hands the slice's linear-forms block to a System; we only need System's
    // name here (the definition lives in slice.cpp, which includes system.hpp).
    class System;

    class Slice
    {
        blocks::LinearFormsBlock block_;  
        VariableGroup sliced_vars_;       
        bool is_homogeneous_ = false;     
 
    public:

        Slice() = default;

        static Slice FromCoefficients(VariableGroup const& v, Mat<complex_mp> const& augmented_coefficients, bool homogeneous = false)
        {
            assert(static_cast<size_t>(augmented_coefficients.cols()) == v.size() + 1 &&
                   "a slice coefficient matrix must have (num_variables + 1) columns");
            Slice s;
            s.sliced_vars_ = v;
            s.is_homogeneous_ = homogeneous;
            s.block_ = blocks::LinearFormsBlock(v.size(), augmented_coefficients);
            return s;
        }

        static Slice ThroughPoint(VariableGroup const& v, Mat<complex_mp> const& coefficients, Vec<complex_mp> const& point)
        {
            if (static_cast<size_t>(coefficients.cols()) != v.size())
                throw std::runtime_error("Slice::ThroughPoint coefficient block must have num_variables columns (it is not augmented)");
            if (static_cast<size_t>(point.size()) != v.size())
                throw std::runtime_error("Slice::ThroughPoint point must have num_variables entries");
 
            return FromCoefficients(v, AugmentAffineThroughPoint(coefficients, point), /*homogeneous=*/false);
        }

        static Slice RandomReal(VariableGroup const& v, unsigned dim, bool homogeneous = false, bool orthogonal = true, Vec<complex_mp> const* through_point = nullptr)
        {
            typedef void (*funtype) (complex_mp&, unsigned); // the type for number generation
            // bounded-modulus draw (away from 0 and infinity), matching patches and the start systems;
            // kept REAL so a real slice stays real.  The orthogonal path is real too (real=true below):
            // it QR-factors a matrix of REAL units, yielding a real orthogonal coefficient block (issue
            // #294 -- previously the orthogonal path hardcoded the complex orthonormal matrix, so a
            // "real" slice came out complex).
            funtype gen = bertini::multiprecision::RandomRealBoundedModulusAssign;
            return Make(v, dim, homogeneous, orthogonal, /*real=*/true, gen, through_point);
        }

        static Slice RandomComplex(VariableGroup const& v, unsigned dim, bool homogeneous = false, bool orthogonal = true, Vec<complex_mp> const* through_point = nullptr)
        {
            typedef void (*funtype) (complex_mp&, unsigned); // the type for number generation
            // bounded-modulus draw (away from 0 and infinity), matching patches and the start systems.
            funtype gen = bertini::multiprecision::RandomComplexBoundedModulusAssign;
            return Make(v, dim, homogeneous, orthogonal, /*real=*/false, gen, through_point);
        }

        static Slice Make(VariableGroup const& v, unsigned dim, bool homogeneous, bool orthogonal, bool real, std::function<void(complex_mp&, unsigned)> gen, Vec<complex_mp> const* through_point = nullptr)
        {
            const unsigned num_vars = static_cast<unsigned>(v.size());
 
            // more forms than variables cannot be independent: the orthonormalization below would
            // hand back dependent rows without a word, and the slice would cut nothing new (b2#380)
            if (dim > num_vars)
                throw std::invalid_argument("Slice: " + std::to_string(dim) + " linear forms on "
                    + std::to_string(num_vars) + " variables; a slice cannot have more forms than variables");
            if (homogeneous && through_point && dim + 1 > num_vars)
                throw std::invalid_argument("Slice: " + std::to_string(dim) + " homogeneous linear forms through a point on "
                    + std::to_string(num_vars) + " variables; the forms live in the point's orthogonal complement, which has "
                    + std::to_string(num_vars - 1) + " dimensions");
 
            if (through_point && static_cast<unsigned>(through_point->size()) != num_vars)
                throw std::runtime_error("Slice: through_point must have num_variables entries");
 
            Mat<complex_mp> coeffs(dim, num_vars); // the variable coefficients (one row per form)
 
            if (orthogonal)
            {
                // conjugate-orthonormal coefficient matrix (orthonormal linear forms), drawn the b1 way
                // via RandomConjugateOrthonormalMatrix (ADR-0041) -- it generates square and truncates,
                // so the old transpose dance is gone.  Built at max precision, like the rest of the slice.
                // A real slice QR-factors a matrix of REAL units (a real orthogonal block); a complex
                // slice uses complex units.  (Issue #294: the real path must stay real.)
                auto prev_precision = DefaultPrecision();
                DefaultPrecision(MaxPrecisionAllowed());
                if (real)
                    coeffs = bertini::RandomConjugateOrthonormalMatrix<real_mp>(dim, num_vars)
                                 .unaryExpr([](real_mp const& r){ return complex_mp(r); });
                else
                    coeffs = bertini::RandomConjugateOrthonormalMatrix<complex_mp>(dim, num_vars);
                DefaultPrecision(prev_precision);
            }
            else
            {
                for (unsigned ii(0); ii < dim; ++ii)
                    for (unsigned jj(0); jj < num_vars; ++jj)
                        gen(coeffs(ii, jj), MaxPrecisionAllowed());
            }
 
            assert(static_cast<unsigned>(coeffs.rows()) == dim);
            assert(static_cast<unsigned>(coeffs.cols()) == num_vars);
 
            // Through a point.  Homogeneous: project the rows into p's orthogonal complement (a.p = 0),
            // zero constant column.  Affine: the constant column is -A*p, so every form vanishes at p.
            if (through_point && homogeneous)
                return FromCoefficients(v, AugmentHomogeneousThroughPoint(coeffs, *through_point, orthogonal), /*homogeneous=*/true);
            if (through_point)
                return FromCoefficients(v, AugmentAffineThroughPoint(coeffs, *through_point), /*homogeneous=*/false);
 
            // Assemble the augmented matrix: [ coeffs | constants ].  A homogeneous slice's constant
            // column is zero; otherwise it is freshly generated.
            Mat<complex_mp> augmented(dim, num_vars + 1);
            augmented.leftCols(num_vars) = coeffs;
            if (homogeneous)
                augmented.col(num_vars).setZero();
            else
                for (unsigned ii(0); ii < dim; ++ii)
                    gen(augmented(ii, num_vars), MaxPrecisionAllowed());
 
            return FromCoefficients(v, augmented, homogeneous);
        }
 

        template<typename NumT>
        void Eval(Vec<NumT> & result, Vec<NumT> const& x) const
        {
            result.resize(Dimension());
            NumT path_value(0); // ignored: linear forms are autonomous
            block_.EvalInPlace<NumT>(result, x, path_value);
        }

        template<typename NumT>
        Vec<NumT> Eval(Vec<NumT> const& x) const
        {
            Vec<NumT> result(Dimension());
            Eval(result, x);
            return result;
        }

        template<typename NumT>
        void Jacobian(Mat<NumT> & result, Vec<NumT> const& x) const
        {
            result.resize(Dimension(), NumVariables());
            NumT path_value(0);
            block_.JacobianInPlace<NumT>(result, x, path_value);
        }

        template<typename NumT>
        Mat<NumT> Jacobian(Vec<NumT> const& x) const
        {
            Mat<NumT> result(Dimension(), NumVariables());
            Jacobian(result, x);
            return result;
        }
 

        Mat<complex_mp> const& Coefficients() const
        {
            return block_.Coefficients();
        }

        blocks::LinearFormsBlock const& AsLinearFormsBlock() const
        {
            return block_;
        }

        void AddTo(System & s) const;

        System AsSystem() const;

        Slice Concatenate(Slice const& other) const
        {
            if (NumVariables() != other.NumVariables())
                throw std::runtime_error("Slice::Concatenate requires both slices to be on the same number of variables");
 
            Mat<complex_mp> const& A = Coefficients();
            Mat<complex_mp> const& B = other.Coefficients();
            Mat<complex_mp> stacked(A.rows() + B.rows(), A.cols());
            stacked.topRows(A.rows()) = A;
            stacked.bottomRows(B.rows()) = B;
 
            return FromCoefficients(sliced_vars_, stacked, is_homogeneous_ && other.is_homogeneous_);
        }
 

        Slice Head(unsigned m) const
        {
            if (m > Dimension())
                throw std::runtime_error("Slice::Head asked for more forms than the slice has");
            return FromCoefficients(sliced_vars_, Coefficients().topRows(m), is_homogeneous_);
        }

        Slice Tail(unsigned m) const
        {
            if (m > Dimension())
                throw std::runtime_error("Slice::Tail asked for more forms than the slice has");
            return FromCoefficients(sliced_vars_, Coefficients().bottomRows(m), is_homogeneous_);
        }

        Slice Rows(std::vector<unsigned> const& indices) const
        {
            Mat<complex_mp> const& C = Coefficients();
            Mat<complex_mp> sub(static_cast<Eigen::Index>(indices.size()), C.cols());
            for (size_t ii = 0; ii < indices.size(); ++ii)
            {
                if (indices[ii] >= Dimension())
                    throw std::runtime_error("Slice::Rows asked for a form index outside the slice");
                sub.row(static_cast<Eigen::Index>(ii)) = C.row(indices[ii]);
            }
            return FromCoefficients(sliced_vars_, sub, is_homogeneous_);
        }
 

        unsigned Dimension() const
        {
            return static_cast<unsigned>(block_.NumFunctions());
        }

        unsigned NumVariables() const
        {
            return static_cast<unsigned>(sliced_vars_.size());
        }

        VariableGroup const& Variables() const
        {
            return sliced_vars_;
        }

        bool IsHomogeneous() const
        {
            return is_homogeneous_;
        }
 

        unsigned Precision() const
        {
            return block_.Precision();
        }

        void Precision(unsigned new_precision) const
        {
            block_.Precision(new_precision);
        }
 
    private:

        static Mat<complex_mp> AugmentAffineThroughPoint(Mat<complex_mp> const& coeffs, Vec<complex_mp> const& point)
        {
            const auto num_vars = coeffs.cols();
 
            auto prev_precision = DefaultPrecision();
            DefaultPrecision(MaxPrecisionAllowed());
 
            Mat<complex_mp> A = coeffs;   // copies, re-materialized at max precision below
            Vec<complex_mp> p = point;
            bertini::Precision(A, MaxPrecisionAllowed());   // qualify: class has member Precision overloads
            bertini::Precision(p, MaxPrecisionAllowed());
 
            Mat<complex_mp> augmented(coeffs.rows(), num_vars + 1);
            augmented.leftCols(num_vars) = A;
            augmented.col(num_vars) = -(A * p);   // b = -A*p, so A*p + b = 0
 
            DefaultPrecision(prev_precision);
            return augmented;
        }

        static Mat<complex_mp> AugmentHomogeneousThroughPoint(Mat<complex_mp> const& coeffs, Vec<complex_mp> const& point, bool orthogonal)
        {
            const auto dim = coeffs.rows();
            const auto num_vars = coeffs.cols();
 
            auto prev_precision = DefaultPrecision();
            DefaultPrecision(MaxPrecisionAllowed());
 
            Mat<complex_mp> A = coeffs;
            Vec<complex_mp> p = point;
            bertini::Precision(A, MaxPrecisionAllowed());   // qualify: class has member Precision overloads
            bertini::Precision(p, MaxPrecisionAllowed());
 
            real_mp p_norm_sq = p.squaredNorm();            // pbar.p = ||p||^2, real-positive for p != 0
            if (p_norm_sq == real_mp(0))
                throw std::runtime_error("Slice: a homogeneous slice through a point needs a nonzero point");
 
            // project every row into p's orthogonal complement (bilinear a.p = 0): subtract the
            // rank-one correction (A*p) pbar^T / ||p||^2.
            Vec<complex_mp> Ap = A * p;                      // bilinear a_i . p, one per row
            A = A - (Ap / complex_mp(p_norm_sq)) * p.conjugate().transpose();
 
            if (orthogonal)
            {
                // orthonormalize the projected rows via a QR of A^T: the columns of Q span the row
                // space of A (which lies in p's complement), so Q^T's rows stay orthogonal to p.
                Eigen::HouseholderQR<Mat<complex_mp>> qr(A.transpose());
                Mat<complex_mp> Q = qr.householderQ() * Mat<complex_mp>::Identity(num_vars, num_vars);
                A = Q.leftCols(dim).transpose();
            }
 
            Mat<complex_mp> augmented(dim, num_vars + 1);
            augmented.leftCols(num_vars) = A;
            augmented.col(num_vars).setZero();              // homogeneous: no constant term
 
            DefaultPrecision(prev_precision);
            return augmented;
        }
 
        friend class boost::serialization::access;
 
        template <typename Archive>
        void serialize(Archive& ar, const unsigned /*version*/) {
            ar & block_;
            ar & sliced_vars_;
            ar & is_homogeneous_;
 
            // a loaded slice's variables join the live node universe, as a loaded System's do
            // (ADR-0068): otherwise its system of linear forms could not be combined with the
            // system it slices
            if (Archive::is_loading::value)
                for (auto& v : sliced_vars_)
                    if (v)
                        v = std::static_pointer_cast<node::Variable>(node::Reintern(v));
        }
 
        friend std::ostream& operator<<(std::ostream&, Slice const&);
    };

    std::ostream& operator<<(std::ostream& out, Slice const& s);
} // re: namespace bertini
 
#endif