diff --git a/examples/life.py b/examples/life.py
index 60fdd9c..6c77008 100644
--- a/examples/life.py
+++ b/examples/life.py
@@ -2,15 +2,15 @@ import cam
 import ruleset as rs
 import neighborhood as nh
 
-def game_of_life(cell, neighbors):
+def game_of_life(coordinate, grid, neighbors):
     """
     Rules of the Game of Life.
 
     Note we ignore the second component of the neighbors tuples since
     life depends on all neighbors
     """
-    total = sum(map(lambda x: int(x[0].value), neighbors))
-    if cell.value:
+    total = sum(map(lambda x: x[1], neighbors))
+    if grid[coordinate]:
         if total < 2 or total > 3:
             return False
         else:
@@ -24,5 +24,5 @@ if __name__ == '__main__':
     c = cam.CAM(1, (100, 100))
     c.randomize()
     r = rs.Ruleset(rs.Rule.SATISFY)
-    n = nh.Neighborhood.moore(c.master.grid, True)
-    c.start_plot(50, r, n, game_of_life)
+    n = nh.Neighborhood.moore(c.master, True)
+    c.start_plot(100, r, n, game_of_life)
diff --git a/src/cam.py b/src/cam.py
index 0168afa..a56d1a3 100644
--- a/src/cam.py
+++ b/src/cam.py
@@ -1,12 +1,13 @@
 """
 
-
+@author: jrpotter
+@date: June 01, 2015
 """
 import time
 import copy
 
+import camtools
 import ruleset as rs
-import cell_plane as cp
 import neighborhood as nh
 
 import numpy as np
@@ -23,14 +24,12 @@ class CAM:
     all methods needed (i.e. supported) to interact/configure the cellular automata as desired.
     """
 
-    def __init__(self, cps=1, dimen=(100,100)):
+    def __init__(self, cps=1, dimen=(100, 100)):
         """
 
         """
-        cps = max(cps, 1)
-        self._dimen = dimen
-        self._planes = [cp.CellPlane(dimen) for i in range(cps)]
-        self.master = self._planes[0]
+        self.planes = np.zeros((max(cps, 1),) + dimen, dtype=np.int32)
+        self.master = self.planes[0]
 
 
     def start_plot(self, clock, ruleset, neighborhood, *args):
@@ -46,11 +45,11 @@ class CAM:
         ax.get_xaxis().set_visible(False)
         ax.get_yaxis().set_visible(False)
 
-        mshown = plt.matshow(self.master.to_binary(), fig.number)
+        mshown = plt.matshow(self.master, fig.number)
 
         def animate(frame):
             self.tick(ruleset, neighborhood, *args)
-            mshown.set_array(self.master.to_binary())
+            mshown.set_array(self.master)
             fig.canvas.draw()
 
         ani.FuncAnimation(fig, animate, interval=clock)
@@ -64,7 +63,7 @@ class CAM:
 
         """
         while True:
-            print(self.to_binary())
+            print(self.master)
             time.sleep(clock / 1000)
             self.tick(ruleset, neighborhood, *args)
 
