diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) 2014 Brent Baumgartner, Harang Ju, Alex Thomas, Joseph Barrow
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/LambdaNet.cabal b/LambdaNet.cabal
new file mode 100644
--- /dev/null
+++ b/LambdaNet.cabal
@@ -0,0 +1,28 @@
+-- Initial LambdaNet.cabal generated by cabal init.  For further
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+name:                LambdaNet
+version:             0.1.0.0
+synopsis:            A configurable and extensible neural network library
+description:         LambdaNet is an artificial neural network library that allows
+                     users to compose their own networks from function primitives.
+license:             MIT
+license-file:        LICENSE
+author:              Brent Baumgartner, Alex Thomas, Harang Ju, Joseph Barrow
+maintainer:          Joseph Barrow <jdb7hw@virginia.edu>
+copyright:           2014
+category:            Machine Learning
+build-type:          Simple
+cabal-version:       >=1.8
+
+library
+  exposed-modules:     Network.Network, Network.Neuron, Network.Layer, Network.Trainer
+  -- other-modules:
+  build-depends:
+    base ==4.5.*,
+    hmatrix == 0.16.0.6,
+    random,
+    random-shuffle >= 0.0.4,
+    split,
+    binary,
+    bytestring
diff --git a/Network/Layer.hs b/Network/Layer.hs
new file mode 100644
--- /dev/null
+++ b/Network/Layer.hs
@@ -0,0 +1,130 @@
+{-# LANGUAGE FlexibleContexts,
+             RecordWildCards #-}
+
+module Network.Layer
+( LayerDefinition(..)
+, Layer(..)
+, ShowableLayer(..)
+, Connectivity
+, RandomTransform
+
+, layerToShowable
+, showableToLayer
+
+, createLayer
+, connectFully
+, randomList
+, boxMuller
+, normals
+, uniforms
+, boundedUniforms
+) where
+
+import Network.Neuron
+import System.Random
+import Numeric.LinearAlgebra
+import Data.Binary (encode, decode, Binary(..))
+
+-- | The LayerDefinition type is an intermediate type initialized by the
+--   library user to define the different layers of the network.
+data LayerDefinition a = LayerDefinition { neuronDef :: (Neuron a)
+                                         , neuronCount :: Int
+                                         , connect :: (Connectivity a)
+                                         }
+
+-- | The Layer type, which stores the weight matrix, the bias matrix, and
+--   a neuron type.
+data Layer a = Layer { weightMatrix :: (Matrix a)
+                     , biasVector :: (Vector a)
+                     , neuron :: (Neuron a)
+                     }
+
+-- | We have to define a new type to be able to serialize and store
+--   networks.
+data ShowableLayer a = ShowableLayer { weights :: (Matrix a)
+                                     , biases :: (Vector a)
+                                     } deriving Show
+
+-- | We want Showable layer to be packable in the binary format, so we
+--   define it as an instance of showable.
+
+instance (Element a, Binary a) => Binary (ShowableLayer a) where
+  put ShowableLayer{..} = do put weights; put biases
+  get = do weights <- get; biases <- get; return ShowableLayer{..}
+
+-- | Connectivity is the type alias for a function that defines the connective
+--   matrix for two layers (fully connected, convolutionally connected, etc.)
+type Connectivity a = Int -> Int -> Matrix a
+
+-- | A random transformation type alias. It is a transformation defined on an
+--   infinite list of uniformly distributed random numbers, and returns a list
+--   distributed on the transforming distribution.
+type RandomTransform a = [a] -> [a]
+
+-- | The createLayer function takes in a random transformation on an infinite
+--   stream of uniformly generated numbers, a source of entropy, and two
+--   layer definitions, one for the previous layer and one for the next layer.
+--   It returns a layer defined by the Layer type -- a weight matrix, a bias
+--   vector, and a neuron type.
+createLayer ::
+  (RandomGen g, Random a, Floating (Vector a), Container Vector a, Floating a)
+  => RandomTransform a -> g -> LayerDefinition a -> LayerDefinition a -> Layer a
+createLayer t g layerDef layerDef' =
+  Layer (randomMatrix * (connectivity i j))
+        (randomVector * bias)
+        (neuronDef layerDef)
+  where randomMatrix = (i >< j) (randomList t g)
+        randomVector = i |> (randomList t g)
+        i = neuronCount layerDef'
+        j = neuronCount layerDef
+        connectivity = connect layerDef'
+        bias = i |> (repeat 1) -- bias connectivity (full)
+
+-- | The connectFully function takes the number of input neurons for a layer, i,
+--   and the number of output neurons of a layer, j, and returns an i x j
+--   connectivity matrix for a fully connected network.
+connectFully :: Int -> Int -> Matrix Float
+connectFully i j = (i >< j) (repeat 1)
+
+-- | We want to be able to convert between layers and showable layers,
+--   and vice-versa
+layerToShowable :: (Floating (Vector a), Container Vector a, Floating a)
+  => Layer a -> ShowableLayer a
+layerToShowable l = ShowableLayer (weightMatrix l) (biasVector l)
+
+-- | To go from a showable to a layer, we also need a neuron type,
+--   which is an unfortunate restriction owed to Haskell's inability to
+--   serialize functions.
+showableToLayer :: (Floating (Vector a), Container Vector a, Floating a)
+  => (ShowableLayer a, LayerDefinition a) -> Layer a
+showableToLayer (s, d) = Layer (weights s) (biases s) (neuronDef d)
+
+-- | Initialize an infinite random list given a random transform and a source
+--   of entroy.
+randomList :: (RandomGen g, Random a, Floating a)
+  => RandomTransform a -> g -> [a]
+randomList transform = transform . randoms
+
+-- | Define a transformation on the uniform distribution to generate
+--   normally distributed numbers in Haskell (the Box-Muller transform)
+boxMuller :: Floating a => a -> a -> (a, a)
+boxMuller x1 x2 = (z1, z2) where z1 = sqrt ((-2) * log x1) * cos (2 * pi * x2)
+                                 z2 = sqrt ((-2) * log x1) * sin (2 * pi * x2)
+
+-- | This is a function of type RandomTransform that transforms a list of
+--   uniformly distributed numbers to a list of normally distributed numbers.
+normals :: Floating a => [a] -> [a]
+normals (x1:x2:xs) = z1:z2:(normals xs) where (z1, z2) = boxMuller x1 x2
+normals _ = []
+
+-- | A non-transformation to return a list of uniformly distributed numbers
+--   from a list of uniformly distributed numbers. It's really a matter of
+--   naming consistency. It generates numbers on the range (0, 1]
+uniforms :: Floating a => [a] -> [a]
+uniforms xs = xs
+
+-- | An affine transformation to return a list of uniforms on the range
+--   (a, b]
+boundedUniforms :: Floating a => (a, a) -> [a] -> [a]
+boundedUniforms (lower, upper) xs = map affine xs
+  where affine x = lower + x * (upper - lower)
diff --git a/Network/Network.hs b/Network/Network.hs
new file mode 100644
--- /dev/null
+++ b/Network/Network.hs
@@ -0,0 +1,67 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+module Network.Network
+( Network(..)
+, TrainingData
+
+, createNetwork
+, loadNetwork
+, predict
+, apply
+, saveNetwork
+) where
+
+import Network.Neuron
+import Network.Layer
+import System.Random
+import Numeric.LinearAlgebra
+import qualified Data.ByteString.Lazy as B
+import System.IO
+import Data.Binary (encode, decode, Binary(..))
+
+-- | Networks are constructed front to back. Start by adding an input layer,
+--   then each hidden layer, and finally an output layer.
+data Network a = Network { layers :: [Layer a] }
+
+-- | A tuple of (input, expected output)
+type TrainingData a = (Vector a, Vector a)
+
+-- | The createNetwork function takes in a random transform used for weight
+--   initialization, a source of entropy, and a list of layer definitions,
+--   and returns a network with the weights initialized per the random transform.
+createNetwork :: (RandomGen g, Random a, Floating a, Floating (Vector a), Container Vector a)
+  => RandomTransform a -> g -> [LayerDefinition a] -> Network a
+createNetwork t g [] = Network []
+createNetwork t g (layerDef : []) = Network []
+createNetwork t g (layerDef : layerDef' : otherLayerDefs) =
+  Network (layer : layers restOfNetwork)
+  where layer = createLayer t g layerDef layerDef'
+        restOfNetwork = createNetwork t g (layerDef' : otherLayerDefs)
+
+-- | Predict folds over each layer of the network using the input vector as the
+--   first value of the accumulator. It operates on whatever network you pass in.
+predict :: (Floating (Vector a), Container Vector a, Product a)
+  => Vector a -> Network a -> Vector a
+predict input network = foldl apply input (layers network)
+
+-- | A function used in the fold in predict that applies the activation
+--   function and pushes the input through a layer of the network.
+apply :: (Floating (Vector a), Container Vector a, Product a)
+  => Vector a -> Layer a -> Vector a
+apply vector layer = mapVector sigma (weights <> vector + bias)
+  where sigma = activation (neuron layer)
+        weights = weightMatrix layer
+        bias = biasVector layer
+
+-- | Given a filename and a network, we want to save the weights and biases
+--   of the network to the file for later use.
+saveNetwork :: (Binary (ShowableLayer a), Floating a, Floating (Vector a), Container Vector a)
+  => FilePath -> Network a -> IO ()
+saveNetwork file n = B.writeFile file (encode $ map layerToShowable (layers n))
+
+-- | Given a filename, and a list of layer definitions, we want to reexpand
+--   the data back into a network.
+loadNetwork :: (Binary (ShowableLayer a), Floating a, Floating (Vector a), Container Vector a)
+  => FilePath -> [LayerDefinition a] -> IO (Network a)
+loadNetwork file defs = B.readFile file >>= \sls ->
+  return $ Network (map showableToLayer (zip (decode sls) defs))
diff --git a/Network/Neuron.hs b/Network/Neuron.hs
new file mode 100644
--- /dev/null
+++ b/Network/Neuron.hs
@@ -0,0 +1,64 @@
+module Network.Neuron
+( Neuron(..)
+
+, ActivationFunction
+, ActivationFunction'
+, sigmoidNeuron
+, tanhNeuron
+, recluNeuron
+
+, sigmoid
+, sigmoid'
+, tanh
+, tanh'
+, reclu
+, reclu'
+) where
+
+-- | Using this structure allows users of the library to create their own
+--   neurons by creating two functions - an activation function and its
+--   derivative - and packaging them up into a neuron type.
+data Neuron a = Neuron { activation :: (ActivationFunction a)
+                       , activation' :: (ActivationFunction' a)
+                       }
+
+type ActivationFunction a = a -> a
+type ActivationFunction' a = a -> a
+
+-- | Our provided neuron types: sigmoid, tanh, reclu
+sigmoidNeuron :: (Floating a) => Neuron a
+sigmoidNeuron = Neuron sigmoid sigmoid'
+
+tanhNeuron :: (Floating a) => Neuron a
+tanhNeuron = Neuron tanh tanh'
+
+recluNeuron :: (Floating a) => Neuron a
+recluNeuron = Neuron reclu reclu'
+
+-- | The sigmoid activation function, a standard activation function defined
+--   on the range (0, 1).
+sigmoid :: (Floating a) => a -> a
+sigmoid t = 1 / (1 + exp (-1 * t))
+
+-- | The derivative of the sigmoid function conveniently can be computed in
+--   terms of the sigmoid function.
+sigmoid' :: (Floating a) => a -> a
+sigmoid' t = s * (1 - s)
+              where s = sigmoid t
+
+-- | The hyperbolic tangent activation function is provided in Prelude. Here
+--   we provide the derivative. As with the sigmoid function, the derivative
+--   of tanh can be computed in terms of tanh.
+tanh' :: (Floating a) => a -> a
+tanh' t = 1 - s ^ 2
+               where s = tanh t
+
+-- | The rectified linear activation function. This is a more "biologically
+--   accurate" activation function that still retains differentiability.
+reclu :: (Floating a) => a -> a
+reclu t = log (1 + exp t)
+
+-- | The derivative of the rectified linear activation function is just the
+--   sigmoid.
+reclu' :: (Floating a) => a -> a
+reclu' t = sigmoid t
diff --git a/Network/Trainer.hs b/Network/Trainer.hs
new file mode 100644
--- /dev/null
+++ b/Network/Trainer.hs
@@ -0,0 +1,145 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+module Network.Trainer
+( BackpropTrainer(..)
+, CostFunction
+, CostFunction'
+, Selection
+
+, quadraticCost
+, quadraticCost'
+, minibatch
+, online
+, backprop
+, inputs
+, outputs
+, deltas
+, fit
+, evaluate
+) where
+
+import Network.Network
+import Network.Neuron
+import Network.Layer
+import System.Random
+import System.Random.Shuffle (shuffle')
+import Data.List.Split (chunksOf)
+import Numeric.LinearAlgebra
+
+-- | Trainer is a typeclass for all trainer types - a trainer will take in
+--   an instance of itself, a network, a list of training data, and return a
+--   new network trained on the data.
+--class Trainer a where
+--  fit :: (Floating b) => a -> Network b -> [TrainingData b] -> Network b
+
+-- | A BackpropTrainer performs simple backpropagation on a neural network.
+--   It can be used as the basis for more complex trainers.
+data BackpropTrainer a = BackpropTrainer { eta :: a
+                                         , cost :: CostFunction a
+                                         , cost' :: CostFunction' a
+                                         }
+
+-- | A CostFunction is used for evaluating a network's performance on a given
+--   input
+type CostFunction a = Vector a -> Vector a -> a
+
+-- | A CostFunction' (derivative) is used in backPropagation
+type CostFunction' a = Vector a -> Vector a -> Vector a
+
+-- | A selection function for performing gradient descent
+type Selection a = [TrainingData a] -> [[TrainingData a]]
+
+-- | The quadratic cost function (1/2) * sum (y - a) ^ 2
+quadraticCost :: (Floating (Vector a), Container Vector a)
+  => Vector a -> Vector a -> a
+quadraticCost y a = sumElements $ 0.5 * (a - y) ** 2
+
+-- | The derivative of the quadratic cost function sum (y - a)
+quadraticCost' :: (Floating (Vector a))
+  => Vector a -> Vector a -> Vector a
+quadraticCost' y a = a - y
+
+-- | The minibatch function becomes a Selection when partially applied
+--   with the minibatch size
+minibatch :: (Floating (Vector a), Container Vector a)
+  => Int -> [TrainingData a] -> [[TrainingData a]]
+minibatch size = chunksOf size
+
+-- | If we want to train the network online
+online :: (Floating (Vector a), Container Vector a)
+  => [TrainingData a] -> [[TrainingData a]]
+online = minibatch 1
+
+-- | Declare the BackpropTrainer to be an instance of Trainer.
+--instance (Floating a) => Trainer (BackpropTrainer a) where
+fit :: (Floating (Vector a), Container Vector a, Product a)
+  => Selection a -> BackpropTrainer a -> Network a -> [TrainingData a] -> Network a
+fit s t n examples = foldl (backprop t) n $
+  s (shuffle' examples (length examples) (mkStdGen 4))
+
+-- | Perform backpropagation on a single training data instance.
+backprop :: (Floating (Vector a), Container Vector a, Product a)
+  => BackpropTrainer a -> Network a -> [TrainingData a] -> Network a
+backprop t n (e:es) = updateNetwork t n
+  (deltas t n e) (outputs (fst e) n)
+
+-- | Update the weights and biases of a network given a list of deltas
+updateNetwork :: (Floating (Vector a), Container Vector a, Product a)
+  => BackpropTrainer a -> Network a -> [Vector a] -> [Vector a] -> Network a
+updateNetwork t n deltas os =
+  Network $ map (updateLayer t) (zip3 (layers n) deltas os)
+
+-- | The mapped function to update the weight and biases in a single layer
+updateLayer :: (Floating (Vector a), Container Vector a, Product a)
+  => BackpropTrainer a -> (Layer a, Vector a, Vector a) -> Layer a
+updateLayer t (l, delta, output) = Layer newWeight newBias n
+  where n = neuron l
+        newWeight = (weightMatrix l) -
+          (eta t) `scale` ((reshape 1 delta) <> (reshape (dim output) output))
+        newBias = (biasVector l) - (eta t) `scale` delta
+
+-- | The outputs function scans over each layer of the network and stores the
+--   activated results
+outputs :: (Floating (Vector a), Container Vector a, Product a)
+  => Vector a -> Network a -> [Vector a]
+outputs input network = scanl apply input (layers network)
+
+-- | The inputs function performs a similar task to outputs, but returns a list
+--   of vectors of unactivated inputs
+inputs :: (Floating (Vector a), Container Vector a, Product a)
+  => Vector a -> Network a -> [Vector a]
+inputs input network = if null (layers network) then []
+  else unactivated : inputs activated (Network (tail $ layers network))
+    where unactivated = weightMatrix layer <> input + biasVector layer
+          layer = head $ layers network
+          activated = mapVector (activation (neuron layer)) unactivated
+
+-- | The deltas function returns a list of layer deltas.
+deltas :: (Floating (Vector a), Container Vector a, Product a)
+  => BackpropTrainer a -> Network a -> TrainingData a -> [Vector a]
+deltas t n example = hiddenDeltas
+  (Network (reverse (layers n))) outputDelta (tail $ reverse is)
+    ++ [outputDelta]
+  where outputDelta = costd (snd example) output *
+          mapVector activationd lastInput
+        costd = cost' t
+        activationd = activation' (neuron (last (layers n)))
+        output = last os
+        lastInput = last is
+        is = inputs (fst example) n
+        os = outputs (fst example) n
+
+-- | Compute the hidden layer deltas
+hiddenDeltas :: (Floating (Vector a), Container Vector a, Product a)
+  => Network a -> Vector a -> [Vector a] -> [Vector a]
+hiddenDeltas n prevDelta is = if length (layers n) <= 1 then []
+  else delta : hiddenDeltas rest delta (tail is)
+  where rest = Network (tail $ layers n)
+        delta = (trans w) <> prevDelta * spv
+        w = weightMatrix (head $ layers n)
+        spv = mapVector (activation' (neuron (head $ layers n))) (head is)
+
+-- | Use the cost function to determine the error of a network
+evaluate :: (Floating (Vector a), Container Vector a, Product a)
+  => BackpropTrainer a -> Network a -> TrainingData a -> a
+evaluate t n example = (cost t) (snd example) (predict (fst example) n)
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain