ConClusion 0.2.1 → 0.2.2
raw patch · 8 files changed
+675/−648 lines, 8 filesdep ~aesondep ~basedep ~formattingPVP: major bump suggested
API removals or changes: PVP suggests a major version bump
Dependency ranges changed: aeson, base, formatting, text
API changes (from Hackage documentation)
- ConClusion.Numeric.Statistics: meanDeviation :: (Source r e, Fractional e, Unbox e, Numeric r e, Manifest r e) => Matrix r e -> Matrix r e
+ ConClusion.Numeric.Statistics: meanDeviation :: (Fractional e, Unbox e, Numeric r e, Manifest r e) => Matrix r e -> Matrix r e
Files
- Changelog.md +10/−0
- ConClusion.cabal +14/−10
- app/ConClusion.hs +5/−25
- src/ConClusion/Array/Conversion.hs +20/−19
- src/ConClusion/Array/Util.hs +48/−47
- src/ConClusion/BinaryTree.hs +36/−33
- src/ConClusion/Chemistry/Topology.hs +85/−81
- src/ConClusion/Numeric/Statistics.hs +457/−433
Changelog.md view
@@ -1,5 +1,15 @@ # Changelog +## 0.2.2+ - Dependency version bound updates+ - GHC 9.6 compatibility++## 0.2.1+ - Exchange some Mutable by Manifest restrictions for Massiv++## 0.2.0+ - Restructuring Modules+ ## 0.1.0 - updates to Massiv 1.0.0.0; changes lots of array types. - infrastructure updates
ConClusion.cabal view
@@ -1,18 +1,18 @@ cabal-version: 1.12 --- This file has been generated from package.yaml by hpack version 0.34.7.+-- This file has been generated from package.yaml by hpack version 0.35.2. -- -- see: https://github.com/sol/hpack name: ConClusion-version: 0.2.1+version: 0.2.2 synopsis: Cluster algorithms, PCA, and chemical conformere analysis description: Please see the README on GitLab at <https://gitlab.com/theoretical-chemistry-jena/quantum-chemistry/ConfoCluster> category: Statistics, Chemistry bug-reports: https://gitlab.com/theoretical-chemistry-jena/quantum-chemistry/ConfoCluster/-/issues author: Phillip Seeber maintainer: phillip.seeber@googlemail.com-copyright: 2022 Phillip Seeber+copyright: 2023 Phillip Seeber license: AGPL-3 license-file: LICENSE.md build-type: Simple@@ -60,13 +60,15 @@ TypeApplications RecordWildCards NamedFieldPuns+ TypeOperators+ LambdaCase ghc-options: -Wall -Wno-unused-top-binds -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wpartial-fields -Wredundant-constraints build-depends:- aeson >=1.5 && <2.2+ aeson >=1.5 && <2.3 , attoparsec >=0.13.0.0 && <0.15- , base >=4.7 && <4.18+ , base >=4.7 && <5 , containers >=0.6.0.0 && <0.7- , formatting >=7.1.0 && <7.2+ , formatting >=7.1.0 && <7.3 , hmatrix >=0.20.0 && <0.21 , massiv >=1.0.0.0 && <1.1 , psqueues >=0.2.7.0 && <0.3@@ -103,19 +105,21 @@ TypeApplications RecordWildCards NamedFieldPuns+ TypeOperators+ LambdaCase ghc-options: -Wall -Wno-unused-top-binds -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N build-depends: ConClusion- , aeson >=1.5 && <2.2+ , aeson >=1.5 && <2.3 , attoparsec >=0.13.0.0 && <0.15- , base >=4.7 && <4.18+ , base >=4.7 && <5 , cmdargs >=0.10.0 && <0.11 , containers >=0.6.0.0 && <0.7- , formatting >=7.1.0 && <7.2+ , formatting >=7.1.0 && <7.3 , hmatrix >=0.20.0 && <0.21 , massiv >=1.0.0.0 && <1.1 , optics >=0.3 && <0.5 , psqueues >=0.2.7.0 && <0.3 , rio >=0.1.13.0 && <0.2- , text >=1.2.0.0 && <2.1+ , text >=1.2.0.0 && <2.2 default-language: Haskell2010
app/ConClusion.hs view
@@ -1,3 +1,4 @@+{-# OPTIONS_GHC -Wno-partial-fields #-} -- | -- Module : Main -- Description : Command Line Interface to analyse CREST results.@@ -353,23 +354,7 @@ clustering :: Maybe Clustering, -- | RIO's 'logFunc'. logFunc :: LogFunc- }---- Lenses-instance (k ~ A_Lens, a ~ FilePath, b ~ a) => LabelOptic "xyz" k App App a b where- labelOptic = lens (xyz :: App -> FilePath) $ \s b -> (s {xyz = b} :: App)--instance (k ~ A_Lens, a ~ [Feature], b ~ a) => LabelOptic "dim" k App App a b where- labelOptic = lens (dim :: App -> [Feature]) $ \s b -> (s {dim = b} :: App)--instance (k ~ A_Lens, a ~ Maybe PrincipalComponentAnalysis, b ~ a) => LabelOptic "pca" k App App a b where- labelOptic = lens (pca :: App -> Maybe PrincipalComponentAnalysis) $ \s b -> (s {pca = b} :: App)--instance (k ~ A_Lens, a ~ Maybe Clustering, b ~ a) => LabelOptic "clustering" k App App a b where- labelOptic = lens clustering $ \s b -> s {clustering = b}--instance (k ~ A_Lens, a ~ LogFunc, b ~ a) => LabelOptic "logFunc" k App App a b where- labelOptic = lens logFunc $ \s b -> s {logFunc = b}+ } deriving (Generic) -- Reader Classes instance HasXYZ App where@@ -391,11 +376,7 @@ newtype PrincipalComponentAnalysis = PrincipalComponentAnalysis { -- | Number of dimensions to keep from PCA. keep :: Int- }---- Lenses-instance (k ~ A_Lens, a ~ Int, b ~ a) => LabelOptic "keep" k PrincipalComponentAnalysis PrincipalComponentAnalysis a b where- labelOptic = lens keep $ \s b -> s {keep = b}+ } deriving (Generic) -- | Settings for clustering. data Clustering@@ -498,8 +479,7 @@ -- Log function setup. logOptions <-- id- . setLogUseTime False+ setLogUseTime False . setLogUseColor True . setLogVerboseFormat False . setLogUseLoc False@@ -523,7 +503,7 @@ _ <- option ',' (char ',') return fType where- eParser = char 'e' *> return Energy+ eParser = char 'e' $> Energy bParser = do _ <- char 'b' skipMany (char ' ')
src/ConClusion/Array/Conversion.hs view
@@ -1,17 +1,18 @@--- |--- Module : ConClusion.Array.Conversion--- Description : Type castings and conversions of array types--- Copyright : Phillip Seeber, 2022--- License : AGPL-3--- Maintainer : phillip.seeber@googlemail.com--- Stability : experimental--- Portability : POSIX, Windows-module ConClusion.Array.Conversion- ( vecH2M,- vecM2H,- matH2M,- matM2H,- )+{- |+Module : ConClusion.Array.Conversion+Description : Type castings and conversions of array types+Copyright : Phillip Seeber, 2023+License : AGPL-3+Maintainer : phillip.seeber@googlemail.com+Stability : experimental+Portability : POSIX, Windows+-}+module ConClusion.Array.Conversion (+ vecH2M,+ vecM2H,+ matH2M,+ matM2H,+) where import Data.Massiv.Array as Massiv hiding (IndexException)@@ -34,13 +35,13 @@ {-# SCC matH2M #-} matH2M :: (Manifest r e, Load r Ix1 e, Element e) => LA.Matrix e -> Massiv.Matrix r e matH2M hMat = Massiv.resize' (Sz $ nRows :. nCols) . vecH2M . LA.flatten $ hMat- where- nRows = LA.rows hMat- nCols = LA.cols hMat+ where+ nRows = LA.rows hMat+ nCols = LA.cols hMat -- | Converts a matrix from Massiv to HMatrix representation. {-# SCC matM2H #-} matM2H :: (Element e, Manifest r e, Load r Ix1 e) => Massiv.Matrix r e -> LA.Matrix e matM2H mMat = LA.reshape nCols . vecM2H . Massiv.flatten $ mMat- where- Sz (_nRows :. nCols) = Massiv.size mMat+ where+ Sz (_nRows :. nCols) = Massiv.size mMat
src/ConClusion/Array/Util.hs view
@@ -1,20 +1,21 @@--- |--- Module : ConClusion.Array.Util--- Description : Additional tools to work with numerical arrays--- Copyright : Phillip Seeber, 2022--- License : AGPL-3--- Maintainer : phillip.seeber@googlemail.com--- Stability : experimental--- Portability : POSIX, Windows-module ConClusion.Array.Util- ( IndexException (..),- magnitude,- normalise,- angle,- minDistAt,- minDistAtVec,- iMinimumM,- )+{- |+Module : ConClusion.Array.Util+Description : Additional tools to work with numerical arrays+Copyright : Phillip Seeber, 2023+License : AGPL-3+Maintainer : phillip.seeber@googlemail.com+Stability : experimental+Portability : POSIX, Windows+-}+module ConClusion.Array.Util (+ IndexException (..),+ magnitude,+ normalise,+ angle,+ minDistAt,+ minDistAtVec,+ iMinimumM,+) where import Data.Massiv.Array as Massiv hiding (IndexException)@@ -41,27 +42,27 @@ -- | Find the minimal distance in a distance matrix, which is not the main diagonal. {-# SCC minDistAt #-} minDistAt ::- ( Manifest r e,- MonadThrow m,- Ord e+ ( Manifest r e+ , MonadThrow m+ , Ord e ) => Massiv.Matrix r e -> m (e, Ix2) minDistAt arr | isEmpty arr = throwM $ SizeEmptyException (Massiv.size arr) | otherwise = return . unsafePerformIO $ ifoldlP minFold start chFold start arr- where- ix0 = pureIndex 0- e0 = arr Massiv.! ix0- start = (e0, ix0)- minFold acc@(eA, _) ix@(m :. n) e = if e < eA && m > n then (e, ix) else acc- chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc+ where+ ix0 = pureIndex 0+ e0 = arr Massiv.! ix0+ start = (e0, ix0)+ minFold acc@(eA, _) ix@(m :. n) e = if e < eA && m > n then (e, ix) else acc+ chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc -- | Find the minimal element of a vector, which is at a larger than the supplied index. minDistAtVec ::- ( Manifest r e,- MonadThrow m,- Ord e+ ( Manifest r e+ , MonadThrow m+ , Ord e ) => Ix1 -> Massiv.Vector r e ->@@ -72,31 +73,31 @@ | otherwise = do let (minE, minIx) = unsafePerformIO $ ifoldlP minFold startAcc chFold startAcc searchVec return (minE, minIx + ixStart + 1)- where- Sz nElems = Massiv.size vec- searchVec = Massiv.drop (Sz $ ixStart + 1) vec- ix0 = 0- e0 = searchVec Massiv.! ix0- startAcc = (e0, ix0)- minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc- chFold acc ch = min acc ch+ where+ Sz nElems = Massiv.size vec+ searchVec = Massiv.drop (Sz $ ixStart + 1) vec+ ix0 = 0+ e0 = searchVec Massiv.! ix0+ startAcc = (e0, ix0)+ minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc+ chFold acc ch = min acc ch -- | Like 'Massiv.minimumM' but also returns the index of the minimal element. iMinimumM ::- ( Manifest r a,- MonadThrow m,- Index ix,- Ord a+ ( Manifest r a+ , MonadThrow m+ , Index ix+ , Ord a ) => Array r ix a -> m (a, ix) iMinimumM arr | isEmpty arr = throwM $ SizeEmptyException (Massiv.size arr) | otherwise = return . unsafePerformIO $ ifoldlP minFold start chFold start arr- where- ix0 = pureIndex 0- e0 = arr Massiv.! ix0- start = (e0, ix0)+ where+ ix0 = pureIndex 0+ e0 = arr Massiv.! ix0+ start = (e0, ix0) - minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc- chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc+ minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc+ chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc
src/ConClusion/BinaryTree.hs view
@@ -1,17 +1,18 @@--- |--- Module : ConClusion.BinaryTree--- Description : Custom binary tree type with some special functions--- Copyright : Phillip Seeber, 2022--- License : AGPL-3--- Maintainer : phillip.seeber@googlemail.com--- Stability : experimental--- Portability : POSIX, Windows-module ConClusion.BinaryTree- ( BinTree (..),- root,- takeBranchesWhile,- takeLeafyBranchesWhile,- )+{- |+Module : ConClusion.BinaryTree+Description : Custom binary tree type with some special functions+Copyright : Phillip Seeber, 2023+License : AGPL-3+Maintainer : phillip.seeber@googlemail.com+Stability : experimental+Portability : POSIX, Windows+-}+module ConClusion.BinaryTree (+ BinTree (..),+ root,+ takeBranchesWhile,+ takeLeafyBranchesWhile,+) where import Data.Aeson hiding (Array)@@ -35,27 +36,29 @@ root (Leaf e) = e root (Node e _ _) = e --- | Steps down each branch of a tree until some criterion is satisfied or the --- end of the branch is reached. Each end of the branch is added to a result.+{- | Steps down each branch of a tree until some criterion is satisfied or the+end of the branch is reached. Each end of the branch is added to a result.+-} takeBranchesWhile :: (a -> Bool) -> BinTree a -> Massiv.Vector DL a takeBranchesWhile chk tree = go tree (Massiv.empty @DL)- where- go (Leaf v) acc = if chk v then acc `snoc` v else acc- go (Node v l r) acc =- let vAcc = if chk v then acc `snoc` v else acc- lAcc = go l vAcc- rAcc = go r lAcc- in if chk v then rAcc else vAcc+ where+ go (Leaf v) acc = if chk v then acc `snoc` v else acc+ go (Node v l r) acc =+ let vAcc = if chk v then acc `snoc` v else acc+ lAcc = go l vAcc+ rAcc = go r lAcc+ in if chk v then rAcc else vAcc --- | Takes the first value in each branch, that does not fullfill the criterion --- anymore and adds it to the result. Terminal leafes of the branches are always--- taken.+{- | Takes the first value in each branch, that does not fullfill the criterion+anymore and adds it to the result. Terminal leafes of the branches are always+taken.+-} takeLeafyBranchesWhile :: (a -> Bool) -> BinTree a -> Massiv.Vector DL a takeLeafyBranchesWhile chk tree = go tree (Massiv.empty @DL)- where- go (Leaf v) acc = acc `snoc` v- go (Node v l r) acc =- let vAcc = if chk v then acc else acc `snoc` v- lAcc = go l vAcc- rAcc = go r lAcc- in if chk v then rAcc else vAcc+ where+ go (Leaf v) acc = acc `snoc` v+ go (Node v l r) acc =+ let vAcc = if chk v then acc else acc `snoc` v+ lAcc = go l vAcc+ rAcc = go r lAcc+ in if chk v then rAcc else vAcc
src/ConClusion/Chemistry/Topology.hs view
@@ -1,29 +1,30 @@--- |--- Module : ConClusion.Chemistry.Topology--- Description : Principal Component Analysis--- Copyright : Phillip Seeber, 2021--- License : AGPL-3--- Maintainer : phillip.seeber@googlemail.com--- Stability : experimental--- Portability : POSIX, Windows------ This module implements routines to work with simple molden style XYZ trajectories. This includes--- parsers as well as functions to obtain structural features in internal coordinates.------ For an introduction into PCA see <https://www.cs.cmu.edu/~elaw/papers/pca.pdf>.------ Diherdrals require a special metric, see <https://onlinelibrary.wiley.com/doi/full/10.1002/prot.20310)>.-module ConClusion.Chemistry.Topology- ( Molecule,- Trajectory,- xyz,- trajectory,- B (..),- A (..),- D (..),- Feature (..),- getFeatures,- )+{- |+Module : ConClusion.Chemistry.Topology+Description : Principal Component Analysis+Copyright : Phillip Seeber, 2023+License : AGPL-3+Maintainer : phillip.seeber@googlemail.com+Stability : experimental+Portability : POSIX, Windows++This module implements routines to work with simple molden style XYZ trajectories. This includes+parsers as well as functions to obtain structural features in internal coordinates.++For an introduction into PCA see <https://www.cs.cmu.edu/~elaw/papers/pca.pdf>.++Diherdrals require a special metric, see <https://onlinelibrary.wiley.com/doi/full/10.1002/prot.20310)>.+-}+module ConClusion.Chemistry.Topology (+ Molecule,+ Trajectory,+ xyz,+ trajectory,+ B (..),+ A (..),+ D (..),+ Feature (..),+ getFeatures,+) where import ConClusion.Array.Conversion@@ -40,12 +41,12 @@ -- | A Molecule in cartesian coordinates. data Molecule = Molecule- { -- | The energy of the molecule.- energy :: !Double,- -- | Chemical symbols or names of the atoms. \(N\) vector.- atoms :: !(VectorB.Vector Text),- -- | Cartesian coordinates. Atoms as rows, xyz as columns. \(N \times 3\) matrix.- coordinates :: !(Massiv.Matrix S Double)+ { energy :: !Double+ -- ^ The energy of the molecule.+ , atoms :: !(VectorB.Vector Text)+ -- ^ Chemical symbols or names of the atoms. \(N\) vector.+ , coordinates :: !(Massiv.Matrix S Double)+ -- ^ Cartesian coordinates. Atoms as rows, xyz as columns. \(N \times 3\) matrix. } type Trajectory = Seq Molecule@@ -86,9 +87,9 @@ return Molecule- { energy = energy,- atoms = VectorB.fromListN nAtoms . fmap fst $ atoms,- coordinates = Massiv.compute . Massiv.concat' (Dim 2) . fmap snd $ atoms+ { energy = energy+ , atoms = VectorB.fromListN nAtoms . fmap fst $ atoms+ , coordinates = Massiv.compute . Massiv.concat' (Dim 2) . fmap snd $ atoms } -- | Parser for trajectories in XYZ format as produced by CREST.@@ -117,61 +118,63 @@ | Dihedral D -- | Calculate a bond distance.-bond :: MonadThrow m => B -> Molecule -> m Double-bond (B a b) Molecule {coordinates}+bond :: (MonadThrow m) => B -> Molecule -> m Double+bond (B a b) Molecule{coordinates} | a == b = throwM . IndexException $ "selected atoms are identical" | otherwise = do- coordA <- compute @U <$> (coordinates !?> a)- coordB <- compute @U <$> (coordinates !?> b)- vecAB <- coordA .-. coordB- return . sqrt . Massiv.sum . Massiv.map (^ (2 :: Int)) $ vecAB+ coordA <- compute @U <$> (coordinates !?> a)+ coordB <- compute @U <$> (coordinates !?> b)+ vecAB <- coordA .-. coordB+ return . sqrt . Massiv.sum . Massiv.map (^ (2 :: Int)) $ vecAB -- | Calculates the sinus of an angle defined by three atoms.-angle :: MonadThrow m => A -> Molecule -> m Double-angle (A a b c) Molecule {coordinates}+angle :: (MonadThrow m) => A -> Molecule -> m Double+angle (A a b c) Molecule{coordinates} | a == b || b == c || a == c = throwM . IndexException $ "selected atoms are identical" | otherwise = do- coordA <- compute @U <$> (coordinates !?> a)- coordB <- compute @U <$> (coordinates !?> b)- coordC <- compute @U <$> (coordinates !?> c)- vecAB <- coordA .-. coordB- vecCB <- coordC .-. coordB- return $ ArrayUtil.angle vecAB vecCB+ coordA <- compute @U <$> (coordinates !?> a)+ coordB <- compute @U <$> (coordinates !?> b)+ coordC <- compute @U <$> (coordinates !?> c)+ vecAB <- coordA .-. coordB+ vecCB <- coordC .-. coordB+ return $ ArrayUtil.angle vecAB vecCB --- | Calculates the dihedral angle defined by four atoms. Respects rotation direction. Obtains the--- result in radian.-dihedral' :: MonadThrow m => D -> Molecule -> m Double-dihedral' (D a b c d) Molecule {coordinates}+{- | Calculates the dihedral angle defined by four atoms. Respects rotation direction. Obtains the+result in radian.+-}+dihedral' :: (MonadThrow m) => D -> Molecule -> m Double+dihedral' (D a b c d) Molecule{coordinates} | a == b || a == c || a == d || b == c || b == d || c == d = throwM . IndexException $ "selected atoms are indentical" | otherwise = do- coordA <- compute @U <$> (coordinates !?> a)- coordB <- compute @U <$> (coordinates !?> b)- coordC <- compute @U <$> (coordinates !?> c)- coordD <- compute @U <$> (coordinates !?> d)- vecAB <- coordA .-. coordB- vecBC <- coordB .-. coordC- vecCD <- coordC .-. coordD- let planeABC = vecH2M $ LA.cross (vecM2H vecAB) (vecM2H vecBC) :: Massiv.Vector U Double- planeBCD = vecH2M $ LA.cross (vecM2H vecBC) (vecM2H vecCD) :: Massiv.Vector U Double- normVecRot = LA.cross (vecM2H vecCD) (vecM2H vecBC) :: LA.Vector Double- rotDir =- if vecH2M normVecRot !.! vecAB < 0- then -1- else 1- return $ rotDir * ArrayUtil.angle planeABC planeBCD+ coordA <- compute @U <$> (coordinates !?> a)+ coordB <- compute @U <$> (coordinates !?> b)+ coordC <- compute @U <$> (coordinates !?> c)+ coordD <- compute @U <$> (coordinates !?> d)+ vecAB <- coordA .-. coordB+ vecBC <- coordB .-. coordC+ vecCD <- coordC .-. coordD+ let planeABC = vecH2M $ LA.cross (vecM2H vecAB) (vecM2H vecBC) :: Massiv.Vector U Double+ planeBCD = vecH2M $ LA.cross (vecM2H vecBC) (vecM2H vecCD) :: Massiv.Vector U Double+ normVecRot = LA.cross (vecM2H vecCD) (vecM2H vecBC) :: LA.Vector Double+ rotDir =+ if vecH2M normVecRot !.! vecAB < 0+ then -1+ else 1+ return $ rotDir * ArrayUtil.angle planeABC planeBCD --- | Calculates a metric value of the dihedral angle defined by four atoms. This must create 2--- values in the feature matrix, instead of one.--- See <https://onlinelibrary.wiley.com/doi/full/10.1002/prot.20310)>-dihedral :: MonadThrow m => D -> Molecule -> m (Double, Double)+{- | Calculates a metric value of the dihedral angle defined by four atoms. This must create 2+values in the feature matrix, instead of one.+See <https://onlinelibrary.wiley.com/doi/full/10.1002/prot.20310)>+-}+dihedral :: (MonadThrow m) => D -> Molecule -> m (Double, Double) dihedral d mol = do dihedRad <- dihedral' d mol return (sin dihedRad, cos dihedRad) -- | Get all selected features from a molecule. getFeature ::- ( Traversable t,- MonadThrow m+ ( Traversable t+ , MonadThrow m ) => -- | Selection of multiple features. t Feature ->@@ -179,7 +182,7 @@ Molecule -> m (Massiv.Vector Massiv.DL Double) getFeature sel mol = do- features <- for sel $ \s -> case s of+ features <- for sel $ \case Energy -> return . Left . energy $ mol Bond b -> Left <$> bond b mol Angle a -> Left <$> angle a mol@@ -196,16 +199,17 @@ return featureVec --- | Obtains the feature matrix \(\mathbf{X}\) for a principal component analysis. Given \(m\)--- features to analyse in \(n\) measurements, \(\mathbf{X}\) will be a \(m \times n\) matrix.+{- | Obtains the feature matrix \(\mathbf{X}\) for a principal component analysis. Given \(m\)+features to analyse in \(n\) measurements, \(\mathbf{X}\) will be a \(m \times n\) matrix.+-} {-# SCC getFeatures #-} getFeatures ::- ( MonadThrow m,- Traversable f+ ( MonadThrow m+ , Traversable f ) => f Feature -> Trajectory -> m (Massiv.Matrix DL Double) getFeatures sel trj = traverse toCols trj >>= concatM (Dim 1)- where- toCols v = expandInner @U @Ix2 (Sz 1) const . compute @U <$> getFeature sel v+ where+ toCols v = expandInner @U @Ix2 (Sz 1) const . compute @U <$> getFeature sel v
src/ConClusion/Numeric/Statistics.hs view
@@ -1,41 +1,42 @@--- |--- Module : ConClusion.Numeric.Statistics--- Description : Statistical Functions--- Copyright : Phillip Seeber, 2021--- License : AGPL-3--- Maintainer : phillip.seeber@googlemail.com--- Stability : experimental--- Portability : POSIX, Windows-module ConClusion.Numeric.Statistics- ( -- * PCA- PCA (..),- pca,+{- |+Module : ConClusion.Numeric.Statistics+Description : Statistical Functions+Copyright : Phillip Seeber, 2023+License : AGPL-3+Maintainer : phillip.seeber@googlemail.com+Stability : experimental+Portability : POSIX, Windows+-}+module ConClusion.Numeric.Statistics (+ -- * PCA+ PCA (..),+ pca, - -- * Variance- normalise,- meanDeviation,- covariance,+ -- * Variance+ normalise,+ meanDeviation,+ covariance, - -- * Distance Metrics- DistFn,- lpNorm,- manhattan,- euclidean,- mahalanobis,+ -- * Distance Metrics+ DistFn,+ lpNorm,+ manhattan,+ euclidean,+ mahalanobis, - -- * Cluster Algorithms- Clusters,+ -- * Cluster Algorithms+ Clusters, - -- ** DBScan- DistanceInvalidException (..),- dbscan,+ -- ** DBScan+ DistanceInvalidException (..),+ dbscan, - -- ** Hierarchical Cluster Analysis- Dendrogram,- JoinStrat (..),- hca,- cutDendroAt,- )+ -- ** Hierarchical Cluster Analysis+ Dendrogram,+ JoinStrat (..),+ hca,+ cutDendroAt,+) where import ConClusion.Array.Conversion@@ -47,6 +48,7 @@ import qualified Data.IntSet as IntSet import Data.Massiv.Array as Massiv import Data.Massiv.Array.Unsafe as Massiv+import Data.Type.Equality import qualified Numeric.LinearAlgebra as LA import RIO hiding (Vector) import System.IO.Unsafe (unsafePerformIO)@@ -59,36 +61,35 @@ eig :: ( -- Manifest r1 Ix1 (Complex Double), -- Manifest r2 Ix1 (Complex Double),- LA.Field e,- Manifest r3 e,- Manifest r1 (Complex Double),- Manifest r2 (Complex Double),- Load r1 Ix1 (Complex Double),- Load r2 Ix1 (Complex Double),- Load r3 Ix1 e,- MonadThrow m+ LA.Field e+ , Manifest r3 e+ , Manifest r1 (Complex Double)+ , Manifest r2 (Complex Double)+ , Load r1 Ix1 (Complex Double)+ , Load r2 Ix1 (Complex Double)+ , Load r3 Ix1 e+ , MonadThrow m ) => Matrix r3 e -> m (Vector r1 (Complex Double), Matrix r2 (Complex Double)) eig covM | m /= n = throwM $ IndexException "eigenvalue problems can only be solved for square matrix" | otherwise = return . bimap vecH2M matH2M . LA.eig $ cov- where- Sz (m :. n) = size covM- cov = matM2H covM+ where+ Sz (m :. n) = size covM+ cov = matM2H covM --- | Sort eigenvalues and eigenvectors by magnitude of the eigenvalues in descending order (largest--- eigenvalues first). Eigenvectors are the columns of the input matrix.+{- | Sort eigenvalues and eigenvectors by magnitude of the eigenvalues in descending order (largest+eigenvalues first). Eigenvectors are the columns of the input matrix.+-} {-# SCC eigSort #-} eigSort ::- ( Load r2 Ix2 e,- MonadThrow m,- Source r1 e,- Source r2 e,- Manifest r1 e,- Manifest r2 e,- Unbox e,- Ord e+ ( Load r2 Ix2 e+ , MonadThrow m+ , Manifest r1 e+ , Manifest r2 e+ , Unbox e+ , Ord e ) => (Vector r1 e, Matrix r2 e) -> m (Vector r1 e, Matrix r2 e)@@ -96,61 +97,62 @@ | m /= n = throwM $ IndexException "matrix of the eigenvectors is not a square matrix" | n /= n' = throwM $ IndexException "different number of eigenvalues and eigenvectors" | otherwise = do- let ixedEigenvalues = Massiv.zip vec ixVec- (eigenValSortAsc, ixSort) = (\a -> (get fst a, get snd a)) . quicksort . compute @U $ ixedEigenvalues- eigenVecSortAsc = backpermute' (Sz $ m :. n) (\(r :. c) -> r :. (ixSort ! c)) mat- eigenValSort = reverse' (Dim 1) eigenValSortAsc- eigenVecSort = reverse' (Dim 1) eigenVecSortAsc- return (compute eigenValSort, compute eigenVecSort)- where- Sz (m :. n) = size mat- Sz n' = size vec- ixVec = makeArrayLinear @D Seq (Sz n') id- get acc = compute @U . Massiv.map acc+ let ixedEigenvalues = Massiv.zip vec ixVec+ (eigenValSortAsc, ixSort) = (\a -> (get fst a, get snd a)) . quicksort . compute @U $ ixedEigenvalues+ eigenVecSortAsc = backpermute' (Sz $ m :. n) (\(r :. c) -> r :. (ixSort ! c)) mat+ eigenValSort = reverse' (Dim 1) eigenValSortAsc+ eigenVecSort = reverse' (Dim 1) eigenVecSortAsc+ return (compute eigenValSort, compute eigenVecSort)+ where+ Sz (m :. n) = size mat+ Sz n' = size vec+ ixVec = makeArrayLinear @D Seq (Sz n') id+ get acc = compute @U . Massiv.map acc -- | Adjust function for priority queues. Updates the priority at a given key if present. pqAdjust :: (Ord k, Hashable k, Ord p) => (p -> p) -> k -> PQ.HashPSQ k p v -> PQ.HashPSQ k p v pqAdjust f k q = snd $ PQ.alter f' k q- where- f' = \op -> case op of- Nothing -> (False, Nothing)- Just (p, v) -> (False, Just (f p, v))+ where+ f' = \case+ Nothing -> (False, Nothing)+ Just (p, v) -> (False, Just (f p, v)) ---------------------------------------------------------------------------------------------------- -- Principal Component Analysis data PCA = PCA- { -- | Original feature matrix.- x :: Matrix U Double,- -- | Feature matrix in mean deviation form.- x' :: Matrix U Double,- -- | Transformed data.- y :: Matrix U Double,- -- | Transformation matrix to transform feature matrix into PCA result matrix.- a :: Matrix U Double,- -- | Mean squared error introduced by PCA.- mse :: Double,- -- | Percentage of the behaviour captured in the remaining dimensions.- remaining :: Double,- -- | All eigenvalues from the diagonalisation of the covariance matrix.- allEigenValues :: Vector U Double,- -- | Eigenvalues that were kept for PCA.- pcaEigenValues :: Vector U Double,- -- | All eigenvectors from the diagonalisation of the covariance matrix.- allEigenVecs :: Matrix U Double,- -- | Eigenvectors that were kept for PCA.- pcaEigenVecs :: Matrix U Double+ { x :: Matrix U Double+ -- ^ Original feature matrix.+ , x' :: Matrix U Double+ -- ^ Feature matrix in mean deviation form.+ , y :: Matrix U Double+ -- ^ Transformed data.+ , a :: Matrix U Double+ -- ^ Transformation matrix to transform feature matrix into PCA result matrix.+ , mse :: Double+ -- ^ Mean squared error introduced by PCA.+ , remaining :: Double+ -- ^ Percentage of the behaviour captured in the remaining dimensions.+ , allEigenValues :: Vector U Double+ -- ^ All eigenvalues from the diagonalisation of the covariance matrix.+ , pcaEigenValues :: Vector U Double+ -- ^ Eigenvalues that were kept for PCA.+ , allEigenVecs :: Matrix U Double+ -- ^ All eigenvectors from the diagonalisation of the covariance matrix.+ , pcaEigenVecs :: Matrix U Double+ -- ^ Eigenvectors that were kept for PCA. } --- | Transform the input values with a transformation matrix \(\mathbf{A}\), where \(\mathbf{A}\) is--- constructed from the eigenvectors associated to the largest eigenvalues.+{- | Transform the input values with a transformation matrix \(\mathbf{A}\), where \(\mathbf{A}\) is+constructed from the eigenvectors associated to the largest eigenvalues.+-} {-# SCC transformToPCABasis #-} transformToPCABasis :: ( -- Source (R r) Ix2 e, -- Extract r Ix2 e,- Manifest r e,- Numeric r e,- MonadThrow m+ Manifest r e+ , Numeric r e+ , MonadThrow m ) => -- | Number of dimensions to keep from PCA. Int ->@@ -168,24 +170,25 @@ | nDim >= nE = throwM $ IndexException "more than the possible amount of dimensions has been selected" | mE /= mF = throwM $ IndexException "eigenvector matrix and feature matrix have mismatching dimensions" | otherwise = do- matA <- compute . transpose <$> extractM (0 :. 0) (Sz $ mE :. nDim) eigenVecMat- pcaData <- matA .><. featureMat- return (pcaData, matA)- where- Sz (mE :. nE) = size eigenVecMat- Sz (mF :. _nF) = size featureMat+ matA <- compute . transpose <$> extractM (0 :. 0) (Sz $ mE :. nDim) eigenVecMat+ pcaData <- matA .><. featureMat+ return (pcaData, matA)+ where+ Sz (mE :. nE) = size eigenVecMat+ Sz (mF :. _nF) = size featureMat --- | Performs a PCA on the feature matrix \(\mathbf{X}\) by solving the eigenproblem of the--- covariance matrix. The function takes the feature matrix directly and perfoms the conversion--- to mean deviation form, the calculation of the covariance matrix and the eigenvalue problem--- automatically.+{- | Performs a PCA on the feature matrix \(\mathbf{X}\) by solving the eigenproblem of the+covariance matrix. The function takes the feature matrix directly and perfoms the conversion+to mean deviation form, the calculation of the covariance matrix and the eigenvalue problem+automatically.+-} {-# SCC pca #-} pca ::- ( Numeric r Double,- Manifest r Double,- Load r Ix1 Double,- Load r Ix2 Double,- MonadThrow m+ ( Numeric r Double+ , Manifest r Double+ , Load r Ix1 Double+ , Load r Ix2 Double+ , MonadThrow m ) => -- | Dimensionalty after PCA transformation. Int ->@@ -221,61 +224,62 @@ return PCA- { x = compute x,- x' = compute x',- y = compute pcaData,- a = compute matA,- mse = mse,- remaining = remaining,- allEigenValues = eValS,- pcaEigenValues = compute pcaEigenValues,- allEigenVecs = compute eVecS,- pcaEigenVecs = compute pcaEigenVecs+ { x = compute x+ , x' = compute x'+ , y = compute pcaData+ , a = compute matA+ , mse = mse+ , remaining = remaining+ , allEigenValues = eValS+ , pcaEigenValues = compute pcaEigenValues+ , allEigenVecs = compute eVecS+ , pcaEigenVecs = compute pcaEigenVecs }- where- Sz (m :. n) = size x+ where+ Sz (m :. n) = size x ---------------------------------------------------------------------------------------------------- -- Variance --- | Subtract the mean value of all columns from the feature matrix. Brings the feature matrix to--- mean deviation form.+{- | Subtract the mean value of all columns from the feature matrix. Brings the feature matrix to+mean deviation form.+-} {-# SCC meanDeviation #-} meanDeviation ::- ( Source r e,- Fractional e,- Unbox e,- Numeric r e,- Manifest r e+ ( Fractional e+ , Unbox e+ , Numeric r e+ , Manifest r e ) => Matrix r e -> Matrix r e meanDeviation mat = mat !-! compute meanMat- where- Sz (_ :. n) = Massiv.size mat- featueMean = Massiv.foldlInner (+) 0 mat .* (1 / fromIntegral n)- meanMat = expandInner (Sz n) const . compute @U $ featueMean+ where+ Sz (_ :. n) = Massiv.size mat+ featueMean = Massiv.foldlInner (+) 0 mat .* (1 / fromIntegral n)+ meanMat = expandInner (Sz n) const . compute @U $ featueMean --- | Obtains the covariance matrix \(\mathbf{C_X}\) from the feature matrix \(\mathbf{X}\).--- \[--- \mathbf{C_X} \equiv \frac{1}{n - 1} \mathbf{X} \mathbf{X}^T--- \]--- where \(n\) is the number of columns in the matrix.------ The feature matrix should be in mean deviation form, see 'meanDeviation'.+{- | Obtains the covariance matrix \(\mathbf{C_X}\) from the feature matrix \(\mathbf{X}\).+\[+ \mathbf{C_X} \equiv \frac{1}{n - 1} \mathbf{X} \mathbf{X}^T+\]+where \(n\) is the number of columns in the matrix.++The feature matrix should be in mean deviation form, see 'meanDeviation'.+-} {-# SCC covariance #-} covariance :: (Numeric r e, Manifest r e, Fractional e) => Matrix r e -> Matrix r e covariance x = (1 / (fromIntegral n - 1)) *. (x !><! (compute . transpose $ x))- where- Sz (_ :. n) = size x+ where+ Sz (_ :. n) = size x -- | Normalise each value so that the maximum absolute value in each row becomes one. normalise ::- ( Ord e,- Unbox e,- Numeric r e,- Fractional e,- Manifest r e+ ( Ord e+ , Unbox e+ , Numeric r e+ , Fractional e+ , Manifest r e ) => Array r Ix2 e -> Array r Ix2 e@@ -284,8 +288,8 @@ maxPerRow = compute @U . foldlInner max 0 $ absMat divMat = compute . Massiv.map (1 /) . expandInner @U @Ix2 (Sz n) const $ maxPerRow in divMat !*! mat- where- Sz (_ :. n) = size mat+ where+ Sz (_ :. n) = size mat ---------------------------------------------------------------------------------------------------- -- Distance Measures@@ -313,41 +317,45 @@ ab = Massiv.zipWith zipFn a b d = foldlInner foldFn acc ab in d- where- Sz (_ :. n) = size mat+ where+ Sz (_ :. n) = size mat --- | The \(L_p\) norm between two vectors. Generalisation of Manhattan and Euclidean distances.--- \[--- d(\mathbf{a}, \mathbf{b}) = \left( \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert ^p \right) ^ \frac{1}{p}--- \]+{- | The \(L_p\) norm between two vectors. Generalisation of Manhattan and Euclidean distances.+\[+ d(\mathbf{a}, \mathbf{b}) = \left( \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert ^p \right) ^ \frac{1}{p}+\]+-} {-# SCC lpNorm #-} lpNorm :: (Manifest r e, Floating e) => Int -> DistFn r e lpNorm p = compute . buildDistMat zipFn foldFn 0- where- zipFn a b = (^ p) . abs $ a - b- foldFn a b = (** (1 / fromIntegral p)) $ a + b+ where+ zipFn a b = (^ p) . abs $ a - b+ foldFn a b = (** (1 / fromIntegral p)) $ a + b --- | The Manhattan distance between two vectors. Specialisation of the \(L_p\) norm for \(p = 1\).--- \[--- d(\mathbf{a}, \mathbf{b}) = \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert--- \]+{- | The Manhattan distance between two vectors. Specialisation of the \(L_p\) norm for \(p = 1\).+\[+ d(\mathbf{a}, \mathbf{b}) = \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert+\]+-} {-# SCC manhattan #-} manhattan :: (Manifest r e, Floating e) => DistFn r e manhattan = lpNorm 1 --- | The Euclidean distance between two vectors. Specialisation of the \(L_p\) norm for \(p = 2\).--- \[--- d(\mathbf{a}, \mathbf{b}) = \sqrt{\sum \limits_{i=1}^n (\mathbf{a}_i - \mathbf{b}_i)^2}--- \]+{- | The Euclidean distance between two vectors. Specialisation of the \(L_p\) norm for \(p = 2\).+\[+ d(\mathbf{a}, \mathbf{b}) = \sqrt{\sum \limits_{i=1}^n (\mathbf{a}_i - \mathbf{b}_i)^2}+\]+-} {-# SCC euclidean #-} euclidean :: (Manifest r e, Floating e) => DistFn r e euclidean = lpNorm 2 --- | Mahalanobis distance between points. Suitable for non correlated axes.--- \[--- d(\mathbf{a}, \mathbf{b}) = \sqrt{(\mathbf{a} - \mathbf{b})^T \mathbf{S}^{-1} (\mathbf{a} - \mathbf{b})}--- \]--- where \(\mathbf{S}\) is the covariance matrix.+{- | Mahalanobis distance between points. Suitable for non correlated axes.+\[+ d(\mathbf{a}, \mathbf{b}) = \sqrt{(\mathbf{a} - \mathbf{b})^T \mathbf{S}^{\-1} (\mathbf{a} - \mathbf{b})}+\]+where \(\mathbf{S}\) is the covariance matrix.+-} {-# SCC mahalanobis #-} mahalanobis :: (Unbox e, Numeric r e, Manifest r e, LA.Field e, Load r Ix1 e) => DistFn r e mahalanobis points =@@ -363,10 +371,10 @@ ) ixArray in compute . Massiv.map sqrt $ distMat- where- Sz (_ :. n) = size points- cov = covariance . meanDeviation $ points- covInv = matH2M . LA.inv . matM2H $ cov+ where+ Sz (_ :. n) = size points+ cov = covariance . meanDeviation $ points+ covInv = matH2M . LA.inv . matM2H $ cov ---------------------------------------------------------------------------------------------------- -- DBScan@@ -382,12 +390,12 @@ -- | DBScan algorithm. {-# SCC dbscan #-} dbscan ::- ( MonadThrow m,- Ord e,- Num e,- Typeable e,- Show e,- Source r e+ ( MonadThrow m+ , Ord e+ , Num e+ , Typeable e+ , Show e+ , Source r e ) => -- | Distance measure to build the distance matrix of all points. DistFn r e ->@@ -404,86 +412,88 @@ | nPoints < 1 = throwM $ SizeNegativeException (Sz nPoints) | epsilon <= 0 = throwM $ DistanceInvalidException epsilon | otherwise =- let pointNeighbours = ifoldlInner collectNeighbours mempty distMat- allClusters = joinOverlapping . compute @B $ pointNeighbours- largeClusters = sfilter (\s -> IntSet.size s >= nPoints) allClusters- in return $ compute largeClusters- where- -- The distance matrix in the measure chosen by the distance function.- distMat = distFn points+ let pointNeighbours = ifoldlInner collectNeighbours mempty distMat+ allClusters = joinOverlapping . compute @B $ pointNeighbours+ largeClusters = sfilter (\s -> IntSet.size s >= nPoints) allClusters+ in return $ compute largeClusters+ where+ -- The distance matrix in the measure chosen by the distance function.+ distMat = distFn points - -- Function to collect the neighbours of a point within the search radius epsilon.- {-# SCC collectNeighbours #-}- collectNeighbours (_ :. n) acc d = if d <= epsilon then IntSet.insert n acc else acc+ -- Function to collect the neighbours of a point within the search radius epsilon.+ {-# SCC collectNeighbours #-}+ collectNeighbours (_ :. n) acc d = if d <= epsilon then IntSet.insert n acc else acc - -- Construct the overlap matrix of all clusters.- compareSets :: (IntSet -> IntSet -> Bool) -> Vector B IntSet -> Matrix D Bool- compareSets fn clVec =- let a = expandOuter sz const clVec- b = transpose a- compareMat = Massiv.zipWith fn a b- in compareMat- where- sz = size clVec+ -- Construct the overlap matrix of all clusters.+ compareSets :: (IntSet -> IntSet -> Bool) -> Vector B IntSet -> Matrix D Bool+ compareSets fn clVec =+ let a = expandOuter sz const clVec+ b = transpose a+ compareMat = Massiv.zipWith fn a b+ in compareMat+ where+ sz = size clVec - -- Overlap matrix. Checks if two sets have any overlap. Sets do overlap with themself.- overlap :: Vector B IntSet -> Matrix D Bool- overlap = compareSets (\a b -> not $ IntSet.disjoint a b)+ -- Overlap matrix. Checks if two sets have any overlap. Sets do overlap with themself.+ overlap :: Vector B IntSet -> Matrix D Bool+ overlap = compareSets (\a b -> not $ IntSet.disjoint a b) - -- Check if any set overlaps wiht **any** other set.- anyOtherOverlap :: Vector B IntSet -> Bool- anyOtherOverlap = Massiv.any (== True) . imap (\(m :. n) v -> if m == n then False else v) . overlap+ -- Check if any set overlaps wiht **any** other set.+ anyOtherOverlap :: Vector B IntSet -> Bool+ anyOtherOverlap = Massiv.any id . imap (\(m :. n) v -> (m /= n) && v) . overlap - -- Check if two sets are identical. Sets are identical to themself.- same :: Vector B IntSet -> Matrix D Bool- same = compareSets (==)+ -- Check if two sets are identical. Sets are identical to themself.+ same :: Vector B IntSet -> Matrix D Bool+ same = compareSets (==) - -- Join all overlapping clusters recursively.- {-# SCC joinOverlapping #-}- joinOverlapping :: Vector B IntSet -> Vector B IntSet- joinOverlapping clVec =- let -- The overlap matrix of the clusters.- ovlpMat = compute @U . overlap $ clVec- anyOvlp = anyOtherOverlap clVec+ -- Join all overlapping clusters recursively.+ {-# SCC joinOverlapping #-}+ joinOverlapping :: Vector B IntSet -> Vector B IntSet+ joinOverlapping clVec =+ let+ -- The overlap matrix of the clusters.+ ovlpMat = compute @U . overlap $ clVec+ anyOvlp = anyOtherOverlap clVec - -- Join all sets that have overlap but keep them redundantly (no reduction of the amount- -- of clusters).- joined =- ifoldlInner- (\(_ :. n) acc ovlp -> if ovlp then (clVec ! n) <> acc else acc)- mempty- ovlpMat+ -- Join all sets that have overlap but keep them redundantly (no reduction of the amount+ -- of clusters).+ joined =+ ifoldlInner+ (\(_ :. n) acc ovlp -> if ovlp then (clVec ! n) <> acc else acc)+ mempty+ ovlpMat - -- Find all sets at different indices that are the same. This is an upper triangular- -- matrix with the main diagonal being False.- sameMat =- compute @U- . imap (\(m :. n) v -> if m >= n then False else v)- . same- . compute @B- $ joined+ -- Find all sets at different indices that are the same. This is an upper triangular+ -- matrix with the main diagonal being False.+ sameMat =+ compute @U+ . imap (\(m :. n) v -> (m < n) && v)+ . same+ . compute @B+ $ joined - -- Remove all sets that are redundant. Redundancy is checked by two criteria:- -- 1. Is this cluster the same set of points as any other cluster? If yes, it is- -- redundant.- -- 2. Is this cluster isolated it is not redundant.- nonRed =- compute @B- . sifilter- ( \ix _ ->- let sameAsAnyOther = Massiv.any (== True) $ sameMat !> ix- in not sameAsAnyOther- )- $ joined- in if anyOvlp then joinOverlapping nonRed else clVec+ -- Remove all sets that are redundant. Redundancy is checked by two criteria:+ -- 1. Is this cluster the same set of points as any other cluster? If yes, it is+ -- redundant.+ -- 2. Is this cluster isolated it is not redundant.+ nonRed =+ compute @B+ . sifilter+ ( \ix _ ->+ let sameAsAnyOther = Massiv.any id $ sameMat !> ix+ in not sameAsAnyOther+ )+ $ joined+ in+ if anyOvlp then joinOverlapping nonRed else clVec ---------------------------------------------------------------------------------------------------- -- Hierarchical Cluster Analysis -- | Nodes of a dendrogram. data DendroNode e = DendroNode- { distance :: e,- cluster :: IntSet+ { distance :: e+ , cluster :: IntSet } deriving (Eq, Show, Generic) @@ -495,21 +505,22 @@ newtype Dendrogram e = Dendrogram {unDendro :: BinTree (DendroNode e)} deriving (Show, Eq, Generic) -instance ToJSON e => ToJSON (Dendrogram e)+instance (ToJSON e) => ToJSON (Dendrogram e) -instance FromJSON e => FromJSON (Dendrogram e)+instance (FromJSON e) => FromJSON (Dendrogram e) --- | An accumulator to finally build a dendrogram by a bottom-up algorithm. Not to be exposed in the--- API.+{- | An accumulator to finally build a dendrogram by a bottom-up algorithm. Not to be exposed in the+API.+-} type DendroAcc e = Vector B (Dendrogram e) -- | Manifest version of the dendrogram accumulator. type DendroAccM m e = MArray (PrimState m) B Ix1 (Dendrogram e) -- | Cut a 'Dendrogram' at a given distance and obtain all clusters from it.-cutDendroAt :: Ord e => Dendrogram e -> e -> Clusters+cutDendroAt :: (Ord e) => Dendrogram e -> e -> Clusters cutDendroAt dendro dist =- let nodes = takeLeafyBranchesWhile (\DendroNode {distance} -> distance >= dist) . unDendro $ dendro+ let nodes = takeLeafyBranchesWhile (\DendroNode{distance} -> distance >= dist) . unDendro $ dendro clusters = compute @B . Massiv.map cluster . compute @B $ nodes in clusters @@ -529,7 +540,7 @@ -- | Lance Williams formula to update distances. {-# SCC lanceWilliams #-} lanceWilliams ::- Fractional e =>+ (Fractional e) => -- | How to calculate distance between clusters of points. JoinStrat e -> -- | Number of points in cluster \(A\).@@ -547,42 +558,44 @@ -- | Updated distance \(D \(A \cup B, C\) e lanceWilliams js nA nB nC dAB dAC dBC = alpha1 * dAC + alpha2 * dBC + beta * dAB + gamma * abs (dAC - dBC)- where- nA' = fromIntegral nA- nB' = fromIntegral nB- nC' = fromIntegral nC- (alpha1, alpha2, beta, gamma) = case js of- SingleLinkage -> (1 / 2, 1 / 2, 0, - 1 / 2)- CompleteLinkage -> (1 / 2, 1 / 2, 0, 1 / 2)- Median -> (1 / 2, 1 / 2, - 1 / 4, 0)- UPGMA -> (nA' / (nA' + nB'), nB' / (nA' + nB'), 0, 0)- WPGMA -> (1 / 2, 1 / 2, 0, 0)- Centroid -> (nA' / (nA' + nB'), nB' / (nA' + nB'), - (nA' * nB') / ((nA' + nB') ^ (2 :: Int)), 0)- Ward -> ((nA' + nC') / (nA' + nB' + nC'), (nA' + nC') / (nA' + nB' + nC'), - (nA' + nC') / (nA' + nB' + nC'), 0)- LWFB b -> ((1 - b) / 2, (1 - b) / 2, b, 0)- LW a b c d -> (a, b, c, d)+ where+ nA' = fromIntegral nA+ nB' = fromIntegral nB+ nC' = fromIntegral nC+ (alpha1, alpha2, beta, gamma) = case js of+ SingleLinkage -> (1 / 2, 1 / 2, 0, -1 / 2)+ CompleteLinkage -> (1 / 2, 1 / 2, 0, 1 / 2)+ Median -> (1 / 2, 1 / 2, -1 / 4, 0)+ UPGMA -> (nA' / (nA' + nB'), nB' / (nA' + nB'), 0, 0)+ WPGMA -> (1 / 2, 1 / 2, 0, 0)+ Centroid -> (nA' / (nA' + nB'), nB' / (nA' + nB'), -(nA' * nB') / ((nA' + nB') ^ (2 :: Int)), 0)+ Ward -> ((nA' + nC') / (nA' + nB' + nC'), (nA' + nC') / (nA' + nB' + nC'), -(nA' + nC') / (nA' + nB' + nC'), 0)+ LWFB b -> ((1 - b) / 2, (1 - b) / 2, b, 0)+ LW a b c d -> (a, b, c, d) ---------------------------------------------------------------------------------------------------- -- Müllner Generic Hierarchical Clustering --- | A neighbourlist. At index @i@ of the vector it contains a tuple with the minimal distance of--- this cluster to any other cluster and the index of the other cluster.+{- | A neighbourlist. At index @i@ of the vector it contains a tuple with the minimal distance of+this cluster to any other cluster and the index of the other cluster.+-} type Neighbourlist r e = Vector r (e, Ix1) -- | A distance matrix. type DistanceMatrix r e = Matrix r e --- | Performance improved hierarchical clustering algorithm. @GENERIC_LINKAGE@ from figure 3,--- <https://arxiv.org/pdf/1109.2378.pdf>.+{- | Performance improved hierarchical clustering algorithm. @GENERIC_LINKAGE@ from figure 3,+<https://arxiv.org/pdf/1109.2378.pdf>.+-} {-# SCC hca #-} hca ::- ( MonadThrow m,- Manifest r e,- Manifest r (e, Ix1),- Load r Ix1 e,- Ord e,- Unbox e,- Fractional e+ ( MonadThrow m+ , Manifest r e+ , Manifest r (e, Ix1)+ , Load r Ix1 e+ , Ord e+ , Unbox e+ , Fractional e ) => DistFn r e -> JoinStrat e ->@@ -591,13 +604,15 @@ hca distFn joinStrat points | Massiv.isEmpty points = throwM $ SizeEmptyException (Sz nPoints) | otherwise = do- let -- The distance matrix from the points.+ let+ -- The distance matrix from the points. distMat = distFn points - -- Initial vector of nearest neighbour to each point.- nNghbr <- nearestNeighbours distMat+ -- Initial vector of nearest neighbour to each point.+ nNghbr <- nearestNeighbours distMat - let -- Initial priority queue of points. Has the minimum distance of all points.+ let+ -- Initial priority queue of points. Has the minimum distance of all points. pq = PQ.fromList . Massiv.toList . Massiv.imap (\k (d, n) -> (k, d, n)) $ nNghbr -- Set of points not joined yet. Initially all points. s = IntSet.fromDistinctAscList [0 .. nPoints - 1]@@ -606,31 +621,32 @@ makeArray @B @Ix1 Par (Sz nPoints)- (\p -> Dendrogram . Leaf $ DendroNode {distance = 0, cluster = IntSet.singleton p})+ (\p -> Dendrogram . Leaf $ DendroNode{distance = 0, cluster = IntSet.singleton p}) - distMatM <- return . unsafePerformIO . thaw $ distMat- nNghbrM <- return . unsafePerformIO . thaw $ nNghbr- dendroAccM <- return . unsafePerformIO . thaw $ dendroAcc+ distMatM <- return . unsafePerformIO . thaw $ distMat+ nNghbrM <- return . unsafePerformIO . thaw $ nNghbr+ dendroAccM <- return . unsafePerformIO . thaw $ dendroAcc - return . unsafePerformIO $ agglomerate joinStrat distMatM nNghbrM pq s dendroAccM- where- Sz (_mFeatures :. nPoints) = size points+ return . unsafePerformIO $ agglomerate joinStrat distMatM nNghbrM pq s dendroAccM+ where+ Sz (_mFeatures :. nPoints) = size points --- | Agglomerative clustering by the improved generic linkage algorithm. This is the main loop--- recursion L 10-43.+{- | Agglomerative clustering by the improved generic linkage algorithm. This is the main loop+recursion L 10-43.+-} {-# SCC agglomerate #-} agglomerate ::- ( MonadThrow m,- PrimMonad m,- MonadUnliftIO m,- PrimState m ~ RealWorld,- Manifest r e,- -- OuterSlice r Ix2 e,+ ( MonadThrow m+ , PrimMonad m+ , MonadUnliftIO m+ , PrimState m ~ RealWorld+ , Manifest r e+ , -- OuterSlice r Ix2 e, -- Manifest (R r) Ix1 e,- Manifest r (e, Ix1),- Shape r Ix1,- Fractional e,- Ord e+ Manifest r (e, Ix1)+ , Shape r Ix1+ , Fractional e+ , Ord e ) => -- | Join strategy for clusters and therefore how to calculate cluster-cluster distances. JoinStrat e ->@@ -649,44 +665,45 @@ agglomerate joinStrat distMat nNghbr pq s dendroAcc | IntSet.null s = throwM $ IndexException "No clusters left. This must never happen." | otherwise = do- -- Obtain candidates for the two clusters to join and the minimal distance in the priority queue.- candidates <- getJoinCandidates nNghbr pq+ -- Obtain candidates for the two clusters to join and the minimal distance in the priority queue.+ candidates <- getJoinCandidates nNghbr pq - -- If the distance between a b is not the minimal distance that the priority queue has found, the- -- neighbour list must be wrong and recalculated.- (a, b, delta, nNghbrU1, pqU1) <- recalculateNghbr candidates s distMat nNghbr pq+ -- If the distance between a b is not the minimal distance that the priority queue has found, the+ -- neighbour list must be wrong and recalculated.+ (a, b, delta, nNghbrU1, pqU1) <- recalculateNghbr candidates s distMat nNghbr pq - -- Remove the minimal element from the priority queue and join clusters a and b. The cluster- -- accumulator is reduced in its size: a is removed and b is updated with the joined cluster.- (newS, pqU2, newAcc) <- joinClusters a b delta s pqU1 dendroAcc+ -- Remove the minimal element from the priority queue and join clusters a and b. The cluster+ -- accumulator is reduced in its size: a is removed and b is updated with the joined cluster.+ (newS, pqU2, newAcc) <- joinClusters a b delta s pqU1 dendroAcc - -- Update the distance matrix in the row and column of b but not at (b,b) and not at (a,b) and- -- (b,a).- newDistMat <- updateDistMat joinStrat a b newS distMat newAcc+ -- Update the distance matrix in the row and column of b but not at (b,b) and not at (a,b) and+ -- (b,a).+ newDistMat <- updateDistMat joinStrat a b newS distMat newAcc - -- Redirect neighbours to b, if they previously pointed to a.- nNghbrU2 <- redirectNeighbours a b newS newDistMat nNghbrU1+ -- Redirect neighbours to b, if they previously pointed to a.+ nNghbrU2 <- redirectNeighbours a b newS newDistMat nNghbrU1 - -- Preserve a lower bound in priority queue and update the nearest neighbour list.- (nNghbrU3, pqU3) <- updateWithNewBDists b newS newDistMat nNghbrU2 pqU2+ -- Preserve a lower bound in priority queue and update the nearest neighbour list.+ (nNghbrU3, pqU3) <- updateWithNewBDists b newS newDistMat nNghbrU2 pqU2 - -- Update the neighbourlist and priority queue with the new distances to b.- (newNNghbr, newPQ) <- updateBNeighbour b s newDistMat nNghbrU3 pqU3+ -- Update the neighbourlist and priority queue with the new distances to b.+ (newNNghbr, newPQ) <- updateBNeighbour b s newDistMat nNghbrU3 pqU3 - -- If the problem has been reduced to a single cluster the algorithm is done and the final- -- dendrogram can be obtained from the accumulator at index b. Otherwise join further.- if IntSet.size newS == 1- then newAcc `readM` b- else agglomerate joinStrat newDistMat newNNghbr newPQ newS newAcc+ -- If the problem has been reduced to a single cluster the algorithm is done and the final+ -- dendrogram can be obtained from the accumulator at index b. Otherwise join further.+ if IntSet.size newS == 1+ then newAcc `readM` b+ else agglomerate joinStrat newDistMat newNNghbr newPQ newS newAcc --- | Obtain candidates for the clusters to join by looking at the minimal distance in the priority--- queue and the neighbourlist. L 11-13+{- | Obtain candidates for the clusters to join by looking at the minimal distance in the priority+queue and the neighbourlist. L 11-13+-} {-# SCC getJoinCandidates #-} getJoinCandidates ::- ( MonadThrow m,- PrimMonad m,- Manifest r (e, Ix1),- Ord e+ ( MonadThrow m+ , PrimMonad m+ , Manifest r (e, Ix1)+ , Ord e ) => MArray (PrimState m) r Ix1 (e, Ix1) -> PQ.HashPSQ Ix1 e Ix1 ->@@ -698,21 +715,22 @@ (_, b) <- nNghbr `readM` a return (a, b, d) --- | If the minimal distance @d@ found is not the distance between @a@ and @b@ recalculate the--- neighbour list, update the priority queue and obtain a new set of a,b and a distance between them.--- L 14-20.+{- | If the minimal distance @d@ found is not the distance between @a@ and @b@ recalculate the+neighbour list, update the priority queue and obtain a new set of a,b and a distance between them.+L 14-20.+-} {-# SCC recalculateNghbr #-} recalculateNghbr ::- ( MonadThrow m,- PrimMonad m,- MonadUnliftIO m,- PrimState m ~ RealWorld,- -- OuterSlice r Ix2 e,+ ( MonadThrow m+ , PrimMonad m+ , MonadUnliftIO m+ , PrimState m ~ RealWorld+ , -- OuterSlice r Ix2 e, -- Manifest (R r) Ix1 e,- Manifest r (e, Ix1),- Manifest r e,- Shape r Ix1,- Ord e+ Manifest r (e, Ix1)+ , Manifest r e+ , Shape r Ix1+ , Ord e ) => (Ix1, Ix1, e) -> IntSet ->@@ -741,14 +759,15 @@ (_, b) <- nNghbr `readM` a recalculateNghbr (a, b, newD) s distMat nNghbr newPQ --- | Joins the selected clusters \(A\) and \(B\) and updates the dendrogram accumulator at index b.--- A will not be removed so that the accumulator never shrinks.--- L 21-24+{- | Joins the selected clusters \(A\) and \(B\) and updates the dendrogram accumulator at index b.+A will not be removed so that the accumulator never shrinks.+L 21-24+-} {-# SCC joinClusters #-} joinClusters ::- ( MonadThrow m,- PrimMonad m,- Ord e+ ( MonadThrow m+ , PrimMonad m+ , Ord e ) => Ix1 -> Ix1 ->@@ -767,8 +786,8 @@ . Dendrogram $ Node ( DendroNode- { distance = d,- cluster = (cluster . root . unDendro $ clA) <> (cluster . root . unDendro $ clB)+ { distance = d+ , cluster = (cluster . root . unDendro $ clA) <> (cluster . root . unDendro $ clB) } ) (unDendro clA)@@ -778,15 +797,16 @@ let newS = IntSet.delete a s return (newS, newPQ, acc) --- | Update the distance matrix with a Lance-Williams update in the rows and columns of cluster b.--- L 25-27+{- | Update the distance matrix with a Lance-Williams update in the rows and columns of cluster b.+L 25-27+-} {-# SCC updateDistMat #-} updateDistMat ::- ( MonadThrow m,- PrimMonad m,- MonadUnliftIO m,- Manifest r e,- Fractional e+ ( MonadThrow m+ , PrimMonad m+ , MonadUnliftIO m+ , Manifest r e+ , Fractional e ) => JoinStrat e -> Ix1 ->@@ -799,31 +819,32 @@ | nDM /= nDM = throwM $ SizeMismatchException (Sz nDM) (Sz nCl) | mDM /= nDM = throwM $ SizeMismatchException (Sz mDM) (Sz nDM) | otherwise = do- dAB <- distMat `readM` (a :. b)- nA <- clSize a- nB <- clSize b- forIO_ ixV $ \ix -> do- dAX <- distMat `readM` (a :. ix)- nX <- clSize ix- modifyM_ distMat (\dBX -> return $ lanceWilliams js nA nB nX dAB dAX dBX) (ix :. b)- modifyM_ distMat (\dBX -> return $ lanceWilliams js nA nB nX dAB dAX dBX) (b :. ix)- return distMat- where- Sz (mDM :. nDM) = sizeOfMArray distMat- Sz nCl = sizeOfMArray dendroAcc- ixV = Massiv.fromList @U Par . IntSet.toAscList . IntSet.delete b $ s- clSize i = IntSet.size . cluster . root . unDendro <$> dendroAcc `readM` i+ dAB <- distMat `readM` (a :. b)+ nA <- clSize a+ nB <- clSize b+ forIO_ ixV $ \ix -> do+ dAX <- distMat `readM` (a :. ix)+ nX <- clSize ix+ modifyM_ distMat (\dBX -> return $ lanceWilliams js nA nB nX dAB dAX dBX) (ix :. b)+ modifyM_ distMat (\dBX -> return $ lanceWilliams js nA nB nX dAB dAX dBX) (b :. ix)+ return distMat+ where+ Sz (mDM :. nDM) = sizeOfMArray distMat+ Sz nCl = sizeOfMArray dendroAcc+ ixV = Massiv.fromList @U Par . IntSet.toAscList . IntSet.delete b $ s+ clSize i = IntSet.size . cluster . root . unDendro <$> dendroAcc `readM` i --- | Updates the neighbourlist. All elements with a smaller index than a, that had a as a nearest--- neighbour are blindly redirected to the union of a and b, now at index b.--- L 28-32+{- | Updates the neighbourlist. All elements with a smaller index than a, that had a as a nearest+neighbour are blindly redirected to the union of a and b, now at index b.+L 28-32+-} {-# SCC redirectNeighbours #-} redirectNeighbours ::- ( MonadThrow m,- PrimMonad m,- MonadUnliftIO m,- Manifest r (e, Ix1),- Manifest r e+ ( MonadThrow m+ , PrimMonad m+ , MonadUnliftIO m+ , Manifest r (e, Ix1)+ , Manifest r e ) => Ix1 -> Ix1 ->@@ -842,20 +863,21 @@ ) ix return nNghbr- where- ixV = compute @U . sfilter (< a) . Massiv.fromList @U Par . IntSet.toAscList $ s+ where+ ixV = compute @U . sfilter (< a) . Massiv.fromList @U Par . IntSet.toAscList $ s --- | Updates the list of nearest neighbours for all combinations that might have changed by--- recalculation with the joined cluster AB at index b.--- L 33-38+{- | Updates the list of nearest neighbours for all combinations that might have changed by+recalculation with the joined cluster AB at index b.+L 33-38+-} {-# SCC updateWithNewBDists #-} updateWithNewBDists ::- ( MonadThrow m,- MonadUnliftIO m,- PrimMonad m,- Manifest r e,- Manifest r (e, Ix1),- Ord e+ ( MonadThrow m+ , MonadUnliftIO m+ , PrimMonad m+ , Manifest r e+ , Manifest r (e, Ix1)+ , Ord e ) => Ix1 -> IntSet ->@@ -879,21 +901,22 @@ newPQ <- readTVarIO pqT return (nNghbr, newPQ)- where- ixV = compute @U . Massiv.sfilter (< b) . Massiv.fromList @U Par . IntSet.toAscList $ s+ where+ ixV = compute @U . Massiv.sfilter (< b) . Massiv.fromList @U Par . IntSet.toAscList $ s --- | Updates the list of nearest neighbours and the priority queue at key b.--- L 39-40+{- | Updates the list of nearest neighbours and the priority queue at key b.+L 39-40+-} {-# SCC updateBNeighbour #-} updateBNeighbour ::- ( MonadThrow m,- PrimMonad m,- RealWorld ~ PrimState m,- MonadUnliftIO m,- Manifest r (e, Ix1),- Manifest r e,- Shape r Ix1,- Ord e+ ( MonadThrow m+ , PrimMonad m+ , RealWorld ~ PrimState m+ , MonadUnliftIO m+ , Manifest r (e, Ix1)+ , Manifest r e+ , Shape r Ix1+ , Ord e ) => Ix1 -> IntSet ->@@ -910,22 +933,23 @@ writeM nNghbr b newNeighbourB let newPQ = pqAdjust (const distB) neighbourB pq return (nNghbr, newPQ)- where- Sz nNeighbours = sizeOfMArray nNghbr+ where+ Sz nNeighbours = sizeOfMArray nNghbr --- | Find the nearest neighbour for each point from a distance matrix. For each point it stores the--- minimum distance and the index of the other point, that is the nearest neighbour but at a higher--- index.+{- | Find the nearest neighbour for each point from a distance matrix. For each point it stores the+minimum distance and the index of the other point, that is the nearest neighbour but at a higher+index.+-} {-# SCC nearestNeighbours #-} nearestNeighbours ::- ( MonadThrow m,- Manifest r e,- Manifest r (e, Ix1),- Load r Ix1 e,- -- OuterSlice r Ix2 e,+ ( MonadThrow m+ , Manifest r e+ , Manifest r (e, Ix1)+ , Load r Ix1 e+ , -- OuterSlice r Ix2 e, -- Source (R r) Ix1 e,- Ord e,- Unbox e+ Ord e+ , Unbox e ) => Matrix r e -> m (Vector r (e, Ix1))@@ -933,24 +957,25 @@ | m /= n = throwM $ IndexException "Distance matrix is not square" | m == 0 = throwM $ IndexException "Distance matrix is empty" | otherwise =- let rows = compute @B . outerSlices $ distMat- minDistIx =- Massiv.imap (\i v -> unsafePerformIO . minDistAtVec i . compute @U $ v) . init $ rows- in return . compute $ minDistIx- where- Sz (m :. n) = size distMat+ let rows = compute @B . outerSlices $ distMat+ minDistIx =+ Massiv.imap (\i v -> unsafePerformIO . minDistAtVec i . compute @U $ v) . init $ rows+ in return . compute $ minDistIx+ where+ Sz (m :. n) = size distMat --- | Make a search row for distances. Takes row x from a distance matrix and zips them with their--- column index. Then keeps only the valid elements of the row, that are still part of the available--- points. A minimum or maximum search can be performed on the resulting vector and a valid pair of--- distance and index can be obtained.+{- | Make a search row for distances. Takes row x from a distance matrix and zips them with their+column index. Then keeps only the valid elements of the row, that are still part of the available+points. A minimum or maximum search can be performed on the resulting vector and a valid pair of+distance and index can be obtained.+-} searchRow ::- ( PrimMonad m,- RealWorld ~ PrimState m,- MonadThrow m,- MonadUnliftIO m,- Manifest r e,- Manifest r (e, Ix1)+ ( PrimMonad m+ , RealWorld ~ PrimState m+ , MonadThrow m+ , MonadUnliftIO m+ , Manifest r e+ , Manifest r (e, Ix1) ) => Ix1 -> IntSet ->@@ -959,8 +984,7 @@ searchRow x s dm = makeMArray Par (size ixV) $ \ix -> do dmIx <- ixV !? ix- val <- (dm `readM` (x :. dmIx)) >>= \dist -> return (dist, dmIx)- return val- where- ixV :: Vector U Ix1- ixV = compute @U . sfilter (> x) . Massiv.fromList @U Par . IntSet.toAscList $ s+ (dm `readM` (x :. dmIx)) >>= \dist -> return (dist, dmIx)+ where+ ixV :: Vector U Ix1+ ixV = compute @U . sfilter (> x) . Massiv.fromList @U Par . IntSet.toAscList $ s