@@ -77,17 +76,15 @@ class CAM:
         is placed into the master grid. Depending on the timing specifications set by the user, this
         may also change secondary cell planes (the master is always updated on each tick).
         """
-        self.master.grid[:] = rs.Ruleset.update(self.master.grid, ruleset, neighborhood, *args)
+        tmp = np.copy(self.master)
+        for i in range(len(self.master.flat)):
+            tmp.flat[i] = ruleset.call(i, self.master, neighborhood, *args)
+        self.master[:] = tmp
 
 
     def randomize(self):
         """
         Set the master grid to a random configuration.
         """
-        @np.vectorize
-        def v_random(cell):
-            cell.value = np.random.random_integers(0, 1)
-            return cell
-
-        v_random(self.master.grid)
+        self.master[:] = np.random.random_integers(0, 1, self.master.shape)
 
diff --git a/src/camtools.py b/src/camtools.py
new file mode 100644
index 0000000..939c7de
--- /dev/null
+++ b/src/camtools.py
@@ -0,0 +1,36 @@
+"""
+
+@author: jrpotter
+@date: June 01, 2015
+"""
+import numpy as np
+
+def flatten(coordinates, grid):
+    """
+    Given the coordinates of a matrix, returns the index of the flat matrix.
+    """
+    index = 0
+    for i in range(len(coordinates)):
+        index += coordinates[i] * np.prod(grid.shape[i+1:], dtype=np.int32)
+
+    return index
+
+
+def unflatten(index, grid):
+    """
+    Given an index of a flat matrix, returns the corresponding coordinates.
+    """
+    coordinates = []
+    for i in range(len(grid.shape)):
+        tmp = np.prod(grid.shape[i+1:], dtype=np.int32)
+        coordinates.append(index // tmp)
+        index -= tmp * coordinates[-1]
+
+    return tuple(coordinates)
+
+
+def comp_add(coor1, coor2):
+    """
+    Adds components of coordinates element-wise.
+    """
+    return tuple(map(sum, zip(coor1, coor2)))
diff --git a/src/cell_plane.py b/src/cell_plane.py
deleted file mode 100644
index e4eb882..0000000
--- a/src/cell_plane.py
+++ /dev/null
@@ -1,69 +0,0 @@
-"""
-
-
-
-
-"""
-import numpy as np
-import matplotlib.pyplot as plt
-
-
-class Cell:
-    """
-    Represents a "cell" in a CellPlane.
-
-    Note we keep track of the index for vectorization purposes. By maintaining each index
-    and batch updating via the given index, we can much more efficiently update the entire
-    cell plane.
-    """
-    def __init__(self, value, *index):
-        self.value = value
-        self.index = index
-
-
-class CellPlane:
-    """
-    A CellPlane represents a layer of the grids that can be placed on top of one another in a 2D CAM.
-
-    The use of multiple cell plane allow for more intricate states of life and death, though there
-    exists only a single master cell plane that controls the others. That is, the master cell plane has
-    a CAM ruleset applied to it, and the other cell planes merely copy the master, though this can
-    be delayed and have different color mappings.
-
-    For example, by setting a delay of two ticks on the second cell plane of a 2-level CAM configuration,
-    one can allow for ECHOing, providing a more intuitive sense of "velocity" based on the master.
-
-    That is not to say one could not have multiple CAM's operating simultaneously though. We can consider
-    a configuration to consist of an arbitrary number of planes, of which one is the master, but multiple
-    masters can exist in separate CAMs that can interact with one another.
-    """
-
-    @staticmethod
-    @np.vectorize
-    def _populate(*indices):
-        """
-        The following joins indices in N-dimensions together.
-
-        This information is stored in a cell (with initial value False) in order for batch processing
-        to be performed when actually updating values and computing whether a cell is on or off. For
-        example, if exploring a 4D array, we want to be able to know which cells we need to check the
-        status of, but this is relative to the current cell, whose position we do not know unless that
-        information is stored with the current cell.
-        """
-        return Cell(False, *indices)
-
-
-    def __init__(self, dimen):
-        """
-
-        """
-        self.grid = CellPlane._populate(*np.indices(dimen))
-
-
-    def to_binary(self):
-        """
-
-        """
-        vfunc = np.vectorize(lambda x: int(x.value))
-        return vfunc(self.grid)
-
diff --git a/src/neighborhood.py b/src/neighborhood.py
index b8f4bfa..356053d 100644
--- a/src/neighborhood.py
+++ b/src/neighborhood.py
@@ -4,7 +4,8 @@
 @author: jrpotter
 @date: May 31st, 2015
 """
-import itertools
+import camtools as ct
+import itertools as it
 
 
 class Neighborhood:
