#  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.
#-------------------------------------------------------------------------------

CREDIT 4ti2
  4ti2 -- A software package for algebraic, geometric and combinatorial problems on linear spaces.
  Copyright by 4ti2 team.
  http://www.4ti2.de/

# path to groebner
custom $groebner;

# path to zsolve
custom $zsolve;

CONFIGURE {
    find_program($groebner, "groebner", { prompt => "the program `groebner' from the 4ti2 package" }) or return;
    # FIXME temporarily excluded 4ti2 version 1.5 from configuration, as its interface is not compatible to version 1.3
    my $ver = `$groebner -h 2>&1 | grep -o \"version 1\..\"`;
    chomp $ver;
    if ($ver eq "version 1.5") {
       die <<'.';
polymake currently does not work with 4ti2 version 1.5 bundled with latte integrale.
Please install version 1.3 available from www.4ti2.de .
.
    }
    $zsolve = $groebner =~ s/groebner$/zsolve/r;
}

# additional command line parameter for 4ti2's "groebner"
# see the manual of 4ti2 for explanations of the various possible options
custom $_4ti2_groebner_param="";

label _4ti2

object_specialization Polytope::Lattice {

    # Writes the Groebner basis of the lattice polytope into BINOMIAL_GENERATORS.
    rule TORIC_IDEAL.GROEBNER.BINOMIAL_BASIS : LATTICE_POINTS_GENERATORS, TORIC_IDEAL.GROEBNER.ORDER_MATRIX | TORIC_IDEAL.GROEBNER.ORDER_NAME {
        my $tempname=new Tempfile;
        my $pts = $this->LATTICE_POINTS_GENERATORS->[0];
        my $npts = $pts->rows;
        my ($term_order, $given) = $this->TORIC_IDEAL->GROEBNER->give_with_name("ORDER_MATRIX | ORDER_NAME");
        groebner_term_order($tempname,$term_order,$given,$npts);
        $this->TORIC_IDEAL->GROEBNER->BINOMIAL_BASIS = groebner_run($pts, $tempname);
    }
    precondition : BOUNDED;
    weight 5.10;


    # Compute Groebner Basis using default ordering.
    # The help text claims to use degrevlex by default but the result seems to be 'ls'
    rule TORIC_IDEAL.GROEBNER(new).BINOMIAL_BASIS, TORIC_IDEAL.GROEBNER(new).ORDER_NAME : LATTICE_POINTS_GENERATORS {
        my $bb = groebner_run($this->LATTICE_POINTS_GENERATORS->[0], new Tempfile);
        $this->TORIC_IDEAL->GROEBNER->ORDER_NAME = "ls";
        $this->TORIC_IDEAL->GROEBNER->BINOMIAL_BASIS= $bb;
    }
    precondition : BOUNDED;
    weight 5.10;


}

object Cone<Rational> {

    rule _4ti2.hilbert, _4ti2.integer_points: HILBERT_BASIS_GENERATORS : CONE_AMBIENT_DIM , FACETS | INEQUALITIES {
        my $tempname=new Tempfile;

        # 4ti2 version 1.6 cannot handle the case of a full-dimensional linear space
        # hence we handle this case separately
        my $n_ineq = $this->give("FACETS | INEQUALITIES")->rows;
        if (defined (my $AH = $this->lookup("LINEAR_SPAN | EQUATIONS"))) {
            $n_ineq += $AH->rows;
        }
        if ( $n_ineq == 0 ) {
            my $cad = $this->CONE_AMBIENT_DIM;
            my $rays = new Matrix<Integer>(0,$cad);
            my $lin = unit_matrix<Integer>($cad);
            $this->HILBERT_BASIS_GENERATORS=[$rays , $lin];
            return;
        }

        zsolve_print($tempname, $this);

        _4ti2_run($zsolve,$tempname);

        zsolve_parse($tempname, $this, hilb=>1);
    }
    weight 5.35;


   # Use 4ti2 to compute the toric ideal of a cone.
   rule _4ti2.toric_ideal : CONE_TORIC_IDEAL.GROEBNER.BINOMIAL_BASIS : HILBERT_BASIS_GENERATORS, CONE_TORIC_IDEAL.GROEBNER.ORDER_MATRIX | CONE_TORIC_IDEAL.GROEBNER.ORDER_NAME {
        my $tempname=new Tempfile;
        my $pts = $this->HILBERT_BASIS_GENERATORS->[0];
        my $npts = $pts->rows();
        my ($term_order, $given) = $this->CONE_TORIC_IDEAL->GROEBNER->give_with_name("ORDER_MATRIX | ORDER_NAME");
        groebner_term_order($tempname,$term_order,$given,$npts);
        $this->CONE_TORIC_IDEAL->GROEBNER->BINOMIAL_BASIS = groebner_run($pts, $tempname);
   }

