#  Copyright (c) 1997-2007
#  Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
#  http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de
#
#  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 PURPOSE.  See the
#  GNU General Public License for more details.
#----------------------------------------------------------------------------
#  $Project: polymake $$Id: flag_vector.rules 7568 2007-01-17 16:25:44Z gawrilow $

object Polytope;

# category: Combinatorics
# Condensed form of the flag vector, containing all entries indexed by sparse sets in `{0,...,@see DIM-1}`
# in the following order: `(1, f_0, f_1, f_2, f_{02}, f_3, f_{03}, f_{13}, f_4, f_{04}, f_{14}, f_{24}, f_{024}, f_5,...)`.
# Use Dehn-Sommerville equations, via user function @see N_FLAGS, to extend.

property FLAG_VECTOR
$type="vector<cardinal>";


# category: Combinatorics
# Coefficients of the cd-index.
# author: Axel Werner

property CD_INDEX_COEFFICIENTS
$type="vector<cardinal>";


# category: Combinatorics
# `h`-vector, defined via recursion on the face lattice of a polytope.
# Coincides for simplicial (resp. simple) polytopes with the combinatorial definition
# of the `h`-vector via shellings (resp. orderings of vertices).

property H_VECTOR
$type="vector<cardinal>";


# category: Combinatorics
# `h`-vector of the bounded subcomplex, defined for not necessarily bounded polyhedra
# which are simple (as polyhedra, i.e., @see VERTEX_DEGREES on the @see FAR_FACE do not matter).
# Coincides with the reverse `h`-vector of the dual simplicial ball.
# Note that this vector will usually start with a number of zero entries.

property BOUNDED_H_VECTOR
$type="vector<cardinal>";


# category: Combinatorics
# Cubical `h`-vector. Defined for cubical polytopes and their duals.

property CUBICAL_H_VECTOR
$type="vector<cardinal>";


# category: Combinatorics
# (Toric) `g`-vector, defined via the (generalized) `h`-vector as `g_i = h_i - h_{i-1}`.
# author: Axel Werner

property G_VECTOR
$type="vector<cardinal>";


H_VECTOR : AMBIENT_DIM, VERTICES, GRAPH
PRECONDITION: BOUNDED, SIMPLE
   $this->BOUNDED && $this->SIMPLE
BODY:
   my $obj= "0 1" . (" 0" x ($this->AMBIENT_DIM - 1));
   my $p= new Apps::polytope::RationalPolytope(VERTICES => $this->VERTICES, LINEAR_OBJECTIVE => $obj);
   client("dgraph", $p, "GRAPH", "VERTICES", "FAR_FACE", "-linear", "LINEAR_OBJECTIVE", "-generic");
   my @h;
   map { ++$h[$_] } (split /\s+/, $p->VERTEX_IN_DEGREES);
   $this->H_VECTOR="@h";

BOUNDED_H_VECTOR : H_VECTOR
PRECONDITION: BOUNDED, SIMPLE
   $this->BOUNDED && $this->SIMPLE
BODY:
   $this->BOUNDED_H_VECTOR= $this->H_VECTOR;

BOUNDED_H_VECTOR : VERTICES, N_VERTICES, FAR_FACE, TOWARDS_FAR_FACE, DIM
PRECONDITION: BOUNDED, SIMPLE_POLYHEDRON
   !$this->BOUNDED && $this->SIMPLE_POLYHEDRON
