diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Alexander Ignatyev (c) 2016
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Alexander Ignatyev nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,82 @@
+## Machine Learning Toolbox
+
+[![Build Status](https://travis-ci.org/Alexander-Ignatyev/mltool.svg?branch=master)](https://travis-ci.org/Alexander-Ignatyev/mltool)
+[![Coverage Status](https://coveralls.io/repos/github/Alexander-Ignatyev/mltool/badge.svg)](https://coveralls.io/github/Alexander-Ignatyev/mltool)
+[![Documentation](https://img.shields.io/badge/mltool-documentation-blue.svg)](https://alexander-ignatyev.github.io/mltool-docs/doc/index.html)
+
+### Supported Methods and Problems
+
+#### Supervised Learning
+
+##### Regression Problem
+
+* Normal Equation;
+
+* Linear Regression using Least Squares approach.
+
+##### Classification Problem
+
+* Softmax Classifier;
+
+* Multi SVM Classifier;
+
+* Logistic Regression;
+
+* Neural Networks, please see the details below.
+
+#### Unsupervised Learning
+
+* Principal Component Analysis (Dimensionality reduction problem);
+
+* K-Means (Clustering).
+
+#### Neural Networks
+
+* Activations: ReLu, Tanh, Sigmoid;
+
+* Loss Functions: Softmax, Multi SVM, Logistic.
+
+### Usage
+
+#### Build the project
+
+    stack build
+
+#### Run samples app
+
+Please run sample app from root dir (because paths to training data sets are hardcoded).
+
+```bash
+cd samples
+stack build
+stack exec linreg      # Linear Regression Sample App
+stack exec logreg      # Logistic Regression (Classification) Sample App
+stack exec digits      # Muticlass Classification Sample App
+                       # (Recognition of Handwritten Digitts
+stack exec digits-pca  # Apply PCA dimensionaly reduction to digits sample app
+stack exec digits-svm  # Support Vector Machines
+stack exec nn          # Neural Network Sample App
+                       # (Recognition of Handwritten Digits)
+stack exec kmeans      # Clustering Sample App
+```
+
+#### Run unit tests
+
+    stack test
+
+
+### Examples
+
+* Linear Regression: [source code](https://github.com/Alexander-Ignatyev/mltool/blob/master/samples/linear_regression/Main.hs);
+
+* Logistic Regression: [source code](https://github.com/Alexander-Ignatyev/mltool/blob/master/samples/logistic_regression/Main.hs);
+
+* Multiclass Logistic Regression: [source code](https://github.com/Alexander-Ignatyev/mltool/blob/master/samples/digits_classification/Main.hs);
+
+* Multiclass Logistic Regression with PCA: [source code](https://github.com/Alexander-Ignatyev/mltool/blob/master/samples/digits_classification_pca/Main.hs);
+
+* Multiclass Support Vector Machine: [source code](https://github.com/Alexander-Ignatyev/mltool/blob/master/samples/digits_classification_svm/Main.hs);
+
+* Neural Networks: [source code](https://github.com/Alexander-Ignatyev/mltool/blob/master/samples/neural_networks/Main.hs);
+
+* K-Means: [source code](https://github.com/Alexander-Ignatyev/mltool/blob/master/samples/kmeans/Main.hs).
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
diff --git a/mltool.cabal b/mltool.cabal
new file mode 100644
--- /dev/null
+++ b/mltool.cabal
@@ -0,0 +1,104 @@
+name:                mltool
+version:             0.1.0.0
+synopsis:            Machine Learning Toolbox
+description:         Please see README.md
+homepage:            https://github.com/alexander-ignatyev/mltool
+license:             BSD3
+license-file:        LICENSE
+author:              Alexander Ignatyev
+maintainer:          ignatyev.alexander@gmail.com
+copyright:           Alexander Ignatyev
+category:            math
+build-type:          Simple
+extra-source-files:  README.md
+cabal-version:       >=1.10
+
+library
+  hs-source-dirs:      src
+  exposed-modules:     MachineLearning
+                     , MachineLearning.Optimization
+                     , MachineLearning.Optimization.GradientDescent
+                     , MachineLearning.Optimization.MinibatchGradientDescent
+                     , MachineLearning.Regression
+                     , MachineLearning.Model
+                     , MachineLearning.LeastSquaresModel
+                     , MachineLearning.LogisticModel
+                     , MachineLearning.MultiSvmClassifier
+                     , MachineLearning.SoftmaxClassifier
+                     , MachineLearning.Classification.Binary
+                     , MachineLearning.Classification.OneVsAll
+                     , MachineLearning.Classification.MultiClass
+                     , MachineLearning.NeuralNetwork
+                     , MachineLearning.NeuralNetwork.Layer
+                     , MachineLearning.NeuralNetwork.Regularization
+                     , MachineLearning.NeuralNetwork.ReluActivation
+                     , MachineLearning.NeuralNetwork.TanhActivation
+                     , MachineLearning.NeuralNetwork.SigmoidActivation
+                     , MachineLearning.NeuralNetwork.MultiSvmLoss
+                     , MachineLearning.NeuralNetwork.SoftmaxLoss
+                     , MachineLearning.NeuralNetwork.LogisticLoss
+                     , MachineLearning.NeuralNetwork.Topology
+                     , MachineLearning.NeuralNetwork.TopologyMaker
+                     , MachineLearning.NeuralNetwork.WeightInitialization
+                     , MachineLearning.PCA
+                     , MachineLearning.Clustering
+                     , MachineLearning.TerminalProgress
+                     , MachineLearning.Regularization
+                     , MachineLearning.Random
+                     , MachineLearning.Types
+                     , MachineLearning.Utils
+  other-modules:       MachineLearning.Classification.Internal
+  build-depends:       base >= 4.7 && < 5
+                     , vector
+                     , hmatrix
+                     , hmatrix-gsl
+                     , hmatrix-gsl-stats
+                     , ascii-progress
+                     , deepseq
+                     , random
+                     , MonadRandom
+  default-language:    Haskell2010
+
+
+test-suite mltool-test
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             Main.hs
+  other-modules:       MachineLearning.Classification.BinaryTest
+                     , MachineLearning.Classification.OneVsAllTest
+                     , MachineLearning.ClusteringTest
+                     , MachineLearning.DataSets
+                     , MachineLearning.LeastSquaresModelTest
+                     , MachineLearning.LogisticModelTest
+                     , MachineLearning.MultiSvmClassifierTest
+                     , MachineLearning.NeuralNetwork.TopologyTest
+                     , MachineLearning.NeuralNetwork.WeightInitializationTest
+                     , MachineLearning.NeuralNetworkTest
+                     , MachineLearning.Optimization.GradientDescentTest
+                     , MachineLearning.Optimization.MinibatchGradientDescentTest
+                     , MachineLearning.PCATest
+                     , MachineLearning.RandomTest
+                     , MachineLearning.RegressionTest
+                     , MachineLearning.SoftmaxClassifierTest
+                     , MachineLearning.UtilsTest
+                     , MachineLearningTest
+                     , Test.HUnit.Approx
+                     , Test.HUnit.Plus
+  build-depends:       base
+                     , mltool
+                     , vector
+                     , hmatrix
+                     , hmatrix-gsl-stats
+                     , random
+                     , MonadRandom
+                     , test-framework
+                     , test-framework-hunit
+                     , test-framework-quickcheck2
+                     , HUnit
+                     , QuickCheck > 2.0
+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N
+  default-language:    Haskell2010
+
+source-repository head
+  type:     git
+  location: https://github.com/alexander-ignatyev/mltool
diff --git a/src/MachineLearning.hs b/src/MachineLearning.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning.hs
@@ -0,0 +1,83 @@
+{-|
+Module: MachineLearning
+Description: Machine Learning
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+-}
+
+
+module MachineLearning
+(
+  addBiasDimension
+  , removeBiasDimension
+  , meanStddev
+  , featureNormalization
+  , mapFeatures
+  , splitToXY
+)
+
+where
+
+import MachineLearning.Types (Vector, Matrix)
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra((|||), (??))
+import qualified Numeric.GSL.Statistics as Stat
+
+import Control.Monad (replicateM, mfilter, MonadPlus)
+import Data.List (sort, group, foldl')
+import qualified Data.Vector as V
+
+
+-- | Add biad dimension to the future matrix
+addBiasDimension :: Matrix -> Matrix
+addBiasDimension x = 1 ||| x
+
+
+-- | Remove biad dimension
+removeBiasDimension :: Matrix -> Matrix
+removeBiasDimension x = x ?? (LA.All, LA.Drop 1)
+
+
+-- | Caclulates mean and stddev values of every feature.
+-- Takes feature matrix X, returns pair of vectors of means and stddevs.
+meanStddev x =
+  let cols = LA.toColumns x
+      means = map Stat.mean cols
+      stddevs = zipWith (\m col -> Stat.stddev_m m col) means cols
+      stddevs' = map (\s -> if s < 2 then 1 else s) stddevs
+  in (LA.row means, LA.row stddevs')
+
+
+featureNormalization (means, stddevs) x = (x - means) / stddevs
+
+-- | Maps the features into all polynomial terms of X up to the degree-th power
+mapFeatures :: Int -> Matrix -> Matrix
+mapFeatures 1 x = x
+mapFeatures degree x = LA.fromColumns $ cols ++ (foldl' (\l d -> (terms d) ++ l) [] [degree, degree-1 .. 2])
+  where cols = LA.toColumns x
+        vv = V.fromList cols
+        ncols = V.length vv
+        makeTerm :: [(Int, Int)] -> Vector
+        makeTerm = foldl' (\c (index, power) -> c * (vv V.! index) ^ power) 1
+        terms :: Int -> [Vector]
+        terms d = foldl' (\l x -> (makeTerm x) : l) [] $ polynomialTerms d [ncols-1, ncols-2 .. 0]
+        
+
+polynomialTerms :: Ord a => Int -> [a] -> [[(a, Int)]]
+polynomialTerms degree terms =
+  map (\x -> map (\y -> (head y, length y)) $ group x)
+  $ combinationsWithReplacement degree terms
+
+
+combinationsWithReplacement :: (MonadPlus m, Ord a) => Int -> m a -> m [a]
+combinationsWithReplacement sample objects = mfilter (\a -> sort a == a) $ replicateM sample objects
+
+
+-- | Splits data matrix to features matrix X and vector of outputs y
+splitToXY m =
+  let x = m ?? (LA.All, LA.DropLast 1)
+      y = LA.flatten $ m ?? (LA.All, LA.TakeLast 1)
+  in (x, y)
diff --git a/src/MachineLearning/Classification/Binary.hs b/src/MachineLearning/Classification/Binary.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Classification/Binary.hs
@@ -0,0 +1,52 @@
+{-|
+Module: MachineLearning.Classification.Binary
+Description: Binary Classification.
+Copyright: (c) Alexander Ignatyev, 2016-2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Binary Classification.
+-}
+
+module MachineLearning.Classification.Binary
+(
+    Opt.MinimizeMethod(..)
+  , module Log
+  , module Model
+  , predict
+  , learn
+  , MLC.calcAccuracy
+  , Regularization(..)
+)
+
+where
+
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Regularization (Regularization(..))
+import qualified MachineLearning.Optimization as Opt
+import qualified MachineLearning.LogisticModel as Log
+import qualified MachineLearning.Model as Model
+import qualified MachineLearning.Classification.Internal as MLC
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+
+-- | Binary Classification prediction function.
+-- Takes a matrix of features X and a vector theta.
+-- Returns predicted class.
+predict :: Matrix -> Vector -> Vector
+predict x theta = V.map (\r -> if r >= 0.5 then 1 else 0) h
+  where h = Model.hypothesis Log.Logistic x theta
+
+
+-- | Learns Binary Classification.
+learn :: Opt.MinimizeMethod -- ^ (e.g. BFGS2 0.1 0.1)
+         -> R                  -- ^ epsilon, desired precision of the solution;
+         -> Int                -- ^ maximum number of iterations allowed;
+         -> Regularization     -- ^ regularization parameter lambda;
+         -> Matrix             -- ^ matrix X;
+         -> Vector             -- ^ binary vector y;
+         -> Vector             -- ^ initial Theta;
+         -> (Vector, Matrix)   -- ^ solution vector and optimization path.
+learn mm eps numIters lambda x y initialTheta = Opt.minimize mm Log.Logistic eps numIters lambda x y initialTheta
diff --git a/src/MachineLearning/Classification/Internal.hs b/src/MachineLearning/Classification/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Classification/Internal.hs
@@ -0,0 +1,40 @@
+{-|
+Module: MachineLearning.Classification.Internal
+Description: Classification Internal module.
+Copyright: (c) Alexander Ignatyev, 2016-2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Defines Internal Classification functions.
+-}
+
+module MachineLearning.Classification.Internal
+(
+  calcAccuracy
+  , processOutputOneVsAll
+)
+
+where
+
+import MachineLearning.Types (R, Vector)
+import qualified Data.Vector.Storable as V
+
+
+-- | Calculates accuracy of Classification predictions.
+-- Takes vector expected y and vector predicted y.
+-- Returns number from 0 to 1, the closer to 1 the better accuracy.
+-- Suitable for both Classification Types: Binary and Multiclass.
+calcAccuracy :: Vector -> Vector -> R
+calcAccuracy yExpected yPredicted = (1 - (V.sum discrepancy) / (fromIntegral $ V.length discrepancy))
+  where discrepancy = V.zipWith f yExpected yPredicted
+        f y1 y2 = if round y1 == round y2 then 0 else 1
+
+
+-- | Process outputs for One-vs-All Classification.
+-- Takes number of labels and output vector y.
+-- Returns list of vectors of binary outputs (One-vs-All Classification).
+-- It is supposed that labels are integerets start at 0.
+processOutputOneVsAll :: Int -> Vector -> [Vector]
+processOutputOneVsAll numLabels y = map f [0 .. numLabels-1]
+  where f sample = V.map (\a -> if round a == sample then 1 else 0) y
diff --git a/src/MachineLearning/Classification/MultiClass.hs b/src/MachineLearning/Classification/MultiClass.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Classification/MultiClass.hs
@@ -0,0 +1,96 @@
+{-|
+Module: MachineLearning.Classification.MultiClass
+Description: MultiClass Classification.
+Copyright: (c) Alexander Ignatyev, 2016-2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+MultiClass Classification.
+-}
+
+module MachineLearning.Classification.MultiClass
+(
+  Classifier(..)
+  , MultiClassModel(..)
+  , processOutput
+  , Regularization(..)
+  , ccostReg
+  , cgradientReg
+)
+
+where
+
+import MachineLearning.Types
+import MachineLearning.Model
+import MachineLearning.Regularization (Regularization(..))
+import qualified MachineLearning as ML
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+
+-- | Classifier type class represents Multi-class classification models.
+class Classifier a where
+  -- | Score function
+  cscore :: a -> Matrix -> Matrix -> Matrix
+
+  -- | Hypothesis function
+  -- Takes X (m x n) and theta (n x k), returns y (m x k).
+  chypothesis :: a -> Matrix -> Matrix -> Vector
+  
+  -- | Cost function J(Theta), a.k.a. loss function.
+  -- It takes regularizarion parameter lambda, matrix X (m x n), vector y (m x 1) and vector theta (n x 1).
+  ccost :: a -> Regularization -> Matrix -> Vector -> Matrix -> R
+
+  -- | Gradient function.
+  -- It takes regularizarion parameter lambda, X (m x n), y (m x 1) and theta (n x 1).
+  -- Returns vector of gradients (n x 1).
+  cgradient :: a -> Regularization -> Matrix -> Vector -> Matrix -> Matrix
+
+  -- | Returns Number of Classes
+  cnumClasses :: a -> Int
+
+
+-- | MultiClassModel is Model wrapper class around Classifier
+data MultiClassModel m = MultiClass m
+
+
+instance (Classifier a) => Model (MultiClassModel a) where
+  hypothesis (MultiClass m) x theta = chypothesis m x theta'
+    where theta' = unflatten (cnumClasses m) theta
+
+  cost (MultiClass m) lambda x y theta = ccost m lambda x y theta'
+    where theta' = unflatten (cnumClasses m) theta
+
+  gradient (MultiClass m) lambda x y theta = LA.flatten $ cgradient m lambda x y theta'
+    where theta' = unflatten (cnumClasses m) theta
+
+
+unflatten :: Int -> Vector -> Matrix
+unflatten nLabels v = LA.reshape cols v
+  where cols = (V.length v) `div` nLabels
+
+
+-- | Process outputs for MultiClass Classification.
+-- Takes Classifier and output vector y.
+-- Returns matrix of binary outputs.
+-- It is supposed that labels are integerets start at 0.
+processOutput :: (Classifier c) => c -> Vector -> Matrix
+processOutput c y = LA.fromColumns $ map f [0 .. (cnumClasses c)-1]
+  where f sample = V.map (\a -> if round a == sample then 1 else 0) y
+
+
+-- | Calculates regularization for Classifier.ccost.
+-- It takes regularization parameter and theta.
+ccostReg :: Regularization -> Matrix -> R
+ccostReg RegNone _ = 0
+ccostReg (L2 lambda) theta = (LA.norm_2 thetaReg) * 0.5 * lambda
+  where thetaReg = ML.removeBiasDimension theta
+
+
+-- | Calculates regularization for Classifier.cgradient.
+-- It takes regularization parameter and theta.
+cgradientReg :: Regularization -> Matrix -> Matrix
+cgradientReg RegNone _ = 0
+cgradientReg (L2 lambda) theta = ((LA.scalar lambda) * thetaReg)
+  where thetaReg = 0 LA.||| (ML.removeBiasDimension theta)
diff --git a/src/MachineLearning/Classification/OneVsAll.hs b/src/MachineLearning/Classification/OneVsAll.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Classification/OneVsAll.hs
@@ -0,0 +1,58 @@
+{-|
+Module: MachineLearning.Classification.OneVsAll
+Description: One-vs-All Classification.
+Copyright: (c) Alexander Ignatyev, 2016-2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+One-vs-All Classification.
+-}
+
+module MachineLearning.Classification.OneVsAll
+(
+    Opt.MinimizeMethod(..)
+  , module Log
+  , module Model
+  , predict
+  , learn
+  , MLC.calcAccuracy
+  , Regularization(..)
+)
+
+where
+
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Regularization (Regularization(..))
+import qualified MachineLearning.Optimization as Opt
+import qualified MachineLearning.LogisticModel as Log
+import qualified MachineLearning.Model as Model
+import qualified MachineLearning.Classification.Internal as MLC
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+
+-- | One-vs-All Classification prediction function.
+-- Takes a matrix of features X and a list of vectors theta,
+-- returns predicted class number assuming that class numbers start at 0.
+predict :: Matrix -> [Vector] -> Vector
+predict x thetas = predictions'
+  where predict = Model.hypothesis Log.Logistic x
+        predictions = LA.toRows . LA.fromColumns $ map predict thetas
+        predictions' = LA.vector $ map (fromIntegral . LA.maxIndex) predictions
+
+
+-- | Learns One-vs-All Classification
+learn :: Opt.MinimizeMethod -- ^ (e.g. BFGS2 0.1 0.1)
+         -> R                  -- ^ epsilon, desired precision of the solution;
+         -> Int                -- ^ maximum number of iterations allowed;
+         -> Regularization     -- ^ regularization parameter lambda;
+         -> Int                -- ^ number of labels
+         -> Matrix             -- ^ matrix X;
+         -> Vector             -- ^ vector y
+         -> [Vector]             -- ^ initial theta list;
+         -> ([Vector], [Matrix])  -- ^ solution vector and optimization path.
+learn mm eps numIters lambda nLabels x y initialThetaList =
+  let ys = MLC.processOutputOneVsAll nLabels y
+      minimize = Opt.minimize mm Log.Logistic eps numIters lambda x
+  in unzip $ zipWith minimize ys initialThetaList
diff --git a/src/MachineLearning/Clustering.hs b/src/MachineLearning/Clustering.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Clustering.hs
@@ -0,0 +1,123 @@
+{-|
+Module: MachineLearning.Clustering
+Description: Clustering
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Cluster Analysis a.k.a. Clustering - grouping data into coherent subsets.
+-}
+
+module MachineLearning.Clustering
+(
+  Cluster(..)
+  , kmeans
+    
+  -- * Exported for testing purposes only.
+  , kmeansIterM
+)
+
+where
+
+import MachineLearning.Types (R, Vector, Matrix)
+import Data.List (sortOn, groupBy, minimumBy)
+import Control.Applicative ((<$>))
+import Control.Monad (forM)
+import qualified Control.Monad.Random as RndM
+import qualified Data.Vector as V
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Random (sampleM)
+
+
+-- | Cluster type (list of samples associtaed with the cluster).
+type Cluster = V.Vector Vector
+
+
+-- | Gets list if the nearest centroid to the sample.
+nearestCentroidIndex :: V.Vector Vector  -- ^ list of cluster cetroids;
+                     -> Vector            -- ^ sample;
+                     -> Int               -- ^ index of the nearest centroid.
+nearestCentroidIndex centroids v =
+  let distances = V.map (\centroid -> LA.norm_2 (v-centroid)) centroids
+  in V.minIndex distances
+
+
+-- | Calculates cost associated with one claster.
+calcClusterCost :: Cluster  -- ^ cluster;
+                -> Vector   -- ^ cluster centroid;
+                -> R        -- ^ cost value.
+calcClusterCost cluster centroid = sum $ fmap (\sample -> LA.norm_2 $ sample-centroid) cluster
+
+
+-- | Calculates cost value for all clusters.
+calcCost :: V.Vector Cluster  -- ^ cluster list;
+         -> V.Vector Vector   -- ^ list of cluster centroids;
+         -> R                  -- ^ cost value.
+calcCost clusters centroids = sum $ V.zipWith calcClusterCost clusters centroids
+
+
+-- | Calculates centroid (middle point) of the given cluster.
+getNewCentroid :: Cluster      -- ^ cluster;
+               -> Vector       -- ^ centroid.
+getNewCentroid cluster =
+  let n = length cluster
+      centroid = (sum cluster) / (fromIntegral n)
+  in centroid
+
+
+-- | Calculates new cluster centroids for each cluster.
+moveCentroids :: V.Vector Cluster    -- ^ list of clusters;
+              -> V.Vector Vector     -- ^ list of cluster centroids.
+moveCentroids clusters = fmap getNewCentroid clusters
+
+
+-- | Build cluster list from list of clusters indices.
+buildClusterList :: V.Vector Vector   -- ^ list of samples;
+                 -> V.Vector Int      -- ^ list of cluster indices (associated cluster index for each sample);
+                 -> V.Vector Cluster  -- ^ list of clusters.
+buildClusterList samples clusterIndicesList = V.fromList $ fmap getClusterSamples clusters''
+  where clusters' = groupBy (\l r -> snd l == snd r) $ sortOn snd $ zip [0..] $ V.toList clusterIndicesList
+        clusters'' = map (map fst) clusters'
+        getClusterSamples clusterIndices = V.fromList $ fmap (samples V.!) clusterIndices
+
+
+-- -- | Run K-Means algorithm once.
+kmeansIter :: V.Vector Vector           -- ^ list of samples;
+              -> Int                    -- ^ number of clusters (`K`);
+              -> V.Vector Vector        -- ^ list of initial centroids;
+              -> (V.Vector Cluster, [R])  -- ^ (list of clusters, cost values).
+kmeansIter samples k initialCentroids =
+  let iter centroids js =
+        let clusterIndicesList = fmap (nearestCentroidIndex centroids) samples
+            clusters = buildClusterList samples clusterIndicesList
+            centroids' = moveCentroids clusters
+            j = calcCost clusters centroids'
+            diff = sum . fmap LA.norm_2 $ V.zipWith (-) centroids centroids'
+        in if diff < 0.001 then (clusters, j:js)
+           else iter centroids' (j:js)
+  in iter initialCentroids []
+
+
+-- | Run K-Means algorithm once inside Random Monad.
+kmeansIterM :: RndM.RandomGen g =>
+               V.Vector Vector  -- ^ list of samples;
+               -> Int           -- ^ number of clusters (`K`);
+               -> Int           -- ^ iteration number;
+               -> RndM.Rand g (V.Vector Cluster, [R])  -- ^ (list of clusters, cost values) inside Random Monad.
+kmeansIterM samples k _ = do
+  centroids <- sampleM k samples
+  return (kmeansIter samples k centroids)
+
+
+-- | Clusters data using K-Means Algorithm inside Random Monad.
+-- Runs K-Means algorithm `N` times, returns the clustering with smaller cost.
+kmeans :: RndM.RandomGen g =>
+           Int                     -- ^ number of K-Means Algorithm runs (`N`);
+           -> Matrix                  -- ^ data to cluster;
+           -> Int                     -- ^ desired number of clusters (`K`);
+           -> RndM.Rand g (V.Vector Cluster)  -- ^ list of clusters inside Random Monad.
+kmeans nIters x k = fst <$>
+    (minimumBy (\(_, js1) (_, js2) -> compare (head js1) (head js2))) <$>
+    forM [1..nIters] (kmeansIterM samples k)
+  where samples = V.fromList $ LA.toRows x
diff --git a/src/MachineLearning/LeastSquaresModel.hs b/src/MachineLearning/LeastSquaresModel.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/LeastSquaresModel.hs
@@ -0,0 +1,40 @@
+{-|
+Module: MachineLearning.LeastSquaresModel
+Description: Least Squares Model
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+-}
+
+module MachineLearning.LeastSquaresModel
+(
+  LeastSquaresModel(..)
+)
+
+where
+
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra((<>), (#>), (<.>))
+import qualified Numeric.LinearAlgebra.Data as LAD
+import qualified Data.Vector.Storable as V
+
+import qualified MachineLearning.Regularization as R
+
+import MachineLearning.Model
+
+data LeastSquaresModel = LeastSquares
+
+instance Model LeastSquaresModel where
+  hypothesis LeastSquares x theta = x #> theta
+
+  cost LeastSquares lambda x y theta = 
+    let m = x #> theta - y
+        nExamples = fromIntegral $ LA.rows x
+        regTerm = R.costReg lambda theta
+    in (LA.sumElements (m * m) * 0.5 + regTerm) / nExamples
+
+  gradient LeastSquares lambda x y theta = ((LA.tr x) #> (x #> theta - y) + regTerm) / nExamples
+    where nExamples = fromIntegral $ LAD.rows x
+          regTerm = R.gradientReg lambda theta
diff --git a/src/MachineLearning/LogisticModel.hs b/src/MachineLearning/LogisticModel.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/LogisticModel.hs
@@ -0,0 +1,57 @@
+{-|
+Module: MachineLearning.LogisticModel
+Description: Logistic Regression Model
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+-}
+
+module MachineLearning.LogisticModel
+(
+  module MachineLearning.Model
+  , LogisticModel(..)
+  , sigmoid
+  , sigmoidGradient
+)
+
+where
+
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra((<>), (#>), (<.>))
+import qualified Data.Vector.Storable as V
+
+import MachineLearning.Model
+import qualified MachineLearning.Regularization as R
+
+data LogisticModel = Logistic
+
+
+-- | Calculates sigmoid
+sigmoid :: Floating a => a -> a
+sigmoid z = 1 / (1+exp(-z))
+
+
+-- | Calculates derivatives of sigmoid
+sigmoidGradient :: Floating a => a -> a
+sigmoidGradient z = s * (1-s)
+  where s = sigmoid z
+
+
+instance Model LogisticModel where
+  hypothesis Logistic x theta = sigmoid (x #> theta)
+
+  cost m lambda x y theta =
+    let h = hypothesis m x theta
+        nExamples = fromIntegral $ LA.rows x
+        tau = 1e-7
+        jPositive = log(tau + h) <.> (-y)
+        jNegative = log((1 + tau) - h) <.> (1-y)
+        regTerm = R.costReg lambda theta
+    in (jPositive - jNegative + regTerm) / nExamples
+
+  gradient m lambda x y theta = (((LA.tr x) #> (h  - y)) + regTerm) / nExamples
+    where h = hypothesis m x theta
+          nExamples = fromIntegral $ LA.rows x
+          regTerm = R.gradientReg lambda theta
diff --git a/src/MachineLearning/Model.hs b/src/MachineLearning/Model.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Model.hs
@@ -0,0 +1,36 @@
+{-|
+Module: MachineLearning.Model
+Description: Regression Model
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Regression Model type class.
+
+-}
+
+module MachineLearning.Model
+(
+  Model(..)
+)
+
+where
+
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Regularization (Regularization)
+
+class Model a where
+  -- | Hypothesis function, a.k.a. score function (for lassifition problem)
+  -- Takes X (m x n) and theta (n x 1), returns y (m x 1).
+  hypothesis :: a -> Matrix -> Vector -> Vector
+  
+  -- | Cost function J(Theta), a.k.a. loss function.
+  -- It takes regularizarion parameter, matrix X (m x n), vector y (m x 1) and vector theta (n x 1).
+  cost :: a -> Regularization -> Matrix -> Vector -> Vector -> R
+
+  -- | Gradient function.
+  -- It takes regularizarion parameter, X (m x n), y (m x 1) and theta (n x 1).
+  -- Returns vector of gradients (n x 1).
+  gradient :: a -> Regularization -> Matrix -> Vector -> Vector -> Vector
+
diff --git a/src/MachineLearning/MultiSvmClassifier.hs b/src/MachineLearning/MultiSvmClassifier.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/MultiSvmClassifier.hs
@@ -0,0 +1,66 @@
+{-|
+Module: MachineLearning.MultiSvmClassifier
+Description: Multiclass Support Vector Machines Classifier.
+Copyright: (c) Alexander Ignatyev, 2017.
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Multicalss Support Vector Machines Classifier.
+-}
+
+module MachineLearning.MultiSvmClassifier
+(
+  module MachineLearning.Model
+  , module MachineLearning.Classification.MultiClass
+  , MultiSvmClassifier(..)
+)
+
+where
+
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra((<>), (<.>), (|||))
+import qualified Data.Vector.Storable as V
+
+import qualified MachineLearning as ML
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Utils (sumByRows, reduceByRowsV)
+import MachineLearning.Model
+import MachineLearning.Classification.MultiClass
+
+
+-- | Multiclass SVM Classifier, takes delta and number of futures. Delta = 1.0 is good for all cases.
+data MultiSvmClassifier = MultiSvm R Int
+
+
+instance Classifier MultiSvmClassifier where
+  cscore (MultiSvm _ _) x theta = x <> (LA.tr theta)
+
+  chypothesis m x theta = predictions
+    where scores = cscore m x theta
+          predictions = reduceByRowsV (fromIntegral . LA.maxIndex) scores
+
+  ccost m@(MultiSvm d _) lambda x y theta =
+    let nSamples = fromIntegral $ LA.rows x
+        scores = cscore m x theta
+        correct_scores = LA.remap (LA.asColumn $ V.fromList [0..(fromIntegral $ LA.rows x)-1]) (LA.toInt $ LA.asColumn y) scores
+        margins = scores - (correct_scores - (LA.scalar d))
+        margins' = margins * LA.step margins
+        loss = LA.sumElements(margins') / nSamples - d
+        reg = (ccostReg lambda theta) / nSamples
+    in loss + reg
+
+  cgradient m@(MultiSvm d _) lambda x y theta =
+    let nSamples = fromIntegral $ LA.rows x
+        ys = processOutput m y
+        scores = cscore m x theta
+        correct_scores = LA.remap (LA.asColumn $ V.fromList [0..(fromIntegral $ LA.rows x)-1]) (LA.toInt $ LA.asColumn y) scores
+        margins = scores - (correct_scores - (LA.scalar d))
+        margins' = (1-ys)*(LA.step margins)  -- step == cmap (\x -> if x>0 then 1 else 0)
+        k = sumByRows margins'
+        margins'' = margins' - (ys * k)
+        dw = ((LA.tr margins'') <> x) / nSamples
+        reg = (cgradientReg lambda theta) / nSamples
+     in dw + reg
+
+  cnumClasses (MultiSvm _ nLabels) = nLabels
diff --git a/src/MachineLearning/NeuralNetwork.hs b/src/MachineLearning/NeuralNetwork.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork.hs
@@ -0,0 +1,68 @@
+{-|
+Module: MachineLearning.NeuralNetwork
+Description: Neural Network
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Simple Neural Networks.
+-}
+
+module MachineLearning.NeuralNetwork
+(
+    Model(..)
+  , NeuralNetworkModel(..)
+  , MLC.calcAccuracy
+  , T.Topology
+  , T.initializeTheta
+  , T.initializeThetaIO
+  , T.initializeThetaM
+  , Regularization(..)
+)
+
+where
+
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Utils (reduceByRowsV)
+import qualified MachineLearning.Classification.Internal as MLC
+import MachineLearning.Model (Model(..))
+import qualified MachineLearning.NeuralNetwork.Topology as T
+import MachineLearning.Regularization (Regularization(..))
+
+
+-- | Neural Network Model.
+-- Takes neural network topology as a constructor argument.
+newtype NeuralNetworkModel = NeuralNetwork T.Topology
+
+
+instance Model NeuralNetworkModel where
+  hypothesis (NeuralNetwork topology) x theta = predictions
+    where thetaList = T.unflatten topology theta
+          scores = calcScores topology x thetaList
+          predictions = reduceByRowsV (fromIntegral . LA.maxIndex) scores
+
+  cost (NeuralNetwork topology) lambda x y theta =
+    let (ys, thetaList) = processParams topology y theta
+        scores = calcScores topology x thetaList
+    in T.loss topology lambda scores thetaList ys
+
+  gradient (NeuralNetwork topology) lambda x y theta =
+    let (ys, thetaList) = processParams topology y theta
+        (scores, cacheList) = T.propagateForward topology x thetaList
+        grad = T.flatten $ T.propagateBackward topology lambda scores cacheList ys
+    in grad
+
+
+-- | Score function. Takes a topology, X and theta list.
+calcScores :: T.Topology -> Matrix -> [(Matrix, Matrix)] -> Matrix
+calcScores topology x thetaList = fst $ T.propagateForward topology x thetaList
+
+
+processParams :: T.Topology -> Vector -> Vector -> (Matrix, [(Matrix, Matrix)])
+processParams topology y theta =
+  let nOutputs = T.numberOutputs topology
+      ys = LA.fromColumns $ MLC.processOutputOneVsAll nOutputs y
+      thetaList = T.unflatten topology theta
+  in (ys, thetaList)
diff --git a/src/MachineLearning/NeuralNetwork/Layer.hs b/src/MachineLearning/NeuralNetwork/Layer.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/Layer.hs
@@ -0,0 +1,57 @@
+{-# LANGUAGE RankNTypes #-}
+{-|
+Module: MachineLearning.NeuralNetwork.Layer
+Description: Neural Network's Layer
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Neural Network's Layer
+-}
+
+
+module MachineLearning.NeuralNetwork.Layer
+(
+  Layer(..)
+  , Cache(..)
+  , affineForward
+  , affineBackward
+)
+
+where
+
+
+import MachineLearning.Types (R, Matrix)
+import MachineLearning.Utils (sumByColumns)
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra ((<>))
+import qualified Control.Monad.Random as RndM
+
+
+data Cache = Cache {
+  cacheZ :: Matrix
+  , cacheX :: Matrix
+  , cacheW :: Matrix
+  };
+
+
+data Layer = Layer {
+  lUnits :: Int
+  , lForward :: Matrix -> Matrix -> Matrix -> Matrix
+  , lBackward :: Matrix -> Cache -> (Matrix, Matrix, Matrix)
+  , lActivation :: Matrix -> Matrix
+  , lActivationGradient :: Matrix -> Matrix -> Matrix
+  , lInitializeThetaM :: forall g. RndM.RandomGen g => (Int, Int) -> RndM.Rand g (Matrix, Matrix)
+  }
+
+
+affineForward :: Matrix -> Matrix -> Matrix -> Matrix
+affineForward x b w = (x <> LA.tr w) + b
+
+
+affineBackward delta (Cache _ x w) = (dx, db, dw)
+  where m = fromIntegral $ LA.rows x
+        dx = delta <> w
+        db = (sumByColumns delta)/m
+        dw = (LA.tr delta <> x)/m
diff --git a/src/MachineLearning/NeuralNetwork/LogisticLoss.hs b/src/MachineLearning/NeuralNetwork/LogisticLoss.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/LogisticLoss.hs
@@ -0,0 +1,39 @@
+{-|
+Module: MachineLearning.NeuralNetwork.LogisticLoss
+Description: Logistic Loss.
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Logistic Loss.
+-}
+
+module MachineLearning.NeuralNetwork.LogisticLoss
+(
+  scores
+  , gradient
+  , loss
+)
+
+where
+
+
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Types (R, Matrix)
+import qualified MachineLearning.LogisticModel as LM
+
+
+scores :: Matrix -> Matrix
+scores = LM.sigmoid
+
+
+gradient :: Matrix -> Matrix -> Matrix
+gradient scores y = scores - y
+
+
+-- Logistic Loss function
+loss :: Matrix -> Matrix -> R
+loss x y = (LA.sumElements $ (-y) * log(tau + x) - (1-y) * log ((1+tau)-x))/m
+  where tau = 1e-7
+        m = fromIntegral $ LA.rows x
diff --git a/src/MachineLearning/NeuralNetwork/MultiSvmLoss.hs b/src/MachineLearning/NeuralNetwork/MultiSvmLoss.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/MultiSvmLoss.hs
@@ -0,0 +1,52 @@
+{-|
+Module: MachineLearning.NeuralNetwork.MultiSvmLoss
+Description: Multi SVM Loss.
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Multi SVM Loss.
+-}
+
+module MachineLearning.NeuralNetwork.MultiSvmLoss
+(
+  scores
+  , gradient
+  , loss
+)
+
+where
+
+
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Types (R, Matrix)
+import MachineLearning.Utils (sumByRows, reduceByRows)
+
+
+-- | SVM Delta
+svmD :: R
+svmD = 1.0
+
+
+scores = id
+
+
+gradient :: Matrix -> Matrix -> Matrix
+gradient scores y =
+    let correct_scores = sumByRows $ scores*(LA.step y)
+        margins = scores - (correct_scores - (LA.scalar svmD))
+        margins' = (1-y)*(LA.step margins)
+        k = sumByRows margins'
+    in margins' - (y * k)
+
+
+loss :: Matrix -> Matrix -> R
+loss scores y = 
+  let nSamples = fromIntegral $ LA.rows scores
+      correct_scores = sumByRows $ scores*(LA.step y)
+      margins = scores - (correct_scores - (LA.scalar svmD))
+      margins' = margins * (LA.step margins)
+      loss = LA.sumElements margins' - nSamples * svmD
+  in loss / nSamples
diff --git a/src/MachineLearning/NeuralNetwork/Regularization.hs b/src/MachineLearning/NeuralNetwork/Regularization.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/Regularization.hs
@@ -0,0 +1,38 @@
+{-|
+Module: MachineLearning.NeuralNetwork.Regularization
+Description: Regularization
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Regularization.
+-}
+
+module MachineLearning.NeuralNetwork.Regularization
+(
+  Regularization(..)
+  , forwardReg
+  , backwardReg
+)
+
+where
+
+
+import MachineLearning.Types (R, Matrix)
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Regularization (Regularization(..))
+
+
+-- | Calcaulates regularization for forward propagation.
+-- It takes regularization parameter and theta list.
+forwardReg :: Regularization -> [(Matrix, Matrix)] -> R
+forwardReg RegNone _ = 0
+forwardReg (L2 lambda) thetaList = 0.5 * lambda * (sum $ map LA.norm_2 $ snd $ unzip thetaList)
+
+
+-- | Calculates regularization for step of backward propagation.
+-- It takes regularization parameter and theta.
+backwardReg :: Regularization -> Matrix -> Matrix
+backwardReg RegNone _ = 0
+backwardReg (L2 lambda) w = w * (LA.scalar lambda)
diff --git a/src/MachineLearning/NeuralNetwork/ReluActivation.hs b/src/MachineLearning/NeuralNetwork/ReluActivation.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/ReluActivation.hs
@@ -0,0 +1,30 @@
+{-|
+Module: MachineLearning.NeuralNetwork.ReluActivation
+Description: ReLu Activation
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+ReLu Activation.
+-}
+
+module MachineLearning.NeuralNetwork.ReluActivation
+(
+  relu
+  , gradient
+)
+
+where
+
+
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Types (Matrix)
+
+
+relu :: Matrix -> Matrix
+relu x = x * (LA.step x)
+
+
+gradient :: Matrix -> Matrix -> Matrix
+gradient x dx = dx * (LA.step x)  -- == dx[x<0] = 0
diff --git a/src/MachineLearning/NeuralNetwork/SigmoidActivation.hs b/src/MachineLearning/NeuralNetwork/SigmoidActivation.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/SigmoidActivation.hs
@@ -0,0 +1,26 @@
+{-|
+Module: MachineLearning.NeuralNetwork.SigmoidActivation
+Description: Sigmoid
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Sigmoid Activation.
+-}
+
+module MachineLearning.NeuralNetwork.SigmoidActivation
+(
+    LM.sigmoid
+    , gradient
+)
+
+where
+
+
+import MachineLearning.Types (R, Matrix)
+import qualified MachineLearning.LogisticModel as LM
+
+
+gradient :: Matrix -> Matrix -> Matrix
+gradient z da = da * LM.sigmoidGradient z
diff --git a/src/MachineLearning/NeuralNetwork/SoftmaxLoss.hs b/src/MachineLearning/NeuralNetwork/SoftmaxLoss.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/SoftmaxLoss.hs
@@ -0,0 +1,42 @@
+{-|
+Module: MachineLearning.NeuralNetwork.SoftmaxLoss
+Description: Softmax Loss.
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Softmax Loss.
+-}
+
+module MachineLearning.NeuralNetwork.SoftmaxLoss
+(
+  scores
+  , gradient
+  , loss
+)
+
+where
+
+
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Types (R, Matrix)
+import MachineLearning.Utils (sumByRows, reduceByRows)
+
+scores x = x - reduceByRows V.maximum x
+
+
+gradient scores y =
+  let sum_probs = sumByRows $ exp scores
+      probs = (exp scores) / sum_probs
+  in probs - y
+
+
+-- Softmax Loss function
+loss :: Matrix -> Matrix -> R
+loss scores y = loss / m
+  where m = fromIntegral $ LA.rows scores
+        sum_probs = sumByRows $ exp scores
+        t = sumByRows $ scores * y
+        loss = LA.sumElements $ (log sum_probs) - t
diff --git a/src/MachineLearning/NeuralNetwork/TanhActivation.hs b/src/MachineLearning/NeuralNetwork/TanhActivation.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/TanhActivation.hs
@@ -0,0 +1,32 @@
+{-|
+Module: MachineLearning.NeuralNetwork.TanhActivation.
+Description: Tanh Activation.
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Tanh Activation.
+-}
+
+module MachineLearning.NeuralNetwork.TanhActivation
+(
+  tanh
+  , gradient
+)
+
+where
+
+
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Types (Matrix)
+
+
+tanhGradient :: Matrix -> Matrix
+tanhGradient x = 1 - tanhx*tanhx
+  where tanhx = tanh x
+
+
+gradient :: Matrix -> Matrix -> Matrix
+gradient x dx = dx * (tanhGradient x)
+
diff --git a/src/MachineLearning/NeuralNetwork/Topology.hs b/src/MachineLearning/NeuralNetwork/Topology.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/Topology.hs
@@ -0,0 +1,151 @@
+{-|
+Module: MachineLearning.NeuralNetwork.Topology
+Description: Neural Network's Topology
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Neural Network's Topology
+-}
+
+module MachineLearning.NeuralNetwork.Topology
+(
+  Topology
+  , LossFunc(..)
+  , makeTopology
+  , loss
+  , propagateForward
+  , propagateBackward
+  , numberOutputs
+  , initializeTheta
+  , initializeThetaIO
+  , initializeThetaM
+  , flatten
+  , unflatten
+)
+
+where
+
+import Control.Monad (zipWithM)
+import Data.List (foldl')
+import qualified Control.Monad.Random as RndM
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Utils (listOfTuplesToList)
+import MachineLearning.NeuralNetwork.Layer (Layer(..), Cache(..))
+import MachineLearning.NeuralNetwork.Regularization (Regularization, forwardReg, backwardReg)
+
+
+-- | Loss function's type.
+-- Takes x, weights and y.
+type LossFunc = Matrix -> Matrix -> R
+
+
+-- | Neural network topology has at least 2 elements: numver of input and number of outputs.
+-- And sizes of hidden layers between 2 elements.
+data Topology = Topology [(Int, Int)] [Layer] LossFunc
+
+
+-- | Makes Neural Network's Topology.
+-- Takes number of inputs, list of hidden layers, output layer and loss function.
+makeTopology :: Int -> [Layer] -> Layer -> LossFunc -> Topology
+makeTopology nInputs hiddenLayers outputLayer lossFunc =
+  let layers = hiddenLayers ++ [outputLayer]
+      layerSizes = nInputs : (map lUnits layers)
+      sizes = getThetaSizes layerSizes
+  in Topology sizes layers lossFunc
+      
+
+-- | Calculates loss for the given topology.
+-- Takes topology, regularization, x, weights, y.
+loss :: Topology -> Regularization -> Matrix -> [(Matrix, Matrix)] -> Matrix -> R
+loss (Topology _ _ lf) reg x weights y =
+  let lossValue = lf x y
+      regValue = forwardReg reg weights
+  in lossValue + regValue
+
+
+-- | Implementation of forward propagation algorithm.
+propagateForward :: Topology -> Matrix -> [(Matrix, Matrix)] -> (Matrix, [Cache])
+propagateForward (Topology _ layers _) x thetaList = foldl' f (x, []) $ zip thetaList layers
+  where f (a, cs) (theta, hl) =
+          let (a', cache) = forwardPass hl a theta
+          in (a', cache:cs)
+
+
+-- | Makes one forward step for the given layer.
+forwardPass :: Layer -> Matrix -> (Matrix, Matrix) -> (Matrix, Cache)
+forwardPass layer a (b, w) = (a', Cache z a w)
+  where z = lForward layer a b w
+        a' = lActivation layer z
+
+
+-- | Implementation of backward propagation algorithm.
+propagateBackward :: Topology -> Regularization -> Matrix -> [Cache] -> Matrix -> [(Matrix, Matrix)]
+propagateBackward (Topology _ layers _) reg scores (cache:cacheList) y = gradientList
+  where cache' = Cache scores (cacheX cache) (cacheW cache)
+        cacheList' = cache':cacheList
+        gradientList = snd $ foldl' f (y, []) $ zip cacheList' $ reverse layers
+        f (da, grads) (cache, hl) =
+          let (da', db, dw) = backwardPass hl reg da cache
+          in (da', (db, dw):grads)
+
+
+-- | Makes one backward step for the given layer.
+backwardPass :: Layer -> Regularization -> Matrix -> Cache -> (Matrix, Matrix, Matrix)
+backwardPass layer reg da cache = (da', db, dw')
+  where delta = lActivationGradient layer (cacheZ cache) da
+        (da', db, dw) = lBackward layer delta cache
+        dw' = dw + (backwardReg reg (cacheW cache))
+
+
+-- | Returns number of outputs of the given topology.
+numberOutputs :: Topology -> Int
+numberOutputs (Topology nnt _ _) = fst $ last nnt
+
+
+-- | Returns dimensions of weight matrices for given neural network topology
+getThetaSizes :: [Int] -> [(Int, Int)]
+getThetaSizes nn = zipWith (\r c -> (r, c)) (tail nn) nn
+
+
+-- | Create and initialize weights vector with random values
+-- for given neural network topology.
+-- Takes a seed to initialize generator of random numbers as a first parameter.
+initializeTheta :: Int -> Topology -> Vector
+initializeTheta seed topology = RndM.evalRand (initializeThetaM topology) gen
+  where gen = RndM.mkStdGen seed
+
+
+-- | Create and initialize weights vector with random values
+-- for given neural network topology inside IO Monad.
+initializeThetaIO :: Topology -> IO Vector
+initializeThetaIO = RndM.evalRandIO . initializeThetaM
+
+
+-- | Create and initialize weights vector with random values
+-- for given neural network topology inside RandomMonad.
+initializeThetaM :: RndM.RandomGen g => Topology -> RndM.Rand g Vector
+initializeThetaM topology = flatten <$> initializeThetaListM topology
+
+
+-- | Create and initialize list of weights matrices with random values
+-- for given neural network topology.
+initializeThetaListM :: RndM.RandomGen g => Topology -> RndM.Rand g [(Matrix, Matrix)]
+initializeThetaListM (Topology sizes layers _) = zipWithM lInitializeThetaM layers sizes
+
+
+-- | Flatten list of matrices into vector.
+flatten :: [(Matrix, Matrix)] -> Vector
+flatten ms = V.concat $ map LA.flatten $ listOfTuplesToList ms
+
+
+-- | Unflatten vector into list of matrices for given neural network topology.
+unflatten :: Topology -> Vector -> [(Matrix, Matrix)]
+unflatten (Topology sizes _ _) v =
+  let offsets = reverse $ foldl' (\os (r, c) -> (r+r*c + head os):os) [0] (init sizes)
+      ms = zipWith (\o (r, c) -> (LA.reshape r (slice o r), LA.reshape c (slice (o+r) (r*c)))) offsets sizes
+      slice o n = V.slice o n v
+  in ms
diff --git a/src/MachineLearning/NeuralNetwork/TopologyMaker.hs b/src/MachineLearning/NeuralNetwork/TopologyMaker.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/TopologyMaker.hs
@@ -0,0 +1,87 @@
+{-|
+Module: MachineLearning.NeuralNetwork.TopologyMaker
+Description: Topology Maker
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Topology Maker.
+-}
+
+module MachineLearning.NeuralNetwork.TopologyMaker
+(
+  Activation(..)
+  , Loss(..)
+  , makeTopology
+)
+
+where
+
+
+import qualified MachineLearning.NeuralNetwork.Topology as T
+import MachineLearning.NeuralNetwork.Layer (Layer(..), affineForward, affineBackward)
+import MachineLearning.NeuralNetwork.WeightInitialization (nguyen)
+import qualified MachineLearning.NeuralNetwork.ReluActivation as Relu
+import qualified MachineLearning.NeuralNetwork.TanhActivation as Tanh
+import qualified MachineLearning.NeuralNetwork.SigmoidActivation as Sigmoid
+import qualified MachineLearning.NeuralNetwork.SoftmaxLoss as Softmax
+import qualified MachineLearning.NeuralNetwork.MultiSvmLoss as MultiSvm
+import qualified MachineLearning.NeuralNetwork.LogisticLoss as Logistic
+
+
+data Activation = ASigmoid | ARelu | ATanh
+
+data Loss = LLogistic | LSoftmax | LMultiSvm
+
+
+-- | Creates toplogy. Takes number of inputs, number of outputs and list of hidden layers.
+makeTopology :: Activation -> Loss -> Int -> Int -> [Int] -> T.Topology
+makeTopology a l nInputs nOutputs hlUnits = T.makeTopology nInputs hiddenLayers outputLayer (loss l)
+  where hiddenLayers = map (mkAffineLayer a) hlUnits
+        outputLayer = mkOutputLayer l nOutputs
+
+
+mkAffineLayer a nUnits = Layer {
+  lUnits = nUnits
+  , lForward = affineForward
+  , lActivation = hiddenActivation a
+  , lBackward = affineBackward
+  , lActivationGradient = hiddenGradient a
+  , lInitializeThetaM = nguyen
+  }
+
+
+mkOutputLayer l nUnits = Layer {
+  lUnits = nUnits
+  , lForward = affineForward
+  , lActivation = outputActivation l
+  , lBackward = affineBackward
+  , lActivationGradient = outputGradient l
+  , lInitializeThetaM = nguyen
+  }
+
+
+hiddenActivation ASigmoid = Sigmoid.sigmoid
+hiddenActivation ARelu = Relu.relu
+hiddenActivation ATanh = Tanh.tanh
+
+
+hiddenGradient ASigmoid = Sigmoid.gradient
+hiddenGradient ARelu = Relu.gradient
+hiddenGradient ATanh = Tanh.gradient
+
+
+outputActivation LLogistic = Logistic.scores
+outputActivation LSoftmax = Softmax.scores
+outputActivation LMultiSvm = MultiSvm.scores
+
+
+outputGradient LLogistic = Logistic.gradient
+outputGradient LSoftmax = Softmax.gradient
+outputGradient LMultiSvm = MultiSvm.gradient
+
+
+loss LLogistic = Logistic.loss
+loss LSoftmax = Softmax.loss
+loss LMultiSvm = MultiSvm.loss
diff --git a/src/MachineLearning/NeuralNetwork/WeightInitialization.hs b/src/MachineLearning/NeuralNetwork/WeightInitialization.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/NeuralNetwork/WeightInitialization.hs
@@ -0,0 +1,47 @@
+{-|
+Module: MachineLearning.NeuralNetwork.WeightInitialization
+Description: Weight Initialization
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Various Weight Initialization algorithms.
+-}
+
+module MachineLearning.NeuralNetwork.WeightInitialization
+(
+  nguyen
+  , he
+)
+
+where
+
+
+import qualified Numeric.LinearAlgebra as LA
+import qualified Control.Monad.Random as RndM
+import MachineLearning.Types (Matrix)
+import MachineLearning.Random (getRandomRMatrixM)
+
+
+
+-- | Weight Initialization Algorithm discussed in Nguyen at al.: https://web.stanford.edu/class/ee373b/nninitialization.pdf
+-- Nguyen, Derrick, Widrow, B. Improving the learning speed of 2-layer neural networks by choosing initial values of adaptive weights.
+-- In Proc. IJCNN, 1990; 3: 21-26.
+nguyen :: RndM.RandomGen g => (Int, Int) -> RndM.Rand g (Matrix, Matrix)
+nguyen (r, c) = do
+  let b = LA.konst 0 (1, r)
+      eps = (sqrt 6) / (sqrt . fromIntegral $ r + c)
+  w <- getRandomRMatrixM r c (-eps, eps)
+  return (b, w)
+
+
+-- | Weight Initialization Algorithm discussed in He at al.: https://arxiv.org/abs/1502.01852
+-- Kaiming He, Xiangyu Zhang, Shaoqing Ren, Jian Sun.
+-- Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification.
+he :: RndM.RandomGen g => (Int, Int) -> RndM.Rand g (Matrix, Matrix)
+he (r, c) = do
+  let b = LA.konst 0 (1, r)
+      eps = sqrt (2/(fromIntegral $ r + c))
+  w <- getRandomRMatrixM r c (-eps, eps)
+  return (b, w)
diff --git a/src/MachineLearning/Optimization.hs b/src/MachineLearning/Optimization.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Optimization.hs
@@ -0,0 +1,84 @@
+{-|
+Module: MachineLearning.Optimization
+Description: Optimization
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Optimization module.
+-}
+
+module MachineLearning.Optimization
+(
+    MinimizeMethod(..)
+  , minimize
+  , checkGradient
+)
+
+where
+
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Model as Model
+import MachineLearning.Regularization (Regularization)
+import qualified MachineLearning.Optimization.GradientDescent as GD
+import qualified MachineLearning.Optimization.MinibatchGradientDescent as MGD
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+
+import qualified Numeric.GSL.Minimization as Min
+
+data MinimizeMethod = GradientDescent R       -- ^ Gradient descent, takes alpha. Requires feature normalization.
+                    | MinibatchGradientDescent Int Int R  -- ^ Minibacth Gradietn Descent, takes seed, batch size and alpha
+                    | ConjugateGradientFR R R -- ^ Fletcher-Reeves conjugate gradient algorithm,
+                                              -- takes size of first trial step (0.1 is fine) and tol (0.1 is fine).
+                    | ConjugateGradientPR R R -- ^ Polak-Ribiere conjugate gradient algorithm.
+                                              -- takes size of first trial step (0.1 is fine) and tol (0.1 is fine).
+                    | BFGS2 R R               -- ^ Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm,
+                                              -- takes size of first trial step (0.1 is fine) and tol (0.1 is fine).
+
+
+-- | Returns solution vector (theta) and optimization path.
+-- Optimization path's row format:
+-- [iter number, cost function value, theta values...]
+minimize :: Model.Model a =>
+            MinimizeMethod
+         -> a       -- ^ model (Least Squares, Logistic Regression etc)
+         -> R   -- ^ epsilon, desired precision of the solution
+         -> Int     -- ^ maximum number of iterations allowed
+         -> Regularization   -- ^ regularization parameter
+         -> Matrix  -- ^ X
+         -> Vector  -- ^ y
+         -> Vector  -- ^ initial solution, theta
+         -> (Vector, Matrix) -- ^ solution vector and optimization path
+
+minimize (BFGS2 firstStepSize tol) = minimizeVD Min.VectorBFGS2 firstStepSize tol
+minimize (ConjugateGradientFR firstStepSize tol) = minimizeVD Min.ConjugateFR firstStepSize tol
+minimize (ConjugateGradientPR firstStepSize tol) = minimizeVD Min.ConjugatePR firstStepSize tol
+minimize (GradientDescent alpha) = GD.gradientDescent alpha
+minimize (MinibatchGradientDescent seed batchSize alpha) = MGD.minibatchGradientDescent seed batchSize alpha
+
+
+minimizeVD method firstStepSize tol model epsilon niters reg x y initialTheta
+  = Min.minimizeVD method epsilon niters firstStepSize tol costf gradientf initialTheta
+  where costf = Model.cost model reg x y
+        gradientf = Model.gradient model reg x y
+
+-- | Gradient checking function.
+-- Approximates the derivates of the Model's cost function
+-- and calculates derivatives using the Model's gradient functions.
+-- Returns norn_2 between 2 derivatives.
+-- Takes model, regularization, X, y, theta and epsilon (used to approximate derivatives, 1e-4 is a good value).
+checkGradient :: Model a => a -> Regularization -> Matrix -> Vector -> Vector -> R -> R
+checkGradient model reg x y theta eps = LA.norm_2 $ grad1 - grad2
+  where theta_m = LA.asColumn theta
+        eps_v = V.replicate (V.length theta) eps
+        eps_m = LA.diag eps_v
+        theta_m1 = theta_m + eps_m
+        theta_m2 = theta_m - eps_m
+        cost1 = LA.vector $ map (cost model reg x y) $ LA.toColumns theta_m1
+        cost2 = LA.vector $ map (cost model reg x y) $ LA.toColumns theta_m2
+        grad1 = (cost1 - cost2) / (LA.scalar $ 2*eps)
+        grad2 = gradient model reg x y theta
+
diff --git a/src/MachineLearning/Optimization/GradientDescent.hs b/src/MachineLearning/Optimization/GradientDescent.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Optimization/GradientDescent.hs
@@ -0,0 +1,38 @@
+{-|
+Module: MachineLearning.Optimization.GradientDescent
+Description: Gradient Descent
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+-}
+
+module MachineLearning.Optimization.GradientDescent
+(
+  gradientDescent
+)
+
+where
+
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Regularization (Regularization)
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+import qualified MachineLearning.Model as Model
+
+-- | Gradient Descent method implementation. See "MachineLearning.Regression" for usage details.
+gradientDescent :: Model.Model a => R-> a -> R -> Int -> Regularization -> Matrix -> Vector -> Vector -> (Vector, Matrix)
+gradientDescent alpha model eps maxIters lambda x y theta = helper theta maxIters []
+  where gradient = Model.gradient model lambda
+        cost = Model.cost model lambda
+        helper theta nIters optPath =
+          let theta' = theta - (LA.scale alpha (gradient x y theta))
+              j = cost x y theta'
+              gradientTest = LA.norm_2 (theta' - theta) < eps
+              optPathRow = V.concat [LA.vector [(fromIntegral $ maxIters - nIters), j], theta']
+              optPath' = optPathRow : optPath
+          in if gradientTest || nIters <= 1
+             then (theta, LA.fromRows $ reverse optPath')
+             else helper theta' (nIters - 1) optPath'
diff --git a/src/MachineLearning/Optimization/MinibatchGradientDescent.hs b/src/MachineLearning/Optimization/MinibatchGradientDescent.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Optimization/MinibatchGradientDescent.hs
@@ -0,0 +1,75 @@
+{-|
+Module: MachineLearning.Optimization.MinibatchGradientDescent
+Description: Gradient Descent
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Minibatch Gradient Descent
+-}
+
+module MachineLearning.Optimization.MinibatchGradientDescent
+(
+  minibatchGradientDescent
+)
+
+where
+
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Regularization (Regularization)
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra ((?))
+import qualified Control.Monad.Random as RndM
+
+import qualified MachineLearning.Model as Model
+
+-- | Minibatch Gradient Descent method implementation. See "MachineLearning.Regression" for usage details.
+minibatchGradientDescent :: Model.Model a
+                            => Int              -- ^ seed
+                            -> Int              -- ^ batch size
+                            -> R                -- ^ learning rate, alpha
+                            -> a                -- ^ model to learn
+                            -> R                -- ^ epsilon
+                            -> Int              -- ^ max number of iters
+                            -> Regularization   -- ^ regularization parameter, lambda
+                            -> Matrix           -- ^ matrix of features, X
+                            -> Vector           -- ^ output vector, y
+                            -> Vector           -- ^ vector of initial weights, theta or w
+                            -> (Vector, Matrix) -- ^ vector of weights and learning path 
+minibatchGradientDescent seed batchSize alpha model eps maxIters lambda x y theta =
+  RndM.evalRand (minibatchGradientDescentM batchSize alpha model eps maxIters lambda x y theta) gen
+  where gen = RndM.mkStdGen seed
+
+
+-- | Minibatch Gradient Descent method implementation. See "MachineLearning.Regression" for usage details.
+minibatchGradientDescentM :: (Model.Model a, RndM.RandomGen g)
+                             => Int              -- ^ batch size
+                             -> R                -- ^ learning rate, alpha
+                             -> a                -- ^ model to learn
+                             -> R                -- ^ epsilon
+                             -> Int              -- ^ max number of iters
+                             -> Regularization   -- ^ regularization parameter, lambda
+                             -> Matrix           -- ^ matrix of features, X
+                             -> Vector           -- ^ output vector, y
+                             -> Vector           -- ^ vector of initial weights, theta or w
+                             -> RndM.Rand g (Vector, Matrix) -- ^ vector of weights and learning path 
+minibatchGradientDescentM batchSize alpha model eps maxIters lambda x y theta = do
+  idxList <- RndM.getRandomRs (0, (LA.rows x) - 1)
+  let gradient = Model.gradient model lambda
+      cost = Model.cost model lambda
+      helper theta nIters optPath =
+          let idx = take batchSize idxList
+              x' = x ? idx
+              y' = LA.flatten $ (LA.asColumn y) ? idx
+              theta' = theta - (LA.scale alpha (gradient x' y' theta))
+              j = cost x' y' theta'
+              gradientTest = LA.norm_2 (theta' - theta) < eps
+              optPathRow = V.concat [LA.vector [(fromIntegral $ maxIters - nIters), j], theta']
+              optPath' = optPathRow : optPath
+          in if gradientTest || nIters <= 1
+             then (theta, LA.fromRows $ reverse optPath')
+             else helper theta' (nIters - 1) optPath'
+  return $ helper theta maxIters []
+
diff --git a/src/MachineLearning/PCA.hs b/src/MachineLearning/PCA.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/PCA.hs
@@ -0,0 +1,65 @@
+{-|
+Module: MachineLearning.PCA
+Description: Principal Component Analysis.
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Principal Component Analysis (PCA) - dimensionality reduction algorithm.
+It is mostly used to speed up supervising learning (Regression, Classification, etc) and visualization of data.
+-}
+
+module MachineLearning.PCA
+(
+  getDimReducer
+  , getDimReducer_rv
+)
+
+where
+
+import Data.Maybe (fromMaybe)
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra ((<>), (??))
+import qualified MachineLearning as ML
+
+import MachineLearning.Types (R, Vector, Matrix)
+
+
+-- | Computes "covariance matrix".
+covarianceMatrix :: Matrix -> Matrix
+covarianceMatrix x = ((LA.tr x) <> x) / (fromIntegral $ LA.rows x)
+
+
+-- | Compute eigenvectors (matrix U) and singular values (matrix S) of the given covariance matrix.
+pca :: Matrix -> (Matrix, Vector)
+pca x = (u, s)
+  where sigma = covarianceMatrix x
+        (u, s, _) = LA.svd sigma
+
+
+-- | Gets dimensionality reduction function, retained variance (0..1) and reduced X
+-- for given matrix X and number of dimensions to retain.
+getDimReducer :: Matrix -> Int -> (Matrix -> Matrix, R, Matrix)
+getDimReducer x k = (reducer, retainedVariance, reducer xNorm)
+  where muSigma = ML.meanStddev x
+        xNorm = ML.featureNormalization muSigma x
+        (u, s) = pca xNorm
+        u' = u ?? (LA.All, LA.Take k)
+        reducer xx = (ML.featureNormalization muSigma xx) <> u'
+        retainedVariance = (V.sum $ V.slice 0 k s) / (V.sum s)
+
+
+-- | Gets dimensionality reduction function, retained number of dimensions and reduced X
+-- for given matrix X and variance to retain (0..1].
+getDimReducer_rv :: Matrix -> R -> (Matrix -> Matrix, Int, Matrix)
+getDimReducer_rv x rv = (reducer, k, reducer xNorm)
+  where muSigma = ML.meanStddev x
+        xNorm = ML.featureNormalization muSigma x
+        (u, s) = pca xNorm
+        sum_s = V.sum s
+        variances = (V.scanl' (+) 0 s) / (LA.scalar sum_s)
+        k = fromMaybe ((V.length s) - 1) $ V.findIndex (\v -> v >= rv) variances
+        u' = u ?? (LA.All, LA.Take k)
+        reducer xx = (ML.featureNormalization muSigma xx) <> u'
diff --git a/src/MachineLearning/Random.hs b/src/MachineLearning/Random.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Random.hs
@@ -0,0 +1,80 @@
+{-|
+Module: MachineLearning.Random
+Description: Random generation utility functions.
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Randon generation uitility functions.
+-}
+
+module MachineLearning.Random
+(
+  sample
+  , sampleM
+  , getRandomRListM
+  , getRandomRVectorM
+  , getRandomRMatrixM
+)
+
+where
+
+import Control.Monad (when)
+import System.Random (RandomGen, Random)
+import qualified Control.Monad.ST as ST
+import qualified Data.Vector as V
+import qualified Data.Vector.Storable as SV
+import Numeric.LinearAlgebra ((><))
+import qualified Data.Vector.Mutable as MV
+import qualified Control.Monad.Random as RndM
+import MachineLearning.Types (R, Vector, Matrix)
+
+
+-- | Samples `n` (given as a second parameter) values from `list` (given as a third parameter).
+sample :: RandomGen g => g -> Int -> V.Vector a -> (V.Vector a, g)
+sample gen n xs = RndM.runRand (sampleM n xs) gen
+
+
+-- | Samples `n` (given as a second parameter) values from `list` (given as a third parameter) inside RandomMonad.
+sampleM :: RandomGen g => Int -> V.Vector a -> RndM.Rand g (V.Vector a)
+sampleM n xs = do  -- Random Monad starts
+  let rangeList = V.fromList $ zip (repeat 0) [n..(length xs)-1]
+  rnds <- randomsInRangesM rangeList
+  let (pre, post) = V.splitAt n xs
+  let ys = ST.runST $ do  -- ST Monad starts
+        mv <- V.thaw pre
+        V.zipWithM_ (\val r -> when (r < n) $ MV.write mv (mod r n) val) post rnds
+        V.unsafeFreeze mv
+  return ys
+
+
+-- | Returns a list of random values distributed in a closed interval `range`
+getRandomRListM :: (RandomGen g, Random a) =>
+                   Int             -- ^ list's lengths
+                   -> (a, a)       -- ^ range
+                   -> RndM.Rand g [a]          -- ^ list of random values inside RandomMonad
+getRandomRListM size range = mapM (\_ -> RndM.getRandomR range) [1..size]
+
+
+-- | Returns a vector of random values distributed in a closed interval `range`
+getRandomRVectorM :: RandomGen g =>
+                    Int                     -- ^ vector's length
+                    -> (R, R)               -- ^ range
+                    -> RndM.Rand g Vector   -- vector of randon values inside RandomMonad
+getRandomRVectorM size range = SV.fromList <$> getRandomRListM size range
+
+
+-- | Returns a matrix of random values distributed in a closed interval `range`
+getRandomRMatrixM :: RandomGen g =>
+                    Int                     -- ^ number of rows
+                    -> Int                  -- ^ number of columns
+                    -> (R, R)               -- ^ range
+                    -> RndM.Rand g Matrix   -- vector of randon values inside RandomMonad
+getRandomRMatrixM r c range = (r><c) <$> getRandomRListM (r*c) range
+
+
+-- | Takes a list of ranges `(lo, hi)`,
+-- returns a list of random values uniformly distributed in the list of closed intervals [(lo, hi)].
+randomsInRangesM :: (RndM.RandomGen g, RndM.Random a) => V.Vector (a, a) -> RndM.Rand g (V.Vector a)
+randomsInRangesM rangeList = mapM RndM.getRandomR rangeList
diff --git a/src/MachineLearning/Regression.hs b/src/MachineLearning/Regression.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Regression.hs
@@ -0,0 +1,46 @@
+{-|
+Module: MachineLearning.Regression
+Description: Regression
+Copyright: (c) Alexander Ignatyev, 2016-2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+-}
+
+module MachineLearning.Regression
+(
+  Model.Model(..)
+  , LeastSquares.LeastSquaresModel(..)
+  , Optimization.MinimizeMethod(..)
+  , Optimization.minimize
+  , normalEquation
+  , normalEquation_p
+  , Regularization(..)
+)
+
+where
+
+import MachineLearning.Types (Vector, Matrix)
+import MachineLearning.Optimization as Optimization
+import MachineLearning.Model as Model
+import MachineLearning.LeastSquaresModel as LeastSquares
+import MachineLearning.Regularization (Regularization(..))
+
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra ((<>), (#>))
+
+
+-- | Normal equation using inverse, does not require feature normalization
+-- It takes X and y, returns theta.
+normalEquation :: Matrix -> Vector -> Vector
+normalEquation x y =
+  let trX = LA.tr x
+  in (LA.inv (trX <> x) <> trX) #> y
+
+
+-- | Normal equation using pseudo inverse, requires feature normalization
+-- It takes X and y, returns theta.
+normalEquation_p :: Matrix -> Vector -> Vector
+normalEquation_p x y =
+  let trX = LA.tr x
+  in (LA.pinv (trX <> x) <> trX) #> y
diff --git a/src/MachineLearning/Regularization.hs b/src/MachineLearning/Regularization.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Regularization.hs
@@ -0,0 +1,45 @@
+{-|
+Module: MachineLearning.Regularization
+Description: Regularization
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Regularization.
+-}
+
+module MachineLearning.Regularization
+(
+  Regularization(..)
+  , costReg
+  , gradientReg
+)
+
+where
+
+import MachineLearning.Types (R, Vector)
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+-- | Regularization
+data Regularization = RegNone -- ^ No regularization
+                    | L2 R    -- ^ L2 Regularization
+
+
+
+-- | Calculates regularization for Model.cost function.
+-- It takes regularization parameter and theta.
+costReg :: Regularization -> Vector -> R
+costReg RegNone _ = 0
+costReg (L2 lambda) theta = (thetaReg LA.<.> thetaReg) * lambda * 0.5
+  where thetaReg = V.tail theta
+
+
+
+-- | Calculates regularization for Model.gradient function.
+-- It takes regularization parameter and theta.
+gradientReg :: Regularization -> Vector -> Vector
+gradientReg RegNone _ = 0
+gradientReg (L2 lambda) theta = (LA.scalar lambda) * thetaReg
+  where thetaReg = theta V.// [(0, 0)]
diff --git a/src/MachineLearning/SoftmaxClassifier.hs b/src/MachineLearning/SoftmaxClassifier.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/SoftmaxClassifier.hs
@@ -0,0 +1,70 @@
+{-|
+Module: MachineLearning.SoftmaxClassifier
+Description: Softmax Classifier.
+Copyright: (c) Alexander Ignatyev, 2017.
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Softmax Classifier (Multiclass Logistic Regression).
+-}
+
+module MachineLearning.SoftmaxClassifier
+(
+  module MachineLearning.Model
+  , module MachineLearning.Classification.MultiClass
+  , SoftmaxClassifier(..)
+)
+
+where
+
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra((<>), (<.>), (|||))
+import qualified Data.Vector.Storable as V
+
+import qualified MachineLearning as ML
+import MachineLearning.Types (R, Vector, Matrix)
+import MachineLearning.Utils (reduceByRows, sumByRows)
+import MachineLearning.Model
+import MachineLearning.Classification.MultiClass
+
+
+-- | Softmax Classifier, takes number of classes.
+data SoftmaxClassifier = Softmax Int
+
+instance Classifier SoftmaxClassifier where
+  cscore (Softmax _) x theta = scores - reduceByRows V.maximum scores
+    where scores = x <> (LA.tr theta)
+
+  chypothesis m x theta = V.fromList predictions
+    where scores = cscore m x theta
+          scores' = LA.toRows scores
+          predictions = map (fromIntegral . LA.maxIndex) scores'
+
+  ccost m lambda x y theta =
+    let nSamples = fromIntegral $ LA.rows x
+        scores = cscore m x theta
+        sum_probs = sumByRows $ exp scores
+        loss = LA.sumElements $ (log sum_probs) - remap scores y
+        reg = ccostReg lambda theta
+    in (loss + reg) / nSamples 
+
+  cgradient m lambda x y theta =
+    let nSamples = fromIntegral $ LA.rows x
+        ys = processOutput m y
+        scores = cscore m x theta
+        sum_probs = sumByRows $ exp scores
+        probs = (exp scores) / sum_probs
+        probs' = probs - ys
+        dw =  (LA.tr probs') <> x
+        reg = cgradientReg lambda theta
+    in (dw + reg)/ nSamples
+
+  cnumClasses (Softmax nLabels) = nLabels
+
+
+remap :: Matrix -> Vector -> Matrix
+remap m v = LA.remap cols rows m
+  where cols = LA.asColumn $ V.fromList [0..(fromIntegral $ LA.rows m)-1]
+        rows = LA.toInt $ LA.asColumn v
+
diff --git a/src/MachineLearning/TerminalProgress.hs b/src/MachineLearning/TerminalProgress.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/TerminalProgress.hs
@@ -0,0 +1,68 @@
+{-|
+Module: TerminalProgress
+Description: Learn function with progress bar for terminal.
+Copyright: (c) Alexander Ignatyev, 2017
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Learn function with progress bar for terminal.
+-}
+
+module MachineLearning.TerminalProgress
+(
+  learnWithProgressBar
+  , learnOneVsAllWithProgressBar
+)
+
+where
+
+import Data.List (transpose)
+import MachineLearning.Types (Vector, Matrix)
+import qualified MachineLearning.Classification.Internal as MLC
+import Control.Monad (foldM, mapAndUnzipM)
+import Control.DeepSeq (deepseq)
+import qualified System.Console.AsciiProgress as AP
+import qualified Numeric.LinearAlgebra as LA
+
+
+-- | Learn the given function displaying progress bar in terminal.
+-- It takes function, initial theta and number of iterations to call the function.
+-- It returns theta and optimization path (see "MachineLearning.Optimization" for details).
+learnWithProgressBar :: (Vector -> (Vector, Matrix)) -> Vector -> Int -> IO (Vector, Matrix)
+learnWithProgressBar func initialTheta nIterations = AP.displayConsoleRegions $ do
+  pg <- newProgressBar nIterations
+  (theta, optPaths) <- foldM (doLoop pg func) (initialTheta, []) [1..nIterations]
+  return (theta, buildOptPathMatrix $ reverse optPaths)
+
+
+-- | Learn the given function displaying progress bar in terminal.
+-- It takes function, list of outputs and list of initial thetas and number of iterations to call the function.
+-- It returns list of thetas and list of optimization paths (see "MachineLearning.Optimization" for details).
+learnOneVsAllWithProgressBar :: (Vector -> Vector -> (Vector, Matrix)) -> Vector -> [Vector] -> Int -> IO ([Vector], [Matrix])
+learnOneVsAllWithProgressBar func y initialThetaList nIterations = AP.displayConsoleRegions $ do
+    let numLabels = length initialThetaList
+        ys = MLC.processOutputOneVsAll numLabels y
+    pg <- newProgressBar $ nIterations * (length ys)
+    mapAndUnzipM (learnOneClass pg func nIterations) $ zip ys initialThetaList
+
+
+newProgressBar nIterations = AP.newProgressBar AP.def {
+  AP.pgTotal = fromIntegral nIterations
+  , AP.pgFormat = "Learning :percent [:bar] (for :elapsed, :eta remaining)"
+  }
+
+doLoop pg func (theta, optPaths) _ = do
+  let (theta', optPath) = func theta
+  theta' `deepseq` AP.tick pg
+  return (theta', (optPath : optPaths))
+
+
+learnOneClass pg func nIterations (y, theta) = do
+  (theta, optPaths) <- foldM (doLoop pg $ func y) (theta, []) [1..nIterations]
+  return (theta, buildOptPathMatrix $ reverse optPaths)
+
+
+-- | Build a single optimazation path matrix from list of optimization path matrices.
+buildOptPathMatrix :: [Matrix] -> Matrix
+buildOptPathMatrix matrices = LA.fromBlocks $ map (\m -> [m]) matrices
diff --git a/src/MachineLearning/Types.hs b/src/MachineLearning/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Types.hs
@@ -0,0 +1,31 @@
+{-|
+Module: MachineLearning.Types
+Description: Common Types
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Common type definitions used in all modules.
+-}
+
+module MachineLearning.Types
+(
+    R
+  , Vector
+  , Matrix
+)
+
+where
+
+import qualified Numeric.LinearAlgebra.Data as LAD
+
+-- | Scalar Type (Double)
+type R = LAD.R
+
+-- | Vector Types (of Doubles)
+type Vector = LAD.Vector R
+
+-- | Matrix Types (of Doubles)
+type Matrix = LAD.Matrix R
+
diff --git a/src/MachineLearning/Utils.hs b/src/MachineLearning/Utils.hs
new file mode 100644
--- /dev/null
+++ b/src/MachineLearning/Utils.hs
@@ -0,0 +1,58 @@
+{-|
+Module: MachineLearning.Utils
+Description: Utils
+Copyright: (c) Alexander Ignatyev, 2016
+License: BSD-3
+Stability: experimental
+Portability: POSIX
+
+Various helpful utilities.
+-}
+
+module MachineLearning.Utils
+(
+  reduceByRowsV
+  , reduceByColumnsV
+  , reduceByRows
+  , reduceByColumns
+  , sumByRows
+  , sumByColumns
+  , listOfTuplesToList
+)
+
+where
+
+  
+import MachineLearning.Types (R, Vector, Matrix)
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+
+reduceByRowsV :: (Vector -> R) -> Matrix -> Vector
+reduceByRowsV f = LA.vector . map f . LA.toRows
+
+
+reduceByColumnsV :: (Vector -> R) -> Matrix -> Vector
+reduceByColumnsV f = LA.vector . map f . LA.toColumns
+
+
+reduceByRows :: (Vector -> R) -> Matrix -> Matrix
+reduceByRows f = LA.asColumn . reduceByRowsV f
+
+
+reduceByColumns :: (Vector -> R) -> Matrix -> Matrix
+reduceByColumns f = LA.asRow . reduceByColumnsV f
+
+
+sumByColumns :: Matrix -> Matrix
+sumByColumns = reduceByColumns V.sum
+
+
+sumByRows :: Matrix -> Matrix
+sumByRows = reduceByRows V.sum
+
+
+-- | Converts list of tuples into list.
+listOfTuplesToList :: [(a, a)] -> [a]
+listOfTuplesToList [] = []
+listOfTuplesToList ((a, b):xs) = a : b : listOfTuplesToList xs
diff --git a/test/MachineLearning/Classification/BinaryTest.hs b/test/MachineLearning/Classification/BinaryTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/Classification/BinaryTest.hs
@@ -0,0 +1,80 @@
+module MachineLearning.Classification.BinaryTest
+(
+  tests
+  , testOptPath
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.Types
+import qualified Data.Vector as V
+import qualified Numeric.LinearAlgebra as LA
+
+import MachineLearning.DataSets (dataset2)
+
+import qualified MachineLearning as ML
+import MachineLearning.Classification.Binary
+
+(x, y) = ML.splitToXY dataset2
+
+
+processX muSigma x = ML.addBiasDimension $ ML.featureNormalization muSigma $ ML.mapFeatures 6 x
+
+muSigma = ML.meanStddev (ML.mapFeatures 6 x)
+x1 = processX muSigma x
+zeroTheta = LA.konst 0 (LA.cols x1)
+
+xPredict = LA.matrix 2 [ -0.5, 0.5
+                       , 0.2, -0.2
+                       , 1, 1
+                       , 1, 0
+                       , 0, 0
+                       , 0, 1]
+xPredict1 = processX muSigma xPredict
+yExpected = LA.vector [1, 1, 0, 0, 1, 0]
+
+eps = 0.0001
+
+
+-- Binary
+
+(thetaCGFR, optPathCGFR) = learn (ConjugateGradientFR 0.1 0.1) eps 50 (L2 0.5) x1 y zeroTheta
+(thetaCGPR, optPathCGPR) = learn (ConjugateGradientPR 0.1 0.1) eps 50 (L2 0.5) x1 y zeroTheta
+(thetaBFGS, optPathBFGS) = learn (BFGS2 0.1 0.1) eps 50 (L2 0.5) x1 y zeroTheta
+
+
+isInDescendingOrder :: V.Vector Double -> Bool
+isInDescendingOrder lst = V.and . snd . V.unzip $ V.scanl (\(prev, _) current -> (current, prev-current > (-0.001))) (1/0, True) lst
+
+testOptPath optPath = do
+  let js = V.convert $ (LA.toColumns optPath) !! 1
+  assertBool ("non-increasing errors: " ++ show js) $ isInDescendingOrder js
+
+testAccuracyBinary theta eps = do
+  let yPredicted = predict x1 theta
+      accuracy = calcAccuracy y yPredicted
+  assertApproxEqual "" eps 1 accuracy
+
+tests = [
+  testGroup "learn" [
+      testCase "Conjugate Gradient FR" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaCGFR)
+      , testCase "Conjugate Gradient PR" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaCGPR)
+      , testCase "BFGS" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaBFGS)
+      ]
+  , testGroup "optPath" [
+      testCase "Conjugate Gradient FR" $ testOptPath optPathCGFR
+      , testCase "Conjugate Gradient PR" $ testOptPath optPathCGPR
+      , testCase "BFGS" $ testOptPath optPathBFGS
+      ]
+    , testGroup "accuracy" [
+        testCase "Conjugate Gradient FR" $ testAccuracyBinary thetaCGFR 0.2
+        , testCase "Conjugate Gradient PR" $ testAccuracyBinary thetaCGPR 0.2
+        , testCase "BFGS" $ testAccuracyBinary thetaBFGS 0.2
+        ]
+  ]
diff --git a/test/MachineLearning/Classification/OneVsAllTest.hs b/test/MachineLearning/Classification/OneVsAllTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/Classification/OneVsAllTest.hs
@@ -0,0 +1,66 @@
+module MachineLearning.Classification.OneVsAllTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.Types
+import qualified Data.Vector as V
+import qualified Numeric.LinearAlgebra as LA
+
+import MachineLearning.DataSets (dataset2)
+
+import qualified MachineLearning as ML
+import MachineLearning.Classification.OneVsAll
+import MachineLearning.Classification.BinaryTest (testOptPath)
+
+(x, y) = ML.splitToXY dataset2
+
+
+processX muSigma x = ML.addBiasDimension $ ML.featureNormalization muSigma $ ML.mapFeatures 6 x
+
+muSigma = ML.meanStddev (ML.mapFeatures 6 x)
+x1 = processX muSigma x
+zeroTheta :: Vector
+zeroTheta = LA.konst 0 (LA.cols x1) 
+
+xPredict = LA.matrix 2 [ -0.5, 0.5
+                       , 0.2, -0.2
+                       , 1, 1
+                       , 1, 0
+                       , 0, 0
+                       , 0, 1]
+xPredict1 = processX muSigma xPredict
+yExpected = LA.vector [1, 1, 0, 0, 1, 0]
+
+eps = 0.000001
+
+
+zeroThetam = replicate 2 zeroTheta
+(thetaGDm, optPathGDm) = learn (GradientDescent 0.0001) eps 50 (L2 1) 2 x1 y zeroThetam
+(thetaCGFRm, optPathCGFRm) = learn (ConjugateGradientFR 0.1 0.1) eps 50 (L2 0.5) 2 x1 y zeroThetam
+(thetaCGPRm, optPathCGPRm) = learn (ConjugateGradientPR 0.1 0.1) eps 50 (L2 0.5) 2 x1 y zeroThetam
+(thetaBFGSm, optPathBFGSm) = learn (BFGS2 0.1 0.1) eps 50 (L2 0.5) 2 x1 y zeroThetam
+
+
+tests = [
+  testGroup "learn" [
+      testCase "Gradient Descent" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaGDm)
+      , testCase "Conjugate Gradient FR" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaCGFRm)
+      , testCase "Conjugate Gradient PR" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaCGPRm)
+      , testCase "BFGS" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaBFGSm)
+      ]
+  , testGroup "optPath" [
+      testCase "Gradient Descent" $ mapM_ testOptPath optPathGDm
+      , testCase "Conjugate Gradient FR" $ mapM_ testOptPath optPathCGFRm
+      , testCase "Conjugate Gradient PR" $ mapM_ testOptPath optPathCGPRm
+      , testCase "BFGS" $ mapM_ testOptPath optPathBFGSm
+      ]
+  ]
diff --git a/test/MachineLearning/ClusteringTest.hs b/test/MachineLearning/ClusteringTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/ClusteringTest.hs
@@ -0,0 +1,63 @@
+{-# LANGUAGE BangPatterns #-}
+module MachineLearning.ClusteringTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.Types
+import qualified Data.Vector as V
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra ((><))
+import qualified Control.Monad.Random as RndM
+import MachineLearning.Clustering
+
+x1 :: Matrix
+x1 = (7><5) [ 1, 2, 3, 4, 5
+            , 1, 2, 3, 4, 5
+            , 1, 2, 3, 4, 5
+            , 7, 4, 3, 2, 1
+            , 7, 4, 3, 2, 1
+            , 1, 2, 3, 4, 5
+            , 1, 2, 3, 4, 5]
+
+x2 :: Matrix
+x2 = (7><5) [ 1.1, 2, 3, 4, 5
+            , 1, 2, 4, 4, 5
+            , 1, 2, 3, 4, 5
+            , 5, 4, 3, 2, 1
+            , 7, 4, 3, 2, 1
+            , 1, 2, 3, 4, 5
+            , 0.5, 2, 3, 4, 5]
+
+
+testKmeans x k expectedK = do
+  let gen = RndM.mkStdGen 10171
+      clusters = RndM.evalRand (kmeans 10 x k) gen
+  assertEqual "number of clusters" expectedK (V.length clusters)
+
+isInDescendingOrder :: [Double] -> Bool
+isInDescendingOrder lst = and . snd . unzip $ scanl (\(prev, _) current -> (current, prev-current > (-0.001))) (1/0, True) lst
+
+testDescOrderOfCostValues = do
+  let gen = RndM.mkStdGen 10171
+      samples = V.fromList $ LA.toRows x1
+      (clusters, js) = RndM.evalRand (kmeansIterM samples 3 1) gen
+  assertBool "" (isInDescendingOrder js)
+
+tests = [testGroup "kmeans" [
+            testCase "perfect clusters, k = 2" $ testKmeans x1 2 2
+            , testCase "perfect clusters, k = 3" $ testKmeans x1 3 2
+            , testCase "good clusters, k = 2" $ testKmeans x2 2 2
+            , testCase "good clusters, k = 3" $ testKmeans x2 3 3
+            , testCase "good clusters, k = 4" $ testKmeans x2 4 4
+            , testCase "descending order" testDescOrderOfCostValues
+            ]
+        ]
diff --git a/test/MachineLearning/DataSets.hs b/test/MachineLearning/DataSets.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/DataSets.hs
@@ -0,0 +1,182 @@
+module MachineLearning.DataSets
+(
+  dataset1
+  , dataset2
+)
+
+where
+
+import Numeric.LinearAlgebra.Data (matrix)
+
+dataset1 = matrix 3
+  [
+    2104, 3, 399900,
+    1600, 3, 329900,
+    2400, 3, 369000,
+    1416, 2, 232000,
+    3000, 4, 539900,
+    1985, 4, 299900,
+    1534, 3, 314900,
+    1427, 3, 198999,
+    1380, 3, 212000,
+    1494, 3, 242500,
+    1940, 4, 239999,
+    2000, 3, 347000,
+    1890, 3, 329999,
+    4478, 5, 699900,
+    1268, 3, 259900,
+    2300, 4, 449900,
+    1320, 2, 299900,
+    1236, 3, 199900,
+    2609, 4, 499998,
+    3031, 4, 599000,
+    1767, 3, 252900,
+    1888, 2, 255000,
+    1604, 3, 242900,
+    1962, 4, 259900,
+    3890, 3, 573900,
+    1100, 3, 249900,
+    1458, 3, 464500,
+    2526, 3, 469000,
+    2200, 3, 475000,
+    2637, 3, 299900,
+    1839, 2, 349900,
+    1000, 1, 169900,
+    2040, 4, 314900,
+    3137, 3, 579900,
+    1811, 4, 285900,
+    1437, 3, 249900,
+    1239, 3, 229900,
+    2132, 4, 345000,
+    4215, 4, 549000,
+    2162, 4, 287000,
+    1664, 2, 368500,
+    2238, 3, 329900,
+    2567, 4, 314000,
+    1200, 3, 299000,
+    852, 2, 179900,
+    1852, 4, 299900,
+    1203, 3, 239500
+  ]
+
+dataset2 = matrix 3
+  [
+    0.051267, 0.69956, 1,
+    -0.092742, 0.68494, 1,
+    -0.21371, 0.69225, 1,
+    -0.375, 0.50219, 1,
+    -0.51325, 0.46564, 1,
+    -0.52477, 0.2098, 1,
+    -0.39804, 0.034357, 1,
+    -0.30588, -0.19225, 1,
+    0.016705, -0.40424, 1,
+    0.13191, -0.51389, 1,
+    0.38537, -0.56506, 1,
+    0.52938, -0.5212, 1,
+    0.63882, -0.24342, 1,
+    0.73675, -0.18494, 1,
+    0.54666, 0.48757, 1,
+    0.322, 0.5826, 1,
+    0.16647, 0.53874, 1,
+    -0.046659, 0.81652, 1,
+    -0.17339, 0.69956, 1,
+    -0.47869, 0.63377, 1,
+    -0.60541, 0.59722, 1,
+    -0.62846, 0.33406, 1,
+    -0.59389, 0.005117, 1,
+    -0.42108, -0.27266, 1,
+    -0.11578, -0.39693, 1,
+    0.20104, -0.60161, 1,
+    0.46601, -0.53582, 1,
+    0.67339, -0.53582, 1,
+    -0.13882, 0.54605, 1,
+    -0.29435, 0.77997, 1,
+    -0.26555, 0.96272, 1,
+    -0.16187, 0.8019, 1,
+    -0.17339, 0.64839, 1,
+    -0.28283, 0.47295, 1,
+    -0.36348, 0.31213, 1,
+    -0.30012, 0.027047, 1,
+    -0.23675, -0.21418, 1,
+    -0.06394, -0.18494, 1,
+    0.062788, -0.16301, 1,
+    0.22984, -0.41155, 1,
+    0.2932, -0.2288, 1,
+    0.48329, -0.18494, 1,
+    0.64459, -0.14108, 1,
+    0.46025, 0.012427, 1,
+    0.6273, 0.15863, 1,
+    0.57546, 0.26827, 1,
+    0.72523, 0.44371, 1,
+    0.22408, 0.52412, 1,
+    0.44297, 0.67032, 1,
+    0.322, 0.69225, 1,
+    0.13767, 0.57529, 1,
+    -0.0063364, 0.39985, 1,
+    -0.092742, 0.55336, 1,
+    -0.20795, 0.35599, 1,
+    -0.20795, 0.17325, 1,
+    -0.43836, 0.21711, 1,
+    -0.21947, -0.016813, 1,
+    -0.13882, -0.27266, 1,
+    0.18376, 0.93348, 0,
+    0.22408, 0.77997, 0,
+    0.29896, 0.61915, 0,
+    0.50634, 0.75804, 0,
+    0.61578, 0.7288, 0,
+    0.60426, 0.59722, 0,
+    0.76555, 0.50219, 0,
+    0.92684, 0.3633, 0,
+    0.82316, 0.27558, 0,
+    0.96141, 0.085526, 0,
+    0.93836, 0.012427, 0,
+    0.86348, -0.082602, 0,
+    0.89804, -0.20687, 0,
+    0.85196, -0.36769, 0,
+    0.82892, -0.5212, 0,
+    0.79435, -0.55775, 0,
+    0.59274, -0.7405, 0,
+    0.51786, -0.5943, 0,
+    0.46601, -0.41886, 0,
+    0.35081, -0.57968, 0,
+    0.28744, -0.76974, 0,
+    0.085829, -0.75512, 0,
+    0.14919, -0.57968, 0,
+    -0.13306, -0.4481, 0,
+    -0.40956, -0.41155, 0,
+    -0.39228, -0.25804, 0,
+    -0.74366, -0.25804, 0,
+    -0.69758, 0.041667, 0,
+    -0.75518, 0.2902, 0,
+    -0.69758, 0.68494, 0,
+    -0.4038, 0.70687, 0,
+    -0.38076, 0.91886, 0,
+    -0.50749, 0.90424, 0,
+    -0.54781, 0.70687, 0,
+    0.10311, 0.77997, 0,
+    0.057028, 0.91886, 0,
+    -0.10426, 0.99196, 0,
+    -0.081221, 1.1089, 0,
+    0.28744, 1.087, 0,
+    0.39689, 0.82383, 0,
+    0.63882, 0.88962, 0,
+    0.82316, 0.66301, 0,
+    0.67339, 0.64108, 0,
+    1.0709, 0.10015, 0,
+    -0.046659, -0.57968, 0,
+    -0.23675, -0.63816, 0,
+    -0.15035, -0.36769, 0,
+    -0.49021, -0.3019, 0,
+    -0.46717, -0.13377, 0,
+    -0.28859, -0.060673, 0,
+    -0.61118, -0.067982, 0,
+    -0.66302, -0.21418, 0,
+    -0.59965, -0.41886, 0,
+    -0.72638, -0.082602, 0,
+    -0.83007, 0.31213, 0,
+    -0.72062, 0.53874, 0,
+    -0.59389, 0.49488, 0,
+    -0.48445, 0.99927, 0,
+    -0.0063364, 0.99927, 0,
+    0.63265, -0.030612, 0
+    ]
diff --git a/test/MachineLearning/LeastSquaresModelTest.hs b/test/MachineLearning/LeastSquaresModelTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/LeastSquaresModelTest.hs
@@ -0,0 +1,52 @@
+module MachineLearning.LeastSquaresModelTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.DataSets (dataset1)
+
+import qualified Numeric.LinearAlgebra as LA
+import qualified MachineLearning as ML
+import MachineLearning.Optimization (checkGradient)
+import MachineLearning.Model
+import MachineLearning.LeastSquaresModel
+import MachineLearning.Regularization (Regularization(..))
+
+(x, y) = ML.splitToXY dataset1
+
+x1 = ML.addBiasDimension x
+initialTheta :: LA.Vector LA.R
+initialTheta = LA.konst 1000 (LA.cols x1)
+zeroTheta :: LA.Vector LA.R
+zeroTheta = LA.konst 0 (LA.cols x1)
+
+tests = [ testGroup "model" [
+            testCase "cost, no reg"      $ assertApproxEqual "" 1e-5 1.6190245331702874e12 (cost LeastSquares RegNone x1 y initialTheta)
+            , testCase "cost, lambda = 0"      $ assertApproxEqual "" 1e-5 1.6190245331702874e12 (cost LeastSquares (L2 0) x1 y initialTheta)
+            , testCase "cost, lambda = 1"    $ assertApproxEqual "" 1e-5 1.619024554446883e12 (cost LeastSquares (L2 1) x1 y initialTheta)
+            , testCase "cost, lambda = 1000" $ assertApproxEqual "" 1e-5 1.619045809766032e12 (cost LeastSquares (L2 1000) x1 y initialTheta)
+            , testCase "gradient, no reg" $ assertVector "" 1e-5 gradient_l0 (gradient LeastSquares RegNone x1 y initialTheta)
+            , testCase "gradient, lambda = 0" $ assertVector "" 1e-5 gradient_l0 (gradient LeastSquares (L2 0) x1 y initialTheta)
+            , testCase "gradient, lambda = 1" $ assertVector "" 1e-5 gradient_l1 (gradient LeastSquares (L2 1) x1 y initialTheta)
+            , testCase "gradient, lambda = 1000" $ assertVector "" 1e-5 gradient_l1000 (gradient LeastSquares (L2 1000) x1 y initialTheta)
+            ]
+          , testGroup "gradient checking" [
+              testCase "non-zero theta, no reg" $ assertBool "" $ (checkGradient LeastSquares RegNone x1 y initialTheta 1e-4) < 10
+              , testCase "non-zero theta, non-zero lambda" $ assertBool "" $ (checkGradient LeastSquares (L2 2) x1 y initialTheta 1e-4) < 10
+              , testCase "zero theta, non-zero lambda" $ assertBool "" $ (checkGradient LeastSquares (L2 2) x1 y zeroTheta 1e-4) < 1
+              , testCase "non-zero theta, zero lambda" $ assertBool "" $ (checkGradient LeastSquares (L2 0) x1 y initialTheta 1e-4) < 5
+              , testCase "zero theta, zero lambda" $ assertBool "" $ (checkGradient LeastSquares (L2 0) x1 y zeroTheta 1e-4) < 1
+              ]
+        ]
+
+gradient_l0    = LA.vector [1664438.4042553192,3.865303999468085e9,5567440.808510638]
+gradient_l1    = LA.vector [1664438.4042553192,3.865304020744681e9,5567462.085106383]
+gradient_l1000 = LA.vector [1664438.4042553192,3.86532527606383e9, 5588717.404255319]
diff --git a/test/MachineLearning/LogisticModelTest.hs b/test/MachineLearning/LogisticModelTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/LogisticModelTest.hs
@@ -0,0 +1,78 @@
+module MachineLearning.LogisticModelTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.DataSets (dataset2)
+
+import qualified Numeric.LinearAlgebra as LA
+import qualified MachineLearning as ML
+import MachineLearning.Optimization (checkGradient)
+import MachineLearning.Model
+import MachineLearning.LogisticModel
+import MachineLearning.Regularization (Regularization(..))
+
+(x, y) = ML.splitToXY dataset2
+
+x1 = ML.addBiasDimension $ ML.mapFeatures 6 x
+onesTheta :: LA.Vector LA.R
+onesTheta = LA.konst 1 (LA.cols x1)
+zeroTheta :: LA.Vector LA.R
+zeroTheta = LA.konst 0 (LA.cols x1)
+
+gradientCheckingEps = 1e-3
+
+checkGradientTest lambda theta = do
+  let diffs = take 5 $ map (\e -> checkGradient Logistic lambda x1 y theta e) [1e-3, 1.1e-3 ..]
+      diff = minimum $ filter (not . isNaN) diffs
+  assertApproxEqual (show theta) gradientCheckingEps 0 diff
+
+
+tests = [ testGroup "sigmoid" [
+            testCase "zero" $ assertApproxEqual "" 1e-10 0.5 (sigmoid 0)
+            , testCase "big positive value" $ assertApproxEqual "" 1e-10 1 (sigmoid 10e10)
+            , testCase "big negative value" $ assertApproxEqual "" 1e-10 0 (sigmoid $ -10e10) 
+            , testCase "quite big positive value" $ assertApproxEqual "" 1e-2 1 (sigmoid 100)
+            , testCase "quite big negative value" $ assertApproxEqual "" 1e-2 0 (sigmoid $ -100)
+            ]
+          , testGroup "sigmoidGradient" [
+              testCase "zero" $ assertApproxEqual "" 1e-10 0.25 (sigmoidGradient 0)
+              , testCase "big positive value" $ assertApproxEqual "" 1e-10 0 (sigmoidGradient 10e10)
+              , testCase "big negative value" $ assertApproxEqual "" 1e-10 0 (sigmoidGradient $ -10e10) 
+              , testCase "quite big positive value" $ assertApproxEqual "" 1e-2 0 (sigmoidGradient 100)
+              , testCase "quite big negative value" $ assertApproxEqual "" 1e-2 0 (sigmoidGradient $ -100)
+              , testCase "small positive value" $ assertApproxEqual "" 1e-2 0.2 (sigmoidGradient 1)
+              , testCase "small negative value" $ assertApproxEqual "" 1e-2 0.2 (sigmoidGradient $ -1) 
+              ]
+          , testGroup "model" [
+              testCase "cost, lambda = 0" $ assertApproxEqual "" 1e-3 2.020 (cost Logistic (L2 0) x1 y onesTheta)
+              , testCase "cost, lambda = 1" $ assertApproxEqual "" 1e-3 2.135 (cost Logistic (L2 1) x1 y onesTheta)
+              , testCase "cost, lambda = 1000" $ assertApproxEqual "" 1e-3 116.427 (cost Logistic (L2 1000) x1 y onesTheta)
+              , testCase "gradient, lambda = 0" $ assertVector "" 1e-5 gradient_l0 (gradient Logistic (L2 0) x1 y onesTheta)
+              , testCase "gradient, lambda = 1" $ assertVector "" 1e-5 gradient_l1 (gradient Logistic (L2 1) x1 y onesTheta)
+              , testCase "gradient, lambda = 1000" $ assertVector "" 1e-5 gradient_l1000 (gradient Logistic (L2 1000) x1 y onesTheta)
+              ]
+          , testGroup "gradient checking" [
+              testCase "non-zero theta, non-zero lambda" $ checkGradientTest (L2 2) onesTheta
+              , testCase "zero theta, non-zero lambda" $ checkGradientTest (L2 2) zeroTheta
+              , testCase "non-zero theta, zero lambda" $ checkGradientTest (L2 0) onesTheta
+              , testCase "zero theta, zero lambda" $ checkGradientTest (L2 0) zeroTheta
+              , testCase "non-zero theta, no reg" $ checkGradientTest RegNone onesTheta
+              , testCase "zero theta, no reg" $ checkGradientTest RegNone zeroTheta
+
+              ]
+        ]
+
+gradient_l0 = LA.vector [0.34604507367924525,7.660615656904722e-2,0.11004999290013262,0.14211701951318526,7.4399123914273965e-3,0.15963981400206023,5.864636026438637e-2,2.369595228045015e-2,1.7568631688383615e-2,9.87226947583087e-2,8.878427090266403e-2,2.509755728155993e-3,3.348199373717313e-2,1.0975419586624715e-3,0.11520318246520422,5.0480764463147594e-2,1.0229511332420786e-2,8.818652179760128e-3,1.5052078581608905e-2,6.655810974747264e-3,9.010665405440292e-2,6.480865498500844e-2,2.039892642614106e-3,1.4231095091583643e-2,5.737477462013599e-4,1.71609000731311e-2,-2.437841009425525e-4,9.753746346629336e-2]
+
+gradient_l1 = LA.vector [0.34604507367924525,8.508073284023365e-2,0.11852456917131905,0.15059159578437167,1.5914488662613836e-2,0.16811439027324665,6.71209365355728e-2,3.217052855163659e-2,2.6043207959570054e-2,0.10719727102949513,9.725884717385046e-2,1.0984331999342433e-2,4.1956570008359576e-2,9.572118229848912e-3,0.12367775873639067,5.895534073433403e-2,1.8704087603607224e-2,1.729322845094657e-2,2.3526654852795346e-2,1.5130387245933704e-2,9.858123032558937e-2,7.328323125619488e-2,1.0514468913800546e-2,2.270567136277008e-2,9.048324017387801e-3,2.563547634431754e-2,8.230792170243889e-3,0.10601203973747979]
+
+gradient_l1000 = LA.vector [0.34604507367924525,8.551182427755489,8.584626264086573,8.616693290699626,8.48201618357787,8.6342160851885,8.533222631450828,8.498272223466891,8.492144902874823,8.57329896594475,8.563360542089104,8.477086026914597,8.508058264923614,8.475673813145104,8.589779453651644,8.525057035649588,8.484805782518862,8.483394923366202,8.489628349768049,8.481232082161188,8.564682925240843,8.539384926171449,8.476616163829055,8.488807366278024,8.475150018932643,8.491737171259572,8.474332487085498,8.572113734652733]
diff --git a/test/MachineLearning/MultiSvmClassifierTest.hs b/test/MachineLearning/MultiSvmClassifierTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/MultiSvmClassifierTest.hs
@@ -0,0 +1,81 @@
+module MachineLearning.MultiSvmClassifierTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.DataSets (dataset2)
+
+import qualified Numeric.LinearAlgebra as LA
+import qualified MachineLearning as ML
+import MachineLearning.Optimization
+import MachineLearning.Model
+import MachineLearning.MultiSvmClassifier
+
+(x, y) = ML.splitToXY dataset2
+
+model = MultiClass (MultiSvm 1 2)
+
+x1 = ML.addBiasDimension x
+onesTheta :: LA.Vector LA.R
+onesTheta = LA.konst 1 (2 * LA.cols x1)
+zeroTheta :: LA.Vector LA.R
+zeroTheta = LA.konst 0 (2 * LA.cols x1)
+
+processX muSigma x = ML.addBiasDimension $ ML.featureNormalization muSigma $ ML.mapFeatures 6 x
+
+muSigma = ML.meanStddev (ML.mapFeatures 6 x)
+x2 = processX muSigma x
+
+
+xPredict = LA.matrix 2 [ -0.5, 0.5
+                       , 0.2, -0.2
+                       , 1, 1
+                       , 1, 0
+                       , 0, 0
+                       , 0, 1]
+xPredict2 = processX muSigma xPredict
+yExpected = LA.vector [1, 1, 0, 0, 1, 0]
+
+
+gradientCheckingEps :: Double
+gradientCheckingEps = 1e-3
+
+eps = 0.0001
+
+zeroTheta2 = LA.konst 0 (2 * LA.cols x2)
+(thetaGD, _) = minimize (GradientDescent 0.001) model eps 150 (L2 1) x2 y zeroTheta2
+(thetaCGFR, _) = minimize (ConjugateGradientFR 0.1 0.1) model eps 30 (L2 0.5) x2 y zeroTheta2
+(thetaCGPR, _) = minimize (ConjugateGradientPR 0.1 0.1) model eps 30 (L2 0.5) x2 y zeroTheta2
+(thetaBFGS, _) = minimize (BFGS2 0.1 0.1) model eps 30 (L2 0.5) x2 y zeroTheta2
+
+
+checkGradientTest lambda theta eps = do
+  let diffs = take 5 $ map (\e -> checkGradient model lambda x1 y theta e) [1e-3, 1.1e-3 ..]
+      diff = minimum $ filter (not . isNaN) diffs
+  assertApproxEqual "" eps 0 diff
+
+
+tests = [  testGroup "gradient checking" [
+            testCase "non-zero theta, non-zero lambda" $ checkGradientTest (L2 2) onesTheta 3e-2
+              , testCase "zero theta, non-zero lambda" $ checkGradientTest (L2 2) zeroTheta gradientCheckingEps
+              , testCase "non-zero theta, zero lambda" $ checkGradientTest (L2 0) onesTheta gradientCheckingEps
+              , testCase "zero theta, zero lambda" $ checkGradientTest (L2 0) zeroTheta gradientCheckingEps
+              , testCase "non-zero theta, no reg" $ checkGradientTest RegNone onesTheta gradientCheckingEps
+              , testCase "zero theta, no reg" $ checkGradientTest RegNone zeroTheta gradientCheckingEps
+              ]
+        , testGroup "learn" [
+            testCase "Gradient Descent" $ assertVector "" 0.01 yExpected (hypothesis model xPredict2 thetaGD)
+            , testCase "Conjugate Gradient FR" $ assertVector "" 0.01 yExpected (hypothesis model xPredict2 thetaCGFR)
+            , testCase "Conjugate Gradient PR" $ assertVector "" 0.01 yExpected (hypothesis model xPredict2 thetaCGPR)
+            , testCase "BFGS" $ assertVector "" 0.01 yExpected (hypothesis model xPredict2 thetaBFGS)
+            ]
+        ]
+
diff --git a/test/MachineLearning/NeuralNetwork/TopologyTest.hs b/test/MachineLearning/NeuralNetwork/TopologyTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/NeuralNetwork/TopologyTest.hs
@@ -0,0 +1,28 @@
+module MachineLearning.NeuralNetwork.TopologyTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+import qualified Numeric.LinearAlgebra as LA
+import MachineLearning.NeuralNetwork.Topology
+import qualified MachineLearning.NeuralNetwork.TopologyMaker as TM
+
+nnt = TM.makeTopology TM.ASigmoid TM.LLogistic 15 2 [10]
+
+flattenTest = do
+  theta <- initializeThetaIO nnt
+  let theta' = flatten $ unflatten nnt theta
+      norm = LA.norm_2 (theta - theta')
+  assertApproxEqual "flatten" 1e-10 0 norm
+
+tests = [ testGroup "flatten" [
+            testCase "flatten" flattenTest
+            ]
+        ]
diff --git a/test/MachineLearning/NeuralNetwork/WeightInitializationTest.hs b/test/MachineLearning/NeuralNetwork/WeightInitializationTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/NeuralNetwork/WeightInitializationTest.hs
@@ -0,0 +1,50 @@
+module MachineLearning.NeuralNetwork.WeightInitializationTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+import Control.Monad (replicateM)
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+import qualified Numeric.GSL.Statistics as Stat
+import qualified Control.Monad.Random as RndM
+import MachineLearning.NeuralNetwork.WeightInitialization
+
+assertDistribution eps x mu e = do
+  let mean = Stat.mean x
+  assertApproxEqual "mean" eps mu mean
+  assertBool "maximum" $ (V.maximum x) <= e
+  assertBool "minimum" $ (V.minimum x) >= (-e)
+
+
+generateData algo n sz = do
+   rndList <- replicateM n $ RndM.evalRandIO (algo sz)
+   let (bs, ws) = unzip rndList
+       b = LA.flatten $ LA.fromBlocks [bs]
+       w = LA.flatten $ LA.fromBlocks [ws]
+   return (b, w)
+
+
+testWeightInitAlgo eps algo n sz mu e = do
+  (b, w) <- generateData algo n sz
+  assertDistribution eps b 0 0
+  assertDistribution eps w mu e
+
+
+heEps (r, c) = sqrt (2/(fromIntegral $ r + c))
+nguyenEps (r, c) = (sqrt 6) / (sqrt . fromIntegral $ r + c)
+
+sz = (7, 5)
+
+tests = [ testGroup "flatten" [
+            testCase "hu" $ testWeightInitAlgo 1e-1 he 100 sz 0 (heEps sz)
+            , testCase "nguyen" $ testWeightInitAlgo 1e-1 nguyen 100 sz 0 (nguyenEps sz)
+            ]
+        ]
diff --git a/test/MachineLearning/NeuralNetworkTest.hs b/test/MachineLearning/NeuralNetworkTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/NeuralNetworkTest.hs
@@ -0,0 +1,72 @@
+module MachineLearning.NeuralNetworkTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.DataSets (dataset2)
+
+import qualified Control.Monad.Random as RndM
+import qualified Numeric.LinearAlgebra as LA
+import qualified MachineLearning as ML
+import qualified MachineLearning.Optimization as Opt
+import MachineLearning.Model
+import MachineLearning.NeuralNetwork
+import qualified MachineLearning.NeuralNetwork.TopologyMaker as TM
+
+(x, y) = ML.splitToXY dataset2
+
+gradientCheckingEps = 0.1
+
+checkGradientTest eps activation loss lambda = do
+  let nnt = TM.makeTopology activation loss (LA.cols x) 2 [10]
+      model = NeuralNetwork nnt
+      thetas = initializeTheta 1511197 nnt
+      diffs = take 5 $ map (\e -> Opt.checkGradient model lambda x y thetas e) [0.005, 0.0051 ..]
+      diff = minimum $ filter (not . isNaN) diffs
+  assertApproxEqual (show thetas) eps 0 diff
+
+
+xPredict = LA.matrix 2 [ -0.5, 0.5
+                       , 0.2, -0.2
+                       , 1, 1
+                       , 1, 0
+                       , 0, 0]
+yExpected = LA.vector [1, 1, 0, 0, 1]
+
+learnTest activation loss minMethod nIters =
+  let lambda = L2 $ 0.5 / (fromIntegral $ LA.rows x)
+      x1 = ML.mapFeatures 2 x
+      nnt = TM.makeTopology activation loss (LA.cols x1) 2 [10]
+      model = NeuralNetwork nnt
+      xPredict1 = ML.mapFeatures 2 xPredict
+      initTheta = initializeTheta 5191711 nnt
+      (theta, optPath) = Opt.minimize minMethod model 1e-7 nIters lambda x1 y initTheta
+      yPredicted = hypothesis model xPredict1 theta
+      js = (LA.toColumns optPath) !! 1
+  in do
+    assertVector (show js) 0.01 yExpected yPredicted
+
+
+tests = [ testGroup "gradient checking" [
+            testCase "Sigmoid - Logistic: non-zero lambda" $ checkGradientTest 0.1 TM.ASigmoid TM.LLogistic (L2 0.01)
+            , testCase "Sigmoid - Logistic: zero lambda" $ checkGradientTest 0.1 TM.ASigmoid TM.LLogistic (L2 0)
+            , testCase "ReLU - Softmax: non-zero lambda" $ checkGradientTest 0.1 TM.ARelu TM.LSoftmax (L2 0.01)
+            , testCase "ReLU - Softmax: zero lambda" $ checkGradientTest 0.1 TM.ARelu TM.LSoftmax (L2 0)
+            , testCase "Tanh - MultiSvm: non-zero lambda" $ checkGradientTest 0.1 TM.ATanh TM.LMultiSvm (L2 0.01)
+            , testCase "Tanh - MultiSvm: zero lambda" $ checkGradientTest 0.1 TM.ATanh TM.LMultiSvm (L2 0)
+            , testCase "Tanh - MultiSvm: no reg" $ checkGradientTest 0.1 TM.ATanh TM.LMultiSvm RegNone
+            ]
+        , testGroup "learn" [
+            testCase "Sigmoid - Logistic: BFGS" $ learnTest TM.ASigmoid TM.LLogistic (Opt.BFGS2 0.01 0.7) 50
+            , testCase "ReLU - Softmax: BFGS" $ learnTest TM.ARelu TM.LSoftmax (Opt.BFGS2 0.1 0.1) 50
+            , testCase "Tanh - MultiSvm: BFGS" $ learnTest TM.ATanh TM.LMultiSvm (Opt.BFGS2 0.1 0.1) 50
+            ]
+        ]
diff --git a/test/MachineLearning/Optimization/GradientDescentTest.hs b/test/MachineLearning/Optimization/GradientDescentTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/Optimization/GradientDescentTest.hs
@@ -0,0 +1,46 @@
+module MachineLearning.Optimization.GradientDescentTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+import qualified MachineLearning as ML
+import MachineLearning.Regularization (Regularization(..))
+import MachineLearning.LeastSquaresModel (LeastSquaresModel(..))
+import MachineLearning.Optimization.GradientDescent
+
+import MachineLearning.DataSets (dataset1)
+
+(x, y) = ML.splitToXY dataset1
+
+muSigma = ML.meanStddev x
+xNorm = ML.featureNormalization muSigma x
+x1 = ML.addBiasDimension xNorm
+initialTheta = LA.konst 0 (LA.cols x1)
+lsExpectedTheta = LA.vector [340412.660, 110630.879, -8737.743]
+eps = 1e-3
+
+
+isInDescendingOrder :: [Double] -> Bool
+isInDescendingOrder lst = and . snd . unzip $ scanl (\(prev, _) current -> (current, prev >= current)) (1/0, True) lst
+
+testGradientDescent model expectedTheta = do
+  let (theta, optPath) = gradientDescent 0.01 model eps 5000 RegNone x1 y initialTheta
+      js = V.toList $ (LA.toColumns optPath) !! 1
+  assertVector "theta" 0.01 expectedTheta theta
+  assertBool "non-increasing errors" $ isInDescendingOrder js
+
+tests = [testGroup "gradientDescent" [
+            testCase "leastSquares" $ testGradientDescent LeastSquares lsExpectedTheta
+            ]
+        ]
diff --git a/test/MachineLearning/Optimization/MinibatchGradientDescentTest.hs b/test/MachineLearning/Optimization/MinibatchGradientDescentTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/Optimization/MinibatchGradientDescentTest.hs
@@ -0,0 +1,48 @@
+module MachineLearning.Optimization.MinibatchGradientDescentTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import qualified Data.Vector.Storable as V
+import qualified Numeric.LinearAlgebra as LA
+
+import MachineLearning.Types (Vector)
+import qualified MachineLearning as ML
+import MachineLearning.Regularization (Regularization(..))
+import MachineLearning.LeastSquaresModel (LeastSquaresModel(..))
+import MachineLearning.Optimization.MinibatchGradientDescent
+
+import MachineLearning.DataSets (dataset1)
+
+(x, y) = ML.splitToXY dataset1
+
+muSigma = ML.meanStddev x
+xNorm = ML.featureNormalization muSigma x
+x1 = ML.addBiasDimension xNorm
+initialTheta :: Vector
+initialTheta = LA.konst 0 (LA.cols x1)
+lsExpectedTheta = LA.vector [325009.354,113890.981,6876.935]
+eps = 1e-3
+
+
+isInDescendingOrder :: [Double] -> Bool
+isInDescendingOrder lst = and . snd . unzip $ scanl (\(prev, _) current -> (current, prev >= current)) (1/0, True) lst
+
+testMinibatchGradientDescent model expectedTheta = do
+  let (theta, optPath) = minibatchGradientDescent 0 16 0.01 model eps 5000 RegNone x1 y initialTheta
+      js = V.toList $ (LA.toColumns optPath) !! 1
+  assertVector "theta" 0.01 expectedTheta theta
+  assertBool "non-increasing errors" $ isInDescendingOrder js
+
+tests = [testGroup "minibatchGradientDescent" [
+            testCase "leastSquares" $ testMinibatchGradientDescent LeastSquares lsExpectedTheta
+            ]
+        ]
diff --git a/test/MachineLearning/PCATest.hs b/test/MachineLearning/PCATest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/PCATest.hs
@@ -0,0 +1,48 @@
+module MachineLearning.PCATest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+
+import qualified Numeric.LinearAlgebra as LA
+
+import MachineLearning.PCA
+
+getDimReducerSmokeTest =
+  let nFeatures = 4
+      nExamples = 7
+      m = LA.matrix nFeatures [1 .. fromIntegral $ nFeatures*nExamples]
+      m10 = m * 10
+      (reduceDims, retainedVariance, mReduced) = getDimReducer m 2
+      m10Reduced = reduceDims m10
+  in do
+    assertEqual "dimension equality (getDimReducer)" (LA.cols mReduced) 2
+    assertEqual "dimension equality (reduceDims)" (LA.cols m10Reduced) 2
+    assertApproxEqual "retained variance" 1e-10 1 retainedVariance
+
+getDimReducer_rvSmokeTest rv=
+  let nFeatures = 4
+      nExamples = 7
+      m = LA.matrix nFeatures [1 .. fromIntegral $ nFeatures*nExamples]
+      m10 = m * 10
+      (reduceDims, k, mReduced) = getDimReducer_rv m rv
+      m10Reduced = reduceDims m10
+  in do
+    assertEqual "dimension equality (getDimReducer_rv)" (LA.cols mReduced) k
+    assertEqual "dimension equality (reduceDims_rv)" (LA.cols m10Reduced) k
+    assertBool "reduced number of dimensions" $ k <= nFeatures
+
+tests = [ testGroup "smoke test" [
+            testCase "getDimReducer" getDimReducerSmokeTest
+            , testCase "getDimReducer_rv, rv = 0" $ getDimReducer_rvSmokeTest 0
+            , testCase "getDimReducer_rv, rv = 0.5" $ getDimReducer_rvSmokeTest 0.5
+            , testCase "getDimReducer_rv, rv = 1" $ getDimReducer_rvSmokeTest 1
+            , testCase "getDimReducer_rv, rv = 2" $ getDimReducer_rvSmokeTest 2
+            ]
+        ]
diff --git a/test/MachineLearning/RandomTest.hs b/test/MachineLearning/RandomTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/RandomTest.hs
@@ -0,0 +1,73 @@
+module MachineLearning.RandomTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.Random
+
+import Data.List (nub)
+import Control.Monad (foldM_)
+import qualified Data.Vector as V
+import qualified Data.Vector.Storable as SV
+import qualified Numeric.LinearAlgebra as LA
+import qualified System.Random as Rnd
+import qualified Control.Monad.Random as RndM
+
+sampleTest = do
+  gen <- Rnd.newStdGen
+  foldM_ sampleTestIter gen [1..25]
+
+sampleTestIter gen i =
+  let xs = V.fromList [1..100]
+      n = 10 + i
+      (ys, gen') = sample gen n xs
+  in do
+    assertEqual "uniqness" (V.length xs) (length . nub $ V.toList xs) 
+    assertEqual "length" n (V.length ys)
+    assertBool "maximum" $ (V.maximum ys) <= (V.maximum xs)
+    assertBool "minimum" $ (V.minimum ys) >= (V.minimum xs)
+    assertEqual "uniqness of elements" n (length . nub $ V.toList ys) 
+    return gen'
+
+
+randomRListTest = mapM_ (randomRListTestIter ((-1000, 1000)::(Int, Int))) [10..30]
+
+randomRListTestIter range@(lo, hi) len = do
+  rndList <- RndM.evalRandIO (getRandomRListM len range)
+  assertEqual "length" len (length rndList)
+  assertBool "minimum" $ lo <= (minimum rndList)
+  assertBool "maximum" $ hi >= (maximum rndList)
+
+
+randomRVectorTest = mapM_ (randomRVectorTestIter ((-2, 2))) [10..30]
+
+randomRVectorTestIter range@(lo, hi) len = do
+  rndVector <- RndM.evalRandIO (getRandomRVectorM len range)
+  assertEqual "length" len (SV.length rndVector)
+  assertBool "minimum" $ lo <= (SV.minimum rndVector)
+  assertBool "maximum" $ hi >= (SV.maximum rndVector)
+
+
+randomRMatrixTest = mapM_ (randomRMatrixTestIter ((-2, 2))) $ zip [10..15] [12..17]
+
+randomRMatrixTestIter range@(lo, hi) (rows, cols) = do
+  rndMatrix <- RndM.evalRandIO (getRandomRMatrixM rows cols range)
+  assertEqual "rows" rows (LA.rows rndMatrix)
+  assertEqual "columns" cols (LA.cols rndMatrix)
+  assertBool "minimum" $ lo <= (LA.minElement rndMatrix)
+  assertBool "maximum" $ hi >= (LA.maxElement rndMatrix)
+
+
+tests = [ testCase "sample" sampleTest
+        , testCase "getRandomRList" randomRListTest
+        , testCase "getRandomRVector" randomRVectorTest
+        , testCase "getRandomRMatrix" randomRMatrixTest
+        ]
diff --git a/test/MachineLearning/RegressionTest.hs b/test/MachineLearning/RegressionTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/RegressionTest.hs
@@ -0,0 +1,53 @@
+module MachineLearning.RegressionTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import qualified Numeric.LinearAlgebra as LA
+
+import MachineLearning.DataSets (dataset1)
+
+import qualified MachineLearning as ML
+import MachineLearning.Regression
+
+(x, y) = ML.splitToXY dataset1
+
+muSigma = ML.meanStddev x
+xNorm = ML.featureNormalization muSigma x
+x1 = ML.addBiasDimension xNorm
+zeroTheta = LA.konst 0 (LA.cols x1)
+
+xPredict = LA.matrix 2 [1650, 3]
+xPredict1 = ML.addBiasDimension $ ML.featureNormalization muSigma xPredict
+
+theta = normalEquation (ML.addBiasDimension x) y
+yExpected = hypothesis LeastSquares (ML.addBiasDimension xPredict) theta
+
+eps = 0.0001
+thetaNE = normalEquation x1 y
+thetaNE_p = normalEquation_p x1 y
+(thetaGD, _) = minimize (GradientDescent 0.01) LeastSquares eps 5000 RegNone x1 y zeroTheta
+(thetaMBGD, _) = minimize (MinibatchGradientDescent 11711 64 0.05) LeastSquares eps 5000 RegNone x1 y zeroTheta
+(thetaCGFR, _) = minimize (ConjugateGradientFR 0.1 0.1) LeastSquares eps 1500 RegNone x1 y zeroTheta
+(thetaCGPR, _) = minimize (ConjugateGradientPR 0.1 0.1) LeastSquares eps 1500 RegNone x1 y zeroTheta
+(thetaBFGS, _) = minimize (BFGS2 0.1 0.1) LeastSquares eps 1500 RegNone x1 y zeroTheta
+
+
+tests = [ testGroup "minimize" [
+            testCase "Normal Equation" $ assertVector "" 0.01 yExpected (hypothesis LeastSquares xPredict1 thetaNE)
+            , testCase "Normal Equation using pseudo inverse" $ assertVector "" 0.01 yExpected (hypothesis LeastSquares xPredict1 thetaNE_p)
+            , testCase "Gradient Descent" $ assertVector "" 0.01 yExpected (hypothesis LeastSquares xPredict1 thetaGD)
+            , testCase "Minibatch Gradient Descent" $ assertVector "" 1100 yExpected (hypothesis LeastSquares xPredict1 thetaMBGD)
+            , testCase "BFGS" $ assertVector "" 0.01 yExpected (hypothesis LeastSquares xPredict1 thetaBFGS)
+            , testCase "Conjugate Gradient FR" $ assertVector "" 0.01 yExpected (hypothesis LeastSquares xPredict1 thetaCGFR)
+            , testCase "Conjugate Gradient PR" $ assertVector "" 0.01 yExpected (hypothesis LeastSquares xPredict1 thetaCGPR)
+            ]
+        ]
diff --git a/test/MachineLearning/SoftmaxClassifierTest.hs b/test/MachineLearning/SoftmaxClassifierTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/SoftmaxClassifierTest.hs
@@ -0,0 +1,79 @@
+module MachineLearning.SoftmaxClassifierTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.DataSets (dataset2)
+
+import qualified Numeric.LinearAlgebra as LA
+import qualified MachineLearning as ML
+import MachineLearning.Optimization
+import MachineLearning.SoftmaxClassifier
+
+(x, y) = ML.splitToXY dataset2
+
+model = MultiClass (Softmax 2)
+
+x1 = ML.addBiasDimension x
+onesTheta :: LA.Vector LA.R
+onesTheta = LA.konst 1 (2 * LA.cols x1)
+zeroTheta :: LA.Vector LA.R
+zeroTheta = LA.konst 0 (2 * LA.cols x1)
+
+processX muSigma x = ML.addBiasDimension $ ML.featureNormalization muSigma $ ML.mapFeatures 6 x
+
+muSigma = ML.meanStddev (ML.mapFeatures 6 x)
+x2 = processX muSigma x
+
+
+xPredict = LA.matrix 2 [ -0.5, 0.5
+                       , 0.2, -0.2
+                       , 1, 1
+                       , 1, 0
+                       , 0, 0
+                       , 0, 1]
+xPredict2 = processX muSigma xPredict
+yExpected = LA.vector [1, 1, 0, 0, 1, 0]
+
+
+checkSoftmaxGradient theta eps lambda = minimum . take 5 . filter (not . isNaN) $ map check [eps, eps+0.001 ..]
+  where check e = checkGradient model lambda x1 y theta e
+  
+
+gradientCheckingEps :: Double
+gradientCheckingEps = 3e-2
+
+eps = 0.000001
+
+initialTheta = LA.konst 0.001 (2 * LA.cols x2)
+(thetaGD, optPathGD) = minimize (GradientDescent 0.0005) model eps 150 (L2 1) x2 y initialTheta
+(thetaCGFR, optPathCGFR) = minimize (ConjugateGradientFR 0.05 0.2) model eps 30 (L2 1) x2 y initialTheta
+(thetaCGPR, optPathCGPR) = minimize (ConjugateGradientPR 0.05 0.3) model eps 30 (L2 1) x2 y initialTheta
+
+showOptPath optPath = show $  (LA.toColumns optPath) !! 1
+
+
+tests = [  testGroup "gradient checking" [
+              testCase "non-zero theta, non-zero lambda" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient onesTheta 1e-3 (L2 2))
+              , testCase "zero theta, non-zero lambda" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient zeroTheta 1e-3 (L2 2))
+              , testCase "non-zero theta, zero lambda" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient onesTheta 1e-3 (L2 0))
+              , testCase "zero theta, zero lambda" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient zeroTheta 1e-3 (L2 0))
+              , testCase "non-zero theta, no reg" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient onesTheta 1e-3 RegNone)
+              , testCase "zero theta, no reg" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient zeroTheta 1e-3 RegNone)
+              ]
+           
+        , testGroup "learn" [
+            testCase "Gradient Descent" $ assertVector (showOptPath optPathGD) 0.01 yExpected (hypothesis model xPredict2 thetaGD)
+            , testCase "Conjugate Gradient FR" $ assertVector (showOptPath optPathCGFR) 0.01 yExpected (hypothesis model xPredict2 thetaCGFR)
+            , testCase "Conjugate Gradient PR" $ assertVector (showOptPath optPathCGPR) 0.01 yExpected (hypothesis model xPredict2 thetaCGPR)
+            ]
+        ]
+
diff --git a/test/MachineLearning/UtilsTest.hs b/test/MachineLearning/UtilsTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearning/UtilsTest.hs
@@ -0,0 +1,75 @@
+module MachineLearning.UtilsTest
+(
+  tests
+)
+
+where
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Approx
+import Test.HUnit.Plus
+
+import MachineLearning.Types (Vector, Matrix)
+import qualified Numeric.LinearAlgebra as LA
+import qualified Data.Vector.Storable as V
+import Numeric.LinearAlgebra ((><))
+
+import MachineLearning.Utils
+
+a :: Matrix
+a = (4><5) [ 1, 7, 9, 11, 2
+           , 77, 4, 6, 9, 0
+           , -11, -55, 3, 11, 55
+           , 7, 9, 11, 13, 15
+           ]
+
+sumRows :: Matrix
+sumRows = (4><1) [ 30
+                 , 96
+                 , 3
+                 , 55
+                 ]
+
+
+sumCols :: Matrix
+sumCols = (1><5) [74, -35, 29, 44, 72]
+
+
+maxRows :: Matrix
+maxRows = (4><1) [ 11
+                 , 77
+                 , 55
+                 , 15
+                 ]
+
+
+minCols :: Matrix
+minCols = (1><5) [-11, -55, 3, 9, 0]
+
+
+sumRowsV = LA.flatten sumRows
+sumColsV = LA.flatten sumCols
+maxRowsV = LA.flatten maxRows
+minColsV = LA.flatten minCols
+
+
+tests = [ testGroup "reduceV" [
+            testCase "reduceByRowsV: sum" $ assertVector "" 1e-10 sumRowsV (reduceByRowsV V.sum a)
+            , testCase "reduceByColumnsV: sum" $ assertVector "" 1e-10 sumColsV (reduceByColumnsV V.sum a)
+            , testCase "reduceByRowsV: max" $ assertVector "" 1e-10 maxRowsV (reduceByRowsV V.maximum a)
+            , testCase "reduceByColumnsV: min" $ assertVector "" 1e-10 minColsV (reduceByColumnsV V.minimum a)
+            ]
+          , testGroup "reduce" [
+              testCase "reduceByRows: sum" $ assertMatrix "" 1e-10 sumRows (reduceByRows V.sum a)
+              , testCase "reduceByColumns: sum" $ assertMatrix "" 1e-10 sumCols (reduceByColumns V.sum a)
+              , testCase "reduceByRows: max" $ assertMatrix "" 1e-10 maxRows (reduceByRows V.maximum a)
+              , testCase "reduceByColumns: min" $ assertMatrix "" 1e-10 minCols (reduceByColumns V.minimum a)
+            ]
+          , testGroup "sum" [
+              testCase "sumByRows" $ assertMatrix "" 1e-10 sumRows (sumByRows a)
+              , testCase "sumColumns" $ assertMatrix "" 1e-10 sumCols (sumByColumns a)
+              ]
+          , testCase "listOfTuplesToList" $ assertEqual "" [11, 9, 27, 3, 43, 11] (listOfTuplesToList [(11, 9), (27, 3), (43, 11)])
+        ]
diff --git a/test/MachineLearningTest.hs b/test/MachineLearningTest.hs
new file mode 100644
--- /dev/null
+++ b/test/MachineLearningTest.hs
@@ -0,0 +1,28 @@
+module MachineLearningTest
+(
+  tests
+)
+
+where
+
+
+import Test.Framework (testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit
+import Test.HUnit.Plus
+
+import MachineLearning
+
+import qualified Numeric.LinearAlgebra as LA
+
+x = LA.matrix 3 [17, 19, 29]
+x2 = LA.matrix 9 [17, 19, 29, 289, 323, 493, 361, 551, 841]
+x3 = LA.matrix 19 [17, 19, 29, 289, 323, 493, 361, 551, 841, 4913, 5491, 8381, 6137, 9367, 14297, 6859, 10469, 15979, 24389]
+eps = 1e-5
+
+tests = [testGroup "mapFeatures" [
+            testCase "mapfeatures 1" $ assertMatrix "" eps x (mapFeatures 1 x)
+            , testCase "mapfeatures 2" $ assertMatrix "" eps x2 (mapFeatures 2 x)
+            , testCase "mapfeatures 3" $ assertMatrix "" eps x3 (mapFeatures 3 x)
+            ]
+        ]
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,42 @@
+import Test.Framework (defaultMain, testGroup)
+
+
+import qualified MachineLearningTest as MachineLearning
+import qualified MachineLearning.RegressionTest as Regression
+import qualified MachineLearning.Classification.BinaryTest as Classification.Binary
+import qualified MachineLearning.Classification.OneVsAllTest as Classification.OneVsAll
+import qualified MachineLearning.LeastSquaresModelTest as LeastSquaresModel
+import qualified MachineLearning.LogisticModelTest as LogisticModel
+import qualified MachineLearning.MultiSvmClassifierTest as MultiSvmClassifier
+import qualified MachineLearning.SoftmaxClassifierTest as SoftmaxClassifier
+import qualified MachineLearning.Optimization.GradientDescentTest as GradientDescent
+import qualified MachineLearning.Optimization.MinibatchGradientDescentTest as MinibatchGradientDescent
+import qualified MachineLearning.NeuralNetworkTest as NeuralNetwork
+import qualified MachineLearning.NeuralNetwork.TopologyTest as NeuralNetwork.Topology
+import qualified MachineLearning.NeuralNetwork.WeightInitializationTest as NeuralNetwork.WeightInitialization
+import qualified MachineLearning.PCATest as PCA
+import qualified MachineLearning.ClusteringTest as Clustering
+import qualified MachineLearning.RandomTest as Random
+import qualified MachineLearning.UtilsTest as Utils
+
+main = defaultMain tests
+
+tests = [
+  testGroup "MachineLearning" MachineLearning.tests
+  , testGroup "MachineLearning.Regression" Regression.tests
+  , testGroup "MachineLearning.Classification.Binary" Classification.Binary.tests
+  , testGroup "MachineLearning.Classification.OneVsAll" Classification.OneVsAll.tests
+  , testGroup "MachineLearning.LeastSquaresModel" LeastSquaresModel.tests
+  , testGroup "MachineLearning.LogisticModel" LogisticModel.tests
+  , testGroup "MachineLearning.MultiSvmClassifier" MultiSvmClassifier.tests
+  , testGroup "MachineLearning.SoftmaxClassifier" SoftmaxClassifier.tests
+  , testGroup "MachineLearning.Optimization.GradientDescent" GradientDescent.tests
+  , testGroup "MachineLearning.Optimization.MinibatchGradientDescent" MinibatchGradientDescent.tests
+  , testGroup "MachineLearning.NeuralNetwork" NeuralNetwork.tests
+  , testGroup "MachineLearning.NeuralNetwork.Topology" NeuralNetwork.Topology.tests
+  , testGroup "MachineLearning.NeuralNetwork.WeightInitialization" NeuralNetwork.WeightInitialization.tests
+  , testGroup "MachineLearning.PCA" PCA.tests
+  , testGroup "MachineLearning.Clustering" Clustering.tests
+  , testGroup "MachineLearning.Random" Random.tests
+  , testGroup "MachineLearning.Utils" Utils.tests
+  ]
diff --git a/test/Test/HUnit/Approx.hs b/test/Test/HUnit/Approx.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/HUnit/Approx.hs
@@ -0,0 +1,85 @@
+{-# LANGUAGE ImplicitParams, CPP #-}
+#if __GLASGOW_HASKELL__ >= 707
+{-# LANGUAGE Safe #-}       -- Test.HUnit is not Safe in 7.6 and below
+#endif
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Test.HUnit.Approx
+-- Copyright   :  (C) 2014 Richard Eisenberg
+-- License     :  BSD-style (see LICENSE)
+-- Maintainer  :  Richard Eisenberg (eir@cis.upenn.edu)
+-- Stability   :  intended to be stable
+-- Portability :  not portable (uses implicit parameters)
+--
+-- This module exports combinators to allow approximate equality of
+-- floating-point values in HUnit tests.
+-----------------------------------------------------------------------------
+
+module Test.HUnit.Approx (
+  -- * Assertions
+  assertApproxEqual, (@~?), (@?~),
+
+  -- * Tests
+  (~~?), (~?~)
+  ) where
+
+import Test.HUnit
+import Control.Monad  ( unless )
+
+-- | Asserts that the specified actual value is approximately equal to the
+-- expected value. The output message will contain the prefix, the expected
+-- value, the actual value, and the maximum margin of error.
+--  
+-- If the prefix is the empty string (i.e., @\"\"@), then the prefix is omitted
+-- and only the expected and actual values are output.
+assertApproxEqual :: (Ord a, Num a, Show a)
+                  => String -- ^ The message prefix
+                  -> a      -- ^ Maximum allowable margin of error
+                  -> a      -- ^ The expected value 
+                  -> a      -- ^ The actual value
+                  -> Assertion
+assertApproxEqual preface epsilon expected actual =
+  unless (abs (actual - expected) <= epsilon) (assertFailure msg)
+  where msg = (if null preface then "" else preface ++ "\n") ++
+              "expected: " ++ show expected ++ "\n but got: " ++ show actual ++
+              "\n (maximum margin of error: " ++ show epsilon ++ ")"
+
+-- | Asserts that the specified actual value is approximately equal to the
+-- expected value (with the expected value on the right-hand side). The margin
+-- of error is specified with the implicit parameter @epsilon@.
+(@?~) :: (Ord a, Num a, Show a, ?epsilon :: a)
+      => a        -- ^ The actual value
+      -> a        -- ^ The expected value
+      -> Assertion
+x @?~ y = assertApproxEqual "" ?epsilon y x
+infix 1 @?~
+
+-- | Asserts that the specified actual value is approximately equal to the
+-- expected value (with the expected value on the left-hand side). The margin
+-- of error is specified with the implicit parameter @epsilon@.
+(@~?) :: (Ord a, Num a, Show a, ?epsilon :: a)
+      => a     -- ^ The expected value
+      -> a     -- ^ The actual value
+      -> Assertion
+x @~? y = assertApproxEqual "" ?epsilon x y
+infix 1 @~?
+
+-- | Shorthand for a test case that asserts approximate equality (with the
+-- expected value on the left-hand side, and the actual value on the
+-- right-hand side).
+(~~?) :: (Ord a, Num a, Show a, ?epsilon :: a)
+      => a     -- ^ The expected value
+      -> a     -- ^ The actual value
+      -> Test
+expected ~~? actual = TestCase (expected @~? actual)
+infix 1 ~~?
+
+-- | Shorthand for a test case that asserts approximate equality (with the
+-- actual value on the left-hand side, and the expected value on the
+-- right-hand side).
+(~?~) :: (Ord a, Num a, Show a, ?epsilon :: a)
+      => a     -- ^ The actual value
+      -> a     -- ^ The expected value 
+      -> Test
+actual ~?~ expected = TestCase (actual @?~ expected)
+infix 1 ~?~
diff --git a/test/Test/HUnit/Plus.hs b/test/Test/HUnit/Plus.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/HUnit/Plus.hs
@@ -0,0 +1,37 @@
+module Test.HUnit.Plus
+(
+    assertMaybeDouble
+  , assertOnFunction
+  , assertVector
+  , assertMatrix
+)
+where
+
+import Control.Monad
+import Test.HUnit
+import Test.HUnit.Approx
+import Numeric.LinearAlgebra
+
+assertMaybeDouble :: Maybe Double -> Maybe Double -> Double -> Assertion
+assertMaybeDouble Nothing Nothing _ = assertString ""
+assertMaybeDouble expected Nothing _ = assertString msg
+  where msg = "expected: " ++ show expected ++ "\nbut got: Nothing"
+assertMaybeDouble Nothing actual _ = assertString msg
+  where msg = "expected: Nothing\nbit got: " ++ show actual
+assertMaybeDouble (Just expected) (Just actual) eps = assertApproxEqual "Maybe Double" eps expected actual
+
+assertOnFunction :: (Eq b, Show b) => (a -> b) -> a -> a -> Assertion
+assertOnFunction func expected actual = func expected @=? func actual
+
+assertVector :: String -> R -> Vector R -> Vector R -> Assertion
+assertVector message eps expected actual =
+  let diff = norm_2 (expected-actual)
+      msg = message ++ "\nexpected: " ++ show expected ++ "\nbut got" ++ show actual
+  in unless (diff < eps) (assertFailure msg)
+
+assertMatrix :: String -> R -> Matrix R -> Matrix R -> Assertion
+assertMatrix message eps expected actual =
+  let diff = norm_2 (expected-actual)
+      msg = message ++ "\nexpected: " ++ show expected ++ "\nbut got" ++ show actual
+  in unless (diff < eps) (assertFailure msg)
+