/* Copyright (c) 1997-2010
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Darmstadt, Germany)
   http://www.polymake.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: perturb_matrix.cc 9831 2010-10-30 18:47:05Z gawrilow $
*/

#include <Matrix.h>
#include <Vector.h>
#include <random_generators.h>
#include <polymake.h>

namespace polymake { namespace polytope {

Matrix<Rational> perturb_matrix(const Matrix<Rational>& input_matrix, double eps, bool not_hom, perl::OptionSet options)
{
   long seed;
   if(!(options["seed"] >> seed)) seed=random_seed();
   
   Matrix<double> perturbedMatrix;
   if(not_hom){
     perturbedMatrix=Matrix<double>(input_matrix);
   } else {
     perturbedMatrix=Matrix<double>(input_matrix.minor(All,~scalar2set(0)));
   }
   int n_rows=perturbedMatrix.rows();
   int n_cols=perturbedMatrix.cols();

   UniformRNG< Vector<double> > random_source(n_cols, seed);

   Matrix<double> randomMatrix(n_rows, n_cols);
   copy(random_source.begin(), entire(rows(randomMatrix)));
   perturbedMatrix+=eps*randomMatrix;

   if(not_hom){
     return Matrix<Rational>(perturbedMatrix);
   } else {
     Vector<double> fstCol(input_matrix.col(0));
     Matrix<double> result(fstCol|perturbedMatrix);
     return Matrix<Rational>(result);
   }
}

UserFunction4perl("# @category Other"
		  "# Perturb a given //input_matrix// by adding a random matrix."
		  "# The random matrix consists of vectors that are uniformly distributed"
		  "# on the unit sphere. Optionally, the random matrix can be scaled by"
		  "# a factor //eps//."
		  "# @param Matrix input_matrix"
		  "# @param Scalar eps the factor by which the random matrix is multiplied"
		  "#   default value: 1"
		  "# @param Bool not_hom if set to 1, the first column will also be perturbed;"
		  "#   otherwise the first columns of the //input_matrix// and the perturbed one"
		  "#   coincide (useful for working with homogenized coordinates);"
		  "#   default value: 0 (homogen. coords)"
		  "# @option Int seed controls the outcome of the random number generator;"
		  "#   fixing a seed number guarantees the same outcome. "
		  "# @return Matrix",
		  &perturb_matrix, "perturb_matrix(Matrix; $=1, $=0, { seed => undef })");
} }
