ConClusion 0.0.2 → 0.1.0

raw patch · 6 files changed

+98/−83 lines, 6 filesdep ~massivPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: massiv

API changes (from Hackage documentation)

Files

• Changelog.md +4/−0
• ConClusion.cabal +3/−3
• app/ConClusion.hs +3/−3
• src/ConClusion/Chemistry/Topology.hs +2/−2
• src/ConClusion/Numeric/Data.hs +12/−11
• src/ConClusion/Numeric/Statistics.hs +74/−64

@@ -1,5 +1,9 @@ # Changelog +## 0.1.0+  - updates to Massiv 1.0.0.0; changes lots of array types.+  - infrastructure updates+ ## 0.0.2   - switching from `PSQueue` to `psqueues` as the former is not maintained   - raise upper bounds of base, allows GHC 9.0

@@ -5,7 +5,7 @@ -- see: https://github.com/sol/hpack  name:           ConClusion-version:        0.0.2+version:        0.1.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@@ -66,7 +66,7 @@     , containers >=0.6.0.0 && <0.7     , formatting >=7.1.0 && <7.2     , hmatrix >=0.20.0 && <0.21-    , massiv >=0.6.0.0 && <0.7+    , massiv >=1.0.0.0 && <1.1     , psqueues >=0.2.7.0 && <0.3     , rio >=0.1.13.0 && <0.2   default-language: Haskell2010@@ -111,7 +111,7 @@     , containers >=0.6.0.0 && <0.7     , formatting >=7.1.0 && <7.2     , hmatrix >=0.20.0 && <0.21-    , massiv >=0.6.0.0 && <0.7+    , 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

@@ -621,7 +621,7 @@   valueLines <- case clusters of     Nothing -> return . Massiv.fold . Massiv.map (writeLine . Left) . innerSlices $ mat     Just cl -> do-      let points = innerSlices mat :: Vector Massiv.D (Vector M Double)+      let points = innerSlices mat           labeledCl = compute . Massiv.imap (,) $ cl :: Vector Massiv.B (Int, IntSet)        -- Anotate all points with their assignment to a cluster number. If a point is from DBScan, it@@ -666,12 +666,12 @@           <> bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm F.% "," F.% iForm F.% "," F.% iForm)) "A cos" a b c d      -- Print a vector of numeric values.-    doubleVecB :: Source r Ix1 Double => Vector r Double -> TB.Builder+    doubleVecB :: Source r Double => Vector r Double -> TB.Builder     doubleVecB v = Massiv.fold . Massiv.map dForm $ v      -- Writer for a value line.     writeLine ::-      (Source r Ix1 Double) => Either (Vector r Double) (Maybe Int, Vector r Double) -> TB.Builder+      (Source r Double) => Either (Vector r Double) (Maybe Int, Vector r Double) -> TB.Builder     writeLine (Left vec) = doubleVecB vec <> "\n"     writeLine (Right (cl, vec)) = bformat iForm (fromMaybe (-1) cl) <> doubleVecB vec <> "\n" 

@@ -159,7 +159,7 @@             else 1     return $ rotDir * ND.angle planeABC planeBCD --- | Calculates the a metric value of the dihedral angle defined by four atoms. This must create 2+-- | 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)@@ -207,4 +207,4 @@   m (Massiv.Matrix DL Double) getFeatures sel trj = traverse toCols trj >>= concatM (Dim 1)   where-    toCols v = expandInner @Ix2 (Sz 1) const . compute @U <$> getFeature sel v+    toCols v = expandInner @U @Ix2 (Sz 1) const . compute @U <$> getFeature sel v