BODY:
   my $p= new Apps::polytope::RationalPolytope(VERTICES => $this->VERTICES, LINEAR_OBJECTIVE => $this->TOWARDS_FAR_FACE);
   client("dgraph", $p, "GRAPH", "VERTICES", "FAR_FACE", "-linear", "LINEAR_OBJECTIVE", "-generic");
   my @not_far_face;
   my $v=0;
   foreach my $w ($this->FAR_FACE =~ /(\d+)/g) {
      push @not_far_face, $v..($w-1) unless $v>=$w;
      $v=$w+1;
   }
   push @not_far_face, $v..($this->N_VERTICES-1) unless $v>=$this->N_VERTICES;
   my @vertex_in_degrees= split /\s+/, $p->VERTEX_IN_DEGREES;
   my @h; map { $h[$_]=0 } (0..$this->DIM);
   print STDERR @vertex_in_degrees[@not_far_face], "\n";
   map { ++$h[$_] } (@vertex_in_degrees[@not_far_face]);
   print STDERR "@h\n";
   $this->BOUNDED_H_VECTOR="@h";

H_VECTOR : F_VECTOR
PRECONDITION: SIMPLICIAL
   $this->SIMPLICIAL
WEIGHT 1.10
   client("h_vector", $this, "SIMPLICIAL");

H_VECTOR : F_VECTOR
PRECONDITION: SIMPLE
   $this->SIMPLE
WEIGHT 1.10
   client("h_vector", $this, "SIMPLE");

CUBICAL_H_VECTOR : F_VECTOR
PRECONDITION: CUBICAL
   $this->CUBICAL
WEIGHT 1.10
   client("cubical_h_vector", $this, "CUBICAL");

CUBICAL_H_VECTOR : F_VECTOR
PRECONDITION: COCUBICAL
   $this->COCUBICAL
BODY:
   client("cubical_h_vector", $this, "COCUBICAL");

G_VECTOR, H_VECTOR : DIM, CD_INDEX_COEFFICIENTS
WEIGHT 4.10
   client("toric_g_vector", $this)

FLAG_VECTOR : HASSE_DIAGRAM
WEIGHT 4.10
   client("flag_vector", $this)

F_VECTOR, F2_VECTOR : FLAG_VECTOR, DIM
   my $d=$this->DIM;
   my @Fib=fibonacci($d+1);
   my @flag_vector= $this->FLAG_VECTOR =~ /\d+/g;
   my @f_vector= map { $flag_vector[$Fib[$_+1]] } (0..$d-2);
   my ($as,$sgn)=alternating_sum_and_sign(\@f_vector); # last f-vector entry determined by Euler's equation
   push @f_vector, $sgn*(1+$sgn-$as);
   $this->F_VECTOR="@f_vector";
   my @f2_vector; # array of array references
   my $one_line_f2_vector= [ $f_vector[0] ];
   push @f2_vector, $one_line_f2_vector;
   for (my $i=1; $i<$d-1; ++$i) {
      $one_line_f2_vector=[];
      for (my $k=0; $k<$i-1; ++$k) {
	 push @{$one_line_f2_vector}, $flag_vector[$Fib[$i+1]+$Fib[$k+1]]; # line i
      }
      ($as,$sgn)=alternating_sum_and_sign($one_line_f2_vector);
      push @{$one_line_f2_vector}, $sgn*($f_vector[$i]*(1+$sgn)-$as);
      push @{$one_line_f2_vector}, $f_vector[$i];
      push @f2_vector, $one_line_f2_vector;
   }
   $one_line_f2_vector=[];
   for (my $k=0; $k<$d-1; ++$k) {
      my @col= map { $f2_vector[$_][$k] } ($k+1..$d-2);
      ($as,$sgn)=alternating_sum_and_sign(\@col);
      push @{$one_line_f2_vector}, $sgn*($f_vector[$k]*(1+$sgn)-$as);
   }
   push @$one_line_f2_vector, $f_vector[$d-1];
   push @f2_vector, $one_line_f2_vector;
   my @F2_VECTOR;
   for (my $i=0; $i<$d; ++$i) {
      my @line= map { $i>=$_ ? $f2_vector[$i][$_] : $f2_vector[$_][$i] } (0..$d-1);
      push @F2_VECTOR, "@line\n";
   }
   $this->F2_VECTOR=\@F2_VECTOR;


sub recursive_n_flags {
   my ($d, $Fib, $flag_vector, @type)=@_;
# DIM, first d+1 Fibonacci numbers, FLAG_VECTOR, type w/ unique entries in descending order
   my $k=0;
   unshift @type, $d; # only to determine the proper $k
   while ($k<$#type && $type[$k]!=$type[$k+1]+1) { ++$k }
   shift(@type);            # here we get rid of it again
   --$k;
   if ($k==$#type) {
      my $idx=0;
      $idx += $Fib->[$_+1] for @type;
      return $flag_vector->[$idx];
   } else {
      push @type, -1, $d;
      my ($low,$high)=($type[$k+2]+1,$type[$k]-2);
      pop @type; pop @type;
      my ($sum,$sgn)=(0,-1);
      foreach ($low..$high) {
	 $type[$k+1]=$_;
	 $sum += $sgn*recursive_n_flags($d,$Fib,$flag_vector,@type);
	 $sgn = -$sgn;
      }
      splice @type, $k+1, 1;
      $sum += (1-$sgn)*recursive_n_flags($d,$Fib,$flag_vector,@type);
      $sum *= -$sgn;
      return $sum;
   }
}

# category: Combinatorics
# args: type, ...
# Determine the number of flags of a given type.
# @a type must belong to `{0,...,@see DIM-1}`.
# Example: "N_FLAGS (0,3,4)" determines the entry `f_{034}` of the flag vector.

user_method N_FLAGS {
   my $self=shift;
   my @set=num_sorted_uniq(sort {$b <=> $a} @_);	# in descending order
   if (@set) {
      if ($set[0] >= $self->DIM  ||  $set[-1] < 0) {
	 die "N_FLAGS: set elements out of range (0..", $self->DIM-1, "\n";
      }
      my @Fib=fibonacci($self->DIM+1);
      my @flag_vector= split /\s+/, $self->FLAG_VECTOR;
      return recursive_n_flags($self->DIM,\@Fib,\@flag_vector,@set);
   }
   return 1;			# f_(emptyset) == 1
}

CD_INDEX_COEFFICIENTS : FLAG_VECTOR, DIM
   client("cd_index", $this)

# construct string from a cd-monomial for user method CD_INDEX
sub monomial_string {
   my ($mon, $d, $Fib) = @_;
   return "1" if $d==0;
   my ($monstr,$last_sym,$this_sym,$occurrences)=("","","",0);
   while ($d > 0) {
     if ($d == 1) {
       --$d;
       $this_sym = "c";
     } elsif ($mon >= ${$Fib}[$d-1]) {
       $mon -= ${$Fib}[$d-1];
       $d -= 2;
       $this_sym = "d";
     } else {
       --$d;
       $this_sym = "c";
     }
     if ($this_sym eq $last_sym || $occurrences==0) {
       ++$occurrences;
     } else {
       $monstr .= $occurrences==1 ? $last_sym : $last_sym . "^" . $occurrences;
       $occurrences = 1;
     }
     $last_sym=$this_sym;
   }
   $monstr .= $occurrences==1 ? $last_sym : $last_sym . "^" . $occurrences;
   return $monstr;
}

# category: Combinatorics
# Prettily print the cd-index
# Author: Axel Werner
user_method CD_INDEX() {
   my ($self)=@_;
   my @coefficients= split /\s+/, $self->CD_INDEX_COEFFICIENTS;
   my @Fib=fibonacci($self->DIM+1);
   my $cdindex="";
   my $k=0;
   foreach (@coefficients) {
      if ($_!=0) {
	 if ($k>0) {
	    $cdindex .= " + ";
	 }
	 if ($_!=1) {
	    $cdindex .= $_;
	 }
	 $cdindex .= monomial_string($k, $self->DIM, \@Fib);
      }
      ++$k;
   }
   return $cdindex;
}


# Local Variables:
# mode: perl
# c-basic-offset:3
# End:
