diff --git a/src/main.nim b/src/main.nim
index 0f56cdf..e7b3f41 100644
--- a/src/main.nim
+++ b/src/main.nim
@@ -1,13 +1,119 @@
 import nn/network
-import nn/util/activations
-import nn/util/losses
+import nn/util/matrix
+import nn/util/tris
 
 
-const InitialSize = 50
+import std/tables
+import std/math
+import std/random
+import std/algorithm
+
+
+## A bunch of activation functions
+
+func step*(input: float): float = (if input < 0.0: 0.0 else: 1.0)
+func sigmoid*(input: float): float = 1 / (1 + exp(-input))
+func silu*(input: float): float = 1 / (1 + exp(-input))
+func relu*(input: float): float = max(0.0, input)
+func htan*(input: float): float = 
+    let temp = exp(2 * input)
+    result = (temp - 1) / (temp + 1)
+
+
+func ind2sub(n: int, shape: tuple[rows, cols: int]): tuple[row, col: int] =
+    ## Converts an absolute index into an x, y pair
+    return (n div shape.rows, n mod shape.cols)    
+
+
+proc loss(params: TableRef[string, float]): float =
+    ## Our loss function for tris
+    result = params["moves"]
+    if int(params["result"]) == GameStatus.Draw.int:
+        result += 6.0
+    elif int(params["result"]) == GameStatus.Lose.int:
+        result += 12.0
+    result = sigmoid(result)
+
+
+proc compareNetworks(a, b: NeuralNetwork): int =
+    return cmp(loss(a.params), loss(b.params))
+
+
+proc crossover(a, b: NeuralNetwork): NeuralNetwork =
+    result = deepCopy(a)
+    for i, layer in a.layers:
+        result.layers[i].weights = layer.weights.copy()
+        result.layers[i].biases = layer.biases.copy()
+    var i = 0
+    while i < a.layers.len():
+        result.layers[i] = new(Layer)
+        result.layers[i].inputSize = a.layers[i].inputSize
+        result.layers[i].outputSize = a.layers[i].outputSize
+        # We inherit 50% of our weights and biases from our first 
+        # parent and the other 50% from the other parent
+        result.layers[i].weights = where(rand[float](a.layers[i].weights.shape) >= 0.5, a.layers[i].weights, b.layers[i].weights)
+        result.layers[i].biases = where(rand[float](a.layers[i].biases.shape) >= 0.5, a.layers[i].biases, b.layers[i].biases)
+        # Now we sprinkle some mutations into the inherited weights
+        # and biases, just to spice things up. If learnRate = 0.02,
+        # then 2% of our weights and biases will randomly change
+        result.layers[i].weights = where(rand[float](result.layers[i].weights.shape) < a.learnRate, rand[float](result.layers[i].weights.shape), 
+                                         result.layers[i].weights)
+        result.layers[i].biases = where(rand[float](result.layers[i].biases.shape) < a.learnRate, rand[float](result.layers[i].biases.shape), 
+                                         result.layers[i].biases)                           
+        inc(i)
+    result.learnRate = a.learnRate
+
+
+## Our training program
+const Population = 2
+const Iterations = 100
+const Epochs = 10
+const Take = 2
 
 
 var networks: seq[NeuralNetwork] = @[]
-for _ in 0..<InitialSize:
+for _ in 0..<Population:
     networks.add(newNeuralNetwork(@[9, 8, 10, 9], activationFunc=newActivation(sigmoid, func (x, y: float): float = 0.0), 
-                          lossFunc=newLoss(mse, func (x, y: float): float = 0.0), weightRange=(-1.0, +1.0), learnRate=0.05))
+                          lossFunc=newLoss(loss, func (x, y: float): float = 0.0), weightRange=(-1.0, +1.0), learnRate=0.02))
 
+var gameOne: TrisGame
+var gameTwo: TrisGame
+var one: NeuralNetwork
+var two: NeuralNetwork
+var pos: tuple[row, col: int]
+var lost: bool = false
+for epoch in 0..<Epochs:
+    for iteration in 0..<Iterations:
+        gameOne = newTrisGame()
+        gameTwo = newTrisGame()
+        one = sample(networks)
+        two = sample(networks)
+        while one == two:
+            two = sample(networks)
+        while gameOne.get() == Playing:
+            pos = ind2sub(one.compute(gameOne.map.flatten().asType(float)).argmax() - 1, gameOne.map.shape)
+            gameOne.place(TileKind.Self, pos.row, pos.col)
+            gameTwo.place(TileKind.Enemy, pos.row, pos.col)
+            pos = ind2sub(two.compute(gameTwo.map.flatten().asType(float)).argmax() - 1, gameTwo.map.shape)
+            if TileKind(gameOne.map[pos.row, pos.col]) != Empty:
+                # We consider this a loss
+                one.params["result"] = float(Lose)
+                two.params["result"] = float(Lose)
+                lost = true
+                break
+            gameTwo.place(TileKind.Self, pos.row, pos.col)
+            gameOne.place(TileKind.Enemy, pos.row, pos.col)
+        if not lost:
+            one.params["result"] = gameOne.get().float()
+            two.params["result"] = gameTwo.get().float()
+        else:
+            lost = false
+        one.params["moves"] = gameOne.moves.float()
+        two.params["moves"] = gameTwo.moves.float()
+        networks.sort(cmp=compareNetworks)
+        networks = networks[0..<Take]
+        one = sample(networks)
+        two = sample(networks)
+        while one == two:
+            two = sample(networks)
+        networks.add(one.crossover(two))
\ No newline at end of file
diff --git a/src/nn/layer.nim b/src/nn/layer.nim
deleted file mode 100644
index 55e8013..0000000
--- a/src/nn/layer.nim
+++ /dev/null
@@ -1,85 +0,0 @@
-# Copyright 2022 Mattia Giambirtone & All Contributors
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#    http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-import util/matrix
-import util/losses
-import util/activations
-
-import std/strformat
-import std/random
-
-
-randomize()
-
-type 
-    Layer* = ref object
-        ## A generic neural network
-        ## layer
-        inputSize*: int                                     # The number of inputs we process
-        outputSize*: int                                    # The number of outputs we produce
-        weights*: Matrix[float]                             # The weights for each connection (2D)
-        biases*: Matrix[float]                              # The biases for each neuron (1D)
-        activation: Activation                              # The activation function along with its derivative
-        loss: Loss                                          # The cost function used in training
-        gradients: tuple[weights, biases: Matrix[float]]    # Gradient coefficients for weights and biases
-        learnRate: float                                    # The speed at which we perform gradient descent
-
-
-proc `$`*(self: Layer): string =
-    ## Returns a string representation
-    ## of the layer
-    result = &"Layer(inputs={self.inputSize}, outputs={self.outputSize})"
-
-
-proc newLayer*(inputSize: int, outputSize: int, activation: Activation, loss: Loss, learnRate: float, weightRange: tuple[start, stop: float]): Layer =
-    ## Creates a new layer with inputSize input 
-    ## parameters and outputSize outgoing outputs.
-    ## Weights are initialized with random values
-    ## in the chosen range
-    new(result)
-    result.inputSize = inputSize
-    result.outputSize = outputSize
-    var biases = newSeqOfCap[float](outputSize)
-    var biasGradients = newSeqOfCap[float](outputSize)
-    for _ in 0..<outputSize:
-        biases.add(0.0)
-        biasGradients.add(0.0)
-    var weights = newSeqOfCap[seq[float]](inputSize * outputSize)
-    var weightGradients = newSeqOfCap[seq[float]](inputSize * outputSize)
-    for _ in 0..<outputSize:
-        weights.add(@[])
-        weightGradients.add(@[])
-        for _ in 0..<inputSize:
-            weights[^1].add(rand(weightRange.start..weightRange.stop))
-            weightGradients[^1].add(0)
-    result.biases = newMatrix[float](biases)
-    result.weights = newMatrix[float](weights)
-    result.activation = activation
-    result.loss = loss
-    result.gradients = (weights: newMatrix[float](weightGradients), biases: newMatrix[float](biasGradients))
-    result.learnRate = learnRate
-
-
-proc compute*(self: Layer, data: Matrix[float]): Matrix[float] =
-    ## Computes the output of a given layer with
-    ## the given input data and returns it as a
-    ## one-dimensional array
-    result = (self.weights.dot(data).sum() + self.biases).apply(self.activation.function, axis= -1)
-
-
-proc cost*(self: Layer, x, y: Matrix[float]): float =
-    ## Returns the total cost of this layer
-    for i in 0..x.shape.cols:
-        result += self.loss.function(x[0, i], y[0, i])
-    result /= float(x.shape.cols)
-    
diff --git a/src/nn/network.nim b/src/nn/network.nim
index f91c791..b41fd53 100644
--- a/src/nn/network.nim
+++ b/src/nn/network.nim
@@ -12,37 +12,122 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import layer
 import util/matrix
