#  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 PURPOSE.  See the
#  GNU General Public License for more details.
#-------------------------------------------------------------------------------

object PermutationAction {

   # Strong generating set with respect to [[BASE]].
   property STRONG_GENERATORS : Array<Array<Int>>;

   # The number of [[STRONG_GENERATORS]].
   property N_STRONG_GENERATORS : Int;

   # A base for [[STRONG_GENERATORS]].
   property BASE : Array<Int>;



# Transversals along the stabilizer chain. 
   property TRANSVERSALS : Array<Array<Int>>;

# The number of group elements per transversal.
   property TRANSVERSAL_SIZES : Array<Int>;


   rule initial : {
      my $g=$this->GENERATORS;
      my $length=$g->size && $g->[0]->size;
      foreach (@$g) {
         if ($_->size != $length) {
            croak( "all generators must have the same length" );
         }
         my %vals;
         foreach my $entry (@$_) {
            if ($entry<0 || $entry>=$length) {
               croak( "each generator must be a permutation of (0,...,DEGREE-1)" );
            }
            $vals{$entry}++;
         }
         if (keys %vals != $length) {
            croak( "each generator must be a permutation of (0,...,DEGREE-1)" );
         }
      }
   }
   precondition : exists(GENERATORS);

   rule DEGREE : GENERATORS {
      $this->DEGREE=$this->GENERATORS->[0]->size;
   }

#counts the number of elements in each transversal that are not null
   rule TRANSVERSAL_SIZES : TRANSVERSALS {
      my $trans_sizes=new Array<Int>($this->TRANSVERSALS->size);
      for(my $i=0;$i<$this->TRANSVERSALS->size;$i++){
         my $trans_size=1;
         my $single_trans=$this->TRANSVERSALS->[$i];
         for(my $j=0;$j<$single_trans->size;$j++){
            if($single_trans->[$j]>=0){
               $trans_size++;
            }
         }
         $trans_sizes->[$i]=$trans_size;
      }
      $this->TRANSVERSAL_SIZES=$trans_sizes;
   }

   rule CHARACTER : CONJUGACY_CLASS_REPRESENTATIVES {
       $this->CHARACTER = new Array<QuadraticExtension>(map { n_fixed_points($_) } @{$this->CONJUGACY_CLASS_REPRESENTATIVES});
   }
   precondition : !exists(EXPLICIT_ORBIT_REPRESENTATIVES);

   rule CHARACTER : GENERATORS, CONJUGACY_CLASS_REPRESENTATIVES, EXPLICIT_ORBIT_REPRESENTATIVES {
       $this->CHARACTER = implicit_character($this);
   }

   rule CONJUGACY_CLASSES : GENERATORS, CONJUGACY_CLASS_REPRESENTATIVES {
      $this->CONJUGACY_CLASSES(temporary) = conjugacy_classes($this->GENERATORS, $this->CONJUGACY_CLASS_REPRESENTATIVES);
   }

   rule CONJUGACY_CLASSES, CONJUGACY_CLASS_REPRESENTATIVES : GENERATORS {
       my $pair = conjugacy_classes_and_reps($this->GENERATORS);
       $this->CONJUGACY_CLASSES = $pair->first;
       $this->CONJUGACY_CLASS_REPRESENTATIVES = $pair->second;
   }

   rule ORBITS : GENERATORS {
       $this->ORBITS = orbits_of_action($this);
   }

   rule N_ORBITS : ORBITS {
       $this->N_ORBITS = $this->ORBITS->size();
   }

   rule N_INPUT_RAYS_GENERATORS : INPUT_RAYS_GENERATORS {
       $this->N_INPUT_RAYS_GENERATORS = $this->INPUT_RAYS_GENERATORS->rows;
   }

   rule N_RAYS_GENERATORS : RAYS_GENERATORS {
       $this->N_RAYS_GENERATORS = $this->RAYS_GENERATORS->rows;
   }

   rule N_INEQUALITIES_GENERATORS : INEQUALITIES_GENERATORS {
       $this->N_INEQUALITIES_GENERATORS = $this->INEQUALITIES_GENERATORS->rows;
   }

   rule N_EQUATIONS_GENERATORS : EQUATIONS_GENERATORS {
       $this->N_EQUATIONS_GENERATORS = $this->EQUATIONS_GENERATORS->rows;
   }

