Generated on: Thu Mar 29 07:46:58 PDT 2012 for custom file set

// doxy/ or-tools/ examples/ cpp/

or-tools/examples/cpp/costas_array.cc

Go to the documentation of this file.

// Copyright 2010-2012 Google
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Costas Array Problem
//
// Finds an NxN matrix of 0s and 1s, with only one 1 per row,
// and one 1 per column, such that all displacement vectors
// between each pair of 1s are distinct.
//
// This example contains two separate implementations. CostasHard()
// uses hard constraints, whereas CostasSoft() uses a minimizer to
// minimize the number of duplicates.
#include <time.h>
#include <set>
#include <utility>
#include "base/callback.h"
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "constraint_solver/constraint_solver.h"
#include "constraint_solver/constraint_solveri.h"
#include "base/random.h"

DEFINE_int32(minsize, 0, "Minimum degree of Costas matrix.");
DEFINE_int32(maxsize, 0, "Maximum degree of Costas matrix.");
DEFINE_int32(freevar, 5, "Number of free variables.");
DEFINE_int32(freeorderedvar, 4, "Number of variables in ordered subset.");
DEFINE_int32(sublimit, 20, "Number of attempts per subtree.");
DEFINE_int32(timelimit, 120000, "Time limit for local search.");
DEFINE_bool(soft_constraints, false, "Use soft constraints.");
DEFINE_string(export_profile, "", "filename to save the profile overview");

namespace operations_research {

// Checks that all pairwise distances are unique and returns all violators
void CheckConstraintViolators(const std::vector<int64>& vars,
                              std::vector<int>* const violators) {
  int dim = vars.size();

  // Check that all indices are unique
  for (int i = 0; i < dim; ++i) {
    for (int k = i + 1; k < dim; ++k) {
      if (vars[i] == vars[k]) {
        violators->push_back(i);
        violators->push_back(k);
      }
    }
  }

  // Check that all differences are unique for each level
  for (int level = 1; level < dim; ++level) {
    for (int i = 0; i < dim - level; ++i) {
      const int difference = vars[i + level] - vars[i];

      for (int k = i + 1; k < dim - level; ++k) {
        if (difference == vars[k + level] - vars[k]) {
          violators->push_back(k + level);
          violators->push_back(k);
          violators->push_back(i + level);
          violators->push_back(i);
        }
      }
    }
  }
}

// Check that all pairwise differences are unique
bool CheckCostas(const std::vector<int64>& vars) {
  std::vector<int> violators;

  CheckConstraintViolators(vars, &violators);

  return violators.empty();
}

// Cycles all possible permutations
class OrderedLNS: public BaseLNS {
 public:
  OrderedLNS(const std::vector<IntVar*>& vars, int free_elements) :
    BaseLNS(vars.data(), vars.size()),
    free_elements_(free_elements) {
    index_ = 0;

    // Start of with the first free_elements_ as a permutations, eg. 0,1,2,3,...
    for (int i = 0; i < free_elements; ++i) {
      index_ += i * pow(static_cast<double>(vars.size()),
                        free_elements - i - 1);
    }
  }

  virtual bool NextFragment(std::vector<int>* const fragment) {
    int dim = Size();
    std::set<int> fragment_set;

    do {
      int work_index = index_;

      fragment->clear();
      fragment_set.clear();

      for (int i = 0; i < free_elements_; ++i) {
        int current_index = work_index % dim;
        work_index = work_index / dim;

        std::pair<std::set<int>::iterator, bool> ret =
            fragment_set.insert(current_index);

        // Check if element has been used before
        if (ret.second) {
          fragment->push_back(current_index);
        } else {
          break;
        }
      }

      // Go to next possible permutation
      ++index_;
    // Try again if a duplicate index is used
    } while (fragment_set.size() < free_elements_);

    return true;
  }

 private:
  int index_;
  const int free_elements_;
};

// RandomLNS is used for the local search and frees the
// number of elements specified in 'free_elements' randomly.
class RandomLNS: public BaseLNS {
 public:
  RandomLNS(const std::vector<IntVar*>& vars, int free_elements) :
    BaseLNS(vars.data(), vars.size()),
    free_elements_(free_elements),
    rand_(ACMRandom::HostnamePidTimeSeed()) {
  }

  virtual bool NextFragment(std::vector<int>* const fragment) {
    std::vector<int> weighted_elements;
    std::vector<int64> values;

    // Create weighted vector for randomizer. Add one of each possible elements
    // to the weighted vector and then add more elements depending on the
    // number of conflicts
    for (int i = 0; i < Size(); ++i) {
      values.push_back(Value(i));

      // Insert each variable at least once.
      weighted_elements.push_back(i);
    }

    CheckConstraintViolators(values, &weighted_elements);
    int size = weighted_elements.size();

    // Randomly insert elements in vector until no more options remain
    while (fragment->size() < std::min(free_elements_, size)) {
      const int index = weighted_elements[rand_.Next() % size];
      fragment->push_back(index);

      // Remove all elements with this index from weighted_elements
      for (std::vector<int>::iterator pos = weighted_elements.begin();
           pos != weighted_elements.end(); ) {
        if (*pos == index) {
          // Try to erase as many elements as possible at the same time
          std::vector<int>::iterator end = pos;

          while ((end + 1) != weighted_elements.end() && *(end + 1) == *pos) {
            ++end;
          }

          pos = weighted_elements.erase(pos, end);
        } else {
          ++pos;
        }
      }

      size = weighted_elements.size();
    }

    return true;
  }

 private:
  const int free_elements_;
  ACMRandom rand_;
};

class Evaluator {
 public:
  explicit Evaluator(const std::vector<IntVar*>& vars) : vars_(vars) {}

  // Prefer the value with the smallest domain
  int64 VarEvaluator(int64 index) {
    return vars_[index]->Size();
  }

