signing.hpp

Go to the documentation of this file.

#pragma once
 
#include "total_unimodularity.hpp"
 
#include <set>
 
#include <boost/type_traits/is_const.hpp>
 
#include "matrix_transposed.hpp"
#include "matrix_permuted.hpp"
#include "matrix_reorder.hpp"
#include "matrix.hpp"
#include "bipartite_graph_bfs.hpp"
 
namespace tu
{
 
  template <typename MatrixType>
  bool find_nonzero_column(MatrixType& matrix, size_t column_first, size_t column_beyond, size_t row_first, size_t row_beyond, size_t target_column)
  {
    for (size_t column = column_first; column < column_beyond; ++column)
    {
      for (size_t row = row_first; row < row_beyond; ++row)
      {
        if (matrix(row, column) != 0)
        {
          matrix_permute2(matrix, target_column, column);
          return true;
        }
 
      }
    }
    return false;
  }
 
  template <typename MatrixType>
  void check_sign(const MatrixType& matrix, const std::vector <bipartite_graph_bfs_node>& spanning_tree, const bipartite_graph_dimensions& dim,
      const std::set <size_t>& nodes, size_t current_index, size_t column, std::map <size_t, bool>& changes)
  {
    if (spanning_tree[current_index].predecessor == current_index)
    {
      changes[dim.index_to_row(current_index)] = false;
      return;
    }
 
    int value = matrix(dim.index_to_row(current_index), column);
    size_t last, index = current_index;
    do
    {
      last = index;
      index = spanning_tree[index].predecessor;
      std::pair <size_t, size_t> coords = dim.indexes_to_coordinates(index, last);
      value += matrix(coords.first, coords.second);
    }
    while (nodes.find(index) == nodes.end());
 
    if (changes.find(dim.index_to_row(index)) == changes.end())
    {
      check_sign(matrix, spanning_tree, dim, nodes, index, column, changes);
    }
 
    value += matrix(dim.index_to_row(index), column);
    if (changes[dim.index_to_row(index)])
    {
      value += 2;
    }
 
    value = (value >= 0 ? value : -value) % 4;
    changes[dim.index_to_row(current_index)] = (value == 2);
 
    if (value != 0 && value != 2)
    {
      throw std::logic_error("Signing procedure: modulo-sum of cycle was neither 0, nor 2!");
    }
  }
 
  template <typename T>
  struct abs_greater
  {
    bool operator()(const T& first, const T& second)
    {
      T abs_first = first >= 0 ? first : -first;
      T abs_second = second >= 0 ? second : -second;
      return abs_first > abs_second;
    }
  };
 
  template <typename M>
  bool sign_matrix(M& matrix, submatrix_indices* violator)
  {
    bool result = true;
    matrix_permuted <M> permuted(matrix);
    size_t handled_rows = 0;
 
    for (size_t handled_columns = 0; handled_columns < permuted.size2(); ++handled_columns)
    {
      if (find_nonzero_column(permuted, handled_columns, permuted.size2(), 0, handled_rows, handled_columns))
      {
 
        std::set <size_t> start_nodes;
        std::set <size_t> end_nodes;
        std::set <size_t> all_nodes;
 
        bipartite_graph_dimensions dim(handled_rows, handled_columns);
        for (size_t row = 0; row < handled_rows; ++row)
        {
          if (permuted(row, handled_columns) != 0)
          {
            size_t index = dim.row_to_index(row);
            if (start_nodes.empty())
              start_nodes.insert(index);
            else
              end_nodes.insert(index);
            all_nodes.insert(index);
          }
        }
 
 
        std::vector <bipartite_graph_bfs_node> bfs_result;
        if (!bipartite_graph_bfs(permuted, dim, start_nodes, end_nodes, true, bfs_result))
          throw std::logic_error("Signing procedure: Did not reach all nodes via bfs!");
 
        std::map <size_t, bool> changes;
        for (typename std::set <size_t>::const_iterator iter = end_nodes.begin(); iter != end_nodes.end(); ++iter)
        {
          check_sign(permuted, bfs_result, dim, all_nodes, *iter, handled_columns, changes);
        }
 
        for (std::map <size_t, bool>::iterator iter = changes.begin(); iter != changes.end(); ++iter)
        {
          if (!iter->second)
            continue;
 
          if (boost::is_const <M>::value)
          {
            if (violator)
            {
              std::set <size_t> violator_rows, violator_columns;
 
              size_t index = iter->first;
              do
              {
                if (dim.is_row(index))
                  violator_rows.insert(permuted.perm1()(dim.index_to_row(index)));
                else
                  violator_columns.insert(permuted.perm2()(dim.index_to_column(index)));
 
                index = bfs_result[index].predecessor;
              }
              while (all_nodes.find(index) == all_nodes.end());
              violator_rows.insert(permuted.perm1()(dim.index_to_row(index)));
              violator_columns.insert(permuted.perm2()(handled_columns));
 
              violator->rows = submatrix_indices::indirect_array_type(violator_rows.size());
              violator->columns = submatrix_indices::indirect_array_type(violator_columns.size());
              size_t i = 0;
              for (std::set <size_t>::const_iterator iter = violator_rows.begin(); iter != violator_rows.end(); ++iter)
                violator->rows[i++] = *iter;
              i = 0;
              for (std::set <size_t>::const_iterator iter = violator_columns.begin(); iter != violator_columns.end(); ++iter)
                violator->columns[i++] = *iter;
            }
            return false;
          }
          else
          {
            size_t real_row = permuted.perm1()(dim.index_to_row(iter->first));
            size_t real_column = permuted.perm2()(handled_columns);
            matrix_set_value(matrix, real_row, real_column, -matrix(real_row, real_column));
 
            result = false;
          }
        }
 
        matrix_reorder_rows(permuted, handled_rows, permuted.size1(), handled_columns, permuted.size2(), abs_greater <int> ());
 
        while (handled_rows < permuted.size1())
        {
          if (permuted(handled_rows, handled_columns) == 0)
            break;
          else
            ++handled_rows;
        }
      }
      else
      {
        for (size_t column = handled_columns; column < permuted.size2(); ++column)
        {
          size_t count = 0;
          for (size_t row = handled_rows; row < permuted.size1(); ++row)
          {
            if (permuted(row, column) != 0)
              ++count;
          }
 
          if (count > 0)
          {
            matrix_reorder_rows(permuted, handled_rows, permuted.size1(), handled_columns, permuted.size2(), abs_greater <int> ());
            while (handled_rows < permuted.size1())
            {
              if (permuted(handled_rows, handled_columns) == 0)
                break;
              else
                ++handled_rows;
            }
 
            break;
          }
        }
      }
    }
 
    return result;
  }
 
} /* namespace tu */