   # Compute Groebner Basis using default ordering.
   # The help text claims to use degrevlex by default but the result seems to be 'ls'
   rule _4ti2.toric_ideal : CONE_TORIC_IDEAL.GROEBNER(new).BINOMIAL_BASIS, CONE_TORIC_IDEAL.GROEBNER(new).ORDER_NAME : HILBERT_BASIS_GENERATORS {
       my $bb = groebner_run($this->HILBERT_BASIS_GENERATORS->[0], new Tempfile);
       $this->CONE_TORIC_IDEAL->GROEBNER->ORDER_NAME = "ls";
       $this->CONE_TORIC_IDEAL->GROEBNER->BINOMIAL_BASIS= $bb;
   }


}

object Polytope<Rational> {

    rule _4ti2.integer_points: LATTICE_POINTS_GENERATORS : CONE_AMBIENT_DIM , FACETS | INEQUALITIES {
        my $tempname=new Tempfile;

        zsolve_print($tempname, $this, polytope=>1);

        _4ti2_run($zsolve,$tempname);

        zsolve_parse($tempname, $this, polytope=>1);
    }
    weight 5.35;

}

sub groebner_term_order {
   my ($tempname, $term_order, $given, $lp) = @_;
   if ($given eq "ORDER_NAME") {
       if ($term_order eq 'lp') {
           # lexicographical ordering
           $term_order=unit_matrix<Int>($lp);
       } elsif ($term_order eq 'Lp') {
           # can't find the meaning of this in singular
           # and since rp exist this is propaly wrong
       } elsif ($term_order eq 'Dp') {
           # degree lexicographical ordering
           $term_order = new Matrix<Int>($lp,$lp);
           $term_order->row(0) = ones_vector<Int>($lp);
           foreach (1..$lp-1) { $term_order->elem($_, $_-1) = 1; }
       } elsif ($term_order eq 'dp') {
           # degree reverse lexicographical ordering
           $term_order = new Matrix<Int>($lp,$lp);
           $term_order->row(0) = ones_vector<Int>($lp);
           foreach (1..$lp-1) { $term_order->elem($_, $lp-$_) = -1; }
       } elsif ($term_order eq 'rp') {
           # reverse lexicographical ordering
           $term_order = new Matrix<Int>($lp,$lp);
           foreach (0..$lp-1) { $term_order->elem($_, $lp-1-$_) = 1; }
       } elsif ($term_order eq 'ls') {
           # negative lexicographical ordering
           $term_order= -unit_matrix<Int>($lp);
       } elsif ($term_order eq 'rs') {
           # negative reverse lexicographical ordering
           $term_order = new Matrix<Int>($lp,$lp);
           foreach (0..$lp-1) { $term_order->elem($_, $lp-1-$_) = -1; }
       } elsif ($term_order eq 'ds') {
           # negative degree reverse lexicographical ordering
           $term_order = new Matrix<Int>($lp,$lp);
           $term_order->row(0) = -ones_vector<Int>($lp);
           foreach (1..$lp-1) { $term_order->elem($_, $lp-$_) = 1; }
       } else {
           # has to be 'Ds' since the initial rull of Groebner
           # allows no other
           # negative degree lexicographical ordering
           $term_order = new Matrix<Int>($lp,$lp);
           $term_order->row(0) = -ones_vector<Int>($lp);
           foreach (1..$lp-1) { $term_order->elem($_, $lp-$_) = -1; }
       }
   }
   # write termorder into $tempname.cost
   groebner_write($tempname, 'cost', $term_order);
}

sub groebner_write {
    my ($filebase, $suffix, $m) = @_;

    open(my $M, ">$filebase.$suffix")
      or die "can't create temporary file $filebase.$suffix: $!";
    #header
    print $M $m->rows, " ", $m->cols, "\n";
    print $M dense($m);
    close $M;
}

sub groebner_run {
    my ($pts, $tempname)=@_;

    # write all lattice points into $tempname.mat
    groebner_write($tempname, 'mat', transpose($pts));

    # run groebner with standard precision
    my $runarb = 0;
    if ($Verbose::external) {
        dbg_print( "running 4ti2's groebner: $groebner $_4ti2_groebner_param $tempname" );
    }
    open(my $P, "$groebner $_4ti2_groebner_param $tempname 2>&1 |")
      or die "couldn't run 4ti2's 'groebner': $!";

    local $_;
    while (<$P>) {
        if (/ERROR: Try running 4ti2 with arbitrary precision/) {
            $runarb = 1;
        }
    }
    close $P;
    open($P, "<$tempname.gro")
      or $runarb = 1;
    close $P;

    # run groebner with arbitrary precision if necessary
    if ($runarb == 1) {
        if ($Verbose::external) {
            dbg_print( "running 4ti2's groebner: $groebner -parb $_4ti2_groebner_param $tempname" );
        }
        open(my $P, "$groebner -parb $_4ti2_groebner_param $tempname 2>&1 |")
          or die "couldn't run 4ti2's 'groebner': $!";

        while(<$P>) {}
        close $P;
    }

    # read groebner basis
    open(my $outfile, "<$tempname.gro")
      or die "can't open output file $tempname.gro: $!";
    my $dimline = <$outfile>;
    my @gromat=();
    while (<$outfile>) {
        my @vec = split;
        push @gromat, \@vec;
    }
    close $outfile;
    return @gromat > 0 ? new Matrix(\@gromat) : new Matrix(0,$pts->rows());
}

