ConClusion 0.1.0 → 0.2.0

raw patch · 7 files changed

+226/−216 lines, 7 filesdep ~aesondep ~basedep ~textPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: aeson, base, text

API changes (from Hackage documentation)

Files

• ConClusion.cabal +10/−8
• src/ConClusion/Array/Conversion.hs +46/−0
• src/ConClusion/Array/Util.hs +102/−0
• src/ConClusion/BinaryTree.hs +60/−0
• src/ConClusion/Chemistry/Topology.hs +5/−4
• src/ConClusion/Numeric/Data.hs +0/−198
• src/ConClusion/Numeric/Statistics.hs +3/−6

@@ -1,11 +1,11 @@ cabal-version: 1.12 --- This file has been generated from package.yaml by hpack version 0.34.4.+-- This file has been generated from package.yaml by hpack version 0.34.6. -- -- see: https://github.com/sol/hpack  name:           ConClusion-version:        0.1.0+version:        0.2.0 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@@ -27,8 +27,10 @@  library   exposed-modules:+      ConClusion.Array.Conversion+      ConClusion.Array.Util+      ConClusion.BinaryTree       ConClusion.Chemistry.Topology-      ConClusion.Numeric.Data       ConClusion.Numeric.Statistics   other-modules:       Paths_ConClusion@@ -60,9 +62,9 @@       NamedFieldPuns   ghc-options: -Wall -Wno-unused-top-binds -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wpartial-fields -Wredundant-constraints   build-depends:-      aeson ==1.5.*+      aeson >=1.5 && <2.1     , attoparsec >=0.13.0.0 && <0.15-    , base >=4.7 && <4.16+    , base >=4.7 && <4.17     , containers >=0.6.0.0 && <0.7     , formatting >=7.1.0 && <7.2     , hmatrix >=0.20.0 && <0.21@@ -104,9 +106,9 @@   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.*+    , aeson >=1.5 && <2.1     , attoparsec >=0.13.0.0 && <0.15-    , base >=4.7 && <4.16+    , base >=4.7 && <4.17     , cmdargs >=0.10.0 && <0.11     , containers >=0.6.0.0 && <0.7     , formatting >=7.1.0 && <7.2@@ -115,5 +117,5 @@     , optics >=0.3 && <0.5     , psqueues >=0.2.7.0 && <0.3     , rio >=0.1.13.0 && <0.2-    , text >=1.2.0.0 && <1.3+    , text >=1.2.0.0 && <2.1   default-language: Haskell2010