-import util/activations
-import util/losses
-
-export layer, matrix, losses, activations
 
 
-import std/sequtils
 import std/strformat
+import std/tables
+import std/random
 
 
-type 
+randomize()
+
+
+type
     NeuralNetwork* = ref object
         ## A generic feed-forward 
         ## neural network
-        layers: seq[Layer]
+        layers*: seq[Layer]
+        activation: Activation                                 # The activation function along with its derivative
+        loss: Loss                                             # The cost function along with its derivative
+        # This parameter has a different meaning depending on
+        # whether we're learning using backpropagation with gradient 
+        # descent (in which case it is the amount by which we increase 
+        # our input for the next epoch) or using a genetic approach 
+        # (where it will be the rate of mutation for each layer)
+        learnRate*: float
+        # Extra parameters
+        params*: TableRef[string, float]
+    # Note: The derivatives of the loss and activation
+    # function are only meaningful when performing gradient
+    # descent!
+    Loss* = ref object
+        ## A loss function
+        function: proc (params: TableRef[string, float]): float
+        derivative: proc (x, y: float): float {.noSideEffect.}
+    Activation* = ref object
+        ## An activation function
+        function: proc (input: float): float {.noSideEffect.}
+        derivative: proc (x, y: float): float {.noSideEffect.}
+    Layer* = ref object
+        ## A generic neural network
+        ## layer
+        inputSize*: int                                     # The number of inputs we process
+        outputSize*: int                                    # The number of outputs we produce
+        weights*: Matrix[float]                              # The weights for each connection (2D)
+        biases*: Matrix[float]                               # The biases for each neuron (1D)
+        gradients: tuple[weights, biases: Matrix[float]]    # Gradient coefficients for weights and biases, if using gradient descent
 
 
-proc newNeuralNetwork*(layers: seq[int], activationFunc: Activation, lossFunc: Loss, learnRate: float, 
-                       weightRange: tuple[start, stop: float]): NeuralNetwork =
+proc `$`*(self: Layer): string =
+    ## Returns a string representation
+    ## of the layer
+    result = &"Layer(inputs={self.inputSize}, outputs={self.outputSize})"
+
+
+proc `$`*(self: NeuralNetwork): string =
+    ## Returns a string representation
+    ## of the network
+    result = &"NeuralNetwork(learnRate={self.learnRate}, layers={self.layers})"
+
+
+proc newLoss*(function: proc (params: TableRef[string, float]): float, derivative: proc (x, y: float): float {.noSideEffect.}): Loss =
+    ## Creates a new Loss object
+    new(result)
+    result.function = function
+    result.derivative = derivative
+
+
+proc newActivation*(function: proc (input: float): float {.noSideEffect.}, derivative: proc (x, y: float): float {.noSideEffect.}): Activation =
+    ## Creates a new Activation object
+    new(result)
+    result.function = function
+    result.derivative = derivative
+
+
+proc newLayer*(inputSize: int, outputSize: int, weightRange: tuple[start, stop: float]): Layer =
+    ## Creates a new layer with inputSize input 
+    ## parameters and outputSize outgoing outputs.
+    ## Weights are initialized with random values
+    ## in the chosen range
+    new(result)
+    result.inputSize = inputSize
+    result.outputSize = outputSize
+    var biases = newSeqOfCap[float](outputSize)
+    var biasGradients = newSeqOfCap[float](outputSize)
+    for _ in 0..<outputSize:
+        biases.add(0.0)
+        biasGradients.add(0.0)
+    var weights = newSeqOfCap[seq[float]](inputSize * outputSize)
+    var weightGradients = newSeqOfCap[seq[float]](inputSize * outputSize)
+    for _ in 0..<outputSize:
+        weights.add(@[])
+        weightGradients.add(@[])
+        for _ in 0..<inputSize:
+            weights[^1].add(rand(weightRange.start..weightRange.stop))
+            weightGradients[^1].add(0)
+    result.biases = newMatrix[float](biases)
+    result.weights = newMatrix[float](weights)
+    result.gradients = (weights: newMatrix[float](weightGradients), biases: newMatrix[float](biasGradients))
+
+    
+
+proc newNeuralNetwork*(layers: seq[int], activationFunc: Activation, lossFunc: Loss,
+                       learnRate: float, weightRange: tuple[start, stop: float]): NeuralNetwork =
     ## Initializes a new neural network
     ## with the given layer layout
     new(result)
     result.layers = newSeqOfCap[Layer](len(layers))
     for i in 0..<layers.high():
-        result.layers.add(newLayer(layers[i], layers[i + 1], activationFunc, lossFunc, learnRate, weightRange))
-    
+        result.layers.add(newLayer(layers[i], layers[i + 1], weightRange))
+    result.activation = activationFunc
+    result.loss = lossFunc
+    result.learnRate = learnRate
+    result.params = newTable[string, float]()
 
-proc predict*(self: NeuralNetwork, data: Matrix[float]): Matrix[float] =
-    ## Performs a prediction and returns a 1D array
+
+proc compute*(self: NeuralNetwork, data: Matrix[float]): Matrix[float] =
+    ## Performs a computation and returns a 1D array
     ## with the output
     when not defined(release):
         if data.shape.rows > 1:
@@ -51,17 +136,5 @@ proc predict*(self: NeuralNetwork, data: Matrix[float]): Matrix[float] =
             raise newException(ValueError, &"input is of the wrong shape (expecting (1, {self.layers[0].inputSize}), got ({data.shape.rows}, {data.shape.cols}) instead)")
     result = data
     for layer in self.layers:
-        result = layer.compute(result)
+        result = (layer.weights.dot(result).sum() + layer.biases).apply(self.activation.function, axis= -1)
 
-
-proc classify*(self: NeuralNetwork, data: Matrix[float]): int =
-    ## Performs a prediction and returns the label
-    ## with the highest likelyhood
-    result = maxIndex(self.predict(data).raw[])
-
-
-proc cost*(self: NeuralNetwork, x, y: Matrix[float]): float =
-    ## Returns the total average cost of the network
-    for layer in self.layers:
-        result += layer.cost(x, y)
-    result /= float(self.layers.len())
\ No newline at end of file
diff --git a/src/nn/util/activations.nim b/src/nn/util/activations.nim
deleted file mode 100644
index 85b25aa..0000000
--- a/src/nn/util/activations.nim
+++ /dev/null
@@ -1,35 +0,0 @@
-# Copyright 2022 Mattia Giambirtone & All Contributors
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#    http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-import std/math
-
-
-func step*(input: float): float = (if input < 0.0: 0.0 else: 1.0)
-func sigmoid*(input: float): float = 1 / (1 + exp(-input))
-func silu*(input: float): float = 1 / (1 + exp(-input))
-func relu*(input: float): float = max(0.0, input)
-func htan*(input: float): float = 
-    let temp = exp(2 * input)
-    result = (temp - 1) / (temp + 1)
-
-
-type Activation* = ref object
-    function*: proc (input: float): float {.noSideEffect.}
-    derivative*: proc (x, y: float): float {.noSideEffect.}
-
-
-proc newActivation*(function: proc (input: float): float {.noSideEffect.}, derivative: proc (x, y: float): float {.noSideEffect.}): Activation =
-    ## Creates a new activation object
-    new(result)
-    result.function = function
-    result.derivative = derivative
\ No newline at end of file
diff --git a/src/nn/util/losses.nim b/src/nn/util/losses.nim
deleted file mode 100644
index 5244e08..0000000
--- a/src/nn/util/losses.nim
+++ /dev/null
@@ -1,29 +0,0 @@
-# Copyright 2022 Mattia Giambirtone & All Contributors
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#    http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-import std/math
-
-# Mean squared error
-func mse*(x, y: float): float = pow(x - y, 2)
-
-
-type Loss* = ref object
-    function*: proc (x, y: float): float
-    derivative*: proc (x, y: float): float
-
-
-proc newLoss*(function: proc (x, y: float): float, derivative: proc (x, y: float): float): Loss =
-    ## Creates a new activation object
-    new(result)
-    result.function = function
-    result.derivative = derivative
\ No newline at end of file
diff --git a/src/nn/util/matrix.nim b/src/nn/util/matrix.nim
index 5591d58..91101f6 100644
--- a/src/nn/util/matrix.nim
+++ b/src/nn/util/matrix.nim
@@ -132,6 +132,16 @@ proc rand*[T: int | float](shape: tuple[rows, cols: int], order: MatrixOrder = R
             result.data[].add(rand(0.0..1.0))
 
 
+proc asType*[T](self: Matrix[T], kind: typedesc): Matrix[kind] =
+    ## Same as np.array.astype(...)
+    new(result)
+    new(result.data)
+    for e in self.data[]:
+        result.data[].add(kind(e))
+    result.shape = self.shape
+    result.order = self.order
+
+
 # Simple one-line helpers and forward declarations
 func len*[T](self: Matrix[T]): int {.inline.} = self.data[].len()
 func len*[T](self: MatrixView[T]): int {.inline.} = self.shape.cols
@@ -814,9 +824,7 @@ proc dot*[T](self, other: Matrix[T]): Matrix[T] =
             if self.shape.rows != other.shape.cols:
                 raise newException(ValueError, &"incompatible argument shapes for dot product")
         result = zeros[T]((self.shape.rows, other.shape.cols))
-        echo self
         var other = other.transpose()
-        echo other
         for i in 0..<result.shape.rows:
             for j in 0..<result.shape.cols:
                 result[i, j] = (self[i] * other[j]).sum()
@@ -828,13 +836,7 @@ proc dot*[T](self, other: Matrix[T]): Matrix[T] =
         for i in 0..<result.shape.cols:
             result[0, i] = (self[i] * other[0]).sum()
     elif other.shape.rows > 1:
-        when not defined(release):
-            if self.shape.cols != other.shape.cols:
-                raise newException(ValueError, &"incompatible argument shapes for dot product")
-        result = zeros[T]((0, self.shape.cols))
-        var other = other.transpose()
-        for i in 0..<result.shape.cols:
-            result[0, i] = (self[0] * other[i]).sum()
+        return other.transpose().dot(self)
     else:
         return self * other
 
@@ -877,8 +879,8 @@ proc max*[T](self: Matrix[T]): T =
     return m
 
 
-proc argmax*[T](self: Matrix[T]): T =
-    ## Returns the index largest element
+proc argmax*[T](self: Matrix[T]): int =
+    ## Returns the index of largest element
     ## into the matrix
     var m: T = self[0, 0]
     var 
@@ -893,7 +895,7 @@ proc argmax*[T](self: Matrix[T]): T =
             if self[0, col] > m:
                 m = self[0, col]
             inc(col)
-    return m
+    return self.getIndex(row, col)
 
 
 proc contains*[T](self: Matrix[T], e: T): bool =
@@ -905,6 +907,7 @@ proc contains*[T](self: Matrix[T], e: T): bool =
                 return true
     return false
 
+
 when isMainModule:
     import math
 
diff --git a/src/nn/util/tris.nim b/src/nn/util/tris.nim
index a3be29a..dcab9a6 100644
--- a/src/nn/util/tris.nim
+++ b/src/nn/util/tris.nim
@@ -15,12 +15,14 @@ type
         Draw
     TrisGame* = ref object
         map*: Matrix[int]
+        moves*: int
 
 
 proc newTrisGame*: TrisGame =
     ## Creates a new TrisGame object
     new(result)
     result.map = zeros[int]((3, 3))
+    result.moves = 0
 
 
 proc get*(self: TrisGame): GameStatus =
@@ -57,7 +59,9 @@ proc `$`*(self: TrisGame): string =
 
 proc place*(self: TrisGame, tile: TileKind, x, y: int) =
     ## Places a tile onto the playing board
-    self.map[x, y] = int(tile)
+    if TileKind(self.map[x, y]) == Empty:
+        self.map[x, y] = int(tile)
+        inc(self.moves)
 
 
 when isMainModule: