pruner.h

Go to the documentation of this file.

/* Copyright (C) 2015-2017 Martin Albrecht, Leo Ducas.
   Copyright (C) 2018 Shi Bai
 
   This file is part of fplll. fplll is free software: you
   can redistribute it and/or modify it under the terms of the GNU Lesser
   General Public License as published by the Free Software Foundation,
   either version 2.1 of the License, or (at your option) any later version.
 
   fplll 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 Lesser General Public License for more details.
 
   You should have received a copy of the GNU Lesser General Public License
   along with fplll. If not, see <http://www.gnu.org/licenses/>. */
 
#ifndef FPLLL_PRUNER_H
#define FPLLL_PRUNER_H
 
#include "fplll/defs.h"
#include "fplll/lll.h"
#include <vector>
 
FPLLL_BEGIN_NAMESPACE
 
#define PRUNER_MAX_N 2047
 
class PruningParams
{
 
public:
  double gh_factor;                  //< radius^2/Gaussian heuristic^2
  std::vector<double> coefficients;  //< pruning coefficients
  double expectation;                //< either expected success probability or number of solutions
  PrunerMetric metric;
  std::vector<double> detailed_cost;  //< Expected nodes per level
 
  PruningParams() : gh_factor(1.), expectation(1.), metric(PRUNER_METRIC_PROBABILITY_OF_SHORTEST){};
 
  static PruningParams LinearPruningParams(int block_size, int level);
};
 
template <class FT>
int run_pruner_f(ZZ_mat<mpz_t> &B, int sel_ft, int precision = 0, int prune_start = 0,
                 int prune_end = 1, double prune_pre_nodes = 1e6, double prune_min_prob = -1,
                 double gh_factor = 1.0);
 
int run_pruner(ZZ_mat<mpz_t> &B, FloatType float_type = FT_DEFAULT, int precision = 0,
               int prune_start = 0, int prune_end = 1, double prune_pre_nodes = 1e6,
               double prune_min_prob = -1, double gh_factor = 1.0);
 
template <class FT>
void prune(/*(input)output*/ PruningParams &pruning,
           /*inputs*/
           const double enumeration_radius, const double preproc_cost, const vector<double> &gso_r,
           const double target       = .9,
           const PrunerMetric metric = PRUNER_METRIC_PROBABILITY_OF_SHORTEST,
           const int flags           = PRUNER_GRADIENT);
 
template <class FT>
void prune(/*output*/ PruningParams &pruning,
           /*inputs*/
           double enumeration_radius, const double preproc_cost,
           const vector<vector<double>> &gso_rs, const double target = .9,
           const PrunerMetric metric = PRUNER_METRIC_PROBABILITY_OF_SHORTEST,
           const int flags           = PRUNER_GRADIENT);
 
template <class FT> FT svp_probability(const PruningParams &pruning);
 
template <class FT> FT svp_probability(const vector<double> &pr);
 
