OKlibrary  0.2.1.6
OrthogonalTriples.cpp
Go to the documentation of this file.
00001 // Oliver Kullmann, 25.11.2006 (Swansea)
00002 /* Copyright 2006 - 2007, 2009 Oliver Kullmann
00003 This file is part of the OKlibrary. OKlibrary is free software; you can redistribute
00004 it and/or modify it under the terms of the GNU General Public License as published by
00005 the Free Software Foundation and included in this library; either version 3 of the
00006 License, or any later version. */
00007 
00021 #include <utility>
00022 #include <vector>
00023 #include <iostream>
00024 #include <limits>
00025 #include <cassert>
00026 #include <cstdlib>
00027 #include <tr1/array>
00028 
00029 #include <boost/rational.hpp>
00030 #include <boost/range/distance.hpp>
00031 #include <boost/range/size_type.hpp>
00032 #include <boost/range/const_iterator.hpp>
00033 #include <boost/range/begin.hpp>
00034 #include <boost/range/end.hpp>
00035 #include <boost/utility.hpp>
00036 #include <boost/numeric/ublas/vector.hpp>
00037 #include <boost/numeric/ublas/io.hpp>
00038 #include <boost/graph/adjacency_list.hpp>
00039 #include <boost/graph/adjacency_matrix.hpp>
00040 #include <boost/property_map/property_map.hpp>
00041 
00042 // Some linear algebra #################################
00043 
00044 template <class Vector>
00045 inline void vector_product(Vector& p, const Vector& a, const Vector& b) {
00046   assert(a.size() == 3);
00047   assert(b.size() == 3);
00048   assert(p.size() == 3);
00049   p(0) = a(1) * b(2) - a(2) * b(1);
00050   p(1) = a(2) * b(0) - a(0) * b(2);
00051   p(2) = a(0) * b(1) - a(1) * b(0);
00052 }
00053 template <class Vector>
00054 inline Vector vector_product(const Vector& a, const Vector& b) {
00055   Vector p(3);
00056   vector_product(p, a, b);
00057   return p;
00058 }
00059 
00060 template <class Vector>
00061 inline bool orthogonal(const Vector& a, const Vector& b) {
00062   typedef typename Vector::value_type value_type;
00063   return std::abs(inner_prod(a, b)) <= 0.000001; // HACK --- only for double
00064 }
00065 
00066 // typedef int integer_type;
00067 // typedef boost::rational<integer_type> rational_number_type;
00068 // typedef boost::numeric::ublas::vector<rational_number_type> vector_type;
00069 typedef boost::numeric::ublas::vector<double> vector_type;
00070 typedef std::vector<vector_type> vector2_type;
00071 
00072 // Input and output ########################################
00073 
00074 bool read(vector2_type& ev) {
00075   unsigned int m;
00076   std::cin >> m;
00077   if (not std::cin) return false;
00078   ev.reserve(m);
00079   for (unsigned int i = 0; i < m; ++i) {
00080     vector_type v(3);
00081     std::cin >> v[0] >> v[1] >> v[2];
00082     if (not std::cin) return false;
00083     ev.push_back(v);
00084   }
00085   return true;
00086 }
00087 
00088 template <class Graph>
00089 void output_graph(const Graph& g) {
00090   std::cout << num_vertices(g) << " " << num_edges(g) << "\n\n";
00091   typedef typename boost::graph_traits<Graph>::edge_iterator edge_iterator;
00092   typedef std::pair<edge_iterator, edge_iterator> edge_range;
00093   const edge_range& r(edges(g));
00094   for (edge_iterator i = r.first; i != r.second; ++i)
00095     std::cout << source(*i, g) << " " << target(*i, g) << "\n";
00096 }
00097 template <class Graph, class EliminatedEdges, class Hypergraph>
00098 void output_remaining_edges(const Graph& g, const EliminatedEdges& map, const Hypergraph& G) {
00099   typedef typename boost::graph_traits<Graph>::edge_iterator edge_iterator;
00100   typedef std::pair<edge_iterator, edge_iterator> edge_range;
00101   const edge_range& r(edges(g));
00102   {
00103     typedef typename boost::graph_traits<Graph>::edges_size_type edges_size_type;
00104     edges_size_type count = 0;
00105     for (edge_iterator i = r.first; i != r.second; ++i)
00106       count += not get(map, *i);
00107     assert(count == num_edges(g) - 3 * G.size());
00108     std::cout << count << "\n";
00109   }
00110   for (edge_iterator i = r.first; i != r.second; ++i)
00111     if (not get(map, *i))
00112       std::cout << source(*i, g) << " " << target(*i, g) << "\n";
00113 }
00114 
00115 template <class Hypergraph>
00116 void output_hypergraph(const Hypergraph& G) {
00117   std::cout << G.size() << "\n";
00118   typedef typename Hypergraph::const_iterator const_iterator;
00119   const const_iterator& end(G.end());
00120   for (const_iterator i = G.begin(); i != end; ++i) {
00121     typedef typename Hypergraph::value_type value_type;
00122     typedef typename value_type::const_iterator const_iterator2;
00123     const const_iterator2& end2(i -> end());
00124     for (const_iterator2 j = i -> begin(); j != end2; ++j)
00125       std::cout << *j << " ";
00126     std::cout << "\n";
00127   }
00128 }
00129 
00135 template <class Graph>
00136 struct SAT_translation {
00137   typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
00138   typedef typename boost::property_map<Graph, boost::vertex_index_t>::const_type property_map_index;
00139   typedef typename boost::graph_traits<Graph>::vertices_size_type vertices_size_type;
00140 
00141   const Graph& g;
00142   const property_map_index index;
00143 
00144   SAT_translation(const Graph& g) : g(g), index(get(boost::vertex_index, g)) {}
00145   
00146   int var(const vertex_descriptor v) const {
00147     const vertices_size_type n(get(index, v));
00148     assert(n < static_cast<unsigned int>(std::numeric_limits<int>::max()));
00149     return n + 1;
00150   }
00151 
00152   template <class Hypergraph, class Map>
00153   void operator()(const Hypergraph& h, const Map& m) const {
00154     std::cout << "c Kochen-Specker paradox; generated by OKlibrary (2006)\n";
00155     std::cout << "c " << num_vertices(g) << " vertices, " << h.size() << " orthogonal triples, " << num_edges(g) - 3 * h.size() << " (remaining) orthogonal pairs\n";
00156     std::cout << "p " << num_vertices(g) << " " << h.size() + num_edges(g) << "\n";
00157     {
00158       typedef typename Hypergraph::const_iterator const_iterator;
00159       const const_iterator& end(h.end());
00160       for (const_iterator i(h.begin()); i != end; ++i) {
00161         typedef typename Hypergraph::value_type value_type;
00162         const value_type& H(*i);
00163         typedef typename value_type::const_iterator const_iterator2;
00164         const const_iterator2& end2(H.end());
00165         for (const_iterator2 j(H.begin()); j != end2; ++j)
00166           std::cout << - var(*j) << " ";
00167         std::cout << "0\n";
00168         for (const_iterator2 j(H.begin()); j != end2; ++j)
00169           for (const_iterator2 k(boost::next(j)); k != end2; ++k)
00170             std::cout << var(*j) << " " << var(*k) << " 0\n";
00171       }
00172     }
00173     {
00174       typedef typename boost::graph_traits<Graph>::edge_iterator edge_iterator;
00175       typedef std::pair<edge_iterator, edge_iterator> edge_range;
00176       const edge_range& r(edges(g));
00177       for (edge_iterator i = r.first; i != r.second; ++i)
00178         if (not get(m, *i))
00179           std::cout << var(source(*i, g)) << " " << var(target(*i, g)) << " 0\n";
00180     }
00181   }
00182 };
00183 
00184 // Graph constructions #################################
00185 
00190 template <class Graph, class Range>
00191   // Graph must be undirected
00192 inline Graph orthogonality_relation(const Range& r) {
00193   typedef std::pair<unsigned int, unsigned int> edge_type;
00194   typedef std::vector<edge_type> edge_vector_type;
00195   edge_vector_type ev;
00196 
00197   typedef typename boost::range_size<Range>::type size_type;
00198   typedef typename boost::range_const_iterator<Range>::type const_iterator;
00199   const const_iterator& end(boost::const_end(r));
00200   size_type i = 0;
00201   for (const_iterator ip(boost::const_begin(r)); ip != end; ++ip, ++i) {
00202     size_type j = i+1;
00203     for (const_iterator jp(boost::next(ip)); jp != end; ++jp, ++j) {
00204       if (orthogonal(*ip, *jp)) ev.push_back(edge_type(i, j));
00205     }
00206   }
00207   Graph g(ev.begin(), ev.end(), boost::distance(r));
00208   assert(num_vertices(g) == boost::distance(r));
00209   assert(num_edges(g) == ev.size());
00210   return g;
00211 }
00212 
00217 template <class AdjacencyGraph, class Hypergraph, class EdgeContained>
00218 void extract_triangles(const AdjacencyGraph& g, Hypergraph& hg, EdgeContained& map) {
00219   typedef typename boost::graph_traits<AdjacencyGraph>::vertex_iterator vertex_iterator;
00220   typedef std::pair<vertex_iterator, vertex_iterator> vertex_range;
00221   const vertex_range& r(vertices(g));
00222   for (vertex_iterator i = r.first; i != r.second; ++i)
00223     for (vertex_iterator j = boost::next(i); j != r.second; ++j)
00224       for (vertex_iterator k = boost::next(j); k != r.second; ++k) {
00225         typedef typename boost::graph_traits<AdjacencyGraph>::edge_descriptor edge_descriptor;
00226         typedef std::pair<edge_descriptor, bool> returned_edge_type;
00227         const returned_edge_type& e1(edge(*i, *j, g));
00228         const returned_edge_type& e2(edge(*i, *k, g));
00229         const returned_edge_type& e3(edge(*j, *k, g));
00230         if (e1.second and e2.second and e3.second) {
00231           typedef typename boost::graph_traits<AdjacencyGraph>::vertex_descriptor vertex_descriptor;
00232           typedef std::tr1::array<vertex_descriptor, 3> Hyperedge;
00233           Hyperedge H = {{*i, *j, *k}};
00234           typedef typename Hypergraph::value_type value_type;
00235           hg.push_back(value_type(H.begin(), H.end()));
00236           put(map, e1.first, true); put(map, e2.first, true); put(map, e3.first, true);
00237         }
00238       }
00239 }
00240 
00245 template <class MatrixGraph, class ListGraph>
00246 MatrixGraph adapt(const ListGraph& g) {
00247   MatrixGraph mg(num_vertices(g)); // BOOST ERROR: boost::adjacency_matrix does NOT fulfill concept IteratorConstructibleGraph
00248   typedef typename boost::graph_traits<ListGraph>::edge_iterator edge_iterator;
00249   typedef std::pair<edge_iterator, edge_iterator> edge_range;
00250   const edge_range& r(edges(g));
00251   for (edge_iterator i(r.first); i != r.second; ++i)
00252     add_edge(source(*i, g), target(*i, g), mg);
00253   return mg;
00254 }
00255 
00256 
00257 // Main type defintions
00258 
00259 typedef std::vector<unsigned int> Hyperedge;
00260 typedef std::vector<Hyperedge> Hypergraph;
00261 
00262 typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS> UndirectedGraph;
00263 typedef boost::adjacency_matrix<boost::undirectedS, boost::no_property, boost::property<boost::edge_index_t, bool> > UndirectedGraphM;
00264 
00265 int main () {
00266 
00267     vector2_type list;
00268     if (not read(list))
00269       return EXIT_FAILURE;
00270 
00271     const UndirectedGraph g(orthogonality_relation<UndirectedGraph>(list));
00272     output_graph(g);
00273     std::cout << "\n";
00274 
00275     UndirectedGraphM g_m(adapt<UndirectedGraphM>(g));
00276     assert(num_vertices(g_m) == num_vertices(g));
00277     assert(num_edges(g_m) == num_edges(g));
00278     Hypergraph G;
00279 
00280     typedef boost::property_map<UndirectedGraphM, boost::edge_index_t>::type edge_property_type;
00281     edge_property_type ep(get(boost::edge_index, g_m)); // because of this operation g_m cannot be const???
00282     {
00283       typedef boost::graph_traits<UndirectedGraphM>::edge_iterator edge_iterator;
00284       typedef std::pair<edge_iterator, edge_iterator> edge_range;
00285       const edge_range& r(edges(g_m));
00286       for (edge_iterator i = r.first; i != r.second; ++i)
00287         put(ep, *i, false);
00288     }
00289     extract_triangles(g_m, G, ep);
00290     output_hypergraph(G);
00291     output_remaining_edges(g_m, ep, G);
00292     std::cout << "\n";
00293 
00294     (SAT_translation<UndirectedGraphM>(g_m))(G, ep);
00295 }