diff --git a/.gitignore b/.gitignore
index 83b48e0..525da29 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,3 +12,4 @@ PyTorch/sparseconvnet/SCN/__init__.py
 sparseconvnet.egg-info
 *.zip
 *.rar
+experiment
diff --git a/setup.py b/setup.py
index 7b189ab..84c0106 100644
--- a/setup.py
+++ b/setup.py
@@ -15,7 +15,7 @@
 torch_dir = os.path.dirname(torch.__file__)
 conda_include_dir = '/'.join(torch_dir.split('/')[:-4]) + '/include'
 
-extra = {'cxx': ['-std=c++11', '-fopenmp'], 'nvcc': ['-std=c++11', '-Xcompiler', '-fopenmp']}
+extra = {'cxx': ['-std=c++14', '-fopenmp'], 'nvcc': ['-std=c++14', '-Xcompiler', '-fopenmp']}
 
 setup(
     name='sparseconvnet',
diff --git a/sparseconvnet/SCN/CPU/Convolution.cpp b/sparseconvnet/SCN/CPU/Convolution.cpp
index b97a197..e3700ef 100644
--- a/sparseconvnet/SCN/CPU/Convolution.cpp
+++ b/sparseconvnet/SCN/CPU/Convolution.cpp
@@ -5,6 +5,7 @@
 // LICENSE file in the root directory of this source tree.
 
 // rows x groups x planes -> groups x rows x planes
+
 template <typename T>
 at::Tensor rule_index_select(at::Tensor &src, Int nRules, const Int *rules,
                               Int groups) {
diff --git a/sparseconvnet/SCN/Metadata/32bits.h b/sparseconvnet/SCN/Metadata/32bits.h
index 8e0b26c..3ce93c4 100644
--- a/sparseconvnet/SCN/Metadata/32bits.h
+++ b/sparseconvnet/SCN/Metadata/32bits.h
@@ -4,16 +4,39 @@
 // This source code is licensed under the BSD-style license found in the
 // LICENSE file in the root directory of this source tree.
 
+#pragma once
+
 #include <array>
 
 // Using 32 bit integers for coordinates and memory calculations.
 
 using Int = int32_t;
 
+// Folly's twang_mix64 hashing function
+inline uint64_t twang_mix64(uint64_t key) noexcept {
+  key = (~key) + (key << 21); // key *= (1 << 21) - 1; key -= 1;
+  key = key ^ (key >> 24);
+  key = key + (key << 3) + (key << 8); // key *= 1 + (1 << 3) + (1 << 8)
+  key = key ^ (key >> 14);
+  key = key + (key << 2) + (key << 4); // key *= 1 + (1 << 2) + (1 << 4)
+  key = key ^ (key >> 28);
+  key = key + (key << 31); // key *= 1 + (1 << 31)
+  return key;
+}
+
 // Point<dimension> is a point in the d-dimensional integer lattice
 // (i.e. square-grid/cubic-grid, ...)
+
 template <Int dimension> using Point = std::array<Int, dimension>;
 
+template<Int dimension>
+Point<dimension> generateEmptyKey() {
+  Point<dimension> empty_key;
+  for (Int i = 0; i < dimension; ++i)
+        empty_key[i] = std::numeric_limits<Int>::min();
+  return empty_key;
+}
+
 template <Int dimension>
 Point<dimension> LongTensorToPoint(/*long*/ at::Tensor &t) {
   Point<dimension> p;
@@ -65,4 +88,20 @@ template <Int dimension> struct IntArrayHash {
   }
 };
 
+// FNV Hash function for Point<dimension>
+template <Int dimension> struct FastHash {
+  std::size_t operator()(Point<dimension> const &p) const {
+    std::size_t seed = 16777619;
+    
+    for (auto x : p) {
+
+      // from boost
+      seed ^= twang_mix64(x) + 0x9e3779b9 + (seed<<6) + (seed>>2);
+    }
+
+    return seed;
+  }
+};
+
+
 #define at_kINT at::kInt
diff --git a/sparseconvnet/SCN/Metadata/ContainerMap.h b/sparseconvnet/SCN/Metadata/ContainerMap.h
new file mode 100644
index 0000000..e593d97
--- /dev/null
+++ b/sparseconvnet/SCN/Metadata/ContainerMap.h
@@ -0,0 +1,366 @@
+#ifndef Container_Map_H
+#define Container_Map_H
+
+#include "32bits.h"
+#include "RectangularRegions.h"
+#include <algorithm>
+#include <cmath>
+#include "sparsehash/dense_hash_map"
+#include <memory>
+
+template <Int dimension> class ContainerMapIterator;
+
+template <Int dimension>
+using inner_map_t = std::vector<std::pair<Point<dimension>, int>>;
+
+template <Int dimension>
+using outer_map_t =
+    google::dense_hash_map<Point<dimension>, inner_map_t<dimension>,
+                           IntArrayHash<dimension>,
+                           std::equal_to<Point<dimension>>>;
+
+template <Int dimension> class ContainerMap {
+  using it_t = ContainerMapIterator<dimension>;
+  using RuleBook = std::vector<std::vector<Int>>;
+  using RegionIndex = Point<dimension>;
+
+public:
+  ContainerMap() {
+    const static int default_size = 3;
+    map = getNewGoogleMap<outer_map_t<dimension>>();
+    filterDimensions.fill(default_size);
+    volume = std::pow(default_size, dimension);
+  }
+
+  void setFilterDimensions(const long *szs) {
+    std::copy_n(szs, dimension, filterDimensions.data());
+    volume = 1;
+    for (long size : filterDimensions) {
+      volume *= size;
+    }
+  }
+
+  /* Uses applyOnOverlappingRegions to only operate on active points of the map.
+    Due to disregarding the rest of the points within the region, the number of
+    map accesses is expected to decrease. */
+  template <typename RulesFunc>
+  void populateActiveSites(const RectangularRegion<dimension> &region,
+                           RulesFunc rulesFunc) {
+    applyOnOverlappingRegions(region, [&](const RegionIndex &regionIndex) {
+      auto *container = getContainerByRegionIndex(regionIndex);
+      if (container) {
+        auto bounds = getRegionBounds(regionIndex);
+        RectangularRegion<dimension> gridRegion =
+            RectangularRegion<dimension>(bounds.first, bounds.second);
+        RectangularRegion<dimension> overlap =
+            getOverlappingRegion(gridRegion, region);
+
+        for (auto &point : overlap) {
+          int vectorIndex = gridRegion.offset(point);
+          if ((*container)[vectorIndex].second != -1) {
+            rulesFunc((*container)[vectorIndex]);
+          }
+        }
+      }
+    });
+  }
+
+  /* Uses applyOnOverlappingRegions to populate overlapping points.
+     Used mainly for convolutions and full convolutions. If the
+     output region is greater than 1x1x1, the number of map accesses is
+     expected to decrease. */
+  template <typename RulesFunc>
+  void populateBlock(RectangularRegion<dimension> &region,
+                     RulesFunc rulesFunc) {
+
+    applyOnOverlappingRegions(region, [&](const RegionIndex &regionIndex) {
+      auto &container = insertContainer(regionIndex);
+      auto regionBounds = getRegionBounds(regionIndex);
+      RectangularRegion<dimension> gridRegion =
+          RectangularRegion<dimension>(regionBounds.first, regionBounds.second);
+      RectangularRegion<dimension> overlap =
+          getOverlappingRegion(gridRegion, region);
+ 
+      for (auto &point : overlap) {
+        int vectorIndex = gridRegion.offset(point);
+        auto &elem = container[vectorIndex];
+
+        if (elem.second == -1) {
+          elem.first = point;
+          ++ctr;
+        }
+
+        rulesFunc(elem);
+      }
+    });
+  }
+
+  /* The following functions are to comply with the standard hash
+     map interface (begin, end, insert, find, operator[]). */
+  it_t begin() { return ContainerMapIterator<dimension>(map); }
+
+  it_t end() {
+    auto container = begin();
+    container.forwardToEnd();
+    return container;
+  }
+
+  std::pair<it_t, bool>
+  insert(const std::pair<Point<dimension>, Int> &mapElem) {
+    const Point<dimension> &p = mapElem.first;
+    auto index = getRegionIndex(p);
+    auto outerMapIt =
+        map.insert(std::make_pair(index, getEmptyInnerMap())).first;
+
+    for (auto it = outerMapIt->second.begin(); it != outerMapIt->second.end();
+         ++it) {
+      if (it->second != -1 && it->first == p) {
+        return std::make_pair(
+            ContainerMapIterator<dimension>(outerMapIt, it, map), false);
+      }
+    }
+
+    ++ctr;
+    auto regionBounds = getRegionBounds(index);
+    int vectorIndex = offset(p, regionBounds.first, regionBounds.second);
+    outerMapIt->second[vectorIndex] = mapElem;
+
+    return std::make_pair(
+        ContainerMapIterator<dimension>(
+            outerMapIt, outerMapIt->second.begin() + vectorIndex, map),
+        true);
+  }
+
+  it_t find(const Point<dimension> &p) {
+    auto index = getRegionIndex(p);
+    auto it = map.find(index);
+
+    if (it == map.end()) {
+      return this->end();
+    }
+
+    auto regionBounds = getRegionBounds(index);
+    int vectorIndex = offset(p, regionBounds.first, regionBounds.second);
+    if (it->second[vectorIndex].second != -1) {
+      return ContainerMapIterator<dimension>(
+          it, it->second.begin() + vectorIndex, map);
+    }
+
+    return end();
+  }
+
+  int count(const Point<dimension> &p) const {
+    auto index = getRegionIndex(p);
+    auto it = map.find(index);
+
+    if (it == map.end()) {
+      return 0;
+    }
+
+    auto regionBounds = getRegionBounds(index);
+    int vectorIndex = offset(p, regionBounds.first, regionBounds.second);
+
+    if (it->second[vectorIndex].second != -1) {
+      return 1;
+    }
+
+    return 0;
+  }
+
+  Int &operator[](const Point<dimension> point) {
+    return insert(make_pair(point, 0)).first->second;
+  }
+
+  size_t size() const { return ctr; }
+
+private:
+  /* Helper to easily initialise google style dense hash maps. */
+  template <typename T> static T getNewGoogleMap() {
+    T map;
+    map.set_empty_key(generateEmptyKey<dimension>());
+    return map;
+  }
+
+  const inner_map_t<dimension> *
+  getContainerByRegionIndex(const RegionIndex &index) const {
+    auto it = map.find(index);
+    if (it == map.end()) {
+      return nullptr;
+    }
+
+    return &it->second;
+  }
+
+  RegionIndex getRegionIndex(const Point<dimension> &p) const {
+    Point<dimension> res;
+
+    for (int i = 0; i < dimension; ++i) {
+      res[i] = p[i] / filterDimensions[i];
+    }
+
+    return res;
+  }
+
+  Int offset(const Point<dimension> &p, const Point<dimension> &lb,
+             const Point<dimension> &ub) const {
+    Int of = 0, m = 1;
+    for (Int i = dimension - 1; i >= 0; --i) {
+      of += m * (p[i] - lb[i]);
+      m *= ub[i] - lb[i] + 1;
+    }
+    return of;
+  }
+
+  inner_map_t<dimension> getEmptyInnerMap() const {
+    auto constant_value = make_pair(generateEmptyKey<dimension>(), -1);
+    auto vec = inner_map_t<dimension>(volume, constant_value);
+
+    return vec;
+  }
+
+  /* Given an arbitrary region, this function iterates over all the grid
+  regions in ContainerMap which overlap with the region. Once it finds one,
+  it applies the regionFunc on it. */
+  template <typename RegionFunc>
+  void applyOnOverlappingRegions(const RectangularRegion<dimension> &region,
+                                 RegionFunc regionFunc) {
+    RegionIndex lowerRegionIndex = getRegionIndex(region.lb);
+    RegionIndex upperRegionIndex = getRegionIndex(region.ub);
+    std::vector<Point<dimension>> queries;
+
+    int size = 1;
+    queries.push_back(lowerRegionIndex);
+    regionFunc(lowerRegionIndex);
+
+    for (int i = 0; i < dimension; ++i) {
+      if (lowerRegionIndex[i] != upperRegionIndex[i]) {
+        for (int query = 0; query < size; ++query) {
+          Point<dimension> next_query = queries[query]; // copy the array
+          ++next_query[i];
+          regionFunc(next_query);
+          queries.push_back(next_query);
+        }
+        size *= 2;
+      }
+    }
+  }
+
+  inner_map_t<dimension> &insertContainer(const RegionIndex &index) {
+    return map.insert(std::make_pair(index, getEmptyInnerMap())).first->second;
+  }
+
+  std::pair<Point<dimension>, Point<dimension>>
+  getRegionBounds(const RegionIndex &index) const {
+    std::pair<Point<dimension>, Point<dimension>> bounds;
+    Point<dimension> &lower = bounds.first;
+    Point<dimension> &upper = bounds.second;
+
+    for (int i = 0; i < dimension; ++i) {
+      lower[i] = index[i] * filterDimensions[i];
+      upper[i] = lower[i] + filterDimensions[i] - 1;
+    }
+
+    return bounds;
+  }
+
+  RectangularRegion<dimension>
+  getOverlappingRegion(const RectangularRegion<dimension> &region1,
+                       const RectangularRegion<dimension> &region2) const {
+    Point<dimension> lower;
+    Point<dimension> upper;
+
+    for (int i = 0; i < dimension; ++i) {
+      lower[i] = std::max(region1.lb[i], region2.lb[i]);
+      upper[i] = std::min(region1.ub[i], region2.ub[i]);
+    }
+
+    return RectangularRegion<dimension>(lower, upper);
+  }
+
+  outer_map_t<dimension> map;
+  std::array<long, dimension> filterDimensions;
+  size_t ctr = 0;
+  long volume;
+};
+
+/* Custom iterator for the ContainerMap class. Captures the state of the
+  iterations via two iterators. One respresenting the outer iterator
+  of container maps, and one representing the inner iterator of the
+  structures within the outer container map. Stores a pointer to
+  the map for the purposes of end() comparisons. */
+template <Int dimension> class ContainerMapIterator {
+public:
+  ContainerMapIterator(outer_map_t<dimension> &p) : map(&p) {
+    outer_it = p.begin();
+
+    if (outer_it != p.end()) {
+      inner_it = outer_it->second.begin();
+      while (inner_it->second == -1) {
+        ++inner_it;
+      } // there must be at least one active point in the vector
+    }
+  }
+
+  ContainerMapIterator(typename outer_map_t<dimension>::iterator o_it,
+                       typename inner_map_t<dimension>::iterator i_it,
+                       outer_map_t<dimension> &p)
+      : outer_it(o_it), inner_it(i_it), map(&p) {}
+
+  void forwardToEnd() { outer_it = map->end(); }
+
+  const ContainerMapIterator<dimension> &operator++() {
+
+    while ((++inner_it != outer_it->second.end()) && (inner_it->second == -1));
+
+    if (inner_it == outer_it->second.end()) {
+      ++outer_it;
+
+      if (outer_it != map->end()) {
+        inner_it = outer_it->second.begin();
+        while (inner_it->second == -1) {
+          ++inner_it;
+        }
+      }
+    }
+
+    return *this;
+  }
+
+  std::pair<Point<dimension>, Int> &operator*() { return *inner_it; }
+
+  const std::pair<Point<dimension>, Int> &operator*() const {
+    return *inner_it;
+  }
+
+  const typename inner_map_t<dimension>::iterator operator->() const {
+    return inner_it;
+  }
+
+  typename inner_map_t<dimension>::iterator operator->() { return inner_it; }
+
+  bool operator==(const ContainerMapIterator<dimension> &it) const {
+    bool end1 = outer_it == map->end();
+    bool end2 = it.outer_it == map->end();
+
+    if (end1 && end2) {
+      return true;
+    }
+
+    if (end1 || end2) {
+      return false;
+    }
+
+    return outer_it == it.outer_it && inner_it == it.inner_it;
+  }
+
+  bool operator!=(const ContainerMapIterator<dimension> &it) const {
+    return !(*this == it);
+  }
+
+private:
+  typename outer_map_t<dimension>::iterator outer_it;
+  typename inner_map_t<dimension>::iterator inner_it;
+  outer_map_t<dimension> *map;
+};
+
+#endif
diff --git a/sparseconvnet/SCN/Metadata/ConvolutionRules.h b/sparseconvnet/SCN/Metadata/ConvolutionRules.h
index f6ab632..935f1c6 100644
--- a/sparseconvnet/SCN/Metadata/ConvolutionRules.h
+++ b/sparseconvnet/SCN/Metadata/ConvolutionRules.h
@@ -18,6 +18,23 @@ void Convolution_InputSgToRulesAndOutputSg(SparseGrid<dimension> &inputGrid,
   for (auto const &inIter : inputGrid.mp) {
     auto outRegion = OutputRegionCalculator<dimension>(
         inIter.first, size, stride, outputSpatialSize);
+
+#if defined(DICT_CONTAINER_HASH)
+    outputGrid.mp.populateBlock(
+        outRegion, [&](std::pair<Point<dimension>, Int> &data) {
+          if (data.second == -1) {
+            data.second = outputGrid.ctr++;
+          }
+
+          auto inRegion =
+              InputRegionCalculator<dimension>(data.first, size, stride);
+          Int rulesOffset = inRegion.offset(inIter.first);
+          rules[rulesOffset].push_back(inIter.second + inputGrid.ctr);
+          rules[rulesOffset].push_back(data.second);
+        });
+#endif
+
+#if !defined(DICT_CONTAINER_HASH)
     for (auto j : outRegion) {
       auto inRegion = InputRegionCalculator<dimension>(j, size, stride);
       Int rulesOffset = inRegion.offset(inIter.first);
@@ -30,6 +47,7 @@ void Convolution_InputSgToRulesAndOutputSg(SparseGrid<dimension> &inputGrid,
       rules[rulesOffset].push_back(inIter.second + inputGrid.ctr);
       rules[rulesOffset].push_back(mapVal.first->second);
     }
+#endif
   }
 }
 
diff --git a/sparseconvnet/SCN/Metadata/FullConvolutionRules.h b/sparseconvnet/SCN/Metadata/FullConvolutionRules.h
index 6055211..50d73dc 100644
--- a/sparseconvnet/SCN/Metadata/FullConvolutionRules.h
+++ b/sparseconvnet/SCN/Metadata/FullConvolutionRules.h
@@ -19,6 +19,20 @@ void FullConvolution_InputSgToRulesAndOutputSg(
   for (auto const &inIter : inputGrid.mp) {
     auto outRegion =
         InputRegionCalculator<dimension>(inIter.first, size, stride);
+#if defined(DICT_CONTAINER_HASH)
+    outputGrid.mp.populateBlock(
+        outRegion, [&](std::pair<Point<dimension>, Int> &data) {
+          if (data.second == -1) {
+            data.second = outputGrid.ctr++;
+          }
+
+          Int rulesOffset = outRegion.offset(data.first);
+          rules[rulesOffset].push_back(inIter.second + inputGrid.ctr);
+          rules[rulesOffset].push_back(data.second);
+        });
+#endif
+
+#if !defined(DICT_CONTAINER_HASH)
     for (auto j : outRegion) {
       Int rulesOffset = outRegion.offset(j);
       auto mapVal = outputGrid.mp.insert(std::make_pair(j, 0));
@@ -26,10 +40,11 @@ void FullConvolution_InputSgToRulesAndOutputSg(
       if (mapVal.second) {
         mapVal.first->second = outputGrid.ctr++;
       }
-      
+
       rules[rulesOffset].push_back(inIter.second + inputGrid.ctr);
       rules[rulesOffset].push_back(mapVal.first->second);
     }
+#endif
   }
 }
 
diff --git a/sparseconvnet/SCN/Metadata/KdTreeAdaptor.h b/sparseconvnet/SCN/Metadata/KdTreeAdaptor.h
new file mode 100644
index 0000000..758ffc6
--- /dev/null
+++ b/sparseconvnet/SCN/Metadata/KdTreeAdaptor.h
@@ -0,0 +1,147 @@
+#ifdef DICT_KD_TREE
+
+#ifndef Kd_Tree_Adaptor_H
+#define Kd_Tree_Adaptor_H
+
+#include "32bits.h"
+#include "nanoflann.hpp"
+#include <google/dense_hash_map>
+#include <map>
+#include <torch/extension.h>
+
+template <Int dimension> struct PointContainer {
+  using VectorIndex = int;
+
+  google::dense_hash_map<Point<dimension>, VectorIndex, IntArrayHash<dimension>,
+                         std::equal_to<Point<dimension>>>
+      data;
+  std::vector<std::pair<Point<dimension>, Int>> vec;
+
+  PointContainer() { data.set_empty_key(generateEmptyKey<dimension>()); }
+
+  /* Following methods are required in order to work with nanoflann.ghpp */
+  inline size_t kdtree_get_point_count() const { return data.size(); }
+
+  inline Int kdtree_get_pt(const size_t idx, const size_t dim) const {
+    return vec[idx].first[dim];
+  }
+
+  template <class BBOX> bool kdtree_get_bbox(BBOX & /* bb */) const {
+    return false;
+  }
+};
+
+template <Int dimension> class KdTreeContainer {
+
+  using it_t = typename std::vector<std::pair<Point<dimension>, Int>>::iterator;
+
+  using index_t = nanoflann::KDTreeSingleIndexAdaptor<
+      nanoflann::L2_Simple_Adaptor<Int, PointContainer<dimension>>,
+      PointContainer<dimension>, dimension>;
+
+public:
+  KdTreeContainer() { points = PointContainer<dimension>(); }
+
+  KdTreeContainer(const KdTreeContainer &k) : points(k.points) {
+    if (k.kdTreeIndex) {
+      init();
+    }
+  }
+
+  void operator=(const KdTreeContainer &k) {
+    points = k.points;
+    if (k.kdTreeIndex) {
+      init();
+    } else {
+      kdTreeIndex.reset();
+    }
+  }
+
+  // Build the kd-tree given the points prestent in the point container
+  // strcuture.
+  void init() {
+    kdTreeIndex = std::make_unique<index_t>(
+        dimension, points, nanoflann::KDTreeSingleIndexAdaptorParams(leafSize));
+    kdTreeIndex->buildIndex();
+  }
+
+  size_t size() const { return points.data.size(); }
+
+  void setLeafSize(long *sz) {
+    if (!sz)
+      throw "Invalid filter size!";
+
+    int acc = 1;
+    long *ptr = sz;
+    for (int i = 0; i < dimension; ++i, ++ptr)
+      acc *= *ptr;
+
+    leafSize = acc;
+  }
+
+  // Given a point and a radius, finds all the points within the radius of the
+  // point. Since we are
+  // using the L2 adaptor, the radius should be squared.
+  const std::vector<std::pair<size_t, Int>>
+  search(const Point<dimension> &point, Int radius) {
+    std::vector<std::pair<size_t, Int>> ret_index;
+    nanoflann::SearchParams params;
+
+    kdTreeIndex->radiusSearch(&point[0], radius, ret_index, params);
+
+    return ret_index;
+  }
+
+  /* Following are the functions to remain compliant with the map like interface
+     for
+      sparse grid maps. */
+
+  std::pair<Point<dimension>, Int> getIndexPointData(int index) {
+    return points.vec[index];
+  }
+
+  std::pair<it_t, bool>
+  insert(const std::pair<Point<dimension>, Int> &mapElem) {
+    auto mapRes =
+        points.data.insert(make_pair(mapElem.first, points.vec.size() - 1));
+
+    if (mapRes.second) {
+      points.vec.push_back(mapElem);
+      return make_pair(std::prev(points.vec.end()), true);
+    }
+
+    // Return the iterator to the element in the vector since it is present in
+    // the map.
+    return make_pair(points.vec.begin() + mapRes.first->second, false);
+  }
+
+  int count(const Point<dimension> &point) const {
+    return points.data.count(point);
+  }
+
+  it_t begin() { return points.vec.begin(); }
+
+  it_t end() { return points.vec.end(); }
+
+  it_t find(const Point<dimension> &p) {
+    auto it = points.data.find(p);
+
+    if (it == points.data.end()) {
+      return points.vec.end();
+    }
+
+    return points.vec.begin() + it->second;
+  }
+
+  Int &operator[](Point<dimension> point) {
+    return insert(make_pair(point, 0)).first->second;
+  }
+
+private:
+  PointContainer<dimension> points;
+  std::unique_ptr<index_t> kdTreeIndex;
+  int leafSize = 100;
+};
+
+#endif
+#endif
diff --git a/sparseconvnet/SCN/Metadata/Metadata.cpp b/sparseconvnet/SCN/Metadata/Metadata.cpp
index d66623c..edb23dd 100644
--- a/sparseconvnet/SCN/Metadata/Metadata.cpp
+++ b/sparseconvnet/SCN/Metadata/Metadata.cpp
@@ -5,7 +5,6 @@
 // LICENSE file in the root directory of this source tree.
 
 #include "Metadata.h"
-
 #include "ActivePoolingRules.h"
 #include "ConvolutionRules.h"
 #include "FullConvolutionRules.h"
@@ -16,10 +15,9 @@
 
 template <Int dimension> SparseGrid<dimension>::SparseGrid() : ctr(0) {
   // Sparsehash needs a key to be set aside and never used
-  Point<dimension> empty_key;
-  for (Int i = 0; i < dimension; ++i)
-    empty_key[i] = std::numeric_limits<Int>::min();
-  mp.set_empty_key(empty_key);
+  #if defined(DICT_GOOGLE_HASH)
+    mp.set_empty_key(generateEmptyKey<dimension>());
+  #endif
 }
 
 template <typename T> T *OptionalTensorData(at::Tensor &tensor) {
@@ -128,6 +126,7 @@ void Metadata<dimension>::setInputSpatialLocations(
       v += nPlanes;
     }
   }
+
   if (locations.size(1) == dimension + 1) {
     // add new samples to batch as necessary
     auto &SGs = *inputSGs;
@@ -379,6 +378,7 @@ template <Int dimension> void Metadata<dimension>::generateRuleBooks2s2() {
   Point<dimension> p1;
   Point<2 * dimension> p2;
   Point<3 * dimension> p3;
+
   for (Int i = 0; i < dimension; ++i) {
     p1[i] = p2[i] = p3[i] = inS[i] = inputSpatialSize[i];
     p2[i + dimension] = s3[i] = 3;
@@ -436,6 +436,7 @@ RuleBook &Metadata<dimension>::getSubmanifoldRuleBook(
   auto &rb = submanifoldRuleBooks[p];
   if (rb.empty()) {
     auto &SGs = grids[LongTensorToPoint<dimension>(spatialSize)];
+
 #if defined(ENABLE_OPENMP)
     openMP ? SubmanifoldConvolution_SgsToRules_OMP(SGs, rb, size.data<long>()) :
 #endif
@@ -593,7 +594,7 @@ Metadata<dimension>::compareSparseHelper(Metadata<dimension> &mR,
     auto &Ls = L[sample];
     auto &Rs = R[sample];
     for (auto const &iter : sgL.mp) {
-      if (sgR.mp.find(iter.first) == sgR.mp.end()) {
+      if (!sgR.mp.count(iter.first)) {
         Ls.push_back(sgL.mp[iter.first] + sgL.ctr);
       } else {
         cLs.push_back(sgL.mp[iter.first] + sgL.ctr);
@@ -601,7 +602,7 @@ Metadata<dimension>::compareSparseHelper(Metadata<dimension> &mR,
       }
     }
     for (auto const &iter : sgR.mp) {
-      if (sgL.mp.find(iter.first) == sgL.mp.end()) {
+      if (!sgR.mp.count(iter.first)) {
         Rs.push_back(sgR.mp[iter.first] + sgR.ctr);
       }
     }
@@ -637,7 +638,7 @@ Metadata<dimension>::copyFeaturesHelper(Metadata<dimension> &mR,
     auto &sgR = sgsR[sample];
     auto &rs = r[sample];
     for (auto const &iter : sgL.mp) {
-      if (sgR.mp.find(iter.first) != sgR.mp.end()) {
+      if (sgR.mp.count(iter.first)) {
         rs.push_back(sgL.mp[iter.first] + sgL.ctr);
         rs.push_back(sgR.mp[iter.first] + sgR.ctr);
       }
diff --git a/sparseconvnet/SCN/Metadata/Metadata.h b/sparseconvnet/SCN/Metadata/Metadata.h
index 0570b83..60f4004 100644
--- a/sparseconvnet/SCN/Metadata/Metadata.h
+++ b/sparseconvnet/SCN/Metadata/Metadata.h
@@ -6,7 +6,8 @@
 
 #ifndef Metadata_H
 #define Metadata_H
-#include "32bits.h"
+
+#include "ContainerMap.h"
 #include <algorithm>
 #include <array>
 #include <cassert>
@@ -15,16 +16,43 @@
 #include "sparsehash/dense_hash_map"
 #include <iostream>
 #include <limits>
+#include <memory>
 #include <numeric>
 #include <random>
 #include <string>
 #include <unordered_map>
 #include <vector>
+#include "KdTreeAdaptor.h"
+
+#define DICT_GOOGLE_HASH
 
+#if defined(DICT_GOOGLE_HASH)
 template <Int dimension>
 using SparseGridMap =
     google::dense_hash_map<Point<dimension>, Int, IntArrayHash<dimension>,
                            std::equal_to<Point<dimension>>>;
+#endif
+
+#if defined(DICT_STANDARD_HASH)
+template <Int dimension>
+using SparseGridMap =
+    std::unordered_map<Point<dimension>, Int, IntArrayHash<dimension>>;
+#endif
+
+#if defined(DICT_CONTAINER_HASH)
+template <Int dimension> using SparseGridMap = ContainerMap<dimension>;
+#endif
+
+#if defined(DICT_BOOST_VECTOR)
+#include <boost/container/flat_map.hpp>
+template <Int dimension>
+using SparseGridMap = boost::container::flat_map<Point<dimension>, Int>;
+#endif
+
+#if defined(DICT_KD_TREE)
+template <Int dimension> using SparseGridMap = KdTreeContainer<dimension>;
+#endif
+
 template <Int dimension> class SparseGrid {
 public:
   Int ctr;
diff --git a/sparseconvnet/SCN/Metadata/RandomizedStrideRules.h b/sparseconvnet/SCN/Metadata/RandomizedStrideRules.h
index b24e676..335a94e 100644
--- a/sparseconvnet/SCN/Metadata/RandomizedStrideRules.h
+++ b/sparseconvnet/SCN/Metadata/RandomizedStrideRules.h
@@ -95,6 +95,10 @@ void RSR_InputSgToRulesAndOutputSg(SparseGrid<dimension> &inputGrid,
       rules[rulesOffset].push_back(outIter->second);
     }
   }
+
+  #if defined(DICT_KD_TREE)
+    outputGrid.mp.init();
+  #endif
 }
 
 template <Int dimension>
diff --git a/sparseconvnet/SCN/Metadata/RectangularRegions.h b/sparseconvnet/SCN/Metadata/RectangularRegions.h
index c5e9b3c..c97ce97 100644
--- a/sparseconvnet/SCN/Metadata/RectangularRegions.h
+++ b/sparseconvnet/SCN/Metadata/RectangularRegions.h
@@ -7,6 +7,7 @@
 #ifndef RECTANGULARREGIONS_H
 #define RECTANGULARREGIONS_H
 
+#include "32bits.h"
 
 // For iterating over the rectangular region with corners lb and ub.
 // The .end() method and operator!= are designed to allow range based for
@@ -17,8 +18,10 @@ template <Int dimension> class RectangularRegion {
 public:
   Point<dimension> lb;
   Point<dimension> ub;
+
   RectangularRegion(Point<dimension> &lb, Point<dimension> &ub)
       : lb(lb), ub(ub) {}
+
   RectangularRegionIterator<dimension> begin() {
     return RectangularRegionIterator<dimension>(*this, lb);
   }
@@ -27,8 +30,8 @@ template <Int dimension> class RectangularRegion {
     // Otherwise it would need to represent a point just outside the region
     return RectangularRegionIterator<dimension>(*this, ub);
   }
-  Int
-  offset(const Point<dimension> &p) { // Enumerate the points inside the region
+  Int offset(const Point<dimension> &p)
+      const { // Enumerate the points inside the region
     Int of = 0, m = 1;
     for (Int i = dimension - 1; i >= 0; i--) {
       of += m * (p[i] - lb[i]);
@@ -36,6 +39,16 @@ template <Int dimension> class RectangularRegion {
     }
     return of;
   }
+
+  bool contains(const Point<dimension> &p) const {
+    for (int i = 0; i < dimension; ++i) {
+      if (p[i] < lb[i] || p[i] > ub[i]) {
+        return false;
+      }
+    }
+
+    return true;
+  }
 };
 
 template <Int dimension> class RectangularRegionIterator {
@@ -69,7 +82,8 @@ template <Int dimension> class RectangularRegionIterator {
 
     return *this;
   }
-  Point<dimension> &operator*() { return point; }
+
+  const Point<dimension> &operator*() { return point; }
 };
 
 // Only to be used for checking the end point of range based for loops.
diff --git a/sparseconvnet/SCN/Metadata/SubmanifoldConvolutionRules.h b/sparseconvnet/SCN/Metadata/SubmanifoldConvolutionRules.h
index b63ca26..ebb5728 100644
--- a/sparseconvnet/SCN/Metadata/SubmanifoldConvolutionRules.h
+++ b/sparseconvnet/SCN/Metadata/SubmanifoldConvolutionRules.h
@@ -7,29 +7,82 @@
 #ifndef SUBMANIFOLDCONVOLUTIONRULES_H
 #define SUBMANIFOLDCONVOLUTIONRULES_H
 
-// Full input region for an output point
+#include <map>
+#include <math.h>
+
+// Compute the input layer lower and upper bounds and return the L2 distance
+// betweem them
 template <Int dimension>
-RectangularRegion<dimension>
-InputRegionCalculator_Submanifold(const Point<dimension> &output, long *size) {
-  Point<dimension> lb, ub;
+void computeInputBounds(const Point<dimension> &output, long *size,
+                        Point<dimension> &lb, Point<dimension> &ub) {
   for (Int i = 0; i < dimension; i++) {
     Int pad = size[i] / 2;
     lb[i] = output[i] - pad;
     ub[i] = output[i] + size[i] - 1 - pad;
   }
+}
+
+// Full input region for an output point. We also store the distance of the
+// bordering points to avoind recomputing this later.
+template <Int dimension>
+RectangularRegion<dimension>
+InputRegionCalculator_Submanifold(const Point<dimension> &output, long *size) {
+  Point<dimension> lb, ub;
+  computeInputBounds<dimension>(output, size, lb, ub);
   return RectangularRegion<dimension>(lb, ub);
 }
 
+template <Int dimension>
+std::pair<Point<dimension>, int>
+computeSearchParams(const RectangularRegion<dimension> &region) {
+  Point<dimension> middle;
+  int L2_radius = 0;
+  for (Int i = 0; i < dimension; i++) {
+    int diff = (region.ub[i] - region.lb[i]) >> 1;
+    middle[i] = region.lb[i] + diff;
+    L2_radius += (diff + 1) * (diff + 1); // accounts for rounding error as well
+  }
+  return make_pair(middle, L2_radius);
+}
+
 // Call for each convolutional / max-pooling layer, once for each batch item.
 // rules is used to carry out the "lowering" whilst carrying out the convolution
-
 template <Int dimension>
 double SubmanifoldConvolution_SgToRules(SparseGrid<dimension> &grid,
                                         RuleBook &rules, long *size) {
   double countActiveInputs = 0;
-  for (auto const &outputIter : grid.mp) {
+#if defined(DICT_KD_TREE)
+  grid.mp.init();
+#endif
+
+  for (const auto &outputIter : grid.mp) {
     auto inRegion =
         InputRegionCalculator_Submanifold<dimension>(outputIter.first, size);
+
+#if defined(DICT_KD_TREE)
+    auto searchParams = computeSearchParams(inRegion);
+    auto results = grid.mp.search(searchParams.first, searchParams.second);
+    for (const auto &inputIndex : results) {
+      auto inputPoint = grid.mp.getIndexPointData(inputIndex.first);
+      if (inRegion.contains(inputPoint.first)) {
+        int offset = inRegion.offset(inputPoint.first);
+        rules[offset].push_back(inputPoint.second + grid.ctr);
+        rules[offset].push_back(outputIter.second + grid.ctr);
+        countActiveInputs++;
+      }
+    }
+#endif
+
+#if defined(DICT_CONTAINER_HASH)
+    grid.mp.populateActiveSites(inRegion, [&](const std::pair<Point<dimension>, Int> &data) {
+      int offset = inRegion.offset(data.first);
+      rules[offset].push_back(data.second + grid.ctr);
+      rules[offset].push_back(outputIter.second + grid.ctr);
+      countActiveInputs++;
+    });
+#endif
+
+#if !defined(DICT_KD_TREE) && !defined(DICT_CONTAINER_HASH)
     Int rulesOffset = 0;
     for (auto inputPoint : inRegion) {
       auto inputIter = grid.mp.find(inputPoint);
@@ -40,6 +93,7 @@ double SubmanifoldConvolution_SgToRules(SparseGrid<dimension> &grid,
       }
       rulesOffset++;
     }
+#endif
   }
   return countActiveInputs;
 }
@@ -51,6 +105,7 @@ Int SubmanifoldConvolution_SgsToRules(SparseGrids<dimension> &SGs,
   Int countActiveInputs = 0;
   rules.clear();
   rules.resize(sd);
+
   for (Int i = 0; i < (Int)SGs.size(); i++)
     countActiveInputs +=
         SubmanifoldConvolution_SgToRules<dimension>(SGs[i], rules, size);