   rule ORBIT_REPRESENTATIVES : GENERATORS {
       $this->ORBIT_REPRESENTATIVES = orbit_representatives($this->GENERATORS);
   }

   rule PERMUTATION_TO_ORBIT_ORDER : GENERATORS, ORBIT_REPRESENTATIVES {
       $this->PERMUTATION_TO_ORBIT_ORDER = to_orbit_order($this->GENERATORS, $this->ORBIT_REPRESENTATIVES);
   }

}

object MatrixActionOnVectors {

    # @category Symmetry
    # orbits of vectors under a matrix group
    property VECTORS_ORBITS : Array<Set<Int>>;

    rule CONJUGACY_CLASSES : GENERATORS, CONJUGACY_CLASS_REPRESENTATIVES {
        $this->CONJUGACY_CLASSES = conjugacy_classes($this->GENERATORS, $this->CONJUGACY_CLASS_REPRESENTATIVES);
    }

    rule CONJUGACY_CLASSES, CONJUGACY_CLASS_REPRESENTATIVES : GENERATORS {
        my $pair = conjugacy_classes_and_reps($this->GENERATORS);
        $this->CONJUGACY_CLASSES = $pair->first;
        $this->CONJUGACY_CLASS_REPRESENTATIVES = $pair->second;
    }

   rule CHARACTER : CONJUGACY_CLASS_REPRESENTATIVES {
       $this->CHARACTER = new Array<QuadraticExtension>(map { trace($_) } @{$this->CONJUGACY_CLASS_REPRESENTATIVES});
   }
   precondition : !exists(EXPLICIT_ORBIT_REPRESENTATIVES);
}

object Group {
    
    rule REGULAR_REPRESENTATION : PERMUTATION_ACTION.GENERATORS {
        $this->REGULAR_REPRESENTATION = regular_representation($this->PERMUTATION_ACTION);
    }
    
}

# @category Symmetry
# Given a permutation action //a// on some indices and an ordered list //domain// of sets containing these indices,
# we ask for the permutation action induced by //a// on this list of sets.
# @param PermutationAction a a permutation action that acts on some indices
# @param Array<Set<Int>> domain a list of sets of indices that //a// should act on
# @return PermutationActionOnSets
# @example Consider the symmetry group of the 3-cube acting on vertices, and induce from it the action on the facets:
# > $a = cube_group(3)->PERMUTATION_ACTION;
# > $f = new Array<Set>([[0,2,4,6],[1,3,5,7],[0,1,4,5],[2,3,6,7],[0,1,2,3],[4,5,6,7]]);
# > print $a->GENERATORS;
# | 1 0 3 2 5 4 7 6
# | 0 2 1 3 4 6 5 7
# | 0 1 4 5 2 3 6 7
# > print induced_action($a,$f)->GENERATORS;
# | 1 0 2 3 4 5
# | 2 3 0 1 4 5
# | 0 1 4 5 2 3
# @example To see what the permutation [0,2,1] induces on the array [{0,1},{0,2}], do the following:
# > $a = new Array<Set<Int>>(2);
# > $a->[0] = new Set<Int>(0,1);
# > $a->[1] = new Set<Int>(0,2);
# > $p = new PermutationAction(GENERATORS=>[[0,2,1]]);
# > print induced_action($p,$a)->properties;
# | type: PermutationActionOnSets
# |
# | GENERATORS
# | 1 0
user_function induced_action(PermutationAction, Array<Set<Int>>) {
    my ($action, $domain) = @_;
    my $iod = index_of($domain);
    my $ia = new PermutationActionOnSets;
    map {
        if (defined($action->lookup($_))) {
            $ia->$_ = induced_permutations($action->$_, $domain, $iod);
        }
    } ('GENERATORS', 'CONJUGACY_CLASS_REPRESENTATIVES');
    return $ia;
}

function nonhomog_container_orbit($$) {
    my ($gens, $elem) = @_;
    return orbit<on_nonhomog_container>($gens, $elem);
}

