tropical/src/augmented_matroid_fan.cc

/*
	This program 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 2
	of the License, or (at your option) any later version.

	This program 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 LicenseRestrictedIncidenceMatrix
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor,
	Boston, MA  02110-1301, USA.

	---
	Copyright (C) 2011 - 2015, Simon Hampe <simon.hampe@googlemail.com>

	---
	Copyright (c) 2016-2021
	Ewgenij Gawrilow, Michael Joswig, and the polymake team
	Technische Universität Berlin, Germany
	https://polymake.org

	Computes an augmented bergman matroid fan by extending the chains of the
    lattice of flats with their respective compatible independent sets.
	*/


#include "polymake/client.h"
#include "polymake/Matrix.h"
#include "polymake/IncidenceMatrix.h"
#include "polymake/Rational.h"
#include "polymake/Vector.h"
#include "polymake/Set.h"
#include "polymake/graph/Lattice.h"
#include "polymake/graph/Decoration.h"
#include "polymake/graph/maximal_chains.h"
#include "polymake/tropical/specialcycles.h"

namespace polymake { namespace tropical {

using graph::Lattice;
using graph::lattice::Sequential;
using graph::lattice::BasicDecoration;

template <typename Addition>
BigObject augmented_matroid_fan(BigObject matroid)
{
    const Int n = matroid.give("N_ELEMENTS");
    const Array<Set<Int>> bases = matroid.give("BASES");
    const Int full_rank = matroid.give("RANK");
    const Set<Int> loops = matroid.give("LOOPS");

    const BigObject lattice_flats = matroid.give("LATTICE_OF_FLATS");
    const int top_index = lattice_flats.give("TOP_NODE");
    const NodeMap<Directed,Set<Int>> flats = lattice_flats.give("FACES");

    const bool ascending = (top_index != 0); // Is FACES ordered rank-ascending, or rank-descending?

    if (loops.size() != 0) {
        throw std::runtime_error("augmented_matroid_fan: The passed matroid is not loopless!");
    }

    // Create the rays.

    // We have two homogenizing coordinates (in the columns).
    // First 0 says "this is a ray"
    // Second 0 comes from homogenization onto the tropical torus
    Matrix<Rational> rays(n + flats.size(), n + 2);

    // First type of ray corresponds to elements of the matroid's ground set.
    // These are the rays generated by all the standard basis vectors.
    // Hence we start the identity matrix (with two homogenizing coordinates).
    for (int i = 0; i < n; i++) {
        rays(i, i + 2) = Addition::orientation();
    }

    // Second type of ray corresponds to nonempty proper flats.
    // These rays are generated by the negative complement of the
    // indicator vector.
    // E.g. if the flat is {0, 2, 3} out of 5 Elements
    // Its indicator vector is    ( 1,  0,  1,  1,  0)
    // And the form we require is ( 0, -1,  0,  0, -1)

    for (int i = 0; i < flats.size(); i++) {
        if (i == top_index) { // Skip the ground set, it's not a proper flat.
            continue;
        }

        Set<Int> flat = flats[i];
        for (int j = 0; j < n; j++) {
            if(!flat.exists(j)) {
                // We already wrote first n rows previously. Hence + n
                // If we have a desecending order, we skip the first flat because its the ground set.
                const Int row_index = ascending ? i + n : i + n - 1;

                const Int col_index = j + 2; // + 2 due to the two homogenizing coords.
                rays(row_index, col_index) = -1 * Addition::orientation();
            }
        }
    }

    // Lastly, we add an artificial "ray" that corresponds to the origin in a fan.
    // This is required because a tropical cycle is in general a polyhedral complex.
    rays(rays.rows() - 1, 0) = 1;


    // Next we must compute the maximal cones of the fan.
    // These maximal cones correspond to tuples (I, F), where
    // I is an independent set in the matroid, and F is a maximal flag of proper
    // flats with the condition that every proper flat in F must contain I.

    // We generate all the maximal flags of flats without considering I
    // We want the empty face in our flags, but not the ground set - false, true.
    const Lattice<BasicDecoration, Sequential> dec_flats(lattice_flats);
    const Array<Set<Int>> maximal_flags = maximal_chains(dec_flats, false, true);

    // Output maximal cones of the fan.
    RestrictedIncidenceMatrix<> maximal_cones(0);
    Int maximal_cones_rows = 0;

    // All the maximal flags of flats give us the maximal cones of the form ({}, F)
    // These are the first maximal cones we create.
    for (int i = 0; i < maximal_flags.size(); i++) {
        const Set<Int> flag = maximal_flags[i];
        Set<Int> new_row;
        for (Int face_index : flag) {
            // + n Since we wish to adress the rays which correspond to faces.
            // - 1 for descending order, since our 0 corresponds to the "first"
            //    face which is not the ground set, but 0 corresponds to the
            //    ground set in a descending order.

            Int new_elem = ascending ? face_index + n : face_index + n - 1;
            new_row += new_elem;
        }
        maximal_cones /= new_row;
        maximal_cones_rows++;
    }


    // We now iteratively create the maximal cones for the intermediate ranks.
    // We iterate on the rank 0 < r < full_rank of the independent set I in (I, F).

    // Assume we are in the situtation that we wish to construct all (I, F)
    // with rank(I) = r, and we know how all the  (I', F') look like with rank(I') < r.

    // Take a (I', F'), with rank(I') = r - 1, we write
    // F' = { F_(r-1), F_(r), F_(r+1), ..., F_(full_rank - 1) }
    // with F_(r-1) ⊂ F_(r) ⊂ F_(r+1) ⊂ ... ⊂ F_(full_rank - 1)

    // Then these (I', F') generate new (I, F) with rank(I) = r.
    // We take F = F' - F_(r-1).
    // And I = I' + e, for any e in F_(r) \ F_(r-1).

    // Notice that we get all possible pairs (I, F) in this way if we take all
    // previous (I', F'). Let I = (x_1, ..., x_r). Then we may arrive at
    // (I,F) by considering the sequence of flats cl(x_1), cl(x_1, x_2), ..., cl(x_1, x_2, x_r)
    // By permuting these x_1, ..., x_r we arrive at dupliacates of (I,F).
    // In fact, these are exactly all the duplicates.
    // So, in order to prevent duplicates, we demand that x_1 > x_2 > ... > x_r. (viewed by their indices)
    // Exactly one permutation will satisify this property.

    // These two variables indicate the rows of the matrix we have
    // written in the previous outermost iteration.
    // I.e. they correspond to the (I', F') from which we will generate
    // the next iteration's (I, F).
    Int last_rows_from = 0;
    Int last_rows_to = maximal_cones_rows; // Notice we begin with the already generated rank 0 cones.

    for (int rank = 1; rank < full_rank; rank++) {
        // Iterate through all (I', F')
        for (Int row_index = last_rows_from; row_index < last_rows_to; row_index++) {

            const Set<Int> the_row = maximal_cones.row(row_index);
            const Set<Int> i_prime = the_row * range(0, n - 1); // Corresponds to I'.
            const Set<Int> f_prime = the_row - i_prime;         // Corresponds to F' (in indices).

            Set<Int> f = f_prime;   // Will correspond to F after if branch.
            Set<Int> smallest_face; // Will correspond to F_(r-1) after if branch.
            Set<Int> second_smallest_face; // Will correspond to F_(r) after if branch.
            if (ascending) { // Smallest face is indexed by smallest integer.
                Int index = f.front(); f.pop_front();
                // Need to correct index.
                // Must subtract by n, since the first n elements correspond to matroid elements in the_row.
                index -= n;
                smallest_face = flats[index];
                index = f.front();
                index -= n; // Ditto with the correction.
                second_smallest_face = flats[index];
            } else { // Smallest face is indexed by largest integer.
                Int index = f.back(); f.pop_back();
                // Need to correct the index.
                // Must subtract by n, since the first n elements correspond to matroid elements in the_row.
                // 0 refers to ground set in flats;
                // But 0 refers to "first" set which is not ground set in maximal_cones after - n.
                // So we add 1 back.
                index = index - n + 1;
                smallest_face = flats[index];
                index = f.back();
                index = index - n + 1; // Ditto with the correction.
                second_smallest_face = flats[index];
            }

            // We now have constructed F, F_(r) and F_(r-1).
            // We can now easily itearte over all e with e in F_(r) \ F_(r-1)
            // These e allow us to construct all possible I = I' + e.
            // This is all we need to generate another row.
            Set<Int> all_es = second_smallest_face - smallest_face;
            // The smallest value of I'. If I' is empty, then effectively no min value.
            const Int min_value = i_prime.size() == 0 ? n + 1 : i_prime.front();

            for (const Int e : all_es) {
                if (e > min_value) {
                    // To create (I, F) with I = (x_1, ..., x_r)
                    // Where x_1 was added first, then x_2, and so on.
                    // We only allow the creation via x_1 > ... > x_r.
                    // This prevents duplicates.
                    break;
                }

                Set<Int> new_row = i_prime + scalar2set(e) + f;
                maximal_cones /= new_row;
                maximal_cones_rows++;
            }
        }

        last_rows_from = last_rows_to;
        last_rows_to = maximal_cones_rows;
    }

    // Now add the maximal cones which correspond to bases
    // I.e. they are of the form (B, {}) where B is a basis.
    IncidenceMatrix<> bases_matrix(bases);
    maximal_cones /= bases_matrix;
    maximal_cones_rows += bases_matrix.rows();

    // Since a cycle is in general a polyhedral complex we add
    // the zero as a vertex of all the maximal cones.
    for (int i = 0; i < maximal_cones_rows; i++) {
        maximal_cones(i, rays.rows() - 1) = 1;
    }

    return BigObject("Cycle", mlist<Min>(),
                     "PROJECTIVE_VERTICES", rays,
                     "MAXIMAL_POLYTOPES", IncidenceMatrix<>(std::move(maximal_cones)),
                     "WEIGHTS", ones_vector<Integer>(maximal_cones_rows));
}


UserFunctionTemplate4perl("# @category Matroids"
                          "# Computes the augmented Bergman fan of a matroid."
                          "# Note that this is potentially very slow for large matroids."
                          "# A definition of the augmented Bergman fan can be found in arXiv:2002.03341. See also the follow up paper arXiv:2010.06088."
                          "# The algorithim used to construct the augemented Bergman fan closely follows its description in the first paper."
                          "# @param matroid::Matroid A matroid. Should be loopfree."
                          "# @tparam Addition Min or max, determines the matroid fan coordinates."
                          "# @example [application matroid]"
                          "# > $m = matroid::fano_matroid;"
                          "# > $f = tropical::augmented_matroid_fan<Min>($m);"
                          "# @return tropical::Cycle<Addition>",
                          "augmented_matroid_fan<Addition>(matroid::Matroid)");

} }