@@ -0,0 +1,46 @@+-- |+-- 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,+  )+where++import Data.Massiv.Array as Massiv hiding (IndexException)+import Data.Massiv.Array.Manifest.Vector as Massiv+import Numeric.LinearAlgebra as LA hiding (magnitude, (<>))+import RIO+import qualified RIO.Vector.Storable as VectorS++-- | Converts a vector from the HMatrix package to the Massiv representation.+{-# SCC vecH2M #-}+vecH2M :: (Element e, Mutable r e, Load r Ix1 e) => VectorS.Vector e -> Massiv.Vector r e+vecH2M hVec = fromVector' Seq (Sz $ VectorS.length hVec) hVec++-- | Converts a vector from the Massiv representation to the HMatrix representation.+{-# SCC vecM2H #-}+vecM2H :: (Manifest r e, Load r Ix1 e, Element e) => Massiv.Vector r e -> LA.Vector e+vecM2H = Massiv.toVector++-- | Converts a matrix from the HMatrix representation to the Massiv representation.+{-# SCC matH2M #-}+matH2M :: (Mutable 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++-- | 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

@@ -0,0 +1,102 @@+-- |+-- 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,+  )+where++import Data.Massiv.Array as Massiv hiding (IndexException)+import RIO+import System.IO.Unsafe (unsafePerformIO)++-- | Exception regarding indexing in some kind of aaray.+newtype IndexException = IndexException String deriving (Show)++instance Exception IndexException++-- | Magnitude of a vector (length).+magnitude :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> e+magnitude v = sqrt $ v !.! v++-- | Normalise a vector.+normalise :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e+normalise v = v .* (1 / magnitude v)++-- | Angle between two vectors.+angle :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e -> e+angle a b = acos $ a !.! b / (magnitude a * magnitude b)++-- | Find the minimal distance in a distance matrix, which is not the main diagonal.+{-# SCC minDistAt #-}+minDistAt ::+  ( 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++-- | Find the minimal element of a vector, which is at a larger than the supplied index.+minDistAtVec ::+  ( Manifest r e,+    MonadThrow m,+    Ord e+  ) =>+  Ix1 ->+  Massiv.Vector r e ->+  m (e, Ix1)+minDistAtVec ixStart vec+  | isEmpty vec = throwM $ SizeEmptyException (Massiv.size vec)+  | ixStart >= nElems = throwM $ IndexOutOfBoundsException (Sz nElems) ixStart+  | 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++-- | Like 'Massiv.minimumM' but also returns the index of the minimal element.+iMinimumM ::+  ( 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)++    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

@@ -0,0 +1,60 @@+-- |+-- 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,+  )+where++import Data.Aeson hiding (Array)+import Data.Massiv.Array as Massiv hiding (IndexException)+import RIO++-- | A binary tree.+data BinTree e = Leaf e | Node e (BinTree e) (BinTree e)+  deriving (Eq, Show, Generic)++instance (FromJSON e) => FromJSON (BinTree e)++instance (ToJSON e) => ToJSON (BinTree e)++instance Functor BinTree where+  fmap f (Leaf a) = Leaf (f a)+  fmap f (Node a l r) = Node (f a) (fmap f l) (fmap f r)++-- | Look at the root of a binary tree.+root :: BinTree e -> e+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.+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++-- | 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

@@ -26,8 +26,9 @@   ) where -import ConClusion.Numeric.Data hiding (angle)-import qualified ConClusion.Numeric.Data as ND+import ConClusion.Array.Conversion+import ConClusion.Array.Util hiding (angle)+import qualified ConClusion.Array.Util as ArrayUtil import Data.Attoparsec.Text hiding (D) import Data.Foldable import Data.Massiv.Array as Massiv hiding (B, D, forM)@@ -135,7 +136,7 @@     coordC <- compute @U <$> (coordinates !?> c)     vecAB <- coordA .-. coordB     vecCB <- coordC .-. coordB-    return $ ND.angle vecAB vecCB+    return $ ArrayUtil.angle vecAB vecCB  -- | Calculates the dihedral angle defined by four atoms. Respects rotation direction. Obtains the -- result in radian.@@ -157,7 +158,7 @@           if vecH2M normVecRot !.! vecAB < 0             then -1             else 1-    return $ rotDir * ND.angle planeABC planeBCD+    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.