template <class FT> class Pruner
{
public:
  class TestPruner;
  friend class TestPruner;
 
  explicit Pruner(const int n) : metric(PRUNER_METRIC_PROBABILITY_OF_SHORTEST), flags(0), n(n)
  {
    verbosity = flags & PRUNER_VERBOSE;
    set_tabulated_consts();
    d = n / 2;
    min_pruning_coefficients.resize(d);
    btmp.resize(d);
    bftmp.resize(n);
    fill(min_pruning_coefficients.begin(), min_pruning_coefficients.end(), 0.);
  }
 
  Pruner(const FT enumeration_radius, const FT preproc_cost, const vector<double> &gso_r,
         const FT target = 0.9, const PrunerMetric metric = PRUNER_METRIC_PROBABILITY_OF_SHORTEST,
         int flags = PRUNER_GRADIENT)
      : enumeration_radius(enumeration_radius), preproc_cost(preproc_cost), target(target),
        metric(metric), flags(flags)
  {
    verbosity = flags & PRUNER_VERBOSE;
 
    n = gso_r.size();
    d = n / 2;
    if (flags & PRUNER_CVP)
    {
      symmetry_factor = 1;
    }
    min_pruning_coefficients.resize(d);
    btmp.resize(d);
    bftmp.resize(n);
    fill(min_pruning_coefficients.begin(), min_pruning_coefficients.end(), 0.);
    set_tabulated_consts();
 
    // optimize overall cost, use
    if (flags & PRUNER_SINGLE)
    {
      opt_single = true;
      if (flags & PRUNER_HALF)
      {
        throw std::invalid_argument(
            "Error: flags PRUNER_HALF and PRUNER_SINGLE are mutually exclusive.");
      }
    }
 
    // need to fix target if possible
    if (metric == PRUNER_METRIC_PROBABILITY_OF_SHORTEST)
    {
      // check if the target is reasonable; otherwise set it to 0.99.
      // note the target may affect the cost in target_function().
      if (this->target >= 1.0 || this->target <= 0.0)
      {
        throw std::invalid_argument("Invalid value for target with metric "
                                    "PRUNER_METRIC_PROBABILITY_OF_SHORTEST (need 0 < target < 1).");
      }
    }
    else if (metric == PRUNER_METRIC_EXPECTED_SOLUTIONS)
    {
      // check if the target is reasonable; otherwise set it to 0.99.
      // note the target is allowed to be larger than 1 in this case.
      // Also the target may affect the cost in target_function().
      if (this->target <= 0.0)
      {
        throw std::invalid_argument("Invalid value for target with metric "
                                    "PRUNER_METRIC_EXPECTED_SOLUTIONS (need 0 < target).");
      }
    }
    else
    {
      throw std::invalid_argument("Pruner was set to an unknown metric");
    }
    // normalization on |b_i^*|
    load_basis_shape(gso_r);
  }
 
  Pruner(const FT enumeration_radius, const FT preproc_cost, const vector<vector<double>> &gso_rs,
         const FT target = 0.9, const PrunerMetric metric = PRUNER_METRIC_PROBABILITY_OF_SHORTEST,
         int flags = PRUNER_GRADIENT)
      : enumeration_radius(enumeration_radius), preproc_cost(preproc_cost), target(target),
        metric(metric), flags(flags)
  {
    verbosity = flags & PRUNER_VERBOSE;
    n         = gso_rs[0].size();
    d         = n / 2;
    if (flags & PRUNER_CVP)
    {
      symmetry_factor = 1;
    }
    min_pruning_coefficients.resize(d);
    btmp.resize(d);
    bftmp.resize(n);
    fill(min_pruning_coefficients.begin(), min_pruning_coefficients.end(), 0.);
    set_tabulated_consts();
 
    // optimize overall cost, use
    if (flags & PRUNER_SINGLE)
    {
      opt_single = true;
      if (flags & PRUNER_HALF)
      {
        throw std::invalid_argument(
            "Error: flags PRUNER_HALF and PRUNER_SINGLE are mutually exclusive.");
      }
    }
 
    // need to fix target if possible
    if (metric == PRUNER_METRIC_PROBABILITY_OF_SHORTEST)
    {
      // check if the target is reasonable; otherwise set it to 0.99.
      // note the target may affect the cost in target_function().
      if (this->target >= 1.0 || this->target <= 0.0)
      {
        throw std::invalid_argument("Invalid value for target with metric "
                                    "PRUNER_METRIC_PROBABILITY_OF_SHORTEST (0 < target < 1).");
      }
    }
    else if (metric == PRUNER_METRIC_EXPECTED_SOLUTIONS)
    {
      // check if the target is reasonable; otherwise set it to 0.99.
      // note the target is allowed to be larger than 1 in this case.
      // Also the target may affect the cost in target_function().
      if (this->target <= 0.0)
      {
        throw std::invalid_argument(
            "Invalid value for target with metric PRUNER_METRIC_EXPECTED_SOLUTIONS (0 < target).");
      }
    }
    else
    {
      throw std::invalid_argument("Pruner was set to an unknown metric");
    }
    load_basis_shapes(gso_rs);
  }
 
  void optimize_coefficients(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_cost_vary_prob(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_cost_fixed_prob(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_evec(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_full(/*io*/ vector<double> &pr);
 
  double single_enum_cost(/*i*/ const vector<double> &pr, vector<double> *detailed_cost = nullptr)
  {
    // cout << "# pr is --> " << pr[0] << endl;
    evec b(d);
    load_coefficients(b, pr);
 
    return single_enum_cost(b, detailed_cost).get_d();
  }
 
  double repeated_enum_cost(/*i*/ const vector<double> &pr)
  {
    vec b(n);
    load_coefficients(b, pr);
    return repeated_enum_cost(b).get_d();
  }
 
  double measure_metric(/*i*/ const vector<double> &pr)
  {
    vec b(n);
    load_coefficients(b, pr);
    return measure_metric(b).get_d();
  }
 
  FT gaussian_heuristic();
 
private:
  FT enumeration_radius;
  FT preproc_cost;
  FT target;
  PrunerMetric metric;
  bool shape_loaded = false;  
  int flags;
  int n;  
  int d;  
  vector<FT> min_pruning_coefficients;
  bool opt_single = false;
 
  double descent_starting_clock;
  static bool tabulated_values_imported;  
  static FT tabulated_factorial[PRUNER_MAX_N];
  static FT tabulated_ball_vol[PRUNER_MAX_N];
 
  FT epsilon  = std::pow(2., -7);  
  FT min_step = std::pow(2., -6);  
  FT min_cf_decrease =
      .995;  //< Maximal ratio of two consectuive cost_factor in the descent before stopping
  FT step_factor     = std::pow(2, .5);  
  FT shell_ratio     = .995;             
  FT symmetry_factor = .5;               
 
  // Following typedefs are given for readability
  using vec  = vector<FT>;  
  using evec = vector<FT>;  
  using poly = vector<FT>;  
 
  /* normalised by 1/det^(1/n) */
  vec r;                    
  vec ipv;                  
  FT normalization_factor;  
  FT normalized_radius;     
  int verbosity = 0;
  vec r_old;  
  FT logvol;
 
  /* temp. variables used internally */
  vector<FT> btmp;
  vector<FT> bftmp;
 
  void set_tabulated_consts();
 
  void load_basis_shape(const vector<double> &gso_r, bool reset_normalization = true);
 
  void load_basis_shapes(const vector<vector<double>> &gso_rs);
 
  void load_coefficients(/*o*/ evec &b, /*i*/ const vector<double> &pr);
 
  void print_coefficients(/*i*/ const vector<double> &pr);
 
  void print_coefficients(/*i*/ const evec &pr);
 
  void save_coefficients(/*o*/ vector<double> &pr, /*i*/ const evec &b);
 
  inline bool enforce(/*io*/ vec &b, /*opt i*/ const int j = 0);
 
  inline FT eval_poly(const int ld, /*i*/ const poly &p, const FT x);
 
  inline void integrate_poly(const int ld, /*io*/ poly &p);
 
  inline FT relative_volume(/*i*/ const int rd, const evec &b);
 
  FT single_enum_cost(/*i*/ const vec &b, vector<double> *detailed_cost = nullptr);
  FT single_enum_cost_evec(/*i*/ const evec &b, vector<double> *detailed_cost = nullptr);
  FT single_enum_cost_lower(/*i*/ const vec &b, vector<double> *detailed_cost = nullptr);
  FT single_enum_cost_upper(/*i*/ const vec &b, vector<double> *detailed_cost = nullptr);
 
  FT svp_probability(/*i*/ const vec &b);
  FT svp_probability_evec(/*i*/ const evec &b);
  FT svp_probability_lower(/*i*/ const vec &b);
  FT svp_probability_upper(/*i*/ const vec &b);
 
  FT expected_solutions(/*i*/ const vec &b);
  FT expected_solutions_evec(/*i*/ const evec &b);
  FT expected_solutions_lower(/*i*/ const vec &b);
  FT expected_solutions_upper(/*i*/ const vec &b);
 
  FT measure_metric(/*i*/ const vec &b);
 
  FT target_function(/*i*/ const evec &b);
 
  void target_function_gradient(/*i*/ const vec &b, /*o*/ vec &res);
 
  FT repeated_enum_cost(/*i*/ const evec &b);
 
  void optimize_coefficients_preparation(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_evec_core(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_full_core(/*io*/ vector<double> &pr);
 
  void greedy(evec &b);
 
  int gradient_descent_step(/*io*/ vec &b);
 
  int gradient_descent(/*io*/ vec &b);
 
  int nelder_mead_step(/*io*/ evec &b);
 
  void optimize_coefficients_local_adjust_decr_single(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_local_adjust_incr_prob(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_local_adjust_smooth(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_incr_prob(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_decr_prob(/*io*/ vector<double> &pr);
 
  void optimize_coefficients_local_adjust_prob(/*io*/ vector<double> &pr);
};
 
template <class FT>
bool Pruner<FT>::tabulated_values_imported = false;  
template <class FT> FT Pruner<FT>::tabulated_factorial[PRUNER_MAX_N];
template <class FT> FT Pruner<FT>::tabulated_ball_vol[PRUNER_MAX_N];
 
template <class FT> inline bool Pruner<FT>::enforce(/*io*/ vec &b, /*opt i*/ const int j)
{
  int dn      = b.size();
  int c       = (dn == d) ? 1 : 2;
  bool status = false;
 
  // the last coefficient being 1
  if ((b[dn - 1] < .999) & (j != dn - 1))
  {
    status    = 1;
    b[dn - 1] = 1.;
  }
 
  for (int i = 0; i < dn; ++i)
  {
    status |= (b[i] > 1.0001);
    b[i] = b[i] > 1 ? 1. : b[i];
 
    // note min_pruning_coefficients always has length d
    if (i / c < d && b[i] <= min_pruning_coefficients[i / c])
      b[i] = min_pruning_coefficients[i / c];
  }
 
  for (int i = j; i < dn - 1; ++i)
  {
    if (b[i + 1] < b[i])
    {
      status |= (b[i + 1] + .000001 < b[i]);
      b[i + 1] = b[i];
    }
  }
 
  for (int i = std::min(j - 1, dn - 2); i >= 0; --i)
  {
    if (b[i + 1] < b[i])
    {
      status |= (b[i + 1] + .000001 < b[i]);
      b[i] = b[i + 1];
    }
  }
  return status;
}
 
// put here to avoid warnings on inlines function should be in the header.
#include "pruner_simplex.h"
 
FPLLL_END_NAMESPACE
 
#endif /* FPLLL_PRUNER_H */