@@ -51,17 +51,17 @@  -- | Converts a vector from the HMatrix package to the Massiv representation. {-# SCC vecH2M #-}-vecH2M :: (Element e, Mutable r Ix1 e) => VectorS.Vector e -> Massiv.Vector r e+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 Ix1 e, Element e) => Massiv.Vector r e -> LA.Vector e+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 Ix1 e, Element e) => LA.Matrix e -> Massiv.Matrix r e+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@@ -69,27 +69,27 @@  -- | Converts a matrix from Massiv to HMatrix representation. {-# SCC matM2H #-}-matM2H :: (Manifest r Ix1 e, Element e, Resize r Ix2, Load r Ix2 e) => Massiv.Matrix r e -> LA.Matrix e+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 Ix1 e, Floating e) => Massiv.Vector r e -> e+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 Ix1 e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e+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 Ix1 e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e -> e+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 Ix2 e,+  ( Manifest r e,     MonadThrow m,     Ord e   ) =>@@ -107,7 +107,7 @@  -- | Find the minimal element of a vector, which is at a larger than the supplied index. minDistAtVec ::-  ( Manifest r Ix1 e,+  ( Manifest r e,     MonadThrow m,     Ord e   ) =>@@ -131,8 +131,9 @@  -- | Like 'Massiv.minimumM' but also returns the index of the minimal element. iMinimumM ::-  ( Manifest r ix a,+  ( Manifest r a,     MonadThrow m,+    Index ix,     Ord a   ) =>   Array r ix a ->@@ -149,7 +150,7 @@     chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc  -- | Quickly print a matrix with numerical values-printMat :: (Source r Ix2 e, Real e) => Massiv.Matrix r e -> Massiv.Matrix D Text+printMat :: (Source r e, Real e) => Massiv.Matrix r e -> Massiv.Matrix D Text printMat mat = Massiv.map (sformat (left 4 ' ' %. fixed 2)) mat  ----------------------------------------------------------------------------------------------------

@@ -55,11 +55,16 @@ -- | Solves eigenvalue problem of a square matrix and obtains its eigenvalues and eigenvectors. {-# SCC eig #-} eig ::-  ( Mutable r1 Ix1 (Complex Double),-    Mutable r2 Ix1 (Complex Double),+  ( -- Mutable r1 Ix1 (Complex Double),+    -- Mutable r2 Ix1 (Complex Double),     LA.Field e,-    Manifest r3 Ix1 e,-    Resize r3 Ix2,+    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,+    -- Resize r3 Ix2,     Load r3 Ix2 e,     MonadThrow m   ) =>@@ -78,10 +83,10 @@ eigSort ::   ( Load r2 Ix2 e,     MonadThrow m,-    Source r1 Ix1 e,-    Source r2 Ix2 e,-    Mutable r1 Ix1 e,-    Mutable r2 Ix2 e,+    Source r1 e,+    Source r2 e,+    Mutable r1 e,+    Mutable r2 e,     Unbox e,     Ord e   ) =>@@ -141,9 +146,9 @@ -- constructed from the eigenvectors associated to the largest eigenvalues. {-# SCC transformToPCABasis #-} transformToPCABasis ::-  ( Source (R r) Ix2 e,-    Extract r Ix2 e,-    Mutable r Ix2 e,+  ( -- Source (R r) Ix2 e,+    -- Extract r Ix2 e,+    Mutable r e,     Numeric r e,     MonadThrow m   ) =>@@ -177,10 +182,9 @@ {-# SCC pca #-} pca ::   ( Numeric r Double,-    Mutable r Ix2 Double,-    Manifest r Ix1 Double,-    Source (R r) Ix2 Double,-    Extract r Ix2 Double,+    Mutable r Double,+    Load r Ix1 Double,+    Load r Ix2 Double,     MonadThrow m   ) =>   -- | Dimensionalty after PCA transformation.@@ -238,11 +242,11 @@ -- mean deviation form. {-# SCC meanDeviation #-} meanDeviation ::-  ( Source r Ix2 e,+  ( Source r e,     Fractional e,     Unbox e,     Numeric r e,-    Mutable r Ix2 e+    Mutable r e   ) =>   Matrix r e ->   Matrix r e@@ -260,7 +264,7 @@ -- -- The feature matrix should be in mean deviation form, see 'meanDeviation'. {-# SCC covariance #-}-covariance :: (Numeric r e, Mutable r Ix2 e, Fractional e) => Matrix r e -> Matrix r e+covariance :: (Numeric r e, Mutable 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@@ -271,15 +275,14 @@     Unbox e,     Numeric r e,     Fractional e,-    Source r Ix2 e,-    Mutable r Ix2 e+    Mutable r e   ) =>   Array r Ix2 e ->   Array r Ix2 e normalise mat =   let absMat = Massiv.map abs mat       maxPerRow = compute @U . foldlInner max 0 $ absMat-      divMat = compute . Massiv.map (1 /) . expandInner @Ix2 (Sz n) const $ maxPerRow+      divMat = compute . Massiv.map (1 /) . expandInner @U @Ix2 (Sz n) const $ maxPerRow    in divMat !*! mat   where     Sz (_ :. n) = size mat@@ -292,7 +295,7 @@  -- | Builds the distance measures in a permutation matrix/distance matrix. buildDistMat ::-  (Mutable r Ix2 e) =>+  (Manifest r e) =>   -- | Zip function to combine the elements of vectors \(\mathbf{a}\) \(\mathbf{b}\). Usually @(-)@.   -- \( f(\mathbf{a}_i, \mathbf{b}_i) = \mathbf{c} \)   (e -> e -> a) ->@@ -305,7 +308,7 @@   -- | Resulting distance matrix.   Matrix D a buildDistMat zipFn foldFn acc mat =-  let a = transposeOuter . expandOuter @Ix3 (Sz n) const $ mat+  let a = transposeOuter @D @Ix3 . expandOuter (Sz n) const $ mat       b = transposeInner a       ab = Massiv.zipWith zipFn a b       d = foldlInner foldFn acc ab@@ -318,7 +321,7 @@ --   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 :: (Mutable r Ix2 e, Floating e) => Int -> DistFn r e+lpNorm :: (Mutable r e, Floating e) => Int -> DistFn r e lpNorm p = compute . buildDistMat zipFn foldFn 0   where     zipFn a b = (^ p) . abs $ a - b@@ -329,7 +332,7 @@ --   d(\mathbf{a}, \mathbf{b}) = \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert -- \] {-# SCC manhattan #-}-manhattan :: (Mutable r Ix2 e, Floating e) => DistFn r e+manhattan :: (Mutable 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\).@@ -337,7 +340,7 @@ --   d(\mathbf{a}, \mathbf{b}) = \sqrt{\sum \limits_{i=1}^n (\mathbf{a}_i - \mathbf{b}_i)^2} -- \] {-# SCC euclidean #-}-euclidean :: (Mutable r Ix2 e, Floating e) => DistFn r e+euclidean :: (Mutable r e, Floating e) => DistFn r e euclidean = lpNorm 2  -- | Mahalanobis distance between points. Suitable for non correlated axes.@@ -346,9 +349,9 @@ -- \] -- where \(\mathbf{S}\) is the covariance matrix. {-# SCC mahalanobis #-}-mahalanobis :: (Unbox e, Numeric r e, Mutable r Ix2 e, Mutable r Ix1 e, LA.Field e) => DistFn r e+mahalanobis :: (Unbox e, Numeric r e, Mutable r e, LA.Field e, Load r Ix1 e) => DistFn r e mahalanobis points =-  let a = transposeOuter . expandOuter @Ix3 (Sz n) const $ points+  let a = transposeOuter @D @Ix3 . expandOuter (Sz n) const $ points       b = transposeInner a       abDiff = compute @U $ a !-! b       ixArray = makeArray @U @Ix2 @Ix2 Par (Sz $ n :. n) id@@ -384,7 +387,7 @@     Num e,     Typeable e,     Show e,-    Source r Ix2 e+    Source r e   ) =>   -- | Distance measure to build the distance matrix of all points.   DistFn r e ->@@ -416,7 +419,7 @@     -- Construct the overlap matrix of all clusters.     compareSets :: (IntSet -> IntSet -> Bool) -> Vector B IntSet -> Matrix D Bool     compareSets fn clVec =-      let a = expandOuter @Ix2 sz const clVec+      let a = expandOuter sz const clVec           b = transpose a           compareMat = Massiv.zipWith fn a b        in compareMat@@ -574,11 +577,12 @@ {-# SCC hca #-} hca ::   ( MonadThrow m,-    Mutable r Ix1 e,-    Mutable r Ix2 e,-    Mutable r Ix1 (e, Ix1),-    Manifest (R r) Ix1 e,-    OuterSlice r Ix2 e,+    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,     Fractional e@@ -588,7 +592,7 @@   Matrix r e ->   m (Dendrogram e) hca distFn joinStrat points-  | Massiv.isEmpty points = throw $ SizeEmptyException (Sz nPoints)+  | Massiv.isEmpty points = throwM $ SizeEmptyException (Sz nPoints)   | otherwise = do     let -- The distance matrix from the points.         distMat = distFn points@@ -623,10 +627,11 @@     PrimMonad m,     MonadUnliftIO m,     PrimState m ~ RealWorld,-    Mutable r Ix2 e,-    OuterSlice r Ix2 e,-    Manifest (R r) Ix1 e,-    Mutable r Ix1 (e, Ix1),+    Mutable r e,+    -- OuterSlice r Ix2 e,+    -- Manifest (R r) Ix1 e,+    Mutable r (e, Ix1),+    Shape r Ix1,     Fractional e,     Ord e   ) =>@@ -645,7 +650,7 @@   -- | The final dendrogram, after all clusters have been joined.   m (Dendrogram e) agglomerate joinStrat distMat nNghbr pq s dendroAcc-  | IntSet.null s = throw $ IndexException "No clusters left. This must never happen."+  | 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@@ -683,7 +688,7 @@ getJoinCandidates ::   ( MonadThrow m,     PrimMonad m,-    Mutable r Ix1 (e, Ix1),+    Mutable r (e, Ix1),     Ord e   ) =>   MArray (PrimState m) r Ix1 (e, Ix1) ->@@ -705,10 +710,11 @@     PrimMonad m,     MonadUnliftIO m,     PrimState m ~ RealWorld,-    OuterSlice r Ix2 e,-    Manifest (R r) Ix1 e,-    Mutable r Ix1 (e, Ix1),-    Mutable r Ix2 e,+    -- OuterSlice r Ix2 e,+    -- Manifest (R r) Ix1 e,+    Mutable r (e, Ix1),+    Mutable r e,+    Shape r Ix1,     Ord e   ) =>   (Ix1, Ix1, e) ->@@ -782,7 +788,7 @@   ( MonadThrow m,     PrimMonad m,     MonadUnliftIO m,-    Mutable r Ix2 e,+    Mutable r e,     Fractional e   ) =>   JoinStrat e ->@@ -806,8 +812,8 @@       modifyM_ distMat (\dBX -> return $ lanceWilliams js nA nB nX dAB dAX dBX) (b :. ix)     return distMat   where-    Sz (mDM :. nDM) = msize distMat-    Sz nCl = msize dendroAcc+    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 @@ -819,8 +825,8 @@   ( MonadThrow m,     PrimMonad m,     MonadUnliftIO m,-    Mutable r Ix1 (e, Ix1),-    Mutable r Ix2 e+    Mutable r (e, Ix1),+    Mutable r e   ) =>   Ix1 ->   Ix1 ->@@ -850,8 +856,8 @@   ( MonadThrow m,     MonadUnliftIO m,     PrimMonad m,-    Mutable r Ix2 e,-    Mutable r Ix1 (e, Ix1),+    Mutable r e,+    Mutable r (e, Ix1),     Ord e   ) =>   Ix1 ->@@ -879,15 +885,17 @@   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.+-- | 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,-    Mutable r Ix1 (e, Ix1),-    Mutable r Ix2 e,+    Mutable r (e, Ix1),+    Mutable r e,+    Shape r Ix1,     Ord e   ) =>   Ix1 ->@@ -906,7 +914,7 @@       let newPQ = pqAdjust (const distB) neighbourB pq       return (nNghbr, newPQ)   where-    Sz nNeighbours = msize nNghbr+    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@@ -914,10 +922,11 @@ {-# SCC nearestNeighbours #-} nearestNeighbours ::   ( MonadThrow m,-    Mutable r Ix1 e,-    Mutable r Ix1 (e, Ix1),-    OuterSlice r Ix2 e,-    Source (R r) Ix1 e,+    Mutable r e,+    Mutable r (e, Ix1),+    Load r Ix1 e,+    -- OuterSlice r Ix2 e,+    -- Source (R r) Ix1 e,     Ord e,     Unbox e   ) =>@@ -934,16 +943,17 @@   where     Sz (m :. n) = size distMat --- Make a search row for distances. Takes row x from a distance matrix and zips them with their+-- | 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,-    Mutable r Ix2 e,-    Mutable r Ix1 (e, Ix1)+    Mutable r e,+    Mutable r (e, Ix1)   ) =>   Ix1 ->   IntSet ->