#  Copyright (c) 1997-2006
#  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: proj-maple 7475 2006-11-22 22:42:02Z gawrilow $

# Set up a quadratic optimization problem for Maple to find
# the nearest point in a V-polytope to a given point.

# usage:  polymake --script proj-maple FILE

application 'polytope';

if (@ARGV != 1 || ! -f $ARGV[0]) {
   die "usage: polymake --script proj FILE\n";
}

my $p=load($ARGV[0]);
my $d=$p->AMBIENT_DIM;
my @vertices=$p->dehomogenize($p->VERTICES);
my $n=@vertices;


my $sum_c="";
foreach my $i (0..$n-2) {
  $sum_c.="c$i";
  $sum_c.="+" unless $i>=$n-2;
}

print << ".";
with(Optimization):
with(LinearAlgebra):
proj := proc(p)
.

print "  local x,result";
foreach my $i (0..$n-2) { print ",c$i"; }
print ";\n";
foreach my $i (0..$n-2) {
  print "  assume(c$i,'real');\n";
}

print "  x:=";
my $i:=0;
foreach my $v (@vertices) {
  chomp $v;
  $v =~ s/\s+/,/g;
  $v = "<".$v.">";
  if ($i<=$n-2) {
    print "c$i*$v+";
  } else {
    print "(1-($sum_c))*$v;\n";
  }
  ++$i;
}

print << ".";
  result:=QPSolve((x-p).(x-p),{$sum_c<=1},assume=nonnegative);
  [result,eval(x,result[2])]
end proc;
.

1;

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