@@ -1,198 +0,0 @@--- |--- Module      : ConClusion.Numeric.Data--- Description : Type castings, conversions and utilities--- Copyright   : Phillip Seeber, 2021--- License     : AGPL-3--- Maintainer  : phillip.seeber@googlemail.com--- Stability   : experimental--- Portability : POSIX, Windows-module ConClusion.Numeric.Data-  ( -- * Conversion of array types.--    -- Conversion between HMatrix and Massiv's array types-    IndexException (..),-    vecH2M,-    vecM2H,-    matH2M,-    matM2H,--    -- * Array Processing-    magnitude,-    normalise,-    angle,-    minDistAt,-    minDistAtVec,-    iMinimumM,--    -- * Utilities-    printMat,--    -- * Binary Trees-    BinTree (..),-    root,-    takeBranchesWhile,-    takeLeafyBranchesWhile,-  )-where--import Data.Aeson hiding (Array)-import Data.Massiv.Array as Massiv hiding (IndexException)-import Data.Massiv.Array.Manifest.Vector as Massiv-import Formatting-import Numeric.LinearAlgebra as LA hiding (magnitude, (<>))-import RIO-import qualified RIO.Vector.Storable as VectorS-import System.IO.Unsafe (unsafePerformIO)---- | Exception regarding indexing in some kind of aaray.-newtype IndexException = IndexException String deriving (Show)--instance Exception IndexException---- | Converts a vector from the HMatrix package to the Massiv representation.-{-# SCC vecH2M #-}-vecH2M :: (Element e, Mutable r e, Load r Ix1 e) => VectorS.Vector e -> Massiv.Vector r e-vecH2M hVec = fromVector' Seq (Sz $ VectorS.length hVec) hVec---- | Converts a vector from the Massiv representation to the HMatrix representation.-{-# SCC vecM2H #-}-vecM2H :: (Manifest r e, Load r Ix1 e, Element e) => Massiv.Vector r e -> LA.Vector e-vecM2H = Massiv.toVector---- | Converts a matrix from the HMatrix representation to the Massiv representation.-{-# SCC matH2M #-}-matH2M :: (Mutable 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---- | 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---- | Magnitude of a vector (length).-magnitude :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> e-magnitude v = sqrt $ v !.! v---- | Normalise a vector.-normalise :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e-normalise v = v .* (1 / magnitude v)---- | Angle between two vectors.-angle :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e -> e-angle a b = acos $ a !.! b / (magnitude a * magnitude b)---- | Find the minimal distance in a distance matrix, which is not the main diagonal.-{-# SCC minDistAt #-}-minDistAt ::-  ( 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---- | Find the minimal element of a vector, which is at a larger than the supplied index.-minDistAtVec ::-  ( Manifest r e,-    MonadThrow m,-    Ord e-  ) =>-  Ix1 ->-  Massiv.Vector r e ->-  m (e, Ix1)-minDistAtVec ixStart vec-  | isEmpty vec = throwM $ SizeEmptyException (Massiv.size vec)-  | ixStart >= nElems = throwM $ IndexOutOfBoundsException (Sz nElems) ixStart-  | 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---- | Like 'Massiv.minimumM' but also returns the index of the minimal element.-iMinimumM ::-  ( 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)--    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---- | Quickly print a matrix with numerical values-printMat :: (Source r e, Real e) => Massiv.Matrix r e -> Massiv.Matrix D Text-printMat mat = Massiv.map (sformat (left 4 ' ' %. fixed 2)) mat--------------------------------------------------------------------------------------------------------- Binary Trees.---- | A binary tree.-data BinTree e = Leaf e | Node e (BinTree e) (BinTree e)-  deriving (Eq, Show, Generic)--instance (FromJSON e) => FromJSON (BinTree e)--instance (ToJSON e) => ToJSON (BinTree e)--instance Functor BinTree where-  fmap f (Leaf a) = Leaf (f a)-  fmap f (Node a l r) = Node (f a) (fmap f l) (fmap f r)---- | Look at the root of a binary tree.-root :: BinTree e -> e-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.-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---- | 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

@@ -38,7 +38,9 @@   ) where -import ConClusion.Numeric.Data hiding (normalise)+import ConClusion.Array.Conversion+import ConClusion.Array.Util hiding (normalise)+import ConClusion.BinaryTree import Data.Aeson hiding (Array) import Data.Complex import qualified Data.HashPSQ as PQ@@ -64,8 +66,6 @@     Load r1 Ix1 (Complex Double),     Load r2 Ix1 (Complex Double),     Load r3 Ix1 e,-    -- Resize r3 Ix2,-    Load r3 Ix2 e,     MonadThrow m   ) =>   Matrix r3 e ->@@ -579,9 +579,6 @@   ( MonadThrow m,     Mutable r e,     Mutable r (e, Ix1),-    -- Manifest (R r) Ix1 e,-    -- OuterSlice r Ix2 e,-    Shape r Ix1,     Load r Ix1 e,     Ord e,     Unbox e,