diff --git a/Changelog.md b/Changelog.md --- a/Changelog.md +++ b/Changelog.md @@ -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 diff --git a/ConClusion.cabal b/ConClusion.cabal --- a/ConClusion.cabal +++ b/ConClusion.cabal @@ -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 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 diff --git a/app/ConClusion.hs b/app/ConClusion.hs --- a/app/ConClusion.hs +++ b/app/ConClusion.hs @@ -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 ' ') diff --git a/src/ConClusion/Array/Conversion.hs b/src/ConClusion/Array/Conversion.hs --- a/src/ConClusion/Array/Conversion.hs +++ b/src/ConClusion/Array/Conversion.hs @@ -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 diff --git a/src/ConClusion/Array/Util.hs b/src/ConClusion/Array/Util.hs --- a/src/ConClusion/Array/Util.hs +++ b/src/ConClusion/Array/Util.hs @@ -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 diff --git a/src/ConClusion/BinaryTree.hs b/src/ConClusion/BinaryTree.hs --- a/src/ConClusion/BinaryTree.hs +++ b/src/ConClusion/BinaryTree.hs @@ -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 diff --git a/src/ConClusion/Chemistry/Topology.hs b/src/ConClusion/Chemistry/Topology.hs --- a/src/ConClusion/Chemistry/Topology.hs +++ b/src/ConClusion/Chemistry/Topology.hs @@ -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 . --- --- Diherdrals require a special metric, see . -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 . + +Diherdrals require a special metric, see . +-} +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 -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 +-} +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 diff --git a/src/ConClusion/Numeric/Statistics.hs b/src/ConClusion/Numeric/Statistics.hs --- a/src/ConClusion/Numeric/Statistics.hs +++ b/src/ConClusion/Numeric/Statistics.hs @@ -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 !> 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, --- . +{- | Performance improved hierarchical clustering algorithm. @GENERIC_LINKAGE@ from figure 3, +. +-} {-# 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