HABQT (empty) → 0.1.0.0
raw patch · 23 files changed
+2137/−0 lines, 23 filesdep +HABQTdep +QuickCheckdep +basesetup-changed
Dependencies added: HABQT, QuickCheck, base, hmatrix, hmatrix-gsl, mtl, mwc-random, newtype-generics, optparse-applicative, streaming, utility-ht, validation, vector
Files
- ChangeLog.md +3/−0
- HABQT.cabal +129/−0
- LICENSE +30/−0
- README.md +63/−0
- Setup.hs +3/−0
- app/Main.hs +181/−0
- libHABQT/ForeignHABQT.hs +54/−0
- libHABQT/LibHABQT.hs +108/−0
- libHABQT/hsinit.c +17/−0
- src/HABQTlib.hs +80/−0
- src/HABQTlib/Data.hs +226/−0
- src/HABQTlib/Data/Particle.hs +252/−0
- src/HABQTlib/MeasurementProcessing.hs +151/−0
- src/HABQTlib/RandomStates.hs +40/−0
- src/HABQTlib/UnsafeAPI.hs +116/−0
- test/FidelityTests.hs +24/−0
- test/MeasurementTests.hs +87/−0
- test/ParticleProcessingTests.hs +259/−0
- test/RankReductionTests.hs +81/−0
- test/StateGenTests.hs +72/−0
- test/SuperpositionSemigroupTests.hs +37/−0
- test/TestHelpers.hs +107/−0
- test/Tests.hs +17/−0
+ ChangeLog.md view
@@ -0,0 +1,3 @@+# Changelog for HABQT++## Unreleased changes
+ HABQT.cabal view
@@ -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
+ LICENSE view
@@ -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.
+ README.md view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,3 @@+import Distribution.Simple++main = defaultMain
+ app/Main.hs view
@@ -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")
+ libHABQT/ForeignHABQT.hs view
@@ -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
+ libHABQT/LibHABQT.hs view
@@ -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
+ libHABQT/hsinit.c view
@@ -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();+}
+ src/HABQTlib.hs view
@@ -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
+ src/HABQTlib/Data.hs view
@@ -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)
+ src/HABQTlib/Data/Particle.hs view
@@ -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
+ src/HABQTlib/MeasurementProcessing.hs view
@@ -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
+ src/HABQTlib/RandomStates.hs view
@@ -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)
+ src/HABQTlib/UnsafeAPI.hs view
@@ -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, [])
+ test/FidelityTests.hs view
@@ -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
+ test/MeasurementTests.hs view
@@ -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
+ test/ParticleProcessingTests.hs view
@@ -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
+ test/RankReductionTests.hs view
@@ -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
+ test/StateGenTests.hs view
@@ -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)
+ test/SuperpositionSemigroupTests.hs view
@@ -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)
+ test/TestHelpers.hs view
@@ -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
+ test/Tests.hs view
@@ -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