diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,10 @@ # Revision history for trexio-hs -## 0.1.0 -- YYYY-mm-dd +## 0.1.1 -- 2025-07-10 + +* Adapt test suite for more restrictive checks for various fields such as determinant_list +* Fix memory corruption bug with determinant IO + +## 0.1.0 -- 2025-01-06 * First version. Released on an unsuspecting world. diff --git a/src-int/TREXIO/Internal/Base.hsc b/src-int/TREXIO/Internal/Base.hsc --- a/src-int/TREXIO/Internal/Base.hsc +++ b/src-int/TREXIO/Internal/Base.hsc @@ -54,6 +54,14 @@ | InvalidState | VersionParsingIssue | PhaseChange + | InvalidMoIndex + | InvalidArg9 + | InvalidArg10 + | InvalidArg11 + | InvalidArg12 + | InvalidArg13 + | InvalidArg14 + | CorruptionAttempt deriving (Show, Eq, Ord, Generic) instance Enum ExitCode where @@ -96,6 +104,14 @@ fromEnum InvalidState = #const TREXIO_INVALID_STATE fromEnum VersionParsingIssue = #const TREXIO_VERSION_PARSING_ISSUE fromEnum PhaseChange = #const TREXIO_PHASE_CHANGE + fromEnum InvalidMoIndex = #const TREXIO_INVALID_MO_INDEX + fromEnum InvalidArg9 = #const TREXIO_INVALID_ARG_9 + fromEnum InvalidArg10 = #const TREXIO_INVALID_ARG_10 + fromEnum InvalidArg11 = #const TREXIO_INVALID_ARG_11 + fromEnum InvalidArg12 = #const TREXIO_INVALID_ARG_12 + fromEnum InvalidArg13 = #const TREXIO_INVALID_ARG_13 + fromEnum InvalidArg14 = #const TREXIO_INVALID_ARG_14 + fromEnum CorruptionAttempt = #const TREXIO_CORRUPTION_ATTEMPT toEnum (#const TREXIO_FAILURE) = Failure toEnum (#const TREXIO_SUCCESS) = Success @@ -136,6 +152,14 @@ toEnum (#const TREXIO_INVALID_STATE) = InvalidState toEnum (#const TREXIO_VERSION_PARSING_ISSUE) = VersionParsingIssue toEnum (#const TREXIO_PHASE_CHANGE) = PhaseChange + toEnum (#const TREXIO_INVALID_MO_INDEX) = InvalidMoIndex + toEnum (#const TREXIO_INVALID_ARG_9) = InvalidArg9 + toEnum (#const TREXIO_INVALID_ARG_10) = InvalidArg10 + toEnum (#const TREXIO_INVALID_ARG_11) = InvalidArg11 + toEnum (#const TREXIO_INVALID_ARG_12) = InvalidArg12 + toEnum (#const TREXIO_INVALID_ARG_13) = InvalidArg13 + toEnum (#const TREXIO_INVALID_ARG_14) = InvalidArg14 + toEnum (#const TREXIO_CORRUPTION_ATTEMPT) = CorruptionAttempt toEnum _ = error "toEnum(ExitCode): invalid argument" instance Exception ExitCode where @@ -234,22 +258,16 @@ close :: (MonadIO m, MonadThrow m) => Trexio -> m () close trexio = checkEC $ close_ trexio -foreign import capi "trexio.h trexio_get_int64_num" intsPerDet_ :: - Trexio -> - Ptr CInt -> - IO ExitCodeC - --- | Get the number of 'Int64's required to store a single determinant. This is a --- convenience function that avoids manual calculation via the size if Int64 and --- the number of MOs. -intsPerDet :: (MonadIO m, MonadThrow m) => Trexio -> m Int -intsPerDet trexio = liftIO . alloca $ \nPtr -> do - checkEC $ intsPerDet_ trexio nPtr - fromIntegral <$> peek nPtr - foreign import capi "trexio.h trexio_mark_safety" markSafety_ :: Trexio -> Int32 -> IO ExitCodeC -- | After a file has been opened unsafely, it can be marked as safe by this function -- forcefully. You are responsible for internal incosistencies that may arise from this. markSafety :: (MonadIO m, MonadThrow m) => Trexio -> m () markSafety trexio = checkEC $ markSafety_ trexio 0 + +-------------------------------------------------------------------------------- + +-- | Bit field type for storing determinant information. +newtype BitFieldT = BitFieldT Int64 + deriving (Eq, Ord, Show) + deriving Storable via (Int64) diff --git a/src/TREXIO.hs b/src/TREXIO.hs --- a/src/TREXIO.hs +++ b/src/TREXIO.hs @@ -78,7 +78,6 @@ -- * High Level Interface scheme, - intsPerDet, withTrexio, module TREXIO.HighLevel, ) where diff --git a/src/TREXIO/HighLevel.hs b/src/TREXIO/HighLevel.hs --- a/src/TREXIO/HighLevel.hs +++ b/src/TREXIO/HighLevel.hs @@ -18,6 +18,9 @@ import Data.Massiv.Array as Massiv hiding (Dim, forM) import Foreign.C.ConstPtr import Foreign.C.Types +import Foreign.Marshal.Alloc +import Foreign.Ptr +import Foreign.Storable import Language.Haskell.TH import TREXIO.Internal.Base import TREXIO.Internal.TH @@ -38,3 +41,12 @@ return $ hasBind <> readBind <> writeBind return $ deleteBind <> fieldBinds ) + +{- | Get the number of 'Int64's required to store a single determinant. This is a +convenience function that avoids manual calculation via the size if Int64 and +the number of MOs. +-} +intsPerDet :: (MonadIO m, MonadThrow m) => Trexio -> m Int +intsPerDet trexio = liftIO . alloca $ \nPtr -> do + checkEC $ trexio_get_int64_num trexio nPtr + fromIntegral <$> peek nPtr diff --git a/src/TREXIO/Internal/Marshaller.hs b/src/TREXIO/Internal/Marshaller.hs --- a/src/TREXIO/Internal/Marshaller.hs +++ b/src/TREXIO/Internal/Marshaller.hs @@ -1,5 +1,6 @@ module TREXIO.Internal.Marshaller where +import Control.Exception import Data.Massiv.Array as Massiv hiding (withMArray) import Data.Massiv.Array.Unsafe import Data.Maybe (fromJust) @@ -22,8 +23,9 @@ mArr <- thaw v withMArray mArr f --- | Get an 'MArray' from C memory. Haskell and C side use the same memory reference. --- Be careful with this function. When the undelying pointer vanishes, the array is nonsense +{- | Get an 'MArray' from C memory. Haskell and C side use the same memory reference. +Be careful with this function. When the undelying pointer vanishes, the array is nonsense +-} unsafeToMArray :: (Index ix, Storable a) => Sz ix -> Ptr a -> IO (MArray s S ix a) unsafeToMArray sz ptr = do fPtr <- newForeignPtr_ ptr @@ -174,3 +176,10 @@ , fromIntegral t , fromIntegral s ] + +-- | Like 'callocArray' but safe with bracket pattern +callocaArray :: (Storable a) => Int -> (Ptr a -> IO b) -> IO b +callocaArray nEl = + bracket + (callocArray nEl) + free diff --git a/src/TREXIO/Internal/TH.hs b/src/TREXIO/Internal/TH.hs --- a/src/TREXIO/Internal/TH.hs +++ b/src/TREXIO/Internal/TH.hs @@ -28,6 +28,7 @@ import Data.Text.IO qualified as T import Data.Vector qualified as V import Foreign hiding (peekArray, void, withArray) +import Foreign qualified as F import Foreign.C.ConstPtr import Foreign.C.String import Foreign.C.Types @@ -946,31 +947,29 @@ let Sz2 nDets _nMos = size dets $(mkWriteSzFn scheme d1) trexio nDets - allocaArray (nDets * nInt64PerDet * 2) $ \detBuf -> do + callocaArray (nDets * nInt64PerDet * 2) $ \(detBuf :: Ptr Int64) -> do -- Write each determinant to the buffer forM_ [0 .. nDets - 1] $ \i -> do + -- Get a single determinant (up and down spin components) and + -- convert to a Storable Vector of Word8 let det = dets !> i - detToByteString bv accFn = - BV.cloneToByteString - . Massiv.toVector - . compute @U - . Massiv.map accFn - $ bv - up = detToByteString det fst - down = detToByteString det snd - nBytes = BS.length up upPtr = detBuf `plusPtr` (i * nInt64PerDet * 2 * sizeOf (undefined :: Int64)) - downPtr = upPtr `plusPtr` (nInt64PerDet * sizeOf (undefined :: Int64)) + dnPtr = upPtr `plusPtr` (nInt64PerDet * sizeOf (undefined :: Int64)) - -- Up spin - BS.unsafeUseAsCString up $ \charPtr -> do - copyBytes (castPtr upPtr) charPtr nBytes + let toOrbList tix = ifoldlS (\acc idx b -> if b == 1 then idx : acc else acc) mempty (Massiv.map tix det) + orbUp = toOrbList fst + orbDn = toOrbList snd + nOccUp = length orbUp + nOccDn = length orbDn - -- Down spin - BS.unsafeUseAsCString down $ \charPtr -> do - copyBytes (castPtr downPtr) charPtr nBytes + -- Write the Bitfields to corresponding parts of the determinant buffer + F.withArray (fromIntegral <$> orbUp) $ \upListPtr -> do + trexio_to_bitfield_list (ConstPtr upListPtr) (fromIntegral nOccUp) (ConstPtr upPtr) (fromIntegral nInt64PerDet) - -- Call the C funciton with the buffer + F.withArray (fromIntegral <$> orbDn) $ \dnListPtr -> do + trexio_to_bitfield_list (ConstPtr dnListPtr) (fromIntegral nOccDn) (ConstPtr dnPtr) (fromIntegral nInt64PerDet) + + -- Call the C function with the buffer checkEC $ $(varE . mkName $ mkCFnName Write groupName dataName) trexio diff --git a/src/TREXIO/LowLevel.hs b/src/TREXIO/LowLevel.hs --- a/src/TREXIO/LowLevel.hs +++ b/src/TREXIO/LowLevel.hs @@ -16,9 +16,11 @@ -} module TREXIO.LowLevel where +import Data.Int import Data.Map qualified as Map import Foreign.C.ConstPtr import Foreign.C.Types +import Foreign.Ptr import TREXIO.Internal.Base import TREXIO.Internal.TH import TREXIO.LowLevel.Scheme @@ -30,3 +32,23 @@ -- Import all C functions for all fields and operations concat <$> traverse (uncurry mkCBindings) groups ) + +-- | 'Int64's required per determinant +foreign import capi "trexio.h trexio_get_int64_num" + trexio_get_int64_num :: + Trexio -> + Ptr CInt -> + IO ExitCodeC + +-- | Take a list of occupied orbital indices and create a bitfield from it +foreign import capi "trexio.h trexio_to_bitfield_list" + trexio_to_bitfield_list :: + -- | @orb_list@ + ConstPtr Int32 -> + -- | @occupied_num@ + Int32 -> + -- | Preallocated, zeroed bit array representing a determinant + ConstPtr Int64 -> + -- | Number of Int64 required per determinant as obtained by 'trexio_get_int64_num' + Int32 -> + IO ExitCodeC diff --git a/test/trexio-test.hs b/test/trexio-test.hs --- a/test/trexio-test.hs +++ b/test/trexio-test.hs @@ -1,9 +1,13 @@ +import Control.Concurrent (threadDelay) import Control.Exception.Safe import Control.Monad +import Data.Bit.ThreadSafe (Bit) import Data.Massiv.Array as Massiv hiding (Size, elem, forM, forM_, mapM, mapM_, take, zip, zipWith) +import Data.Massiv.Array qualified as Massiv import Data.Maybe (catMaybes, fromJust) import Data.Set qualified as Set import Data.Text (Text) +import Debug.Trace import Hedgehog (Gen, MonadGen, Size, forAll, property, (===)) import Hedgehog.Gen qualified as Gen import Hedgehog.Range qualified as Range @@ -26,11 +30,11 @@ "0D" [ testGroup "Integers" . appFn $ [ ("nucleus.num", genDim, deleteNucleus, hasNucleusNum, readNucleusNum, writeNucleusNum) - , ("grid.max_ang_num", genInt, deleteGrid, hasGridMaxAngNum, readGridMaxAngNum, writeGridMaxAngNum) + , ("grid.max_ang_num", genPosInt, deleteGrid, hasGridMaxAngNum, readGridMaxAngNum, writeGridMaxAngNum) , ("state.id", genIndex, deleteState, hasStateId, readStateId, writeStateId) ] , testGroup "Floats" . appFn $ - [ ("nucleus.repulsion", genFloat, deleteNucleus, hasNucleusRepulsion, readNucleusRepulsion, writeNucleusRepulsion) + [ ("nucleus.repulsion", genPosFloat, deleteNucleus, hasNucleusRepulsion, readNucleusRepulsion, writeNucleusRepulsion) ] , testGroup "Strings" . appFn $ [ ("metadata.description", genIdentifier, deleteMetadata, hasMetadataDescription, readMetadataDescription, writeMetadataDescription) @@ -47,35 +51,52 @@ , testGroup "Strings" . appFn $ [ ("metadata.author", genVector genIdentifier, deleteMetadata, hasMetadataAuthor, readMetadataAuthor, writeMetadataAuthor) ] - , testCase "Determinant IO" . withSystemTempFile "trexio.dat" $ \fp _ -> - withTrexio fp FileWrite Hdf5 $ \trexio -> do - writeMoNum trexio 3 + , testProperty "Determinant IO" . property $ do + -- Generate occupation numbers + nMo <- forAll $ Gen.int (Range.linear 1 1000) + nUp <- forAll $ Gen.int (Range.linear 2 (nMo - 1)) + nDn <- forAll $ Gen.int (Range.linear 2 (nMo - 1)) - has1 <- hasDeterminantList trexio - has1 @?= False + -- Generate determinants + dets <- forAll $ genDet nMo (nUp, nDn) + let Sz2 nDets _ = Massiv.size dets - ingoreExcp [AttrMissing] (readDeterminantList trexio) + -- Generate coefficients + let + coeffs <- forAll . fmap (Massiv.fromList Par) $ Gen.list (Range.singleton nDets) (Gen.double $ Range.linearFrac (-10) 10) - let detList = - Massiv.fromLists' - Seq - [ [(1, 1), (1, 0), (0, 0)] - , [(1, 0), (1, 1), (0, 0)] - ] - detCoeffs = Massiv.fromList Seq [1.0, 2.0] + liftIO $ withSystemTempFile "trexio.dat" $ \fp _ -> + withTrexio fp FileWrite Hdf5 $ \trexio -> do + -- Write MO numbers and electron occupation numbers to the file + writeMoNum trexio nMo + writeElectronUpNum trexio nUp + writeElectronDnNum trexio nDn - writeDeterminantList trexio detList + nMo' <- readMoNum trexio + nUp' <- readElectronUpNum trexio + nDn' <- readElectronDnNum trexio - has2 <- hasDeterminantList trexio - has2 @?= True + -- Check + nMo @?= nMo' + nUp @?= nUp' + nDn @?= nDn' - detList' <- readDeterminantList trexio - detList' @?= detList + -- Write determinants + writeDeterminantList trexio dets - writeDeterminantCoefficient trexio detCoeffs + -- Read back determinants + readDets <- readDeterminantList trexio - detCoeffs' <- readDeterminantCoefficient trexio - detCoeffs' @?= detCoeffs + -- Check for equality + dets @?= readDets + + -- Write coefficients + writeDeterminantCoefficient trexio coeffs + + -- Read back coefficients + coeffs' <- readDeterminantCoefficient trexio + + coeffs @?= coeffs' ] , testGroup "2D" @@ -181,6 +202,9 @@ genInt :: Gen Int genInt = Gen.integral (Range.linearFrom 0 (-1_000_000) 1_000_000) +genPosInt :: Gen Int +genPosInt = Gen.integral (Range.linearFrom 0 0 100) + -- | Generate a @dim@ value, which is a positive integer genDim :: Gen Int genDim = Gen.integral (Range.linear 1 100) @@ -189,9 +213,15 @@ genIndex :: Gen Int genIndex = Gen.integral (Range.linear 0 1000) +genMoIndex :: Gen Word +genMoIndex = Gen.word (Range.linear 0 300) + genFloat :: Gen Double genFloat = Gen.realFloat (Range.linearFrac (-1_000_000) 1_000_000) +genPosFloat :: Gen Double +genPosFloat = Gen.realFloat (Range.linearFrac 0 1_000_000) + -- | Generate a identifier, that is a single word without spaces or stuff genIdentifier :: Gen Text genIdentifier = Gen.text (Range.linear 1 10) Gen.alphaNum @@ -444,3 +474,34 @@ where sz2zs :: Size -> Size sz2zs x = round $ fromIntegral x * (0.25 :: Double) + +-- | Generate multiple determinants valid for given system +genDet :: + -- | Number of MOs in the system + Int -> + -- | Number of Up and Down electrons + (Int, Int) -> + -- | List of determinants + Gen (Matrix U (Bit, Bit)) +genDet nMo (nUp, nDn) = do + detsL <- Gen.set (Range.linear 1 100) detGen + dets <- case Massiv.stackOuterSlicesM . Set.toList $ detsL of + Nothing -> error "Failed to stack slices" + Just dets' -> return dets' + return . compute $ dets + where + occGen :: (MonadGen m) => Int -> m (Set.Set Int) + occGen nOcc = Gen.set (Range.singleton nOcc) (Gen.int (Range.linear 0 (nMo - 1))) + + detGen :: (MonadGen m) => m (Massiv.Vector U (Bit, Bit)) + detGen = do + -- Generate indices of occupied orbitals + occUp <- occGen nUp + occDn <- occGen nDn + + -- Generate a single determinant + let detUp = makeArray @U Par (Sz nMo) $ \i -> if i `Set.member` occUp then 1 else 0 + detDn = makeArray @U Par (Sz nMo) $ \i -> if i `Set.member` occDn then 1 else 0 + det = Massiv.zip detUp detDn + + return . compute $ det diff --git a/trexio-hs.cabal b/trexio-hs.cabal --- a/trexio-hs.cabal +++ b/trexio-hs.cabal @@ -7,7 +7,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 0.1.0 +version: 0.1.1 synopsis: Bindings to the TREXIO library for wave function data homepage: https://github.com/TREX-CoE/trexio-hs license: BSD-3-Clause @@ -16,7 +16,7 @@ category: Data build-type: Simple extra-doc-files: CHANGELOG.md -tested-with: GHC == {9.6, 9.8, 9.10} +tested-with: GHC == {9.6, 9.8, 9.10, 9.12} description: This package provides low- and high-level Haskell bindings for [TREXIO, a portable file format for storing wave function data](https://trex-coe.github.io/trexio/). The vast majority of the bindings in this package is generated via TemplateHaskell from the TREXIO JSON specification, that then defines the C-API. @@ -75,7 +75,7 @@ bitvec >= 1.1.5.0 && < 1.2, bytestring >= 0.10 && < 0.13, casing >= 0.1.4 && < 0.2, - containers >= 0.6 && < 0.8, + containers >= 0.6 && < 0.9, filepath >= 1.4 && < 1.6, massiv >= 1.0.0.0 && < 1.1, safe-exceptions >= 0.1.7 && < 0.2,