packages feed

Synapse (empty) → 0.1.0.0

raw patch · 28 files changed

+3554/−0 lines, 28 filesdep +HUnitdep +Synapsedep +base

Dependencies added: HUnit, Synapse, base, easyplot, hashable, random, terminal-progress-bar, unordered-containers, vector

Files

+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Revision history for Synapse
+
+## 0.1.0.0 -- 2025-03-20
+
+* First version.
+ LICENSE view
@@ -0,0 +1,20 @@+Copyright (c) 2025 JktuJQ
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be included
+in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ README.md view
@@ -0,0 +1,380 @@+# ![Synapse logo](SynapseLogo.png) Synapse ![Synapse logo](SynapseLogo.png)
+
+`Synapse` is a machine learning library written in pure Haskell,
+that makes creating and training neural networks an easy job.
+
+## 🔨 Design 🔨
+
+Goals of `Synapse` library are to provide interface that is:
+
+* easily extensible
+* simple
+* powerful
+
+Haskell ecosystem only offers a few of machine learning libraries,
+but all of them have very complicated interface -
+[neural](https://hackage.haskell.org/package/neural) does not have too much extensibility and its typing is pretty hard,
+[grenade](https://hackage.haskell.org/package/grenade), although very powerful, uses a lot of type trickery that is impossible to reason about for beginners
+and [hasktorch](https://hackage.haskell.org/package/hasktorch) is a wrapper that is not as convenient as one might want.
+
+`Synapse` tries to resemble Python's [Keras](https://keras.io/api/) API,
+which is a unified interface for backends such as [Pytorch](https://pytorch.org/) and [Tensorflow](https://www.tensorflow.org/).
+
+Even though Python code practices usually are not what we want to see in Haskell code,
+mantaining the same level of accessibility and flexibility is what `Synapse` is focused on.
+
+## 💻 Usage 💻
+
+`Synapse` comes batteries-included with its own matrices,
+autodifferentiation system,
+and neural networks building blocks.
+
+### `Vec`s and `Mat`s
+
+Clever usage of Haskell's typeclasses allows operations to be written 'as is' for different types.
+
+Following example emulates dense layer forward propagation on plain matrices.
+
+```Haskell
+input   = M.rowVec $ V.fromList [1.0, 2.0, 3.0]
+
+weights = M.replicate (3, 3) 1.0
+bias    = M.rowVec $ V.replicate 3 0.0
+
+output  = tanh (input `matMul` weights `addMatRows` bias)
+```
+
+### Symbolic operations and autograd system
+
+`Synapse` autograd system implements reverse-mode dynamic automatic differentiation, where graph of operations is created on the fly.
+Its usage is much simpler than [ad](https://hackage.haskell.org/package/ad),
+and it is easily extensible - all you need is to define a function that produces `Symbol` which will hold gradients of
+that function and that is it!
+There is a [backprop](https://hackage.haskell.org/package/backprop) library that has a very similar implementation to `Synapse`'s one,
+so if you are familiar with [backprop](https://hackage.haskell.org/package/backprop) you will have no problems using `Synapse.Autograd` module.
+
+Speaking of the previous example - you might want to record gradients of those operations, and that is just as easy:
+
+```Haskell
+input   = symbol (SymbolIdentifier "input") $ M.rowVec $ V.fromList [1.0, 2.0, 3.0]
+
+weights = symbol (SymbolIdentifier "weights") $ M.replicate (3, 3) 1.0
+bias    = symbol (SymbolIdentifier "bias") $ M.rowVec $ V.replicate 3 0.0
+
+output  = tanh (input `matMul` weights `addMatRows` bias)
+```
+
+You just need to set the names for your symbolic variables, and you are good to go - `Synapse` will take care of the rest.
+
+Look at the `Synapse` implementation of some common operations:
+
+```Haskell
+(+) a b = symbolicBinaryOp (+) a b [(a, id), (b, id)]
+(*) a b = symbolicBinaryOp (*) a b [(a, (* b)), (b, (a *))]
+
+exp x = symbolicUnaryOp exp x [(x, (* exp x))]
+sin x = symbolicUnaryOp sin x [(x, (* cos x))]
+```
+
+Provided functions (`symbolicUnaryOp`, `symbolicBinaryOp`) expect an operation that will be performed on *values* of those symbols,
+symbols themselves and a list of tuples,
+where each tuple represents a gradient (the symbol which gradient is taken and function that implements chain rule - multiplies already calculated gradient (*symbol*) by the gradient of the function (another *symbol*)).
+Using those, defining new symbolic operations is very easy,
+and you should note that any composition of symbolic operations is itself a symbolic operation.
+
+This implementation is really easy to use:
+
+```Haskell
+a = symbol (SymbolIdentifier "a") 4
+b = symbol (SymbolIdentifier "b") 3
+c = renameSymbol (SymbolIdentifier "c") ((a * a) + (b * b))
+
+getGradientsOf (a * (b + b) * a) `wrt` a  -- == 4 * 4 * 3
+getGradientsOf (c * c) `wrt` c            -- == 2 * 25
+
+nthGradient 2 (a * (b + b) * a) a         -- == 4 * 3
+nthGradient 4 (sin c) c                   -- == sin c
+```
+
+`Synapse` does not care what is the type of your symbols: it might be `Int`, `Double`, `Vec`, `Mat`, anything that instantiates `Symbolic` typeclass -
+types just need to match with each other and with types of operations.
+
+### Neural networks
+
+`Synapse` takes as much of [Keras](https://keras.io/api/) API as it is possible, but is also provides additional abstractions leveraging Haskell's type system.
+
+#### Functions
+
+Everything that is a function in a common sense, is a function too in `Synapse`.
+
+`Activation`, `Loss`, `Metric`, `LearningRate`, `Constraint`, `Initializer`, `Regularizer` newtypes
+just wrap any plain Haskell function with needed type!
+
+That means that to create new activation function, loss function, etc. you just need to create Haskell function with appropriate constrains.
+
+```Haskell
+sigmoid :: (Symbolic a, Floating a) => a -> a
+sigmoid x = recip $ exp (negate x) + symbolicOne x
+
+sigmoidActivation :: (Symbolic a, Floating a) => Activation a
+sigmoidActivation = Activation sigmoid
+```
+
+`symbolicOne` function represents constant that corresponds to identity element for given `x`.
+Usage of const literal `1.0` does not work, because that identity element also needs to know `x`s dimensions.
+If `x` is matrix, you need to get `M.replicate (M.size x) 1.0`,
+not the singleton `1.0`.
+Writing additional constraints like `Symbolic` and having to create constants using `symbolicOne` might seem tedious,
+but that ensures type safety.
+
+You can also specialize the function:
+
+```Haskell
+type ActivationFn a = SymbolMat a -> SymbolMat a
+
+sigmoid :: (Symbolic a, Floating a) => ActivationFn a
+sigmoid x = recip $ exp (negate x) +. 1.0
+
+sigmoidActivation :: (Symbolic a, Floating a) => Activation a
+sigmoidActivation = Activation sigmoid
+```
+
+Even with all those limitations, it is still easy to create your own functions for any task.
+
+#### Layer system
+
+`AbstractLayer` typeclass is the most low-leveled abstraction of entire `Synapse` library.
+3 functions (`getParameters`, `updateParameters` and `symbolicForward`) are the backbone of entire neural networks interface.
+Docs on that typeclass as well as docs on those functions extensively describe invariants that `Synapse` expects from their implementations.
+
+With the help of `Layer` existential datatype `Synapse` is able to
+build networks from any types that are instances of `AbstractLayer` typeclass,
+which means that this system is easily extendable.
+
+`Dense` layer, for example, supports regularizers, constraints, recording gradients on its forward operations, and that is built upon those 3 functions.
+
+#### Training
+
+Here is the `train` function signature:
+
+```Haskell
+train
+    :: (Symbolic a, Floating a, Ord a, Show a, RandomGen g, AbstractLayer model, Optimizer optimizer)
+    => model a
+    -> optimizer a
+    -> Hyperparameters a
+    -> Callbacks model optimizer a
+    -> g
+    -> IO (model a, [OptimizerParameters optimizer a], Vec (RecordedMetric a), g)
+```
+
+Let's break it down into pieces and examine that function.
+
+#### Models
+
+`model` represents any `AbstractLayer` instance, so it is the model with parameters that are going to be trained.
+
+Any layer can be a model, but more commonly you would use `SequentialModel`.
+`SequentialModel` is a newtype that wraps list of `Layer`s.
+`buildSequentialModel` function builds the model, ensuring that dimensions of layers match.
+
+That is achieved by `LayerConfiguration` type alias and corresponding functions like `layerDense` and `layerActivation`.
+`LayerConfiguration` represents functions that are able to build a new layer upon the other layer, using information about its output dimension.
+
+```Haskell
+layers = [ Layer . layerDense 1
+         , Layer . layerActivation (Activation tanh)
+         , Layer . layerDense 1
+         ] :: [LayerConfiguration (Layer Double)]
+```
+
+You just write your layers like that and let `buildSequentialModel` to figure out how to compose them.
+
+It would look like this:
+
+```Haskell
+model = buildSequentialModel
+            (InputSize 1)
+            [ Layer . layerDense 1
+            , Layer . layerActivation (Activation tanh)
+            , Layer . layerDense 1
+            ] :: SequentialModel Double
+```
+
+`InputSize` indicates size of input that will be supported by this model.
+Model can take any matrix of size `(n, i)`, where `i` was supplied as `InputSize i` when the model was built.
+
+Since `AbstractLayer` instance is a trainable model and `SequentialModel` is a model, it means that it is also an instance `AbstractLayer`.
+Some layers are inherently a composition of other layers (LSTM layers are the example) and `Synapse` supports this automatically.
+
+#### Optimizers
+
+`optimizer` represents any `Optimizer` instance.
+Any optimizer has its parameters (`OptimizerParameters`) which it uses to update parameters of a model.
+
+Update is done with the functions `optimizerUpdateStep` and `optimizerUpdateParameters`.
+The second function is a mass update, so it needs gradients on all parameters of the model which are represented by symbolic matrices, while the first updates only one parameters which does not need to be symbolic, due to supplied exact gradient value.
+
+It is pretty easy to implement your own optimizer.
+See how `Synapse` implements `SGD`:
+
+```Haskell
+data SGD a = SGD
+    { sgdMomentum :: a
+    , sgdNesterov :: Bool
+    } deriving (Eq, Show)
+
+instance Optimizer SGD where
+    type OptimizerParameters SGD a = Mat a
+
+    optimizerInitialParameters _ parameter = zeroes (M.size parameter)
+
+    optimizerUpdateStep (SGD momentum nesterov) (lr, gradient) (parameter, velocity) = (parameter', velocity')
+      where
+        velocity' = velocity *. momentum - gradient *. lr
+
+        parameter' = if nesterov
+                     then parameter + velocity' *. momentum - gradient *. lr
+                     else parameter + velocity'
+```
+
+#### Hyperparameters
+
+Any training has some hyperparameters that configure that training.
+
+```Haskell
+data Hyperparameters a = Hyperparameters
+    { hyperparametersEpochs       :: Int
+    , hyperparametersBatchSize    :: Int
+
+    , hyperparametersDataset      :: VecDataset a
+
+    , hyperparametersLearningRate :: LearningRate a
+    , hyperparametersLoss         :: Loss a
+
+    , hyperparametersMetrics      :: Vec (Metric a)
+    }
+```
+
+Those hyperparameters include the number of epochs, batch size,
+dataset of vector samples (vector input and vector output),
+learning rate function, loss function
+and metrics that will be recorded during training.
+
+#### Callbacks
+
+`Synapse` allows 'peeking' in the training process using callbacks system.
+
+Several type aliases
+(`CallbackFnOnTrainBegin`,
+`CallbackFnOnEpochBegin`,
+`CallbackFnOnBatchBegin`,
+`CallbackFnOnBatchEnd`,
+`CallbackFnOnEpochEnd`,
+`CallbackFnOnTrainEnd`)
+represent functions that take mutable references to training parameters and do something with them (read, print/save, modify, etc.).
+
+Callback system interface should be used with caution, because some changes might break the training completely,
+but nonetheless, it is a very powerful instrument.
+
+#### Training process
+
+Training itself consists of following steps:
+
+1. Training beginning (setting up initial parameters of model and optimizer)
+2. Epoch training (shuffling, batching and processing of the dataset)
+3. Batch training (update of parameters based on the result of batch processing, recording of metrics)
+4. Training end (collecting results of training)
+
+All of that is handled by the `train` function.
+
+Here is an example of sine wave approximator which you could find at [tests](./test/) directory:
+
+```Haskell
+let sinFn x = (-3.0) * sin (x + 5.0)
+let model = buildSequentialModel (InputSize 1) [ Layer . layerDense 1
+                                               , Layer . layerActivation (Activation cos)
+                                               , Layer . layerDense 1
+                                               ] :: SequentialModel Double
+let dataset = Dataset $ V.fromList $ [Sample (singleton x) (sinFn $ singleton x) | x <- [-pi, -pi+0.2 .. pi]]
+(trainedModel, _, losses, _) <- train model
+                                   (SGD 0.2 False)
+                                   (Hyperparameters 500 16 dataset (LearningRate $ const 0.01) (Loss mse) V.empty)
+                                   emptyCallbacks
+                                   (mkStdGen 1)
+_ <- plot (PNG "test/plots/sin.png")
+          [ Data2D [Title "predicted sin", Style Lines, Color Red] [Range (-pi) pi] [(x, unSingleton $ forward (singleton x) trainedModel) | x <- [-pi, -pi+0.05 .. pi]]
+          , Data2D [Title "true sin", Style Lines, Color Green] [Range (-pi) pi] [(x, sinFn x) | x <- [-pi, -pi+0.05 .. pi]]
+          ]
+let unpackedLosses = unRecordedMetric (unsafeIndex losses 0)
+let lastLoss = unsafeIndex unpackedLosses (V.size unpackedLosses - 1)
+assertBool "trained well enough" (lastLoss < 0.01)
+```
+
+##### Prefix system
+
+`Synapse` manages gradients and parameters for layers with erased type information using prefix system.
+
+`SymbolIdentifier` is a prefix for name of symbolic parameters that are used in calculation.
+Every used parameter should have unique name to be recognised by the autograd -
+it must start with given prefix and end with the numerical index of said parameter.
+For example 3rd layer with 2 parameters (weights and bias) should
+name its weights symbol `"ml3w1"` and name its bias symbol `"ml3w2"` (`"ml3w"` prefix will be supplied).
+
+Prefix system along with layer system require to carefully ensure all the invariants that `Synapse` imposes if you are willing to extend them (write your own layers, training loops, etc.).
+But the user of this library should not worry about those getting in the way, because they are hidden behind an abstraction.
+
+## 📖 Future plans 📖
+
+`Synapse` library is still under development and there is work to be done on:
+
+* **Performance**
+
+    `Synapse` 'brings its own guns' and, although it makes the library independent,
+it could make `Synapse` to miss on some things that are battle-tested and tuned to performance.
+
+    That is especially true for `Synapse.Tensors` implementations of `Vec` and `Mat`.
+Those are built upon [vector](https://hackage.haskell.org/package/vector) library,
+which is good, but it is not suitable for heavy numerical calculations.
+
+    [hmatrix](https://hackage.haskell.org/package/hmatrix) which offers numerical computations based on BLAS and LAPACK is way more efficient.
+It would be great if `Synapse` library could work with any matrix backend.
+
+* **Tensors**
+
+    It is really a shame that `Synapse.Tensors` does not have actual tensors.
+`Tensor` datatype would allow to get rid of `Vec` and `Mat` datatypes in favour of more powerful abstraction.
+Tensor broadcasting could also be created to address issues that `Symbolic` typeclass is trying to solve.
+`Tensor` datatype could even present a unified frontend for any matrix backend that `Synapse` could use.
+
+* **GPU support**
+
+    This clause addresses all of the issues above.
+It would severely increase performance of `Synapse` library,
+and it would greatly work with backend-independent tensors.
+Haskell ecosystem offers great [accelerate](https://hackage.haskell.org/package/accelerate) library which could help with all those problems.
+
+* **More out-of-the-box solutions**
+
+    At this point, `Synapse` does not offer a wide variety of layers, activations, models, optimizers out-of-the-box.
+Suppling more of them would definitely help.
+
+* **Developer abstractions**
+
+    Currently, `Synapse`'s backbones are `SymbolIdentifier`s and implementations of `AbstractLayer`,
+which might be cumbersome for developers and be a bit unreliable.
+If those systems could work naturally, without some hardcoding of values, it would be great.
+
+* **More monads!**
+
+    `Synapse` library does not use a lot of Haskell instruments, like optics, monad transformers, etc.
+Although it makes the library easy for beginners,
+I am sure that some of those instruments can offer richer expressiveness for the code
+and also address the issue of 'Developer abstractions'.
+
+## ❤️ Contributions ❤️
+
+`Synapse` library would benefit from every contribution to it: docs, code, small advertisement - you name it.
+
+If you want to help me in the development, you could always contact me -
+my [GitHub profile](https://github.com/JktuJQ) has links to my socials.
+ Synapse.cabal view
@@ -0,0 +1,110 @@+cabal-version:      3.0
+name:               Synapse
+version:            0.1.0.0
+build-type:         Simple
+
+
+author:             JktuJQ
+maintainer:         odmamontov@gmail.com
+
+
+homepage:           https://github.com/JktuJQ/Synapse
+extra-doc-files:    README.md
+                  , CHANGELOG.md
+
+license:            MIT
+license-file:       LICENSE
+
+synopsis:
+    Synapse is a machine learning library written in pure Haskell.
+description:
+    Synapse is a machine learning library written in pure Haskell, that makes creating and training neural networks an easy job.
+category:           Math
+ 
+
+source-repository head
+  type:     git
+  location: https://github.com/JktuJQ/Synapse
+  subdir:   
+
+
+common warnings
+    ghc-options: -Wall
+
+library
+    import:           warnings
+    default-language: Haskell2010
+
+    build-depends:    base ^>= 4.17.2.1
+
+                    , vector >= 0.13.2 && < 0.14
+
+                    , hashable >= 1.4.7 && < 1.5
+                    , unordered-containers >= 0.2.20 && < 0.3
+
+                    , random >= 1.3.0 && < 1.4
+
+                    , terminal-progress-bar >= 0.4.2 && < 0.5
+
+    hs-source-dirs:   src
+    exposed-modules:  Synapse
+                    
+                      -- Tensors
+                    , Synapse.Tensors
+                    , Synapse.Tensors.Mat
+                    , Synapse.Tensors.Vec
+
+
+                      -- Autograd
+                    , Synapse.Autograd
+                    
+
+                      -- NN
+                    , Synapse.NN
+
+                    , Synapse.NN.Layers
+                    , Synapse.NN.Layers.Layer
+
+                    , Synapse.NN.Layers.Initializers
+                    , Synapse.NN.Layers.Constraints
+                    , Synapse.NN.Layers.Regularizers
+
+                    , Synapse.NN.Layers.Activations
+                    , Synapse.NN.Layers.Dense
+
+                    , Synapse.NN.Models
+
+                    , Synapse.NN.Training
+
+                    , Synapse.NN.Batching
+                    
+                    , Synapse.NN.Losses
+                    , Synapse.NN.Metrics
+
+                    , Synapse.NN.LearningRates
+
+                    , Synapse.NN.Optimizers
+                    
+    -- other-modules:
+    -- other-extensions:
+
+test-suite Synapse-test
+    import:           warnings
+    default-language: Haskell2010
+    type:             exitcode-stdio-1.0
+
+    build-depends:    base ^>=4.17.2.1
+                    , Synapse
+
+                    , random >= 1.3.0 && < 1.4
+
+                    , HUnit >= 1.6.2 && < 1.7
+
+                    , easyplot >= 1.0 && < 1.1
+
+    hs-source-dirs:   test
+    main-is:          Main.hs
+    other-modules:    TensorsTest
+                    , AutogradTest
+                    , NNTest
+    -- other-extensions:
+ src/Synapse.hs view
@@ -0,0 +1,11 @@+{- | Synapse
+
+"Synapse" is a machine learning library written in pure Haskell,
+that makes creating and training neural networks an easy job.
+
+There are examples of "Synapse" usage at <https://github.com/JktuJQ/Synapse>,
+you are encouraged to look at library's repository to get yourself
+a grasp on basic concepts.
+-}
+
+module Synapse () where
+ src/Synapse/Autograd.hs view
@@ -0,0 +1,347 @@+{- | This module implements reverse-mode automatic differentiation.
+
+Machine learning and training of models are based on calculating gradients of operations.
+This can be done symbolically by dynamically creating a graph of all operations,
+which is then traversed to obtain the gradient.
+
+"Synapse" provides several operations that support automatic differentiation,
+but you could easily extend list of those: you just need to define function
+that returns 'Symbol' with correct local gradients.
+You can check out implementations in the source to give yourself a reference
+and read more about it in 'Symbol' datatype docs.
+-}
+
+
+-- 'FlexibleInstances' and 'TypeFamilies' are needed to instantiate 'Indexable', 'ElementwiseScalarOps', 'SingletonOps', 'VecOps', 'MatOps' typeclasses.
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeFamilies      #-}
+
+
+module Synapse.Autograd
+    ( -- * 'Symbolic' and 'Symbol'
+
+      Symbolic (symbolicZero, symbolicOne, symbolicN)
+
+    , SymbolIdentifier (SymbolIdentifier, unSymbolIdentifier)
+
+    , Symbol (Symbol, symbolIdentifier, unSymbol, symbolGradients)
+    , SymbolVec
+    , SymbolMat
+
+    , symbol
+    , constSymbol
+    , renameSymbol
+
+    , symbolicUnaryOp
+    , symbolicBinaryOp
+
+      -- * 'Gradients' calculation
+
+    , Gradients (unGradients)
+    , allGradients
+    , getGradientsOf
+    , wrt
+    , nthPartialGradient
+    , nthGradient
+    ) where
+
+
+import Synapse.Tensors (DType, Indexable(..), ElementwiseScalarOps(..), ElementwiseScalarOps(..), SingletonOps(..), VecOps(..), MatOps(..))
+
+import Synapse.Tensors.Vec (Vec)
+import qualified Synapse.Tensors.Vec as V
+
+import Synapse.Tensors.Mat (Mat)
+import qualified Synapse.Tensors.Mat as M
+
+import Data.Foldable (foldl')
+import Data.String (IsString(..))
+
+import Data.Hashable (Hashable(..))
+
+import qualified Data.HashMap.Lazy as HM
+
+
+{- | 'Symbolic' typeclass describes types with few properties that are needed for autogradient.
+
+Members of this typeclass could have default implementation due to 'Num', but such implementation is not always correct.
+'Synapse.Tensors.Vec.Vec's and 'Synapse.Tensors.Mat.Mat's do not have only one zero or identity element, and so numerical literal is not enough.
+'symbolicZero' and 'symbolicOne' function additionally take reference value to consider dimensions.
+Absence of default implementations forces to manually ensure correctness of those functions.
+
+"Synapse" provides implementations for primitive types ('Int', 'Float', 'Double'),
+and for containers types ('Synapse.Tensors.Vec.Vec', 'Synapse.Tensors.Vec.Vec').
+
+Detailed laws of 'Symbolic' properties are in the docs for associated functions.
+-}
+class (Eq a, Num a) => Symbolic a where
+    -- | Returns additive and multiplicative (elementwise) zero element. Argument is passed for the reference of the dimension.
+    symbolicZero :: a -> a
+
+    -- | Returns multiplicative (elementwise) identity element. Argument is passed for the reference of the dimension.
+    symbolicOne :: a -> a
+
+    -- | Returns what could be considered @N@ constant (sum of @N@ 'symbolicOne's). Argument is passed for the reference of the dimension.
+    symbolicN :: Int -> a -> a
+    symbolicN n c
+        | n < 0     = negate $ go (abs n) c
+        | n == 0    = symbolicZero c
+        | otherwise = go n c
+      where
+        go 0 x = x
+        go i x = go (i - 1) (x + symbolicOne x)
+
+-- Instances
+
+instance Symbolic Int where
+    symbolicZero = const 0
+    symbolicOne = const 1
+    symbolicN = const
+
+instance Symbolic Float where
+    symbolicZero = const 0.0
+    symbolicOne = const 1.0
+    symbolicN n = const (fromIntegral n)
+
+instance Symbolic Double where
+    symbolicZero = const 0.0
+    symbolicOne = const 1.0
+    symbolicN n = const (fromIntegral n)
+
+
+instance Symbolic a => Symbolic (Vec a) where
+    symbolicZero reference = V.replicate (V.size reference) 0
+    symbolicOne reference = V.replicate (V.size reference) 1
+    symbolicN n reference = V.replicate (V.size reference) (fromIntegral n)
+
+instance Symbolic a => Symbolic (Mat a) where
+    symbolicZero reference = M.replicate (M.size reference) 0
+    symbolicOne reference = M.replicate (M.size reference) 1
+    symbolicN n reference = M.replicate (M.size reference) (fromIntegral n)
+
+
+-- | 'SymbolIdentifier' is a newtype that wraps string, which needs to uniquely represent symbol.
+newtype SymbolIdentifier = SymbolIdentifier 
+    { unSymbolIdentifier :: String  -- ^ Identifier of a symbol.
+    } deriving (Eq, Show)
+
+instance IsString SymbolIdentifier where
+    fromString = SymbolIdentifier
+
+instance Semigroup SymbolIdentifier where
+    (<>) (SymbolIdentifier a) (SymbolIdentifier b) = SymbolIdentifier $ a <> b
+
+instance Monoid SymbolIdentifier where
+    mempty = SymbolIdentifier ""
+
+{- | Datatype that represents symbol variable (variable which operations are recorded to symbolically obtain derivatives).
+
+Any operation returning @Symbol a@ where @a@ is 'Symbolic' could be autogradiented - returned 'Symbol' has 'symbolGradients' list,
+which allows "Synapse" to build a graph of computation and obtain needed gradients.
+'symbolGradients' list contains pairs: first element in that pair is symbol wrt which you can take gradient and
+the second element is closure that represents chain rule - it takes incoming local gradient of said symbol and multiplies it by local derivative.
+You can check out implementations of those operations in the source to give yourself a reference.
+-}
+data Symbol a = Symbol
+    { symbolIdentifier :: SymbolIdentifier                    -- ^ Name of a symbol (identifier for differentiation).
+    , unSymbol         :: a                                   -- ^ Value of a symbol.
+    , symbolGradients  :: [(Symbol a, Symbol a -> Symbol a)]  -- ^ List of gradients (wrt to what Symbol and closure to calculate gradient). 
+    }
+
+-- | Creates new symbol that refers to a variable (so it must have a name to be able to be differentiated wrt).
+symbol :: SymbolIdentifier -> a -> Symbol a
+symbol name value = Symbol name value []
+
+-- | Creates new symbol that refers to constant (so it does not have name and thus its gradients are not saved).
+constSymbol :: a -> Symbol a
+constSymbol = symbol mempty
+
+-- | Renames symbol which allows differentiating wrt it. Note: renaming practically creates new symbol for the gradient calculation.
+renameSymbol :: SymbolIdentifier -> Symbol a -> Symbol a
+renameSymbol name (Symbol _ value localGradients) = Symbol name value localGradients
+
+
+-- | @SymbolVec a@ type alias stands for @Symbol (Vec a)@.
+type SymbolVec a = Symbol (Vec a)
+
+-- | @SymbolMat a@ type alias stands for @Symbol (Mat a)@.
+type SymbolMat a = Symbol (Mat a)
+
+
+-- Typeclasses
+
+instance Show a => Show (Symbol a) where
+    show (Symbol name value _) = "Symbol " ++ show name ++ ": " ++ show value
+
+
+instance Eq (Symbol a) where
+    (==) (Symbol name1 _ _) (Symbol name2 _ _) = name1 == name2
+
+
+instance Hashable (Symbol a) where
+    hashWithSalt salt (Symbol (SymbolIdentifier name) _ _) = hashWithSalt salt name
+
+
+-- Symbolic symbols
+
+instance Symbolic a => Symbolic (Symbol a) where
+    symbolicZero x = constSymbol $ symbolicZero $ unSymbol x
+    symbolicOne x = constSymbol $ symbolicOne $ unSymbol x
+
+
+type instance DType (SymbolVec a) = a
+
+type instance DType (SymbolMat a) = a
+
+
+-- | Converts unary operation into symbolic one.
+symbolicUnaryOp :: (a -> a) -> Symbol a -> [(Symbol a, Symbol a -> Symbol a)] -> Symbol a
+symbolicUnaryOp op x = Symbol mempty (op (unSymbol x))
+
+-- | Converts binary operation into symbolic one.
+symbolicBinaryOp :: (a -> a -> a) -> Symbol a -> Symbol a -> [(Symbol a, Symbol a -> Symbol a)] -> Symbol a
+symbolicBinaryOp op a b = Symbol mempty (op (unSymbol a) (unSymbol b))
+
+instance Symbolic a => Num (Symbol a) where
+    (+) a b = symbolicBinaryOp (+) a b [(a, id), (b, id)]
+    (-) a b = symbolicBinaryOp (-) a b [(a, id), (b, negate)]
+    negate x = symbolicUnaryOp negate x [(x, negate)]
+    (*) a b = symbolicBinaryOp (*) a b [(a, (* b)), (b, (a *))]
+    abs x = symbolicUnaryOp abs x [(x, signum)]
+    signum x = symbolicUnaryOp signum x [(x, id)]
+    fromInteger = constSymbol . fromInteger
+
+instance (Symbolic a, Fractional a) => Fractional (Symbol a) where
+    (/) a b = symbolicBinaryOp (/) a b [(a, (/ b)), (b, (* (negate a / (b * b))))]
+    recip x = symbolicUnaryOp recip x [(x, (* (negate (symbolicOne x) / (x * x))))]
+    fromRational = constSymbol . fromRational
+
+instance (Symbolic a, Floating a) => Floating (Symbol a) where
+    pi = constSymbol pi
+    (**) a b = symbolicBinaryOp (**) a b [(a, (* (b * a ** (b - symbolicOne b)))), (b, (* (a ** b * log a)))]
+    sqrt x = symbolicUnaryOp sqrt x [(x, (* (recip $ symbolicN 2 x * sqrt x)))]
+    exp x = symbolicUnaryOp exp x [(x, (* exp x))]
+    log x = symbolicUnaryOp log x [(x, (* recip x))]
+    sin x = symbolicUnaryOp sin x [(x, (* cos x))]
+    cos x = symbolicUnaryOp cos x [(x, (* negate (sin x)))]
+    asin x = symbolicUnaryOp asin x [(x, (* recip (sqrt (symbolicOne x - x * x))))]
+    acos x = symbolicUnaryOp acos x [(x, (* negate (recip (sqrt (symbolicOne x - x * x)))))]
+    atan x = symbolicUnaryOp atan x [(x, (* recip (symbolicOne x + x * x)))]
+    sinh x = symbolicUnaryOp sinh x [(x, (* cosh x))]
+    cosh x = symbolicUnaryOp cosh x [(x, (* sinh x))]
+    asinh x = symbolicUnaryOp asinh x [(x, (* recip (sqrt (symbolicOne x + x * x))))]
+    acosh x = symbolicUnaryOp acosh x [(x, (* recip (sqrt (x * x - symbolicOne x))))]
+    atanh x = symbolicUnaryOp atanh x [(x, (* recip (symbolicOne x - x * x)))]
+
+instance Symbolic a => ElementwiseScalarOps (Symbol (Vec a)) where
+    (+.) x n = x + constSymbol (V.replicate (V.size $ unSymbol x) n)
+    (-.) x n = x - constSymbol (V.replicate (V.size $ unSymbol x) n)
+    (*.) x n = x * constSymbol (V.replicate (V.size $ unSymbol x) n)
+    (/.) x n = x / constSymbol (V.replicate (V.size $ unSymbol x) n)
+    (**.) x n = x ** constSymbol (V.replicate (V.size $ unSymbol x) n)
+
+    elementsMin x n = symbolicUnaryOp (`elementsMin` n) x
+                      [(x, (* constSymbol (V.generate (V.size $ unSymbol x) $ \i -> if unsafeIndex (unSymbol x) i <= n then 1 else 0)))]
+    elementsMax x n = symbolicUnaryOp (`elementsMax` n) x
+                      [(x, (* constSymbol (V.generate (V.size $ unSymbol x) $ \i -> if unsafeIndex (unSymbol x) i >= n then 1 else 0)))]
+
+instance Symbolic a => ElementwiseScalarOps (SymbolMat a) where
+    (+.) x n = x + constSymbol (M.replicate (M.size $ unSymbol x) n)
+    (-.) x n = x - constSymbol (M.replicate (M.size $ unSymbol x) n)
+    (*.) x n = x * constSymbol (M.replicate (M.size $ unSymbol x) n)
+    (/.) x n = x / constSymbol (M.replicate (M.size $ unSymbol x) n)
+    (**.) x n = x ** constSymbol (M.replicate (M.size $ unSymbol x) n)
+
+    elementsMin x n = symbolicUnaryOp (`elementsMin` n) x
+                      [(x, (* constSymbol (M.generate (M.size $ unSymbol x) $ \i -> if unsafeIndex (unSymbol x) i <= n then 1 else 0)))]
+    elementsMax x n = symbolicUnaryOp (`elementsMax` n) x
+                      [(x, (* constSymbol (M.generate (M.size $ unSymbol x) $ \i -> if unsafeIndex (unSymbol x) i >= n then 1 else 0)))]
+
+instance Symbolic a => SingletonOps (SymbolVec a) where
+    singleton = constSymbol . singleton
+    isSingleton = isSingleton . unSymbol
+    unSingleton = unSingleton . unSymbol
+
+    extendSingleton vec reference = constSymbol $ extendSingleton (unSymbol vec) (unSymbol reference)
+
+    elementsSum x = symbolicUnaryOp elementsSum x [(x, (`extendSingleton` x))]
+    elementsProduct x = let innerProduct = unSingleton $ elementsProduct $ unSymbol x
+                        in symbolicUnaryOp elementsProduct x [(x, \path -> extendSingleton path x * constSymbol (V.generate (V.size $ unSymbol x) $ \i -> innerProduct / unsafeIndex (unSymbol x) i))]
+
+    mean x = symbolicUnaryOp mean x [(x, \path -> extendSingleton path x /. fromIntegral (V.size (unSymbol x)))]
+
+    norm x = symbolicUnaryOp norm x [(x, \path -> extendSingleton path x * (x /. unSingleton (norm $ unSymbol x)))]
+
+instance Symbolic a => SingletonOps (SymbolMat a) where
+    singleton = constSymbol . singleton
+    isSingleton = isSingleton . unSymbol
+    unSingleton = unSingleton . unSymbol
+
+    extendSingleton mat reference = constSymbol $ extendSingleton (unSymbol mat) (unSymbol reference)
+
+    elementsSum x = symbolicUnaryOp elementsSum x [(x, (`extendSingleton` x))]
+    elementsProduct x = let innerProduct = unSingleton $ elementsProduct $ unSymbol x
+                        in symbolicUnaryOp elementsProduct x [(x, \path -> extendSingleton path x * constSymbol (M.generate (M.size $ unSymbol x) $ \i -> innerProduct / unsafeIndex (unSymbol x) i))]
+
+    mean x = symbolicUnaryOp mean x [(x, \path -> extendSingleton path x /. fromIntegral (M.nElements $ unSymbol x))]
+
+    norm x = symbolicUnaryOp norm x [(x, \path -> extendSingleton path x * (x /. unSingleton (norm $ unSymbol x)))]
+
+instance Symbolic a => VecOps (SymbolVec a) where
+    dot a b = elementsSum $ a * b
+
+instance Symbolic a => MatOps (SymbolMat a) where
+    transpose x = symbolicUnaryOp M.transpose x [(x, (* transpose x))]
+
+    addMatRow mat row = symbolicBinaryOp addMatRow mat row [(mat, id), (row, constSymbol . M.rowVec . flip M.indexRow 0 . unSymbol)]
+
+    matMul a b = symbolicBinaryOp M.matMul a b [(a, (`matMul` transpose b)), (b, (transpose a `matMul`))]
+
+
+-- Gradients calculation
+
+-- | 'Gradients' datatype holds all gradients of one symbol with respect to other symbols.
+newtype Gradients a = Gradients
+    { unGradients :: HM.HashMap (Symbol a) (Symbol a)  -- ^ Map of gradients.
+    }
+
+-- | Returns key-value pairs of all gradients of symbol.
+allGradients :: Gradients a -> [(Symbol a, Symbol a)]
+allGradients = HM.toList . unGradients
+
+
+-- Typeclasses
+
+instance Show a => Show (Gradients a) where
+    show gradients = show $ allGradients gradients
+
+
+-- | Generates 'Gradients' for given symbol.
+getGradientsOf :: Symbolic a => Symbol a -> Gradients a
+getGradientsOf differentiatedSymbol = Gradients $ HM.insert differentiatedSymbol wrtItself $
+                                                  HM.delete (Symbol mempty undefined []) $
+                                                  go HM.empty differentiatedSymbol wrtItself
+  where
+    wrtItself = symbolicOne differentiatedSymbol
+
+    go grads s pathValue =
+        foldr (\(child, mulPath) grad -> let pathValue' = mulPath pathValue
+                                             grad' = HM.alter (\e -> Just $ case e of
+                                                                                Nothing -> pathValue'
+                                                                                Just x  -> x + pathValue'
+                                                              ) child grad
+                                         in go grad' child pathValue') grads (symbolGradients s)
+
+-- | Chooses gradient with respect to given symbol.
+wrt :: Symbolic a => Gradients a -> Symbol a -> Symbol a
+wrt gradients x = HM.findWithDefault (symbolicZero x) x $ unGradients gradients
+
+-- | Takes partial gradients wrt to all symbols in a list sequentially, returning last result.
+nthPartialGradient :: Symbolic a => Symbol a -> [Symbol a] -> Symbol a
+nthPartialGradient = foldl' $ \y x -> getGradientsOf y `wrt` x
+
+-- | Takes nth order gradient of one symbol wrt other symbol. If n is negative number, an error is returned.
+nthGradient :: Symbolic a => Int -> Symbol a -> Symbol a -> Symbol a
+nthGradient n y x
+    | n < 0 = error "Cannot take negative order gradient"
+    | otherwise = nthPartialGradient y (replicate n x)
+ src/Synapse/NN.hs view
@@ -0,0 +1,37 @@+{- | This module provides functions and datatypes that are needed in work with neural networks.
+
+"Synapse" tries to supply as much of the common interface of any neural network library,
+such as the collection of different models, layers, optimisers, training algorithms and etc.
+-}
+
+
+module Synapse.NN 
+    ( -- * Re-exports
+
+      module Synapse.NN.Layers
+
+    , module Synapse.NN.Models
+    , module Synapse.NN.Training
+    , module Synapse.NN.Batching
+
+    , module Synapse.NN.Losses
+    , module Synapse.NN.Metrics
+
+    , module Synapse.NN.LearningRates
+    
+    , module Synapse.NN.Optimizers
+    ) where
+
+
+import Synapse.NN.Layers
+
+import Synapse.NN.Models
+import Synapse.NN.Training
+import Synapse.NN.Batching
+
+import Synapse.NN.Losses
+import Synapse.NN.Metrics
+
+import Synapse.NN.LearningRates
+
+import Synapse.NN.Optimizers
+ src/Synapse/NN/Batching.hs view
@@ -0,0 +1,104 @@+{- | Implements batching - technology that allows packing and processing multiple samples at once.
+-}
+
+
+-- 'TypeFamilies' are needed to instantiate 'DType'.
+{-# LANGUAGE TypeFamilies #-}
+
+
+module Synapse.NN.Batching
+    ( -- * 'Sample' datatype
+    
+      Sample (Sample, sampleInput, sampleOutput)
+
+      -- * 'Dataset' datatype
+
+    , Dataset (Dataset, unDataset)
+
+    , datasetSize
+
+    , shuffleDataset
+    , splitDataset
+
+    , VecDataset
+    , BatchedDataset
+    , batchVectors
+    ) where
+
+
+import Synapse.Tensors (DType, Indexable(unsafeIndex))
+
+import Synapse.Tensors.Vec (Vec(Vec))
+import qualified Synapse.Tensors.Vec as V
+
+import Synapse.Tensors.Mat (Mat)
+import qualified Synapse.Tensors.Mat as M
+
+import Control.Monad.ST (runST)
+
+import System.Random (RandomGen, uniformR)
+
+import Data.Vector (thaw, unsafeFreeze)
+import Data.Vector.Mutable (swap)
+
+
+-- | 'Sample' datatype represents known pair of inputs and outputs of function that is unknown.
+data Sample a = Sample
+    { sampleInput  :: a  -- ^ Sample input.
+    , sampleOutput :: a  -- ^ Sample output.
+    } deriving (Eq, Show)
+
+type instance DType (Sample a) = DType a
+
+
+-- | 'Dataset' newtype wraps vector of 'Sample's - it represents known information about unknown function.
+newtype Dataset a = Dataset 
+    { unDataset :: Vec (Sample a)  -- ^ Unwraps 'Dataset' newtype.
+    } deriving (Eq, Show)
+
+type instance DType (Dataset a) = DType a
+
+-- | Returns size of dataset.
+datasetSize :: Dataset a -> Int
+datasetSize = V.size . unDataset
+
+
+-- | Shuffles any 'Dataset' using Fisher-Yates algorithm.
+shuffleDataset :: RandomGen g => Dataset a -> g -> (Dataset a, g)
+shuffleDataset (Dataset dataset) gen
+    | V.size dataset <= 1 = (Dataset dataset, gen)
+    | otherwise           = runST $ do
+        mutableVector <- thaw $ V.unVec dataset
+        gen' <- go mutableVector (V.size dataset - 1) gen
+        shuffledVector <- unsafeFreeze mutableVector
+        return (Dataset $ Vec shuffledVector, gen')
+  where
+    go _ 0 seed = return seed
+    go v lastIndex seed = let (swapIndex, seed') = uniformR (0, lastIndex) seed
+                          in swap v swapIndex lastIndex >> go v (lastIndex - 1) seed'
+
+-- | Splits dataset such that size of left dataset divided on size of right dataset will be equal to given ratio.
+splitDataset :: Dataset a -> Float -> (Dataset a, Dataset a)
+splitDataset (Dataset dataset) ratio = let (left, right) = V.splitAt (round $ fromIntegral (V.size dataset) * ratio) dataset
+                                       in (Dataset left, Dataset right)
+
+
+-- | 'VecDataset' type alias represents 'Dataset's with samples of vector functions.
+type VecDataset a = Dataset (Vec a)
+-- | 'BatchedDataset' type alias represents 'Dataset's with samples of vector functions where multiple samples were batched together.
+type BatchedDataset a = Dataset (Mat a)
+
+-- | Batches 'VecDataset' by grouping a given amount of samples into batches.
+batchVectors :: Int -> VecDataset a -> BatchedDataset a
+batchVectors batchSize (Dataset dataset) = Dataset $ V.fromList $ map groupBatch $ split dataset
+  where
+    split vector
+        | V.size vector <= batchSize = [vector]
+        | otherwise                  = let (current, remainder) = V.splitAt batchSize vector
+                                       in current : split remainder
+    
+    groupBatch vector = let (rows, inputCols) = (V.size vector, V.size $ sampleInput $ unsafeIndex vector 0)
+                            group (r, c) = unsafeIndex ((if c < inputCols then sampleInput else sampleOutput) (unsafeIndex vector r)) (c `mod` inputCols)
+                            fullBatch = M.generate (rows, inputCols + V.size (sampleOutput $ unsafeIndex vector 0)) group
+                            (batchInput, batchOutput, _, _) = M.split fullBatch (rows, inputCols)
+                        in Sample batchInput batchOutput
+ src/Synapse/NN/Layers.hs view
@@ -0,0 +1,26 @@+{- | Provides several types of layers and functions to control their construction and usage.
+-}
+
+
+module Synapse.NN.Layers
+    ( -- * Re-exports
+
+      module Synapse.NN.Layers.Layer
+    
+    , module Synapse.NN.Layers.Initializers
+    , module Synapse.NN.Layers.Constraints
+    , module Synapse.NN.Layers.Regularizers
+
+    , module Synapse.NN.Layers.Activations
+    , module Synapse.NN.Layers.Dense
+    ) where
+
+
+import Synapse.NN.Layers.Layer
+
+import Synapse.NN.Layers.Initializers
+import Synapse.NN.Layers.Constraints
+import Synapse.NN.Layers.Regularizers
+
+import Synapse.NN.Layers.Activations
+import Synapse.NN.Layers.Dense
+ src/Synapse/NN/Layers/Activations.hs view
@@ -0,0 +1,77 @@+{- | Provides activation functions - unary functions that are differentiable almost everywhere and so they can be used in backward loss propagation.
+-}
+
+
+module Synapse.NN.Layers.Activations
+    ( -- * 'ActivationFn' type alias and 'Activation' newtype
+
+      ActivationFn
+    , activateScalar
+    , activateMat
+
+    , Activation (Activation, unActivation)
+    , layerActivation
+
+      -- * Activation functions
+
+    , relu
+    , sigmoid
+    ) where
+
+
+import Synapse.NN.Layers.Layer (AbstractLayer(..), LayerConfiguration)
+
+import Synapse.Tensors (ElementwiseScalarOps((+.), (/.)), SingletonOps(unSingleton))
+
+import Synapse.Tensors.Mat (Mat)
+import qualified Synapse.Tensors.Mat as M
+
+import Synapse.Autograd (Symbol(unSymbol), SymbolMat, Symbolic, constSymbol)
+
+
+-- | 'ActivationFn' is a type alias that represents unary functions that differentiable almost everywhere.
+type ActivationFn a = SymbolMat a -> SymbolMat a
+
+
+-- | Applies activation function to a scalar to produce new scalar.
+activateScalar :: Symbolic a => ActivationFn a -> a -> a
+activateScalar fn = unSingleton . unSymbol . fn . constSymbol . M.singleton
+
+-- | Applies activation function to a scalar to produce new scalar.
+activateMat :: Symbolic a => ActivationFn a -> Mat a -> Mat a
+activateMat fn = unSymbol . fn . constSymbol
+
+
+{- | 'Activation' newtype wraps 'ActivationFn's - unary functions that can be thought of as activation functions for neural network layers.
+
+Any activation function must be differentiable almost everywhere and so
+it must be function that operates on 'Synapse.Autograd.Symbol's, which is allows for function to be differentiated when needed.
+-}
+newtype Activation a = Activation
+    { unActivation :: ActivationFn a  -- ^ Unwraps 'Activation' newtype.
+    }
+
+instance AbstractLayer Activation where
+    inputSize _ = Nothing
+    outputSize _ = Nothing
+
+    nParameters _ = 0
+    getParameters _ _ = []
+    updateParameters = const
+
+    symbolicForward _ input (Activation fn) = (fn input, M.singleton 0)
+
+-- | Creates configuration for activation layer.
+layerActivation :: Activation a -> LayerConfiguration (Activation a)
+layerActivation = const
+
+
+-- Activation functions
+
+-- | ReLU function.
+relu :: (Symbolic a, Fractional a) => ActivationFn a
+relu x = (x + abs x) /. 2.0
+
+-- | Sigmoid function.
+sigmoid :: (Symbolic a, Floating a) => ActivationFn a
+sigmoid x = recip $ exp (negate x) +. 1.0
+ src/Synapse/NN/Layers/Constraints.hs view
@@ -0,0 +1,70 @@+{- | Allows to constraint values of layers parameters.
+
+'ConstraintFn' type alias represents functions that are able to constrain the values of matrix
+and 'Constraint' newtype wraps 'ConstraintFn's.
+
+'ConstraintFn's should be applied on matrices from the 'Synapse.NN.Layers.Layer.updateParameters' function.
+-}
+
+
+module Synapse.NN.Layers.Constraints
+    ( -- * 'ConstraintFn' type alias and 'Constraint' newtype
+      
+      ConstraintFn
+
+    , Constraint (Constraint, unConstraint)
+
+      -- * Value constraints
+    
+    , nonNegative
+    , clampMin
+    , clampMax
+    , clampMinMax
+
+      -- * Matrix-specific constraints
+    
+    , centralize
+    ) where
+
+
+import Synapse.Tensors (ElementwiseScalarOps((+.), (-.)), SingletonOps(unSingleton, mean))
+
+import Synapse.Tensors.Mat (Mat)
+
+import Data.Ord (clamp)
+
+
+-- | 'ConstraintFn' type alias represents functions that are able to constrain the values of matrix.
+type ConstraintFn a = Mat a -> Mat a
+
+
+-- | 'Constraint' newtype wraps 'ConstraintFn's - functions that are able to constrain the values of matrix.
+newtype Constraint a = Constraint
+    { unConstraint :: ConstraintFn a  -- ^ Unwraps 'Constraint' newtype.
+    }
+
+
+-- Value constraints
+
+-- | Ensures that matrix values are non-negative.
+nonNegative :: (Num a, Ord a) => ConstraintFn a
+nonNegative = fmap (max 0)
+
+-- | Ensures that matrix values are more or equal than given value.
+clampMin :: Ord a => a -> ConstraintFn a
+clampMin minimal = fmap (max minimal)
+
+-- | Ensures that matrix values are less or equal than given value.
+clampMax :: Ord a => a -> ConstraintFn a
+clampMax maximal = fmap (min maximal)
+
+-- | Ensures that matrix values are clamped between given values.
+clampMinMax :: Ord a => (a, a) -> ConstraintFn a
+clampMinMax = fmap . clamp
+
+
+-- Matrix-specific constraints
+
+-- | Ensures that matrix values are centralized by mean around given value.
+centralize :: Fractional a => a -> ConstraintFn a
+centralize center mat = mat -. unSingleton (mean mat) +. center
+ src/Synapse/NN/Layers/Dense.hs view
@@ -0,0 +1,94 @@+{- | Provides dense layer implementation.
+
+'Dense' datatype represents densely-connected neural network layer and
+it performs following operation: @x `matMul` w + b@, where @w@ is weights and @b@ is bias (if present) of a layer.
+-}
+
+
+-- 'TypeFamilies' are needed to instantiate 'DType'.
+{-# LANGUAGE TypeFamilies #-}
+
+
+module Synapse.NN.Layers.Dense
+    ( -- * 'Dense' datatype
+
+      Dense (Dense, denseWeights, denseBias, denseConstraints, denseRegularizers)
+    , layerDenseWith
+    , layerDense
+    ) where
+
+
+
+import Synapse.Tensors (DType, SingletonOps(singleton), MatOps(addMatRow, matMul))
+
+import Synapse.Tensors.Vec (Vec)
+
+import Synapse.Tensors.Mat (Mat)
+import qualified Synapse.Tensors.Mat as M
+
+import Synapse.Autograd (Symbolic, SymbolIdentifier(SymbolIdentifier), symbol, SymbolMat)
+
+import Synapse.NN.Layers.Layer (AbstractLayer(..), LayerConfiguration)
+import Synapse.NN.Layers.Initializers (Initializer(Initializer), zeroes, ones)
+import Synapse.NN.Layers.Constraints (Constraint(Constraint))
+import Synapse.NN.Layers.Regularizers (Regularizer(Regularizer))
+
+
+{- | 'Dense' datatype represents densely-connected neural network layer.
+
+'Dense' performs following operation: @x `matMul` w + b@, where @w@ is weights and @b@ is bias (if present) of a layer.
+-}
+data Dense a = Dense
+    { denseWeights      :: Mat a                           -- ^ Matrix that represents weights of dense layer.
+    , denseBias         :: Vec a                           -- ^ Vector that represents bias of dense layer.
+
+    , denseConstraints  :: (Constraint a, Constraint a)    -- ^ Constraints on weights and bias of dense layer.
+    , denseRegularizers :: (Regularizer a, Regularizer a)  -- ^ Regularizers on weights and bias of dense layer.
+    }
+
+-- | Creates symbol for weights.
+weightsSymbol :: SymbolIdentifier -> Mat a -> SymbolMat a
+weightsSymbol prefix = symbol (prefix <> SymbolIdentifier "1")
+
+
+-- | Creates symbol for bias.
+biasSymbol :: SymbolIdentifier -> Vec a -> SymbolMat a
+biasSymbol prefix = symbol (prefix <> SymbolIdentifier "2") . M.rowVec
+
+type instance DType (Dense a) = a
+
+instance AbstractLayer Dense where
+    inputSize = Just . M.nRows . denseWeights
+    outputSize = Just . M.nCols . denseWeights
+
+    nParameters _ = 2
+    getParameters prefix (Dense weights bias _ _) = [weightsSymbol prefix weights, biasSymbol prefix bias]
+    updateParameters (Dense _ _ constraints@(Constraint weightsConstraintFn, Constraint biasConstraintFn) regularizers) [weights', biasMat'] =
+        Dense (weightsConstraintFn weights') (M.indexRow (biasConstraintFn biasMat') 0) constraints regularizers
+    updateParameters _ _ = error "Parameters update failed - wrong amount of parameters was given"
+
+    symbolicForward prefix input (Dense weights bias _ (Regularizer weightsRegularizerFn, Regularizer biasRegularizerFn)) =
+        let symbolWeights = weightsSymbol prefix weights
+            symbolBias = biasSymbol prefix bias
+        in
+        ( (input `matMul` symbolWeights) `addMatRow` symbolBias
+        , weightsRegularizerFn symbolWeights + biasRegularizerFn symbolBias
+        )
+
+-- | Creates configuration of dense layer.
+layerDenseWith
+    :: Symbolic a
+    => (Initializer a, Constraint a, Regularizer a)  -- ^ Weights initializer, constraint and regularizer.
+    -> (Initializer a, Constraint a, Regularizer a)  -- ^ Bias initializer, constraint and regularizer.
+    -> Int                                           -- ^ Amount of neurons.
+    -> LayerConfiguration (Dense a)
+layerDenseWith (Initializer weightsInitializer, weightsConstraints, weightsRegularizer)
+               (Initializer biasInitializer, biasConstraints, biasRegularizer)
+               neurons input =
+    Dense (weightsInitializer (input, neurons)) (M.indexRow (biasInitializer (1, neurons)) 0) 
+          (weightsConstraints, biasConstraints) (weightsRegularizer, biasRegularizer)
+
+-- | Creates default configuration of dense layer - no constraints and weight are initialized with ones, bias is initialized with zeroes.
+layerDense :: Symbolic a => Int -> LayerConfiguration (Dense a)
+layerDense = layerDenseWith (Initializer ones, Constraint id, Regularizer (const $ singleton 0))
+                            (Initializer zeroes, Constraint id, Regularizer (const $ singleton 0))
+ src/Synapse/NN/Layers/Initializers.hs view
@@ -0,0 +1,181 @@+{- | Allows to initialize values of layers parameters.
+
+'InitializerFn' type alias represents functions that are able to initialize matrix with given size
+and 'Initializer' newtype wraps 'InitializerFn's.
+
+"Synapse" provides 4 types of initializers:
+* Non-random constant initializers
+* Random uniform distribution initializers
+* Random normal distribution initializers
+* Matrix-specific initializers
+-}
+
+
+module Synapse.NN.Layers.Initializers
+    ( -- * 'InitializerFn' type alias and 'Initializer' newtype
+
+      InitializerFn
+
+    , Initializer (Initializer, unInitializer)
+        
+      -- * Non-random constant initializers
+
+    , constants
+    , zeroes
+    , ones
+
+      -- * Random uniform distribution initializers
+
+    , randomUniform
+    , lecunUniform
+    , heUniform
+    , glorotUniform
+
+      -- * Random normal distribution initializers
+    
+    , randomNormal
+    , lecunNormal
+    , heNormal
+    , glorotNormal
+
+      -- * Matrix-like initializers
+    
+    , identity
+    , orthogonal
+    ) where
+
+
+import Synapse.Tensors.Mat (Mat)
+import qualified Synapse.Tensors.Mat as M
+
+import System.Random (uniformListR, uniformRs, UniformRange, RandomGen)
+
+
+-- | 'InitializerFn' type alias represents functions that are able to initialize matrix with given size.
+type InitializerFn a = (Int, Int) -> Mat a
+
+
+-- | 'Initializer' newtype wraps 'InitializerFn's - functions that are able to initialize matrix with given size.
+newtype Initializer a = Initializer 
+    { unInitializer :: InitializerFn a  -- ^ Unwraps 'Initializer' newtype.
+    }
+
+
+-- Non-random constant initializers
+
+-- | Initializes list that is filled with given constant.
+constants :: Num a => a -> InitializerFn a
+constants c (input, output) = M.replicate (input, output) c
+
+-- | Initializes list that is filled with zeroes.
+zeroes :: Num a => InitializerFn a
+zeroes = constants 0
+
+-- | Initializes list that is filled with ones.
+ones :: Num a => InitializerFn a
+ones = constants 1
+
+
+-- Random uniform distribution initializers
+
+{- | Initializes list with samples from random uniform distribution in range.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+randomUniform :: (UniformRange a, RandomGen g) => (a, a) -> g -> InitializerFn a
+randomUniform range gen sizes@(input, output) = M.fromList sizes $ fst $ uniformListR (input * output) range gen
+
+{- | Initializes list with samples from random LeCun uniform distribution in range.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+lecunUniform :: (UniformRange a, Floating a, RandomGen g) => g -> InitializerFn a
+lecunUniform gen sizes@(input, _) = let limit = sqrt $ 3.0 / fromIntegral input
+                                    in randomUniform (-limit, limit) gen sizes
+
+{- | Initializes list with samples from random He uniform distribution in range.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+heUniform :: (UniformRange a, Floating a, RandomGen g) => g -> InitializerFn a
+heUniform gen sizes@(input, _) = let limit = sqrt $ 6.0 / fromIntegral input
+                                 in randomUniform (-limit, limit) gen sizes
+
+{- | Initializes list with samples from random Glorot uniform distribution in range.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+glorotUniform :: (UniformRange a, Floating a, RandomGen g) => g -> InitializerFn a
+glorotUniform gen sizes@(input, output) = let limit = sqrt $ 6.0 / fromIntegral (input + output)
+                                          in randomUniform (-limit, limit) gen sizes
+
+
+-- Random normal distribution initializers
+
+{- | Initializes list with samples from random normal distribution in range which could be truncated.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+randomNormal :: (UniformRange a, Floating a, Ord a, RandomGen g) => Maybe a -> a -> a -> g -> InitializerFn a
+randomNormal truncated mean stdDev gen sizes@(input, output) = let us = pairs $ uniformRs (0.0, 1.0) gen
+                                                                   ns = concatMap ((\(n1, n2) -> [n1, n2]) . transformBoxMuller) us
+                                                                   ns' = map ((+ mean) . (* stdDev)) ns
+                                                                   ns'' = case truncated of
+                                                                               Nothing  -> ns'
+                                                                               Just eps -> filter (\x -> abs (x - mean) < eps) ns'
+                                                               in M.fromList sizes $ take (input * output) ns''
+  where
+    pairs [] = []
+    pairs [x] = [(x, 1)]
+    pairs (a:b:xs) = (a, b) : pairs xs
+
+    transformBoxMuller (u1, u2) = let r = sqrt $ (-2.0) * log u1
+                                      theta = 2.0 * pi * u2
+                                  in (r * cos theta, r * sin theta)
+
+{- | Initializes list with samples from random LeCun normal distribution in range
+which is truncated for values more than two standard deviations from mean.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+lecunNormal :: (UniformRange a, Floating a, Ord a, RandomGen g) => g -> InitializerFn a
+lecunNormal gen sizes@(input, _) = let mean = 0
+                                       stdDev = sqrt $ 1.0 / fromIntegral input
+                                   in randomNormal (Just $ 2.0 * stdDev) mean stdDev gen sizes
+
+{- | Initializes list with samples from random He normal distribution in range
+which is truncated for values more than two standard deviations from mean.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+heNormal :: (UniformRange a, Floating a, Ord a, RandomGen g) => g -> InitializerFn a
+heNormal gen sizes@(input, _) = let mean = 0
+                                    stdDev = sqrt $ 2.0 / fromIntegral input
+                                in randomNormal (Just $ 2.0 * stdDev) mean stdDev gen sizes
+
+{- | Initializes list with samples from random Glorot normal distribution in range
+which is truncated for values more than two standard deviations from mean.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+glorotNormal :: (UniformRange a, Floating a, Ord a, RandomGen g) => g -> InitializerFn a
+glorotNormal gen sizes@(input, output) = let mean = 0
+                                             stdDev = sqrt $ 2.0 / fromIntegral (input + output)
+                                         in randomNormal (Just $ 2.0 * stdDev) mean stdDev gen sizes
+
+
+-- Matrix-like initializers
+
+-- | Initializes flat identity matrix. If dimensions do not represent square matrix, an error is thrown.
+identity :: Num a => InitializerFn a
+identity (input, output)
+    | input /= output = error "Given dimensions do not represent square matrix"
+    | otherwise       = M.identity input
+
+{- | Initializes float orthogonal matrix obtained from a random normal distribution
+that is truncated for values more than two standard deviations from mean.
+
+This function does not preserve seed generator - split generator before calling this function.
+-}
+orthogonal :: (UniformRange a, Floating a, Ord a, RandomGen g) => g -> InitializerFn a
+orthogonal gen sizes = M.orthogonalized $ randomNormal Nothing 0.0 1.0 gen sizes
+ src/Synapse/NN/Layers/Layer.hs view
@@ -0,0 +1,101 @@+{- | This module provides nessesary abstraction over layers of neural networks.
+
+'AbstractLayer' typeclass defines interface of all layers of neural network model.
+Its implementation is probably the most low-leveled abstraction of the "Synapse" library.
+Notes on how to correctly implement that typeclass are in the docs for it.
+
+'Layer' is the existential datatype that wraps any 'AbstractLayer' instance.
+That is the building block of any neural network.
+-}
+
+-- 'TypeFamilies' are needed to instantiate 'DType'.
+{-# LANGUAGE TypeFamilies #-}
+
+-- 'ExistentialQuantification' is needed to define 'Layer' datatype.
+{-# LANGUAGE ExistentialQuantification #-}
+
+
+module Synapse.NN.Layers.Layer
+    ( -- * 'AbstractLayer' typeclass
+
+      AbstractLayer (inputSize, outputSize, nParameters, getParameters, updateParameters, symbolicForward)
+    , forward
+
+      -- * 'Layer' existential datatype
+
+    , Layer (Layer)
+
+      -- * 'LayerConfiguration' type alias
+    , LayerConfiguration
+    ) where
+
+
+import Synapse.Tensors (DType)
+import Synapse.Tensors.Mat (Mat)
+
+import Synapse.Autograd (Symbolic, SymbolIdentifier, Symbol(unSymbol), SymbolMat, constSymbol)
+
+
+{- | 'AbstractLayer' typeclass defines basic interface of all layers of neural network model.
+
+Every layer should be able to pass 'Synapse.Autograd.SymbolMat' through itself to produce new 'Synapse.Autograd.SymbolMat'
+(make prediction based on its parameters) using 'symbolicForward' function,
+which allows for gradients to be calculated after predictions, which in turn makes training possible.
+
+'nParameters', 'getParameters' and 'updateParameters' functions allow training of parameters of the layer.
+Their implementations should match - that is 'getParameters' function should return list of length 'nParameters'
+and 'updateParameters' should expect a list of the same length with the matrices in the same order as were they in 'getParameters'.
+
+"Synapse" manages gradients and parameters for layers with erased type information using prefix system.
+'Synapse.Autograd.SymbolIdentifier' is a prefix for name of symbolic parameters that are used in calculation.
+Every used parameter should have unique name to be recognised by the autograd -
+it must start with given prefix and end with the numerical index of said parameter.
+For example 3rd layer with 2 parameters (weights and bias) should
+name its weights symbol \"ml3w1\" and name its bias symbol \"ml3w2\" (\"ml3w\" prefix will be supplied).
+
+Important: this typeclass correct implementation is very important (as it is the \'heart\' of "Synapse" library)
+for work of the neural network and training, read the docs thoroughly to ensure that all the invariants are met.
+-}
+class AbstractLayer l where
+    -- | Returns the size of the input that is supported for 'forward' and 'symbolicForward' functions. 'Nothing' means size independence (activation functions are the example).
+    inputSize :: l a -> Maybe Int
+    -- | Returns the size of the output that is supported for 'forward' and 'symbolicForward' functions. 'Nothing' means size independence (activation functions are the example).
+    outputSize :: l a -> Maybe Int
+    -- | Returns the number of parameters of this layer.
+    nParameters :: l a -> Int
+    
+    -- | Returns a list of all parameters (those must be of the exact same order as they are named by the layer (check 'symbolicForward' docs)).
+    getParameters :: SymbolIdentifier -> l a -> [SymbolMat a]
+    -- | Updates parameters based on supplied list (length of that list, the order and the form of parameters is EXACTLY the same as those from 'getParameters')
+    updateParameters :: l a -> [Mat a] -> l a
+
+    {- | Passes symbolic matrix through to produce new symbolic matrix, while retaining gradients graph.
+    Second matrix is a result of application of regularizer on a layer.
+    -}
+    symbolicForward :: (Symbolic a, Floating a, Ord a) => SymbolIdentifier -> SymbolMat a -> l a -> (SymbolMat a, SymbolMat a)
+
+-- | Passes matrix through to produce new matrix.
+forward :: (AbstractLayer l, Symbolic a, Floating a, Ord a) => Mat a -> l a -> Mat a
+forward input = unSymbol . fst . symbolicForward mempty (constSymbol input)
+
+
+-- | 'Layer' existential datatype wraps anything that implements 'AbstractLayer'.
+data Layer a = forall l. (AbstractLayer l) => Layer (l a)
+
+type instance DType (Layer a) = a
+
+instance AbstractLayer Layer where
+    inputSize (Layer l) = inputSize l
+    outputSize (Layer l) = outputSize l
+
+    nParameters (Layer l) = nParameters l
+    getParameters prefix (Layer l) = getParameters prefix l
+    updateParameters (Layer l) = Layer . updateParameters l
+
+    symbolicForward prefix input (Layer l) = symbolicForward prefix input l
+
+
+-- | 'LayerConfiguration' type alias represents functions that are able to build layers.
+type LayerConfiguration l
+    =  Int  -- ^ Output size of previous layer.
+    -> l    -- ^ Resulting layer.
+ src/Synapse/NN/Layers/Regularizers.hs view
@@ -0,0 +1,44 @@+{- | Provides collection of functions that impose penalties on parameters which is done by adding result to loss value.
+-}
+
+
+module Synapse.NN.Layers.Regularizers
+    ( -- * 'RegularizerFn' type alias and 'Regularizer' newtype
+
+      RegularizerFn
+    , Regularizer (Regularizer, unRegularizer)
+
+      -- * Regularizers
+
+    , l1
+    , l2
+    ) where
+
+
+import Synapse.Tensors (ElementwiseScalarOps((*.)), SingletonOps(elementsSum))
+
+import Synapse.Autograd (SymbolMat, Symbolic)
+
+
+-- | 'RegularizerFn' type alias represents functions that impose penalties on parameters which is done by adding result of regularization to loss value.
+type RegularizerFn a = SymbolMat a -> SymbolMat a
+
+
+{- | 'Regularizer' newtype wraps 'RegularizerFn's - functions that impose penalties on parameters.
+
+Every regularization function must return symbol of singleton matrix.
+-}
+newtype Regularizer a = Regularizer
+    { unRegularizer :: RegularizerFn a  -- ^ Unwraps 'Regularizer' newtype.
+    }
+
+
+-- Regularizers
+
+-- | Applies a L1 regularization penalty (sum of absolute values of parameter multiplied by a coefficient).
+l1 :: (Symbolic a, Num a) => a -> RegularizerFn a
+l1 k mat = elementsSum (abs mat) *. k
+
+-- | Applies a L1 regularization penalty (sum of squared values of parameter multiplied by a coefficient).
+l2 :: (Symbolic a, Num a) => a -> RegularizerFn a
+l2 k mat = elementsSum (mat * mat) *. k
+ src/Synapse/NN/LearningRates.hs view
@@ -0,0 +1,69 @@+{- | Provides learning rate functions - functions that return coefficient which modulates how big are updates of parameters in training.
+-}
+
+
+module Synapse.NN.LearningRates
+    ( -- * 'LearningRateFn' type alias and 'LearningRate' newtype
+
+      LearningRateFn
+    , LearningRate (LearningRate, unLearningRate)
+
+      -- * Learning rate decay functions
+    
+    , exponentialDecay
+    , inverseTimeDecay
+    , polynomialDecay
+    , cosineDecay
+    , piecewiseConstantDecay
+    ) where
+
+
+-- | 'LearningRateFn' type alias represents functions that return coefficient which modulates how big are updates of parameters in training.
+type LearningRateFn a = Int -> a
+
+
+-- | 'LearningRate' newtype wraps 'LearningRateFn's - functions that modulate how big are updates of parameters in training.
+newtype LearningRate a = LearningRate
+    { unLearningRate :: LearningRateFn a  -- ^ Unwraps 'LearningRate' newtype.
+    }
+
+
+-- Learning rate decay functions
+
+-- | Takes initial learning rate, decay steps and decay rate and calculates exponential decay learning rate (@initial * decay_rate ^ (step / decay_steps)@).
+exponentialDecay :: Num a => a -> Int -> a -> LearningRateFn a
+exponentialDecay initial steps rate step = initial * rate ^ (step `div` steps)
+
+-- | Takes initial learning rate, decay steps and decay rate and calculates inverse time decay learning rate (@initial / (1 + rate * step / steps))@).
+inverseTimeDecay :: Fractional a => a -> Int -> a -> LearningRateFn a
+inverseTimeDecay initial steps rate step = initial / (1.0 + rate * fromIntegral (step `div` steps))
+
+{- | Takes initial learning rate, decay steps, polynomial power and end decay and calculates polynomial decay learning rate
+(@if step < steps then initial * (1 - step / steps) ** power else end@).
+-}
+polynomialDecay :: Floating a => a -> Int -> a -> a -> LearningRateFn a
+polynomialDecay initial steps power end step
+    | step < steps = initial * fromIntegral (1 - min step steps `div` steps) ** power
+    | otherwise    = end
+
+{- | Takes initial learning rate, decay steps, alpha coefficient and warmup steps and target (optional) and calculates cosine decay learning rate
+(@(1 - alpha) * (0.5 * (1.0 + cos (pi * step / steps))) + alpha) *
+(if warmup then (if step < warmupSteps then (warmupLR - initial) * step / warmupSteps else warmupLR) else initial)@).
+-}
+cosineDecay :: Floating a => a -> Int -> a -> Maybe (Int, a) -> LearningRateFn a
+cosineDecay initial steps alpha warmup step =
+    ((1.0 - alpha) * (0.5 * (1.0 + cos (pi * fromIntegral (step `div` steps)))) + alpha) *
+    case warmup of
+        Nothing -> initial
+        Just (warmupSteps, warmupLR) -> if step < warmupSteps
+                                        then (warmupLR - initial) * fromIntegral (step `div` warmupSteps)
+                                        else warmupLR
+
+{- | Takes list of boundaries and learning rate values and last rate value for those boundaries and calculates piecewise constant decay learning rate
+(@if step < bound1 then value1 else if step < bound2 then value2 else lastRate@ for @[(bound1, value1), (bound2, value2)]@).
+-}
+piecewiseConstantDecay :: [(Int, a)] -> a -> LearningRateFn a
+piecewiseConstantDecay [] lastRate _ = lastRate
+piecewiseConstantDecay ((stepBound, rateValue):xs) lastRate step
+    | step < stepBound = rateValue
+    | otherwise        = piecewiseConstantDecay xs lastRate step
+ src/Synapse/NN/Losses.hs view
@@ -0,0 +1,62 @@+{- | Provides collection of functions that are used to as a reference of what needs to be minimised during training.
+
+'LossFn' type alias represents those functions, and "Synapse" offers a variety of them.
+-}
+
+
+module Synapse.NN.Losses
+    ( -- * 'LossFn' type alias and 'Loss' newtype
+
+      LossFn
+    
+    , Loss (Loss, unLoss)
+
+      -- * Regression losses
+    
+    , mse
+    , msle
+    , mae
+    , mape
+    , logcosh
+    ) where
+
+
+import Synapse.Tensors (ElementwiseScalarOps((+.), (*.), (**.)), SingletonOps(mean))
+
+import Synapse.Autograd (SymbolMat, Symbolic)
+
+
+-- | 'LossFn' type alias represents functions that are able to provide a reference of what relation between matrices needs to be minimised.
+type LossFn a = SymbolMat a -> SymbolMat a -> SymbolMat a
+
+
+{- | 'Loss' newtype wraps 'LossFn's - differentiable functions that are able to provide a reference of what relation between matrices needs to be minimised.
+
+Every loss function must return symbol of singleton matrix.
+-}
+newtype Loss a = Loss 
+    { unLoss :: LossFn a  -- ^ Unwraps 'Loss' newtype.
+    }
+
+
+-- Regression losses
+
+-- | Computes the mean of squares of errors.
+mse :: (Symbolic a, Floating a) => LossFn a
+mse true predicted = mean $ (true - predicted) **. 2.0
+
+-- | Computes the mean squared logarithmic error.
+msle :: (Symbolic a, Floating a) => LossFn a
+msle true predicted = mean $ (log (true +. 1) - log (predicted +. 1)) **. 2.0
+
+-- | Computes the mean of absolute error.
+mae :: (Symbolic a, Floating a) => LossFn a
+mae true predicted = mean $ abs (true - predicted)
+
+-- | Computes the mean absolute percentage error.
+mape :: (Symbolic a, Floating a) => LossFn a
+mape true predicted = mean (abs (true - predicted) / true) *. 100
+
+-- | Computes the logarithm of the hyperbolic cosine of the error.
+logcosh :: (Symbolic a, Floating a) => LossFn a
+logcosh true predicted = mean $ log $ cosh (true - predicted)
+ src/Synapse/NN/Metrics.hs view
@@ -0,0 +1,38 @@+{- | Provides collection of functions that are used to judge the performance of models.
+
+'MetricFn' type alias represents those functions, and "Synapse" offers a variety of them.
+-}
+
+
+module Synapse.NN.Metrics
+    ( -- 'MetricFn' type alias and 'Metric' newtype
+
+      MetricFn
+    , lossFnToMetricFn
+    
+    , Metric (Metric, unMetric)
+    ) where
+
+
+import Synapse.NN.Losses (LossFn)
+
+import Synapse.Tensors.Mat (Mat)
+
+import Synapse.Autograd (Symbol(unSymbol), constSymbol)
+
+
+-- | 'MetricFn' type alias represents functions that are able to provide a reference of performance of neural network model.
+type MetricFn a = Mat a -> Mat a -> Mat a
+
+-- | Converts any loss function to a metric function (because the same constraint is imposed on both 'MetricFn' and 'LossFn').
+lossFnToMetricFn :: LossFn a -> MetricFn a
+lossFnToMetricFn loss true predicted = unSymbol $ loss (constSymbol true) (constSymbol predicted)
+
+
+{- | 'Metric' newtype wraps 'MetricFn's - functions that are able to provide a reference of performance of neural network model.
+
+Every metric function must return symbol of singleton matrix.
+-}
+newtype Metric a = Metric 
+    { unMetric :: MetricFn a  -- ^ Unwraps 'Metric' newtype.
+    }
+ src/Synapse/NN/Models.hs view
@@ -0,0 +1,80 @@+{- | Provides interface for creating and using neural network models.
+-}
+
+
+-- 'TypeFamilies' are needed to instantiate 'Synapse.Tensors.DType'.
+{-# LANGUAGE TypeFamilies #-}
+
+
+module Synapse.NN.Models
+    ( -- * Common for models
+      InputSize (InputSize)
+
+      -- * 'SequentialModel' datatype
+      
+    , SequentialModel (SequentialModel, unSequentialModel)
+    , buildSequentialModel
+
+    , layerPrefix
+    ) where
+
+
+import Synapse.Tensors (DType, SingletonOps(singleton))
+
+import Synapse.Autograd (SymbolIdentifier (SymbolIdentifier))
+
+import Synapse.NN.Layers.Layer(AbstractLayer(..), Layer, LayerConfiguration)
+
+import Data.Maybe (fromMaybe)
+import Data.Foldable (foldl')
+
+
+-- | 'InputSize' newtype wraps 'Int' - amount of features of input that the model should support (@InputSize 3@ means that model supports any matrix with size (x, 3)).
+newtype InputSize = InputSize Int
+
+
+-- | 'SequentialModel' datatype represents any model grouping layers linearly.
+newtype SequentialModel a = SequentialModel
+    { unSequentialModel :: [Layer a]  -- ^ Returns layers of 'SequentialModel'.
+    }
+
+instance Show a => Show (SequentialModel a) where
+    show (SequentialModel layers) = go 1 layers
+      where
+        go _ [] = ""
+        go i (x:xs) = "Layer " ++ show i ++ " parameters: " ++ show (getParameters (layerPrefix mempty i) x) ++ ";\n" ++ go (i + 1) xs
+
+-- | Builds sequential model using input size and layer configurations to ensure that layers are compatible with each other.
+buildSequentialModel :: InputSize -> [LayerConfiguration (Layer a)] -> SequentialModel a
+buildSequentialModel (InputSize i) layerConfigs = SequentialModel $ go i layerConfigs
+  where
+    go _ [] = []
+    go prevSize (l:ls) = let layer = l prevSize
+                             outputMaybe = outputSize layer
+                             output = fromMaybe prevSize outputMaybe
+                         in layer : go output ls
+
+type instance DType (SequentialModel a) = a
+
+
+-- | Forms prefix for layers according to 'Synapse.NN.Layers.Layer.AbstractLayer' requirements.
+layerPrefix :: SymbolIdentifier -> Int -> SymbolIdentifier
+layerPrefix prefix i = mconcat [prefix, SymbolIdentifier "l", SymbolIdentifier (show i), SymbolIdentifier "w"]
+
+instance AbstractLayer SequentialModel where
+    inputSize = inputSize . head . unSequentialModel
+    outputSize = outputSize . head . unSequentialModel
+
+    nParameters = foldl' (\parameters layer -> parameters + nParameters layer) 0 . unSequentialModel
+    getParameters prefix =
+        snd . foldl' (\(i, acc) layer -> (i + 1, acc ++ getParameters (layerPrefix prefix i) layer)) (1, []) . unSequentialModel
+    updateParameters (SequentialModel model) = SequentialModel . go model
+      where
+        go [] _ = []
+        go (layer:layers) parameters = let (x, parameters') = splitAt (nParameters layer) parameters
+                                       in updateParameters layer x : go layers parameters'
+
+    symbolicForward prefix input =
+        snd . foldl' (\(i, (mat, loss)) layer -> let (mat', newLoss) = symbolicForward (layerPrefix prefix i) mat layer
+                                                 in (i + 1, (mat', loss + newLoss)))
+              (1, (input, singleton 0)) . unSequentialModel
+ src/Synapse/NN/Optimizers.hs view
@@ -0,0 +1,82 @@+{- | This module implements several optimizers that are used in training.
+-}
+
+
+-- 'TypeFamilies' are needed to use 'DType' and define 'Optimizer' typeclass.
+{-# LANGUAGE TypeFamilies #-}
+
+
+module Synapse.NN.Optimizers
+    ( -- * 'Optimizer' typeclass
+      
+      Optimizer (OptimizerParameters, optimizerInitialParameters, optimizerUpdateStep)
+
+    , optimizerUpdateParameters
+    
+      -- * Optimizers
+
+    , SGD (SGD, sgdMomentum, sgdNesterov)
+    ) where
+
+
+import Synapse.Tensors (DType, ElementwiseScalarOps((*.)))
+
+import Synapse.Tensors.Mat (Mat)
+import qualified Synapse.Tensors.Mat as M
+
+import Synapse.Autograd (Symbolic, Symbol(unSymbol), SymbolMat, Gradients, wrt)
+
+import Synapse.NN.Layers.Initializers (zeroes)
+
+import Data.Kind (Type)
+
+
+-- | 'Optimizer' typeclass represents optimizer - algorithm that defines an update rule of neural network parameters.
+class Optimizer optimizer where
+    -- | 'OptimizerParameters' represent optimizer-specific parameters that it needs to implement update rule.
+    type OptimizerParameters optimizer a :: Type
+
+    -- | Returns initial state of optimizer-specific parameters for given variable.
+    optimizerInitialParameters :: Num a => optimizer a -> Mat a -> OptimizerParameters optimizer a
+    
+    -- | Performs the update step of optimizer.
+    optimizerUpdateStep
+        :: Num a
+        => optimizer a                               -- ^ Optimizer itself.
+        -> (a, Mat a)                                -- ^ Learning rate and gradient of given parameter.
+        -> (Mat a, OptimizerParameters optimizer a)  -- ^ Given parameter and current state of optimizer-specific parameters.
+        -> (Mat a, OptimizerParameters optimizer a)  -- ^ Updated parameter and a new state of optimizer-specific parameters.
+
+-- | 'optimizerUpdateParameters' function updates whole model using optimizer by performing 'optimizerUpdateStep' for every parameter.
+optimizerUpdateParameters
+    :: (Symbolic a, Optimizer optimizer)
+    => optimizer a                                       -- ^ Optimizer itself.
+    -> (a, Gradients (Mat a))                            -- ^ Learning rate and gradients of all parameters. 
+    -> [(SymbolMat a, OptimizerParameters optimizer a)]  -- ^ Given parameters and current state of optimizer-specific parameters.
+    -> [(Mat a, OptimizerParameters optimizer a)]        -- ^ Updated parameters and a new state of optimizer-specific parameters.
+optimizerUpdateParameters _ _ [] = []
+optimizerUpdateParameters optimizer (lrValue, gradients) ((parameter, optimizerParameter):xs) =
+    optimizerUpdateStep optimizer (lrValue, unSymbol $ gradients `wrt` parameter) (unSymbol parameter, optimizerParameter)
+    : optimizerUpdateParameters optimizer (lrValue, gradients) xs
+
+
+-- | 'SGD' is a optimizer that implements stochastic gradient-descent algorithm.
+data SGD a = SGD
+    { sgdMomentum :: a     -- ^ Momentum coefficient.
+    , sgdNesterov :: Bool  -- ^ Nesterov update rule.
+    } deriving (Eq, Show)
+
+type instance DType (SGD a) = a
+
+instance Optimizer SGD where
+    type OptimizerParameters SGD a = Mat a
+
+    optimizerInitialParameters _ parameter = zeroes (M.size parameter)
+
+    optimizerUpdateStep (SGD momentum nesterov) (lr, gradient) (parameter, velocity) = (parameter', velocity')
+      where
+        velocity' = velocity *. momentum - gradient *. lr
+
+        parameter' = if nesterov
+                     then parameter + velocity' *. momentum - gradient *. lr
+                     else parameter + velocity'
+ src/Synapse/NN/Training.hs view
@@ -0,0 +1,265 @@+{- | This module provides datatypes and functions that implement neural networks training.
+-}
+
+
+-- | 'OverloadedStrings' are needed to use strings in progress bar.
+{-# LANGUAGE OverloadedStrings #-}
+
+
+module Synapse.NN.Training
+    ( -- * 'Callbacks' datatype and associated type aliases
+
+      CallbackFnOnTrainBegin
+    , CallbackFnOnEpochBegin
+    , CallbackFnOnBatchBegin
+    , CallbackFnOnBatchEnd
+    , CallbackFnOnEpochEnd
+    , CallbackFnOnTrainEnd
+
+    , Callbacks
+          ( Callbacks
+          , callbacksOnTrainBegin
+          , callbacksOnEpochBegin
+          , callbacksOnBatchBegin
+          , callbacksOnBatchEnd
+          , callbacksOnEpochEnd
+          , callbacksOnTrainEnd
+          )
+    , emptyCallbacks
+
+    -- * 'Hyperparameters' datatype
+
+    , Hyperparameters
+          ( Hyperparameters
+          , hyperparametersEpochs
+          , hyperparametersBatchSize
+          , hyperparametersDataset
+          , hyperparametersLearningRate
+          , hyperparametersLoss
+          , hyperparametersMetrics
+          )
+
+    , RecordedMetric (RecordedMetric, unRecordedMetric)
+
+      -- * Training
+
+    , train
+    ) where
+
+
+import Synapse.Tensors (SingletonOps(unSingleton))
+import Synapse.Tensors.Vec (Vec(Vec), unVec)
+import Synapse.Tensors.Mat (Mat)
+
+import Synapse.Autograd (Symbolic, SymbolIdentifier(SymbolIdentifier), unSymbol, symbol, constSymbol, getGradientsOf)
+
+import Synapse.NN.Layers.Layer (AbstractLayer(..))
+import Synapse.NN.Optimizers (Optimizer(..), optimizerUpdateParameters)
+import Synapse.NN.Batching (Sample(Sample), unDataset, shuffleDataset, VecDataset, batchVectors, BatchedDataset)
+import Synapse.NN.LearningRates (LearningRate(LearningRate))
+import Synapse.NN.Losses (Loss(Loss))
+import Synapse.NN.Metrics (Metric(Metric))
+
+import Data.Functor ((<&>))
+import Control.Monad (forM_, when)
+import Data.IORef (IORef, newIORef, writeIORef, readIORef, modifyIORef')
+
+import System.Random (RandomGen)
+
+import System.ProgressBar
+
+import qualified Data.Vector as V
+import qualified Data.Vector.Mutable as MV
+
+
+-- | Type of callback that is called at the beginning of training.
+type CallbackFnOnTrainBegin model optimizer a
+    =  IORef (model a)                          -- ^ Initial model state.
+    -> IORef [OptimizerParameters optimizer a]  -- ^ Initial optimizer parameters.
+    -> IO ()
+
+-- | Type of callback that is called at the beginning of training epoch.
+type CallbackFnOnEpochBegin model optimizer a
+    =  IORef Int                                -- ^ Current epoch.
+    -> IORef (model a)                          -- ^ Model state at the beginning of the epoch processing.
+    -> IORef [OptimizerParameters optimizer a]  -- ^ Optimizer parameters at the beginning of the epoch processing.
+    -> IORef (BatchedDataset a)                 -- ^ Batched shuffled dataset.
+    -> IO ()
+
+-- | Type of callback that is called at the beginning of training batch.
+type CallbackFnOnBatchBegin model optimizer a
+    =  IORef Int                                -- ^ Current epoch.
+    -> IORef Int                                -- ^ Current batch number.
+    -> IORef (model a)                          -- ^ Model state at the beginning of the batch processing.
+    -> IORef [OptimizerParameters optimizer a]  -- ^ Optimizer parameters at the beginning of the batch processing.
+    -> IORef (Sample (Mat a))                   -- ^ Batch that is being processed.
+    -> IORef a                                  -- ^ Learning rate value.
+    -> IO ()
+
+-- | Type of callback that is called at the end of training batch.
+type CallbackFnOnBatchEnd model optimizer a
+    =  IORef Int                                -- ^ Current epoch.
+    -> IORef Int                                -- ^ Current batch number.
+    -> IORef (model a)                          -- ^ Model state at the end of the batch processing.
+    -> IORef [OptimizerParameters optimizer a]  -- ^ Optimizer parameters at the end of the batch processing.
+    -> IORef (Vec a)                            -- ^ Metrics that were recorded on this batch.
+    -> IO ()
+
+-- | Type of callback that is called at the end of training epoch.
+type CallbackFnOnEpochEnd model optimizer a
+    =  IORef Int                                -- ^ Current epoch.
+    -> IORef (model a)                          -- ^ Model state at the end of the epoch processing.
+    -> IORef [OptimizerParameters optimizer a]  -- ^ Optimizer parameters at the end of the epoch processing.
+    -> IO ()
+
+-- | Type of callback that is called at the end of training.
+type CallbackFnOnTrainEnd model optimizer a
+    =  IORef (model a)                            -- ^ Model state at the end of the training.
+    -> IORef [OptimizerParameters optimizer a]    -- ^ Optimizer parameters at the end of the training.
+    -> IORef (Vec (RecordedMetric a))             -- ^ Recorded metrics.
+    -> IO ()
+
+{- | 'Callbacks' record datatype holds all callbacks for the training.
+
+All callbacks take 'IORef's to various training parameters,
+which allows to affect training in any way possible.
+
+This interface should be used with caution, because some changes might break the training completely.
+-}
+data Callbacks model optimizer a = Callbacks
+    { callbacksOnTrainBegin :: [CallbackFnOnTrainBegin model optimizer a]  -- ^ Callbacks that will be called at the beginning of training.
+    , callbacksOnEpochBegin :: [CallbackFnOnEpochBegin model optimizer a]  -- ^ Callbacks that will be called at the beginning of training epoch processing.
+    , callbacksOnBatchBegin :: [CallbackFnOnBatchBegin model optimizer a]  -- ^ Callbacks that will be called at the beginning of training batch processing.
+    , callbacksOnBatchEnd   :: [CallbackFnOnBatchEnd   model optimizer a]  -- ^ Callbacks that will be called at the end of training batch processing.
+    , callbacksOnEpochEnd   :: [CallbackFnOnEpochEnd   model optimizer a]  -- ^ Callbacks that will be called at the end of training epoch processing.
+    , callbacksOnTrainEnd   :: [CallbackFnOnTrainEnd   model optimizer a]  -- ^ Callbacks that will be called at the end of training.
+    }
+
+-- | Returns empty 'Callbacks' record. It could also be used to build your own callbacks upon.
+emptyCallbacks :: Callbacks model optimizer a
+emptyCallbacks = Callbacks [] [] [] [] [] []
+
+
+-- | 'Hyperparameters' datatype represents configuration of a training.
+data Hyperparameters a = Hyperparameters
+    { hyperparametersEpochs       :: Int             -- ^ Number of epochs in the training.
+    , hyperparametersBatchSize    :: Int             -- ^ Size of batches that will be used in the training.
+
+    , hyperparametersDataset      :: VecDataset a    -- ^ Dataset with samples of vector functions.
+
+    , hyperparametersLearningRate :: LearningRate a  -- ^ 'LearningRate' that will be used in the training.
+    , hyperparametersLoss         :: Loss a          -- ^ 'Loss' that will be used in the training.
+
+    , hyperparametersMetrics      :: Vec (Metric a)  -- ^ 'Metric's that will be recorded during training.
+    }
+
+-- | 'RecordedMetric' newtype wraps vector of results of metrics.
+newtype RecordedMetric a = RecordedMetric
+    { unRecordedMetric :: Vec a  -- ^ Results of metric recording.
+    }
+
+
+-- | 'whileM_' function implements a monadic @while@ loop which can be @break@ed if the condition becomes false.
+whileM_ :: (Monad m) => m Bool -> m a -> m ()
+whileM_ p f = go
+  where
+    go = p >>= flip when (f >> go)
+
+-- | 'train' function allows training neural networks on datasets with specified parameters.
+train
+    :: (Symbolic a, Floating a, Ord a, Show a, RandomGen g, AbstractLayer model, Optimizer optimizer)
+    => model a                                                                     -- ^ Trained model.
+    -> optimizer a                                                                 -- ^ Optimizer that will be during training.
+    -> Hyperparameters a                                                           -- ^ Hyperparameters of training.
+    -> Callbacks model optimizer a                                                 -- ^ Callbacks that will be used during training.
+    -> g                                                                           -- ^ Generator of random values that will be used to shuffle dataset.
+    -> IO (model a, [OptimizerParameters optimizer a], Vec (RecordedMetric a), g)  -- ^ Updated model, optimizer parameters at the end of training, vector of recorded metrics (loss is also recorded and is the first in vector), updated generator of random values.
+train model optimizer (Hyperparameters epochs batchSize dataset (LearningRate lr) (Loss loss) (Vec metrics)) callbacks gen0 = 
+    let modelIdentifier = SymbolIdentifier "m"
+        inputIdentifier = SymbolIdentifier "input"
+    in do
+    let totalIterations = epochs * ((V.length (unVec $ unDataset dataset) + batchSize - 1) `div` batchSize)
+    progressBar <- newProgressBar
+        (defStyle { stylePrefix = exact <> msg " iterations"
+                  , stylePostfix = msg "\nElapsed time: " <> elapsedTime renderDuration
+                                <> msg "; Remaining time: " <> remainingTime renderDuration ""
+                                <> msg "; Total time: " <> totalTime renderDuration ""
+                  }
+        ) 10 (Progress 0 totalIterations ())
+
+    modelRef <- newIORef model
+
+    optimizerParametersRef <- newIORef $ optimizerInitialParameters optimizer . unSymbol <$> getParameters modelIdentifier model
+
+    mapM_ (\fn -> fn modelRef optimizerParametersRef) (callbacksOnTrainBegin callbacks)
+
+    allMetrics <- V.generateM (1 + V.length metrics) (const $ MV.new totalIterations)
+
+    gen <- newIORef gen0
+
+
+    epochRef <- newIORef 1
+    whileM_ (readIORef epochRef <&> (<= epochs)) $ do
+
+        currentGen <- readIORef gen
+        let (shuffledDataset, gen') = shuffleDataset dataset currentGen
+        _ <- writeIORef gen gen'
+
+        batchedDatasetRef <- newIORef $ batchVectors batchSize shuffledDataset
+
+        mapM_ (\fn -> fn epochRef modelRef optimizerParametersRef batchedDatasetRef) (callbacksOnEpochBegin callbacks)
+
+        epoch <- readIORef epochRef
+
+        batchedDataset <- readIORef batchedDatasetRef
+        batchIRef <- newIORef 1
+        whileM_ (readIORef batchIRef <&> (<= V.length (unVec $ unDataset batchedDataset))) $ do
+            incProgress progressBar 1
+
+            batchI <- readIORef batchIRef
+            let currentIteration = epoch * batchI
+
+            batchRef <- newIORef $ V.unsafeIndex (unVec $ unDataset batchedDataset) (batchI - 1)
+            lrValueRef <- newIORef $ lr currentIteration
+
+            mapM_ (\fn -> fn epochRef batchIRef modelRef optimizerParametersRef batchRef lrValueRef) (callbacksOnBatchBegin callbacks)
+
+            (Sample batchInput batchOutput) <- readIORef batchRef
+
+            (prediction, regularizersLoss) <- readIORef modelRef <&> symbolicForward modelIdentifier (symbol inputIdentifier batchInput)
+
+            lrValue <- readIORef lrValueRef
+            let lossValue = loss (constSymbol batchOutput) prediction
+
+            parameters <- readIORef modelRef <&> getParameters modelIdentifier
+            optimizerParameters <- readIORef optimizerParametersRef
+            let (parameters', optimizerParameters') = unzip $ optimizerUpdateParameters optimizer (lrValue, getGradientsOf $ lossValue + regularizersLoss)
+                                                                                        (zip parameters optimizerParameters)
+
+            _ <- modifyIORef' modelRef (`updateParameters` parameters')
+            _ <- writeIORef optimizerParametersRef optimizerParameters'
+
+            metricsValuesRef <- newIORef $ Vec $ V.cons (unSingleton lossValue) $ V.map (\(Metric metric) -> unSingleton $ metric batchOutput (unSymbol prediction)) metrics
+
+            mapM_ (\fn -> fn epochRef batchIRef modelRef optimizerParametersRef metricsValuesRef) (callbacksOnBatchEnd callbacks)
+
+            (Vec metricsValues) <- readIORef metricsValuesRef
+            forM_ [0 .. (V.length metricsValues - 1)] $ \metricI -> MV.write (V.unsafeIndex allMetrics metricI) (currentIteration - 1) (V.unsafeIndex metricsValues metricI)
+
+            modifyIORef' batchIRef (+ 1)
+
+        mapM_ (\fn -> fn epochRef modelRef optimizerParametersRef) (callbacksOnEpochEnd callbacks)
+
+        modifyIORef' epochRef (+ 1)
+
+
+    recordedMetricsRef <- V.mapM (fmap (RecordedMetric . Vec) . V.unsafeFreeze) allMetrics >>= newIORef . Vec
+
+    mapM_ (\fn -> fn modelRef optimizerParametersRef recordedMetricsRef) (callbacksOnTrainEnd callbacks)
+
+    trainedModel <- readIORef modelRef
+    trainedOptimizerParameters <- readIORef optimizerParametersRef
+    recordedMetrics <- readIORef recordedMetricsRef
+    gen'' <- readIORef gen
+
+    return (trainedModel, trainedOptimizerParameters, recordedMetrics, gen'')
+ src/Synapse/Tensors.hs view
@@ -0,0 +1,140 @@+{- | Module that provides mathematical base for neural networks.
+
+This module implements 'Synapse.Tensors.Vec.Vec' and 'Synapse.Tensors.Mat.Mat' datatypes and provides
+several useful function to work with them.
+
+Most of typeclasses of this module are working with 'DType' type family.
+That is to permit instances on types that are not exactly containers, but rather wrappers of containers,
+and it allows imposing additional constraints on inner type.
+The best example is 'Synapse.Autograd.Symbol' from "Synapse.Autograd".
+-}
+
+
+{- 'ConstrainedClassMethods', 'FlexibleContexts', 'TypeFamilies', 
+are needed to define 'DType' typefamily, 'Indexable', 'ElementwiseScalarOps', 'SingletonOps', 'VecOps', 'MatOps' typeclasses.
+-}
+{-# LANGUAGE ConstrainedClassMethods #-}
+{-# LANGUAGE FlexibleContexts        #-}
+{-# LANGUAGE TypeFamilies            #-}  
+
+
+module Synapse.Tensors
+    ( -- * 'DType' type family
+      
+      DType
+        
+      -- * 'Indexable' typeclass
+
+    , Indexable (Index, unsafeIndex, (!), (!?))
+
+      -- * Container-scalar operations
+    , ElementwiseScalarOps ((+.), (-.), (*.), (/.), (**.), elementsMin, elementsMax)
+    , SingletonOps (singleton, isSingleton, unSingleton, extendSingleton, elementsSum, elementsProduct, mean, norm)
+
+      -- * Specific container operations
+    , VecOps (dot)
+    , MatOps (transpose, addMatRow, matMul)
+    ) where
+
+
+import Data.Kind (Type)
+
+
+-- | 'DType' type family allows to get type of element for any container of that type - even for nested ones!
+type family DType f :: Type
+
+
+infixl 9 !, !?
+-- | 'Indexable' typeclass provides indexing interface for datatypes. 
+class Indexable f where
+    -- | Type of index for 'Indexable' container.
+    type Index f :: Type
+
+    -- | Unsafe indexing.
+    unsafeIndex :: f -> Index f -> DType f
+
+    -- | Indexing with bounds checking.
+    (!) :: f -> Index f -> DType f
+
+    -- | Safe indexing.
+    (!?) :: f -> Index f -> Maybe (DType f)
+
+
+infixl 6 +., -.
+infixl 7 *., /.
+infixr 8 **.
+{- | 'ElementwiseScalarOps' typeclass allows containers over numerical values easily work with scalars by using elementwise operations.
+
+This typeclass operates on 'DType' to permit instances on types that are not exactly containers, but rather wrappers of containers.
+The best example is 'Synapse.Autograd.Symbol' from "Synapse.Autograd".
+-}
+class ElementwiseScalarOps f where
+    -- | Adds given value to every element of the container.
+    (+.) :: Num (DType f) => f -> DType f -> f
+    -- | Subtracts given value from every element of the functor.
+    (-.) :: Num (DType f) => f -> DType f -> f
+    -- | Multiplies every element of the functor by given value.
+    (*.) :: Num (DType f) => f -> DType f -> f
+    -- | Divides every element of the functor by given value.
+    (/.) :: Fractional (DType f) => f -> DType f -> f
+    -- | Exponentiates every element of the functor by given value.
+    (**.) :: Floating (DType f) => f -> DType f -> f
+
+    -- | Applies 'min' operation with given value to every element.
+    elementsMin :: Ord (DType f) => f -> DType f -> f
+    -- | Applies 'max' operation with given value to every element.
+    elementsMax :: Ord (DType f) => f -> DType f -> f
+
+{- | 'SingletonOps' typeclass provides operations that relate to singleton containers (scalars that are wrapped in said container).
+
+All functions of that typeclass must return singletons (scalars that are wrapped in container).
+
+This typeclass operates on 'DType' to permit instances on types that are not exactly containers, but rather wrappers of containers.
+The best example is 'Synapse.Autograd.Symbol' from "Synapse.Autograd".
+-}
+class SingletonOps f where
+    -- | Initializes singleton container.
+    singleton :: DType f -> f
+    -- | Return true if container represents a singleton.
+    isSingleton :: f -> Bool
+    -- | Unwraps singleton container.
+    unSingleton :: f -> DType f
+
+    -- Extends singleton to needed size using second argument as a reference.
+    extendSingleton :: f -> f -> f
+
+    -- | Sums all elements of container.
+    elementsSum :: Num (DType f) => f -> f
+    -- | Multiplies all elements of container ('Fractional' constraint is needed for efficient gradient calculation, although it may be overly restrictive in some situations).
+    elementsProduct :: Fractional (DType f) => f -> f
+    
+    -- | Calculates the mean of all elements of container.
+    mean :: Fractional (DType f) => f -> f
+
+    -- | Calculates the Frobenius norm of all elements of container.
+    norm :: Floating (DType f) => f -> f
+
+
+{- | 'VecOps' typeclass provides vector-specific operations.
+
+This typeclass operates on 'DType' to permit instances on types that are not exactly containers, but rather wrappers of containers.
+The best example is 'Synapse.Autograd.Symbol' from "Synapse.Autograd".
+-}
+class VecOps f where
+    -- | Calculates dot product of two vectors.
+    dot :: Num (DType f) => f -> f -> f
+
+{- | 'MatOps' typeclass provides matrix-specific operations.
+
+This typeclass operates on 'DType' to permit instances on types that are not exactly containers, but rather wrappers of containers.
+The best example is 'Synapse.Autograd.Symbol' from "Synapse.Autograd".
+-}
+class MatOps f where
+    -- | Transposes matrix.
+    transpose :: f -> f
+
+    -- | Add matrix that represents row to every row of given matrix.
+    addMatRow :: Num (DType f) => f -> f -> f
+
+    -- | Mutiplies two matrices.
+    matMul :: Num (DType f) => f -> f -> f
+ src/Synapse/Tensors/Mat.hs view
@@ -0,0 +1,558 @@+{- | Implementation of matrix.
+
+Matrices are a backbone of any machine learning library,
+since most of the operations are implemented by the matrices
+combinations (matrix multiplication, elementwise operations).
+
+'Mat' datatype provides interface for all of those operations.
+-}
+
+
+-- 'TypeFamilies' are needed to instantiate 'Indexable', 'ElementwiseScalarOps', 'SingletonOps', 'MatOps' typeclasses.
+{-# LANGUAGE TypeFamilies #-}
+
+
+module Synapse.Tensors.Mat
+    ( --  * 'Mat' datatype and simple getters.
+
+      Mat (nRows, nCols)
+
+    , nElements
+    , size
+
+    , isTransposed
+    , isSubmatrix
+
+      -- * Utility
+
+    , force
+    , toLists
+
+      -- * Constructors
+
+    , empty
+    , singleton
+    , fromList
+    , fromLists
+    , generate
+    , replicate
+
+      -- * Vec operations
+
+    , rowVec
+    , colVec
+    , fromVec
+
+    , indexRow
+    , indexCol
+    , safeIndexRow
+    , safeIndexCol
+
+    , diagonal
+    , flatten
+
+     -- * Combining
+
+    , map
+    , mapRow
+    , mapCol
+    , for
+    , imap
+    , elementwise
+
+      -- * Operations with matrices
+
+    , setSize
+    , extend
+    , shrink
+
+    , swapRows
+    , swapCols
+    , transpose
+
+      -- * Submatrices
+
+    , minor
+    , submatrix
+    , split
+    , join
+    , (<|>)
+    , (<->)
+
+      -- * Mathematics
+
+    , zeroes
+    , ones
+    , identity
+
+    , adamarMul
+    , matMul
+
+    , det
+    , rref
+    , inverse
+    , orthogonalized
+    ) where
+
+
+import Synapse.Tensors (DType, Indexable(..), (!), ElementwiseScalarOps(..), SingletonOps(..), MatOps(..))
+
+import Synapse.Tensors.Vec (Vec(Vec))
+import qualified Synapse.Tensors.Vec as SV
+
+import Prelude hiding (map, replicate, zip)
+import Data.Foldable (Foldable(..))
+import Data.List (find)
+import Data.Tuple (swap)
+
+import qualified Data.Vector as V
+
+
+{- | Mathematical matrix (collection of elements).
+
+This implementation focuses on sharing parts of matrices and clever indexation
+to reduce overhead of several essential operations.
+Those include splitting matrices into submatrices, transposing - their asymptotical complexity becomes O(1).
+However there are few downsides:
+the first is that severely splitted matrix is hard to garbage collect and is not cache-friendly
+and the second is that mass traversal operations on those sparse matrices might not fuse and combine well.
+'force' function and some functions that by their nature are forced address those issues,
+but most of the time those problems are not significant enough and you are just better using convenient functions instead of workarounds.
+-}
+data Mat a = Mat
+    { nRows     :: {-# UNPACK #-} !Int         -- ^ Number of rows. 
+    , nCols     :: {-# UNPACK #-} !Int         -- ^ Number of columns.
+    , rowStride :: {-# UNPACK #-} !Int         -- ^ How much increasing row index affects true indexing.
+    , colStride :: {-# UNPACK #-} !Int         -- ^ How much increasing column index affects true indexing.
+    , rowOffset :: {-# UNPACK #-} !Int         -- ^ Row offset (from which row index does the matrix actually start).
+    , colOffset :: {-# UNPACK #-} !Int         -- ^ Column offset (from which column index does the matrix actually start).
+    , storage   ::                 V.Vector a  -- ^ Internal storage (elements are stored in a vector using row-major ordering).
+    }
+
+-- | Number of elements in a matrix.
+nElements :: Mat a -> Int
+nElements mat = nRows mat * nCols mat
+
+-- | Size of matrix.
+size :: Mat a -> (Int, Int)
+size mat = (nRows mat, nCols mat)
+
+
+-- | Returns whether the matrix is transposed. If the matrix consists of only one element, it is considered never transposed.
+isTransposed :: Mat a -> Bool
+isTransposed mat = colStride mat /= 1 && rowStride mat == 1
+
+-- | Returns whether the matrix is a submatrix from another matrix.
+isSubmatrix :: Mat a -> Bool
+isSubmatrix mat = any (0 /=) [rowOffset mat, colOffset mat]
+
+
+-- | Converts two dimensional matrix index to one dimensional vector index.
+indexMatToVec :: Mat a -> (Int, Int) -> Int
+indexMatToVec (Mat _ _ rk ck r0 c0 _) (r, c) = (r0 + r) * rk + (c0 + c) * ck
+
+-- | Converts one dimensional vector index to two dimensional matrix index.
+indexVecToMat :: Mat a -> Int -> (Int, Int)
+indexVecToMat mat@(Mat _ _ rk ck r0 c0 _) i = let t = isTransposed mat
+                                                  (r', c') = (if t then swap else id) $ quotRem i (if t then ck else rk)
+                                              in (r' - r0, c' - c0)
+
+
+
+-- | Copies matrix data dropping any extra memory that may be held if given matrix is a submatrix.
+force :: Mat a -> Mat a
+force x = Mat (nRows x) (nCols x) (nCols x) 1 0 0 (V.fromList [unsafeIndex x (r, c) | r <- [0 .. nRows x - 1], c <- [0 .. nCols x - 1]])
+
+-- | Converts matrix to list of lists.
+toLists :: Mat a -> [[a]]
+toLists x = [[unsafeIndex x (r, c) | c <- [0 .. nCols x - 1] ] | r <- [0 .. nRows x - 1]]
+
+
+-- Typeclasses
+
+instance Show a => Show (Mat a) where
+    show mat = "(" ++ show (size mat) ++ "): " ++ show (toLists mat)
+
+
+type instance DType (Mat a) = a
+
+
+instance Indexable (Mat a) where
+    type Index (Mat a) = (Int, Int)
+
+    unsafeIndex x (r, c) = V.unsafeIndex (storage x) (indexMatToVec x (r, c))
+    (!) x (r, c)
+        | r < 0 || r >= nRows x || c < 0 || c >= nCols x = error $ "Index " ++ show (r, c) ++ " is out of bounds for matrix with size " ++ show (size x)
+        | otherwise                                      = unsafeIndex x (r, c)
+    (!?) x (r, c)
+        | r < 0 || r >= nRows x || c < 0 || c >= nCols x = Nothing
+        | otherwise                                      = Just $ unsafeIndex x (r, c)
+
+
+instance Num a => Num (Mat a) where
+    (+) = elementwise (+)
+    (-) = elementwise (-)
+    negate = fmap (0 -)
+    (*) = adamarMul
+    abs = fmap abs
+    signum = fmap signum
+    fromInteger = singleton . fromInteger
+
+instance Fractional a => Fractional (Mat a) where
+    (/) = elementwise (/)
+    recip = fmap (1/)
+    fromRational = singleton . fromRational
+
+instance Floating a => Floating (Mat a) where
+    pi = singleton pi
+    (**) = elementwise (**)
+    sqrt = fmap sqrt
+    exp = fmap exp
+    log = fmap log
+    sin = fmap sin
+    cos = fmap cos
+    asin = fmap asin
+    acos = fmap acos
+    atan = fmap atan
+    sinh = fmap sinh
+    cosh = fmap cosh
+    asinh = fmap asinh
+    acosh = fmap acosh
+    atanh = fmap atanh
+
+
+instance ElementwiseScalarOps (Mat a) where
+    (+.) x n = fmap (+ n) x
+    (-.) x n = fmap (subtract n) x
+    (*.) x n = fmap (* n) x
+    (/.) x n = fmap (/ n) x
+    (**.) x n = fmap (** n) x
+
+    elementsMin x n = fmap (min n) x
+    elementsMax x n = fmap (max n) x
+
+instance SingletonOps (Mat a) where
+    singleton x =  Mat 1 1 1 1 0 0 (V.singleton x)
+    isSingleton mat = nElements mat == 1
+    unSingleton mat
+        | not $ isSingleton mat = error "Matrix is not a singleton"
+        | otherwise             = unsafeIndex mat (0, 0)
+
+    extendSingleton mat reference = extend mat (unSingleton mat) (size reference)
+
+    elementsSum mat@(Mat rows cols _ _ _ _ _) = singleton $ sum [unsafeIndex mat (r, c) | r <- [0 .. rows - 1], c <- [0 .. cols - 1]]
+    elementsProduct mat@(Mat rows cols _ _ _ _ _) = singleton $ product [unsafeIndex mat (r, c) | r <- [0 .. rows - 1], c <- [0 .. cols - 1]]
+    mean x = elementsSum x /. fromIntegral (nElements x)
+    norm x = sqrt $ elementsSum $ x * x
+
+
+instance Eq a => Eq (Mat a) where
+    (==) a@(Mat rows1 cols1 _ _ _ _ _) b@(Mat rows2 cols2 _ _ _ _ _)
+        | rows1 /= rows2 || cols1 /= cols2 = False
+        | otherwise                        = and [unsafeIndex a (r, c) == unsafeIndex b (r, c) | r <- [0 .. rows1 - 1], c <- [0 .. cols1 - 1]]
+
+
+instance Functor Mat where
+    fmap f (Mat rows cols rk ck r0 c0 x) = Mat rows cols rk ck r0 c0 (fmap f x)
+    (<$) = fmap . const
+
+instance Applicative Mat where
+    pure x = Mat 1 1 1 1 0 0 (V.singleton x)
+    (<*>) = elementwise (\f x -> f x)
+
+instance Foldable Mat where
+    foldr f x = V.foldr f x . storage
+    foldl f x = V.foldl f x . storage
+
+    foldr' f x = V.foldr' f x . storage
+    foldl' f x = V.foldl' f x . storage
+
+    foldr1 f = V.foldr1 f . storage
+    foldl1 f = V.foldl1 f . storage
+
+    toList = V.toList . storage
+
+    null x = nElements x == 0
+
+    length = nElements
+
+instance Traversable Mat where
+    sequenceA mat@(Mat rows cols _ _ _ _ _) = fmap (Mat rows cols cols 1 0 0) . sequenceA . storage $ force mat
+
+
+-- Constructors
+
+-- | Creates empty 'Mat'.
+empty :: Mat a
+empty = Mat 0 0 0 0 0 0 V.empty
+
+-- | Creates 'Mat' from list (will throw an error, if elements of that list do not form a matrix of given size).
+fromList :: (Int, Int) -> [a] -> Mat a
+fromList (rows, cols) xs
+    | V.length m /= rows * cols = error "Given dimensions do not match with list length"
+    | otherwise                 = Mat rows cols cols 1 0 0 m
+  where
+    m = V.fromList xs
+
+-- | Creates 'Mat' from list of lists (alias for @fromLists (rows, cols) (concat xs)@).
+fromLists :: (Int, Int) -> [[a]] -> Mat a
+fromLists (rows, cols) xs = fromList (rows, cols) (concat xs)
+
+-- | Creates 'Mat' of given size using generating function.
+generate :: (Int, Int) -> ((Int, Int) -> a) -> Mat a
+generate (rows, cols) f = Mat rows cols cols 1 0 0 (V.generate (rows * cols) (f . flip quotRem cols))
+
+-- | Creates 'Mat' of given size filled with given element.
+replicate :: (Int, Int) -> a -> Mat a
+replicate (rows, cols) x = Mat rows cols cols 1 0 0 (V.replicate (rows * cols) x)
+
+
+-- Vec operations
+
+-- | Converts 'Synapse.Tensors.Vec.Vec' to a one row 'Mat'.
+rowVec :: Vec a -> Mat a
+rowVec (Vec x) = Mat 1 (V.length x) (V.length x) 1 0 0 x
+
+-- | Converts 'Synapse.Tensors.Vec.Vec' to a one column 'Mat'.
+colVec :: Vec a -> Mat a
+colVec (Vec x) = Mat (V.length x) 1 1 1 0 0 x
+
+-- | Initializes 'Mat' from given 'Synapse.Tensors.Vec.Vec'.
+fromVec :: (Int, Int) -> Vec a -> Mat a
+fromVec (rows, cols) (Vec x)
+    | V.length x /= rows * cols = error "Given dimensions do not match with vector length"
+    | otherwise                 = Mat rows cols cols 1 0 0 x
+
+
+-- | Extracts row from 'Mat'. If row is not present, an error is thrown.
+indexRow :: Mat a -> Int -> Vec a
+indexRow mat@(Mat rows cols _ _ _ _ x) r
+    | r < 0 || r >= rows = error "Given row is not present in the matrix"
+    | isTransposed mat   = Vec $ V.fromList [unsafeIndex mat (r, c) | c <- [0 .. cols - 1]]
+    | otherwise          = Vec $ V.slice (indexMatToVec mat (r, 0)) cols x
+
+-- | Extracts column from 'Mat'. If column is not present, an error is thrown.
+indexCol :: Mat a -> Int -> Vec a
+indexCol mat@(Mat rows cols _ _ _ _ x) c
+    | c < 0 || c >= cols = error "Given column is not present in the matrix"
+    | isTransposed mat   = Vec $ V.slice (indexMatToVec mat (0, c)) rows x
+    | otherwise          = Vec $ V.fromList [unsafeIndex mat (r, c) | r <- [0 .. rows - 1]]
+
+-- | Extracts row from 'Mat'.
+safeIndexRow :: Mat a -> Int -> Maybe (Vec a)
+safeIndexRow mat@(Mat rows cols _ _ _ _ x) r
+    | r < 0 || r >= rows = Nothing
+    | isTransposed mat   = Just $ Vec $ V.fromList [unsafeIndex mat (r, c) | c <- [0 .. cols - 1]]
+    | otherwise          = Just $ Vec $ V.slice (indexMatToVec mat (r, 0)) cols x
+
+-- | Extracts column from 'Mat'.
+safeIndexCol :: Mat a -> Int -> Maybe (Vec a)
+safeIndexCol mat@(Mat rows cols _ _ _ _ x) c
+    | c < 0 || c >= cols = Nothing
+    | isTransposed mat   = Just $ Vec $ V.slice (indexMatToVec mat (0, c)) rows x
+    | otherwise          = Just $ Vec $ V.fromList [unsafeIndex mat (r, c) | r <- [0 .. rows - 1]]
+
+
+-- | Extracts diagonal from 'Mat'.
+diagonal :: Mat a -> Vec a
+diagonal x = Vec $ V.fromList [unsafeIndex x (n, n) | n <- [0 .. min (nRows x) (nCols x) - 1]]
+
+-- | Flattens 'Mat' to a 'Vec'.
+flatten :: Mat a -> Vec a
+flatten = Vec . storage . force
+
+
+-- Combining
+
+-- | Applies function to every element of 'Mat'.
+map :: (a -> b) -> Mat a -> Mat b
+map = fmap
+
+-- | Applies function to a given row. If new 'Synapse.Tensors.Vec.Vec' is longer then кщц, it is truncated.
+mapRow :: Int -> (Vec a -> Vec a) -> Mat a -> Mat a
+mapRow row f mat = let newRow = f $ indexRow mat row
+                   in imap (\(r, c) x -> if r == row then unsafeIndex newRow c else x) mat
+
+-- | Applies function to a given column. If new 'Synapse.Tensors.Vec.Vec' is longer then column, it is truncated.
+mapCol :: Int -> (Vec a -> Vec a) -> Mat a -> Mat a
+mapCol col f mat = let newCol = f $ indexCol mat col
+                   in imap (\(r, c) x -> if c == col then unsafeIndex newCol r else x) mat
+
+-- | Flipped 'map'.
+for :: Mat a -> (a -> b) -> Mat b
+for = flip map
+
+-- | Applies function to every element and its position of 'Mat'.
+imap :: ((Int, Int) -> a -> b) -> Mat a -> Mat b
+imap f mat@(Mat rows cols rk ck r0 c0 x) = Mat rows cols rk ck r0 c0 (V.imap (f . indexVecToMat mat) x)
+
+-- | Zips two 'Mat's together using given function.
+elementwise :: (a -> b -> c) -> Mat a -> Mat b -> Mat c
+elementwise f a b
+    | size a /= size b = error "Two matrices have different sizes"
+    | otherwise        = let (rows, cols) = size a
+                         in Mat rows cols cols 1 0 0 $
+                            V.fromList [f (unsafeIndex a (r, c)) (unsafeIndex b (r, c)) | r <- [0 .. rows - 1], c <- [0 .. cols - 1]]
+
+
+-- Operations with matrices
+
+-- | Sets new size for a matrix relative to top left corner and uses given element for new entries if the matrix is extended.
+setSize :: Mat a -> a -> (Int, Int) -> Mat a
+setSize mat x = flip generate $ \(r, c) -> if r < nRows mat && c < nCols mat then unsafeIndex mat (r, c) else x
+
+-- | Extends matrix size relative to top left corner using given element for new entries. The matrix is never reduced in size.
+extend :: Mat a -> a -> (Int, Int) -> Mat a
+extend mat x (rows, cols) = setSize mat x (max (nRows mat) rows, max (nCols mat) cols)
+
+-- | Shrinks matrix size relative to top left corner. The matrix is never extended in size.
+shrink :: Mat a -> (Int, Int) -> Mat a
+shrink mat (rows, cols) = setSize mat undefined (min (nRows mat) rows, min (nCols mat) cols)
+
+
+-- | Swaps two rows.
+swapRows :: Mat a -> Int -> Int -> Mat a
+swapRows mat@(Mat rows cols _ _ _ _ _) row1 row2
+    | row1 < 0 || row2 < 0 || row1 >= rows || row2 >= rows = error "Given row indices are out of bounds"
+    | otherwise                                            = generate (rows, cols) $ \(r, c) -> if r == row1 then unsafeIndex mat (row2, c)
+                                                                                                else if r == row2 then unsafeIndex mat (row1, c)
+                                                                                                else unsafeIndex mat (r, c)
+-- | Swaps two columns.
+swapCols :: Mat a -> Int -> Int -> Mat a
+swapCols mat@(Mat rows cols _ _ _ _ _) col1 col2
+    | col1 < 0 || col2 < 0 || col1 >= cols || col2 >= cols = error "Given column indices are out of bounds"
+    | otherwise                                            = generate (rows, cols) $ \(r, c) -> if c == col1 then unsafeIndex mat (r, col2)
+                                                                                                else if c == col2 then unsafeIndex mat (r, col1)
+                                                                                                else unsafeIndex mat (r, c)
+
+
+-- Submatrices
+
+-- | Extacts minor matrix, skipping given row and column.
+minor :: Mat a -> (Int, Int) -> Mat a
+minor mat@(Mat rows cols _ _ _ _ _) (r', c') = Mat (rows - 1) (cols - 1) (cols - 1) 1 0 0 $
+                                             V.fromList [unsafeIndex mat (r, c) | r <- [0 .. rows - 1], c <- [0 .. cols - 1], r /= r' && c /= c']
+
+-- | Extracts submatrix, that is located between given two positions.
+submatrix :: Mat a -> ((Int, Int), (Int, Int)) -> Mat a
+submatrix (Mat rows cols rk ck r0 c0 x) ((r1, c1), (r2, c2))
+    | r1 < 0 || c1 < 0 || r2 < 0 || c2 < 0 ||
+      r2 < r1 || r2 > rows || c2 < c2 || c2 > cols = error "Given row and column limits are incorrect"
+    | otherwise                                    = Mat (r2 - r1) (c2 - c1) rk ck (r0 + r1) (c0 + c1) x
+
+-- | Splits matrix into 4 parts, given position is a pivot, that corresponds to first element of bottom-right subpart.
+split :: Mat a -> (Int, Int) -> (Mat a, Mat a, Mat a, Mat a)
+split x (r, c) = (submatrix x ((0, 0), (r, c)),       submatrix x ((0, c), (r, nCols x)),
+                  submatrix x ((r, 0), (nRows x, c)), submatrix x ((r, c), (nRows x, nCols x)))
+
+-- | Joins 4 blocks of matrices.
+join :: (Mat a, Mat a, Mat a, Mat a) -> Mat a
+join (tl@(Mat rowsTL colsTL _ _ _ _ _), tr@(Mat rowsTR colsTR _ _ _ _ _),
+      bl@(Mat rowsBL colsBL _ _ _ _ _), br@(Mat rowsBR colsBR _ _ _ _ _))
+    | rowsTL /= rowsTR || rowsBL /= rowsBR ||
+      colsTL /= colsBL || colsTR /= colsBR = error "Matrices dimensions do not match"
+    | otherwise                            = generate (rowsTL + rowsBL, colsTL + colsTR) $
+                                             \(r, c) -> uncurry unsafeIndex $
+                                                        if r >= rowsTL && c >= colsTL then (br, (r - rowsTL, c - colsTL))
+                                                        else if r >= rowsTL           then (bl, (r - rowsTL, c))
+                                                        else if c >= colsTL           then (tr, (r, c - colsTL))
+                                                        else                               (tl, (r, c))
+
+infixl 9 <|>
+-- | Joins two matrices horizontally.
+(<|>) :: Mat a -> Mat a -> Mat a
+(<|>) a@(Mat rows1 cols1 _ _ _ _ _) b@(Mat rows2 cols2 _ _ _ _ _)
+    | rows1 /= rows2 = error "Given matrices must have the same number of rows"
+    | otherwise      = generate (rows1, cols1 + cols2) $ \(r, c) -> if c < cols1
+                                                                    then unsafeIndex a (r, c)
+                                                                    else unsafeIndex b (r, c - cols1)
+
+infixl 9 <->
+-- | Joins two matrices vertically.
+(<->) :: Mat a -> Mat a -> Mat a
+(<->) a@(Mat rows1 cols1 _ _ _ _ _) b@(Mat rows2 cols2 _ _ _ _ _)
+    | cols1 /= cols2 = error "Given matrices must have the same number of columns"
+    | otherwise      = generate (rows1 + rows2, cols1) $ \(r, c) -> if r < rows1
+                                                                    then unsafeIndex a (r, c)
+                                                                    else unsafeIndex b (r - rows1, c)
+
+
+-- Functions that work on mathematical matrix (type constraint refers to a number)
+
+-- | Creates 'Mat' that is filled with zeroes.
+zeroes :: Num a => (Int, Int) -> Mat a
+zeroes = flip replicate 0
+
+-- | Creates 'Mat' that is filled with ones.
+ones :: Num a => (Int, Int) -> Mat a
+ones = flip replicate 1
+
+-- | Creates identity matrix.
+identity :: Num a => Int -> Mat a
+identity n = generate (n, n) $ \(r, c) -> if r == c then 1 else 0
+
+
+-- | Adamar multiplication (elementwise multiplication).
+adamarMul :: Num a => Mat a -> Mat a -> Mat a
+adamarMul = elementwise (*)
+
+instance Num a => MatOps (Mat a) where
+    transpose (Mat rows cols rk ck r0 c0 x) = Mat cols rows ck rk c0 r0 x
+    addMatRow mat row
+        | nRows row /= 1         = error "Given row matrix is not a row"
+        | nCols row /= nCols mat = error "Number of columns does not match"
+        | otherwise              = imap (\(_, c) -> (+ unsafeIndex row (0, c))) mat
+    matMul a@(Mat rows1 cols1 _ _ _ _ _) b@(Mat rows2 cols2 _ _ _ _ _)
+        | cols1 /= rows2 = error "Matrices dimensions do not match"
+        | otherwise      = generate (rows1, cols2) $ \(r, c) -> unSingleton $ indexRow a r `SV.dot` indexCol b c
+
+-- | Determinant of a square matrix. If matrix is empty, zero is returned.
+det :: Num a => Mat a -> a
+det mat@(Mat rows cols _ _ _ _ _)
+    | rows /= cols       = error "Matrix is not square"
+    | nElements mat == 0 = 0
+    | nElements mat == 1 = unsafeIndex mat (0, 0)
+    | otherwise          = let mat' = if isTransposed mat then transpose mat else mat
+
+                               determinant :: Num a => Mat a -> a
+                               determinant x@(Mat 2 _ _ _ _ _ _) = unsafeIndex x (0, 0) * unsafeIndex x (1, 1) -
+                                                                   unsafeIndex x (1, 0) * unsafeIndex x (0, 1)
+                               determinant x                     = sum [((-1) ^ i) * (indexRow x 0 ! i) * determinant (minor x (0, i))
+                                                                      | i <- [0 .. nRows x - 1]]
+                           in determinant mat'
+
+-- | Row reduced echelon form of matrix.
+rref :: (Eq a, Fractional a) => Mat a -> Mat a
+rref mat@(Mat rows cols _ _ _ _ _) = go mat 0 [0 .. rows - 1]
+  where
+    go m _ [] = m
+    go m lead (r:rs) = case find ((0 /=) . unsafeIndex m) [(i, j) | j <- [lead .. cols - 1], i <- [r .. rows - 1]] of
+                           Nothing                -> m
+                           Just (pivotRow, lead') -> let newRow = SV.map (/ unsafeIndex m (pivotRow, lead')) (indexRow m pivotRow)
+                                                         m'   = swapRows m pivotRow r
+                                                         m''  = mapRow r (const newRow) m'
+                                                         m''' = imap (\(row, c) -> if row == r
+                                                                                   then id
+                                                                                   else subtract (newRow ! c * unsafeIndex m'' (row, lead'))
+                                                                     ) m''
+                                                     in go m''' (lead' + 1) rs
+
+-- | Inverse of a square matrix. If given matrix is empty, empty matrix is returned.
+inverse :: (Eq a, Fractional a) => Mat a -> Maybe (Mat a)
+inverse mat@(Mat rows cols _ _ _ _ _)
+    | rows /= cols       = error "Matrix is not square"
+    | nElements mat == 0 = Just empty
+    | otherwise          = let mat' = mat <|> identity rows
+                               reduced = rref mat'
+                               (left, right, _, _) = split reduced (rows, cols)
+                           in case V.find (== 0) (SV.unVec $ diagonal left) of
+                                  Nothing -> Just right
+                                  Just _  -> Nothing
+
+-- | Orthogonalizes matrix by rows using Gram-Schmidt algorithm.
+orthogonalized :: Floating a => Mat a -> Mat a
+orthogonalized mat = foldl' (\mat' row -> mapRow row SV.normalized $ iterationGramSchmidt mat' row) mat [0 .. nRows mat]
+  where
+    iterationGramSchmidt mat' row = foldl' (\mat'' r ->
+                                            mapRow row (subtract $ indexRow mat'' r *. unSingleton (indexRow mat'' r `SV.dot` indexRow mat'' row)) mat''
+                                           ) mat' [0 .. row]
+ src/Synapse/Tensors/Vec.hs view
@@ -0,0 +1,295 @@+{- | Implementation of mathematical vector.
+
+'Vec' is only a newtype wrapper around 'Data.Vector', which
+implements several mathematical operations on itself.
+
+'Vec' offers meaningful abstraction and easy interface
+(you can unwrap it to perform more complex tasks).
+-}
+
+
+-- 'TypeFamilies' are needed to instantiate 'Indexable', 'ElementwiseScalarOps', 'SingletonOps', 'VecOps' typeclasses.
+{-# LANGUAGE TypeFamilies #-}
+
+
+module Synapse.Tensors.Vec
+    ( -- * 'Vec' datatype and simple getters.
+
+      Vec (Vec, unVec)
+
+    , size
+
+      -- * Constructors
+
+    , empty
+    , singleton
+    , fromList
+    , generate
+    , replicate
+
+      -- * Concatenation and splitting
+
+    , cons
+    , snoc
+    , (++)
+    , concat
+
+    , splitAt
+
+      -- * Combining
+
+    , map
+    , imap
+    , for
+    , zipWith
+    , zip
+
+    -- * Mathematics
+
+    , zeroes
+    , ones
+
+    , squaredMagnitude
+    , magnitude
+    , clampMagnitude
+    , normalized
+
+    , linearCombination
+    , dot
+
+    , angleBetween
+    , lerp
+    ) where
+
+
+import Synapse.Tensors (DType, Indexable(..), ElementwiseScalarOps(..), SingletonOps(..), VecOps(..))
+
+import Prelude hiding ((++), concat, splitAt, map, replicate, zip, zipWith)
+import Data.Foldable (Foldable(..))
+import Data.Ord (clamp)
+
+import qualified Data.Vector as V
+
+
+-- | Mathematical vector (collection of elements).
+newtype Vec a = Vec
+    { unVec :: V.Vector a  -- ^ Internal representation.
+    } deriving (Eq, Read)
+
+-- | Size of a vector - number of elements.
+size :: Vec a -> Int
+size = V.length . unVec
+
+
+-- Typeclasses
+
+instance Show a => Show (Vec a) where
+    show (Vec x) = show x
+
+
+type instance DType (Vec a) = a
+
+
+instance Indexable (Vec a) where
+    type Index (Vec a) = Int
+
+    unsafeIndex (Vec x) = V.unsafeIndex x
+    (!) (Vec x) = (V.!) x
+    (!?) (Vec x) = (V.!?) x
+
+
+instance Num a => Num (Vec a) where
+    (+) = zipWith (+)
+    (-) = zipWith (-)
+    negate = fmap (0 -)
+    (*) = zipWith (*)
+    abs = fmap abs
+    signum = fmap signum
+    fromInteger = singleton . fromInteger
+
+instance Fractional a => Fractional (Vec a) where
+    (/) = zipWith (/)
+    recip = fmap (1/)
+    fromRational = singleton . fromRational
+
+instance Floating a => Floating (Vec a) where
+    pi = singleton pi
+    (**) = zipWith (**)
+    sqrt = fmap sqrt
+    exp = fmap exp
+    log = fmap log
+    sin = fmap sin
+    cos = fmap cos
+    asin = fmap asin
+    acos = fmap acos
+    atan = fmap atan
+    sinh = fmap sinh
+    cosh = fmap cosh
+    asinh = fmap asinh
+    acosh = fmap acosh
+    atanh = fmap atanh
+
+instance ElementwiseScalarOps (Vec a) where
+    (+.) x n = fmap (+ n) x
+    (-.) x n = fmap (subtract n) x
+    (*.) x n = fmap (* n) x
+    (/.) x n = fmap (/ n) x
+    (**.) x n = fmap (** n) x
+
+    elementsMin x n = fmap (min n) x
+    elementsMax x n = fmap (max n) x
+
+instance SingletonOps (Vec a) where
+    singleton = pure
+    isSingleton vec = size vec == 1
+    unSingleton vec
+        | not $ isSingleton vec = error "Vector is not a singleton"
+        | otherwise             = unsafeIndex vec 0
+
+    extendSingleton vec reference = replicate (size reference) (unSingleton vec)
+
+    elementsSum = singleton . V.sum . unVec
+    elementsProduct = singleton . V.product . unVec
+    mean x = elementsSum x /. fromIntegral (size x)
+    norm x = sqrt $ elementsSum $ x * x
+
+instance Functor Vec where
+    fmap f = Vec . V.map f . unVec
+    (<$) = fmap . const
+
+instance Applicative Vec where
+    pure = Vec . V.singleton
+    (<*>) = zipWith (\f x -> f x)
+
+instance Foldable Vec where
+    foldr f x = V.foldr f x . unVec
+    foldl f x = V.foldl f x . unVec
+
+    foldr' f x = V.foldr' f x . unVec
+    foldl' f x = V.foldl' f x . unVec
+
+    foldr1 f = V.foldr1 f . unVec
+    foldl1 f = V.foldl1 f . unVec
+
+    toList = V.toList . unVec
+
+    null x = size x == 0
+
+    length = size
+
+instance Traversable Vec where
+    traverse f (Vec x) = Vec <$> traverse f x
+
+
+-- Constructors
+
+-- | Creates empty 'Vec'.
+empty :: Vec a
+empty = Vec V.empty
+
+-- | Creates 'Vec' from list.
+fromList :: [a] -> Vec a
+fromList = Vec . V.fromList
+
+-- | Creates 'Vec' of given length using generating function.
+generate :: Int -> (Int -> a) -> Vec a
+generate n = Vec . V.generate n
+
+-- | Creates 'Vec' of given length filled with given element.
+replicate :: Int -> a -> Vec a
+replicate n = generate n . const
+
+
+-- Concatenation and splitting
+
+-- | Prepend 'Vec' with given element.
+cons :: a -> Vec a -> Vec a
+cons x = Vec . V.cons x . unVec
+
+-- | Append 'Vec' with given element.
+snoc :: Vec a -> a -> Vec a
+snoc (Vec vec) x = Vec $ V.snoc vec x
+
+-- | Concatenate two 'Vec's.
+infixr 5 ++
+(++) :: Vec a -> Vec a -> Vec a
+(++) (Vec x) (Vec y) = Vec $ x V.++ y
+
+-- | Concatenate all 'Vec's.
+concat :: [Vec a] -> Vec a
+concat = foldr1 (++)
+
+-- | Splits 'Vec' into two 'Vec's at a given index.
+splitAt :: Int -> Vec a -> (Vec a, Vec a)
+splitAt i (Vec v) = let (v1, v2) = V.splitAt i v
+                    in (Vec v1, Vec v2)
+
+
+-- Combining
+
+-- | Map a function over a 'Vec'.
+map :: (a -> b) -> Vec a -> Vec b
+map = fmap
+
+-- | Apply a function to every element of a 'Vec' and its index.
+imap :: (Int -> a -> b) -> Vec a -> Vec b
+imap f = Vec . V.imap f . unVec
+
+-- | 'map' with its arguments flipped.
+for :: Vec a -> (a -> b) -> Vec b
+for = flip fmap
+
+-- | Zips two 'Vec's with the given function.
+zipWith :: (a -> b -> c) -> Vec a -> Vec b -> Vec c
+zipWith f (Vec a) (Vec b) = Vec $ V.zipWith f a b
+
+-- | Zips two 'Vec's.
+zip :: Vec a -> Vec b -> Vec (a, b)
+zip = zipWith (,)
+
+
+-- Functions that work on mathematical vector (type constraint refers to a number)
+
+-- | Creates 'Vec' of given length filled with zeroes.
+zeroes :: Num a => Int -> Vec a
+zeroes = flip generate (const 0)
+
+-- | Creates 'Vec' of given length filled with ones.
+ones :: Num a => Int -> Vec a
+ones = flip generate (const 1)
+
+
+-- | Squared magnitude of a 'Vec'.
+squaredMagnitude :: Num a => Vec a -> a
+squaredMagnitude x = sum (fmap (^ (2 :: Int)) x)
+
+-- | Magnitude of a 'Vec'.
+magnitude :: Floating a => Vec a -> a
+magnitude = sqrt . squaredMagnitude
+
+-- | Clamps 'Vec' magnitude.
+clampMagnitude :: (Floating a, Ord a) => a -> Vec a -> Vec a
+clampMagnitude m x = x *. (min (magnitude x) m / magnitude x)
+
+-- | Normalizes 'Vec' by dividing each component by 'Vec' magnitude.
+normalized :: Floating a => Vec a -> Vec a
+normalized x = x /. magnitude x
+
+
+-- | Computes linear combination of 'Vec's. Returns empty 'Vec' if empty list was passed to this function.
+linearCombination :: Num a => [(a, Vec a)] -> Vec a
+linearCombination [] = empty
+linearCombination (x:xs) = foldl' (\acc (a, vec) -> acc + vec *. a) (snd x *. fst x) xs
+
+
+instance Num a => VecOps (Vec a) where
+    dot a b = elementsSum $ a * b
+
+
+-- | Calculates an angle between two 'Vec's.
+angleBetween :: Floating a => Vec a -> Vec a -> a
+angleBetween a b = acos $ unSingleton (a `dot` b) / (magnitude a * magnitude b)
+
+-- | Linearly interpolates between two 'Vec's. Given parameter will be clamped between [0.0, 1.0].
+lerp :: (Floating a, Ord a) => a -> Vec a -> Vec a -> Vec a
+lerp k a b = b - (b - a) *. clamp (0.0, 1.0) k
+ test/AutogradTest.hs view
@@ -0,0 +1,78 @@+-- | Tests @Synapse.Autograd@ module and its submodules.
+
+
+module AutogradTest
+    ( tests
+    ) where
+
+
+import Synapse.Autograd
+
+import Synapse.Tensors (MatOps(..))
+import qualified Synapse.Tensors.Vec as V
+import qualified Synapse.Tensors.Mat as M
+
+import Test.HUnit
+
+
+testNumOps :: Test
+testNumOps = TestLabel "testIntOps" $ TestList
+    [ TestCase $ assertEqual "zero gradient" 0 (unSymbol $ getGradientsOf a `wrt` b)
+    , TestCase $ assertEqual "identity gradient" 1 (unSymbol $ getGradientsOf a `wrt` a)
+    , TestCase $ assertEqual "composited identity gradient" 1 (unSymbol $ getGradientsOf c `wrt` c)
+    , TestCase $ assertEqual "addition gradient" 1 (unSymbol $ getGradientsOf (a + b) `wrt` a)
+    , TestCase $ assertEqual "subtraction gradient" (-1) (unSymbol $ getGradientsOf (a - b) `wrt` b)
+    , TestCase $ assertEqual "multiplication gradient" 3 (unSymbol $ getGradientsOf (a * b) `wrt` a)
+    , TestCase $ assertEqual "composed operations gradient" (4 * 4 * 3) (unSymbol $ getGradientsOf (a * (b + b) * a) `wrt` a)
+    , TestCase $ assertEqual "renamed symbol gradient" (2 * 25) (unSymbol $ getGradientsOf (c * c) `wrt` c)
+    ]
+  where
+    a = symbol (SymbolIdentifier "a") 4 :: Symbol Int
+    b = symbol (SymbolIdentifier "b") 3 :: Symbol Int
+    c = renameSymbol (SymbolIdentifier "c") ((a * a) + (b * b)) :: Symbol Int
+
+testNthOrderGradients :: Test
+testNthOrderGradients = TestLabel "testNthOrderGradients" $ TestList
+    [ TestCase $ assertEqual "identity 2nd gradient" 0 (unSymbol $ nthGradient 2 a a)
+    , TestCase $ assertEqual "composed operations 2nd gradient" (4 * 3) (unSymbol $ nthGradient 2 (a * (b + b) * a) a)
+    , TestCase $ assertEqual "4th gradient" (unSymbol $ sin c) (unSymbol $ nthGradient 4 (sin c) c)
+    ]
+  where
+    a = symbol (SymbolIdentifier "a") 4 :: Symbol Int
+    b = symbol (SymbolIdentifier "b") 3 :: Symbol Int
+    c = symbol (SymbolIdentifier "c") 0.5 :: Symbol Float
+
+testVecOps :: Test
+testVecOps = TestLabel "testVecOps" $ TestList
+    [ TestCase $ assertEqual "vector addition gradient" (V.replicate 3 1.0) (unSymbol $ getGradientsOf (a + b) `wrt` a)
+    , TestCase $ assertEqual "vector elementwise multiplication gradient" (unSymbol b) (unSymbol $ getGradientsOf (a * b) `wrt` a)
+    , TestCase $ assertEqual "vector op gradient" (V.map cos $ unSymbol a) (unSymbol $ getGradientsOf (sin a) `wrt` a)
+    ]
+  where
+    a = symbol (SymbolIdentifier "a") (V.fromList [1.0, 2.0, 3.0]) :: SymbolVec Float
+    b = symbol (SymbolIdentifier "b") (V.fromList [3.0, 2.0, 1.0]) :: SymbolVec Float
+
+testMatOps :: Test
+testMatOps = TestLabel "testMatOps" $ TestList
+    [ TestCase $ assertEqual "matrix addition gradient" (M.replicate (3, 3) 1.0) (unSymbol $ getGradientsOf (a + b) `wrt` a)
+    , TestCase $ assertEqual "matrix elementwise multiplication gradient" (unSymbol b) (unSymbol $ getGradientsOf (a * b) `wrt` a)
+    , TestCase $ assertEqual "matrix op gradient" (M.map cos $ unSymbol a) (unSymbol $ getGradientsOf (sin a) `wrt` a)
+    , TestCase $ assertEqual "matrix transposing gradient" (M.transpose $ unSymbol c) (unSymbol $ getGradientsOf (transpose c) `wrt` c)
+    , TestCase $ assertEqual "matrix multiplication+addition gradient" (M.transpose $ unSymbol c) (unSymbol $ getGradientsOf (c `matMul` d + e) `wrt` d)
+    , TestCase $ assertEqual "matrix composed multiplication gradient" (M.transpose $ unSymbol c) (unSymbol $ getGradientsOf (c `matMul` d `matMul` e) `wrt` d)
+    ]
+  where
+    a = symbol (SymbolIdentifier "a") (M.replicate (3, 3) 3.0) :: SymbolMat Float
+    b = symbol (SymbolIdentifier "b") (M.replicate (3, 3) (-3.0)) :: SymbolMat Float
+    c = symbol (SymbolIdentifier "c") (M.fromLists (1, 2) [[5.0, 3.0]]) :: SymbolMat Float
+    d = symbol (SymbolIdentifier "d") (M.fromLists (2, 1) [[1.0], [-2.0]]) :: SymbolMat Float
+    e = symbol (SymbolIdentifier "e") (M.fromLists (1, 1) [[1.0]]) :: SymbolMat Float
+
+
+tests :: Test
+tests = TestLabel "AutogradTest" $ TestList
+    [ testNumOps
+    , testNthOrderGradients
+    , testVecOps
+    , testMatOps
+    ]
+ test/Main.hs view
@@ -0,0 +1,19 @@+-- | Test suite for "Synapse" library.
+
+
+module Main (main) where
+
+
+import qualified TensorsTest
+import qualified AutogradTest
+import qualified NNTest
+
+import Test.HUnit (Test(TestList), runTestTTAndExit)
+
+
+main :: IO ()
+main = runTestTTAndExit $ TestList
+    [ TensorsTest.tests
+    , AutogradTest.tests
+    , NNTest.tests
+    ]
+ test/NNTest.hs view
@@ -0,0 +1,111 @@+-- | Tests "Synapse.NN" module and its submodules.
+
+
+module NNTest
+    ( tests
+    )
+    where
+
+
+import Synapse.Tensors
+import qualified Synapse.Tensors.Vec as V
+
+import Synapse.NN.Models
+import Synapse.NN.Layers
+import Synapse.NN.Optimizers
+import Synapse.NN.LearningRates
+import Synapse.NN.Losses
+import Synapse.NN.Batching
+import Synapse.NN.Training
+
+import System.Random
+
+import Graphics.EasyPlot
+
+import Test.HUnit
+
+
+testSin :: Test  -- -3 sin (x + 5)
+testSin = TestLabel "testSin" $ TestCase $ do
+    let sinFn x = (-3.0) * sin (x + 5.0)
+    let model = buildSequentialModel (InputSize 1) [ Layer . layerDense 1
+                                                   , Layer . layerActivation (Activation cos)
+                                                   , Layer . layerDense 1
+                                                   ] :: SequentialModel Double
+
+    let dataset = Dataset $ V.fromList $ [Sample (singleton x) (sinFn $ singleton x) | x <- [-pi, -pi+0.2 .. pi]]
+    (trainedModel, _, losses, _) <- train model
+                                       (SGD 0.2 False)
+                                       (Hyperparameters 500 16 dataset (LearningRate $ const 0.01) (Loss mse) V.empty)
+                                       emptyCallbacks
+                                       (mkStdGen 1)
+    _ <- plot (PNG "test/plots/sin.png")
+              [ Data2D [Title "predicted sin", Style Lines, Color Red] [Range (-pi) pi] [(x, unSingleton $ forward (singleton x) trainedModel) | x <- [-pi, -pi+0.05 .. pi]]
+              , Data2D [Title "true sin", Style Lines, Color Green] [Range (-pi) pi] [(x, sinFn x) | x <- [-pi, -pi+0.05 .. pi]]
+              ]
+
+    let unpackedLosses = unRecordedMetric (unsafeIndex losses 0)
+    let lastLoss = unsafeIndex unpackedLosses (V.size unpackedLosses - 1)
+
+    assertBool "trained well enough" (lastLoss < 0.01)
+
+testSqrt :: Test  -- sqrt(x)
+testSqrt = TestLabel "testSqrt" $ TestCase $ do
+    let sqrtFn x = sqrt x
+    let model = buildSequentialModel (InputSize 1) [ Layer . layerDense 1
+                                                   , Layer . layerActivation (Activation tanh)
+                                                   , Layer . layerDense 1
+                                                   ] :: SequentialModel Double
+
+    let dataset = Dataset $ V.fromList $ [Sample (singleton x) (sqrtFn $ singleton x) | x <- [0.0, 0.2 .. 4.0]]
+    (trainedModel, _, losses, _) <- train model
+                                       (SGD 0.2 True)
+                                       (Hyperparameters 500 16 dataset (LearningRate $ const 0.01) (Loss mse) V.empty)
+                                       emptyCallbacks
+                                       (mkStdGen 1)
+    _ <- plot (PNG "test/plots/sqrt.png")
+               [ Data2D [Title "predicted sqrt", Style Lines, Color Red] [Range 0.0 6.0] $ [(x, unSingleton $ forward (singleton x) trainedModel) | x <- [0.0, 0.05 .. 4.0]]
+               , Data2D [Title "true sqrt", Style Lines, Color Green] [Range 0.0 6.0] $ [(x, sqrtFn x) | x <- [0.0, 0.05 .. 4.0]]
+               ]
+
+    let unpackedLosses = unRecordedMetric (unsafeIndex losses 0)
+    let lastLoss = unsafeIndex unpackedLosses (V.size unpackedLosses - 1)
+
+    assertBool "trained well enough" (lastLoss < 0.01)
+
+testTrigonometry :: Test  -- sin(2.0 * cos(x) + 3.0) + 2.5
+testTrigonometry = TestLabel "testTrigonometry" $ TestCase $ do
+    let trigonometryFn x = sin (2.0 * cos x + 3.0) + 2.5
+    let model = buildSequentialModel (InputSize 1) [ Layer . layerDense 1
+                                                   , Layer . layerActivation (Activation sin)
+                                                   , Layer . layerDense 1
+                                                   , Layer . layerActivation (Activation sin)
+                                                   , Layer . layerDense 1
+                                                   ] :: SequentialModel Double
+
+    let dataset = Dataset $ V.fromList $ [Sample (singleton x) (trigonometryFn $ singleton x) | x <- [-(2.0 * pi),((-(2.0 * pi)) + 0.1)..(2.0 * pi)]]
+                                         
+    (trainedModel, _, losses, _) <- train model
+                                       (SGD 0.3 True)
+                                       (Hyperparameters 1000 1 dataset (LearningRate $ const 0.001) (Loss mse) V.empty)
+                                       emptyCallbacks
+                                       (mkStdGen 1)
+
+    _ <- plot (PNG "test/plots/trigonometry.png")
+               [ Data2D [Title "predicted trigonometry", Style Lines, Color Red] [Range ((-2.0) * pi) (2.0 * pi)] $ [(x, unSingleton $ forward (singleton x) trainedModel) | x <- [-(2.0 * pi),((-(2.0 * pi)) + 0.1)..(2.0 * pi)]]
+               , Data2D [Title "true trigonometry", Style Lines, Color Green] [Range ((-2.0) * pi) (2.0 * pi)] $ [(x, trigonometryFn x) | x <- [-(2.0 * pi),((-(2.0 * pi)) + 0.1)..(2.0 * pi)]]
+               ]
+
+    let unpackedLosses = unRecordedMetric (unsafeIndex losses 0)
+    let lastLoss = unsafeIndex unpackedLosses (V.size unpackedLosses - 1)
+
+    assertBool "trained well enough" (lastLoss < 0.01)
+
+
+tests :: Test
+tests = TestLabel "NNTest" $ TestList
+    [ testSin
+    , testSqrt
+    , testTrigonometry
+    ]
+
+ test/TensorsTest.hs view
@@ -0,0 +1,150 @@+-- | Tests "Synapse.Tensors" module and its submodules.
+
+
+module TensorsTest
+    ( tests
+    ) where
+
+
+import Synapse.Tensors ((!), ElementwiseScalarOps(..))
+
+import Synapse.Tensors.Vec (Vec)
+import qualified Synapse.Tensors.Vec as V
+
+import Synapse.Tensors.Mat (Mat)
+import qualified Synapse.Tensors.Mat as M
+
+import Test.HUnit
+
+
+testVecOps :: Test
+testVecOps = TestLabel "testVecOps" $ TestList
+    [ TestCase $ assertBool "true result == manual operations == linear combination" $
+                 all (V.fromList [-1, 2, 5] ==)
+                     [vec1 *. 2 - vec2, V.linearCombination [(2, vec1), (-1, vec2)]]
+    , TestCase $ assertEqual "addition" (V.replicate 3 4) (vec1 + vec2)
+    , TestCase $ assertEqual "dot multiplication" 10 (vec1 `V.dot` vec2)
+    ]
+  where
+    vec1 = V.fromList [1, 2, 3] :: Vec Int
+    vec2 = V.fromList [3, 2, 1] :: Vec Int
+
+testVecMagnitude :: Test
+testVecMagnitude = TestLabel "testVecMagnitude" $ TestList
+    [ TestCase $ assertEqual "magnitude value" 5.0 (V.magnitude vec)
+    , TestCase $ assertEqual "normalized" (V.fromList [0.6, 0.8]) (V.normalized vec)
+    , TestCase $ assertEqual "clamped magnitude" 0.5 (V.magnitude $ V.clampMagnitude 0.5 vec)
+    ]
+  where
+    vec = V.fromList [3.0, 4.0] :: Vec Float
+
+testVecAngle :: Test
+testVecAngle = TestLabel "testVecAngle" $ TestList
+    [ TestCase $ assertEqual "90-degree angle" (pi / 2.0) (V.angleBetween vec1 vec2)
+    , TestCase $ assertEqual "magnitude independent angle" (pi / 2.0) (V.angleBetween (vec1 *. 3.0) (vec2 *. 4.0))
+    , TestCase $ assertEqual "180-degree angle" pi (V.angleBetween vec3 (negate vec3))
+    , TestCase $ assertEqual "45-degree angle" (pi / 4.0) (V.angleBetween vec4 vec5)
+
+    , TestCase $ assertEqual "lerp" (V.fromList [0.5, 0.5]) (V.lerp 0.5 vec1 vec2)
+    ]
+  where
+    vec1 = V.fromList [1.0, 0.0] :: Vec Float
+    vec2 = V.fromList [0.0, 1.0] :: Vec Float
+    vec3 = V.fromList [3.0, 4.0] :: Vec Float
+    vec4 = V.fromList [1.0, 0.0] :: Vec Float
+    vec5 = V.fromList [1.0, 1.0] :: Vec Float
+
+
+testMatExtracting :: Test
+testMatExtracting = TestLabel "testMatExtracting" $ TestList
+    [ TestCase $ assertBool "transposed indexation" $ and [mat ! (r, c) == matT ! (c, r) | r <- [0 .. 2], c <- [0 .. 2]]
+    , TestCase $ assertBool "get rows"
+                 (  (V.fromList [1, 2, 3] == M.indexRow mat 0)
+                 && (V.fromList [4, 5, 6] == M.indexRow mat 1)
+                 && (V.fromList [7, 8, 9] == M.indexRow mat 2)
+                 )
+    , TestCase $ assertBool "get columns"
+                 (  (V.fromList [1, 4, 7] == M.indexCol mat 0)
+                 && (V.fromList [2, 5, 8] == M.indexCol mat 1)
+                 && (V.fromList [3, 6, 9] == M.indexCol mat 2)
+                 )
+    , TestCase $ assertEqual "transposed transposed" mat (M.transpose matT)
+    , TestCase $ assertEqual "diagonal" (V.fromList [1, 5, 9]) (M.diagonal mat)
+    ]
+  where
+    mat = M.fromLists (3, 3) [[1, 2, 3], [4, 5, 6], [7, 8, 9]] :: Mat Int
+    matT = M.transpose mat :: Mat Int
+
+testMatOps :: Test
+testMatOps = TestLabel "testMatOps" $ TestList
+    [ TestCase $ assertEqual "num addition" (M.fromLists (3, 3) [[2, 3, 4], [5, 6, 7], [8, 9, 10]]) (mat1 +. 1)
+    , TestCase $ assertEqual "num multiplication" (M.fromLists (3, 3) [[2, 4, 6], [8, 10, 12], [14, 16, 18]]) (mat1 *. 2)
+
+    , TestCase $ assertEqual "addition" (M.replicate (3, 3) 10) (mat1 + mat2)
+    , TestCase $ assertEqual "transposed addition" (M.fromLists (3, 3) [[10, 8, 6], [12, 10, 8], [14, 12, 10]]) (mat1 + M.transpose mat2)
+
+    , TestCase $ assertEqual "adamar multiplication" (M.fromLists (3, 3) [[9, 16, 21], [24, 25, 24], [21, 16, 9]]) (mat1 * mat2)
+    , TestCase $ assertEqual "matrix multiplication" (M.fromLists (3, 3) [[30, 24, 18], [84, 69, 54], [138, 114, 90]]) (mat1 `M.matMul` mat2)
+
+    , TestCase $ assertEqual "map" (M.replicate (3, 3) 5) (M.map (const (5 :: Int)) mat1)
+    , TestCase $ assertEqual "imap" (M.fromLists (3, 3) [[1, 3, 5], [5, 7, 9], [9, 11, 13]]) (M.imap (\(r, c) x -> r + c + x) mat1)
+
+    , TestCase $ assertEqual "swap rows" (M.fromLists (3, 3) [[7, 8, 9], [4, 5, 6], [1, 2, 3]]) (M.swapRows mat1 0 2)
+    , TestCase $ assertEqual "swap cols" (M.fromLists (3, 3) [[2, 1, 3], [5, 4, 6], [8, 7, 9]]) (M.swapCols mat1 0 1)
+
+    , TestCase $ assertEqual "extend" (M.fromLists (4, 4) [[1, 2, 3, 0], [4, 5, 6, 0], [7, 8, 9, 0], [0, 0, 0, 0]]) (M.extend mat1 0 (4, 4))
+    , TestCase $ assertEqual "failed extend" mat1 (M.extend mat1 0 (2, 2))
+    , TestCase $ assertEqual "shrink" (M.fromLists (2, 2) [[1, 2], [4, 5]]) (M.shrink mat1 (2, 2))
+    , TestCase $ assertEqual "failed shrink" mat1 (M.shrink mat1 (4, 4))
+    ]
+  where
+    mat1 = M.fromLists (3, 3) [[1, 2, 3], [4, 5, 6], [7, 8, 9]] :: Mat Int
+    mat2 = M.fromLists (3, 3) [[9, 8, 7], [6, 5, 4], [3, 2, 1]] :: Mat Int
+
+testMatSubmatrices :: Test
+testMatSubmatrices = TestLabel "testMatSubmatrices" $ TestList
+    [ TestCase $ assertEqual "minor" (M.fromLists (2, 2) [[1, 3], [7, 9]]) (M.minor mat1 (1, 1))
+    , TestCase $ assertEqual "submatrix" (M.fromLists (2, 2) [[2, 3], [5, 6]]) (M.submatrix mat1 ((0, 1), (2, 3)))
+    , TestCase $ assertEqual "split" (M.replicate (4, 4) 1, M.replicate (4, 4) 2, M.replicate (4, 4) 3, M.replicate (4, 4) 4) (M.split mat2 (4, 4))
+    , TestCase $ assertEqual "join" (M.fromLists (2, 2) [[1, 2], [3, 4]]) (M.join (matTL, matTR, matBL, matBR))
+    , TestCase $ assertEqual "join horizontal" (M.fromLists (1, 2) [[1, 2]]) (matTL M.<|> matTR)
+    , TestCase $ assertEqual "join vertical" (M.fromLists (2, 1) [[1], [3]]) (matTL M.<-> matBL)
+    ]
+  where
+    mat1 = M.fromLists (3, 3) [[1, 2, 3], [4, 5, 6], [7, 8, 9]] :: Mat Int
+    mat2 = M.generate (8, 8) (\(r, c) -> if r >= 4 && c >= 4 then 4 else if r >= 4 then 3 else if c >= 4 then 2 else 1) :: Mat Int
+    matTL = M.singleton 1 :: Mat Int
+    matTR = M.singleton 2 :: Mat Int
+    matBL = M.singleton 3 :: Mat Int
+    matBR = M.singleton 4 :: Mat Int
+
+testMatComplicatedOps :: Test
+testMatComplicatedOps = TestLabel "testMatComplicatedOps" $ TestList
+    [ TestCase $ assertEqual "det == 0" 0 (M.det mat1)
+    , TestCase $ assertEqual "det" (-230) (M.det mat2)
+    , TestCase $ assertEqual "rref" (M.fromLists (2, 4) [[1.0, 0.0, -3.0, -4.0], [0.0, 1.0, 1.0, 1.0]]) (M.rref mat3)
+    , TestCase $ assertEqual "inverse" (Just $ M.fromLists (2, 2) [[-3.0, -4.0], [1.0, 1.0]]) (M.inverse mat4)
+    , TestCase $ assertBool "orthogonal" $ M.fromLists (2, 2) [[-(sqrt 2.0 / 2.0), -(sqrt 2.0 / 2.0)], [sqrt 2.0 / 2.0, -(sqrt 2.0 / 2.0)]] ==
+                                           M.orthogonalized mat5
+    ]
+  where
+    mat1 = M.fromLists (3, 3) [[1, 2, 3], [4, 5, 6], [7, 8, 9]] :: Mat Int
+    mat2 = M.fromLists (3, 3) [[2, 5, 5], [4, -10, 0], [-3, -2, 1]] :: Mat Int
+    mat3 = M.fromLists (2, 4) [[1.0, 4.0, 1.0, 0.0], [-1.0, -3.0, 0.0, 1.0]] :: Mat Float
+    mat4 = M.fromLists (2, 2) [[1.0, 4.0], [-1.0, -3.0]] :: Mat Float
+    mat5 = M.fromLists (2, 2) [[1.0, 1.0], [-2.0, 2.0]] :: Mat Float
+
+
+tests :: Test
+tests = TestLabel "TensorsTest" $ TestList
+    [ -- 'Vec'
+      testVecOps
+    , testVecMagnitude
+    , testVecAngle
+
+      -- 'Mat'
+    , testMatExtracting
+    , testMatOps
+    , testMatSubmatrices
+    , testMatComplicatedOps
+    ]