sub _4ti2_run {
    my $prog = shift;
    # run $prog
    if ($Verbose::external) {
        dbg_print( "running 4ti2: $prog -q @_" );
    }
    system("$prog -q @_".(!$DebugLevel && " >/dev/null 2>&1"))
        and die "couldn't run 4ti2: $prog -q @_\n";
}

sub zsolve_print {
    my ($tempname, $this, %flags) = @_;

    # prepare data: ineq / eq
    my $dim = $this->CONE_AMBIENT_DIM;
    if ($flags{polytope}) { --$dim; }

    my $F=$this->give("FACETS | INEQUALITIES");
    my $n_ineq = $F->rows;
    my $facets = dense(eliminate_denominators_in_rows($F));

    if (!check_int_limit($facets)) {
        die("input values exceed int limits");
    }

    my $n_lines = $n_ineq;
    my $ah;
    if (defined (my $AH = $this->lookup("LINEAR_SPAN | EQUATIONS"))) {
        $n_lines += $AH->rows;
        $ah=dense(eliminate_denominators_in_rows($AH));
        if (!check_int_limit($ah)) {
            die("input values exceed int limits");
        }
    }

    # matrix file
    open(my $M, ">$tempname.mat")
        or die "can't create temporary file $tempname.mat: $!";
    if ($flags{polytope}) {
        print $M "$n_lines $dim\n", $facets->minor(All,range(1,$dim));
    } else {
        print $M "$n_lines $dim\n", $facets;
    }

    # right hand side
    open(my $R, ">$tempname.rhs")
        or die "can't create temporary file $tempname.rhs: $!";
    if ($flags{polytope}) {
        print $R "1 $n_lines\n", -($facets->col(0)), " ";
    } else {
        my $rhs= zero_vector<Integer>($n_lines);
        print $R "1 $n_lines\n", $rhs, " ";
    }

    # relations file
    open(my $E, ">$tempname.rel")
        or die "can't create temporary file $tempname.rel: $!";
    print $E "1 $n_lines\n", "> " x $n_ineq;

    # equalities if necessary
    if ($n_lines > $n_ineq) {
        if ($flags{polytope}) {
            print $M $ah->minor(All,range(1,$dim));
            print $R -($ah->col(0));
        } else {
            print $M $ah;
            # right hand side has been dealt with already above.
        }
        print $E "= " x ($n_lines-$n_ineq);
    }

    close $E;
    close $M;
    close $R;
}

sub zsolve_parse {
    my ($tempname, $this, %flags) = @_;

    my (@rays, @points, @lineality);
    my $is_polytope = $flags{polytope};

    # parse output
    open(my $outfile, "<$tempname.zhom")
        or die "can't open output file $tempname.zhom: $!";
    local $_;
    my $dimline = <$outfile>;

    while (<$outfile>) {
        next if /^$/;
        # read one point and append to output
        push @rays, [ $is_polytope ? (0, split) : split ];
    }
    close $outfile;

    if ($is_polytope) {
        open(my $outfile, "<$tempname.zinhom")
            or die "can't open output file $tempname.zinhom: $!";
        $dimline = <$outfile>;
        while (<$outfile>) {
            next if /^$/;
            # read one point and append to output including proj 1
            push @points, [ 1, split ];
        }
        close $outfile;
    }

    if ($is_polytope || $flags{hilb}) {
        if (-e "$tempname.zfree") {
            open(my $outfile, "<$tempname.zfree")
                or die "can't open output file $tempname.zfree: $!";
            $dimline = <$outfile>;
            while (<$outfile>) {
                next if /^$/;
                # read one point and append to output including proj 0
                push @lineality, [ $is_polytope ? (0, split) : split ];
            }
            close $outfile;
        }
    }

    $this->LATTICE_POINTS_GENERATORS = [ \@points, \@rays , \@lineality ] if $is_polytope;
    $this->HILBERT_BASIS_GENERATORS = [ \@rays , \@lineality ] if $flags{hilb};
}


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