diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,3 @@
+# Changelog for HABQT
+
+## Unreleased changes
diff --git a/HABQT.cabal b/HABQT.cabal
new file mode 100644
--- /dev/null
+++ b/HABQT.cabal
@@ -0,0 +1,129 @@
+name:           HABQT
+version:        0.1.0.0
+synopsis:       Hierarchical adaptive Bayesian quantum tomography for quantum bits
+homepage:       https://github.com/Belinsky-L-V/HABQT#readme
+bug-reports:    https://github.com/Belinsky-L-V/HABQT/issues
+author:         Leonid Belinsky
+maintainer:     belinsky.leonid@gmail.com
+copyright:      Copyright (c) 2018 Leonid Belinsky
+license:        BSD3
+license-file:   LICENSE
+build-type:     Simple
+cabal-version:  >= 1.10
+category:       Math, Quantum
+description:
+    Extends adaptive Bayesian quantum tomography as described in
+    <https://doi.org/10.1103/PhysRevA.85.052120> by using a hierarchical
+    distribution over density matrices of all possible ranks.
+    .
+    \Includes:
+    .
+    * a Haskell library
+    .
+    * a shared library which provides a C
+      interface to the tomography function
+    .
+    * an executable that simulates
+    tomography of random states and outputs infidelity between true states and
+    mean Bayesian estimates to a file
+    .
+    Please refer to @HABQT-simulation \-\-help@ for executable usage
+    instructions,
+    <https://github.com/Belinsky-L-V/HABQT#readme README on Github>
+    for installation instructions and shared library C ABI description,
+    accompanying Haddock documentation for Haskell API.
+
+
+extra-source-files:
+    ChangeLog.md
+    README.md
+
+source-repository head
+  type: git
+  location: https://github.com/Belinsky-L-V/HABQT
+
+library
+  hs-source-dirs:
+      src
+  ghc-options:  -fPIC
+  build-depends:
+      base >=4.10 && <4.11
+    , hmatrix >=0.18.2 && <0.19
+    , hmatrix-gsl >=0.18 && <0.20
+    , mtl >=2.2.2 && <2.3
+    , mwc-random >=0.13.6 && <0.14
+    , newtype-generics >= 0.5.3 && <0.6
+    , streaming >=0.2.1 && <0.3
+    , utility-ht >=0.0.14 && <0.1
+    , vector >=0.12.0 && <0.13
+    , validation >= 1 && <1.1
+  exposed-modules:
+      HABQTlib
+      HABQTlib.UnsafeAPI
+      HABQTlib.Data
+      HABQTlib.Data.Particle
+      HABQTlib.MeasurementProcessing
+      HABQTlib.RandomStates
+  default-language: Haskell2010
+
+executable HABQT-simulation
+  main-is: Main.hs
+  hs-source-dirs:
+      app
+  ghc-options:  -threaded -rtsopts -with-rtsopts=-N
+  build-depends:
+      HABQT
+    , base >=4.10 && <4.11
+    , optparse-applicative >= 0.14.2 && <0.15
+    , streaming >=0.2.1 && <0.3
+  default-language: Haskell2010
+
+foreign-library HABQT
+  type:
+      native-shared
+  lib-version-info:
+      1:0:0
+  if os(Windows)
+    options: standalone
+  hs-source-dirs:
+      libHABQT
+  c-sources:
+      libHABQT/hsinit.c
+  build-depends:
+      HABQT
+    , base >=4.10 && <4.11
+    , mtl >=2.2.2 && <2.3
+    , mwc-random >=0.13.6 && <0.14
+    , validation >= 1 && <1.1
+    , hmatrix >=0.18.2 && <0.19
+    , vector >=0.12.0 && <0.13
+  other-modules:
+      LibHABQT
+      ForeignHABQT
+  default-language: Haskell2010
+
+test-suite HABQT-test
+  type: exitcode-stdio-1.0
+  main-is: Tests.hs
+  hs-source-dirs:
+      test
+  ghc-options:  -threaded -rtsopts -with-rtsopts=-N
+  build-depends:
+      QuickCheck >=2.10.1 && <2.11
+    , HABQT
+    , base >=4.10 && <4.11
+    , hmatrix >=0.18.2 && <0.19
+    , mwc-random >=0.13.6 && <0.14
+    , streaming >=0.2.1 && <0.3
+    , utility-ht >=0.0.14 && <0.1
+    , vector >=0.12.0 && <0.13
+    , newtype-generics >= 0.5.3 && <0.6
+  other-modules:
+      TestHelpers
+      FidelityTests
+      MeasurementTests
+      ParticleProcessingTests
+      RankReductionTests
+      StateGenTests
+      SuperpositionSemigroupTests
+  default-language: Haskell2010
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Leonid Belinsky (c) 2018
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Leonid Belinsky nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,63 @@
+# HABQT
+
+Aim of the project is to extend adaptive Bayesian quantum tomography as
+described in [2012 paper by Huszár and
+Houlsby](https://doi.org/10.1103/PhysRevA.85.052120) by using a hierarchical
+distribution over density matrices of all possible ranks.
+
+Includes:
+
+* a Haskell library
+* a shared library which provides a C interface to the tomography function
+* an executable that simulates tomography of random states and outputs
+  infidelity between true states and mean Bayesian estimates to a file
+
+Please refer to `HABQT-simulation --help` for executable usage instructions,
+accompanying Haddock documentation for Haskell API, and [libHABQT header
+file](./libHABQT.h) for shared library C ABI description.
+
+### Installation instructions
+
+#### Linux
+
+No special setup should be necessary. Simply make sure GSL, BLAS and LAPACK are
+installed on your system and install normally using stack, cabal or Setup.hs.
+
+#### Windows
+
+Making the necessary shared libraries and tools available on windows can be a
+bit tricky. One way to do this is install them inside MSYS2 that comes with
+stack on windows (on x86\_64 can be found under
+`%LOCALAPPDATA%\Programs\stack\x86_64-windows`):
+
+0. MSYS may not be present, in which case execute `stack build` in HABQT
+   directory. The build will fail due to missing libraries/tools, but MSYS
+   should be installed after it, and you will be able to add them.
+1. Launch msys2.exe
+2. (Optional) Update MSYS2 with `pacman -Syu`. It may be necessary to restart
+   the shell: follow the instructions displayed in it.
+3. Install the appropriate mingw toolchain, which includes necessary tools like
+   pkg-config. E.g.  
+   ```pacman -S mingw-w64-x86_64-toolchain```
+4. Install GSL for the appropriate mingw toolchain (the 64-bit one this case):  
+   ```pacman -S mingw64/mingw-w64-x86_64-gsl```
+5. Install openblas for appropriate mingw toolchain:  
+   ```pacman -S mingw64/mingw-w64-x86_64-openblas```
+6. The versions/naming conventions hmatrix expects differ from what is used in
+   modern MSYS2, so it's necessary to either link of create renamed copies of
+   two libraries from the appropriate mingw toolchain (in my case found under
+   `%LOCALAPPDATA%\Programs\stack\x86_64-windows\msys2-20150512\mingw64\bin`):
+   `libgfortran-4.dll` to `libgfortran.dll` and `libgslcblas-0.dll` to
+   `gsl-0.dll`. I recommend placing links/copies in some directory that isn't
+   normally on PATH (neither windows nor MSYS) and explicitly pointing stack to
+   them during installation with `--extra-lib-dirs ` (example follows in next
+   step).
+7. Outside MSYS2 open a normal windows shell, navigate to HABQT folder and
+   build/install with stack, passing appropriate flags and library dirs:  
+   ```stack build --flag hmatrix:openblas --extra-lib-dirs=D:\lib```  
+   where `D:\lib` contains ` libgfortran.dll` and `gsl-0.dll`.
+
+To use the shared library or executable, you’ll need to have several mingw
+libraries on PATH or in the same directory: libgcc_s_seh-1.dll,
+libgfortran-4.dll, libgsl-23.dll, libgslcblas-0.dll, libopenblas.dll,
+libquadmath-0.dll and libwinpthread-1.dll.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,3 @@
+import Distribution.Simple
+
+main = defaultMain
diff --git a/app/Main.hs b/app/Main.hs
new file mode 100644
--- /dev/null
+++ b/app/Main.hs
@@ -0,0 +1,181 @@
+{-# LANGUAGE RecordWildCards #-}
+
+module Main where
+
+import Data.Semigroup ((<>))
+import HABQTlib.Data
+import HABQTlib.UnsafeAPI
+import Options.Applicative
+import Streaming (Of, Stream)
+import qualified Streaming.Prelude as S
+import qualified System.IO as IO
+
+data CLIargs = CLIargs
+  { cliQbNum :: QBitNum
+  , cliRank :: Rank
+  , cliMeasNum :: Int
+  , cliExpNum :: Int
+  , cliPtNum :: NumberOfParticles
+  , cliVerb :: Int
+  , cliFilePath :: IO.FilePath
+  , cliMHMCiter :: MHMCiter
+  , cliOptIter :: OptIter
+  }
+
+msgPositive :: String
+msgPositive = "must be a positive integer"
+
+readInt :: (Int -> Bool) -> String -> String -> Either String Int
+readInt cond msg s =
+  let pstr :: [(Int, String)]
+      pstr = reads s
+      go [(i, "")] =
+        if cond i
+          then Right i
+          else Left msg
+      go _ = Left msg
+   in go pstr
+
+readPostiveInt :: String -> Either String Int
+readPostiveInt = readInt (> 0) msgPositive
+
+positive :: ReadM Int
+positive = eitherReader readPostiveInt
+
+cliparse :: Parser CLIargs
+cliparse =
+  CLIargs <$>
+  option
+    positive
+    (long "quantum-bit-number" <> short 'q' <>
+     help "Number of quantum bits under tomography" <>
+     showDefault <>
+     value 1 <>
+     metavar "QBNUM") <*>
+  option
+    positive
+    (long "rank" <> short 'r' <> help "Rank of true states" <> showDefault <>
+     value 1 <>
+     metavar "RANK") <*>
+  option
+    positive
+    (long "measurements-per-experiment" <> short 'm' <>
+     help "Number of measurements to make per-experiment" <>
+     showDefault <>
+     value 100 <>
+     metavar "MNUM") <*>
+  option
+    positive
+    (long "experiment-number" <> short 'e' <>
+     help "Number of tomography experiments to run" <>
+     showDefault <>
+     value 10 <>
+     metavar "EXPNUM") <*>
+  option
+    positive
+    (long "particle-number" <> short 'p' <>
+     help
+       "Number of particles to use for approximating the distribution over states (per rank)" <>
+     showDefault <>
+     value 1000 <>
+     metavar "PNUM") <*>
+  option
+    auto
+    (long "verbosity" <> short 'v' <>
+     help
+       "Verbosity level of output to stdout from 0 (no output) to 2 (full output)" <>
+     showDefault <>
+     value 2 <>
+     metavar "VERB") <*>
+  strOption
+    (long "output-file-path" <> short 'o' <>
+     help
+       "Path to file which infidelities of tomographic estimates will be appended to" <>
+     showDefault <>
+     value "./output.txt" <>
+     metavar "OUTPATH") <*>
+  option
+    positive
+    (long "mhmc-iterations" <>
+     help
+       "Number of Metropolis–Hastings steps to perform when resampling (after adjusting for acceptance rate)" <>
+     showDefault <>
+     value 50 <>
+     metavar "MHMCITER") <*>
+  option
+    positive
+    (long "opt-iterations" <>
+     help
+       "Number of optimisation steps to perform when searching for optimal measurment" <>
+     showDefault <>
+     value 50 <>
+     metavar "POVMITER")
+
+writeCsInfid ::
+     Int
+  -> Int
+  -> QBitNum
+  -> Rank
+  -> NumberOfParticles
+  -> MHMCiter
+  -> OptIter
+  -> OutputVerb
+  -> IO.FilePath
+  -> IO ()
+writeCsInfid numRuns numMeas qbn rank pn mi oi v fp = do
+  let tomStr = streamTo numMeas (streamResults' qbn rank pn mi oi v)
+      multStr = S.replicateM numRuns . tomStr
+      dupIO fh s = do
+        IO.hPutStr fh s
+        IO.hFlush IO.stdout
+  IO.withFile fp IO.AppendMode (S.effects . S.effects . multStr . dupIO)
+
+streamTo ::
+     Int
+  -> Stream (Of Double) IO ()
+  -> (String -> IO ())
+  -> Stream (Of ()) IO ()
+streamTo n is sf = do
+  let rs = S.map show is
+      fs = S.take 1 rs
+      rest = S.drop 1 rs
+  S.mapM sf . S.yield $ "\n"
+  S.mapM sf fs
+  S.mapM (\x -> sf (", " ++ x)) . S.take n $ rest
+  S.mapM sf . S.yield $ "\n"
+
+main :: IO ()
+main = do
+  CLIargs {..} <- execParser opts
+  r' <-
+    if cliRank > 2 ^ cliQbNum
+      then do
+        putStrLn "WARNING: RANK > 2 ^ QBNUM, truncating to 2 ^ QBNUM"
+        return (2 ^ cliQbNum)
+      else return cliRank
+  v <-
+    case cliVerb of
+      0 -> return NoOutput
+      1 -> return FidOutput
+      2 -> return FullOutput
+      _ -> do
+        putStrLn
+          "WARNING: VERB isn't equal to 0,1 or 2, defaulting to full output"
+        return FullOutput
+  writeCsInfid
+    cliExpNum
+    cliMeasNum
+    cliQbNum
+    r'
+    cliPtNum
+    cliMHMCiter
+    cliOptIter
+    v
+    cliFilePath
+  where
+    opts =
+      info
+        (cliparse <**> helper)
+        (fullDesc <>
+         progDesc "Simulate HABQT and write infidelities of estimates to file" <>
+         header "Hierarchical Adaptive Bayesian Quantum Tomography simulation")
diff --git a/libHABQT/ForeignHABQT.hs b/libHABQT/ForeignHABQT.hs
new file mode 100644
--- /dev/null
+++ b/libHABQT/ForeignHABQT.hs
@@ -0,0 +1,54 @@
+module ForeignHABQT where
+
+import Data.Bifunctor (bimap)
+import Data.Complex (Complex(..))
+import Data.List (unzip)
+import qualified Data.Vector as V
+import Foreign
+import Foreign.C
+import HABQTlib.Data
+import HABQTlib.MeasurementProcessing (PurePOVM)
+import qualified Numeric.LinearAlgebra as LA
+
+tdm :: DensityMatrix
+tdm = DensityMatrix $ LA.real $ (2 LA.>< 2) [1 ..]
+
+type CMatrix = [[CDouble]]
+
+convertDM :: DensityMatrix -> (CMatrix, CMatrix)
+convertDM (DensityMatrix dm) =
+  let (rp, ip) = LA.fromComplex dm
+      conv = (map . map) realToFrac . LA.toLists
+   in (conv rp, conv ip)
+
+unmarshallSV :: Dim -> Ptr CDouble -> Ptr CDouble -> IO PureStateVector
+unmarshallSV dim rPtr iPtr = do
+  rl <- peekArray dim rPtr
+  il <- peekArray dim iPtr
+  let toHM = (dim LA.>< 1) . map realToFrac
+  return . PureStateVector $ LA.toComplex (toHM rl, toHM il)
+
+convertWPSV :: WeighedPureStateVector -> ([CDouble], [CDouble])
+convertWPSV (WeighedPureStateVector (w, PureStateVector sv)) =
+  let sv' = LA.scale (sqrt w :+ 0) sv
+      (rp, ip) = LA.fromComplex sv'
+      conv = map realToFrac . LA.toList . head . LA.toColumns
+   in (conv rp, conv ip)
+
+convertPOVM :: PurePOVM -> ([CDouble], [CDouble])
+convertPOVM povm =
+  let svs = V.toList povm
+      csvs = map convertWPSV svs
+   in bimap concat concat . unzip $ csvs
+
+marshallPOVM :: PurePOVM -> Ptr CDouble -> Ptr CDouble -> IO ()
+marshallPOVM povm rPtr iPtr = do
+  let (r, i) = convertPOVM povm
+  pokeArray rPtr r
+  pokeArray iPtr i
+
+marshallDM :: DensityMatrix -> Ptr CDouble -> Ptr CDouble -> IO ()
+marshallDM dm rPtr iPtr = do
+  let (r, i) = convertDM dm
+  pokeArray rPtr $ concat r
+  pokeArray iPtr $ concat i
diff --git a/libHABQT/LibHABQT.hs b/libHABQT/LibHABQT.hs
new file mode 100644
--- /dev/null
+++ b/libHABQT/LibHABQT.hs
@@ -0,0 +1,108 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+module LibHABQT where
+
+import Control.Monad.State.Lazy
+import Data.IORef
+import Data.Validation
+import Foreign
+import Foreign.C
+import ForeignHABQT
+import HABQTlib.Data
+import HABQTlib.Data.Particle
+import HABQTlib.UnsafeAPI
+import qualified System.IO as IO
+import qualified System.Random.MWC as MWC
+
+type TomForegin
+   = Ptr CDouble -- measurement SV real
+      -> Ptr CDouble -- measurement SV imag
+          -> Ptr CDouble -- estimate dm real
+              -> Ptr CDouble -- estimate dm imag
+                  -> Ptr CDouble -- next POVM real
+                      -> Ptr CDouble -- next POVM imag
+                          -> IO ()
+
+type ArgStorage = (QBitNum, OutputVerb, MHMCiter, OptIter, MWC.GenIO)
+
+type Storage = IORef (ArgStorage, (ParticleHierarchy, [PureStateVector]))
+
+type InitFun = CInt -> CInt -> CInt -> CInt -> CInt -> IO (StablePtr Storage)
+
+foreign export ccall "tomography_init" tomInit :: InitFun
+
+tomInit :: InitFun
+tomInit qn' pc' mi' oi' verb' = do
+  let inp = validateInputs qn' pc' mi' oi' verb'
+  gen <- MWC.createSystemRandom
+  case inp of
+    Failure errs -> do
+      IO.hSetBuffering IO.stderr IO.LineBuffering
+      IO.hPutStrLn IO.stderr (unlines errs)
+      IO.hPutStrLn
+        IO.stderr
+        "Free the resulting pointer using tomography_free and re-initialise with valid inputs"
+      IO.hFlush IO.stderr
+      ph <- initialiseParticleHierarchy 2 1
+      mem <- newIORef ((1, NoOutput, 1, 1, gen), (ph, []))
+      newStablePtr mem
+    Success (qn, pc, mi, oi, verb) -> do
+      let dim = 2 ^ qn
+          out =
+            case verb of
+              0 -> NoOutput
+              1 -> FullOutput
+      ph <- initialiseParticleHierarchy dim pc
+      mem <- newIORef ((qn, out, mi, oi, gen), (ph, []))
+      newStablePtr mem
+
+foreign export ccall "tomography" foreignTomFun
+  :: StablePtr Storage -> TomForegin
+
+foreignTomFun :: StablePtr Storage -> TomForegin
+foreignTomFun strPtr svrPtr sviPtr dmrPtr dmiPtr povmrPtr povmiPtr = do
+  mem <- deRefStablePtr strPtr
+  ((qn, out, mi, oi, gen), s) <- readIORef mem
+  let dim = 2 ^ qn
+  sv <- unmarshallSV dim svrPtr sviPtr
+  case validSV sv of
+    Success _ -> do
+      ((dm, nextPOVM), sn) <- runStateT (tomographyFun' qn mi oi out gen sv) s
+      writeIORef mem ((qn, out, mi, oi, gen), sn)
+      marshallDM dm dmrPtr dmiPtr
+      marshallPOVM nextPOVM povmrPtr povmiPtr
+    Failure [msg] -> do
+      IO.hSetBuffering IO.stderr IO.LineBuffering
+      IO.hPutStrLn IO.stderr $ msg ++ " Doing nothing."
+      IO.hFlush IO.stderr
+
+foreign export ccall "tomography_free" tomFree :: StablePtr Storage -> IO ()
+
+tomFree :: StablePtr Storage -> IO ()
+tomFree = freeStablePtr
+
+validateInputs ::
+     CInt
+  -> CInt
+  -> CInt
+  -> CInt
+  -> CInt
+  -> Validation [String] (QBitNum, NumberOfParticles, MHMCiter, OptIter, Int)
+validateInputs qn' pn' mi' oi' verb' =
+  let qn = fromIntegral qn'
+      pn = fromIntegral pn'
+      verb = fromIntegral verb'
+      mi = fromIntegral mi'
+      oi = fromIntegral oi'
+      qnM = ["Number of quantum bits must be a positive integer."]
+      pnM = ["Number of particles per rank must be a positive integer."]
+      vM =
+        [ "Verbosity level must be an integer equal to 0 (no output) or 1 (full output)."
+        ]
+      miM = ["Number of MHMC iterations must be a positive integer."]
+      vQn = validate qnM (> 0) qn
+      vPn = validate pnM (> 0) pn
+      vVerb = validate vM ((||) <$> (== 0) <*> (== 1)) verb
+      vMi = validate miM (> 0) mi
+      vOi = validOptIter oi
+   in (,,,,) <$> vQn <*> vPn <*> vMi <*> vOi <*> vVerb
diff --git a/libHABQT/hsinit.c b/libHABQT/hsinit.c
new file mode 100644
--- /dev/null
+++ b/libHABQT/hsinit.c
@@ -0,0 +1,17 @@
+#include <stdlib.h>
+#include "HsFFI.h"
+
+static void HsStart(void) __attribute__((constructor));
+static void HsStart(void)
+{
+  static char *argv[] = { "+RTS", "-A32m", NULL }, **argv_ = argv;
+  static int argc = 2;
+  hs_init(NULL,NULL);
+}
+
+
+static void HsEnd(void) __attribute__((destructor));
+extern void HsEnd(void)
+{
+  hs_exit();
+}
diff --git a/src/HABQTlib.hs b/src/HABQTlib.hs
new file mode 100644
--- /dev/null
+++ b/src/HABQTlib.hs
@@ -0,0 +1,80 @@
+{-|
+Module      : HABQTlib
+
+This module contains functions for performing and simulating HABQT in Haskell.
+
+Note: functions in this module simply call API from "HABQTlib.UnsafeAPI" after validating inputs.
+-}
+module HABQTlib
+  ( TomState
+  , TomFun
+  , tomographyFun
+  , simulatedTomography
+  , streamResults
+  ) where
+
+import Control.Monad.State.Lazy
+import Data.Validation
+import HABQTlib.Data
+import HABQTlib.MeasurementProcessing
+import HABQTlib.UnsafeAPI
+import Streaming (Of, Stream)
+import qualified System.Random.MWC as MWC
+
+validQMO ::
+     QBitNum
+  -> MHMCiter
+  -> OptIter
+  -> ( Validation [String] QBitNum
+     , Validation [String] MHMCiter
+     , Validation [String] OptIter)
+validQMO nq' mi' oi' =
+  let nq = validQBitN nq'
+      mi = validMHMCiter mi'
+      oi = validOptIter oi'
+   in (nq, mi, oi)
+
+-- | Given parameters such as output verbosity level and number of quantum
+-- bits, set up the tomography function.
+tomographyFun ::
+     QBitNum -- ^ Number of quantum bits under tomography
+  -> MHMCiter -- ^ Number of MHMC iterations to perform when resampling
+  -> OptIter -- ^ Number of POVM optimisation steps to perform
+  -> OutputVerb -- ^ Verbosity of stdout output
+  -> MWC.GenIO -- ^ IO generator for variates from "System.Random.MWC"
+  -> Validation [String] TomFun
+tomographyFun nq' mi' oi' outv gen =
+  let (nq, mi, oi) = validQMO nq' mi' oi'
+   in tomographyFun' <$> nq <*> mi <*> oi <*> Success outv <*> Success gen
+
+-- | Given a true state's density matrix and parameters, set up a simulation of
+-- quantum tomography that outputs infidelity between mean estimates and true
+-- state.
+simulatedTomography ::
+     DensityMatrix -- ^ True state's density matrix
+  -> QBitNum -- ^ Number of quantum bits under tomography
+  -> MHMCiter -- ^ Number of MHMC iterations to perform when resampling
+  -> OptIter -- ^ Number of POVM optimisation steps to perform
+  -> OutputVerb -- ^ Verbosity of stdout output
+  -> MWC.GenIO -- ^ IO generator for variates from "System.Random.MWC"
+  -> Validation [String] (StateT PurePOVM TomState Double)
+simulatedTomography trueDM nq' mi' oi' outv gen =
+  let dm = validDM trueDM
+      (nq, mi, oi) = validQMO nq' mi' oi'
+   in simulatedTomography' <$> dm <*> nq <*> mi <*> oi <*> Success outv <*>
+      Success gen
+
+-- | Stream simulated tomography results.
+streamResults ::
+     QBitNum -- ^ Number of quantum bits under tomography
+  -> Rank -- ^ Rank of true state
+  -> NumberOfParticles -- ^ Number of particles (per rank) to use for tomography
+  -> MHMCiter -- ^ Number of MHMC iterations to perform when resampling
+  -> OptIter -- ^ Number of POVM optimisation steps to perform
+  -> OutputVerb -- ^ Verbosity of stdout output
+  -> Validation [String] (Stream (Of Double) IO ())
+streamResults nq' rank' pn' mi' oi' outv =
+  let (nq, mi, oi) = validQMO nq' mi' oi'
+      rank = validRank rank'
+      pn = validPartNum pn'
+   in streamResults' <$> nq <*> rank <*> pn <*> mi <*> oi <*> Success outv
diff --git a/src/HABQTlib/Data.hs b/src/HABQTlib/Data.hs
new file mode 100644
--- /dev/null
+++ b/src/HABQTlib/Data.hs
@@ -0,0 +1,226 @@
+{-# LANGUAGE DeriveGeneric #-}
+
+{-|
+Module      : HABQTlib.Data
+
+This module contains data types and helper functions for working with quantum
+state vectors and density matrices.
+-}
+module HABQTlib.Data
+  ( Dim
+  , Rank
+  , NumberOfParticles
+  , QBitNum
+  , Weight
+  , MHMCiter
+  , OptIter
+  , OutputVerb(..)
+  , DensityMatrix(..)
+  , truncateRank
+  , PureStateVector(..)
+  , pureStateLikelihood
+  , svToDM
+  , WeighedDensityMatrix(..)
+  , mkWDM
+  , mkWDM1
+  , WeighedPureStateVector(..)
+  , (<+>)
+  , fidelity
+  , fidelityDM
+  , validSV
+  , validDM
+  , validPartNum
+  , validRank
+  , validMHMCiter
+  , validOptIter
+  , validQBitN
+  ) where
+
+import Control.Newtype.Generics (Newtype, unpack)
+import Data.Bool.HT (select)
+import Data.Complex
+import Data.Validation
+import GHC.Generics (Generic)
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra (Matrix)
+
+-- | Dimension of Hilbert space.
+type Dim = Int
+
+-- | Rank of mixed state.
+type Rank = Int
+
+-- | Number of particles per rank.
+type NumberOfParticles = Int
+
+-- | Number of quantum bits.
+type QBitNum = Int
+
+-- | Number of MHMC iterations to perform when resampling.
+type MHMCiter = Int
+
+-- | Number of optimisation steps to perform when searching for optimal
+-- measurment.
+type OptIter = Int
+
+-- | Weight associated with a particle.
+type Weight = Double
+
+-- | Output verbosity settings.
+data OutputVerb
+  = NoOutput -- ^ No stdout output
+  | FidOutput -- ^ Only output fidelities and weights of hierarchical mean estimates
+  | FullOutput -- ^ Full output, including resampling diagnostic information
+  deriving (Eq, Show, Ord)
+
+-- | Density matrix are stored as hmatrix matrices of complex doubles.
+newtype DensityMatrix = DensityMatrix
+  { getDensityMatrix :: Matrix (Complex Double)
+  } deriving (Eq, Show, Generic)
+
+-- | Pure state vectors are stored as hmatrix matrices of complex doubles.
+-- Such matrices only have one column.
+newtype PureStateVector = PureStateVector
+  { getStateVector :: Matrix (Complex Double)
+  } deriving (Eq, Show, Generic)
+
+instance Newtype DensityMatrix
+
+instance Newtype PureStateVector
+
+-- | Check whether a pure state vector is properly normed.
+validSV :: PureStateVector -> Validation [String] PureStateVector
+validSV =
+  validate
+    ["State vector must have unit norm."]
+    (\x -> abs (1 - LA.norm_2 (unpack x)) < 1e-12)
+
+-- | Verify that density matrix is Hermitian and has trace 1.
+validDM :: DensityMatrix -> Validation [String] DensityMatrix
+validDM =
+  let traceU :: LA.Matrix (Complex Double) -> Bool
+      traceU dm =
+        (abs (1 - (magnitude . LA.sumElements . LA.takeDiag) dm) < 1e-12)
+      hermU dm = (LA.norm_2 (dm - LA.tr dm) < 1e-6)
+      both = ((&&) <$> traceU <*> hermU) . unpack
+      dmM = ["Density matrix must be Hermitian and have trace of 1."]
+   in validate dmM both
+
+qnM :: [String]
+qnM = ["Number of quantum bits must be a positive integer."]
+
+-- | Verify that number of quantum bits is positive.
+validQBitN :: QBitNum -> Validation [String] QBitNum
+validQBitN = validate qnM (> 0)
+
+miM :: [String]
+miM = ["Number of MHMC iterations must be a positive integer."]
+
+-- | Verify that number of MHMC iterations is a positive integer.
+validMHMCiter :: MHMCiter -> Validation [String] MHMCiter
+validMHMCiter = validate miM (> 0)
+
+pnM :: [String]
+pnM = ["Number of particles per rank must be a positive integer."]
+
+-- | Verify that particle number is a positive integer.
+validPartNum :: NumberOfParticles -> Validation [String] NumberOfParticles
+validPartNum = validate pnM (> 0)
+
+rM :: [String]
+rM = ["Rank must be a positive integer."]
+
+-- | Verify that rank is a positive integer. Setting rank to be higher than
+-- the dimension of space creates poinless performance overhead, but isn't
+-- prevented by validation.
+validRank :: Rank -> Validation [String] Rank
+validRank = validate rM (> 0)
+
+oiM :: [String]
+oiM = ["Number of POVM optimisation iterations must be a positive integer."]
+
+-- | Verify that number of POVM optimisation iterations is positive.
+validOptIter :: OptIter -> Validation [String] OptIter
+validOptIter = validate oiM (> 0)
+
+-- | Weighed density matrix where weight is stored separately as first
+-- coordinate of a tuple.
+newtype WeighedDensityMatrix = WeighedDensityMatrix
+  { getWDM :: (Weight, DensityMatrix)
+  } deriving (Eq, Show, Generic)
+
+-- | A shorter alias for curried 'WeighedDensityMatrix' constructor.
+mkWDM :: Weight -> DensityMatrix -> WeighedDensityMatrix
+mkWDM w dm = WeighedDensityMatrix (w, dm)
+
+-- | Alias for @'mkWDM' 1@.
+mkWDM1 :: DensityMatrix -> WeighedDensityMatrix
+mkWDM1 = mkWDM 1
+
+-- | Weighed state vector where weight is stored separately as first coordinate
+-- of a tuple.
+newtype WeighedPureStateVector = WeighedPureStateVector
+  { getWSV :: (Weight, PureStateVector)
+  } deriving (Eq, Show, Generic)
+
+instance Newtype WeighedDensityMatrix
+
+instance Newtype WeighedPureStateVector
+
+-- | Fidelity (probability of measurement) between pure states.
+fidelity :: PureStateVector -> PureStateVector -> Double
+fidelity (PureStateVector sv1) (PureStateVector sv2) =
+  let ips = magnitude (LA.atIndex (LA.tr sv1 LA.<> sv2) (0, 0)) ^ (2 :: Int)
+   in select ips [(ips < 0, 0), (ips > 1, 1)]
+
+-- | Fidelity (probability of measurement) between mixed states.
+fidelityDM :: DensityMatrix -> DensityMatrix -> Double
+fidelityDM (DensityMatrix dm1) (DensityMatrix dm2) =
+  let (u1, s1) = LA.leftSV dm1
+      (u2, s2) = LA.leftSV dm2
+      ss1 = LA.cmap (\x -> sqrt x :+ 0) s1
+      ss2 = LA.cmap (\x -> sqrt x :+ 0) s2
+   in LA.norm_nuclear
+        (u1 LA.<> LA.diag ss1 LA.<> LA.tr u1 LA.<> u2 LA.<> LA.diag ss2 LA.<>
+         LA.tr u2) ^
+      (2 :: Int)
+
+-- | Calculate the density matrix of a given pure state.
+svToDM :: PureStateVector -> DensityMatrix
+svToDM (PureStateVector sv) = DensityMatrix $ sv LA.<> LA.tr sv
+
+infix 8 <+>
+
+-- | Given two weighed density matrixes, compute their mixture. Associative
+-- operation.
+(<+>) :: WeighedDensityMatrix -> WeighedDensityMatrix -> WeighedDensityMatrix
+wdm0 <+> wdm1 =
+  WeighedDensityMatrix
+    (w0 + w1, DensityMatrix $ LA.scale c0 dm0 + LA.scale c1 dm1)
+  where
+    up = fmap unpack . unpack
+    (w0, dm0) = up wdm0
+    (w1, dm1) = up wdm1
+    cs = LA.fromList [w0 :+ 0, w1 :+ 0]
+    csn = LA.scale (1 / LA.sumElements cs) cs
+    c0 = csn LA.! 0
+    c1 = csn LA.! 1
+
+-- | Probability of obtaining a measurement result when projecting a system in
+-- mixed state determined by a density matrix onto a pure state.
+pureStateLikelihood :: PureStateVector -> DensityMatrix -> Double
+pureStateLikelihood (PureStateVector sv) (DensityMatrix dm) =
+  let p =
+        LA.magnitude . LA.sumElements . LA.takeDiag $ LA.tr sv LA.<> dm LA.<> sv
+   in select p [(p < 0, 0), (p > 1, 1)]
+
+-- | Set smallest eigenvalues of a weighed density matrix to zero until
+-- specified rank is reached.
+truncateRank :: Rank -> WeighedDensityMatrix -> WeighedDensityMatrix
+truncateRank targetRank (WeighedDensityMatrix (w, DensityMatrix dm)) =
+  let (u, s, _) = LA.svd dm
+      st = LA.real $ LA.subVector 0 targetRank s
+      ut = u LA.¿ [0 .. (targetRank - 1)]
+      stn = LA.scale (1 / LA.sumElements st) st
+   in WeighedDensityMatrix
+        (w, DensityMatrix $ ut LA.<> LA.diag stn LA.<> LA.tr ut)
diff --git a/src/HABQTlib/Data/Particle.hs b/src/HABQTlib/Data/Particle.hs
new file mode 100644
--- /dev/null
+++ b/src/HABQTlib/Data/Particle.hs
@@ -0,0 +1,252 @@
+{-# LANGUAGE RecordWildCards #-}
+
+{-|
+Module      :  HABQTlib.Data.Particle
+
+Data structures and functions that deal with storing and processing particle
+hierarchies.
+
+/Warning/: functions in this module assume that the 'ptsParticles' is non-empty
+and 'NumberOfParticles', 'Dim', and 'Rank' are positive, no validation is
+performed. If you use them directly, instead of employing API from
+"HABQTlib", you must ensure those assumptions hold.
+-}
+module HABQTlib.Data.Particle
+  ( Particles(..)
+  , genParticles
+  , updateParticles
+  , ParticleHierarchy
+  , initialiseParticleHierarchy
+  , updateParticleHierarchy
+  , getMixedEstimate
+  , foldOverPts
+  , reduceParticlesToMean
+  , effectiveSize
+  , ResampleArgs(..)
+  , resampleMultinom
+  , resample
+  , ecdf
+  , icdf
+  , nudgeParticle
+  ) where
+
+import Control.Monad (when)
+import Control.Newtype.Generics (over)
+import Data.Bool.HT (if', select)
+import qualified Data.Vector as V
+import HABQTlib.Data
+import HABQTlib.RandomStates
+import Numeric.LinearAlgebra (Complex(..))
+import qualified Numeric.LinearAlgebra as LA
+import qualified System.Random.MWC as MWC
+import Text.Printf (printf)
+
+-- | A vector of weighed density matrices is stored along with their rank and
+-- number. 'ptsWeight' corresponds to the collective weight of particles of
+-- rank 'ptsRank' in the hierarchical model, it is not the sum of individual
+-- weights of particles (that is normalised to unity after every update).
+data Particles = Particles
+  { ptsRank :: Rank
+  , ptsWeight :: Weight
+  , ptsNumber :: NumberOfParticles
+  , ptsParticles :: V.Vector WeighedDensityMatrix
+  } deriving (Show)
+
+-- | Particle hierarchy is described by a vector of 'Particles'.
+type ParticleHierarchy = V.Vector Particles
+
+-- | Generates random particles (according to induced measure).
+genParticles :: Dim -> Rank -> NumberOfParticles -> IO Particles
+genParticles d r n =
+  let w = 1 / fromIntegral n
+   in Particles r 1 n . fmap (mkWDM w) <$> V.replicateM n (genDM d r)
+
+-- | Summarise particles to a mean estimate (weighed by the corresponding
+-- hierarchical weight of the rank).
+reduceParticlesToMean :: Particles -> WeighedDensityMatrix
+reduceParticlesToMean Particles {..} =
+  let wdm = V.foldl1' (<+>) ptsParticles
+      wdmw = over WeighedDensityMatrix (\(w, dm) -> (w * ptsWeight, dm)) wdm
+   in truncateRank ptsRank wdmw
+
+-- | Map density matrices, combine them with their weights, and then perform a
+-- (strict left) fold of results.
+foldOverPts ::
+     (DensityMatrix -> a) -- ^ function to map over density matrices
+  -> (Weight -> a -> b) -- ^ function to combine weights with results of mapping
+  -> (c -> b -> c) -- ^ fold funciton
+  -> c -- ^ seed value for folding
+  -> Particles
+  -> c
+foldOverPts f wf fld z Particles {..} =
+  let wm (WeighedDensityMatrix (w, dm)) = wf w (f dm)
+   in V.foldl' (\l r -> fld l (wm r)) z ptsParticles
+
+fullDataLogLikelihood :: [PureStateVector] -> DensityMatrix -> Double
+fullDataLogLikelihood vs dm =
+  let lps = map (log . (`pureStateLikelihood` dm)) vs
+   in sum lps
+
+-- | Given a measurement result, perform a Bayesian update over the particles.
+updateParticles :: PureStateVector -> Particles -> Particles
+updateParticles sv pts@Particles {..} =
+  let updateF :: WeighedDensityMatrix -> WeighedDensityMatrix
+      updateF (WeighedDensityMatrix (w, dm)) = WeighedDensityMatrix (wnew, dm)
+        where
+          wnew = w * pureStateLikelihood sv dm
+      upts = V.map updateF ptsParticles
+      uw = V.foldl' (\acc (WeighedDensityMatrix (w, _)) -> acc + w) 0 upts
+      npts = V.map (over WeighedDensityMatrix (\(w, dm) -> (w / uw, dm))) upts
+   in pts {ptsWeight = ptsWeight * uw, ptsParticles = npts}
+
+-- | Helper function that generates random particles of each applicable rank.
+initialiseParticleHierarchy :: Dim -> NumberOfParticles -> IO ParticleHierarchy
+initialiseParticleHierarchy d n = V.generateM d (\r -> genParticles d (r + 1) n)
+
+-- | Given a measurement result, update all particles, then normalise resulting
+-- hierarchical weights to sum to unity.
+updateParticleHierarchy ::
+     PureStateVector -> ParticleHierarchy -> ParticleHierarchy
+updateParticleHierarchy sv ph =
+  let uph = V.map (updateParticles sv) ph
+      wgts = V.map ptsWeight uph
+      nwgts = V.map (/ V.sum wgts) wgts
+   in V.zipWith (\x w -> x {ptsWeight = w}) uph nwgts
+
+-- | Summarise the whole particle hierarchy into one mean Bayesian estimate.
+getMixedEstimate :: ParticleHierarchy -> DensityMatrix
+getMixedEstimate ph =
+  let rankEstimates = V.map reduceParticlesToMean ph
+      WeighedDensityMatrix (_, result) = V.foldl1' (<+>) rankEstimates
+   in result
+
+-- | Calculate the effective sample size of particles (weights don’t
+-- necessarily have to be normalised).
+effectiveSize :: Particles -> Double
+effectiveSize Particles {..} =
+  let ss = V.sum . V.map ((^ (2 :: Int)) . fst . getWDM) $ ptsParticles
+      wa = V.foldl' (flip ((+) . fst . getWDM)) 0 ptsParticles
+   in wa ^ (2 :: Int) / ss
+
+-- | Nudges a particle by mixing the state together with some randomly
+-- generated pure state. Relative weight of the random component determines how
+-- “close” a nudged particle is to the original one.
+nudgeParticle ::
+     Dim
+  -> Weight -- ^ Relative weight (from 0 to 1) of random component
+  -> WeighedDensityMatrix
+  -> IO WeighedDensityMatrix
+nudgeParticle dim weightFraction (WeighedDensityMatrix (w, dm)) = do
+  DensityMatrix nudgeDM <- svToDM <$> genPureSV dim
+  let dmw = LA.scale (1 - (weightFraction :+ 0)) (getDensityMatrix dm)
+      dmwn = LA.scale (weightFraction :+ 0) nudgeDM
+  return $ WeighedDensityMatrix (w, DensityMatrix $ dmw + dmwn)
+
+-- | Calculates values of empirical distribution function at data points.
+ecdf :: V.Vector WeighedDensityMatrix -> V.Vector Double
+ecdf = V.postscanl' (+) 0 . V.map (fst . getWDM)
+
+-- | /O(log n)/ Given a non-empty sorted vector (typically an empirical cdf
+-- evaluated at data points returned by ecdf) and a real number return the
+-- (0-based) index of the least element of vector which is greater or equal to
+-- the given real number (or the index of the last element, in case there is no
+-- element smaller than the argument).
+icdf :: V.Vector Double -> Double -> Int
+icdf cdf x =
+  let tIdx = V.length cdf - 1
+      go (lIdx, hIdx) =
+        let mIdx =
+              truncate $ ((fromIntegral lIdx :: Double) + fromIntegral hIdx) / 2
+         in select
+              (go (lIdx, mIdx))
+              [ (lIdx == hIdx, lIdx)
+              , (lIdx + 1 == hIdx, if' (cdf V.! lIdx > x) lIdx hIdx)
+              , (cdf V.! mIdx < x, go (mIdx, hIdx))
+              ]
+   in select (go (0, tIdx)) [(x <= V.head cdf, 0), (x > V.last cdf, tIdx)]
+
+-- | Multinomial resampling of particle vector, which equalises weights of
+-- particles.
+resampleMultinom :: MWC.GenIO -> Particles -> IO Particles
+resampleMultinom gen pts@Particles {..} = do
+  us <- MWC.uniformVector gen ptsNumber
+  let cdf = ecdf ptsParticles
+      idxs = V.map (icdf cdf) us
+      w = 1 / fromIntegral ptsNumber
+      pointR = over WeighedDensityMatrix (\(_, dm) -> (w, dm))
+      resampled = V.map (ptsParticles V.!) idxs
+      normed = V.map pointR resampled
+  return pts {ptsParticles = normed}
+
+mhmcStep ::
+     MWC.GenIO
+  -> Dim
+  -> Double
+  -> [PureStateVector]
+  -> Particles
+  -> IO (Double, Particles)
+mhmcStep gen dim rw ms pts@Particles {..} = do
+  let cr wdm wdm' =
+        exp
+          (fullDataLogLikelihood ms (snd . getWDM $ wdm') -
+           fullDataLogLikelihood ms (snd . getWDM $ wdm))
+  newParticles <-
+    V.mapM (fmap (truncateRank ptsRank) . nudgeParticle dim rw) ptsParticles
+  us <- V.replicateM ptsNumber (MWC.uniform gen :: IO Double)
+  let crs = V.zipWith cr ptsParticles newParticles
+      change = V.zipWith (<=) us crs
+      accRate =
+        (fromIntegral . V.length . V.filter id) change / fromIntegral ptsNumber
+      rwdms = V.zipWith3 if' change newParticles ptsParticles
+      final = pts {ptsParticles = rwdms}
+  return (accRate, final)
+
+resampleMHMC ::
+     ResampleArgs
+  -> DensityMatrix
+  -> Double
+  -> Int
+  -> [PureStateVector]
+  -> Particles
+  -> IO Particles
+resampleMHMC ra@ResampleArgs {..} estimate wr iter mts pts = do
+  (accRate, resampled) <- mhmcStep argGen argDim wr mts pts
+  when (argOut == FullOutput) $
+    printf
+      "(Weight of new particle: %8.3g, MHMC acceptance rate: %8.3g)\n"
+      wr
+      accRate
+  let (iter', wr') =
+        select
+          (iter + 1, wr)
+          [ (accRate < 1e-2, (0, wr * 0.25))
+          , (accRate < 1e-1, (0, wr * 0.5))
+          , (iter < argMinIter, (iter + 1, wr))
+          , (accRate < 0.33, (0, wr * 0.5))
+          ]
+  if iter > argMinIter
+    then return resampled
+    else resampleMHMC ra estimate wr' iter' mts resampled
+
+-- | Arguments for the resampling function.
+data ResampleArgs = ResampleArgs
+  { argOut :: OutputVerb
+  , argGen :: MWC.GenIO
+  , argDim :: Dim
+  , argMinIter :: MHMCiter
+  }
+
+-- | Resample particles. First does one multinomial step that equalises the
+-- weights, then performs MHMC iterations adaptively refining the proposal
+-- distribution based on acceptance rate. 'argMinIter' iterations are performed
+-- for proposal distributions with significant acceptance rates.
+resample :: ResampleArgs -> [PureStateVector] -> Particles -> IO Particles
+resample ra@ResampleArgs {..} mts pts@Particles {..} = do
+  let estimate = getMixedEstimate . V.singleton $ pts
+      nudgeW = 0.95
+  when (argOut == FullOutput) $ do
+    putStrLn ""
+    putStrLn $ "resampling rank " ++ show ptsRank
+    putStrLn ""
+  rm <- resampleMultinom argGen pts
+  resampleMHMC ra estimate nudgeW 0 mts rm
diff --git a/src/HABQTlib/MeasurementProcessing.hs b/src/HABQTlib/MeasurementProcessing.hs
new file mode 100644
--- /dev/null
+++ b/src/HABQTlib/MeasurementProcessing.hs
@@ -0,0 +1,151 @@
+{-|
+Module      : HABQTlib.MeasurementProcessing
+
+Functions that deal with optimising and simulating measurements that take place
+during tomography.
+-}
+module HABQTlib.MeasurementProcessing
+  ( PurePOVM
+  , SingleQbParam
+  , measurementProbs
+  , svToAngles
+  , blochAnglesToSV
+  , mkAntipodalPOVM
+  , productPOVM
+  , simulateMeasuremet
+  , optimiseSingleQbPOVM
+  ) where
+
+import Control.Applicative (liftA2)
+import Control.Newtype.Generics (unpack)
+import Data.Complex (mkPolar, polar)
+import Data.List (unfoldr)
+import Data.Maybe (fromJust)
+import qualified Data.Vector as V
+import HABQTlib.Data
+import HABQTlib.Data.Particle
+import qualified Numeric.GSL as GSL
+import Numeric.LinearAlgebra (Complex((:+)))
+import qualified Numeric.LinearAlgebra as LA
+import qualified System.Random.MWC as MWC
+import System.Random.MWC.Distributions (categorical)
+
+-- | A POVM consisting of projections onto pure states.
+type PurePOVM = V.Vector WeighedPureStateVector
+
+-- | Spherical coordinates of a pure single qubit state on Bloch sphere.
+type SingleQbParam = (Double, Double)
+
+-- | Probabilities to measure elements of POVM when performing the measurement
+-- over a mixed state.
+measurementProbs :: PurePOVM -> DensityMatrix -> V.Vector Double
+measurementProbs povm dm = V.map (weighedProb dm) povm
+
+weighedProb :: DensityMatrix -> WeighedPureStateVector -> Double
+weighedProb dm (WeighedPureStateVector (w, sv)) = w * pureStateLikelihood sv dm
+
+entropy :: V.Vector Double -> Double
+entropy = negate . V.sum . V.map (\p -> p * log p)
+
+povmPointEntropy :: PurePOVM -> DensityMatrix -> Double
+povmPointEntropy povm dm = entropy (measurementProbs povm dm)
+
+povmMeanEntropy :: PurePOVM -> ParticleHierarchy -> Double
+povmMeanEntropy povm ph =
+  let meanEnt = foldOverPts (povmPointEntropy povm) (*) (+) 0
+      totalWeight = V.sum . V.map ptsWeight $ ph
+      rawEntropy = V.sum . V.map (liftA2 (*) ptsWeight meanEnt) $ ph
+   in rawEntropy / totalWeight
+
+-- | Recover a state vector from spherical coordinates.
+blochAnglesToSV :: SingleQbParam -> PureStateVector
+blochAnglesToSV (th, phi) =
+  let z = LA.fromList [1, 0]
+      o = LA.fromList [0, 1]
+   in PureStateVector . LA.fromColumns . pure $
+      LA.scalar (cos (th / 2) :+ 0) * z +
+      LA.scalar (mkPolar (sin (th / 2)) phi) * o
+
+-- | Return spherical coordinates of a single qubit pure state (on Bloch
+-- sphere).
+svToAngles :: PureStateVector -> Maybe SingleQbParam
+svToAngles (PureStateVector sv) =
+  let s = LA.size sv
+      (m0, ph0) = polar $ LA.atIndex sv (0, 0)
+      (_, ph1) = polar $ LA.atIndex sv (1, 0)
+      ph = ph1 - ph0
+      th = 2 * acos m0
+   in if s == (2, 1)
+        then Just (th, ph)
+        else Nothing
+
+-- | Given spherical coordinates, construct a POVM from the given vector and
+-- one orthogonal to it.
+mkAntipodalPOVM :: SingleQbParam -> PurePOVM
+mkAntipodalPOVM c@(th, phi) =
+  let sv = blochAnglesToSV c
+      sv' = blochAnglesToSV (pi - th, phi + pi)
+   in V.fromList $ WeighedPureStateVector <$> [(1, sv), (1, sv')]
+
+reshapeList :: Int -> [a] -> [[a]]
+reshapeList n =
+  unfoldr
+    (\b ->
+       if length b < n
+         then Nothing
+         else Just (splitAt n b))
+
+listToPairs :: [a] -> [(a, a)]
+listToPairs = fmap (\(a:[b]) -> (a, b)) . reshapeList 2
+
+pairToList :: (a, a) -> [a]
+pairToList (a, b) = [a, b]
+
+-- | Given a list of POVM measurements on sub-systems, construct a POVM over
+-- the composite system that includes all of them.
+productPOVM :: [PurePOVM] -> PurePOVM
+productPOVM sqbPovms =
+  let pr ::
+           WeighedPureStateVector
+        -> WeighedPureStateVector
+        -> WeighedPureStateVector
+      pr wsv0 wsv1 =
+        WeighedPureStateVector (w0 * w1, PureStateVector $ LA.kronecker sv0 sv1)
+        where
+          up = fmap unpack . unpack
+          (w0, sv0) = up wsv0
+          (w1, sv1) = up wsv1
+   in V.foldl1' (liftA2 pr) $ V.fromList sqbPovms
+
+-- | Approximate most informative separable POVM over a composite system of
+-- quantum bits, given a list of single-qubit starting points and a particle
+-- distribution.
+optimiseSingleQbPOVM ::
+     OptIter -- ^ Number of optimisation steps to perform
+  -> [PureStateVector] -- ^ single qubit initial states
+  -> ParticleHierarchy
+  -> PurePOVM
+optimiseSingleQbPOVM iter sv0s ph =
+  let nq = length sv0s
+      method = GSL.NMSimplex2
+      precision = 1e-6
+      iterations = iter
+      initialBox = replicate (2 * nq) (pi / 2)
+      estimateDM = getMixedEstimate ph
+      obj params =
+        negate $ povmPointEntropy povm estimateDM - povmMeanEntropy povm ph
+        where
+          povm = productPOVM . fmap mkAntipodalPOVM . listToPairs $ params
+      start = concatMap (pairToList . fromJust . svToAngles) sv0s
+      result = GSL.minimize method precision iterations initialBox obj start
+   in productPOVM . fmap mkAntipodalPOVM . listToPairs . fst $ result
+
+-- | Simulate a POVM over a mixed state and return the state vector on which
+-- the projection was obtained.
+simulateMeasuremet ::
+     DensityMatrix -> PurePOVM -> MWC.GenIO -> IO PureStateVector
+simulateMeasuremet dm povm gen = do
+  let probs = measurementProbs povm dm
+      svs = V.map (\(WeighedPureStateVector (_, sv)) -> sv) povm
+  idx <- categorical probs gen
+  return $ svs V.! idx
diff --git a/src/HABQTlib/RandomStates.hs b/src/HABQTlib/RandomStates.hs
new file mode 100644
--- /dev/null
+++ b/src/HABQTlib/RandomStates.hs
@@ -0,0 +1,40 @@
+{-|
+Module      : HABQTlib.RandomStates
+
+Generation of random pure (from Haar measure) and mixed states (from measure
+induced by partial tracing of purified states).
+-}
+module HABQTlib.RandomStates
+  ( genPureSV
+  , genDM
+  ) where
+
+import HABQTlib.Data
+import Numeric.LinearAlgebra
+  ( norm_2
+  , randn
+  , scalar
+  , sumElements
+  , takeDiag
+  , toComplex
+  , tr
+  )
+import qualified Numeric.LinearAlgebra as LA
+
+-- | Generate a random mixed state of specified rank.
+genDM :: Dim -> Rank -> IO DensityMatrix
+genDM dim r = do
+  r1 <- randn dim r
+  r2 <- randn dim r
+  let a = toComplex (r1, r2)
+  let h = a LA.<> tr a
+  let hTr = sumElements $ takeDiag h
+  return . DensityMatrix $ h / scalar hTr
+
+-- | Generate a random pure state from Hilbert space of given dimension.
+genPureSV :: Dim -> IO PureStateVector
+genPureSV dim = do
+  r1 <- randn dim 1
+  r2 <- randn dim 1
+  let sv = toComplex (r1, r2)
+  return . PureStateVector $ sv / toComplex (scalar (norm_2 sv), 0)
diff --git a/src/HABQTlib/UnsafeAPI.hs b/src/HABQTlib/UnsafeAPI.hs
new file mode 100644
--- /dev/null
+++ b/src/HABQTlib/UnsafeAPI.hs
@@ -0,0 +1,116 @@
+{-|
+Module      : HABQTlib.UnsafeAPI
+
+This module contains functions for performing and simulating HABQT in Haskell.
+
+__Caution__: functions in this module perform no input validation and are partial. For a safe API refer to "HABQTlib".
+-}
+module HABQTlib.UnsafeAPI where
+
+import Control.Applicative (liftA2)
+import Control.Monad (replicateM)
+import Control.Monad.State.Lazy
+import Data.Maybe (fromJust)
+import qualified Data.Vector as V
+import HABQTlib.Data
+import HABQTlib.Data.Particle
+import HABQTlib.MeasurementProcessing
+import HABQTlib.RandomStates
+import qualified Numeric.LinearAlgebra as LA
+import Streaming
+import qualified Streaming.Prelude as S
+import qualified System.Random.MWC as MWC
+import Text.Printf (printf)
+
+-- | Tomography keeps track of the particle hierarchy and list of previous
+-- measurement results, IO is used for verbose output and assorted random state
+-- generation.
+type TomState = StateT (ParticleHierarchy, [PureStateVector]) IO
+
+-- | Tomography function takes a measurement result and returns state-dependent
+-- Bayesian mean estimate of state and the optimal next POVM to perform.
+type TomFun = PureStateVector -> TomState (DensityMatrix, PurePOVM)
+
+-- | Given parameters such as output verbosity level and number of quantum
+-- bits, set up the tomography function.
+tomographyFun' ::
+     QBitNum -- ^ Number of quantum bits under tomography
+  -> MHMCiter -- ^ Number of MHMC iterations to perform when resampling
+  -> OptIter -- ^ Number of POVM optimisation steps to perform
+  -> OutputVerb -- ^ Verbosity of stdout output
+  -> MWC.GenIO -- ^ IO generator for variates from "System.Random.MWC"
+  -> TomFun
+tomographyFun' nq mi oi outv gen nextResult = do
+  (ph, ms) <- get
+  let nextPH = updateParticleHierarchy nextResult ph
+      dim = LA.rows . getStateVector $ nextResult
+      effectiveSizes =
+        V.map (liftA2 (/) effectiveSize (fromIntegral . ptsNumber)) nextPH
+      ra = ResampleArgs outv gen dim mi
+      resampleC es pts =
+        if es < 0.5
+          then resample ra (nextResult : ms) pts
+          else return pts
+  nextPH' <- lift $ V.zipWithM resampleC effectiveSizes nextPH
+  sv0s <- liftIO $ replicateM nq (genPureSV 2)
+  let nextEstimate = getMixedEstimate nextPH'
+      nextPOVM = optimiseSingleQbPOVM oi sv0s nextPH'
+  put (nextPH', nextResult : ms)
+  return (nextEstimate, nextPOVM)
+
+-- | Given a true state's density matrix and parameters, set up a simulation of
+-- quantum tomography that outputs infidelity between mean estimates and true
+-- state.
+simulatedTomography' ::
+     DensityMatrix -- ^ True state's density matrix
+  -> QBitNum -- ^ Number of quantum bits under tomography
+  -> MHMCiter -- ^ Number of MHMC iterations to perform when resampling
+  -> OptIter -- ^ Number of POVM optimisation steps to perform
+  -> OutputVerb -- ^ Verbosity of stdout output
+  -> MWC.GenIO -- ^ IO generator for variates from "System.Random.MWC"
+  -> StateT PurePOVM TomState Double
+simulatedTomography' trueDM nq mi oi outv gen = do
+  povm <- get
+  nextResult <- liftIO $ simulateMeasuremet trueDM povm gen
+  (nextEstimate, nextPOVM) <- lift $ tomographyFun' nq mi oi outv gen nextResult
+  (nextPH, _) <- lift get
+  put nextPOVM
+  let fid = fidelityDM trueDM nextEstimate
+  when (outv > NoOutput) . liftIO $ do
+    let dim = 2 ^ nq
+        rankFids =
+          V.map
+            (fidelityDM trueDM . snd . getWDM . reduceParticlesToMean)
+            nextPH
+        weightsAndFids =
+          V.zip3 (V.enumFromN (1 :: Rank) dim) (V.map ptsWeight nextPH) rankFids
+    putStrLn ""
+    V.mapM_
+      (\(a, b, c) ->
+         printf "(Rank: %4d, Weight: %10.9f, Fidelity: %10.9f)\n" a b c)
+      weightsAndFids
+  return $ 1 - fid
+
+-- | Stream simulated tomography results.
+streamResults' ::
+     QBitNum -- ^ Number of quantum bits under tomography
+  -> Rank -- ^ Rank of true state
+  -> NumberOfParticles -- ^ Number of particles (per rank) to use for tomography
+  -> MHMCiter -- ^ Number of MHMC iterations to perform when resampling
+  -> OptIter -- ^ Number of POVM optimisation steps to perform
+  -> OutputVerb -- ^ Verbosity of stdout output
+  -> Stream (Of Double) IO ()
+streamResults' nq rank pn mi oi outv = do
+  let dim = 2 ^ nq
+  trueDM <- liftIO $ genDM dim rank
+  ph <- liftIO $ initialiseParticleHierarchy dim pn
+  gen <- liftIO MWC.createSystemRandom
+  rPOVM <-
+    liftIO $
+    productPOVM <$>
+    replicateM nq (mkAntipodalPOVM . fromJust . svToAngles <$> genPureSV 2)
+  let tomS = S.repeatM (simulatedTomography' trueDM nq mi oi outv gen)
+      tomS' = evalStateT (distribute tomS) rPOVM
+      initInfid = 1 - fidelityDM trueDM (getMixedEstimate ph)
+  S.yield initInfid
+  evalStateT (distribute tomS') (ph, [])
diff --git a/test/FidelityTests.hs b/test/FidelityTests.hs
new file mode 100644
--- /dev/null
+++ b/test/FidelityTests.hs
@@ -0,0 +1,24 @@
+module FidelityTests
+  ( testFidelity
+  ) where
+
+import HABQTlib.Data
+import qualified Test.QuickCheck as QC
+import TestHelpers
+
+fidProp :: QC.Positive Dim -> QC.Property
+fidProp (QC.Positive dim) =
+  let gen = do
+        v1 <- QC.resize dim QC.arbitrary
+        v2 <- QC.resize dim QC.arbitrary
+        return (v1, v2)
+      fidMatch :: (PureStateVector, PureStateVector) -> QC.Property
+      fidMatch (sv1, sv2) =
+        QC.property $
+        abs (fidelity sv1 sv2 - fidelityDM (svToDM sv1) (svToDM sv2)) <= 1e-12
+   in QC.forAll gen fidMatch
+
+testFidelity :: IO ()
+testFidelity = do
+  putStrLn "Testing density matrix fidelity:"
+  QC.quickCheckWith QC.stdArgs {QC.maxSuccess = 1000} fidProp
diff --git a/test/MeasurementTests.hs b/test/MeasurementTests.hs
new file mode 100644
--- /dev/null
+++ b/test/MeasurementTests.hs
@@ -0,0 +1,87 @@
+module MeasurementTests
+  ( testMeasurements
+  ) where
+
+import Data.Maybe (fromJust)
+import qualified Data.Vector as V
+import HABQTlib.Data
+import HABQTlib.MeasurementProcessing
+import Numeric.LinearAlgebra (Complex(..), (><))
+import qualified Numeric.LinearAlgebra as LA
+import StateGenTests (arbPropPure)
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck (Property, (.&&.))
+import TestHelpers
+
+povmValidProp :: PurePOVM -> Property
+povmValidProp wsts = normedContents .&&. sumsToUnity
+  where
+    dim :: Int
+    dim =
+      (\(WeighedPureStateVector (_, PureStateVector sv)) -> LA.rows sv) $
+      V.head wsts
+    normedContents =
+      QC.conjoin . V.toList $
+      V.map (\(WeighedPureStateVector (_, sv)) -> arbPropPure sv) wsts
+    msum :: LA.Matrix (LA.Complex Double)
+    msum =
+      V.foldl'
+        (\acc (WeighedPureStateVector (w, PureStateVector sv)) ->
+           acc + LA.scalar (w :+ 0) * (sv LA.<> LA.tr sv))
+        ((dim >< dim) (repeat (0 :+ 0)))
+        wsts
+    sumsToUnity = LA.norm_Frob (msum - LA.ident dim) <= 1e-6
+
+povmProbabilityNorm :: PurePOVM -> DensityMatrix -> Property
+povmProbabilityNorm povm dm =
+  QC.property $ abs (1 - V.sum (measurementProbs povm dm)) <= 1e-8
+
+blochHelper :: Double -> Double -> SingleQbParam
+blochHelper = (,)
+
+coordProp :: Property
+coordProp =
+  let prop sv = fidelity sv sv' > 0.99999
+        where
+          sv' = blochAnglesToSV . fromJust . svToAngles $ sv
+   in QC.forAll (QC.resize 2 QC.arbitrary) prop
+
+testBlochCoords :: IO ()
+testBlochCoords = do
+  putStrLn "Testing Bloch coordinate transofrmations:"
+  QC.quickCheckWith QC.stdArgs {QC.maxSuccess = 1000} coordProp
+
+testAntipodal1QbPOVM :: IO ()
+testAntipodal1QbPOVM = do
+  putStrLn "Testing antipodal single qubit POVM:"
+  QC.quickCheckWith QC.stdArgs {QC.maxSuccess = 1000} $
+    ((povmValidProp . mkAntipodalPOVM) .) . blochHelper
+
+genProdPOVM :: Int -> QC.Gen PurePOVM
+genProdPOVM qbnum = do
+  paramVs' <- QC.vectorOf qbnum QC.arbitrary
+  let paramVs = fmap (uncurry blochHelper) paramVs'
+      svs = map mkAntipodalPOVM paramVs
+  return $ productPOVM svs
+
+testProductPOVM :: IO ()
+testProductPOVM = do
+  putStrLn "Testing productPOVMs"
+  let gen = do
+        qbnum <- QC.suchThat QC.arbitrary (\n -> n >= 1 && n <= 6)
+        let dim :: Dim
+            dim = 2 ^ qbnum
+        rank <- QC.suchThat QC.arbitrary (\n -> n <= dim && n > 0)
+        povm <- genProdPOVM qbnum
+        dm <- arbDM dim rank
+        return (povm, dm)
+  QC.quickCheckWith QC.stdArgs {QC.maxSuccess = 10000} $
+    QC.forAll
+      gen
+      (\(povm, dm) -> povmValidProp povm .&&. povmProbabilityNorm povm dm)
+
+testMeasurements :: IO ()
+testMeasurements = do
+  testBlochCoords
+  testAntipodal1QbPOVM
+  testProductPOVM
diff --git a/test/ParticleProcessingTests.hs b/test/ParticleProcessingTests.hs
new file mode 100644
--- /dev/null
+++ b/test/ParticleProcessingTests.hs
@@ -0,0 +1,259 @@
+{-# LANGUAGE NamedFieldPuns #-}
+
+module ParticleProcessingTests
+  ( testParticleHandling
+  ) where
+
+import Control.Applicative (liftA2)
+import Data.Complex (Complex(..))
+import Data.List (sort)
+import qualified Data.Vector as V
+import HABQTlib.Data
+import HABQTlib.Data.Particle
+import HABQTlib.RandomStates
+import qualified Numeric.LinearAlgebra as LA
+import StateGenTests (arbPropDensityMatrix)
+import qualified Streaming.Prelude as S
+import qualified System.Random.MWC as MWC
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck
+  ( Args(..)
+  , Positive(..)
+  , Property
+  , (.&&.)
+  , arbitrary
+  , quickCheckWith
+  , stdArgs
+  , suchThat
+  )
+import TestHelpers
+
+lik :: Double -> Double
+lik p = p / (1 - p)
+
+loglik :: Double -> Double
+loglik = log . lik
+
+likProp :: PureStateVector -> DensityMatrix -> Property
+likProp sv dm =
+  let l = pureStateLikelihood sv dm
+   in QC.property (l >= 0) .&&. QC.property (l <= 1)
+
+testLikelihood :: IO ()
+testLikelihood = do
+  putStrLn "Testing likelihood computation for density matrices:"
+  quickCheckWith stdArgs {maxSuccess = 1000} likProp
+
+expProp :: QC.Blind Particles -> Property
+expProp (QC.Blind pts) =
+  let dm = snd . getWDM . V.foldl1' (<+>) . ptsParticles $ pts
+      dim = LA.rows . getDensityMatrix $ dm
+      z = (dim LA.>< dim) $ repeat 0
+      scaleDM :: Weight -> DensityMatrix -> LA.Matrix (Complex Double)
+      scaleDM w (DensityMatrix dm'') = LA.scale (w :+ 0) dm''
+      dm' = DensityMatrix $ foldOverPts id scaleDM (+) z pts
+      matchProp = QC.property $ dm <==> dm'
+      validProp = arbPropDensityMatrix dm'
+   in matchProp .&&. validProp
+
+testExpectation :: IO ()
+testExpectation = do
+  putStrLn "Testing particle expectation function:"
+  quickCheckWith stdArgs {maxSuccess = 500} expProp
+
+particleUpdateFidR :: PureStateVector -> Particles -> Double
+particleUpdateFidR testState@(PureStateVector sv) particles =
+  let testPure = DensityMatrix $ sv LA.<> LA.tr sv
+      originalEstimate = snd . getWDM $ reduceParticlesToMean particles
+      originalFid = fidelityDM testPure originalEstimate
+      updatedParticles = updateParticles testState particles
+      updatedEstimate = snd . getWDM $ reduceParticlesToMean updatedParticles
+      updatedFid = fidelityDM testPure updatedEstimate
+   in loglik updatedFid - loglik originalFid
+
+updateProp :: Property
+updateProp =
+  let gen = do
+        (dim, rank) <- arbDimRank
+        Positive num <- arbitrary :: QC.Gen (Positive NumberOfParticles)
+        stateVector <- QC.resize dim arbitrary
+        particles <- arbParticles dim rank (fromIntegral num)
+        return (stateVector, particles)
+      prop (sv, pts) =
+        let wtsf Particles {ptsParticles} =
+              V.foldl1' (+) . V.map (fst . getWDM) $ ptsParticles
+            ws = wtsf pts
+            pts' = updateParticles sv pts
+            ws' = wtsf pts'
+            wsProp = abs (ws - 1) < 1e-10
+            ws'Prop = abs (ws' - 1) < 1e-10
+         in wsProp .&&. ws'Prop
+      msg = "of updates move the estimate closer"
+   in QC.forAll
+        gen
+        (\(sv, pts) ->
+           QC.classify (particleUpdateFidR sv pts >= 0) msg (prop (sv, pts)))
+
+ecdfProp :: Particles -> Property
+ecdfProp Particles {ptsParticles} =
+  let cdf = ecdf ptsParticles
+      hp = V.head cdf > 0
+      tp = abs (V.last cdf - 1) < 1e-10
+      cdfSorted = V.fromList . sort . V.toList $ cdf
+      sorted = cdf == cdfSorted
+   in hp .&&. tp .&&. sorted
+
+icdfProp :: Particles -> Property
+icdfProp Particles {ptsParticles} =
+  QC.ioProperty $ do
+    gen <- MWC.createSystemRandom
+    x <- MWC.uniform gen
+    let cdf = ecdf ptsParticles
+        idx = icdf cdf x
+        Just idx' = V.findIndex (> x) cdf
+    return (idx == idx')
+
+hDistProp :: ParticleHierarchy -> PureStateVector -> Double
+hDistProp ph sv =
+  let dm = DensityMatrix $ getStateVector sv LA.<> LA.tr (getStateVector sv)
+      originalEstimate = getMixedEstimate ph
+      updatedEstimate = getMixedEstimate $ updateParticleHierarchy sv ph
+      originalFid = fidelityDM dm originalEstimate
+      updatedFid = fidelityDM dm updatedEstimate
+   in loglik updatedFid - loglik originalFid
+
+hierarchyBatchProp ::
+     Int -> Positive Dim -> Positive NumberOfParticles -> IO Property
+hierarchyBatchProp batchSize (Positive dim) (Positive num) = do
+  ph <- initialiseParticleHierarchy dim num
+  passed <-
+    S.length_ . S.filter (>= 0) . S.replicateM batchSize $
+    hDistProp ph <$> genPureSV dim
+  return . QC.property $
+    fromIntegral passed / (fromIntegral batchSize :: Double) >= 0.95
+
+hierarchyUpdateProp :: Int -> Property
+hierarchyUpdateProp batchSize =
+  let gen = do
+        dim <- suchThat arbitrary (\x -> x > 1 && x < 20)
+        num <- suchThat arbitrary (> 100)
+        return (Positive dim, Positive num)
+   in QC.forAll
+        gen
+        (\(dim, num) -> QC.ioProperty $ hierarchyBatchProp batchSize dim num)
+
+effectiveSampleSizeProp :: Property
+effectiveSampleSizeProp =
+  let genP = do
+        (dim, rank) <- arbDimRank
+        Positive num <- arbitrary :: QC.Gen (Positive NumberOfParticles)
+        arbParticles dim rank (fromIntegral num)
+      sizeP ps@(Particles _ _ n _) =
+        QC.classify
+          (es <= fromIntegral n / 2)
+          "of effective sizes are less than half of total particle number"
+          ((es <= fromIntegral n) && (es >= 0))
+        where
+          es = effectiveSize ps
+   in QC.forAll genP sizeP
+
+nudgeProp :: Property
+nudgeProp =
+  let gen = do
+        (dim, rank) <- arbDimRank
+        dm <- arbDM dim rank
+        return (dim, rank, dm)
+      nudged (dim, rank, dm) =
+        truncateRank rank <$> nudgeParticle dim 1e-2 (mkWDM1 dm)
+      prop (dim, rank, dm) =
+        QC.ioProperty $ do
+          WeighedDensityMatrix (w, ndm) <- nudged (dim, rank, dm)
+          let validDM = arbPropDensityMatrix ndm
+              preservesRank = getRank dm == getRank ndm
+              preservesWeight = w == 1
+          return $
+            QC.classify
+              (fidelityDM ndm dm > 0.99)
+              "fidelities exceed 99%"
+              (validDM .&&. preservesRank .&&. preservesWeight)
+   in QC.forAll gen prop
+
+type ResamplingFunIO = MWC.GenIO -> Particles -> IO Particles
+
+resampleProp :: MWC.GenIO -> ResamplingFunIO -> Property
+resampleProp gen resamplingFun =
+  let genA = do
+        dim <- suchThat (arbitrary :: QC.Gen Dim) (liftA2 (&&) (< 10) (> 1))
+        Positive rank <-
+          suchThat (arbitrary :: QC.Gen (Positive Rank)) (< Positive dim)
+        num <- suchThat arbitrary (> 1000)
+        QC.Blind <$> arbParticles dim rank num
+      preservesMeanIO :: QC.Blind Particles -> IO Property
+      preservesMeanIO (QC.Blind pts) = do
+        rpts <- resamplingFun gen pts
+        let WeighedDensityMatrix (w1, dm1) = reduceParticlesToMean pts
+            WeighedDensityMatrix (w2, dm2) = reduceParticlesToMean rpts
+            sumw wacc (WeighedDensityMatrix (w', _)) = wacc + w'
+            getpw Particles {ptsParticles = ps} = V.foldl' sumw 0 ps
+            w1' = getpw pts
+            w2' = getpw rpts
+            propWeight = QC.property $ 2 * abs (w1 - w2) / (w1 + w2) < 1e-12
+            propWeight' =
+              QC.property $ 2 * abs (w1' - w2') / (w1' + w2') < 1e-12
+            Particles _ _ num _ = rpts
+            propESabs =
+              QC.property $ abs (effectiveSize rpts - fromIntegral num) < 1
+            propESincrease =
+              QC.property $ effectiveSize rpts > effectiveSize pts
+            propAll =
+              propWeight .&&. propWeight' .&&. propESincrease .&&. propESabs
+        return $
+          QC.classify
+            (fidelityDM dm1 dm2 > 0.99)
+            "fidelities between pre- and post-resampling means exceed 99%"
+            propAll
+   in QC.forAll genA (QC.ioProperty . preservesMeanIO)
+
+testParticleResampling :: IO ()
+testParticleResampling = do
+  gen <- MWC.createSystemRandom
+  putStrLn "Testing multinomial particle resampling without checks:"
+  quickCheckWith stdArgs {maxSuccess = 100} $ resampleProp gen resampleMultinom
+
+testParticleNudging :: IO ()
+testParticleNudging = do
+  putStrLn "Testing particle nudging:"
+  quickCheckWith stdArgs {maxSuccess = 1000} nudgeProp
+
+testEffectiveSampleSize :: IO ()
+testEffectiveSampleSize = do
+  putStrLn "Testing calculation of effective sample size:"
+  quickCheckWith stdArgs {maxSuccess = 1000} effectiveSampleSizeProp
+
+testParticleUpdate :: IO ()
+testParticleUpdate = do
+  putStrLn "Testing update of particles:"
+  quickCheckWith stdArgs {maxSuccess = 1000} updateProp
+
+testEcdf :: IO ()
+testEcdf = do
+  putStrLn "Testing ecdf calculation:"
+  quickCheckWith stdArgs {maxSuccess = 10000} ecdfProp
+  putStrLn "Testing iecdf calculation:"
+  quickCheckWith stdArgs {maxSuccess = 10000} icdfProp
+
+testHierarchyUpdate :: IO ()
+testHierarchyUpdate = do
+  putStrLn "Testing update of particle hierarchies:"
+  quickCheckWith stdArgs {maxSuccess = 50} (hierarchyUpdateProp 100)
+
+testParticleHandling :: IO ()
+testParticleHandling = do
+  testLikelihood
+  testExpectation
+  testParticleUpdate
+  testEcdf
+  testHierarchyUpdate
+  testEffectiveSampleSize
+  testParticleNudging
+  testParticleResampling
diff --git a/test/RankReductionTests.hs b/test/RankReductionTests.hs
new file mode 100644
--- /dev/null
+++ b/test/RankReductionTests.hs
@@ -0,0 +1,81 @@
+module RankReductionTests
+  ( testRankReduction
+  ) where
+
+import Data.List (sort)
+import HABQTlib.Data
+import HABQTlib.Data.Particle
+import qualified Numeric.LinearAlgebra as LA
+import StateGenTests (arbPropDensityMatrix)
+import qualified System.Random.MWC as MWC
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck
+  ( Args(..)
+  , Positive(..)
+  , Property
+  , (.&&.)
+  , (===)
+  , ioProperty
+  , quickCheckWith
+  , stdArgs
+  )
+import TestHelpers
+
+sortSVDprop :: Positive Int -> Positive Int -> Property
+sortSVDprop (Positive r) (Positive c) =
+  ioProperty $ do
+    r1 <- LA.randn r c
+    r2 <- LA.randn r c
+    let rm = LA.toComplex (r1, r2)
+        (u, s) = LA.leftSV rm
+        sorted = sort (LA.toList s) == (reverse . LA.toList) s
+        normU = (LA.norm_2 . head . LA.toColumns) u
+    return $ sorted && abs (normU - 1) < 1e-6
+
+truncateProp' :: Rank -> WeighedDensityMatrix -> Property
+truncateProp' rank wdm@(WeighedDensityMatrix (w, _)) =
+  let WeighedDensityMatrix (wt, dmt) = truncateRank rank wdm
+      preservesWeight = abs (2 * (w - wt) / (w + wt)) <= 1e-12
+      (_, s, _) = LA.compactSVD . getDensityMatrix $ dmt
+      setsRank = rank === LA.size s
+      validDM = arbPropDensityMatrix dmt
+   in preservesWeight .&&. setsRank .&&. validDM
+
+truncateProp :: Property
+truncateProp =
+  let gen = do
+        (dim, rank) <- arbDimRank
+        wdm <- arbWDM dim rank
+        newRank <- QC.suchThat QC.arbitrary (\x -> x > 0 && x <= rank)
+        return (newRank, wdm)
+   in QC.forAll gen (uncurry truncateProp')
+
+testTruncation :: IO ()
+testTruncation = do
+  putStrLn "Testing density matrix rank reduction:"
+  quickCheckWith stdArgs {maxSuccess = 1000} truncateProp
+
+particleReductionProp :: MWC.GenIO -> Positive Int -> Positive Int -> Property
+particleReductionProp gen (Positive dim) (Positive num) =
+  ioProperty $ do
+    r <- MWC.uniformR (1, dim) gen
+    particles <- genParticles dim r num
+    let WeighedDensityMatrix (w, dm) = reduceParticlesToMean particles
+        dimProp = dim === LA.rows (getDensityMatrix dm)
+        rankProp = r === getRank dm
+        validProp = arbPropDensityMatrix dm
+        weightProp = abs (w - 1) < 1e-12
+    return $ dimProp .&&. rankProp .&&. validProp .&&. weightProp
+
+testSVDsort :: IO ()
+testSVDsort = do
+  putStrLn "Testing singular value sort:"
+  quickCheckWith stdArgs {maxSuccess = 1000} sortSVDprop
+  putStrLn "Testing particle vector reduction:"
+  gen <- MWC.createSystemRandom
+  quickCheckWith stdArgs {maxSuccess = 1000} (particleReductionProp gen)
+
+testRankReduction :: IO ()
+testRankReduction = do
+  testSVDsort
+  testTruncation
diff --git a/test/StateGenTests.hs b/test/StateGenTests.hs
new file mode 100644
--- /dev/null
+++ b/test/StateGenTests.hs
@@ -0,0 +1,72 @@
+module StateGenTests
+  ( testStateGen
+  , testStateArb
+  , arbPropPure
+  , arbPropDensityMatrix
+  ) where
+
+import Data.Complex (magnitude)
+import HABQTlib.Data
+import HABQTlib.RandomStates
+import Numeric.LinearAlgebra (norm_2, sumElements, takeDiag, tr)
+import qualified System.Random.MWC as MWC
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck (Property)
+import TestHelpers
+
+class StateProp a where
+  traceHermProp :: MWC.GenIO -> a -> Property
+
+dimRankIO :: MWC.GenIO -> Dim -> IO (Dim, Rank)
+dimRankIO gen ub = do
+  dim <- MWC.uniformR (1, ub) gen
+  r <- MWC.uniformR (1, dim) gen
+  return (dim, r)
+
+instance StateProp DensityMatrix where
+  traceHermProp gen _ =
+    QC.ioProperty $ do
+      (dim, r) <- dimRankIO gen 100
+      dm <- getDensityMatrix <$> genDM dim r
+      return $
+        (abs (1 - (magnitude . sumElements . takeDiag) dm) < 1e-12) &&
+        (norm_2 (dm - tr dm) < 1e-12)
+
+instance StateProp PureStateVector where
+  traceHermProp gen _ =
+    QC.ioProperty $ do
+      (dim, _) <- dimRankIO gen 100
+      sv <- getStateVector <$> genPureSV dim
+      return $ abs (1 - norm_2 sv) < 1e-12
+
+testStateGen :: IO ()
+testStateGen = do
+  gen <- MWC.createSystemRandom
+  putStrLn "Testing density matrix generation:"
+  QC.quickCheckWith
+    QC.stdArgs {QC.maxSuccess = 10000}
+    (traceHermProp gen :: DensityMatrix -> Property)
+  putStrLn "Testing pure state generation:"
+  QC.quickCheckWith
+    QC.stdArgs {QC.maxSuccess = 10000}
+    (traceHermProp gen :: PureStateVector -> Property)
+
+arbPropPure :: PureStateVector -> Property
+arbPropPure (PureStateVector sv) = QC.property $ abs (1 - norm_2 sv) < 1e-12
+
+arbPropDensityMatrix :: DensityMatrix -> Property
+arbPropDensityMatrix (DensityMatrix dm) =
+  QC.property $
+  (abs (1 - (magnitude . sumElements . takeDiag) dm) < 1e-12) &&
+  (norm_2 (dm - tr dm) < 1e-12)
+
+testStateArb :: IO ()
+testStateArb = do
+  putStrLn "Testing density matrix arbitrary:"
+  QC.quickCheckWith
+    QC.stdArgs {QC.maxSuccess = 1000}
+    (arbPropDensityMatrix :: DensityMatrix -> Property)
+  putStrLn "Testing pure state arbitrary:"
+  QC.quickCheckWith
+    QC.stdArgs {QC.maxSuccess = 10000}
+    (arbPropPure :: PureStateVector -> Property)
diff --git a/test/SuperpositionSemigroupTests.hs b/test/SuperpositionSemigroupTests.hs
new file mode 100644
--- /dev/null
+++ b/test/SuperpositionSemigroupTests.hs
@@ -0,0 +1,37 @@
+module SuperpositionSemigroupTests
+  ( testWeighedDensityMatrixSemigroup
+  ) where
+
+import HABQTlib.Data
+import Test.QuickCheck
+  ( Arbitrary(..)
+  , Args(..)
+  , Gen
+  , Property
+  , forAll
+  , getSize
+  , quickCheckWith
+  , resize
+  , stdArgs
+  , suchThat
+  )
+import TestHelpers
+
+type WDM = WeighedDensityMatrix
+
+testAssocWithGen :: Gen (WDM, WDM, WDM) -> Property
+testAssocWithGen gen =
+  forAll gen (\(x, y, z) -> ((x <+> y) <+> z) <==> (x <+> (y <+> z)))
+
+gen3DM :: Gen (WDM, WDM, WDM)
+gen3DM = do
+  n <- suchThat getSize (> 0)
+  dm1 <- resize n arbitrary
+  dm2 <- resize n arbitrary
+  dm3 <- resize n arbitrary
+  return (dm1, dm2, dm3)
+
+testWeighedDensityMatrixSemigroup :: IO ()
+testWeighedDensityMatrixSemigroup = do
+  putStrLn "Testing weighed density matrix semigroup:"
+  quickCheckWith stdArgs {maxSuccess = 1000} (testAssocWithGen gen3DM)
diff --git a/test/TestHelpers.hs b/test/TestHelpers.hs
new file mode 100644
--- /dev/null
+++ b/test/TestHelpers.hs
@@ -0,0 +1,107 @@
+module TestHelpers where
+
+import Control.Newtype.Generics (over)
+import Data.Complex
+import qualified Data.Vector as V
+import HABQTlib.Data
+import HABQTlib.Data.Particle
+import qualified Numeric.LinearAlgebra as LA
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck (Arbitrary(..), Gen, Positive(..))
+
+getRank :: DensityMatrix -> Rank
+getRank (DensityMatrix dm) =
+  let (_, s, _) = LA.compactSVD dm
+   in LA.size s
+
+arbDimRank :: QC.Gen (Dim, Rank)
+arbDimRank = do
+  dim <- QC.suchThat arbitrary (> 1)
+  Positive rank <- QC.suchThat arbitrary (\(Positive v) -> v <= dim)
+  return (dim, rank)
+
+arbDM :: Dim -> Rank -> QC.Gen DensityMatrix
+arbDM dim rank = do
+  let genXs =
+        QC.vectorOf
+          dim
+          (QC.vectorOf rank (QC.arbitrary :: QC.Gen (Complex Double)))
+  xs <- QC.suchThat genXs (\ll -> any (/= 0.0 :+ 0.0) (zipWith (!!) ll [0 ..]))
+  let a = LA.fromLists xs
+      h = a LA.<> LA.tr a
+      hTr = LA.sumElements $ LA.takeDiag h
+      dm = h / LA.scalar hTr
+  return $ DensityMatrix dm
+
+arbWDM :: Dim -> Rank -> QC.Gen WeighedDensityMatrix
+arbWDM dim rank = do
+  Positive w <- QC.arbitrary
+  dm <- arbDM dim rank
+  return $ WeighedDensityMatrix (w, dm)
+
+arbNormedVecDim :: Int -> Gen (LA.Matrix (Complex Double))
+arbNormedVecDim dim = do
+  let genXs = QC.vectorOf dim (arbitrary :: Gen (Complex Double))
+  xs <- QC.suchThat genXs (any (/= 0))
+  let sv = LA.asColumn . LA.fromList $ xs
+  return $ sv / LA.toComplex (LA.scalar (LA.norm_2 sv), 0)
+
+arbitraryWeighed :: Arbitrary x => Gen (Weight, x)
+arbitraryWeighed = do
+  dim <- QC.suchThat QC.getSize (> 0)
+  w <- getPositive <$> arbitrary
+  dm <- QC.resize dim arbitrary
+  return (w, dm)
+
+arbParticles :: Dim -> Rank -> NumberOfParticles -> QC.Gen Particles
+arbParticles dim rank num = do
+  vdms <- V.replicateM num (arbWDM dim rank)
+  let wr = V.foldl' (\acc (WeighedDensityMatrix (wi, _)) -> acc + wi) 0 vdms
+      vdmsn = V.map (over WeighedDensityMatrix (\(w, dm) -> (w / wr, dm))) vdms
+  QC.Positive w <- QC.arbitrary
+  return $ Particles rank w num vdmsn
+
+instance Arbitrary DensityMatrix where
+  arbitrary = do
+    dim <- QC.suchThat QC.getSize (> 0)
+    arbDM dim dim
+
+instance Arbitrary PureStateVector where
+  arbitrary = do
+    dim <- QC.suchThat QC.getSize (> 0)
+    sv <- arbNormedVecDim dim
+    return . PureStateVector $ sv
+
+instance Arbitrary WeighedDensityMatrix where
+  arbitrary = do
+    (w, dm) <- arbitraryWeighed
+    return . WeighedDensityMatrix $ (w, dm)
+
+instance Arbitrary WeighedPureStateVector where
+  arbitrary = do
+    (w, sv) <- arbitraryWeighed
+    return . WeighedPureStateVector $ (w, sv)
+
+instance QC.Arbitrary Particles where
+  arbitrary = do
+    (dim, rank) <- arbDimRank
+    num <- QC.suchThat QC.arbitrary (> 100)
+    arbParticles dim rank num
+
+class MEq a where
+  infix 4 <==>
+  (<==>) :: a -> a -> Bool
+
+instance MEq DensityMatrix where
+  DensityMatrix dm1 <==> DensityMatrix dm2 = LA.norm_Frob (dm1 - dm2) < 1e-6
+
+instance MEq WeighedDensityMatrix where
+  WeighedDensityMatrix (w1, DensityMatrix dm1) <==> WeighedDensityMatrix (w2, DensityMatrix dm2) =
+    LA.norm_Frob (LA.scalar (w1 :+ 0) * dm1 - LA.scalar (w2 :+ 0) * dm2) < 1e-6
+
+instance MEq PureStateVector where
+  sv1 <==> sv2 = 1 - fidelity sv1 sv2 < 1e-6
+
+instance MEq WeighedPureStateVector where
+  WeighedPureStateVector (w1, sv1) <==> WeighedPureStateVector (w2, sv2) =
+    abs (w1 - w2) <= 1e-6 * 0.5 * (w1 + w2) && sv1 <==> sv2
diff --git a/test/Tests.hs b/test/Tests.hs
new file mode 100644
--- /dev/null
+++ b/test/Tests.hs
@@ -0,0 +1,17 @@
+import FidelityTests
+import MeasurementTests
+import ParticleProcessingTests
+import RankReductionTests
+import StateGenTests
+import SuperpositionSemigroupTests
+
+main :: IO ()
+main = do
+  putStrLn ""
+  testStateGen
+  testStateArb
+  testFidelity
+  testWeighedDensityMatrixSemigroup
+  testRankReduction
+  testParticleHandling
+  testMeasurements