#  Copyright (c) 1997-2024
#  Ewgenij Gawrilow, Michael Joswig, and the polymake team
#  Technische Universität Berlin, Germany
#  https://polymake.org
#
#  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, or (at your option) any
#  later version: http://www.gnu.org/licenses/gpl.txt.
#
#  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 PURPSE.  See the
#  GNU General Public License for more details.
#-------------------------------------------------------------------------------

object Matroid {

rule RANK : BASES {
  $this->RANK = $this->BASES->[0]->size();
}
weight 0.1;

rule BASES, N_BASES, RANK, N_ELEMENTS : VECTORS {
  bases_from_points($this)
}
weight 4.50;

rule BASES, N_BASES, RANK : LATTICE_OF_FLATS, N_ELEMENTS {
  bases_from_lof($this)
}
weight 4.50;

rule POLYTOPE.CONE_AMBIENT_DIM : N_ELEMENTS {
   $this->POLYTOPE->CONE_AMBIENT_DIM = $this->N_ELEMENTS+1;
}
weight 0.10;

rule POLYTOPE.FEASIBLE : N_BASES {
   $this->POLYTOPE->FEASIBLE = $this->N_BASES > 0;
}
weight 0.10;

rule POLYTOPE.BOUNDED : {
   $this->POLYTOPE->BOUNDED = true;
}
weight 0.10;

rule POLYTOPE.VERTICES, POLYTOPE.CONE_AMBIENT_DIM, POLYTOPE.FEASIBLE, POLYTOPE.BOUNDED : BASES, N_ELEMENTS {
   polytope::matroid_polytope($this)
}
weight 1.12;

rule POLYTOPE.VertexPerm.PERMUTATION = BasesPerm.PERMUTATION;

rule POLYTOPE.VERTICES, POLYTOPE.CONE_AMBIENT_DIM, POLYTOPE.INEQUALITIES, POLYTOPE.EQUATIONS, POLYTOPE.FEASIBLE, POLYTOPE.BOUNDED: BASES, N_ELEMENTS, RANK, CONNECTED, LATTICE_OF_FLATS {
   polytope::matroid_polytope($this, inequalities=>true)
}
precondition : CONNECTED;
weight 1.10;

rule POLYTOPE.EHRHART_POLYNOMIAL : SPARSE_PAVING, N_ELEMENTS, RANK, N_BASES {
	my $n = $this->N_ELEMENTS;
	my $r = $this->RANK;
	my $lambda = binomial($n,$r)-$this->N_BASES;
	if($lambda == 0){
		$this->POLYTOPE->EHRHART_POLYNOMIAL = polytope::ehrhart_polynomial_hypersimplex($r,$n);
	}else{
		$this->POLYTOPE->EHRHART_POLYNOMIAL = polytope::ehrhart_polynomial_hypersimplex($r,$n) + $lambda*(polytope::ehrhart_polynomial_product_simplicies($r,$n-$r) - polytope::ehrhart_polynomial_minimal_matroid($r,$n));
	}
}
precondition : CONNECTED && SPARSE_PAVING;
weight 1.20;

rule POLYTOPE.EHRHART_POLYNOMIAL : SPLIT_FLACETS, N_ELEMENTS, RANK {
	my $n = $this->N_ELEMENTS;
	my $r = $this->RANK;
	my $f = polytope::ehrhart_polynomial_hypersimplex($r,$n);
	for my $i (1..$this->SPLIT_FLACETS->size-1){
		for my $flat (@{$this->SPLIT_FLACETS->[$i]}){
			$f-= polytope::ehrhart_polynomial_cuspidal_matroid($r,$n,$flat->size,$r-$i)-polytope::ehrhart_polynomial_hypersimplex($i,$flat->size)*polytope::ehrhart_polynomial_hypersimplex($r-$i,$n-$flat->size);
		}
	}
	$this->POLYTOPE->EHRHART_POLYNOMIAL = $f;
}
precondition : CONNECTED && SPLIT;
weight 1.50;


rule RANK : NON_BASES { # should be more expensive than the computation from BASES
   $this->RANK=$this->NON_BASES->[0]->size();
}
precondition : NON_BASES {
   $this->NON_BASES->size() > 0;
}
weight 0.1;

rule BASES : RANK, N_ELEMENTS {
   $this->BASES = all_subsets_of_k(sequence(0, $this->N_ELEMENTS), $this->RANK);
}
precondition : NON_BASES {
   $this->NON_BASES->size() == 0;
}
weight 0.10;

rule BASES : NON_BASES, N_ELEMENTS {
   $this->BASES = invert_bases($this->NON_BASES, $this->N_ELEMENTS);
}
precondition : NON_BASES {
   $this->NON_BASES->size() > 0;
}
weight 3.10;

rule NON_BASES : BASES, N_ELEMENTS {
  $this->NON_BASES = invert_bases($this->BASES, $this->N_ELEMENTS)
}
weight 3.10;

rule CIRCUITS : BASES, N_ELEMENTS {
   $this->CIRCUITS = bases_to_circuits($this->BASES, $this->N_ELEMENTS);
}
weight 4.10;

rule BASES : CIRCUITS, N_ELEMENTS {
   $this->BASES = circuits_to_bases($this->CIRCUITS, $this->N_ELEMENTS);
}
weight 4.10;

rule BASES : CIRCUITS, N_ELEMENTS, RANK {
   $this->BASES = circuits_to_bases_rank($this->CIRCUITS, $this->N_ELEMENTS, $this->RANK);
}
weight 3.10;

rule MATROID_HYPERPLANES : CIRCUITS, N_ELEMENTS, RANK {
   $this->MATROID_HYPERPLANES = circuits_to_hyperplanes($this->CIRCUITS, $this->N_ELEMENTS, $this->RANK);
}
weight 3.10;

rule N_MATROID_HYPERPLANES : MATROID_HYPERPLANES {
    $this->N_MATROID_HYPERPLANES = $this->MATROID_HYPERPLANES->size();
}
weight 0.10;

rule LOOPS : BASES, N_ELEMENTS {
    loops($this);
}
weight 1.10;

rule N_BASES : BASES {
   $this->N_BASES = $this->BASES->size();
}
weight 0.1;

rule N_CIRCUITS : CIRCUITS {
   $this->N_CIRCUITS = $this->CIRCUITS->size();
}
weight 0.1;

rule N_LOOPS : LOOPS {
   $this->N_LOOPS = $this->LOOPS->size();
}
weight 0.1;


rule BASES, MAXIMAL_TRANSVERSAL_PRESENTATION : N_ELEMENTS, TRANSVERSAL_PRESENTATION {
   #Create bipartite graph from edges
   # Elements from the ground set are nodes 0,..,n-1.
   # Sets from the set system are nodes n,..,n+k
   my @edges = ();
   my $n = $this->N_ELEMENTS;
   my $tp = $this->TRANSVERSAL_PRESENTATION;
   for my $i (0 .. $tp->size()-1) {
      for my $x (@{$tp->[$i]}) {
         push @edges, new Set<Int>($x, $n + $i);
      }
   }
   my $g = graph_from_edges(\@edges);
   my $p = polytope::fractional_matching_polytope($g->ADJACENCY);
   $p->LP = new polytope::LinearProgram(LINEAR_OBJECTIVE=>ones_vector<Rational>(scalar(@edges)+1));
   my $maxverts = new Matrix<Rational>($p->VERTICES->minor($p->LP->MAXIMAL_FACE,~[0]));
   my @basislist = map { new Set<Int>(map {  min(@{$g->EDGES->[$_]}) } @{support($_)}) } @{rows($maxverts)};
   #Need to make the list of bases irredundant, some might occur multiple times
   $this->BASES = new Array<Set<Int>>( new Set<Set<Int>>(\@basislist));
   if($maxverts->rows() > 0) {
      my $one_matching = new Set<Int>( map { max(@{$g->EDGES->[$_]}) - $n} @{support($maxverts->row(0))});
      $this->MAXIMAL_TRANSVERSAL_PRESENTATION =
      maximal_transversal_presentation($n,$this->BASES,$tp, $one_matching);
   }
   else {
      $this->MAXIMAL_TRANSVERSAL_PRESENTATION = new IncidenceMatrix(0,$this->N_ELEMENTS);
   }
}
weight 2.10;
precondition : defined(TRANSVERSAL_PRESENTATION);

rule TRANSVERSAL : {
   $this->TRANSVERSAL = 1;
}
weight 0.10;
precondition : defined(TRANSVERSAL_PRESENTATION | MAXIMAL_TRANSVERSAL_PRESENTATION);

rule TRANSVERSAL, TRANSVERSAL_PRESENTATION : N_ELEMENTS, RANK, LOOPS, LATTICE_OF_FLATS, LATTICE_OF_CYCLIC_FLATS {
   my @a = check_transversality($this);
   $this->TRANSVERSAL = $a[0];
   if($a[0]) {
      $this->TRANSVERSAL_PRESENTATION = $a[1];
   }
   else {
      $this->TRANSVERSAL_PRESENTATION = new Array<Set<Int>>();
   }
}
weight 2.10;

rule TRANSVERSAL, MAXIMAL_TRANSVERSAL_PRESENTATION : LATTICE_OF_CYCLIC_FLATS {
   $this->TRANSVERSAL = 1;
   $this->MAXIMAL_TRANSVERSAL_PRESENTATION = nested_presentation($this);
}
weight 1.10;
precondition : NESTED;

rule NESTED : LATTICE_OF_CYCLIC_FLATS {
   $this->NESTED = is_nested_matroid($this);
}
weight 1.10;

rule SIMPLE : CIRCUITS {
    $this->SIMPLE = 1;
    foreach (@{$this->CIRCUITS}){
        if($_->size() < 3){
            $this->SIMPLE = 0;
            last;
        }
    }
}
weight 1.10;

rule PAVING : CIRCUITS, RANK {
    $this->PAVING=( $this->CIRCUITS->size() == 0 or $this->RANK<=minimum([map{$_->size()}@{$this->CIRCUITS}]));
}
weight 1.10;



rule BINARY, BINARY_VECTORS : BASES, N_ELEMENTS, RANK {
    binary_representation($this);
}
weight 2.10;

rule TERNARY, TERNARY_VECTORS : BASES, N_ELEMENTS, RANK {
    ternary_representation($this);
}
weight 5.10;

rule REGULAR : BINARY, TERNARY {
    $this->REGULAR=$this->BINARY && $this->TERNARY; # See [Oxley:Matroid theory (2nd ed.) Thm. 6.6.3]
}
weight 0.1;

rule VECTORS : BINARY_VECTORS, TERNARY_VECTORS {
   #Normalize such that entries 0,1(,2) are used for binary/ternary rep
   my $normalized_ternary = new Matrix( [ map { [map {$_ % 3} @{$_}]} @{rows($this->TERNARY_VECTORS)}]);
   my $normalized_binary  = new Matrix( [ map { [map {$_ % 2} @{$_}]} @{rows($this->BINARY_VECTORS)}]);
   $this->VECTORS=3*($normalized_binary)-2*($normalized_ternary);
}
precondition : REGULAR;
weight 1.10;

rule TUTTE_POLYNOMIAL: CIRCUITS, N_ELEMENTS {
    $this->TUTTE_POLYNOMIAL = tutte_polynomial_from_circuits( $this->N_ELEMENTS,
                                                              $this->CIRCUITS);
}
weight 5.10;


rule BETA_INVARIANT : TUTTE_POLYNOMIAL {
    my $beta=0;
    for(my $i=0;$i<$this->TUTTE_POLYNOMIAL->monomials_as_matrix->rows();++$i){
        if($this->TUTTE_POLYNOMIAL->monomials_as_matrix->row($i)==new Vector("1 0")){
           $beta = new Integer($this->TUTTE_POLYNOMIAL->coefficients_as_vector->[$i]);
           last;
        }
    }
    $this->BETA_INVARIANT=$beta;
}
weight 1.10;

rule SERIES_PARALLEL, CONNECTED : BETA_INVARIANT { # See [David Speyer:tropical linear spaces prop. 6.3]
   if ($this->BETA_INVARIANT==0) {
      $this->CONNECTED=0;
      $this->SERIES_PARALLEL=0;
   } else {
      $this->CONNECTED=1;
      if ($this->BETA_INVARIANT==1) {
         $this->SERIES_PARALLEL=1;
      } else {
         $this->SERIES_PARALLEL=0;
      }
   }
}
weight 0.1;
precondition : N_ELEMENTS {
   $this->N_ELEMENTS >= 2;
}

rule SERIES_PARALLEL : {
   $this->SERIES_PARALLEL = $this->N_ELEMENTS;
}
weight 0.10;
precondition : N_ELEMENTS {
   $this->N_ELEMENTS <= 1;
}

rule UNIFORM : N_ELEMENTS, RANK, N_BASES {
   $this->UNIFORM = ( $this->N_BASES == binomial( $this->N_ELEMENTS, $this->RANK));
}
weight 0.10;

rule LATTICE_OF_FLATS.ADJACENCY, LATTICE_OF_FLATS.DECORATION, LATTICE_OF_FLATS.INVERSE_RANK_MAP, LATTICE_OF_FLATS.TOP_NODE, LATTICE_OF_FLATS.BOTTOM_NODE : MATROID_HYPERPLANES, N_ELEMENTS, RANK {
	my $hy;
	#A rank-0-matroid has no hyperplanes, so we need to treat this separately
	if($this->MATROID_HYPERPLANES->size() == 0) {
		$hy = new IncidenceMatrix(0,$this->N_ELEMENTS);
	}
	else {
		#If the empty set is a hyperplane, we need to tell the incidence matrix its dimension.
		$hy=new IncidenceMatrix<NonSymmetric>(@{$this->MATROID_HYPERPLANES}, $this->N_ELEMENTS);
	}
	$this->LATTICE_OF_FLATS = lattice_of_flats($hy, $this->RANK);
}
weight 6.20;

rule MATROID_HYPERPLANES, N_ELEMENTS, RANK : LATTICE_OF_FLATS {
   my $lof = $this->LATTICE_OF_FLATS;
   $this->N_ELEMENTS = $lof->FACES->[$lof->TOP_NODE]->size();
   my $rank = $this->LATTICE_OF_FLATS->rank();
   $this->RANK = max(0,$rank);
   if($rank  <= 0 ) {
      $this->MATROID_HYPERPLANES = new Array<Set<Int>>();
   }
   else {
      $this->MATROID_HYPERPLANES = new Array<Set<Int>> ([
         map {$lof->FACES->[$_]} @{$lof->nodes_of_rank($rank-1)}
         ]);
   }
}
weight 1.10;
incurs HyperplanePerm;

rule N_FLATS : LATTICE_OF_FLATS {
   $this->N_FLATS = $this->LATTICE_OF_FLATS->N_NODES;
}
weight 0.10;

rule LATTICE_OF_CYCLIC_FLATS.ADJACENCY, LATTICE_OF_CYCLIC_FLATS.DECORATION, LATTICE_OF_CYCLIC_FLATS.INVERSE_RANK_MAP, LATTICE_OF_CYCLIC_FLATS.TOP_NODE, LATTICE_OF_CYCLIC_FLATS.BOTTOM_NODE : LATTICE_OF_FLATS, CIRCUITS, RANK {
   $this->LATTICE_OF_CYCLIC_FLATS = lattice_of_cyclic_flats($this);
}
weight 6.20;

rule N_CYCLIC_FLATS : LATTICE_OF_CYCLIC_FLATS {
   $this->N_CYCLIC_FLATS = $this->LATTICE_OF_CYCLIC_FLATS->N_NODES;
}
weight 0.10;

rule BASES, RANK : LATTICE_OF_CYCLIC_FLATS, N_ELEMENTS {
   my $loc = $this->LATTICE_OF_CYCLIC_FLATS;
   $this->RANK = $loc->rank() + $this->N_ELEMENTS - $loc->DECORATION->[$loc->TOP_NODE]->face->size();
   $this->BASES = bases_from_cyclic_flats($this->N_ELEMENTS, $this->RANK, $loc);
}
weight 6.10;

rule SPLIT_FLACETS : N_ELEMENTS, RANK, POLYTOPE.FACETS, POLYTOPE.AFFINE_HULL, CONNECTED_COMPONENTS {
    split_flacets($this);
}
weight 2.10;

rule SPLIT : SPLIT_FLACETS, CONNECTED_COMPONENTS, CIRCUITS {
    $this->SPLIT=split_compatibility_check($this);
}
weight 2.10;

rule POSITROID : SPLIT_FLACETS, CONNECTED_COMPONENTS {
       my $bool = 1;
       foreach my $S (@{$this->SPLIT_FLACETS}) {
               foreach my $flat (@{$S}){
                       foreach my $comp (@{$this->CONNECTED_COMPONENTS}) {
                               my $set = $comp*$flat;
                               next if( $set->size==0 );
                               my $count = -1;
                               my $pre = -1;
                               for my $i (@{$comp}){
                                       if( $count==-1 && !$set->contains($i) ){
                                               $pre = $i;
                                               $count = 0;
                                               next;
                                       }
                                       if( $set->contains($i) && !$set->contains($pre) ){
                                               ++$count;
                                       }
                                       if( $i==$comp->back && !$set->contains($i) && $set->contains($comp->front)){
                                               ++$count;
                                       }
                                       if($count==2 ){
                                               $bool = 0;
                                               last;
                                       }
                                       $pre  =$i;
                               }
                               last if(!$bool);
                       }
                       last if(!$bool);
               }
       }

       $this->POSITROID = $bool;
}



rule CONNECTED_COMPONENTS : CIRCUITS, N_ELEMENTS{
    $this->CONNECTED_COMPONENTS = connected_components_from_circuits($this->CIRCUITS, $this->N_ELEMENTS);
}
weight 1.10;



rule CONNECTED : N_CONNECTED_COMPONENTS {
   $this->CONNECTED = ($this->N_CONNECTED_COMPONENTS <= 1);
}
weight 0.1;

rule N_CONNECTED_COMPONENTS : CONNECTED_COMPONENTS {
   $this->N_CONNECTED_COMPONENTS = $this->CONNECTED_COMPONENTS->size();
}
weight 0.1;

rule LAMINAR : CIRCUITS, LATTICE_OF_FLATS {
   $this->LAMINAR = is_laminar_matroid($this);
}
weight 2.10;

rule H_VECTOR : BASES, RANK {
   my $bases = $this->BASES;
   my $r = $this->RANK;
   my $h_vector = new Vector<Integer>($r+1);
   foreach (@{$bases}) {
      $h_vector->[$r-internal_activity($_, $bases)]++;
   }
   $this->H_VECTOR=$h_vector;
}

rule F_VECTOR : H_VECTOR {
    my $h = $this->H_VECTOR;
    my $d = $h->dim() - 1;
    my $f = new Vector<Integer>($d);
    for (my $j=1; $j<=$d; $j++) {
        for (my $i=0; $i<=$j; $i++) {
            $f->[$j-1] += binomial($d-$i, $j-$i) * $h->[$i];
        }
    }
    $this->F_VECTOR = $f;
}

rule H_VECTOR : F_VECTOR {
    my $f = $this->F_VECTOR;
    my $d = $f->dim();
    my $h = new Vector<Integer>($d+1);
    $h->[0] = 1;
    for (my $j=1; $j<=$d; $j++) {
        for (my $i=0; $i<=$j; $i++) {
            my $c = binomial($d-$i, $j-$i);
            if ($i>0) {
                $c *= $f->[$i-1];
            }
            if (($j-$i)%2 == 0) {
                $h->[$j] += $c;
            } else {
                $h->[$j] -= $c;
            }
        }
    }
    $this->H_VECTOR = $h;
}

rule REVLEX_BASIS_ENCODING : BASES, RANK, N_ELEMENTS {
    $this->REVLEX_BASIS_ENCODING = bases_to_revlex_encoding($this->BASES, $this->RANK, $this->N_ELEMENTS);
}

rule BASES : REVLEX_BASIS_ENCODING, RANK, N_ELEMENTS {
    $this->BASES = bases_from_revlex_encoding($this->REVLEX_BASIS_ENCODING, $this->RANK, $this->N_ELEMENTS);
}

rule N_AUTOMORPHISMS = AUTOMORPHISM_GROUP.ORDER;

rule CATENARY_G_INVARIANT : LATTICE_OF_FLATS, RANK {
   $this->CATENARY_G_INVARIANT = catenary_g_invariant($this);
}
weight 3.10;

rule G_INVARIANT : N_ELEMENTS, CATENARY_G_INVARIANT {
   $this->G_INVARIANT = g_invariant_from_catenary($this->N_ELEMENTS, $this->CATENARY_G_INVARIANT);
}
weight 3.10;

## USER METHODS ###

# The following four are mainly for backwards compat and convenience, since cocircuits and coloops
# are common terms. However, we want to avoid duplicating data.

# @category Axiom systems
user_method COCIRCUITS {
   return shift->DUAL->CIRCUITS;
}

# @category Enumerative properties
user_method N_COCIRCUITS {
   return shift->DUAL->N_CIRCUITS;
}

# @category Axiom systems
user_method COLOOPS {
   return shift->DUAL->LOOPS;
}

# @category Enumerative properties
user_method N_COLOOPS {
   return shift->DUAL->N_LOOPS;
}


# @category Axiom systems
# calculate the rank of a set with respect to a given matroid
# @return Int
user_method rank(Set) {
   my ($self,$S) = @_;
   my $rk=0;
   for (my $i=0; $i<$self->N_BASES; ++$i) {
       my $inter = $S * $self->BASES->[$i];
       $rk = max($rk, $inter->size());
   }
   return $rk;
}

# @category Advanced properties
# @param Matroid M
# @return Bool Whether this matroid is isomorphic to M
user_method is_isomorphic_to(Matroid) {
   my ($self, $M) = @_;
   return defined(find_row_col_permutation(new IncidenceMatrix($self->BASES), new IncidenceMatrix($M->BASES)));
}

}