@@ -23,45 +24,7 @@ class Neighborhood:
     at this point isn't possible.
     """
 
-    class NeighborhoodKey:
-        """
-        Allows proper sorting of neighborhoods.
-
-        Lists should be returned in order, where cell's with smaller indices (in most significant axis first)
-        are listed before cell's with larger ones. For example, in a 3D grid, the neighbors corresponding to:
-
-            offsets = (-1, -1, -1), (-1, 1, 0), (-1, 0, -1), and (1, 0, -1)
-
-        are returned in the following order:
-
-            offsets = (-1, -1, -1), (-1, 0, -1), (1, 0, -1), (-1, 1, 0)
-
-        since the z-axis is most significant, followed by the y-axis, and lastly the x-axis.
-        """
-        def __init__(self, obj, *args):
-            self.obj = obj
-        def __lt__(self, other):
-             return self.compare(self.obj, other.obj) < 0
-        def __gt__(self, other):
-            return self.compare(self.obj, other.obj) > 0
-        def __eq__(self, other):
-            return self.compare(self.obj, other.obj) == 0
-        def __le__(self, other):
-            return self.compare(self.obj, other.obj) <= 0
-        def __ge__(self, other):
-            return self.compare(self.obj, other.obj) >= 0
-        def __ne__(self, other):
-            return self.compare(self.obj, other.obj) != 0
-        def compare(self, other):
-            for i in reversed(range(len(a))):
-                if a[i] < b[i]:
-                    return -1
-                elif a[i] > b[i]:
-                    return 1
-            return 0
-
-
-    def __init__(self, grid, wrap_around=True):
+    def __init__(self, wrap_around=True):
         """
         Sets up an empty neighborhood.
 
@@ -69,42 +32,41 @@ class Neighborhood:
         Note the offsets have a tuple as a key representing the position being offsetted by, and as a value,
         the current state the given cell at the offset is checked to be.
         """
-        self.grid = grid
         self.offsets = {}
         self.wrap_around = wrap_around
 
 
-    def neighbors(self, cell):
+    def neighbors(self, index, grid):
         """
         Returns all cells in the given neighborhood.
 
-        The returned cells are grouped with the value the cell is checked to be (a 2-tuple (Cell, value) pair).
-        These are sorted based on the NeighborhoodKey comparison class defined above.
+        The returned list of indices represent the index in question, the value at the given index, and
+        the expected value as defined in the offsets.
         """
-        cells = []
-        for k in sorted(self.offsets.keys()):
-            position = [sum(x) for x in zip(cell.index, k)]
-            for i in range(len(position)):
-                if self.wrap_around:
-                    position[i] = position[i] % self.grid.shape[i]
-                elif i < 0 or i >= self.grid.shape[i]:
-                    break
+        indices = []
+        for key in sorted(self.offsets.keys()):
+            if self.wrap_around:
+                f_index = (key + index) % len(grid.flat)
+                indices.append((f_index, grid.flat[f_index], self.offsets[key]))
             else:
-                cells.append((self.grid[tuple(position)], self.offsets[k]))
+                pass
 
-        return cells
+        return indices
 
 
-    def extend(self, offsets, strict=False):
+    def extend(self, offsets, grid, strict=False):
         """
         Adds new offsets to the instance member offsets.
 
         We complain if the strict flag is set to True and an offset has already been declared with a different value.
+        Note also that all offsets are indices of the flattened matrix. This allows for quick row indexing as opposed
+        to individual coordinates.
         """
+        f_offsets = {ct.flatten(k, grid): v for k, v in offsets.items()}
         if not strict:
-            self.offsets.update(offsets)
+            self.offsets.update(f_offsets)
         else:
-            for k in offsets.keys():
+            for k in f_offsets.keys():
                 value = self.offsets.get(k, None)
                 if value is None:
                     self.offsets[k] = offsets[k]
@@ -125,12 +87,12 @@ class Neighborhood:
         """
         offsets = {}
         variants = ([-1, 0, 1],) * len(grid.shape)
-        for current in itertools.product(*variants):
+        for current in it.product(*variants):
             if any(current):
                 offsets[current] = value
 
-        m_neighborhood = cls(grid, wrap_around)
-        m_neighborhood.extend(offsets)
+        m_neighborhood = cls(wrap_around)
+        m_neighborhood.extend(offsets, grid)
 
         return m_neighborhood
 