  // Penalize for each time the value appears in a different domain,
  // as values have to be unique
  int64 ValueEvaluator(int64 id, int64 value) {
    int appearance = 0;

    for (int i = 0; i < vars_.size(); ++i) {
      if (i != id && vars_[i]->Contains(value)) {
        ++appearance;
      }
    }

    return appearance;
  }
 private:
  std::vector<IntVar*> vars_;
};

// Computes a Costas Array using soft constraints.
// Instead of enforcing that all distance vectors are unique, we
// minimize the number of duplicate distance vectors.
void CostasSoft(const int dim) {
  Solver solver("Costas");

  // For the matrix and for the count of occurrences of each possible distance
  // for each stage
  const int num_elements = dim + (2 * dim + 1) * (dim);

  // create the variables
  std::vector<IntVar*> vars;
  solver.MakeIntVarArray(num_elements, -dim, dim, "var_", &vars);

  // the costas matrix
  std::vector<IntVar*> matrix(dim);
  // number of occurrences per stage
  std::vector<IntVar*> occurences;

  // All possible values of the distance vector
  // used to count the occurrence of all these values
  std::vector<int64> possible_values(2 * dim + 1);

  for (int64 i = -dim; i <= dim; ++i) {
    possible_values[i + dim] = i;
  }

  int index = 0;

  // Initialize the matrix that contains the coordinates of all '1' per row
  for (; index < dim; ++index) {
    matrix[index] = vars[index];
    vars[index]->SetMin(1);
  }

  // First constraint for the elements in the Costas Matrix. We want
  // them to be unique.
  std::vector<IntVar*> matrix_count;
  solver.MakeIntVarArray(2 * dim + 1, 0, dim, "matrix_count_", &matrix_count);
  solver.AddConstraint(solver.MakeDistribute(matrix, possible_values,
                                             matrix_count));

  // Here we only consider the elements from 1 to dim.
  for (int64 j = dim + 1; j <= 2 * dim; ++j) {
    // Penalize if an element occurs more than once.
    vars[index]
         = solver.MakeSemiContinuousExpr(solver.MakeSum(matrix_count[j], -1),
                                         0, 1)->Var();

    occurences.push_back(vars[index++]);
  }

  // Count the number of duplicates for each stage
  for (int i = 1; i < dim; ++i) {
    std::vector<IntVar*> subset(dim - i);

    // Initialize each stage
    for (int j = 0; j < dim - i; ++j) {
      IntVar* const diff = solver.MakeDifference(vars[j + i], vars[j])->Var();
      subset[j] = diff;
    }

    // Count the number of occurrences for all possible values
    std::vector<IntVar*> domain_count;
    solver.MakeIntVarArray(2 * dim + 1, 0, dim, "domain_count_", &domain_count);
    solver.AddConstraint(solver.MakeDistribute(subset,
                                               possible_values,
                                               domain_count));

    // Penalize occurrences of more than one
    for (int64 j = 0; j <= 2 * dim; ++j) {
      vars[index] =
          solver.MakeSemiContinuousExpr(solver.MakeSum(domain_count[j], -1),
                                        0, dim - i)->Var();

      occurences.push_back(vars[index++]);
    }
  }

  // We would like to minimize the penalties that we just computed
  IntVar* const objective_var = solver.MakeSum(occurences)->Var();
  OptimizeVar* const total_duplicates = solver.MakeMinimize(objective_var, 1);

  SearchMonitor* const log = solver.MakeSearchLog(1000, objective_var);

  // Out of all solutions, we just want to store the last one.
  SolutionCollector* const collector = solver.MakeLastSolutionCollector();
  collector->Add(vars);

  // The first solution that the local optimization is based on
  Evaluator evaluator(matrix);
  DecisionBuilder* const first_solution =
      solver.MakePhase(matrix,
                       NewPermanentCallback(&evaluator,
                                            &Evaluator::VarEvaluator),
                       NewPermanentCallback(&evaluator,
                                            &Evaluator::ValueEvaluator));

  SearchLimit* const search_time_limit =
      solver.MakeLimit(FLAGS_timelimit, kint64max, kint64max, kint64max);

  // Locally optimize solutions for LNS
  SearchLimit* const fail_limit = solver.MakeLimit(kint64max, kint64max,
                                                   FLAGS_sublimit, kint64max);

  DecisionBuilder* const subdecision_builder =
      solver.MakeSolveOnce(first_solution, fail_limit);

  std::vector<LocalSearchOperator*> localSearchOperators;

  // Apply RandomLNS to free FLAGS_freevar variables at each stage
  localSearchOperators.push_back(
      solver.RevAlloc(new OrderedLNS(matrix, FLAGS_freevar)));

  // Go through all possible permutations one by one
  localSearchOperators.push_back(
      solver.RevAlloc(new OrderedLNS(matrix, FLAGS_freeorderedvar)));

  LocalSearchPhaseParameters* const ls_params =
      solver.MakeLocalSearchPhaseParameters(
          solver.ConcatenateOperators(localSearchOperators),
          subdecision_builder);

  DecisionBuilder* const second_phase =
      solver.MakeLocalSearchPhase(matrix.data(),
                                  matrix.size(),
                                  first_solution,
                                  ls_params);

  // Try to find a solution
  solver.Solve(second_phase,
               collector,
               log,
               total_duplicates,
               search_time_limit);

  if (collector->solution_count() > 0) {
    std::vector<int64> costas_matrix;
    string output;

    for (int n = 0; n < dim; ++n) {
      const int64 v = collector->Value(0, vars[n]);
      costas_matrix.push_back(v);
      StringAppendF(&output, "%3lld", v);
    }

    if (!CheckCostas(costas_matrix)) {
      LOG(INFO) << "No Costas Matrix found, closest solution displayed.";
    }

    LOG(INFO) << output;
  } else {
    LOG(INFO) << "No solution found";
  }
}

// Computes a Costas Array.
void CostasHard(const int dim) {
  SolverParameters parameters;
  if (!FLAGS_export_profile.empty()) {
    parameters.profile_level = SolverParameters::NORMAL_PROFILING;
  }
  Solver solver("costas", parameters);

  // create the variables
  std::vector<IntVar*> vars;
  solver.MakeIntVarArray(dim, -dim, dim, "var", &vars);

  // Initialize the matrix that contains the coordinates of all '1' per row
  for (int m = 0; m < dim; ++m) {
    vars[m]->SetMin(1);
  }

  solver.AddConstraint(solver.MakeAllDifferent(vars));

  // Check that the pairwise difference is unique
  for (int i = 1; i < dim; ++i) {
    std::vector<IntVar*> subset(dim - i);

    for (int j = 0; j < dim - i; ++j) {
      subset[j] = solver.MakeDifference(vars[j + i], vars[j])->Var();
    }

    solver.AddConstraint(solver.MakeAllDifferent(subset));
  }

  DecisionBuilder* const db = solver.MakePhase(vars,
                                               Solver::CHOOSE_FIRST_UNBOUND,
                                               Solver::ASSIGN_MIN_VALUE);
  solver.NewSearch(db);

  if (solver.NextSolution()) {
    std::vector<int64> costas_matrix;
    string output;

    for (int n = 0; n < dim; ++n) {
      const int64 v = vars[n]->Value();
      costas_matrix.push_back(v);
      StringAppendF(&output, "%3lld", v);
    }

    LOG(INFO) << output << " (" << solver.wall_time() << "ms)";

    CHECK(CheckCostas(costas_matrix)) <<
        ": Solution is not a valid Costas Matrix.";
  } else {
    LOG(INFO) << "No solution found.";
  }
  if (!FLAGS_export_profile.empty()) {
    solver.ExportProfilingOverview(FLAGS_export_profile);
  }
}
}  // namespace operations_research

int main(int argc, char **argv) {
  google::ParseCommandLineFlags(&argc, &argv, true);
  int min = 1;
  int max = 10;

  if (FLAGS_minsize != 0) {
    min = FLAGS_minsize;

    if (FLAGS_maxsize != 0) {
      max = FLAGS_maxsize;
    } else {
      max = min;
    }
  }

  for (int i = min; i <= max; ++i) {
    LOG(INFO) << "Computing Costas Array for dim = " << i;
    if (FLAGS_soft_constraints) {
      operations_research::CostasSoft(i);
    } else {
      operations_research::CostasHard(i);
    }
  }

  return 0;
}