{-# LANGUAGE DeriveGeneric #-}
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)
type Dim = Int
type Rank = Int
type NumberOfParticles = Int
type QBitNum = Int
type MHMCiter = Int
type OptIter = Int
type Weight = Double
data OutputVerb
= NoOutput
| FidOutput
| FullOutput
deriving (Eq, Show, Ord)
newtype DensityMatrix = DensityMatrix
{ getDensityMatrix :: Matrix (Complex Double)
} deriving (Eq, Show, Generic)
newtype PureStateVector = PureStateVector
{ getStateVector :: Matrix (Complex Double)
} deriving (Eq, Show, Generic)
instance Newtype DensityMatrix
instance Newtype PureStateVector
validSV :: PureStateVector -> Validation [String] PureStateVector
validSV =
validate
["State vector must have unit norm."]
(\x -> abs (1 - LA.norm_2 (unpack x)) < 1e-12)
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."]
validQBitN :: QBitNum -> Validation [String] QBitNum
validQBitN = validate qnM (> 0)
miM :: [String]
miM = ["Number of MHMC iterations must be 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."]
validPartNum :: NumberOfParticles -> Validation [String] NumberOfParticles
validPartNum = validate pnM (> 0)
rM :: [String]
rM = ["Rank must be a positive integer."]
validRank :: Rank -> Validation [String] Rank
validRank = validate rM (> 0)
oiM :: [String]
oiM = ["Number of POVM optimisation iterations must be a positive integer."]
validOptIter :: OptIter -> Validation [String] OptIter
validOptIter = validate oiM (> 0)
newtype WeighedDensityMatrix = WeighedDensityMatrix
{ getWDM :: (Weight, DensityMatrix)
} deriving (Eq, Show, Generic)
mkWDM :: Weight -> DensityMatrix -> WeighedDensityMatrix
mkWDM w dm = WeighedDensityMatrix (w, dm)
mkWDM1 :: DensityMatrix -> WeighedDensityMatrix
mkWDM1 = mkWDM 1
newtype WeighedPureStateVector = WeighedPureStateVector
{ getWSV :: (Weight, PureStateVector)
} deriving (Eq, Show, Generic)
instance Newtype WeighedDensityMatrix
instance Newtype WeighedPureStateVector
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)]
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)
svToDM :: PureStateVector -> DensityMatrix
svToDM (PureStateVector sv) = DensityMatrix $ sv LA.<> LA.tr sv
infix 8 <+>
(<+>) :: 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
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)]
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)