module Mcmc.Proposal.Hamiltonian.Masses
( Mu,
MassesI,
toGMatrix,
cleanMatrix,
getMassesI,
getMus,
Dimension,
vectorToMasses,
massesToVector,
tuneDiagonalMassesOnly,
tuneAllMasses,
)
where
import Data.Maybe
import qualified Data.Vector as VB
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Unboxed as VU
import Mcmc.Proposal.Hamiltonian.Common
import qualified Numeric.LinearAlgebra as L
import qualified Statistics.Covariance as S
import qualified Statistics.Function as S
import qualified Statistics.Sample as S
type Mu = L.Vector Double
type MassesI = L.GMatrix
precision :: Double
precision :: Double
precision = Double
1e-8
isDiag :: L.Matrix Double -> Bool
isDiag :: Matrix Double -> Bool
isDiag Matrix Double
xs = Double -> Double
forall a. Num a => a -> a
abs (Double
sumDiag Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sumFull) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
precision
where
xsAbs :: Matrix Double
xsAbs = (Double -> Double) -> Matrix Double -> Matrix Double
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
L.cmap Double -> Double
forall a. Num a => a -> a
abs Matrix Double
xs
sumDiag :: Double
sumDiag = Vector Double -> Double
forall (c :: * -> *) e. Container c e => c e -> e
L.sumElements (Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
L.takeDiag Matrix Double
xsAbs)
sumFull :: Double
sumFull = Matrix Double -> Double
forall (c :: * -> *) e. Container c e => c e -> e
L.sumElements Matrix Double
xsAbs
isSparse :: L.Matrix Double -> Bool
isSparse :: Matrix Double -> Bool
isSparse Matrix Double
xs = Double
nNonZero Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nMax
where
n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows Matrix Double
xs
m :: Int
m = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
5 Int
n
nMax :: Int
nMax = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
m
f :: Double -> Double
f Double
x = if Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
precision then Double
1 else Double
0 :: Double
xsInd :: Matrix Double
xsInd = (Double -> Double) -> Matrix Double -> Matrix Double
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
L.cmap Double -> Double
f Matrix Double
xs
nNonZero :: Double
nNonZero = Matrix Double -> Double
forall (c :: * -> *) e. Container c e => c e -> e
L.sumElements Matrix Double
xsInd
toAssocMatrix :: L.Matrix Double -> L.AssocMatrix
toAssocMatrix :: Matrix Double -> AssocMatrix
toAssocMatrix Matrix Double
xs
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m = [Char] -> [((Int, Int), Double)]
forall a. HasCallStack => [Char] -> a
error [Char]
"toAssocMatrix: Matrix not square."
| Bool
otherwise =
[ ((Int
i, Int
j), Double
e)
| Int
i <- [Int
0 .. (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)],
Int
j <- [Int
0 .. (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)],
let e :: Double
e = Matrix Double
xs Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`L.atIndex` (Int
i, Int
j),
Double -> Double
forall a. Num a => a -> a
abs Double
e Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
precision
]
where
n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows Matrix Double
xs
m :: Int
m = Matrix Double -> Int
forall t. Matrix t -> Int
L.cols Matrix Double
xs
toGMatrix :: L.Matrix Double -> L.GMatrix
toGMatrix :: Matrix Double -> GMatrix
toGMatrix Matrix Double
xs
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [Char] -> GMatrix
forall a. HasCallStack => [Char] -> a
error [Char]
"toGMatrix: Matrix empty."
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m = [Char] -> GMatrix
forall a. HasCallStack => [Char] -> a
error [Char]
"toGMatrix: Matrix not square."
| Matrix Double -> Bool
isDiag Matrix Double
xs = Int -> Int -> Vector Double -> GMatrix
L.mkDiagR Int
n Int
m (Vector Double -> GMatrix) -> Vector Double -> GMatrix
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
L.takeDiag Matrix Double
xs
| Matrix Double -> Bool
isSparse Matrix Double
xs = AssocMatrix -> GMatrix
L.mkSparse (AssocMatrix -> GMatrix) -> AssocMatrix -> GMatrix
forall a b. (a -> b) -> a -> b
$ Matrix Double -> AssocMatrix
toAssocMatrix Matrix Double
xs
| Bool
otherwise = Matrix Double -> GMatrix
L.mkDense Matrix Double
xs
where
n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows Matrix Double
xs
m :: Int
m = Matrix Double -> Int
forall t. Matrix t -> Int
L.cols Matrix Double
xs
cleanMatrix :: L.Matrix Double -> L.Matrix Double
cleanMatrix :: Matrix Double -> Matrix Double
cleanMatrix Matrix Double
xs =
(Vector Double -> Matrix Double
forall a. (Num a, Element a) => Vector a -> Matrix a
L.diag (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double) -> Vector Double -> Vector Double
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
L.cmap Double -> Double
cleanDiag Vector Double
xsDiag) Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ (Double -> Double) -> Matrix Double -> Matrix Double
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
L.cmap Double -> Double
cleanOffDiag Matrix Double
xsOffDiag
where
xsDiag :: Vector Double
xsDiag = Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
L.takeDiag Matrix Double
xs
cleanDiag :: Double -> Double
cleanDiag Double
x
| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x = Double
1
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = Double
1
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
precision = Double
precision
| Bool
otherwise = Double
x
xsOffDiag :: Matrix Double
xsOffDiag = Matrix Double
xs Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
- Vector Double -> Matrix Double
forall a. (Num a, Element a) => Vector a -> Matrix a
L.diag Vector Double
xsDiag
cleanOffDiag :: Double -> Double
cleanOffDiag Double
x
| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x = Double
0
| Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
precision = Double
0
| Bool
otherwise = Double
x
getMassesI :: L.Herm Double -> L.GMatrix
getMassesI :: Herm Double -> GMatrix
getMassesI Herm Double
xs
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [Char] -> GMatrix
forall a. HasCallStack => [Char] -> a
error [Char]
"getMassesI: Matrix empty."
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m = [Char] -> GMatrix
forall a. HasCallStack => [Char] -> a
error [Char]
"getMassesI: Matrix not square."
| Double
sign Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
1.0 = [Char] -> GMatrix
forall a. HasCallStack => [Char] -> a
error [Char]
"getMassesI: Determinant of matrix is negative."
| Bool
otherwise = Matrix Double -> GMatrix
toGMatrix (Matrix Double -> GMatrix) -> Matrix Double -> GMatrix
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
cleanMatrix (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Matrix Double
xsI
where
xs' :: Matrix Double
xs' = Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym (Herm Double -> Matrix Double) -> Herm Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Herm Double
xs
n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows Matrix Double
xs'
m :: Int
m = Matrix Double -> Int
forall t. Matrix t -> Int
L.cols Matrix Double
xs'
(Matrix Double
xsI, (Double
_, Double
sign)) = Matrix Double -> (Matrix Double, (Double, Double))
forall t. Field t => Matrix t -> (Matrix t, (t, t))
L.invlndet Matrix Double
xs'
getMus :: Masses -> L.Vector Double
getMus :: Herm Double -> Vector Double
getMus Herm Double
xs
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [Char] -> Vector Double
forall a. HasCallStack => [Char] -> a
error [Char]
"getMu: Matrix empty."
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m = [Char] -> Vector Double
forall a. HasCallStack => [Char] -> a
error [Char]
"getMu: Matrix not square."
| Bool
otherwise = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
L.fromList ([Double] -> Vector Double) -> [Double] -> Vector Double
forall a b. (a -> b) -> a -> b
$ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
0.0
where
xs' :: Matrix Double
xs' = Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym Herm Double
xs
n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows Matrix Double
xs'
m :: Int
m = Matrix Double -> Int
forall t. Matrix t -> Int
L.cols Matrix Double
xs'
type Dimension = Int
massesToVector :: Masses -> VU.Vector Double
massesToVector :: Herm Double -> Vector Double
massesToVector = Vector Double -> Vector Double
forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
VU.convert (Vector Double -> Vector Double)
-> (Herm Double -> Vector Double) -> Herm Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
L.flatten (Matrix Double -> Vector Double)
-> (Herm Double -> Matrix Double) -> Herm Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym
vectorToMasses :: Dimension -> VU.Vector Double -> Masses
vectorToMasses :: Int -> Vector Double -> Herm Double
vectorToMasses Int
d = Matrix Double -> Herm Double
forall t. Matrix t -> Herm t
L.trustSym (Matrix Double -> Herm Double)
-> (Vector Double -> Matrix Double) -> Vector Double -> Herm Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Vector Double -> Matrix Double
forall t. Storable t => Int -> Vector t -> Matrix t
L.reshape Int
d (Vector Double -> Matrix Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Vector Double
forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
VU.convert
massMin :: Double
massMin :: Double
massMin = Double
precision
massMax :: Double
massMax :: Double
massMax = Double -> Double
forall a. Fractional a => a -> a
recip Double
precision
samplesDiagonalMin :: Int
samplesDiagonalMin :: Int
samplesDiagonalMin = Int
61
samplesAllMinWith :: Dimension -> Int
samplesAllMinWith :: Int -> Int
samplesAllMinWith Int
d = Int
samplesDiagonalMin Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
samplesDiagonalMin Int
d
getSampleSize :: VS.Vector Double -> Int
getSampleSize :: Vector Double -> Int
getSampleSize = Vector Double -> Int
forall a. Storable a => Vector a -> Int
VS.length (Vector Double -> Int)
-> (Vector Double -> Vector Double) -> Vector Double -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Vector Double
forall a. (Storable a, Eq a) => Vector a -> Vector a
VS.uniq (Vector Double -> Vector Double)
-> (Vector Double -> Vector Double)
-> Vector Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Vector Double
forall e (v :: * -> *). (Ord e, Vector v e) => v e -> v e
S.gsort
getNewMassDiagonalWithRescue :: Int -> Double -> Double -> Double
getNewMassDiagonalWithRescue :: Int -> Double -> Double -> Double
getNewMassDiagonalWithRescue Int
sampleSize Double
massOld Double
massEstimate
| Int
sampleSize Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
samplesDiagonalMin = Double
massOld
| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
massEstimate = Double
massOld
| Double
massEstimate Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 = Double
massOld
| Double
massMin Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
massNew = Double
massMin
| Double
massNew Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
massMax = Double
massMax
| Bool
otherwise = Double
massNew
where
massNewSqrt :: Double
massNewSqrt = Double -> Double
forall a. Fractional a => a -> a
recip Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
sqrt Double
massOld Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
massEstimate)
massNew :: Double
massNew = Double
massNewSqrt Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
2
findClosestPositiveDefiniteMatrix :: L.Matrix Double -> L.Matrix Double
findClosestPositiveDefiniteMatrix :: Matrix Double -> Matrix Double
findClosestPositiveDefiniteMatrix Matrix Double
a
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = [Char] -> Matrix Double
forall a. HasCallStack => [Char] -> a
error [Char]
"findClosestPositiveDefiniteMatrix: Matrix empty."
| Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m = [Char] -> Matrix Double
forall a. HasCallStack => [Char] -> a
error [Char]
"findClosestPositiveDefiniteMatrix: Matrix not square."
| Matrix Double -> Bool
isPositiveDefinite Matrix Double
a = Matrix Double
a
| Bool
otherwise = Matrix Double -> Double -> Matrix Double
go Matrix Double
a3 Double
1
where
n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows Matrix Double
a
m :: Int
m = Matrix Double -> Int
forall t. Matrix t -> Int
L.cols Matrix Double
a
b :: Matrix Double
b = Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym (Herm Double -> Matrix Double) -> Herm Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Herm Double
forall t. Field t => Matrix t -> Herm t
L.sym Matrix Double
a
(Matrix Double
_, Vector Double
s, Matrix Double
v) = Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
forall t.
Field t =>
Matrix t -> (Matrix t, Vector Double, Matrix t)
L.svd Matrix Double
b
h :: Matrix Double
h = (Matrix Double -> Matrix Double
forall m mt. Transposable m mt => m -> mt
L.tr Matrix Double
v) Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
L.<> (Vector Double -> Matrix Double
forall a. (Num a, Element a) => Vector a -> Matrix a
L.diag Vector Double
s Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
L.<> Matrix Double
v)
a2 :: Matrix Double
a2 = Double -> Matrix Double -> Matrix Double
forall t (c :: * -> *). Linear t c => t -> c t -> c t
L.scale Double
0.5 (Matrix Double
b Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ Matrix Double
h)
a3 :: Matrix Double
a3 = Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym (Herm Double -> Matrix Double) -> Herm Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Herm Double
forall t. Field t => Matrix t -> Herm t
L.sym Matrix Double
a2
isPositiveDefinite :: Matrix Double -> Bool
isPositiveDefinite = Maybe (Matrix Double) -> Bool
forall a. Maybe a -> Bool
isJust (Maybe (Matrix Double) -> Bool)
-> (Matrix Double -> Maybe (Matrix Double))
-> Matrix Double
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Herm Double -> Maybe (Matrix Double)
forall t. Field t => Herm t -> Maybe (Matrix t)
L.mbChol (Herm Double -> Maybe (Matrix Double))
-> (Matrix Double -> Herm Double)
-> Matrix Double
-> Maybe (Matrix Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> Herm Double
forall t. Matrix t -> Herm t
L.trustSym
i :: Matrix Double
i = Int -> Matrix Double
forall a. (Num a, Element a) => Int -> Matrix a
L.ident Int
n
eps :: Double
eps = Double
2.2204460492503131e-16
go :: Matrix Double -> Double -> Matrix Double
go Matrix Double
x Double
k
| Matrix Double -> Bool
isPositiveDefinite Matrix Double
x = Matrix Double
x
| Bool
otherwise =
let minEig :: Double
minEig = Vector Double -> Double
forall (c :: * -> *) e. Container c e => c e -> e
L.minElement (Vector Double -> Double) -> Vector Double -> Double
forall a b. (a -> b) -> a -> b
$ (Complex Double -> Double)
-> Vector (Complex Double) -> Vector Double
forall b (c :: * -> *) e.
(Element b, Container c e) =>
(e -> b) -> c e -> c b
L.cmap Complex Double -> Double
forall a. Complex a -> a
L.realPart (Vector (Complex Double) -> Vector Double)
-> Vector (Complex Double) -> Vector Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Vector (Complex Double)
forall t. Field t => Matrix t -> Vector (Complex Double)
L.eigenvalues Matrix Double
x
nu :: Double
nu = Double -> Double
forall a. Num a => a -> a
negate Double
minEig Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eps
x' :: Matrix Double
x' = Matrix Double
x Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ Double -> Matrix Double -> Matrix Double
forall t (c :: * -> *). Linear t c => t -> c t -> c t
L.scale Double
nu Matrix Double
i
in Matrix Double -> Double -> Matrix Double
go Matrix Double
x' (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
tuneDiagonalMassesOnly ::
(a -> Positions) ->
VB.Vector a ->
(Masses, MassesI) ->
(Masses, MassesI)
tuneDiagonalMassesOnly :: (a -> Vector Double)
-> Vector a -> (Herm Double, GMatrix) -> (Herm Double, GMatrix)
tuneDiagonalMassesOnly a -> Vector Double
toVec Vector a
xs (Herm Double
ms, GMatrix
msI)
| Vector a -> Int
forall a. Vector a -> Int
VB.length Vector a
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
samplesDiagonalMin = (Herm Double
ms, GMatrix
msI)
| Int
dimState Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
dimMs = [Char] -> (Herm Double, GMatrix)
forall a. HasCallStack => [Char] -> a
error [Char]
"tuneDiagonalMassesOnly: Dimension mismatch."
| Bool
otherwise =
let msDirty :: Matrix Double
msDirty = Matrix Double
msOld Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
- Vector Double -> Matrix Double
forall a. (Num a, Element a) => Vector a -> Matrix a
L.diag Vector Double
msDiagonalOld Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ Vector Double -> Matrix Double
forall a. (Num a, Element a) => Vector a -> Matrix a
L.diag Vector Double
msDiagonalNew
ms' :: Herm Double
ms' = Matrix Double -> Herm Double
forall t. Matrix t -> Herm t
L.trustSym (Matrix Double -> Herm Double) -> Matrix Double -> Herm Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
findClosestPositiveDefiniteMatrix (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
cleanMatrix Matrix Double
msDirty
msI' :: GMatrix
msI' = Herm Double -> GMatrix
getMassesI Herm Double
ms'
in (Herm Double
ms', GMatrix
msI')
where
xs' :: Vector (Vector Double)
xs' = (a -> Vector Double) -> Vector a -> Vector (Vector Double)
forall a b. (a -> b) -> Vector a -> Vector b
VB.map a -> Vector Double
toVec Vector a
xs
xs'' :: Matrix Double
xs'' = [Vector Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
L.fromRows ([Vector Double] -> Matrix Double)
-> [Vector Double] -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Vector (Vector Double) -> [Vector Double]
forall a. Vector a -> [a]
VB.toList Vector (Vector Double)
xs'
dimState :: Int
dimState = Vector Double -> Int
forall a. Storable a => Vector a -> Int
VS.length (Vector Double -> Int) -> Vector Double -> Int
forall a b. (a -> b) -> a -> b
$ Vector (Vector Double) -> Vector Double
forall a. Vector a -> a
VB.head Vector (Vector Double)
xs'
sampleSizes :: Vector Int
sampleSizes = [Int] -> Vector Int
forall a. Storable a => [a] -> Vector a
VS.fromList ([Int] -> Vector Int) -> [Int] -> Vector Int
forall a b. (a -> b) -> a -> b
$ (Vector Double -> Int) -> [Vector Double] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Vector Double -> Int
getSampleSize ([Vector Double] -> [Int]) -> [Vector Double] -> [Int]
forall a b. (a -> b) -> a -> b
$ Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
L.toColumns Matrix Double
xs''
msOld :: Matrix Double
msOld = Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym Herm Double
ms
dimMs :: Int
dimMs = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows Matrix Double
msOld
msDiagonalOld :: Vector Double
msDiagonalOld = Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
L.takeDiag Matrix Double
msOld
msDiagonalEstimate :: Vector Double
msDiagonalEstimate = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
VS.fromList ([Double] -> Vector Double) -> [Double] -> Vector Double
forall a b. (a -> b) -> a -> b
$ (Vector Double -> Double) -> [Vector Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double
forall a. Fractional a => a -> a
recip (Double -> Double)
-> (Vector Double -> Double) -> Vector Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> Double
forall (v :: * -> *). Vector v Double => v Double -> Double
S.variance) ([Vector Double] -> [Double]) -> [Vector Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
L.toColumns Matrix Double
xs''
msDiagonalNew :: Vector Double
msDiagonalNew =
(Int -> Double -> Double -> Double)
-> Vector Int -> Vector Double -> Vector Double -> Vector Double
forall a b c d.
(Storable a, Storable b, Storable c, Storable d) =>
(a -> b -> c -> d) -> Vector a -> Vector b -> Vector c -> Vector d
VS.zipWith3
Int -> Double -> Double -> Double
getNewMassDiagonalWithRescue
Vector Int
sampleSizes
Vector Double
msDiagonalOld
Vector Double
msDiagonalEstimate
defaultGraphicalLassoPenalty :: Double
defaultGraphicalLassoPenalty :: Double
defaultGraphicalLassoPenalty = Double
0.4
tuneAllMasses ::
(a -> Positions) ->
VB.Vector a ->
(Masses, MassesI) ->
(Masses, MassesI)
tuneAllMasses :: (a -> Vector Double)
-> Vector a -> (Herm Double, GMatrix) -> (Herm Double, GMatrix)
tuneAllMasses a -> Vector Double
toVec Vector a
xs (Herm Double
ms, GMatrix
msI)
| Vector a -> Int
forall a. Vector a -> Int
VB.length Vector a
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
samplesDiagonalMin = (Herm Double
ms, GMatrix
msI)
| Vector a -> Int
forall a. Vector a -> Int
VB.length Vector a
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Int
samplesAllMinWith Int
dimMs = (Herm Double, GMatrix)
fallbackDiagonal
| Matrix Double -> Int
forall t. Field t => Matrix t -> Int
L.rank Matrix Double
xs'' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
dimState = (Herm Double, GMatrix)
fallbackDiagonal
| Int
dimState Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
dimMs = [Char] -> (Herm Double, GMatrix)
forall a. HasCallStack => [Char] -> a
error [Char]
"tuneAllMasses: Dimension mismatch."
| Bool
otherwise = (Herm Double
msPD', GMatrix
msPDI')
where
fallbackDiagonal :: (Herm Double, GMatrix)
fallbackDiagonal = (a -> Vector Double)
-> Vector a -> (Herm Double, GMatrix) -> (Herm Double, GMatrix)
forall a.
(a -> Vector Double)
-> Vector a -> (Herm Double, GMatrix) -> (Herm Double, GMatrix)
tuneDiagonalMassesOnly a -> Vector Double
toVec Vector a
xs (Herm Double
ms, GMatrix
msI)
xs' :: Vector (Vector Double)
xs' = (a -> Vector Double) -> Vector a -> Vector (Vector Double)
forall a b. (a -> b) -> Vector a -> Vector b
VB.map a -> Vector Double
toVec Vector a
xs
xs'' :: Matrix Double
xs'' = [Vector Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
L.fromRows ([Vector Double] -> Matrix Double)
-> [Vector Double] -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Vector (Vector Double) -> [Vector Double]
forall a. Vector a -> [a]
VB.toList Vector (Vector Double)
xs'
dimState :: Int
dimState = Vector Double -> Int
forall a. Storable a => Vector a -> Int
VS.length (Vector Double -> Int) -> Vector Double -> Int
forall a b. (a -> b) -> a -> b
$ Vector (Vector Double) -> Vector Double
forall a. Vector a -> a
VB.head Vector (Vector Double)
xs'
dimMs :: Int
dimMs = Matrix Double -> Int
forall t. Matrix t -> Int
L.rows (Matrix Double -> Int) -> Matrix Double -> Int
forall a b. (a -> b) -> a -> b
$ Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym Herm Double
ms
(Vector Double
_, Vector Double
ss, Matrix Double
xsNormalized) = Matrix Double -> (Vector Double, Vector Double, Matrix Double)
S.scale Matrix Double
xs''
sigmaNormalized :: Matrix Double
sigmaNormalized =
Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym (Herm Double -> Matrix Double) -> Herm Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$
([Char] -> Herm Double)
-> ((Herm Double, Herm Double) -> Herm Double)
-> Either [Char] (Herm Double, Herm Double)
-> Herm Double
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either [Char] -> Herm Double
forall a. HasCallStack => [Char] -> a
error (Herm Double, Herm Double) -> Herm Double
forall a b. (a, b) -> a
fst (Either [Char] (Herm Double, Herm Double) -> Herm Double)
-> Either [Char] (Herm Double, Herm Double) -> Herm Double
forall a b. (a -> b) -> a -> b
$
Double -> Matrix Double -> Either [Char] (Herm Double, Herm Double)
S.graphicalLasso Double
defaultGraphicalLassoPenalty Matrix Double
xsNormalized
sigma :: Matrix Double
sigma = Vector Double -> Matrix Double -> Matrix Double
S.rescaleSWith Vector Double
ss Matrix Double
sigmaNormalized
msI' :: GMatrix
msI' = Matrix Double -> GMatrix
toGMatrix Matrix Double
sigma
ms' :: Herm Double
ms' = Matrix Double -> Herm Double
forall t. Field t => Matrix t -> Herm t
L.sym (Matrix Double -> Herm Double) -> Matrix Double -> Herm Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
cleanMatrix (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
forall t. Field t => Matrix t -> Matrix t
L.inv Matrix Double
sigma
msPD' :: Herm Double
msPD' = Matrix Double -> Herm Double
forall t. Matrix t -> Herm t
L.trustSym (Matrix Double -> Herm Double) -> Matrix Double -> Herm Double
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
findClosestPositiveDefiniteMatrix (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym Herm Double
ms'
msPDI' :: GMatrix
msPDI' = if Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym Herm Double
ms' Matrix Double -> Matrix Double -> Bool
forall a. Eq a => a -> a -> Bool
== Herm Double -> Matrix Double
forall t. Herm t -> Matrix t
L.unSym Herm Double
msPD' then GMatrix
msI' else Herm Double -> GMatrix
getMassesI Herm Double
msPD'