function homog_container_orbit($$) {
    my ($gens, $elem) = @_;
    return orbit<on_container>($gens, $elem);
}

                                             
# @category Symmetry
# Given a matrix action //a// of a group //G// on some //n//-dimensional vector space and a total degree //d//,
# calculate the //G//-invariant polynomials of total degree 0 < //deg// ≤ //d// in //n// variables.
# This is done by calculating the //a//-orbit of all monomials of total degree at most //d//.
# By a theorem of Noether, for //d// = the order of //G// the output is guaranteed to generate the entire ring of //G//-invariant polynomials.
# No effort is made to calculate a basis of the ideal generated by these invariants. 
# @param MatrixActionOnVectors a the matrix action
# @param Int d the maximal degree of the sought invariants
# @option Bool action_is_affine is the action //a// affine, ie., ignore the first row and column of the generating matrices? Default yes
# @example [application polytope] To calculate the invariants of degree at most six of the matrix action of the symmetry group of the 3-cube, type
# > $c = cube(3, group=>1);
# > print join "\n", @{group::invariant_polynomials($c->GROUP->MATRIX_ACTION, 6)};
# | x_0^2 + x_1^2 + x_2^2
# | x_0^2*x_1^2 + x_0^2*x_2^2 + x_1^2*x_2^2
# | x_0^2*x_1^2*x_2^2
# | x_0^4 + x_1^4 + x_2^4
# | x_0^4*x_1^2 + x_0^4*x_2^2 + x_0^2*x_1^4 + x_0^2*x_2^4 + x_1^4*x_2^2 + x_1^2*x_2^4
# | x_0^6 + x_1^6 + x_2^6

user_function invariant_polynomials<Scalar>(MatrixActionOnVectors<Scalar>, $, { action_is_affine => 1 }) {
    my ($a, $d, $options) = @_;
    my $n = $a->GENERATORS->[0]->rows();

    my $is_affine = $options->{action_is_affine};
    if ($is_affine) {
        $n--;
    }
    
    my $zero = new Scalar(0);
    my $one  = new Scalar(1);
    my $zero_poly = new Polynomial<Scalar>($zero, zero_vector<Int>($n));
    my $one_poly  = new Polynomial<Scalar>($one,  zero_vector<Int>($n));

    # turn the matrices in G into matrices of constant polynomials; take is_affine into account
    my @G_as_poly;
    foreach my $g(@{$a->ALL_GROUP_ELEMENTS}) {
        my $m = new Matrix<Polynomial<Scalar>>($n, $n);
        foreach my $i(0..$n-1) {
            foreach my $j(0..$n-1) {
                $m->elem($i,$j) = $is_affine
                    ? $g->elem($i+1,$j+1) * $one_poly
                    : $g->elem($i,  $j)   * $one_poly;
            }
        }
        push @G_as_poly, $m;
    }

    # make a vector of monomials (x_0, x_1, ..., x_{n-1})
    my @a_monoms;
    foreach(0..$n-1) {
        push @a_monoms, new Polynomial<Scalar>($one, unit_vector<Int>($n, $_));
    }
    my $v_monoms = new Vector<Polynomial<Scalar>>(\@a_monoms);

    # generate the exponent vectors of monomials of degree at most $d
    my $deg_ineq = new Vector(-ones_vector($n+1));
    $deg_ineq->[0] = $d;
    my $exponents = new Matrix<Int>(new polytope::Polytope(INEQUALITIES=>(unit_matrix($n+1)/$deg_ineq))->LATTICE_POINTS->minor(All, ~[0]));

    # now generate the invariants
    my $invariants = new Set<Polynomial<Scalar>>;
    foreach my $alpha(@{$exponents}) {
        # for each monomial x^alpha we implement inv_alpha = sum_{g in G} (g.x0)^alpha_0 ... (g.x_{n-1})^alpha_{n-1}
        my $inv_alpha = $zero_poly;
        foreach my $g(@G_as_poly) {
            my $img_of_monoms = $g * $v_monoms;
            my $new_monom = $one_poly;
            foreach my $i (0..$n-1) {
                if ($alpha->[$i] > 0) {
                    $new_monom *= ($img_of_monoms->[$i])^(int($alpha->[$i])); # omitting the int() miscompiles on gcc 5.3.0
                }
            }
            $inv_alpha += $new_monom;
        }
        if ($inv_alpha != $zero_poly && $inv_alpha->deg() > 0) { # don't remember constants or constant coefficients
            $inv_alpha -= $inv_alpha->constant_coefficient();
            $invariants += $inv_alpha / $inv_alpha->lc();
        }
    }
    return $invariants;
}


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