# @category Producing a matroid from matroids
# Produces the __dual__ of a given matroid //m//.
# Not quite the same as calling m->[[DUAL]]. The latter returns a subobject and properties of this subobject
# are only computed upon demand.
# This function returns the dual as an independent object by computing its [[BASES]].
# @param Matroid m"
# @return Matroid"
user_function dual(Matroid) {
   my $m = shift;
   return new Matroid(N_ELEMENTS=>$m->N_ELEMENTS, BASES=>$m->DUAL->BASES);
}

# @category Other
# calculate the internal activity of a base with respect to a given ordering of all bases.
# Let the given base B = B_j be the j-th one in the given ordering B_1, B_2, ...
# The internal activity of B_j is the number of "necessary" elements in B_j,
# i.e., the number of elements i in B_j such that B_j - {i} is not a subset of any B_k, k<j.
# @return Int
user_function internal_activity(Set<Int>, Array<Set<Int>>) {
   my ($B, $bases) = @_;
   my $activity = $B->size();
   foreach my $ridge (@{all_subsets_of_k($B, $B->size()-1)}) {
      foreach my $basis (@{$bases}) {
         if ($basis == $B) {
            last;
         }
         if (incl($ridge,$basis) < 0) {
            $activity--;
            last;
         }
      }
   }
   return $activity;
}


# Local Variables:
# mode: perl
# cperl-indent-level: 3
# indent-tabs-mode:nil
# End: