Add lattice and geometric network generation

src/graph/generation/Makefile.am

@@ -22,7 +22,10 @@ libgraph_tool_generation_la_SOURCES = \
     graph_union.cc \
     graph_union_vprop.cc \
     graph_union_eprop.cc \
-    graph_triangulation.cc
+    graph_triangulation.cc \
+    graph_lattice.cc \
+    graph_geometric.cc
+
 
 libgraph_tool_generation_la_include_HEADERS = \
     graph_generation.hh \
@@ -30,4 +33,6 @@ libgraph_tool_generation_la_include_HEADERS = \
     graph_predecessor.hh \
     graph_union.hh \
     graph_triangulation.hh \
+    graph_lattice.hh \
+    graph_geometric.hh \
     sampler.hh

src/graph/generation/graph_generation.cc

@@ -97,6 +97,9 @@ void edge_property_union(GraphInterface& ugi, GraphInterface& gi,
                          boost::any uprop, boost::any prop);
 void triangulation(GraphInterface& gi, python::object points, boost::any pos,
                    string type, bool periodic);
+void lattice(GraphInterface& gi, python::object oshape, bool periodic);
+void geometric(GraphInterface& gi, python::object opoints, double r,
+               python::object orange, bool periodic, boost::any pos);
 
 using namespace boost::python;
 
@@ -110,4 +113,6 @@ BOOST_PYTHON_MODULE(libgraph_tool_generation)
     def("vertex_property_union", &vertex_property_union);
     def("edge_property_union", &edge_property_union);
     def("triangulation", &triangulation);
+    def("lattice", &lattice);
+    def("geometric", &geometric);
 }

src/graph/generation/graph_geometric.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007-2010 Tiago de Paula Peixoto <tiago@forked.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 3
+// of the License, or (at your option) any later version.
+//
+// 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+
+#include "graph_geometric.hh"
+
+#include <boost/python.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+
+typedef graph_tool::detail::get_all_graph_views
+    ::apply<graph_tool::detail::scalar_pairs, mpl::bool_<false>,
+            mpl::bool_<true>,mpl::bool_<false>,
+            mpl::bool_<true>,mpl::bool_<true> >::type graph_views;
+
+typedef property_map_types::apply<mpl::vector<vector<double> >,
+                                  GraphInterface::vertex_index_map_t,
+                                  mpl::bool_<false> >::type prop_types;
+
+
+void geometric(GraphInterface& gi, python::object opoints, double r,
+               python::object orange, bool periodic, boost::any pos)
+{
+    python::object shape = opoints.attr("shape");
+    size_t size = python::extract<size_t>(shape[0]);
+    vector<vector<double> > points(size);
+    vector<pair<double, double> > range(python::len(orange));
+
+    size = python::extract<size_t>(shape[1]);
+    for(size_t i = 0; i < points.size(); ++i)
+    {
+        points[i].resize(size);
+        for (size_t j = 0; j < points[i].size(); ++j)
+            points[i][j] = python::extract<double>(opoints[i][j]);
+    }
+
+    for(size_t i = 0; i < range.size(); ++i)
+    {
+        range[i].first = python::extract<double>(orange[i][0]);
+        range[i].second = python::extract<double>(orange[i][1]);
+    }
+
+    run_action<graph_views>()(gi, bind<void>(get_geometric(), _1, _2,
+                                             ref(points), ref(range), r,
+                                             periodic),
+                              prop_types())(pos);
+    gi.ReIndexEdges();
+}

src/graph/generation/graph_geometric.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007-2010 Tiago de Paula Peixoto <tiago@forked.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 3
+// of the License, or (at your option) any later version.
+//
+// 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef GRAPH_GEOMETRIC_HH
+#define GRAPH_GEOMETRIC_HH
+
+#include <iostream>
+#include <tr1/unordered_map>
+#include <boost/functional/hash.hpp>
+#include "graph_util.hh"
+
+namespace graph_tool
+{
+using namespace std;
+using namespace boost;
+
+
+void get_box(const vector<double>& p, double w, vector<int>& box)
+{
+    if (box.size() != p.size())
+        box.resize(p.size());
+    for (size_t i = 0; i < p.size(); ++i)
+        box[i] = int(p[i] / w);
+}
+
+template <class Point>
+double get_dist(const Point& p1, const Point& p2)
+{
+    double r = 0;
+    for (size_t i = 0; i < p1.size(); ++i)
+        r += pow(p1[i] - p2[i], 2);
+    return sqrt(r);
+}
+
+void periodic(int& x, int size)
+{
+    if (x >= int(size))
+        x -= size;
+    if (x < 0)
+        x += size;
+}
+
+
+bool is_boundary(vector<double>& point, double r,
+                 const vector<pair<double, double> > ranges)
+{
+    bool boundary = 0;
+    for (size_t j = 0; j < ranges.size(); ++j)
+    {
+        if (point[j] - r < ranges[j].first)
+        {
+            point[j] += ranges[j].second - ranges[j].first;
+            boundary = true;
+        }
+        if (point[j] + r >= ranges[j].second)
+        {
+            point[j] -= ranges[j].second - ranges[j].first;
+            boundary = true;
+        }
+    }
+    return boundary;
+}
+
+template <class Graph>
+bool is_adjacent(typename graph_traits<Graph>::vertex_descriptor v1,
+                 typename graph_traits<Graph>::vertex_descriptor v2,
+                 Graph& g)
+{
+    typename graph_traits<Graph>::out_edge_iterator e, e_end;
+    for (tie(e, e_end) = out_edges(v1, g); e != e_end; ++e)
+        if (target(*e, g) == v2)
+            return true;
+    return false;
+}
+
+struct get_geometric
+{
+    template <class Graph, class Pos>
+    void operator()(Graph& g, Pos upos, vector<vector<double> >& points,
+                    vector<pair<double, double> >& ranges,
+                    double r, bool periodic_boundary) const
+    {
+        int N = points.size();
+        typename Pos::checked_t pos = upos.get_checked();
+        vector<int> box;
+
+        tr1::unordered_multimap<vector<int>,
+                                typename graph_traits<Graph>::vertex_descriptor,
+                                boost::hash<vector<int> > > boxes;
+
+        for (int i = 0; i < N; ++i)
+        {
+            typename graph_traits<Graph>::vertex_descriptor v = add_vertex(g);
+            pos[v].resize(points[i].size());
+            copy(points[i].begin(), points[i].end(), pos[v].begin());
+
+            get_box(points[i], r, box);
+            boxes.insert(make_pair(box, v));
+
+            if (periodic_boundary)
+            {
+                //insert the same vertex at the opposite box
+                vector<double> point = points[i];
+                if (is_boundary(point, r, ranges))
+                {
+                    get_box(point, r, box);
+                    boxes.insert(make_pair(box, v));
+                }
+            }
+        }
+
+        #pragma omp parallel for default(shared) private(i, box) \
+            schedule(dynamic)
+        for (int i = 0; i < N; ++i)
+        {
+            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
+            get_box(points[i], r, box);
+
+            for (size_t k = 0; k < pow(3, box.size()); ++k)
+            {
+                for (size_t j = 0; j < box.size(); ++j)
+                    box[j] += (int(k / pow(3, j)) % 3) - 1;
+
+                typeof(boxes.begin()) iter, end;
+                for (tie(iter, end) = boxes.equal_range(box);
+                     iter != end; ++iter)
+                {
+                    typename graph_traits<Graph>::vertex_descriptor w =
+                        iter->second;
+                    double d = get_dist(pos[v], pos[w]);
+
+                    if (periodic_boundary)
+                    {
+                        if (is_adjacent(v, w, g))
+                            continue;
+                        //get correct distance for boundary
+                        vector<double> point = pos[w];
+                        if (is_boundary(point, r, ranges))
+                            d = min(d, get_dist(pos[v], point));
+                    }
+
+                    if (w > v && d <= r)
+                    {
+                        #pragma omp critical
+                        add_edge(v, w, g);
+                    }
+                }
+
+                for (size_t j = 0; j < box.size(); ++j)
+                    box[j] -= (int(k / pow(3, j)) % 3) - 1;
+            }
+        }
+    }
+};
+
+} // namespace graph_tool
+
+#endif // GRAPH_GEOMETRIC_HH

src/graph/generation/graph_lattice.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007-2010 Tiago de Paula Peixoto <tiago@forked.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 3
+// of the License, or (at your option) any later version.
+//
+// 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+
+#include "graph_lattice.hh"
+
+#include <boost/python.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void lattice(GraphInterface& gi, python::object oshape, bool periodic)
+{
+    vector<size_t> shape(python::len(oshape));
+    for(size_t i = 0; i < shape.size(); ++i)
+        shape[i] = python::extract<size_t>(oshape[i]);
+    get_lattice()(gi.GetGraph(), shape, periodic);
+    gi.ReIndexEdges();
+}

src/graph/generation/graph_lattice.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007-2010 Tiago de Paula Peixoto <tiago@forked.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 3
+// of the License, or (at your option) any later version.
+//
+// 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef GRAPH_LATTICE_HH
+#define GRAPH_LATTICE_HH
+
+#include <iostream>
+#include "graph_util.hh"
+
+namespace graph_tool
+{
+using namespace std;
+using namespace boost;
+
+
+void get_pos(size_t i, const vector<size_t>& shape, vector<int>& pos)
+{
+    size_t offset = 1;
+    for (size_t j = 0; j < shape.size(); ++j)
+    {
+        size_t L = shape[j];
+        pos[j] = ((i / offset) % L);
+        offset *= L;
+    }
+}
+
+size_t get_idx(vector<int>& pos, const vector<size_t>& shape)
+{
+    size_t offset = 1, idx = 0;
+    for (size_t j = 0; j < shape.size(); ++j)
+    {
+        idx += pos[j] * offset;
+        offset *= shape[j];
+    }
+    return idx;
+}
+
+void periodic(int& x, size_t size)
+{
+    if (x >= int(size))
+        x -= size;
+    if (x < 0)
+        x += size;
+}
+
+struct get_lattice
+{
+    template <class Graph>
+    void operator()(Graph& g, vector<size_t>& shape,
+                    bool periodic_boundary) const
+    {
+        int N = 1;
+        for (int i = 0; i < shape.size(); ++i)
+            N *= shape[i];
+        for (int i = 0; i < N; ++i)
+            add_vertex(g);
+
+        vector<int> pos(shape.size());
+        #pragma omp parallel for default(shared) private(i) \
+            firstprivate(pos) schedule(dynamic)
+        for (int i = 0; i < N; ++i)
+        {
+            get_pos(i, shape, pos);
+            for (size_t j = 0; j < shape.size(); ++j)
+            {
+                for (int k = -1; k <= 1; k += 2)
+                {
+                    pos[j] += k;
+                    if (periodic_boundary)
+                        periodic(pos[j], shape[j]);
+                    if (pos[j] > 0 && size_t(pos[j]) < shape[j])
+                    {
+                        int m = get_idx(pos, shape);
+                        if (m > i)
+                        {
+                            #pragma omp critical
+                            add_edge(vertex(i, g), vertex(m, g), g);
+                        }
+                    }
+                    pos[j] -= k;
+                    if (periodic_boundary)
+                        periodic(pos[j], shape[j]);
+                }
+            }
+        }
+    }
+};
+
+} // namespace graph_tool
+
+#endif // GRAPH_LATTICE_HH

src/graph_tool/generation/__init__.py

@@ -34,6 +34,8 @@ Summary
    line_graph
    graph_union
    triangulation
+   lattice
+   geometric_graph
 
 Contents
 ++++++++
@@ -47,7 +49,7 @@ from .. stats import label_parallel_edges, label_self_loops
 import sys, numpy, numpy.random
 
 __all__ = ["random_graph", "random_rewire", "predecessor_tree", "line_graph",
-           "graph_union", "triangulation"]
+           "graph_union", "triangulation", "lattice", "geometric_graph"]
 
 
 def random_graph(N, deg_sampler, deg_corr=None, directed=True,
@@ -642,7 +644,7 @@ def triangulation(points, type="simple", periodic=False):
     --------
     >>> from numpy.random import seed, random
     >>> seed(42)
-    >>> points = random((500,2))*4
+    >>> points = random((500, 2)) * 4
     >>> g, pos = gt.triangulation(points)
     >>> weight = g.new_edge_property("double") # Edge weights corresponding to
     ...                                        # Euclidean distances
@@ -693,3 +695,138 @@ def triangulation(points, type="simple", periodic=False):
     libgraph_tool_generation.triangulation(g._Graph__graph, points,
                                            _prop("v", g, pos), type, periodic)
     return g, pos
+
+
+def lattice(shape, periodic=False):
+    r"""
+    Generate a N-dimensional square lattice.
+
+    Parameters
+    ----------
+    shape : list or :class:`~numpy.ndarray`
+        List of sizes in each dimension.
+    periodic : bool (optional, default: False)
+        If ``True``, periodic boundary conditions will be used.
+
+    Returns
+    -------
+    lattice_graph : :class:`~graph_tool.Graph`
+        The generated graph.
+
+    See Also
+    --------
+    triangulation: 2D or 3D triangulation
+    random_graph: random graph generation
+
+    Examples
+    --------
+    >>> g = gt.lattice([10,10])
+    >>> gt.graph_draw(g, size=(8,8), output="lattice.png")
+    <...>
+    >>> g = gt.lattice([10,20], periodic=True)
+    >>> gt.graph_draw(g, size=(8,8), output="lattice_periodic.png")
+    <...>
+    >>> g = gt.lattice([10,10,10])
+    >>> gt.graph_draw(g, size=(8,8), output="lattice_3d.png")
+    <...>
+
+    .. image:: lattice.png
+    .. image:: lattice_periodic.png
+    .. image:: lattice_3d.png
+
+    *Left:* 10x10 2D lattice. *Middle:* 10x20 2D periodic lattice (torus).
+    *Right:* 10x10x10 3D lattice.
+
+    References
+    ----------
+    .. [lattice] http://en.wikipedia.org/wiki/Square_lattice
+
+    """
+
+    g = Graph(directed=False)
+    libgraph_tool_generation.lattice(g._Graph__graph, shape, periodic)
+    return g
+
+
+def geometric_graph(points, radius, ranges=None):
+    r"""
+    Generate a geometric network form a set of N-dimensional points.
+
+    Parameters
+    ----------
+    points : list or :class:`~numpy.ndarray`
+        List of points. This must be a two-dimensional array, where the rows are
+        coordinates in a N-dimensional space.
+    radius : float
+        Pairs of points with an euclidean distance lower than this parameters
+        will be connected.
+    ranges : list or :class:`~numpy.ndarray` (optional, default: None)
+        If provided, periodic boundary conditions will be assumed, and the
+        values of this parameter it will be used as the ranges in all
+        dimensions. It must be a two-dimensional array, where each row will
+        cointain the lower and upper bound of each dimension.
+
+    Returns
+    -------
+    geometric_graph : :class:`~graph_tool.Graph`
+        The generated graph.
+    pos : :class:`~graph_tool.PropertyMap`
+        A vertex property map with the position of each vertex.
+
+    Notes
+    -----
+    A geometric graph [geometric-graph]_ is generated by connecting points
+    embedded in a N-dimensional euclidean space which are at a distance equal to
+    or smaller than a given radius.
+
+    See Also
+    --------
+    triangulation: 2D or 3D triangulation
+    random_graph: random graph generation
+    lattice : N-dimensional square lattice
+
+    Examples
+    --------
+    >>> from numpy.random import seed, random
+    >>> seed(42)
+    >>> points = random((500, 2)) * 4
+    >>> g, pos = gt.geometric_graph(points, 0.3)
+    >>> gt.graph_draw(g, pos=pos, pin=True, size=(8,8), output="geometric.png")
+    <...>
+    >>> g, pos = gt.geometric_graph(points, 0.3, [(0,4), (0,4)])
+    >>> gt.graph_draw(g, size=(8,8), output="geometric_periodic.png")
+    <...>
+
+    .. image:: geometric.png
+    .. image:: geometric_periodic.png
+
+    *Left:* Geometric network with random points. *Right:* Same network, but
+     with periodic boundary conditions.
+
+    References
+    ----------
+    .. [geometric-graph] Jesper Dall and Michael Christensen, "Random geometric
+       graphs", Phys. Rev. E 66, 016121 (2002), DOI: 10.1103/PhysRevE.66.016121
+
+    """
+
+    g = Graph(directed=False)
+    pos = g.new_vertex_property("vector<double>")
+    if type(points) != numpy.ndarray:
+        points = numpy.array(points)
+    if len(points.shape) < 2:
+        raise ValueError("points list must be a two-dimensional array!")
+    if ranges is not None:
+        periodic = True
+        if type(ranges) != numpy.ndarray:
+            ranges = numpy.array(ranges, dtype="float")
+        else:
+            ranges = array(ranges, dtype="float")
+    else:
+        periodic = False
+        ranges = ()
+
+    libgraph_tool_generation.geometric(g._Graph__graph, points, float(radius),
+                                       ranges, periodic,
+                                       _prop("v", g, pos))
+    return g, pos