@@ -154,8 +116,8 @@ class Neighborhood:
                 offsets[tuple(variant)] = value
             variant[i] = 0
 
-        n_neighborhood = cls(grid, wrap_around)
-        n_neighborhood.extend(offsets)
+        n_neighborhood = cls(wrap_around)
+        n_neighborhood.extend(offsets, grid)
 
         return n_neighborhood
 
diff --git a/src/ruleset.py b/src/ruleset.py
index f0667b2..6828cb5 100644
--- a/src/ruleset.py
+++ b/src/ruleset.py
@@ -5,9 +5,8 @@
 @author: jrpotter
 @date: May 31st, 2015
 """
-import copy
 import enum
-import numpy as np
+import camtools
 
 
 class Rule(enum.Enum):
@@ -29,26 +28,6 @@ class Ruleset:
     a neighborhood instance's offsets member.
     """
 
-    @staticmethod
-    @np.vectorize
-    def update(cell, rules, neighborhood, *args):
-        """
-        Allow for batch processing of rules.
-
-        We choose our processing function based on the specified rule and update every cell in the grid simultaneously
-        via a vectorization.
-        """
-        tmp = copy.deepcopy(cell)
-        if rules.method == Rule.MATCH:
-            tmp.value = rules.matches(cell, neighborhood, *args)
-        elif rules.method == Rule.TOLERATE:
-            tmp.value = rules.tolerate(cell, neighborhood, *args)
-        elif rules.method == Rule.SATISFY:
-            tmp.value = rules.satisfies(cell, neighborhood, *args)
-
-        return tmp
-
-
     def __init__(self, method):
         """
 
@@ -56,22 +35,22 @@ class Ruleset:
         self.method = method
 
 
-    def matches(self, cell, neighborhood):
+    def matches(self, index, grid, neighborhood):
         """
         Determines that neighborhood matches expectation exactly.
 
         Note this is just like the tolerate method with a tolerance of 1, but
         recoding allows for short circuiting.
         """
-        residents = neighborhood.neighbors(cell)
+        residents = neighborhood.neighbors(index, grid)
         for resident in residents:
-            if resident[0].value != resident[1]:
+            if grid[resident[0]] != resident[1]:
                 return False
 
         return True
 
 
-    def tolerate(self, cell, neighborhood, tolerance):
+    def tolerate(self, index, grid, neighborhood, tolerance):
         """
         Determines that neighborhood matches expectation within tolerance.
 
@@ -79,23 +58,40 @@ class Ruleset:
         consider this cell to be alive. Note tolerance must be a value 0 <= t <= 1.
         """
         matches = 0
-        residents = neighborhood.neighbors(cell)
+        residents = neighborhood.neighbors(index, grid)
         for resident in residents:
-            if resident[0].value == resident[1]:
+            if grid[resident[0]] == resident[1]:
                 matches += 1
 
         return (matches / len(residents)) >= tolerance
 
 
-    def satisfies(self, cell, neighborhood, valid_func):
+    def satisfies(self, index, grid, neighborhood, valid_func):
         """
         Allows custom function to relay next state of given cell.
 
         The passed function is supplied the list of 2-tuple elements, of which the first is a Cell and the second is
         the expected state as declared in the Neighborhood, as well as the grid and cell in question.
         """
-        residents = neighborhood.neighbors(cell)
+        residents = neighborhood.neighbors(index, grid)
+        coordinate = camtools.unflatten(index, grid)
 
-        return valid_func(cell, residents)
+        return valid_func(coordinate, grid, residents)
 
 
+    def call(self, index, grid, neighborhood, *args):
+        """
+        Allow for batch processing of rules.
+
+        We choose our processing function based on the specified rule and update every cell in the grid simultaneously
+        via a vectorization.
+        """
+        if self.method == Rule.MATCH:
+            func = self.matches
+        elif self.method == Rule.TOLERATE:
+            func = self.tolerate
+        elif self.method == Rule.SATISFY:
+            func = self.satisfies
+
+        return int(func(index, grid, neighborhood, *args))
+