diff options
Diffstat (limited to 'contrib/haskell/src')
58 files changed, 6196 insertions, 0 deletions
diff --git a/contrib/haskell/src/Hkl.hs b/contrib/haskell/src/Hkl.hs new file mode 100644 index 0000000..d52a69a --- /dev/null +++ b/contrib/haskell/src/Hkl.hs @@ -0,0 +1,16 @@ +module Hkl (module X) where + +import Hkl.C as X +import Hkl.DataSource as X +import Hkl.Detector as X +import Hkl.Engine as X +import Hkl.Flat as X +import Hkl.H5 as X +import Hkl.Lattice as X +import Hkl.MyMatrix as X +import Hkl.Nxs as X +import Hkl.PyFAI as X +import Hkl.Script as X +import Hkl.Tiff as X +import Hkl.Types as X +import Hkl.Xrd as X diff --git a/contrib/haskell/src/Hkl/C.hsc b/contrib/haskell/src/Hkl/C.hsc new file mode 100644 index 0000000..6066d51 --- /dev/null +++ b/contrib/haskell/src/Hkl/C.hsc @@ -0,0 +1,160 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE CPP #-} + +module Hkl.C + ( compute + , computePipe + , solve + , solveTraj + , solveTrajPipe + , module X + ) where + +import Prelude hiding (min, max) + +import Control.Monad (forever) +import Control.Monad.Trans.State.Strict +import Foreign ( ForeignPtr + , FunPtr + , Ptr + , nullPtr + , newForeignPtr + , withForeignPtr + , withArray) +import Foreign.C ( CInt(..), CSize(..), CString + , withCString) + +import Pipes (Pipe, await, lift, yield) + +import Hkl.C.Detector +import Hkl.C.Engine +import Hkl.C.EngineList +import Hkl.C.Geometry as X +import Hkl.C.GeometryList as X +import Hkl.C.Sample +import Hkl.Detector +import Hkl.Types + +#include "hkl.h" + +-- Engine + +solve' :: Ptr HklEngine -> Engine -> IO (ForeignPtr HklGeometryList) +solve' engine (Engine _ ps _) = do + let positions = [v | (Parameter _ v _) <- ps] + let n = toEnum (length positions) + withArray positions $ \values -> + c_hkl_engine_pseudo_axis_values_set engine values n unit nullPtr + >>= newForeignPtr c_hkl_geometry_list_free + +solve :: Geometry -> Detector a -> Sample b -> Engine -> IO [Geometry] +solve g@(Geometry f _ _ _) d s e@(Engine name _ _) = do + withSample s $ \sample -> + withDetector d $ \detector -> + withGeometry g $ \geometry -> + withEngineList f $ \engines -> + withCString name $ \cname -> do + c_hkl_engine_list_init engines geometry detector sample + engine <- c_hkl_engine_list_engine_get_by_name engines cname nullPtr + solve' engine e >>= peekHklGeometryList + +solveTraj :: Geometry -> Detector a -> Sample b -> [Engine] -> IO [Geometry] +solveTraj g@(Geometry f _ _ _) d s es = do + let name = engineName (head es) + withSample s $ \sample -> + withDetector d $ \detector -> + withGeometry g $ \geometry -> + withEngineList f $ \engines -> + withCString name $ \cname -> do + c_hkl_engine_list_init engines geometry detector sample + engine <- c_hkl_engine_list_engine_get_by_name engines cname nullPtr + mapM (\e -> solve' engine e >>= getSolution0) es + +-- Pipe + +data Diffractometer = Diffractometer { difEngineList :: (ForeignPtr HklEngineList) + , difGeometry :: (ForeignPtr Geometry) + , difDetector :: (ForeignPtr HklDetector) + , difSample :: (ForeignPtr HklSample) + } + deriving (Show) + +withDiffractometer :: Diffractometer -> (Ptr HklEngineList -> IO b) -> IO b +withDiffractometer d fun = do + let f_engines = difEngineList d + withForeignPtr f_engines fun + +newDiffractometer :: Geometry -> Detector a -> Sample b -> IO Diffractometer +newDiffractometer g@(Geometry f _ _ _) d s = do + f_engines <- newEngineList f + f_geometry <- newGeometry g + f_detector <- newDetector d + f_sample <- newSample s + withForeignPtr f_sample $ \sample -> + withForeignPtr f_detector $ \detector -> + withForeignPtr f_geometry $ \geometry -> + withForeignPtr f_engines $ \engines -> do + c_hkl_engine_list_init engines geometry detector sample + return $ Diffractometer { difEngineList = f_engines + , difGeometry = f_geometry + , difDetector = f_detector + , difSample = f_sample + } + +computePipe :: Detector a -> Sample b -> Pipe Geometry [Engine] IO () +computePipe d s = forever $ do + g <- await + e <- lift $ compute g d s + yield e + +solveTrajPipe :: Geometry -> Detector a -> Sample b -> Pipe Engine Geometry IO () +solveTrajPipe g d s = do + dif <- lift $ newDiffractometer g d s + solveTrajPipe' dif + +solveTrajPipe' :: Diffractometer -> Pipe Engine Geometry IO () +solveTrajPipe' dif = flip evalStateT dif $ forever $ do + -- Inside here we are using `StateT Diffractometer (Pipe Engine Geometry IO ()) r` + e <- lift await + dif_ <- get + let name = engineName e + solutions <- lift . lift $ withDiffractometer dif_ $ \engines -> + withCString name $ \cname -> do + engine <- c_hkl_engine_list_engine_get_by_name engines cname nullPtr + solve' engine e >>= getSolution0 + put dif_ + lift $ yield solutions + +foreign import ccall unsafe "hkl.h hkl_engine_list_engine_get_by_name" + c_hkl_engine_list_engine_get_by_name :: Ptr HklEngineList --engine list + -> CString -- engine name + -> Ptr () -- gerror need to deal about this + -> IO (Ptr HklEngine) -- the returned HklEngine + +foreign import ccall unsafe "hkl.h hkl_engine_pseudo_axis_values_set" + c_hkl_engine_pseudo_axis_values_set :: Ptr HklEngine + -> Ptr Double --values + -> CSize -- n_values + -> CInt -- unit_type + -> Ptr () -- GError **error + -> IO (Ptr HklGeometryList) + +foreign import ccall unsafe "hkl.h &hkl_geometry_list_free" + c_hkl_geometry_list_free :: FunPtr (Ptr HklGeometryList -> IO ()) + +compute :: Geometry -> Detector a -> Sample b -> IO [Engine] +compute g@(Geometry f _ _ _) d s = do + withSample s $ \sample -> + withDetector d $ \detector -> + withGeometry g $ \geometry -> + withEngineList f $ \engines -> do + c_hkl_engine_list_init engines geometry detector sample + c_hkl_engine_list_get engines + engineListEnginesGet engines + +foreign import ccall unsafe "hkl.h hkl_engine_list_init" + c_hkl_engine_list_init:: Ptr HklEngineList -> Ptr Geometry -> Ptr HklDetector -> Ptr HklSample -> IO () + +foreign import ccall unsafe "hkl.h hkl_engine_list_get" + c_hkl_engine_list_get:: Ptr HklEngineList -> IO () diff --git a/contrib/haskell/src/Hkl/C/DArray.hsc b/contrib/haskell/src/Hkl/C/DArray.hsc new file mode 100644 index 0000000..82520ee --- /dev/null +++ b/contrib/haskell/src/Hkl/C/DArray.hsc @@ -0,0 +1,25 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE ForeignFunctionInterface #-} + +module Hkl.C.DArray + (DArray(..)) where + +import Foreign (peekArray) +import Foreign.C (CSize, CString) +import Foreign.Storable (Storable(..)) + +#include "hkl.h" + +#let alignment t = "%lu", (unsigned long)offsetof(struct {char x__; t (y__); }, y__) + +data DArray a = DArray CSize [a] deriving Show + +instance Storable (DArray CString) where + alignment _ = #{alignment darray_string} + sizeOf _ = #{size darray_string} + peek ptr = do + n <- (#{peek darray_string, size} ptr) + items <- #{peek darray_string ,item} ptr + ss <- peekArray (fromEnum n) items + return $ DArray n ss diff --git a/contrib/haskell/src/Hkl/C/Detector.hsc b/contrib/haskell/src/Hkl/C/Detector.hsc new file mode 100644 index 0000000..73c6b1d --- /dev/null +++ b/contrib/haskell/src/Hkl/C/Detector.hsc @@ -0,0 +1,41 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE GADTs #-} + +module Hkl.C.Detector + ( HklDetector + , newDetector + , withDetector + ) where + +import Prelude hiding (min, max) + +import Foreign ( ForeignPtr + , FunPtr + , Ptr + , newForeignPtr + , withForeignPtr) +import Foreign.C (CInt(..)) + +import Hkl.Detector + +#include "hkl.h" + +data HklDetector + +-- Detector + +withDetector :: Detector a -> (Ptr HklDetector -> IO b) -> IO b +withDetector d func = do + fptr <- newDetector d + withForeignPtr fptr func + +newDetector :: Detector a -> IO (ForeignPtr HklDetector) +newDetector ZeroD = c_hkl_detector_new 0 >>= newForeignPtr c_hkl_detector_free +newDetector _ = error "Can not use 2D detector with the hkl library" + +foreign import ccall unsafe "hkl.h hkl_detector_new" + c_hkl_detector_new:: CInt -> IO (Ptr HklDetector) + +foreign import ccall unsafe "hkl.h &hkl_detector_free" + c_hkl_detector_free :: FunPtr (Ptr HklDetector -> IO ()) diff --git a/contrib/haskell/src/Hkl/C/Engine.hsc b/contrib/haskell/src/Hkl/C/Engine.hsc new file mode 100644 index 0000000..9d5eced --- /dev/null +++ b/contrib/haskell/src/Hkl/C/Engine.hsc @@ -0,0 +1,81 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE GADTs #-} + +module Hkl.C.Engine + ( HklEngine + , engineName + , peekEngine + ) where + +import Prelude hiding (min, max) + +import Foreign (Ptr, nullPtr) +import Foreign.C (CString, peekCString) +import Foreign.Storable + +import Hkl.C.DArray +import Hkl.Types + +#include "hkl.h" + +-- private types + +data HklEngine + +-- Engine + +engineName :: Engine -> String +engineName (Engine name _ _) = name + +-- Engine + +peekMode :: Ptr HklEngine -> IO Mode +peekMode e = do + name <- c_hkl_engine_current_mode_get e >>= peekCString + (DArray _ ns) <- peek =<< c_hkl_engine_parameters_names_get e + parameters <- mapM f ns + return (Mode name parameters) + where + f n = (c_hkl_engine_parameter_get e n nullPtr >>= peek) + +foreign import ccall unsafe "hkl.h hkl_engine_current_mode_get" + c_hkl_engine_current_mode_get :: Ptr HklEngine -> IO CString + +foreign import ccall unsafe "hkl.h hkl_engine_parameters_names_get" + c_hkl_engine_parameters_names_get:: Ptr HklEngine -> IO (Ptr (DArray CString)) + +foreign import ccall unsafe "hkl.h hkl_engine_parameter_get" + c_hkl_engine_parameter_get:: Ptr HklEngine -> CString -> Ptr () -> IO (Ptr Parameter) -- darray_string + + +peekEngine :: Ptr HklEngine -> IO Engine +peekEngine e = do + name <- peekCString =<< c_hkl_engine_name_get e + ps <- enginePseudoAxesGet e + mode <- peekMode e + return (Engine name ps mode) + +-- engineNameGet :: Ptr HklEngine -> IO String +-- engineNameGet engine = c_hkl_engine_name_get engine >>= peekCString + +foreign import ccall unsafe "hkl.h hkl_engine_name_get" + c_hkl_engine_name_get :: Ptr HklEngine -> IO CString + +foreign import ccall unsafe "hkl.h hkl_engine_pseudo_axis_names_get" + c_hkl_engine_pseudo_axis_names_get:: Ptr HklEngine -> IO (Ptr (DArray CString)) + +-- enginePseudoAxisNamesGet :: Ptr HklEngine -> IO [String] +-- enginePseudoAxisNamesGet e = enginePseudoAxisNamesGet' e >>= mapM peekCString + +enginePseudoAxisGet :: Ptr HklEngine -> CString -> IO Parameter +enginePseudoAxisGet e n = c_hkl_engine_pseudo_axis_get e n nullPtr >>= peek + +foreign import ccall unsafe "hkl.h hkl_engine_pseudo_axis_get" + c_hkl_engine_pseudo_axis_get:: Ptr HklEngine -> CString -> Ptr () -> IO (Ptr Parameter) + +enginePseudoAxesGet :: Ptr HklEngine -> IO [Parameter] +enginePseudoAxesGet ptr = do + (DArray _ ns) <- peek =<< c_hkl_engine_pseudo_axis_names_get ptr + mapM (enginePseudoAxisGet ptr) ns + diff --git a/contrib/haskell/src/Hkl/C/EngineList.hsc b/contrib/haskell/src/Hkl/C/EngineList.hsc new file mode 100644 index 0000000..08232f6 --- /dev/null +++ b/contrib/haskell/src/Hkl/C/EngineList.hsc @@ -0,0 +1,60 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE GADTs #-} + +module Hkl.C.EngineList + ( HklEngineList + , engineListEnginesGet + , newEngineList + , withEngineList + ) where + +import Prelude hiding (min, max) + +import Foreign ( ForeignPtr + , FunPtr + , Ptr + , newForeignPtr + , withForeignPtr + , peekArray) +import Foreign.C ( CSize(..) ) +import Foreign.Storable + +import Hkl.C.Engine +import Hkl.C.Geometry +import Hkl.Types + +#include "hkl.h" + +-- private types + +data HklEngineList + +-- EngineList + +withEngineList :: Factory -> (Ptr HklEngineList -> IO b) -> IO b +withEngineList f func = do + fptr <- newEngineList f + withForeignPtr fptr func + +newEngineList :: Factory -> IO (ForeignPtr HklEngineList) +newEngineList f = newFactory f + >>= c_hkl_factory_create_new_engine_list + >>= newForeignPtr c_hkl_engine_list_free + +foreign import ccall unsafe "hkl.h hkl_factory_create_new_engine_list" + c_hkl_factory_create_new_engine_list:: Ptr HklFactory -> IO (Ptr HklEngineList) + +foreign import ccall unsafe "hkl.h &hkl_engine_list_free" + c_hkl_engine_list_free :: FunPtr (Ptr HklEngineList -> IO ()) + +engineListEnginesGet :: Ptr HklEngineList -> IO [Engine] +engineListEnginesGet e = do + pdarray <- c_hkl_engine_list_engines_get e + n <- (#{peek darray_engine, size} pdarray) :: IO CSize + engines <- #{peek darray_engine ,item} pdarray :: IO (Ptr (Ptr HklEngine)) + enginess <- peekArray (fromEnum n) engines + mapM peekEngine enginess + +foreign import ccall unsafe "hkl.h hkl_engine_list_engines_get" + c_hkl_engine_list_engines_get:: Ptr HklEngineList -> IO (Ptr ()) diff --git a/contrib/haskell/src/Hkl/C/Geometry.hsc b/contrib/haskell/src/Hkl/C/Geometry.hsc new file mode 100644 index 0000000..406c65d --- /dev/null +++ b/contrib/haskell/src/Hkl/C/Geometry.hsc @@ -0,0 +1,188 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE GADTs #-} + +module Hkl.C.Geometry + ( Geometry(..) + , Factory(..) + , HklFactory + , HklMatrix + , HklQuaternion + , factoryFromString + , newFactory + , newGeometry + , withGeometry + ) where + +import Prelude hiding (min, max) + +import Numeric.LinearAlgebra +import Foreign ( ForeignPtr + , FunPtr + , Ptr + , nullPtr + , newForeignPtr + , withForeignPtr) +import Foreign.C (CInt(..), CDouble(..), CSize(..), CString, + peekCString, withCString) +import Foreign.Storable + +import Numeric.Units.Dimensional.Prelude ( meter, nano + , (*~), (/~)) + +import qualified Data.Vector.Storable as V +import qualified Data.Vector.Storable.Mutable as MV + +import Hkl.Types +import Hkl.C.DArray + +#include "hkl.h" + +-- | Factory + +data Factory = K6c | Uhv | MedH | MedV | SoleilSiriusKappa + +instance Show Factory where + show K6c = "K6C" + show Uhv = "ZAXIS" + show MedH = "todo" + show MedV = "todo" + show SoleilSiriusKappa = "SOLEIL SIRIUS KAPPA" + +factoryFromString :: String -> Factory +factoryFromString s + | s == "K6C" = K6c + | s == "ZAXIS" = Uhv + | s == "todo" = MedH + | s == "todo" = MedV + | s == "SOLEIL SIRIUS KAPPA" = SoleilSiriusKappa + | otherwise = error $ "unknown diffractometer type:" ++ s + +-- | Geometry + +data Geometry = Geometry + Factory -- ^ the type of diffractometer + Source -- ^ source + (Vector Double) -- ^ axes position + (Maybe [Parameter]) -- ^ axes configuration + deriving (Show) + + +-- private types + +data HklFactory +data HklMatrix +data HklQuaternion + +#if __GLASGOW_HASKELL__ <= 710 +#let alignment t = "%lu", (unsigned long)offsetof(struct {char x__; t (y__); }, y__) +#endif + +-- Factory + +newFactory :: Factory -> IO (Ptr HklFactory) +newFactory f = withCString (show f) $ \cname -> c_hkl_factory_get_by_name cname nullPtr + +foreign import ccall unsafe "hkl.h hkl_factory_get_by_name" + c_hkl_factory_get_by_name :: CString -- ^ name + -> Ptr () -- ^ GError (null for now) + -> IO (Ptr HklFactory) +-- Geometry + +peekSource :: Ptr Geometry -> IO (Source) +peekSource ptr = do + (CDouble w) <- c_hkl_geometry_wavelength_get ptr unit + return (Source (w *~ nano meter)) + +foreign import ccall unsafe "hkl.h hkl_geometry_wavelength_set" + c_hkl_geometry_wavelength_set :: Ptr Geometry -- geometry + -> CDouble -- wavelength + -> CInt -- unit + -> Ptr () -- *gerror + -> IO () -- IO CInt but for now do not deal with the errors + +pokeSource :: Ptr Geometry -> Source -> IO () +pokeSource ptr (Source lw) = do + let wavelength = CDouble (lw /~ nano meter) + c_hkl_geometry_wavelength_set ptr wavelength unit nullPtr + +foreign import ccall unsafe "hkl.h hkl_geometry_wavelength_get" + c_hkl_geometry_wavelength_get :: Ptr Geometry -- geometry + -> CInt -- unit + -> IO CDouble -- wavelength + +peekAxis :: Ptr Geometry -> CString -> IO Parameter +peekAxis ptr s = c_hkl_geometry_axis_get ptr s nullPtr >>= peek + +instance Storable Geometry where + alignment _ = #{alignment int} + + sizeOf _ = #{size int} + + peek ptr = do + f_name <- c_hkl_geometry_name_get ptr >>= peekCString + let factory = factoryFromString f_name + + source <- peekSource ptr + + (DArray n axis_names) <- peek =<< c_hkl_geometry_axis_names_get ptr + v <- MV.new (fromEnum n) + MV.unsafeWith v $ \values -> + c_hkl_geometry_axis_values_get ptr values n unit + vs <- V.freeze v + + ps <- mapM (peekAxis ptr) axis_names + + return $ Geometry factory source vs (Just ps) + + poke ptr (Geometry _ s vs _) = do + pokeSource ptr s + (DArray n _) <- peek =<< c_hkl_geometry_axis_names_get ptr + V.unsafeWith vs $ \values -> + c_hkl_geometry_axis_values_set ptr values n unit nullPtr + +foreign import ccall unsafe "hkl.h hkl_geometry_axis_values_get" + c_hkl_geometry_axis_values_get :: Ptr Geometry -- geometry + -> Ptr Double -- axis values + -> CSize -- size of axis values + -> CInt -- unit + -> IO () -- IO CInt but for now do not deal with the errors + +foreign import ccall unsafe "hkl.h hkl_geometry_axis_names_get" + c_hkl_geometry_axis_names_get :: Ptr Geometry -- goemetry + -> IO (Ptr (DArray CString)) -- darray_string + +foreign import ccall unsafe "hkl.h hkl_geometry_axis_get" + c_hkl_geometry_axis_get :: Ptr Geometry -- geometry + -> CString -- axis name + -> Ptr () -- gerror + -> IO (Ptr Parameter) -- parameter or nullPtr + +foreign import ccall unsafe "hkl.h hkl_geometry_name_get" + c_hkl_geometry_name_get :: Ptr Geometry -> IO CString + +foreign import ccall unsafe "hkl.h hkl_geometry_axis_values_set" + c_hkl_geometry_axis_values_set :: Ptr Geometry -- geometry + -> Ptr Double -- axis values + -> CSize -- size of axis values + -> CInt -- unit + -> Ptr () -- gerror + -> IO () -- IO CInt but for now do not deal with the errors + +withGeometry :: Geometry -> (Ptr Geometry -> IO b) -> IO b +withGeometry g fun = do + fptr <- newGeometry g + withForeignPtr fptr fun + +newGeometry :: Geometry -> IO (ForeignPtr Geometry) +newGeometry g@(Geometry f _ _ _) = do + ptr <- c_hkl_factory_create_new_geometry =<< newFactory f + poke ptr g + newForeignPtr c_hkl_geometry_free ptr + +foreign import ccall unsafe "hkl.h hkl_factory_create_new_geometry" + c_hkl_factory_create_new_geometry :: Ptr HklFactory -> IO (Ptr Geometry) + +foreign import ccall unsafe "hkl.h &hkl_geometry_free" + c_hkl_geometry_free :: FunPtr (Ptr Geometry -> IO ()) diff --git a/contrib/haskell/src/Hkl/C/GeometryList.hsc b/contrib/haskell/src/Hkl/C/GeometryList.hsc new file mode 100644 index 0000000..a51067c --- /dev/null +++ b/contrib/haskell/src/Hkl/C/GeometryList.hsc @@ -0,0 +1,120 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE CPP #-} + +module Hkl.C.GeometryList + ( HklGeometryList + , geometryDetectorRotationGet + , getSolution0 + , peekHklGeometryList + ) where + +import Prelude hiding (min, max) + +import Control.Monad.Loops (unfoldrM) +import Numeric.LinearAlgebra +import Foreign ( ForeignPtr + , FunPtr + , Ptr + , nullPtr + , newForeignPtr + , withForeignPtr) +import Foreign.C (CInt(..), CDouble(..)) +import Foreign.Storable + +import Hkl.C.Detector +import Hkl.C.Geometry +import Hkl.Detector + +#include "hkl.h" + +-- private types + +data HklGeometryList +data HklGeometryListItem + +-- | HklGeometryList + +getSolution0 :: ForeignPtr HklGeometryList -> IO Geometry +getSolution0 gl = withForeignPtr gl $ \solutions -> + c_hkl_geometry_list_items_first_get solutions + >>= c_hkl_geometry_list_item_geometry_get + >>= peek + +buildMatrix' :: Element a => CInt -> CInt -> ((CInt, CInt) -> IO a) -> IO (Matrix a) +buildMatrix' rc cc f = do + let coordinates' = map (\ ri -> map (\ ci -> (ri, ci)) [0 .. (cc - 1)]) [0 .. (rc - 1)] + l <- mapM (mapM f) coordinates' + return $ fromLists l + + + -- fromLists $ map (map f) + -- $ map (\ ri -> map (\ ci -> (ri, ci)) [0 .. (cc - 1)]) [0 .. (rc - 1)] + +geometryDetectorRotationGet :: Geometry -> Detector a -> IO (Matrix Double) +geometryDetectorRotationGet g d = do + f_geometry <- newGeometry g + f_detector <- newDetector d + withForeignPtr f_detector $ \detector -> + withForeignPtr f_geometry $ \geometry -> do + f_q <- newForeignPtr c_hkl_quaternion_free =<< c_hkl_geometry_detector_rotation_get_binding geometry detector + withForeignPtr f_q $ \quaternion -> do + f_m <- newForeignPtr c_hkl_matrix_free =<< c_hkl_quaternion_to_matrix_binding quaternion + withForeignPtr f_m $ \matrix' -> + buildMatrix' 3 3 (getV matrix') + where + getV :: Ptr HklMatrix -> (CInt, CInt) -> IO Double + getV m (i', j') = do + (CDouble v) <- c_hkl_matrix_get m i' j' + return v + +foreign import ccall unsafe "hkl.h hkl_geometry_detector_rotation_get_binding" + c_hkl_geometry_detector_rotation_get_binding :: Ptr Geometry + -> Ptr HklDetector + -> IO (Ptr HklQuaternion) + +foreign import ccall unsafe "hkl.h hkl_quaternion_to_matrix_binding" + c_hkl_quaternion_to_matrix_binding :: Ptr HklQuaternion + -> IO (Ptr HklMatrix) + +foreign import ccall unsafe "hkl.h &hkl_quaternion_free" + c_hkl_quaternion_free :: FunPtr (Ptr HklQuaternion -> IO ()) + +foreign import ccall unsafe "hkl.h &hkl_matrix_free" + c_hkl_matrix_free :: FunPtr (Ptr HklMatrix -> IO ()) + +foreign import ccall unsafe "hkl.h hkl_matrix_get" + c_hkl_matrix_get :: Ptr HklMatrix + -> CInt + -> CInt + -> IO CDouble + + +peekItems :: Ptr HklGeometryList -> IO [Ptr HklGeometryListItem] +peekItems l = c_hkl_geometry_list_items_first_get l >>= unfoldrM go + where + go e + | e == nullPtr = return Nothing + | otherwise = do + next <- c_hkl_geometry_list_items_next_get l e + return (Just (e, next)) + +peekHklGeometryList :: ForeignPtr HklGeometryList -> IO [Geometry] +peekHklGeometryList l = withForeignPtr l $ \ls -> do + items <- peekItems ls + mapM extract items + where + extract it = c_hkl_geometry_list_item_geometry_get it >>= peek + +foreign import ccall unsafe "hkl.h hkl_geometry_list_items_first_get" + c_hkl_geometry_list_items_first_get :: Ptr HklGeometryList + -> IO (Ptr HklGeometryListItem) + +foreign import ccall unsafe "hkl.h hkl_geometry_list_items_next_get" + c_hkl_geometry_list_items_next_get :: Ptr HklGeometryList + -> Ptr HklGeometryListItem + -> IO (Ptr HklGeometryListItem) + +foreign import ccall unsafe "hkl.h hkl_geometry_list_item_geometry_get" + c_hkl_geometry_list_item_geometry_get :: Ptr HklGeometryListItem + -> IO (Ptr Geometry) diff --git a/contrib/haskell/src/Hkl/C/Lattice.hsc b/contrib/haskell/src/Hkl/C/Lattice.hsc new file mode 100644 index 0000000..5cb1d30 --- /dev/null +++ b/contrib/haskell/src/Hkl/C/Lattice.hsc @@ -0,0 +1,106 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE GADTs #-} + +module Hkl.C.Lattice + ( HklLattice + , newLattice + , withLattice + ) where + +import Prelude hiding (min, max) + +import Foreign ( ForeignPtr + , FunPtr + , Ptr + , nullPtr + , newForeignPtr + , withForeignPtr) +import Foreign.C (CDouble(..)) + +import Numeric.Units.Dimensional.Prelude ( meter, degree, radian, nano + , (*~), (/~)) +import Hkl.Lattice + +#include "hkl.h" + +#if __GLASGOW_HASKELL__ <= 710 +#let alignment t = "%lu", (unsigned long)offsetof(struct {char x__; t (y__); }, y__) +#endif + +-- private types + +data HklLattice + +-- Lattice + +withLattice :: Lattice a -> (Ptr HklLattice -> IO r) -> IO r +withLattice l func = do + fptr <- newLattice l + withForeignPtr fptr func + +newLattice' :: CDouble + -> CDouble + -> CDouble + -> CDouble + -> CDouble + -> CDouble + -> IO (ForeignPtr HklLattice) +newLattice' a b c alpha beta gamma = do + lattice <- c_hkl_lattice_new a b c alpha beta gamma nullPtr + newForeignPtr c_hkl_lattice_free lattice + +newLattice :: Lattice a -> IO (ForeignPtr HklLattice) +newLattice (Cubic la) = do + let a = CDouble (la /~ nano meter) + let alpha = CDouble ((90 *~ degree) /~ radian) + newLattice' a a a alpha alpha alpha +newLattice (Tetragonal la lc) = do + let a = CDouble (la /~ nano meter) + let c = CDouble (lc /~ nano meter) + let alpha = CDouble ((90 *~ degree) /~ radian) + newLattice' a a c alpha alpha alpha +newLattice (Orthorhombic la lb lc) = do + let a = CDouble (la /~ nano meter) + let b = CDouble (lb /~ nano meter) + let c = CDouble (lc /~ nano meter) + let alpha = CDouble ((90 *~ degree) /~ radian) + newLattice' a b c alpha alpha alpha +newLattice (Rhombohedral la aalpha) = do + let a = CDouble (la /~ nano meter) + let alpha = CDouble (aalpha /~ radian) + newLattice' a a a alpha alpha alpha +newLattice (Hexagonal la lc) = do + let a = CDouble (la /~ nano meter) + let c = CDouble (lc /~ nano meter) + let alpha = CDouble ((90 *~ degree) /~ radian) + let gamma = CDouble ((120 *~ degree) /~ radian) + newLattice' a a c alpha alpha gamma +newLattice (Monoclinic la lb lc abeta) = do + let a = CDouble (la /~ nano meter) + let b = CDouble (lb /~ nano meter) + let c = CDouble (lc /~ nano meter) + let alpha = CDouble ((90 *~ degree) /~ radian) + let beta = CDouble (abeta /~ radian) + newLattice' a b c alpha beta alpha +newLattice (Triclinic la lb lc aalpha abeta agamma) = do + let a = CDouble (la /~ nano meter) + let b = CDouble (lb /~ nano meter) + let c = CDouble (lc /~ nano meter) + let alpha = CDouble (aalpha /~ radian) + let beta = CDouble (abeta /~ radian) + let gamma = CDouble (agamma /~ radian) + newLattice' a b c alpha beta gamma + +foreign import ccall unsafe "hkl.h hkl_lattice_new" + c_hkl_lattice_new :: CDouble -- a + -> CDouble -- b + -> CDouble -- c + -> CDouble -- alpha + -> CDouble -- beta + -> CDouble -- gamma + -> Ptr () -- *gerror + -> IO (Ptr HklLattice) + +foreign import ccall unsafe "hkl.h &hkl_lattice_free" + c_hkl_lattice_free :: FunPtr (Ptr HklLattice -> IO ()) diff --git a/contrib/haskell/src/Hkl/C/Sample.hsc b/contrib/haskell/src/Hkl/C/Sample.hsc new file mode 100644 index 0000000..d9c106c --- /dev/null +++ b/contrib/haskell/src/Hkl/C/Sample.hsc @@ -0,0 +1,91 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE CPP #-} + +module Hkl.C.Sample + ( HklSample + , newSample + , withSample + ) where + +import Control.Monad (void) +import Foreign ( ForeignPtr + , FunPtr + , Ptr + , nullPtr + , newForeignPtr + , withForeignPtr) +import Foreign.C (CInt(..), CString, withCString) +import Foreign.Storable + +import Hkl.C.Lattice +import Hkl.Types + +#include "hkl.h" + +-- private types + +data HklSample + +-- Sample + +withSample :: Sample a -> (Ptr HklSample -> IO r) -> IO r +withSample s fun = do + fptr <- newSample s + withForeignPtr fptr fun + +newSample :: Sample a -> IO (ForeignPtr HklSample) +newSample (Sample name l ux uy uz) = + withCString name $ \cname -> do + sample <- c_hkl_sample_new cname + withLattice l $ \lattice -> do + c_hkl_sample_lattice_set sample lattice + go sample ux c_hkl_sample_ux_get c_hkl_sample_ux_set + go sample uy c_hkl_sample_uy_get c_hkl_sample_uy_set + go sample uz c_hkl_sample_uz_get c_hkl_sample_uz_set + newForeignPtr c_hkl_sample_free sample + where + go s p getter setter = do + fptr <- copyParameter =<< (getter s) + withForeignPtr fptr $ \ptr -> do + poke ptr p + void $ setter s ptr nullPtr + +foreign import ccall unsafe "hkl.h hkl_sample_new" + c_hkl_sample_new:: CString -> IO (Ptr HklSample) + +foreign import ccall unsafe "hkl.h hkl_sample_lattice_set" + c_hkl_sample_lattice_set :: Ptr HklSample -> Ptr HklLattice -> IO () + +foreign import ccall unsafe "hkl.h &hkl_sample_free" + c_hkl_sample_free :: FunPtr (Ptr HklSample -> IO ()) + +foreign import ccall unsafe "hkl.h hkl_sample_ux_get" + c_hkl_sample_ux_get :: Ptr HklSample + -> IO (Ptr Parameter) + +foreign import ccall unsafe "hkl.h hkl_sample_uy_get" + c_hkl_sample_uy_get :: Ptr HklSample + -> IO (Ptr Parameter) + +foreign import ccall unsafe "hkl.h hkl_sample_uz_get" + c_hkl_sample_uz_get :: Ptr HklSample + -> IO (Ptr Parameter) + +foreign import ccall unsafe "hkl.h hkl_sample_ux_set" + c_hkl_sample_ux_set :: Ptr HklSample + -> Ptr Parameter + -> Ptr () + -> IO CInt + +foreign import ccall unsafe "hkl.h hkl_sample_uy_set" + c_hkl_sample_uy_set :: Ptr HklSample + -> Ptr Parameter + -> Ptr () + -> IO CInt + +foreign import ccall unsafe "hkl.h hkl_sample_uz_set" + c_hkl_sample_uz_set :: Ptr HklSample + -> Ptr Parameter + -> Ptr () + -> IO CInt diff --git a/contrib/haskell/src/Hkl/DataSource.hs b/contrib/haskell/src/Hkl/DataSource.hs new file mode 100644 index 0000000..87a4b16 --- /dev/null +++ b/contrib/haskell/src/Hkl/DataSource.hs @@ -0,0 +1,51 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE StandaloneDeriving #-} + +module Hkl.DataSource ( ExtendDims(..) + , DataItem(..) + , DataSource(..) + , atIndex' + , openDataSource + , closeDataSource + ) where + +#if __GLASGOW_HASKELL__ < 710 +import Control.Applicative ((<$>), (<*>)) +#endif + +import Control.Monad.Trans.Maybe (MaybeT) +import Data.Array.Repa (Shape) +import Data.ByteString.Char8 (pack) +import Data.Vector.Storable (Vector, any, singleton) +import Pipes (lift) +import Prelude hiding ( any ) + +import Hkl.H5 + +data ExtendDims = ExtendDims | StrictDims deriving (Show) + +data DataItem a where + DataItemH5 :: H5Path -> ExtendDims -> DataItem H5 + DataItemConst :: Double -> DataItem Double +deriving instance Show (DataItem a) + +data DataSource a where + DataSourceH5 :: DataItem H5 -> Dataset -> DataSource H5 + DataSourceConst :: Double -> DataSource Double + +openDataSource :: File -> DataItem a -> IO (DataSource a) +openDataSource hid di@(DataItemH5 name _) = DataSourceH5 + <$> return di + <*> openDataset hid (pack name) Nothing +openDataSource _ (DataItemConst v) = return $ DataSourceConst v + +closeDataSource :: DataSource a -> IO () +closeDataSource (DataSourceH5 _ d) = closeDataset d +closeDataSource (DataSourceConst _) = return () + +atIndex' :: Shape sh => DataSource a -> sh -> MaybeT IO (Vector Double) +atIndex' (DataSourceH5 _ a ) b = lift $ do + v <- get_position_new a b + if any isNaN v then fail "File contains Nan" else return v +atIndex' (DataSourceConst v) _ = lift $ return $ singleton v diff --git a/contrib/haskell/src/Hkl/Detector.hs b/contrib/haskell/src/Hkl/Detector.hs new file mode 100644 index 0000000..f5ffaf4 --- /dev/null +++ b/contrib/haskell/src/Hkl/Detector.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE StandaloneDeriving #-} + +module Hkl.Detector + ( Detector(..) + , ImXpadS140 + , Xpad32 + , ZeroD + , coordinates + ) where + +import Data.Vector.Storable ( Vector + , fromList + ) + +import Hkl.PyFAI.Npt ( NptPoint ( NptPoint ) ) + + +data ImXpadS140 +data Xpad32 +data ZeroD + +data Detector a where + ImXpadS140 :: Detector ImXpadS140 + Xpad32 :: Detector Xpad32 + ZeroD :: Detector ZeroD + +deriving instance Show (Detector a) + +-- | Xpad Family + +type Gap = Double +type Width = Int +type Index = Int + +-- an xpad line is like this (pixel size, index) +-- | s 0 | s 1 | s 2 | ... | 5/2 s (w - 1) || 5/2 s w | s (w + 1) | ... +xpadLine :: Width -> Index -> Double +xpadLine w i' + | i' == 0 = s / 2 + | i' == 1 = s * 3 / 2 + | idx == 0 = s * (fromIntegral i' + 3 * fromIntegral c - 1 / 4) + | idx <= (w - 2) = s * (fromIntegral i' + 3 * fromIntegral c + 1 / 2) + | idx == (w - 1) = s * (fromIntegral i' + 3 * fromIntegral c + 5 / 4) + | otherwise = error $ "wront coordinates" ++ show i' + where + s = 130e-6 + (c, idx) = divMod i' w + +xpadLineWithGap :: Width -> Gap -> Index -> Double +xpadLineWithGap w g i' = s / 2 + (s * fromIntegral i') + g * fromIntegral (div i' w) + where + s = 130e-6 + +interp :: (Int -> Double) -> Double -> Double +interp f p + | p0 == p1 = f p0 + | otherwise = (p - fromIntegral p0) * (f p1 - f p0) + f p0 + where + p0 :: Int + p0 = floor p + + p1 :: Int + p1 = ceiling p + +-- compute the coordinated at a given point + +coordinates :: Detector a -> NptPoint -> Vector Double +coordinates ZeroD (NptPoint 0 0) = fromList [0, 0, 0] +coordinates ZeroD _ = error "No coordinates in a ZeroD detecteor" + +coordinates ImXpadS140 (NptPoint x y) = + fromList [ interp (xpadLine 120) y + , interp (xpadLine 80) x + , 0 + ] + +coordinates Xpad32 (NptPoint x y) = + fromList [ interp (xpadLineWithGap 120 3.57e-3) y + , interp (xpadLine 80) x + , 0] diff --git a/contrib/haskell/src/Hkl/Edf.hs b/contrib/haskell/src/Hkl/Edf.hs new file mode 100644 index 0000000..4c33739 --- /dev/null +++ b/contrib/haskell/src/Hkl/Edf.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Edf + ( Edf(..) + , ExtractEdf(..) + , edfP + , edfFromFile + ) where + +import Data.Attoparsec.Text ( Parser + , (<?>) + , anyChar + , double + , many1 + , manyTill + , parseOnly + , skipSpace + , string + , takeTill + , try + ) +import Data.ByteString.Char8 (readFile, split) +import Data.Text (Text, words) +import Data.Text.Encoding (decodeUtf8) +import Numeric.Units.Dimensional.Prelude (Length, (*~), nano, meter) + +data Edf = Edf { edf'Lambda :: Length Double + , edf'Motors :: [(Text, Double)] + } + deriving (Show) + +class ExtractEdf a where + extractEdf ∷ a → IO () + + +edf'LambdaP :: Parser (Length Double) +edf'LambdaP = do + _ <- manyTill anyChar (try $ string "Lambda = ") + value <- double + pure $ value *~ nano meter + +edf'MotorsP :: Parser [(Text, Double)] +edf'MotorsP = do + _ <- manyTill anyChar (try $ string "motor_pos = ") + vs <- many1 (skipSpace *> double) + _ <- manyTill anyChar (try $ string "motor_mne = ") + ns <- takeTill (\c -> c == ';') + return $ zip (Data.Text.words ns) vs + +edfP :: Parser Edf +edfP = Edf + <$> edf'LambdaP + <*> edf'MotorsP + <?> "edfP" + +edfFromFile :: FilePath -> IO Edf +edfFromFile filename = do + content <- Data.ByteString.Char8.readFile filename + let header = head (split '}' content) + return $ case parseOnly edfP (decodeUtf8 header) of + Left _ -> error $ "Can not parse the " ++ filename ++ " edf file" + Right a -> a + +-- main :: IO () +-- main = do +-- edf <- edfFromFile "/home/picca/test.edf" +-- print edf +-- return () diff --git a/contrib/haskell/src/Hkl/Engine.hs b/contrib/haskell/src/Hkl/Engine.hs new file mode 100644 index 0000000..56cb3c9 --- /dev/null +++ b/contrib/haskell/src/Hkl/Engine.hs @@ -0,0 +1,27 @@ +module Hkl.Engine + ( enginesTrajectoryPipe + , fromToPipe + ) where + +import Control.Monad ( forever, forM_ ) +import Numeric.LinearAlgebra ( Vector, toList ) +import Pipes ( Pipe, Producer, await, yield ) + +import Hkl.Types ( Engine ( Engine ) + , Parameter ( Parameter ) + ) + +engineSetValues :: Engine -> Vector Double -> Engine +engineSetValues (Engine name ps mode) vs = Engine name nps mode + where + nps = zipWith set ps (toList vs) + set (Parameter n _ range) newValue = Parameter n newValue range + +fromToPipe :: Int -> Vector Double -> Vector Double -> Producer (Vector Double) IO () +fromToPipe n from to = forM_ [0..n-1] $ \i -> yield $ vs i + where + vs i = from + step * fromIntegral i + step = (to - from) / (fromIntegral n - 1) + +enginesTrajectoryPipe :: Engine -> Pipe (Vector Double) Engine IO () +enginesTrajectoryPipe e = forever $ await >>= yield . engineSetValues e diff --git a/contrib/haskell/src/Hkl/Flat.hs b/contrib/haskell/src/Hkl/Flat.hs new file mode 100644 index 0000000..62746e4 --- /dev/null +++ b/contrib/haskell/src/Hkl/Flat.hs @@ -0,0 +1,81 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Flat + ( Flat(..) + , Npy + , computeFlat + ) + where + +import Data.Text ( unlines, pack ) +import System.Exit ( ExitCode( ExitSuccess ) ) +import System.FilePath.Posix ( replaceExtension ) + +import Hkl.DataSource ( DataItem ( DataItemH5 ) ) +import Hkl.Nxs ( Nxs ( Nxs ) + , XrdFlat + , DataFrameH5Path ( XrdFlatH5Path ) + ) +import Hkl.Python ( PyVal + , toPyVal + ) +import Hkl.Script ( Py2 + , Script ( Py2Script ) + , run + ) + +data Npy + +data Flat a where + FlatNpy ∷ FilePath → Flat Npy +deriving instance (Show) (Flat a) + +scriptPy2Flat ∷ [Nxs XrdFlat] → FilePath → Script Py2 +scriptPy2Flat ns output = Py2Script (script, scriptName) + where + script = Data.Text.unlines $ + map pack ["#!/bin/env python" + , "" + , "import numpy" + , "from h5py import File" + , "" + , "NEXUSFILES = " ++ toPyVal nxs' + , "IMAGEPATHS = " ++ toPyVal hpaths + , "OUTPUT = " ++ toPyVal output + , "" + , "flat = None" + , "n = None" + , "with File(NEXUSFILES[0], mode='r') as f:" + , " imgs = f[IMAGEPATHS[0]]" + , " flat = numpy.sum(imgs[:], axis=0)" + , " n = imgs.shape[0]" + , "for idx, (nxs, h5path) in enumerate(zip(NEXUSFILES[1:], IMAGEPATHS[1:])):" + , " with File(nxs, mode='r') as f:" + , " imgs = f[h5path]" + , " flat += numpy.sum(imgs[:], axis=0)" + , " n += imgs.shape[0]" + , "numpy.save(OUTPUT, flat.astype('f') / n)" + ] + nxs' ∷ [String] + nxs' = [f | (Nxs f _) ← ns] + + hpaths ∷ [String] + hpaths = [h | (Nxs _ (XrdFlatH5Path (DataItemH5 h _))) ← ns] + + scriptName ∷ FilePath + scriptName = output `replaceExtension` "py" + +computeFlat ∷ [Nxs XrdFlat] → FilePath → IO (Flat Npy) +computeFlat ns o = do + -- create the python script. + let script = scriptPy2Flat ns o + -- execute this script. + ExitSuccess ← run script False + -- return the filepath of the generated file. + return (FlatNpy o) + +instance PyVal (Flat a) where + toPyVal (FlatNpy v) = "numpy.load(" ++ show v ++ ")" diff --git a/contrib/haskell/src/Hkl/H5.hs b/contrib/haskell/src/Hkl/H5.hs new file mode 100644 index 0000000..5858bfa --- /dev/null +++ b/contrib/haskell/src/Hkl/H5.hs @@ -0,0 +1,194 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE UnicodeSyntax #-} +module Hkl.H5 + ( Dataset + , File + , H5 + , H5Path + , check_ndims + , closeDataset + , closeFile + , get_position + , get_position_new + , get_ub + , lenH5Dataspace + , nxEntries + , openDataset + , openH5 + , withH5File + ) + where + +import Bindings.HDF5.Core ( HSize(HSize) + , IndexType(ByName) + , IterOrder(Native) + , hid + , hSize + , indexTypeCode + , iterOrderCode + ) +import Bindings.HDF5.File ( File + , AccFlags(ReadOnly) + , openFile + , closeFile + ) +import Bindings.HDF5.Dataset ( Dataset + , openDataset + , closeDataset + , getDatasetSpace + , readDataset + , readDatasetInto + ) +import Bindings.HDF5.Dataspace ( Dataspace + , SelectionOperator(Set) + , closeDataspace + , createSimpleDataspace + , getSimpleDataspaceExtentNDims + , getSimpleDataspaceExtentNPoints + , selectHyperslab + ) +import Bindings.HDF5.Raw ( HErr_t(HErr_t) + , HId_t(HId_t) + , H5L_info_t + , h5l_iterate + ) +import Control.Exception (bracket) +import Data.Array.Repa (Shape, listOfShape) +import Data.ByteString.Char8 ( pack ) +import Data.IORef ( newIORef, readIORef, writeIORef ) +import Data.Vector.Storable (Vector, freeze) +import Data.Vector.Storable.Mutable (replicate) +import Foreign.StablePtr ( castPtrToStablePtr + , castStablePtrToPtr + , deRefStablePtr + , freeStablePtr + , newStablePtr + ) +import Foreign.Ptr ( FunPtr, freeHaskellFunPtr ) +import Foreign.Ptr.Conventions ( In(In) + , InOut(InOut) + , castWrappedPtr + , withInOut_ + ) +import Foreign.C.String ( CString, peekCString ) +import Foreign.C.Types (CInt(CInt)) +import Numeric.LinearAlgebra (Matrix, reshape) + +{-# ANN module "HLint: ignore Use camelCase" #-} + +data H5 + +type H5Path = String + + +check_ndims :: Dataset -> Int -> IO Bool +check_ndims d expected = do + space_id <- getDatasetSpace d + (CInt ndims) <- getSimpleDataspaceExtentNDims space_id + return $ expected == fromEnum ndims + +toHyperslab :: Shape sh => sh -> [(HSize, Maybe HSize, HSize, Maybe HSize)] +toHyperslab s = [(HSize (fromIntegral n), Just (HSize 1), HSize 1, Just (HSize 1)) | n <- listOfShape s] + +get_position_new :: Shape sh => Dataset -> sh -> IO (Vector Double) +get_position_new dataset s = + withDataspace dataset $ \dataspace -> do + selectHyperslab dataspace Set (toHyperslab s) + withDataspace' $ \memspace -> do + data_out <- Data.Vector.Storable.Mutable.replicate 1 (0.0 :: Double) + readDatasetInto dataset (Just memspace) (Just dataspace) Nothing data_out + freeze data_out + +get_position :: Dataset -> Int -> IO (Vector Double) +get_position dataset n = + withDataspace dataset $ \dataspace -> do + let start = HSize (fromIntegral n) + let stride = Just (HSize 1) + let count = HSize 1 + let block = Just (HSize 1) + selectHyperslab dataspace Set [(start, stride, count, block)] + withDataspace' $ \memspace -> do + data_out <- Data.Vector.Storable.Mutable.replicate 1 (0.0 :: Double) + readDatasetInto dataset (Just memspace) (Just dataspace) Nothing data_out + freeze data_out + +get_ub :: Dataset -> IO (Matrix Double) +get_ub dataset = do + v <- readDataset dataset Nothing Nothing + return $ reshape 3 v + +-- | File + +withH5File :: FilePath -> (File -> IO r) -> IO r +withH5File fp = bracket acquire release + where + acquire = openFile (pack fp) [ReadOnly] Nothing + release = closeFile + +openH5 ∷ FilePath → IO File +openH5 f = openFile (pack f) [ReadOnly] Nothing + +-- | Dataspace + +-- check how to merge both methods + +withDataspace' :: (Dataspace -> IO r) -> IO r +withDataspace' = bracket acquire release + where + acquire = createSimpleDataspace [HSize 1] + release = closeDataspace + +withDataspace :: Dataset -> (Dataspace -> IO r) -> IO r +withDataspace d = bracket acquire release + where + acquire = getDatasetSpace d + release = closeDataspace + +lenH5Dataspace :: Dataset -> IO (Maybe Int) +lenH5Dataspace = withDataspace'' len + where + withDataspace'' f d = withDataspace d f + len space_id = do + (HSize n) <- getSimpleDataspaceExtentNPoints space_id + return $ if n < 0 then Nothing else Just (fromIntegral n) + + +-- | WIP until I have decided what is the right way to go + +type H5Iterate a = HId_t -> CString -> In H5L_info_t -> InOut a -> IO HErr_t + +foreign import ccall "wrapper" mkOp :: H5Iterate a -> IO (FunPtr (H5Iterate a)) + +nxEntries ∷ FilePath → IO [String] +nxEntries f = withH5File f $ \h → do + state <- newIORef [] + statePtr <- newStablePtr state + let opData = InOut $ castStablePtrToPtr statePtr + let startIndex = Nothing + let indexType = ByName + let order = Native + iop <- mkOp callback + _ <- withInOut_ (maybe 0 hSize startIndex) $ \ioStartIndex -> + h5l_iterate (hid h) (indexTypeCode indexType) (iterOrderCode order) ioStartIndex iop opData + + freeHaskellFunPtr iop + freeStablePtr statePtr + + -- retrieve the final state + readIORef state + where + callback ∷ H5Iterate a + callback _g n _i (InOut dataptr) = + do + let opData = castWrappedPtr dataptr + -- get the state + stRef <- deRefStablePtr (castPtrToStablePtr opData) + st <- readIORef stRef + + -- compute the new state + name <- peekCString n + let newSt = st ++ [name] + + -- store the new state + writeIORef stRef newSt + return $ HErr_t 0 diff --git a/contrib/haskell/src/Hkl/Lattice.hs b/contrib/haskell/src/Hkl/Lattice.hs new file mode 100644 index 0000000..9578402 --- /dev/null +++ b/contrib/haskell/src/Hkl/Lattice.hs @@ -0,0 +1,63 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Lattice ( Lattice(..) + , Cubic + , Tetragonal + , Orthorhombic + , Rhombohedral + , Hexagonal + , Monoclinic + , Triclinic + ) where + +import Numeric.Units.Dimensional.Prelude (Length, Angle) + +-- | Lattice + +data Cubic +data Tetragonal +data Orthorhombic +data Rhombohedral +data Hexagonal +data Monoclinic +data Triclinic + +data Lattice a where + -- ^ a = b = c, alpha = beta = gamma = 90 + Cubic ∷ Length Double + → Lattice Cubic -- ^ a = b = c, alpha = beta = gamma = 90 + -- a = b != c, alpha = beta = gamma = 90 + Tetragonal ∷ (Length Double) -- ^ a, b + → (Length Double) -- ^ c + → Lattice Tetragonal + -- ^ a != b != c, alpha = beta = gamma = 90 + Orthorhombic ∷ (Length Double) -- ^ a + → (Length Double) -- ^ b + → (Length Double) -- ^ c + → Lattice Orthorhombic + -- ^ a = b = c, alpha = beta = gamma != 90 + Rhombohedral ∷ (Length Double) -- ^ a, b, c + → (Angle Double) -- ^ alpha, beta, gamma + → Lattice Rhombohedral + -- ^ a = b != c, alpha = beta = 90, gamma = 120 + Hexagonal ∷ (Length Double) -- ^ a, b + → (Length Double) -- ^ c + → Lattice Hexagonal + -- a != b != c, alpha = gamma = 90, beta != 90 + Monoclinic ∷ (Length Double) -- ^ a + → (Length Double) -- ^ b + → (Length Double) -- ^ c + → (Angle Double) -- ^ beta + → Lattice Monoclinic + -- a != b != c, alpha != beta != gamma != 90 + Triclinic ∷ (Length Double) + → (Length Double) -- ^ b + → (Length Double) -- ^ c + → (Angle Double) -- ^ alpha + → (Angle Double) -- ^ beta + → (Angle Double) -- ^ gamma + → Lattice Triclinic + +deriving instance Show (Lattice a) diff --git a/contrib/haskell/src/Hkl/MyMatrix.hs b/contrib/haskell/src/Hkl/MyMatrix.hs new file mode 100644 index 0000000..57877d9 --- /dev/null +++ b/contrib/haskell/src/Hkl/MyMatrix.hs @@ -0,0 +1,50 @@ +module Hkl.MyMatrix + ( Basis(..) + , MyMatrix(..) + , changeBase + , toEulerians + ) where + +import Numeric.LinearAlgebra (Matrix, atIndex, fromLists, inv, (<>)) +import Numeric.Units.Dimensional.Prelude (Angle, (*~), radian) + +data Basis = PyFAIB -- the pyFAI (1, 2, 3) detector coordinates + | HklB -- the hkl coordinates + deriving (Show) + +data MyMatrix a = MyMatrix Basis (Matrix a) deriving (Show) + +changeBase :: MyMatrix Double -> Basis -> MyMatrix Double +changeBase (MyMatrix PyFAIB m) HklB = MyMatrix HklB (passage m p2) +changeBase (MyMatrix HklB m) PyFAIB = MyMatrix PyFAIB (passage m p1) +changeBase m@(MyMatrix PyFAIB _) PyFAIB = m +changeBase m@(MyMatrix HklB _) HklB = m + +passage :: Matrix Double -> Matrix Double -> Matrix Double +passage r p = inv p <> r <> p + +p1 :: Matrix Double -- hkl -> pyFAI +p1 = fromLists [ [0, 0, 1] + , [0, -1, 0] + , [1, 0, 0]] + +p2 :: Matrix Double -- pyFAI -> hkl: +p2 = fromLists [ [ 0, 0, 1] + , [ 0, -1, 0] + , [ 1, 0, 0]] + +toEulerians :: Matrix Double -> (Angle Double, Angle Double, Angle Double) +toEulerians m + | abs c > epsilon = ( atan2 ((m `atIndex` (2, 1)) / c) ((m `atIndex` (2, 2)) / c) *~ radian + , rot2 *~ radian + , atan2 ((m `atIndex` (1, 0)) / c) ((m `atIndex` (0, 0)) / c) *~ radian + ) + | otherwise = ( 0 *~ radian + , rot2 *~ radian + , atan2 (-(m `atIndex` (0, 1))) (m `atIndex` (1, 1)) *~ radian + ) + where + epsilon = 1e-10 + rot2 = asin (-(m `atIndex` (2, 0))) + c = cos rot2 + diff --git a/contrib/haskell/src/Hkl/Nxs.hs b/contrib/haskell/src/Hkl/Nxs.hs new file mode 100644 index 0000000..a7934cc --- /dev/null +++ b/contrib/haskell/src/Hkl/Nxs.hs @@ -0,0 +1,237 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Nxs + ( DataFrameH5(..) + , DataFrameH5Path(..) + , NxEntry + , Nxs(..) + , PoniGenerator + , XrdFlat + , XrdOneD + , XrdMesh + , XrdZeroD + , mkNxs + , withDataFrameH5 + , withDataSource + ) where + +import Bindings.HDF5.Dataset ( readDataset + , getDatasetSpace ) +import Bindings.HDF5.Dataspace ( getSimpleDataspaceExtent ) +import Codec.Picture ( DynamicImage( ImageY16 ) + , Image ( Image ) + ) +#if __GLASGOW_HASKELL__ < 710 +import Control.Applicative ((<$>), (<*>), pure) +#endif +import Control.Exception.Base (bracket) +import Control.Monad.IO.Class (liftIO) +import Pipes.Safe ( MonadSafe, bracket ) + +import Hkl.DataSource ( DataItem + , DataSource ( DataSourceH5 ) + , closeDataSource + , openDataSource + ) +import Hkl.H5 ( File, H5 + , closeFile + , openH5 + ) +import Hkl.PyFAI ( Pose, PoniExt ) +import Hkl.Tiff ( ToTiff + , toTiff + ) + +type NxEntry = String + +-- to remove an put directly into OneD +type PoniGenerator = Pose -> Int -> IO PoniExt + +data XrdFlat +data XrdOneD +data XrdMesh +data XrdZeroD + +data DataFrameH5Path a where + XrdFlatH5Path ∷ (DataItem H5) -- ^ image + → DataFrameH5Path XrdFlat + XrdOneDH5Path ∷ (DataItem H5) -- ^ image + → (DataItem H5) -- ^ gamma + → (DataItem H5) -- ^ delta + → (DataItem H5) -- ^ wavelength + → DataFrameH5Path XrdOneD + XrdMeshH5Path ∷ (DataItem H5) -- ^ Image + → (DataItem H5) -- ^ meshx + → (DataItem H5) -- ^ meshy + → (DataItem H5) -- ^ gamma + → (DataItem H5) -- ^ delta + → (DataItem H5) -- ^ wavelength + → DataFrameH5Path XrdMesh + XrdMeshFlyH5Path ∷ (DataItem H5) -- ^ Image + → (DataItem H5) -- ^ meshx + → (DataItem H5) -- ^ meshy + → (DataItem Double) -- ^ gamma + → (DataItem Double) -- ^ delta + → (DataItem Double) -- ^ wavelength + → DataFrameH5Path XrdMesh + XrdZeroDH5Path ∷ (DataItem H5) -- ^ image + → (DataItem Double) -- ^ wavelength + → DataFrameH5Path XrdZeroD -- used to integrate one static image + +deriving instance Show (DataFrameH5Path a) + +data Nxs a where + Nxs ∷ FilePath → DataFrameH5Path a → Nxs a + +deriving instance Show (Nxs a) + +data DataFrameH5 a where + XrdFlatH5 ∷ (Nxs XrdFlat) -- Nexus Source file + → File -- h5file handler + → (DataSource H5) --images + → DataFrameH5 XrdFlat + DataFrameH5 ∷ (Nxs XrdOneD) -- Nexus file + → File -- h5file handler + → (DataSource H5) -- gamma + → (DataSource H5) -- delta + → (DataSource H5) -- wavelength + → PoniGenerator -- ponie generator + → DataFrameH5 XrdOneD + XrdMeshH5 ∷ (Nxs XrdMesh) -- NexusFile Source File + → File -- h5file handler + → (DataSource H5) -- image + → (DataSource H5) -- meshx + → (DataSource H5) -- meshy + → (DataSource H5) -- gamma + → (DataSource H5) -- delta + → (DataSource H5) -- wavelength + → DataFrameH5 XrdMesh + XrdMeshFlyH5 ∷ (Nxs XrdMesh) -- NexusFile Source File + → File -- h5file handler + → (DataSource H5) -- image + → (DataSource H5) -- meshx + → (DataSource H5) -- meshy + → (DataSource Double) -- gamma + → (DataSource Double) -- delta + → (DataSource Double) -- wavelength + → DataFrameH5 XrdMesh + XrdZeroDH5 ∷ (Nxs XrdZeroD) -- NexusFile Source File + → File -- h5file handler + → (DataSource H5) -- image + → (DataSource Double) -- wavelength + → DataFrameH5 XrdZeroD + +mkNxs ∷ FilePath → NxEntry → (NxEntry → DataFrameH5Path a) → Nxs a +mkNxs f e h = Nxs f (h e) + +-- | Instanciate a DataFrameH5 from a DataFrameH5Path +-- acquire and release the resources + +after ∷ DataFrameH5 a → IO () +after (XrdFlatH5 _ f i) = do + closeDataSource i + closeFile f +after (DataFrameH5 _ f g d w _) = do + closeDataSource g + closeDataSource d + closeDataSource w + closeFile f +after (XrdMeshH5 _ f i x y g d w) = do + closeDataSource i + closeDataSource x + closeDataSource y + closeDataSource g + closeDataSource d + closeDataSource w + closeFile f +after (XrdMeshFlyH5 _ f i x y g d w) = do + closeDataSource i + closeDataSource x + closeDataSource y + closeDataSource g + closeDataSource d + closeDataSource w + closeFile f +after (XrdZeroDH5 _ f i w) = do + closeDataSource i + closeDataSource w + closeFile f + +before :: Nxs a → IO (DataFrameH5 a) +before nxs'@(Nxs f (XrdFlatH5Path i)) = do + h ← openH5 f + XrdFlatH5 + <$> return nxs' + <*> return h + <*> openDataSource h i +-- before nxs'@(Nxs f (XrdOneDH5Path i g d w)) = do +-- h ← openH5 f +-- DataFrameH5 +-- <$> return nxs' +-- <*> return h +-- <*> openDataSource h g +-- <*> openDataSource h d +-- <*> openDataSource h w +-- <*> return gen +before nxs'@(Nxs f (XrdMeshH5Path i x y g d w)) = do + h ← openH5 f + XrdMeshH5 + <$> return nxs' + <*> return h + <*> openDataSource h i + <*> openDataSource h x + <*> openDataSource h y + <*> openDataSource h g + <*> openDataSource h d + <*> openDataSource h w +before nxs'@(Nxs f (XrdMeshFlyH5Path i x y g d w))= do + h ← openH5 f + XrdMeshFlyH5 + <$> return nxs' + <*> return h + <*> openDataSource h i + <*> openDataSource h x + <*> openDataSource h y + <*> openDataSource h g + <*> openDataSource h d + <*> openDataSource h w +before nxs'@(Nxs f (XrdZeroDH5Path i w)) = do + h ← openH5 f + XrdZeroDH5 + <$> return nxs' + <*> return h + <*> openDataSource h i + <*> openDataSource h w + +withDataSource :: Nxs a -> (DataFrameH5 a -> IO r) -> IO r +withDataSource s = Control.Exception.Base.bracket (before s) after + +-- | Pipe + +withDataFrameH5 :: (MonadSafe m) => Nxs XrdOneD -> PoniGenerator -> (DataFrameH5 XrdOneD -> m r) -> m r +withDataFrameH5 nxs'@(Nxs f (XrdOneDH5Path _ g d w)) gen = Pipes.Safe.bracket (liftIO before') (liftIO . after) + where + -- before :: File -> DataFrameH5Path -> m DataFrameH5 + before' :: IO (DataFrameH5 XrdOneD) + before' = do + h ← openH5 f + DataFrameH5 + <$> return nxs' + <*> return h + <*> openDataSource h g + <*> openDataSource h d + <*> openDataSource h w + <*> return gen + +instance ToTiff (Nxs XrdFlat) where + toTiff n = withDataSource n $ + \(XrdFlatH5 _ _ (DataSourceH5 _ i)) → do + ([w, h], _) ← getSimpleDataspaceExtent =<< (getDatasetSpace i) + ImageY16 <$> ( Image + <$> pure (fromIntegral w) + <*> pure (fromIntegral h) + <*> readDataset i Nothing Nothing ) diff --git a/contrib/haskell/src/Hkl/Projects.hs b/contrib/haskell/src/Hkl/Projects.hs new file mode 100644 index 0000000..0a69776 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects.hs @@ -0,0 +1,6 @@ +module Hkl.Projects ( module X ) where + +import Hkl.Projects.D2AM as X +import Hkl.Projects.Diffabs as X +import Hkl.Projects.Mars as X +import Hkl.Projects.Sixs as X diff --git a/contrib/haskell/src/Hkl/Projects/D2AM.hs b/contrib/haskell/src/Hkl/Projects/D2AM.hs new file mode 100644 index 0000000..1a71b06 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/D2AM.hs @@ -0,0 +1,3 @@ +module Hkl.Projects.D2AM (module X) where + +import Hkl.Projects.D2AM.XRD as X diff --git a/contrib/haskell/src/Hkl/Projects/D2AM/XRD.hs b/contrib/haskell/src/Hkl/Projects/D2AM/XRD.hs new file mode 100644 index 0000000..0b431af --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/D2AM/XRD.hs @@ -0,0 +1,105 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} + +module Hkl.Projects.D2AM.XRD + ( d2am ) where + +import Data.Array.Repa (DIM1, ix1) +-- import Data.Char (toUpper) +import Numeric.LinearAlgebra (ident) +import System.FilePath ((</>)) +import Text.Printf (printf) + +import Prelude hiding (concat, lookup, readFile, writeFile) + +import Hkl.MyMatrix +import Hkl.PyFAI +import Hkl.Xrd +import Hkl.Detector + +-- | Samples + +project :: FilePath +project = "/home/experiences/instrumentation/picca/data/d2am" +-- project = "/nfs/ruche-diffabs/diffabs-soleil/com-diffabs/" + +published :: FilePath +published = project </> "published-data" + +sampleRef :: XRDRef +sampleRef = XRDRef "reference" + (published </> "calibration") + (XrdRefEdf + (project </> "16Dec08D5_0268-rsz.edf") + (project </> "16Dec08D5_0268-rsz.poni") + ) + +sampleCalibration :: XRDCalibration Xpad32 +sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" + , xrdCalibrationOutputDir = published </> "calibration" -- TODO pourquoi ce output + , xrdCalibrationDetector = Xpad32 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + + idxs :: [Int] + idxs = [268, 271, 285, 295] + + entry :: Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryEdf + { xrdCalibrationEntryEdf'Edf = project </> printf "16Dec08D5_%04d-rsz.edf" idx + , xrdCalibrationEntryEdf'NptPath = project </> printf "16Dec08D5_%04d-rsz.npt" idx + } + + entries :: [XRDCalibrationEntry] + entries = [ entry idx | idx <- idxs] + +bins :: DIM1 +bins = ix1 1000 + +multibins :: DIM1 +multibins = ix1 10000 + +threshold :: Maybe Threshold +threshold = Just (Threshold 5000) + +skipedFrames :: [Int] +skipedFrames = [] + +lab6 :: XRDSample +lab6 = XRDSample "test" + (published </> "test") + [XrdNxs bins multibins threshold skipedFrames entries] + where + idxs :: [Int] + idxs = [268, 271, 285, 295] + + entry :: Int -> FilePath + entry idx = project </> printf "16Dec08D5_%04d-rsz.edf" idx + + entries :: XrdSource + entries = XrdSourceEdf [entry idx | idx <- idxs] + +-- | Main + +d2am :: IO () +d2am = do + let samples = [lab6] + + p <- getPoniExtRef sampleRef + + -- let poniextref = setPose (Hkl.PyFAI.PoniExt.flip p) (MyMatrix HklB (ident 3)) + let poniextref = move p (Pose (MyMatrix HklB (ident 3))) + + -- full calibration + poniextref' <- calibrate sampleCalibration poniextref + + print poniextref + print poniextref' + + -- integrate each step of the scan + let params = XrdOneDParams poniextref' Nothing Csr -- waiting for PyFAI to manage method in multi geometry + integrateMulti params samples + + return () diff --git a/contrib/haskell/src/Hkl/Projects/Diffabs.hs b/contrib/haskell/src/Hkl/Projects/Diffabs.hs new file mode 100644 index 0000000..1a7b753 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Diffabs.hs @@ -0,0 +1,9 @@ +module Hkl.Projects.Diffabs (module X) where + +import Hkl.Projects.Diffabs.Charlier as X +import Hkl.Projects.Diffabs.Hamon as X +import Hkl.Projects.Diffabs.Hercules as X +import Hkl.Projects.Diffabs.IRDRx as X +import Hkl.Projects.Diffabs.Laure as X +import Hkl.Projects.Diffabs.Martinetto as X +import Hkl.Projects.Diffabs.Melle as X diff --git a/contrib/haskell/src/Hkl/Projects/Diffabs/Charlier.hs b/contrib/haskell/src/Hkl/Projects/Diffabs/Charlier.hs new file mode 100644 index 0000000..49d28a6 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Diffabs/Charlier.hs @@ -0,0 +1,164 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} + +module Hkl.Projects.Diffabs.Charlier + ( charlier ) where + +import Data.Array.Repa (DIM1, ix1) +import Numeric.LinearAlgebra (ident) +import System.FilePath ((</>)) +import Text.Printf (printf) + +import Prelude hiding (concat, lookup, readFile, writeFile) + +import Hkl + +-- | TODO +-- ⋅ gerer le dummy correctement en focntion du type de données des images uint32, int16 +-- ∘ couper la fin du spectre qui nous embète. +-- | Samples + +project :: FilePath +project = "/nfs/ruche-diffabs/diffabs-users/20151386/" + +published :: FilePath +published = project </> "published-data" </> "xrd" + +-- | Calibration part + +project' :: FilePath +project' = "/nfs/ruche-diffabs/diffabs-users/99160066/" + +published':: FilePath +published' = project' </> "published-data" + +h5path' :: NxEntry -> DataFrameH5Path XrdOneD +h5path' nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_53" + gamma = "d13-1-cx1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/actuator_1_1" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +sampleCalibration :: XRDCalibration Xpad32 +sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" + , xrdCalibrationOutputDir = published' </> "calibration" + , xrdCalibrationDetector = Xpad32 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + idxs :: [Int] + idxs = [3, 6, 9, 15, 18, 21, 24, 27, 30, 33, 36, 39, 43] + + entry :: Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryNxs + { xrdCalibrationEntryNxs'Nxs = mkNxs (published' </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path' + , xrdCalibrationEntryNxs'Idx = idx + , xrdCalibrationEntryNxs'NptPath = published' </> "calibration" </> printf "XRD18keV_26.nxs_%02d.npt" idx + } + + entries :: [XRDCalibrationEntry] + entries = [ entry idx | idx <- idxs] + + +sampleRef :: XRDRef +sampleRef = XRDRef "reference" + (published' </> "calibration") + (XrdRefNxs + (mkNxs (published' </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path') + 6 -- BEWARE only the 6th poni was generated with the right Xpad_flat geometry. + ) + +bins :: DIM1 +bins = ix1 8000 + +multibins :: DIM1 +multibins = ix1 25000 + +threshold :: Maybe Threshold +threshold = Just (Threshold 1200) + +h5path :: NxEntry -> DataFrameH5Path XrdMesh +h5path nxentry = XrdMeshH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> meshx) StrictDims) + (DataItemH5 (nxentry </> meshy) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) StrictDims) + (DataItemH5 (nxentry </> beamline </> delta) StrictDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_53" + meshx = "scan_data/actuator_1_1" + meshy = "scan_data/actuator_2_1" + gamma = "d13-1-cx1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "d13-1-cx1__EX__DIF.1-DELTA__#1/raw_value" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +charlemagne :: XrdMeshSample +charlemagne = XrdMeshSample "Charlemagne" + (published </> "Charlemagne") + [ XrdMesh bins multibins threshold (XrdMeshSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-23" </> "XRD18keV_31.nxs") "scan_31" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-23" </> "XRD18keV_32.nxs") "scan_32" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-23" </> "XRD18keV_33.nxs") "scan_33" h5path + ] + ] + +charlesLeChauve :: XrdMeshSample +charlesLeChauve = XrdMeshSample "Charles le Chauve" + (published </> "Charles le Chauve") + [ XrdMesh bins multibins threshold (XrdMeshSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-24" </> "XRD18keV_34.nxs") "scan_34" h5path ] + ] + +louisLePieux :: XrdMeshSample +louisLePieux = XrdMeshSample "Louis le Pieux" + (published </> "Louis Le Pieux") + [ XrdMesh bins multibins threshold (XrdMeshSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-24" </> "XRD18keV_35.nxs") "scan_35" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-24" </> "XRD18keV_36.nxs") "scan_36" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-24" </> "XRD18keV_37.nxs") "scan_37" h5path + ] + ] + +-- | Main + +charlier :: IO () +charlier = do + let samples = [ charlemagne, charlesLeChauve, louisLePieux] + -- let samples = [ louisLePieux ] + -- # need to run f30 by itself because of a segfault in the hkl library + -- for now f30 whcih is an incomplet scan stop the script so put it at the end. + -- let samples = [f30, ceo2] + -- let samples = [ceo2] + let mflat = Nothing + let method = CsrOcl + + p <- getPoniExtRef sampleRef + + -- flip the ref poni in order to fit the reality + -- let poniextref = p + let poniextref = move p (Pose (MyMatrix HklB (ident 3))) + -- let poniextref = setPose (Hkl.PyFAI.PoniExt.flip p) (MyMatrix HklB (ident 3)) + + -- full calibration + poniextref' <- calibrate sampleCalibration poniextref + -- print p + print poniextref + print poniextref' + + -- integrate each step of the scan + integrateMesh (XrdMeshParams poniextref' mflat method) samples + return () diff --git a/contrib/haskell/src/Hkl/Projects/Diffabs/Hamon.hs b/contrib/haskell/src/Hkl/Projects/Diffabs/Hamon.hs new file mode 100644 index 0000000..8923efc --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Diffabs/Hamon.hs @@ -0,0 +1,134 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} + +module Hkl.Projects.Diffabs.Hamon + ( hamon ) where + +import Data.Array.Repa (DIM1, ix1) +import Numeric.LinearAlgebra (ident) +import System.FilePath ((</>)) +import Text.Printf (printf) + +import Prelude hiding (concat, lookup, readFile, writeFile) + +import Hkl + +-- | TODO +-- * take into account a non-centered sample. +-- * find a way to use integrateMulti with a small amount of memory. +-- * better mask for each detector. + +-- | Samples + +project :: FilePath +project = "/nfs/ruche-diffabs/diffabs-soleil/com-diffabs/" + +published :: FilePath +published = project </> "2016" </> "Run4B" </> "OutilsMetallo_CarolineHamon" + +sampleRef :: XRDRef +sampleRef = XRDRef "reference" + (published </> "xrd" </> "calibration") + (XrdRefNxs + (mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_30.nxs") "scan_30" h5path') + 33 + ) + +h5path' :: NxEntry -> DataFrameH5Path XrdOneD +h5path' nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_02" + gamma = "D13-1-CX1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/actuator_1_1" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +sampleCalibration :: XRDCalibration Xpad32 +sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" + , xrdCalibrationOutputDir = published </> "xrd" </> "calibration" -- TODO pourquoi ce output + , xrdCalibrationDetector = Xpad32 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + + idxs :: [Int] + idxs = [5, 33, 100, 246, 300, 436] + + entry :: Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryNxs + { xrdCalibrationEntryNxs'Nxs = mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_30.nxs") "scan_30" h5path' + , xrdCalibrationEntryNxs'Idx = idx + , xrdCalibrationEntryNxs'NptPath = published </> "xrd" </> "calibration" </> printf "IHR_30.nxs_%02d.npt" idx + } + + entries :: [XRDCalibrationEntry] + entries = [ entry idx | idx <- idxs] + + +bins :: DIM1 +bins = ix1 1000 + +multibins :: DIM1 +multibins = ix1 10000 + +threshold :: Maybe Threshold +threshold = Just (Threshold 5000) + +skipedFrames :: [Int] +skipedFrames = [] + +ceo2 :: XRDSample +ceo2 = XRDSample "CeO2" + (published </> "xrd" </> "CeO2") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_29.nxs") "scan_29" h5path' + , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_30.nxs") "scan_30" h5path' + , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_56.nxs") "scan_56" h5path' + , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_58.nxs") "scan_58" h5path' + ] + ] + +-- | Main + +hamon :: IO () +hamon = do + -- | pre-calibrate (extract from nexus to edf in order to do the + -- calibration) + extractEdf sampleCalibration + + p <- getPoniExtRef sampleRef + + let poniextref = move p (Pose (MyMatrix HklB (ident 3))) + + -- full calibration + poniextref' <- calibrate sampleCalibration poniextref + + print poniextref + print poniextref' + + -- Integrate the flyscan mesh + -- 4.680504680504681e-3 per images (2*60+18) / 29484 this contain + -- read/write and computation + -- integrateMesh (XrdMeshParams poniextref' mflat method) [fly] + + -- | set the integration parameters + let mflat = Nothing + let aiMethod = Csr + let params = XrdOneDParams poniextref' mflat aiMethod + + -- integrate each step of the scan + integrate params [ceo2] + + -- this code doesn not work because there is not enought memory on + -- the computer. + -- integrateMulti params [ceo2] + + return () diff --git a/contrib/haskell/src/Hkl/Projects/Diffabs/Hercules.hs b/contrib/haskell/src/Hkl/Projects/Diffabs/Hercules.hs new file mode 100644 index 0000000..7c43650 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Diffabs/Hercules.hs @@ -0,0 +1,168 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Projects.Diffabs.Hercules + ( hercules ) where + +import Data.Array.Repa (DIM1, ix1) +import Numeric.Units.Dimensional.Prelude ((*~), centi, degree, meter) +import System.FilePath ((</>)) +import Text.Printf (printf) + +import Prelude hiding (lookup, readFile, writeFile) + +import Hkl + +-- | Samples + +project ∷ FilePath +project = "/nfs/ruche-diffabs/diffabs-soleil/com-diffabs/" + +published ∷ FilePath +published = "/nfs/ruche-diffabs/diffabs-soleil/com-diffabs/2017/Run2B/TPHercules" + +-- | Calibration part + +mkNxs' ∷ FilePath → Int → (NxEntry → DataFrameH5Path a ) → Nxs a +mkNxs' d idx h = mkNxs f' e h + where + f ∷ FilePath → Int → (FilePath, NxEntry) + f d' i' = (d' </> printf "scan_%d.nxs" i', printf "scan_%d" i') + + (f', e) = f d idx + +h5path ∷ NxEntry → DataFrameH5Path XrdOneD +h5path nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_03" + gamma = "D13-1-CX1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/actuator_1_1" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +sampleRef ∷ XRDRef +sampleRef = XRDRef "reference" + (published </> "calibration") + (XrdRefNxs + (mkNxs' (project </> "2017" </> "Run2" </> "2017-03-21") 91 h5path) + 15 -- BEWARE only the 6th poni was generated with the right Xpad_flat geometry. + ) + +sampleCalibration ∷ XRDCalibration ImXpadS140 +sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" + , xrdCalibrationOutputDir = published </> "calibration" + , xrdCalibrationDetector = ImXpadS140 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + idxs ∷ [Int] + idxs = [15, 16, 17, 18, 19] + + entry ∷ Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryNxs + { xrdCalibrationEntryNxs'Nxs = mkNxs' (project </> "2017" </> "Run2" </> "2017-03-21") 91 h5path + , xrdCalibrationEntryNxs'Idx = idx + , xrdCalibrationEntryNxs'NptPath = published </> "calibration" </> printf "scan_91.nxs_%02d.npt" idx + } + + entries ∷ [XRDCalibrationEntry] + entries = map entry idxs + +-- | Data treatment + +bins ∷ DIM1 +bins = ix1 3000 + +multibins ∷ DIM1 +multibins = ix1 25000 + +threshold ∷ Maybe Threshold +threshold = Just (Threshold 800) + +skipedFrames ∷ [Int] +skipedFrames = [] + +-- Flat + +-- flat ∷ [Nxs XrdFlat] +-- flat = [mkNxs' (project </> "2017" </> "Run1" </> "2017-02-15") idx h5path | idx ← [57, 60 ∷ Int]] -- skip 58 59 for now (problème de droits d'accès) +-- where +-- h5path :: NxEntry -> DataFrameH5Path XrdFlat +-- h5path nxentry = XrdFlatH5Path (DataItemH5 (nxentry </> "scan_data/data_02") StrictDims) + +-- Scan en delta + +mkXRDSample ∷ String → [(FilePath, [Int])] -> XRDSample +mkXRDSample n ps = XRDSample n + (published </> "xrd" </> n) + [ XrdNxs bins multibins threshold skipedFrames n' | n' ← concatMap nxs''' ps ] + where + nxs''' ∷ (FilePath, [Int]) → [XrdSource] + nxs''' (p, idxs) = [XrdSourceNxs (mkNxs' p idx h5path) | idx ← idxs] + + +samples :: [XRDSample] +samples = map (uncurry mkXRDSample) + [ ("CeO2", [ ((project </> "2017" </> "Run2" </> "2017-03-21"), [91 :: Int]) ]) + , ("zgso4_room", [ ((project </> "2017" </> "Run2" </> "2017-03-21"), [96 :: Int]) ]) + , ("zgso4_450C", [ ((project </> "2017" </> "Run2" </> "2017-03-21"), [192 :: Int]) ]) + , ("zgso4_heating", [ ((project </> "2017" </> "Run2" </> "2017-03-21"), [100..190 :: Int]) ]) + , ("zgso4_cooling", [ ((project </> "2017" </> "Run2" </> "2017-03-21"), [199..214 :: Int]) ]) + ] + +-- | Main + +hercules ∷ IO () +hercules = do + + -- | pre-calibrate (extract from nexus to edf in order to do the + -- calibration) + extractEdf sampleCalibration + + -- | compute the flat + -- flat' ← computeFlat flat (published </> "flat" </> "flat.npy") + + -- | get a first ref poniExt + p ← getPoniExtRef sampleRef + -- set the initial position of the poni (pyFAI calibration is not + -- accurate with only one ring) + let poniextref = set p + (63 *~ centi meter) -- distance + (0 *~ meter) -- poni1 + (0 *~ meter) -- poni2 + (0 *~ degree) -- rot1 + (0 *~ degree) -- rot2 + (0 *~ degree) -- rot3 + print poniextref + + -- | full calibration + poniextref' ← calibrate sampleCalibration poniextref + print poniextref' + + -- | set the integration parameters + let mflat = Nothing + let aiMethod = Csr + let params = XrdOneDParams poniextref' mflat aiMethod + + -- -- integrate scan with multi geometry + -- -- splitPixel (the only available now) → 17m47.825s + integrateMulti params samples + + -- -- Integrate each image of the scans + -- -- Lut → 21.52 minutes + -- -- Csr → 21.9 minutes + -- integrate params samples + + -- -- substrack the air from all samples + -- substract params air samples + -- substractMulti params air samples + + return () diff --git a/contrib/haskell/src/Hkl/Projects/Diffabs/IRDRx.hs b/contrib/haskell/src/Hkl/Projects/Diffabs/IRDRx.hs new file mode 100644 index 0000000..86b7a68 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Diffabs/IRDRx.hs @@ -0,0 +1,158 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} + +module Hkl.Projects.Diffabs.IRDRx + ( irdrx ) where + +import Data.Array.Repa (DIM1, ix1) +import Numeric.LinearAlgebra (ident) +import System.FilePath ((</>)) +import Text.Printf (printf) + +import Prelude hiding (concat, lookup, readFile, writeFile) + +import Hkl + +-- | Samples + +project :: FilePath +project = "/nfs/ruche-diffabs/diffabs-soleil/com-diffabs/" + +published :: FilePath +published = project </> "2016" </> "Run5B" </> "irdrx" + +sampleRef :: XRDRef +sampleRef = XRDRef "reference" + (published </> "calibration") + (XrdRefNxs + (mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_39.nxs") "scan_39" h5path') + 10 + ) + +h5path' :: NxEntry -> DataFrameH5Path XrdOneD +h5path' nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_05" + gamma = "D13-1-CX1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/data_03" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +sampleCalibration :: XRDCalibration ImXpadS140 +sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" + , xrdCalibrationOutputDir = published </> "calibration" -- TODO pourquoi ce output + , xrdCalibrationDetector = ImXpadS140 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + + idxs :: [Int] + idxs = [0, 1, 10, 30] + + entry :: Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryNxs + { xrdCalibrationEntryNxs'Nxs = mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_39.nxs") "scan_39" h5path' + , xrdCalibrationEntryNxs'Idx = idx + , xrdCalibrationEntryNxs'NptPath = published </> "calibration" </> printf "scan_39.nxs_%02d.npt" idx + } + + entries :: [XRDCalibrationEntry] + entries = [ entry idx | idx <- idxs] + + +bins :: DIM1 +bins = ix1 1000 + +multibins :: DIM1 +multibins = ix1 10000 + +threshold :: Maybe Threshold +threshold = Just (Threshold 5000) + +skipedFrames :: [Int] +skipedFrames = [] + +lab6 :: XRDSample +lab6 = XRDSample "LaB6" + (published </> "LaB6") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_39.nxs") "scan_39" h5path' + , mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_40.nxs") "scan_40" h5path' + , mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_41.nxs") "scan_41" h5path' + , mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_42.nxs") "scan_42" h5path' + , mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_43.nxs") "scan_43" h5path' + , mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_44.nxs") "scan_44" h5path' + , mkNxs (project </> "2016" </> "Run5" </> "2016-11-09" </> "scan_45.nxs") "scan_45" h5path' + ] + ] + + + +-- meshSample :: String +-- meshSample :: project </> "2016" </> Run5 </> "2016-11-fly" </> "scan5 </> "*" +-- h5path nxentry = exptest_01368 +-- scan_data, sxpos szpos xpad_image 12x273 x 10 (fichiers) +-- delta = -6.2 +-- gamma = 0.0 +-- nrj 18.2 keV +fly :: XrdMeshSample +fly = XrdMeshSample "scan5" + (published </> "scan5") + [ XrdMesh bins multibins threshold + ( XrdMeshSourceNxsFly [mkNxs (project </> "2016" </> "Run5" </> "2016-11-fly" </> "scan5" </> printf "flyscan_%05d.nxs" n) "exptest_01368" h5path | + n <- [7087, 7088, 7089, 7090, 7091, 7092, 7093, 7094, 7095] :: [Int] + ] + ) + ] + where + h5path :: NxEntry -> (DataFrameH5Path XrdMesh) + h5path nxentry = XrdMeshFlyH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> meshx) StrictDims) + (DataItemH5 (nxentry </> meshy) StrictDims) + (DataItemConst gamma) + (DataItemConst delta) + (DataItemConst wavelength) + + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/xpad_image" + meshx = "scan_data/sxpos" + meshy = "scan_data/szpos" + gamma = 0.0 / 180.0 * 3.14159 + delta = -6.2 / 180.0 * 3.14159 + wavelength = 1.54 -- TODO vérifier + +-- | Main + +irdrx :: IO () +irdrx = do + let mflat = Nothing + let method = CsrOcl + + p <- getPoniExtRef sampleRef + + let poniextref = move (Hkl.flip p) (Pose (MyMatrix HklB (ident 3))) + + -- full calibration + poniextref' <- calibrate sampleCalibration poniextref + + print poniextref' + + -- Integrate the flyscan mesh + -- 4.680504680504681e-3 per images (2*60+18) / 29484 this contain + -- read/write and computation + integrateMesh (XrdMeshParams poniextref' mflat method) [fly] + + -- integrate each step of the scan + -- _ <- mapConcurrently (integrate poniextref') [lab6] + return () diff --git a/contrib/haskell/src/Hkl/Projects/Diffabs/Laure.hs b/contrib/haskell/src/Hkl/Projects/Diffabs/Laure.hs new file mode 100644 index 0000000..05706c6 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Diffabs/Laure.hs @@ -0,0 +1,206 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Projects.Diffabs.Laure + ( laure ) where + +import Data.Array.Repa (DIM1, ix1) +import Numeric.LinearAlgebra (ident) +import System.FilePath ((</>)) +import Text.Printf (printf) + +import Prelude hiding (lookup, readFile, writeFile) + +import Hkl + +-- | TODO +-- * Livre 45 p159 +-- * simplify with the list of nxs using list comprehension. +-- * add the flyscan mesh +-- * add possibility to sum a bunch of pixel coordinates from a mesh. on a mask + +-- | Samples + +project ∷ FilePath +project = "/nfs/ruche-diffabs/diffabs-users/20160370/" + +published ∷ FilePath +published = project </> "published-data" + +h5path ∷ NxEntry → DataFrameH5Path XrdOneD +h5path nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_02" + gamma = "D13-1-CX1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/actuator_1_1" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +mkNxs' ∷ FilePath → Int → (NxEntry → DataFrameH5Path a ) → Nxs a +mkNxs' d idx h = mkNxs f' e h + where + f ∷ FilePath → Int → (FilePath, NxEntry) + f d' i' = (d' </> printf "scan_%d.nxs" i', printf "scan_%d" (i' - 1)) + + (f', e) = f d idx + +-- | Calibration part + +sampleRef ∷ XRDRef +sampleRef = XRDRef "reference" + (published </> "calibration") + (XrdRefNxs + (mkNxs' (published </> "calibration") 45 h5path) + 10 -- BEWARE only the 6th poni was generated with the right Xpad_flat geometry. + ) + +sampleCalibration ∷ XRDCalibration ImXpadS140 +sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" + , xrdCalibrationOutputDir = published </> "calibration" + , xrdCalibrationDetector = ImXpadS140 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + idxs ∷ [Int] + idxs = [00, 01, 02, 03, 04, 09, 10, 11, 12, 14, 15, 18, 19, 22, 23, 26, 29, 33, 38, 42, 49, 53] + + entry ∷ Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryNxs + { xrdCalibrationEntryNxs'Nxs = mkNxs' (published </> "calibration") 45 h5path + , xrdCalibrationEntryNxs'Idx = idx + , xrdCalibrationEntryNxs'NptPath = published </> "calibration" </> printf "scan_45.nxs_%02d.npt" idx + } + + entries ∷ [XRDCalibrationEntry] + entries = map entry idxs + +-- | Data treatment + +bins ∷ DIM1 +bins = ix1 3000 + +multibins ∷ DIM1 +multibins = ix1 25000 + +threshold ∷ Maybe Threshold +threshold = Just (Threshold 800) + +skipedFrames ∷ [Int] +skipedFrames = [4] + +-- Flat + +flat ∷ [Nxs XrdFlat] +flat = [mkNxs' (project </> "2017" </> "Run1" </> "2017-02-15") idx h5path' | idx ← [57, 60 ∷ Int]] -- skip 58 59 for now (problème de droits d'accès) + where + h5path' :: NxEntry -> DataFrameH5Path XrdFlat + h5path' nxentry = XrdFlatH5Path (DataItemH5 (nxentry </> "scan_data/data_02") StrictDims) + +-- Scan en delta + +mkXRDSample ∷ String → [(FilePath, [Int])] -> XRDSample +mkXRDSample n ps = XRDSample n + (published </> "xrd" </> n) + [ XrdNxs bins multibins threshold skipedFrames n' | n' ← concatMap nxs''' ps ] + where + nxs''' ∷ (FilePath, [Int]) → [XrdSource] + nxs''' (p, idxs) = [XrdSourceNxs (mkNxs' p idx h5path) | idx ← idxs] + + +air ∷ XRDSample +air = mkXRDSample "air" [ ((project </> "2017" </> "Run1" </> "2017-02-17"), [198 :: Int]) ] + +samples :: [XRDSample] +samples = air : map (uncurry mkXRDSample) + [ ("CeO2", [ ((project </> "2017" </> "Run1" </> "2017-02-15"), [45 :: Int]) ]) + , ("kapton", [ ((project </> "2017" </> "Run1" </> "2017-02-17"), [197 :: Int]) ]) + , ("chlorite", [ ((project </> "2017" </> "Run1" </> "2017-02-15"), [53 :: Int]) ]) + , ("dMnO2", [ ((project </> "2017" </> "Run1" </> "2017-02-16"), [135 :: Int]) ]) + , ("bulk_L2", [ ((project </> "2017" </> "Shutdown1-2" </> "2017-02-19"), [315..316 :: Int]) ]) + , ("L1-H_3", [ ((project </> "2017" </> "Run1" </> "2017-02-15"), concat [ [62..63 :: Int] + , [65..70 :: Int] + , [74, 75 :: Int] + ]) + , ((project </> "2017" </> "Run1" </> "2017-02-16"), [76..89 :: Int]) + ]) + , ("L1-H_4", [ ((project </> "2017" </> "Run1" </> "2017-02-15"), [71..73 :: Int]) + , ((project </> "2017" </> "Run1" </> "2017-02-16"), concat [ [90..94 :: Int] + , [96..103 :: Int] + , [119..127 :: Int] + ]) + ]) + , ("L1-H_5", [ ((project </> "2017" </> "Run1" </> "2017-02-16"), [104..118 :: Int]) ]) + , ("L1-Patine_1", [ ((project </> "2017" </> "Run1" </> "2017-02-16"), [136..151 :: Int]) + , ((project </> "2017" </> "Run1" </> "2017-02-17"), concat [ [152..184 :: Int] + , [186 :: Int] + ]) + ]) + , ("L1-Patine_2", [ ((project </> "2017" </> "Run1" </> "2017-02-17"), [187..196 :: Int]) ]) + , ("L2-H_1", [ ((project </> "2017" </> "Run1" </> "2017-02-17"), [199..213 :: Int]) ]) + , ("L2-H_2", [ ((project </> "2017" </> "Run1" </> "2017-02-17"), [214..220 :: Int]) + , ((project </> "2017" </> "Run1" </> "2017-02-18"), concat [ [221..228 :: Int] + , [259..262 :: Int] + ]) + ]) + , ("L2-H_3", [ ((project </> "2017" </> "Run1" </> "2017-02-18"), [229..248 :: Int]) ]) + , ("L2-PatineFoncee", [ ((project </> "2017" </> "Run1" </> "2017-02-18"), [249..258 :: Int]) ]) + , ("L2-PatineFonceeNew", [ ((project </> "2017" </> "Run1" </> "2017-02-18"), concat [ [263, 264, 266, 267 :: Int] + , [269..273 :: Int]]) + ]) + , ("L2-patineLabo_1", [ ((project </> "2017" </> "Shutdown1-2" </> "2017-02-19"),[295..313 :: Int]) ]) + , ("L2-PatineClaire_1", [ ((project </> "2017" </> "Shutdown1-2" </> "2017-02-19"), [317..324 :: Int]) + , ((project </> "2017" </> "Shutdown1-2" </> "2017-02-20"), [325..356 :: Int]) + ]) + , ("L3-patine_1", [ ((project </> "2017" </> "Run1" </> "2017-02-19"), [274..293 :: Int]) + , ((project </> "2017" </> "Shutdown1-2" </> "2017-02-19"), [294, 295 :: Int]) + ]) + ] + +-- | Main + +laure ∷ IO () +laure = do + + -- | compute the flat + flat' ← computeFlat flat (published </> "flat" </> "flat.npy") + + -- | get a first ref poniExt + p ← getPoniExtRef sampleRef + -- flip the ref poni in order to fit the reality + -- let poniextref = p + let poniextref = move p (Pose (MyMatrix HklB (ident 3))) + -- let poniextref = setPose (Hkl.PyFAI.PoniExt.flip p) (MyMatrix HklB (ident 3)) + print poniextref + + -- | full calibration + poniextref' ← calibrate sampleCalibration poniextref + print poniextref' + + -- | set the integration parameters + let mflat = Just flat' + let aiMethod = Csr + let params = XrdOneDParams poniextref' mflat aiMethod + + -- integrate scan with multi geometry + -- splitPixel (the only available now) → 17m47.825s + integrateMulti params samples + + -- Integrate each image of the scans + -- Lut → 21.52 minutes + -- Csr → 21.9 minutes + integrate params samples + + -- substrack the air from all samples + substract params air samples + substractMulti params air samples + + return () diff --git a/contrib/haskell/src/Hkl/Projects/Diffabs/Martinetto.hs b/contrib/haskell/src/Hkl/Projects/Diffabs/Martinetto.hs new file mode 100644 index 0000000..977d9b3 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Diffabs/Martinetto.hs @@ -0,0 +1,294 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} + +module Hkl.Projects.Diffabs.Martinetto + ( martinetto + , martinetto' + ) where + +import Data.Array.Repa (DIM1, ix1) +import Numeric.LinearAlgebra (ident) +import System.FilePath ((</>)) +import Text.Printf (printf) + +import Prelude hiding (concat, lookup, readFile, writeFile) + +import Hkl + +-- | Samples + +project :: FilePath +project = "/nfs/ruche-diffabs/diffabs-users/99160066/" + +published :: FilePath +published = project </> "published-data" + +h5path' :: NxEntry -> DataFrameH5Path XrdOneD +h5path' nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_53" + gamma = "d13-1-cx1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/actuator_1_1" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +sampleCalibration :: XRDCalibration Xpad32 +sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" + , xrdCalibrationOutputDir = published </> "calibration" + , xrdCalibrationDetector = Xpad32 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + + idxs :: [Int] + idxs = [3, 6, 9, 15, 18, 21, 24, 27, 30, 33, 36, 39, 43] + + entry :: Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryNxs + { xrdCalibrationEntryNxs'Nxs = mkNxs (published </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path' + , xrdCalibrationEntryNxs'Idx = idx + , xrdCalibrationEntryNxs'NptPath = published </> "calibration" </> printf "XRD18keV_26.nxs_%02d.npt" idx + } + + entries :: [XRDCalibrationEntry] + entries = [ entry idx | idx <- idxs] + + +sampleRef :: XRDRef +sampleRef = XRDRef "reference" + (published </> "calibration") + (XrdRefNxs + (mkNxs (published </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path') + 6 -- BEWARE only the 6th poni was generated with the right Xpad_flat geometry. + ) + +h5path :: NxEntry -> DataFrameH5Path XrdOneD +h5path nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_58" + gamma = "D13-1-CX1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/actuator_1_1" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +bins :: DIM1 +bins = ix1 8000 + +multibins :: DIM1 +multibins = ix1 25000 + +threshold :: Maybe Threshold +threshold = Just (Threshold 800) + +skipedFrames :: [Int] +skipedFrames = [] + +ceo2 :: XRDSample +ceo2 = XRDSample "CeO2" + (published </> "CeO2") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (published </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path' ] + ] + +n27t2 :: XRDSample +n27t2 = XRDSample "N27T2" + (published </> "N27T2") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "N27T2_14.nxs") "scan_14" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "N27T2_17.nxs") "scan_17" h5path + ] + ] + +r23 :: XRDSample +r23 = XRDSample "R23" + (published </> "R23") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R23_6.nxs") "scan_6" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R23_12.nxs") "scan_12" h5path + ] + ] + +r18 :: XRDSample +r18 = XRDSample "R18" + (published </> "R18") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R18_20.nxs") "scan_20" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R18_24.nxs") "scan_24" h5path + ] + ] + +a3 :: XRDSample +a3 = XRDSample "A3" + (published </> "A3") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "A3_13.nxs") "scan_13" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "A3_14.nxs") "scan_14" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "A3_15.nxs") "scan_15" h5path + ] + ] + +a2 :: XRDSample +a2 = XRDSample "A2" + (published </> "A2") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "A2_14.nxs") "scan_14" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "A2_17.nxs") "scan_17" h5path + ] + ] + +a26 :: XRDSample +a26 = XRDSample "A26" + (published </> "A26") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_50.nxs") "scan_50" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_51.nxs") "scan_51" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_52.nxs") "scan_52" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_53.nxs") "scan_53" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_54.nxs") "scan_54" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_55.nxs") "scan_55" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_56.nxs") "scan_56" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_57.nxs") "scan_57" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_58.nxs") "scan_58" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "A26_59.nxs") "scan_59" h5path + ] + ] + +d2 :: XRDSample +d2 = XRDSample "D2" + (published </> "D2") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "D2_16.nxs") "scan_16" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "D2_17.nxs") "scan_17" h5path + ] + ] + +d3 :: XRDSample +d3 = XRDSample "D3" + (published </> "D3") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "D3_14.nxs") "scan_14" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "D3_15.nxs") "scan_15" h5path + ] + ] + +f30 :: XRDSample +f30 = XRDSample "F30" + (published </> "F30") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "F30_11.nxs") "scan_11" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "F30_12.nxs") "scan_12" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "F30_13.nxs") "scan_13" h5path + ] + ] + +r11 :: XRDSample +r11 = XRDSample "R11" + (published </> "R11") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R11_5.nxs") "scan_5" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R11_6.nxs") "scan_6" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R11_7.nxs") "scan_7" h5path + ] + ] + +d16 :: XRDSample +d16 = XRDSample "D16" + (published </> "D16") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "D16_12.nxs") "scan_12" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "D16_15.nxs") "scan_15" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "D16_17.nxs") "scan_17" h5path + ] + ] + +k9a2 :: XRDSample +k9a2 = XRDSample "K9A2" + (published </> "K9A2") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "K9A2_1_31.nxs") "scan_31" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "K9A2_1_32.nxs") "scan_32" h5path + ] + ] + +r34n1 :: XRDSample +r34n1 = XRDSample "R34N1" + (published </> "R34N1") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R34N1_28.nxs") "scan_28" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-27" </> "R34N1_37.nxs") "scan_37" h5path + ] + ] + +r35n1 :: XRDSample +r35n1 = XRDSample "R35N1" + (published </> "R35N1") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "R35N1_25.nxs") "scan_19" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "R35N1_26.nxs") "scan_20" h5path + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-26" </> "R35N1_27.nxs") "scan_21" h5path + ] + ] + +-- meshSample :: String +-- meshSample = project </> "2016" "Run2" "2016-03-28" "MELLE_29.nxs" +-- scan_29 scan_data actuator_1_1 actuator_2_1 data_58 (images) + +-- | Main + +martinetto :: IO () +martinetto = do + -- lire le ou les ponis de référence ainsi que leur géométrie + -- associée. + + -- let samples = [ceo2, a2, a3, a26, d2, d3, d16, f30, k9a2, n27t2, r11, r18, r23, r34n1, r35n1] + let samples = [ceo2] + + p <- getPoniExtRef sampleRef + + -- flip the ref poni in order to fit the reality + -- let poniextref = Hkl.PyFAI.PoniExt.flip p + let poniextref = p + -- integrate each step of the scan + let params = XrdOneDParams poniextref Nothing Lut + integrate params samples + + -- plot de la figure. (script python ou autre ?) + return () + +martinetto' :: IO () +martinetto' = do + let samples = [ceo2, a2, a3, a26, d2, d3, d16, f30, k9a2, n27t2, r11, r18, r23, r34n1, r35n1] + let mflat = Nothing + + p <- getPoniExtRef sampleRef + + -- flip the ref poni in order to fit the reality + -- let poniextref = p + let poniextref = move p (Pose (MyMatrix HklB (ident 3))) + -- let poniextref = setPose (Hkl.PyFAI.PoniExt.flip p) (MyMatrix HklB (ident 3)) + + -- full calibration + poniextref' <- calibrate sampleCalibration poniextref + -- print p + print poniextref + print poniextref' + + -- integrate each step of the scan + integrateMulti (XrdOneDParams poniextref' mflat Csr) samples + + return () diff --git a/contrib/haskell/src/Hkl/Projects/Diffabs/Melle.hs b/contrib/haskell/src/Hkl/Projects/Diffabs/Melle.hs new file mode 100644 index 0000000..de837f0 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Diffabs/Melle.hs @@ -0,0 +1,439 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Projects.Diffabs.Melle + ( melle ) where + +-- import Control.Concurrent (setNumCapabilities) +-- import Control.Concurrent.Async (mapConcurrently) +import Data.Array.Repa (DIM1, ix1) +-- import Data.Char (toUpper) +-- import Numeric.LinearAlgebra (ident) +import System.FilePath ((</>)) +import Text.Printf (printf) + +-- import Prelude hiding (concat, lookup, readFile, writeFile) + +import Hkl + +published ∷ FilePath +published = "/nfs/ruche-diffabs/diffabs-soleil/com-diffabs/Reguer/USERSexperiences/melle" + +-- | TODO + +-- | MELLE / VIALAS +-- Session 1 MACRO - 16-17 février 2016 (Logbook n° 42 p 169) +-- Session 2 MICRO KB --28 mars 2016 (Logbook 42 + Logbook 43 p3) +-- Session 3 MICRO pinhole - 22-24 juillet 2016 (Logbook 44 p33) +-- Session 4 MACRO - septembre 2016 (Logbook 44 p63) + +-- | Session 1 + +-- macrofaisceau +-- 16keV +-- Λ = 0,775 +-- detection : XPAD S140 / image = data 54 +-- sample : ω = 5 et χ = 70 + +-- calibration = beam direct + +-- - 3 MESH pour 3 positions du détecteur de diffraction (delta = -4, 3, 10), + +-- macro python: +-- for i in range (10): +-- myx = -12+i*0,5 +-- mv(samplex, myx) +-- ascan(sampley, -8, 12, 100, 10) + +-- scan_26 à 55.nxs +-- diffabs-soleil/com-diffabs/2016/Run1/2016-02-16 ou 02-17 + +-- 2THETA = 1 DELTA SCAN +-- scan_56 = ascan(delta, -4, 70, 18, 3) +-- scan_58 = ascan(delta, -4, 70, 18, 3) + + +-- | Session 2 + +-- microbeam +-- 18keV, ?= 0,6888Å +-- detection : XPAD 3.2 / image = data 58 +-- sample : ? = 5° et ? = 80°. +-- calibration CeO2 +-- data dans le dossier du proposal de Philippe Charlier 2015 1386 +-- voir aussi script Martinetto proposal IHR 99160066 +-- scan_25 = ascan(delta, -14.5, 60.5, 75, 0.5) +-- scan_26 = ascan(delta, -14, 60, 75, 1) +-- scan_27 = ascan(delta, -14, 60, 46, 1) + +-- MESH : MELLE_29.nxs +-- dossier: diffabs-soleil/com-diffabs/2016/Run2/2016-03-28 + +-- calibration + +project2 :: FilePath +project2 = "/nfs/ruche-diffabs/diffabs-users/99160066/" + +published2:: FilePath +published2 = project2 </> "published-data" + +h5path2 :: NxEntry -> DataFrameH5Path XrdOneD +h5path2 nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_53" + gamma = "d13-1-cx1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/actuator_1_1" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +sampleCalibration2 :: XRDCalibration Xpad32 +sampleCalibration2 = XRDCalibration { xrdCalibrationName = "calibration2" + , xrdCalibrationOutputDir = published </> "calibration2" + , xrdCalibrationDetector = Xpad32 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + idxs :: [Int] + idxs = [3, 6, 9, 15, 18, 21, 24, 27, 30, 33, 36, 39, 43] + + entry :: Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryNxs + { xrdCalibrationEntryNxs'Nxs = mkNxs (published2 </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path2 + , xrdCalibrationEntryNxs'Idx = idx + , xrdCalibrationEntryNxs'NptPath = published2 </> "calibration" </> printf "XRD18keV_26.nxs_%02d.npt" idx + } + + entries :: [XRDCalibrationEntry] + entries = [ entry idx | idx <- idxs] + +sampleRef2 :: XRDRef +sampleRef2 = XRDRef "reference" + (published2 </> "calibration") + (XrdRefNxs + (mkNxs (published2 </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path2) + 6 -- BEWARE only the 6th poni was generated with the right Xpad_flat geometry. + ) + +bins :: DIM1 +bins = ix1 8000 + +multibins :: DIM1 +multibins = ix1 25000 + +threshold :: Maybe Threshold +threshold = Just (Threshold 800) + +skipedFrames :: [Int] +skipedFrames = [] + +melleScan :: XRDSample +melleScan = XRDSample "CeO2" + (published </> "xrd" </> "session2" </> "oned") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run2" </> "2016-03-23" </> "XRD18keV_25.nxs") "scan_25" h5path2 + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-23" </> "XRD18keV_26.nxs") "scan_26" h5path2 + , mkNxs (project </> "2016" </> "Run2" </> "2016-03-23" </> "XRD18keV_27.nxs") "scan_27" h5path2 + ] + ] + where + project ∷ FilePath + project = "/nfs/ruche-diffabs/diffabs-users/20151386/" + + +melleMesh :: XrdMeshSample +melleMesh = XrdMeshSample "MELLE_29" + (published </> "xrd" </> "session2" </> "mesh") + [ XrdMesh bins multibins threshold (XrdMeshSourceNxs n) | n <- + [ mkNxs (project2' </> "2016" </> "Run2" </> "2016-03-28" </> "MELLE_29.nxs") "scan_29" h5path2' + ] + ] + where + project2' :: FilePath + project2' = "/nfs/ruche-diffabs/diffabs-users/99160066/" + + h5path2' :: NxEntry -> DataFrameH5Path XrdMesh + h5path2' nxentry = + XrdMeshH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> meshX) StrictDims) + (DataItemH5 (nxentry </> meshY) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> beamline </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_58" + meshX = "scan_data/actuator_1_1" + meshY = "scan_data/actuator_2_1" + gamma = "D13-1-CX1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "D13-1-CX1__EX__DIF.1-DELTA__#1/raw_value" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + + +session2 :: IO () +session2 = do + -- compute the ref poni + p ← getPoniExtRef sampleRef2 + poniextref <- calibrate sampleCalibration2 p + + -- integrate the mesh + let mflat = Nothing + integrateMesh (XrdMeshParams poniextref mflat CsrOcl) [melleMesh] + + -- integrate the scan parts + let params = XrdOneDParams poniextref mflat Csr + integrate params [melleScan] + integrateMulti params [melleScan] + + return () + +-- | session 4 +-- macro +-- 18keV, ?= 0,6888Å +-- detection : XPAD 3.2 + +session4 ∷ IO () +session4 = do + -- calibration + p ← getPoniExtRef sampleRef + poniextref <- calibrate sampleCalibration p + +-- calibration : CeO2 +-- On peut utiliser la calib de IHR_30, mais il faut prendre en compte le décentrage. +-- IHR_56 +-- IHR_58 +-- sont deux autres possibilité de calibration. +-- diffabs-soleil\com-diffabs\2016\Run4\2016-09-07 + + -- | set the integration parameters + let mflat = Nothing + let params = XrdOneDParams poniextref mflat Csr + + -- integrate each step of the scan + integrate params [ceo2] + +-- 1 seul "MESH"(20, 49) à partir d'une serie 2THETA +-- IHR_63 à 95 +-- diffabs-soleil\com-diffabs\2016\Run4\2016-09-07 +-- IHR_96 à 190 +-- diffabs-soleil\com-diffabs\2016\Run4\2016-09-08 +-- obtenu via la macro suivante. +-- for i in range(20): +-- myx = -11 + i +-- mv(txs, myx) # exhantillon à 45 degree donc ce double déplacement correspond au vrai x +-- mv(tys, myx) +-- for j in range(29): +-- myy = 12 + j +-- mv(tabV, myy) +-- ascan(δ, -13.6, 30, 109, 5) + + return () + + where + + project :: FilePath + project = "/nfs/ruche-diffabs/diffabs-soleil/com-diffabs/" + + published' :: FilePath + published' = project </> "2016" </> "Run4B" </> "OutilsMetallo_CarolineHamon" + + sampleRef :: XRDRef + sampleRef = XRDRef "reference" + (published' </> "xrd" </> "calibration") + (XrdRefNxs + (mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_30.nxs") "scan_30" h5path') + 33 + ) + + h5path' :: NxEntry -> DataFrameH5Path XrdOneD + h5path' nxentry = + XrdOneDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemH5 (nxentry </> beamline </> gamma) ExtendDims) + (DataItemH5 (nxentry </> delta) ExtendDims) + (DataItemH5 (nxentry </> beamline </> wavelength) StrictDims) + where + beamline :: String + beamline = beamlineUpper Diffabs + + image = "scan_data/data_02" + gamma = "D13-1-CX1__EX__DIF.1-GAMMA__#1/raw_value" + delta = "scan_data/actuator_1_1" + wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + + sampleCalibration :: XRDCalibration Xpad32 + sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" + , xrdCalibrationOutputDir = published' </> "xrd" </> "calibration" -- TODO pourquoi ce output + , xrdCalibrationDetector = Xpad32 + , xrdCalibrationCalibrant = CeO2 + , xrdCalibrationEntries = entries + } + where + + idxs :: [Int] + idxs = [5, 33, 100, 246, 300, 436] + + entry :: Int -> XRDCalibrationEntry + entry idx = XRDCalibrationEntryNxs + { xrdCalibrationEntryNxs'Nxs = mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_30.nxs") "scan_30" h5path' + , xrdCalibrationEntryNxs'Idx = idx + , xrdCalibrationEntryNxs'NptPath = published' </> "xrd" </> "calibration" </> printf "IHR_30.nxs_%02d.npt" idx + } + + entries :: [XRDCalibrationEntry] + entries = [ entry idx | idx <- idxs] + + bins :: DIM1 + bins = ix1 1000 + + multibins :: DIM1 + multibins = ix1 10000 + + threshold :: Maybe Threshold + threshold = Just (Threshold 5000) + + skipedFrames :: [Int] + skipedFrames = [] + + ceo2 :: XRDSample + ceo2 = XRDSample "CeO2" + (published </> "session4" </> "xrd" </> "CeO2") + [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- + [ mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_29.nxs") "scan_29" h5path' + , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_30.nxs") "scan_30" h5path' + , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_56.nxs") "scan_56" h5path' + , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_58.nxs") "scan_58" h5path' + ] + ] + +-- ** session 5 +-- micro +-- 18.05keV +-- detection XPAD S140 + +-- calibration CeO2 +-- gam = 9 phi = 170 +-- 18p1kev_71 +-- gam = 9 phi = 175 +-- 18p1kev_73 +-- gam = 0 phi = 205 +-- 18p1kev_74 +-- gam = 0.3 phi = 205 +-- 18p1kev_75 +-- ruche-diffabs\diffabs-users\99170085\2017\Run3\2017-05-14 + +-- FLAT (à verifier si suffisant) (faire la somme des trois fichiers) +-- 18p1kev_82 +-- 18p1kev_83 +-- 18p1kev_84 +-- ruche-diffabs\diffabs-users\99170085\2017\Run3\2017-05-14 + +-- FLY -- ???? +-- flyscan_16602 +-- diffabs-soleil\com-diffabs\2017\Run3\fly_IHRSol + +-- 2THETA = 1 DELTA SCAN +-- 18p1kev_85 +-- 18p1kev_86 +-- ruche-diffabs\diffabs-users\99170085\2017\Run3\2017-05-14 + +-- | Samples + +-- published :: FilePath +-- published = project </> "published-data" + +-- beamlineUpper :: Beamline -> String +-- beamlineUpper b = [Data.Char.toUpper x | x <- show b] + +-- nxs :: FilePath -> NxEntry -> (NxEntry -> DataFrameH5Path) -> Nxs +-- nxs f e h = Nxs f e (h e) + +-- nxs' :: FilePath -> NxEntry -> (NxEntry -> a) -> Nxs' a +-- nxs' f e h = Nxs' f e (h e) + +-- h5path :: NxEntry -> DataFrameH5Path +-- h5path nxentry = +-- DataFrameH5Path { h5pImage = DataItem (nxentry </> image) StrictDims +-- , h5pGamma = DataItem (nxentry </> beamline </> gamma) ExtendDims +-- , h5pDelta = DataItem (nxentry </> delta) ExtendDims +-- , h5pWavelength = DataItem (nxentry </> beamline </> wavelength) StrictDims +-- } +-- where +-- beamline :: String +-- beamline = beamlineUpper Diffabs + +-- image = "scan_data/data_53" +-- gamma = "d13-1-cx1__EX__DIF.1-GAMMA__#1/raw_value" +-- delta = "scan_data/actuator_1_1" +-- wavelength = "D13-1-C03__OP__MONO__#1/wavelength" + +-- sampleCalibration :: XRDCalibration +-- sampleCalibration = XRDCalibration { xrdCalibrationName = "calibration" +-- , xrdCalibrationOutputDir = published </> "calibration" +-- , xrdCalibrationEntries = entries +-- } +-- where + +-- idxs :: [Int] +-- idxs = [3, 6, 9, 15, 18, 21, 24, 27, 30, 33, 36, 39, 43] + +-- entry :: Int -> XRDCalibrationEntry +-- entry idx = XRDCalibrationEntryNxs +-- { xrdCalibrationEntryNxs'Nxs = nxs (published </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path +-- , xrdCalibrationEntryNxs'Idx = idx +-- , xrdCalibrationEntryNxs'NptPath = published </> "calibration" </> printf "XRD18keV_26.nxs_%02d.npt" idx +-- } + +-- entries :: [XRDCalibrationEntry] +-- entries = [ entry idx | idx <- idxs] + + +-- sampleRef :: XRDRef +-- sampleRef = XRDRef "reference" +-- (published </> "calibration") +-- (nxs (published </> "calibration" </> "XRD18keV_26.nxs") "scan_26" h5path) +-- 6 -- BEWARE only the 6th poni was generated with the right Xpad_flat geometry. + +-- bins :: DIM1 +-- bins = ix1 8000 + +-- multibins :: DIM1 +-- multibins = ix1 25000 + +-- threshold :: Threshold +-- threshold = Threshold 800 + + +-- p <- getPoniExtRef sampleRef + +-- -- flip the ref poni in order to fit the reality +-- -- let poniextref = p +-- let poniextref = setPose p (MyMatrix HklB (ident 3)) +-- -- let poniextref = setPose (Hkl.PyFAI.PoniExt.flip p) (MyMatrix HklB (ident 3)) + +-- -- full calibration +-- poniextref' <- calibrate sampleCalibration poniextref Xpad32 +-- -- print p +-- print poniextref +-- print poniextref' + +-- -- integrate each step of the scan +-- _ <- mapM_ (integrateMesh poniextref') samples + +-- return () + +melle ∷ IO () +melle = do + session2 + session4 diff --git a/contrib/haskell/src/Hkl/Projects/Mars.hs b/contrib/haskell/src/Hkl/Projects/Mars.hs new file mode 100644 index 0000000..75b46d6 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Mars.hs @@ -0,0 +1,4 @@ +module Hkl.Projects.Mars (module X) where + +import Hkl.Projects.Mars.Schlegel as X +import Hkl.Projects.Mars.Romeden as X diff --git a/contrib/haskell/src/Hkl/Projects/Mars/Romeden.hs b/contrib/haskell/src/Hkl/Projects/Mars/Romeden.hs new file mode 100644 index 0000000..f89589d --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Mars/Romeden.hs @@ -0,0 +1,47 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Projects.Mars.Romeden + ( romeden ) where + +import Codec.Picture ( saveTiffImage ) +import Control.Arrow ( (&&&) ) +import System.FilePath ((</>)) +import System.FilePath.Glob ( compile, globDir1 ) + +import Prelude hiding (concat, lookup, readFile, writeFile) + +import Hkl + +-- | TODO +-- ne pas planter lorsque l'image est manquante dans une nx entry. + +project ∷ FilePath +-- project = "/nfs/ruche-mars/mars-soleil/com-mars/2017_Run2/comisioning_microfaisceau" +-- project = "/home/experiences/instrumentation/picca" +project = "/media/picca/Transcend/ROMEDENNE" + +h5path ∷ NxEntry → DataFrameH5Path XrdFlat +h5path nxentry = + XrdFlatH5Path + (DataItemH5 (nxentry </> image) StrictDims) + where + image ∷ H5Path + image = "image#0/data" + +saveAsTiff' ∷ (Nxs XrdFlat, FilePath) → IO () +saveAsTiff' (n, o) = saveTiffImage o =<< toTiff n + +saveAsTiff ∷ (NxEntry -> DataFrameH5Path XrdFlat) → FilePath → IO () +saveAsTiff h5path' n = mapM_ (saveAsTiff' . (nxs &&& out)) =<< nxEntries n + where + nxs ∷ FilePath → Nxs XrdFlat + nxs nx = mkNxs (project </> n) nx h5path' + + out ∷ FilePath → FilePath + out nx = (project </> n) ++ nx ++ ".tiff" + +-- | Main + +romeden ∷ IO () +romeden = mapM_ (saveAsTiff h5path) =<< globDir1 (compile "*.nxs") project diff --git a/contrib/haskell/src/Hkl/Projects/Mars/Schlegel.hs b/contrib/haskell/src/Hkl/Projects/Mars/Schlegel.hs new file mode 100644 index 0000000..6ada48a --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Mars/Schlegel.hs @@ -0,0 +1,110 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Projects.Mars.Schlegel + ( schlegel ) where + +import Numeric.LinearAlgebra (ident) +import System.FilePath ((</>)) +import Text.Printf (printf) + +import Prelude hiding (concat, lookup, readFile, writeFile) + +import Hkl + +-- | TODO +-- * check if the +-- * find a way to use integrateMulti with a small amount of memory. +-- * better mask for each detector. + +-- | Samples + +project :: FilePath +project = "/nfs/share-temp/picca/20160800" + +published :: FilePath +published = project </> "published-data" + +h5path :: NxEntry -> DataFrameH5Path XrdZeroD +h5path nxentry = + XrdZeroDH5Path + (DataItemH5 (nxentry </> image) StrictDims) + (DataItemConst 0.0485945) + where + image ∷ H5Path + image = "scan_data/data_01" + +sampleCalibration ∷ XrdZeroDCalibration Xpad32 +sampleCalibration = XrdZeroDCalibration (XrdZeroDSample name outputdir entries) Xpad32 LaB6 + where + name ∷ String + name = "lab6" + + outputdir ∷ AbsDirPath + outputdir = published </> "xrd" </> "calibration" + + entries :: [XrdZeroDSource] + entries = [ XrdZeroDSourceNxs $ + mkNxs (project </> "2017" </> "Run3" </> "scan_5_01.nxs") "_5" h5path + ] + + +-- bins :: DIM1 +-- bins = ix1 1000 + +-- multibins :: DIM1 +-- multibins = ix1 10000 + +-- threshold :: Maybe Threshold +-- threshold = Just (Threshold 5000) + +-- skipedFrames :: [Int] +-- skipedFrames = [] + +-- ceo2 :: XRDSample +-- ceo2 = XRDSample "CeO2" +-- (published </> "xrd" </> "CeO2") +-- [ XrdNxs bins multibins threshold skipedFrames (XrdSourceNxs n) | n <- +-- [ mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_29.nxs") "scan_29" h5path' +-- , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_30.nxs") "scan_30" h5path' +-- , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_56.nxs") "scan_56" h5path' +-- , mkNxs (project </> "2016" </> "Run4" </> "2016-09-07" </> "IHR_58.nxs") "scan_58" h5path' +-- ] +-- ] + +-- | Main + +schlegel :: IO () +schlegel = do + -- | pre-calibrate (extract from nexus to edf in order to do the + -- calibration) + extractEdf sampleCalibration + + -- p <- getPoniExtRef sampleRef + + -- let poniextref = move p (Pose (MyMatrix HklB (ident 3))) + + -- -- full calibration + -- poniextref' <- calibrate sampleCalibration poniextref + + -- print poniextref + -- print poniextref' + + -- -- Integrate the flyscan mesh + -- -- 4.680504680504681e-3 per images (2*60+18) / 29484 this contain + -- -- read/write and computation + -- -- integrateMesh (XrdMeshParams poniextref' mflat method) [fly] + + -- -- | set the integration parameters + -- let mflat = Nothing + -- let aiMethod = Csr + -- let params = XrdOneDParams poniextref' mflat aiMethod + + -- -- integrate each step of the scan + -- integrate params [ceo2] + + -- -- this code doesn not work because there is not enought memory on + -- -- the computer. + -- -- integrateMulti params [ceo2] + + return () diff --git a/contrib/haskell/src/Hkl/Projects/Sixs.hs b/contrib/haskell/src/Hkl/Projects/Sixs.hs new file mode 100644 index 0000000..1c6cdb5 --- /dev/null +++ b/contrib/haskell/src/Hkl/Projects/Sixs.hs @@ -0,0 +1,141 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE GADTs #-} +module Hkl.Projects.Sixs + ( main_sixs ) + where + +#if __GLASGOW_HASKELL__ < 710 +import Control.Applicative ((<$>), (<*>)) +#endif + +import Data.ByteString.Char8 (pack) +import Data.Vector.Storable (concat, head) +import Control.Exception (bracket) +import Control.Monad (forM_) +import Numeric.LinearAlgebra (Matrix) +import Numeric.Units.Dimensional.Prelude (meter, nano, (*~)) +import Pipes (Producer, runEffect, (>->), lift, yield) +import Pipes.Prelude (print) +import System.FilePath.Posix ((</>)) + +import Hkl ( DataItem ( DataItemH5 ) + , Dataset + , ExtendDims ( ExtendDims, StrictDims ) + , Factory(Uhv) + , File + , Geometry(Geometry) + , H5 + , Source(Source) + , check_ndims + , closeDataset + , get_position + , get_ub + , lenH5Dataspace + , openDataset + , withH5File + ) + +{-# ANN module "HLint: ignore Use camelCase" #-} + +data DataFrameHklH5Path + = DataFrameHklH5Path + (DataItem H5) -- Image + (DataItem H5) -- Mu + (DataItem H5) -- Omega + (DataItem H5) -- delta + (DataItem H5) -- gamma + (DataItem H5) -- UB + (DataItem H5) -- Wavelength + (DataItem H5) -- DiffractometerType + deriving (Show) + +data DataFrameHklH5 + = DataFrameHklH5 + Dataset -- image + Dataset -- mu + Dataset -- omega + Dataset -- delta + Dataset -- gamma + Dataset -- ub + Dataset -- wavelength + Dataset -- dtype + +data DataFrame + = DataFrame + Int -- n + Geometry -- geometry + (Matrix Double) -- ub + deriving (Show) + +withDataframeH5 :: File -> DataFrameHklH5Path -> (DataFrameHklH5 -> IO r) -> IO r +withDataframeH5 h5file dfp = bracket (hkl_h5_open h5file dfp) hkl_h5_close + +hkl_h5_open :: File -> DataFrameHklH5Path -> IO DataFrameHklH5 +hkl_h5_open h5file (DataFrameHklH5Path i m o d g u w t) = DataFrameHklH5 + <$> openDataset' h5file i + <*> openDataset' h5file m + <*> openDataset' h5file o + <*> openDataset' h5file d + <*> openDataset' h5file g + <*> openDataset' h5file u + <*> openDataset' h5file w + <*> openDataset' h5file t + where + openDataset' :: File -> DataItem H5 -> IO Dataset + openDataset' hid (DataItemH5 name _) = openDataset hid (pack name) Nothing + +hkl_h5_is_valid :: DataFrameHklH5 -> IO Bool +hkl_h5_is_valid (DataFrameHklH5 _ m o d g _ _ _) = do + True <- check_ndims m 1 + True <- check_ndims o 1 + True <- check_ndims d 1 + True <- check_ndims g 1 + return True + +hkl_h5_close :: DataFrameHklH5 -> IO () +hkl_h5_close (DataFrameHklH5 i m o d g u w t) = do + closeDataset i + closeDataset m + closeDataset o + closeDataset d + closeDataset g + closeDataset u + closeDataset w + closeDataset t + +getDataFrame' :: DataFrameHklH5 -> Int -> IO DataFrame +getDataFrame' (DataFrameHklH5 _ m o d g u w _) i = do + mu <- get_position m i + omega <- get_position o i + delta <- get_position d i + gamma <- get_position g i + wavelength <- get_position w 0 + ub <- get_ub u + let positions = Data.Vector.Storable.concat [mu, omega, delta, gamma] + let source = Source (Data.Vector.Storable.head wavelength *~ nano meter) + return $ DataFrame i (Geometry Uhv source positions Nothing) ub + +getDataFrame :: DataFrameHklH5 -> Producer DataFrame IO () +getDataFrame d@(DataFrameHklH5 _ m _ _ _ _ _ _) = do + (Just n) <- lift $ lenH5Dataspace m + forM_ [0..n-1] (\i -> lift (getDataFrame' d i) >>= yield) + +main_sixs :: IO () +main_sixs = do + let root = "/nfs/ruche-sixs/sixs-soleil/com-sixs/2015/Shutdown4-5/XpadAu111/" + let filename = "align_FLY2_omega_00045.nxs" + let dataframe_h5p = DataFrameHklH5Path + (DataItemH5 "com_113934/scan_data/xpad_image" StrictDims) + (DataItemH5 "com_113934/scan_data/UHV_MU" ExtendDims) + (DataItemH5 "com_113934/scan_data/UHV_OMEGA" ExtendDims) + (DataItemH5 "com_113934/scan_data/UHV_DELTA" ExtendDims) + (DataItemH5 "com_113934/scan_data/UHV_GAMMA" ExtendDims) + (DataItemH5 "com_113934/SIXS/I14-C-CX2__EX__DIFF-UHV__#1/UB" StrictDims) + (DataItemH5 "com_113934/SIXS/Monochromator/wavelength" StrictDims) + (DataItemH5 "com_113934/SIXS/I14-C-CX2__EX__DIFF-UHV__#1/type" StrictDims) + + withH5File (root </> filename) $ \h5file -> + withDataframeH5 h5file dataframe_h5p $ \dataframe_h5 -> do + True <- hkl_h5_is_valid dataframe_h5 + runEffect $ getDataFrame dataframe_h5 + >-> Pipes.Prelude.print diff --git a/contrib/haskell/src/Hkl/PyFAI.hs b/contrib/haskell/src/Hkl/PyFAI.hs new file mode 100644 index 0000000..eeed0e9 --- /dev/null +++ b/contrib/haskell/src/Hkl/PyFAI.hs @@ -0,0 +1,9 @@ +module Hkl.PyFAI (module X) where + +import Hkl.PyFAI.AzimuthalIntegrator as X +import Hkl.PyFAI.Calib as X +import Hkl.PyFAI.Calibrant as X +import Hkl.PyFAI.Detector as X +import Hkl.PyFAI.Poni as X +import Hkl.PyFAI.PoniExt as X +import Hkl.PyFAI.Npt as X diff --git a/contrib/haskell/src/Hkl/PyFAI/AzimuthalIntegrator.hs b/contrib/haskell/src/Hkl/PyFAI/AzimuthalIntegrator.hs new file mode 100644 index 0000000..e29df7d --- /dev/null +++ b/contrib/haskell/src/Hkl/PyFAI/AzimuthalIntegrator.hs @@ -0,0 +1,18 @@ +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.PyFAI.AzimuthalIntegrator + ( AIMethod(..) + ) where + +data AIMethod = Numpy | Cython | SplitPixel | Lut | Csr | NoSplitCsr | FullCsr | LutOcl | CsrOcl + +instance Show AIMethod where + show Numpy = "numpy" + show Cython = "cython" + show SplitPixel = "splitpixel" + show Lut = "lut" + show Csr = "csr" + show NoSplitCsr = "nosplit_csr" + show FullCsr = "full_csr" + show LutOcl = "lut_ocl" + show CsrOcl = "csr_ocl" diff --git a/contrib/haskell/src/Hkl/PyFAI/Calib.hs b/contrib/haskell/src/Hkl/PyFAI/Calib.hs new file mode 100644 index 0000000..1c41a09 --- /dev/null +++ b/contrib/haskell/src/Hkl/PyFAI/Calib.hs @@ -0,0 +1,29 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.PyFAI.Calib + ( ToPyFAICalibArg(..) ) where + +import Data.Text (unpack) +import Numeric.Units.Dimensional.Prelude ((/~), nano, meter) + +import Hkl.Types ( WaveLength ) +import Hkl.Detector ( Detector ) +import Hkl.PyFAI.Calibrant ( Calibrant ) +import Hkl.PyFAI.Detector ( toPyFAI ) + +class ToPyFAICalibArg a where + toPyFAICalibArg ∷ a → String + +instance ToPyFAICalibArg FilePath where + toPyFAICalibArg f = f + +instance ToPyFAICalibArg (Detector a) where + toPyFAICalibArg d = "-D" ++ unpack (toPyFAI d) + +instance ToPyFAICalibArg Calibrant where + toPyFAICalibArg c = "-c " ++ show c + +instance ToPyFAICalibArg WaveLength where + toPyFAICalibArg w = "-w " ++ show ((w /~ nano meter) * 10) diff --git a/contrib/haskell/src/Hkl/PyFAI/Calibrant.hs b/contrib/haskell/src/Hkl/PyFAI/Calibrant.hs new file mode 100644 index 0000000..f6bd110 --- /dev/null +++ b/contrib/haskell/src/Hkl/PyFAI/Calibrant.hs @@ -0,0 +1,10 @@ +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.PyFAI.Calibrant + ( Calibrant(..) ) where + +data Calibrant = CeO2 | LaB6 + +instance Show Calibrant where + show CeO2 = "CeO2" + show LaB6 = "LaB6" diff --git a/contrib/haskell/src/Hkl/PyFAI/Detector.hs b/contrib/haskell/src/Hkl/PyFAI/Detector.hs new file mode 100644 index 0000000..9c7e172 --- /dev/null +++ b/contrib/haskell/src/Hkl/PyFAI/Detector.hs @@ -0,0 +1,19 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.PyFAI.Detector + ( ToPyFAI(..) + ) where + +import Data.Text (Text) + +import Hkl.Detector ( Detector ( Xpad32, ImXpadS140, ZeroD ) ) + +class ToPyFAI a where + toPyFAI ∷ a → Text + +instance ToPyFAI (Detector a) where + toPyFAI Xpad32 = "Xpad_flat" + toPyFAI ImXpadS140 = "imxpad_s140" + toPyFAI ZeroD = error "Unsupported Detector" diff --git a/contrib/haskell/src/Hkl/PyFAI/Npt.hs b/contrib/haskell/src/Hkl/PyFAI/Npt.hs new file mode 100644 index 0000000..a7567cf --- /dev/null +++ b/contrib/haskell/src/Hkl/PyFAI/Npt.hs @@ -0,0 +1,99 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Hkl.PyFAI.Npt + ( Npt(..) + , NptEntry(..) + , NptPoint(..) + , nptP + , nptFromFile + ) where + +import Control.Applicative +import Data.Attoparsec.Text +import Data.Text +import Data.Text.IO (readFile) +import Numeric.Units.Dimensional.Prelude (Angle, Length, (*~), meter, radian) + +type Calibrant = Text + +data NptPoint = NptPoint { nptPointX :: Double + , nptPointY :: Double + } + deriving (Show) + +data NptEntry = NptEntry { nptEntryId :: Int + , nptEntryTth :: Angle Double + , nptEntryRing :: Int + , nptPoints :: [NptPoint] + } + deriving (Show) + +data Npt = Npt { nptComment :: [Text] + , nptCalibrant :: Calibrant + , nptWavelength :: Length Double + , npdDSpacing :: [Length Double] + , nptEntries :: [NptEntry] + } + deriving (Show) + +commentP :: Parser Text +commentP = "#" *> takeTill isEndOfLine <* endOfLine <?> "commentP" + +headerP :: Parser [Text] +headerP = many1 commentP <?> "headerP" + +calibrantP :: Parser Text +calibrantP = "calibrant: " *> takeTill isEndOfLine <* endOfLine <?> "calibrantP" + +dspacingP :: Parser [Length Double] +dspacingP = "dspacing:" *> many1 lengthP' <* endOfLine <?> "dspasingP" + +doubleP :: Text -> Parser Double +doubleP key = string key *> double <* endOfLine <?> "doubleP" + +lengthP' :: Parser (Length Double) +lengthP' = do + skipSpace + value <- double + pure $ value *~ meter + +lengthP :: Text -> Parser (Length Double) +lengthP key = do + value <- doubleP key + pure $ value *~ meter + +angleP :: Text -> Parser (Angle Double) +angleP key = do + value <-doubleP key + pure $ value *~ radian + +intP :: Text -> Parser Int +intP key = string key *> decimal <* endOfLine <?> "intP" + +nptPointP :: Parser NptPoint +nptPointP = NptPoint + <$> ("point: x=" *> double) + <*> (" y=" *> double <* endOfLine) + +nptEntryP :: Parser NptEntry +nptEntryP = NptEntry + <$> (skipSpace *> intP "New group of points: ") + <*> angleP "2theta: " + <*> intP "ring: " + <*> many nptPointP + +nptP :: Parser Npt +nptP = Npt + <$> headerP + <*> calibrantP + <*> lengthP "wavelength: " + <*> dspacingP + <*> many1 nptEntryP + <?> "nptP" + +nptFromFile :: FilePath -> IO Npt +nptFromFile filename = do + content <- Data.Text.IO.readFile filename + return $ case parseOnly nptP content of + Left _ -> error $ "Can not parse the " ++ filename ++ " npt file" + Right a -> a diff --git a/contrib/haskell/src/Hkl/PyFAI/Poni.hs b/contrib/haskell/src/Hkl/PyFAI/Poni.hs new file mode 100644 index 0000000..f8ec7eb --- /dev/null +++ b/contrib/haskell/src/Hkl/PyFAI/Poni.hs @@ -0,0 +1,257 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ExistentialQuantification #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.PyFAI.Poni + ( Pose(..) + -- Poni + , Poni + , PoniPath + , poniP + , poniToText + -- PoniEntry + , PoniEntry + , poniEntryFlip + , poniEntryFromList + , poniEntryRotation + , poniEntryTranslation + , poniEntryToList + , poniEntrySet + , poniEntryMove + -- other + , fromAxisAndAngle + ) where + +import Control.Applicative ((<$>), (<|>), (<*>), (*>), (<*), many, optional, pure) +import Data.Attoparsec.Text (Parser, (<?>), endOfLine, isEndOfLine, many1, double, string, takeTill) +import Data.Text (Text, append, intercalate, pack) +import Data.Vector.Storable (Vector, fromList) +import Numeric.LinearAlgebra (Matrix, (<>), atIndex, fromLists, ident, scalar) +import Numeric.Units.Dimensional.Prelude (Angle, Length, (+), (*~), (/~), (/~~), one, meter, radian, degree) + +import Hkl.Detector +import Hkl.MyMatrix +import Hkl.PyFAI.Detector +import Hkl.Types + +#if !MIN_VERSION_hmatrix(0, 17, 0) +import Numeric.LinearAlgebra (trans) +tr:: Matrix t -> Matrix t +tr = trans +#else +import Numeric.LinearAlgebra (tr) +#endif + +type PoniPath = FilePath + +-- | Pose + +data Pose = Pose (MyMatrix Double) deriving (Show) + +-- | ADetector + +data ADetector = forall a. ADetector (Detector a) + +instance Show ADetector where + show (ADetector v) = show v + +instance ToPyFAI ADetector where + toPyFAI (ADetector v) = toPyFAI v + +-- | Poni + +data PoniEntry = PoniEntry { poniEntryHeader :: [Text] + , poniEntryDetector :: (Maybe ADetector) -- ^ Detector Name + , poniEntryPixelSize1 :: (Length Double) -- ^ pixels size 1 + , poniEntryPixelSize2 :: (Length Double) -- ^ pixels size 1 + , poniEntryDistance :: (Length Double) -- ^ pixels size 2 + , poniEntryPoni1 :: (Length Double) -- ^ poni1 + , poniEntryPoni2 :: (Length Double) -- ^ poni2 + , poniEntryRot1 :: (Angle Double) -- ^ rot1 + , poniEntryRot2 :: (Angle Double) -- ^ rot2 + , poniEntryRot3 :: (Angle Double) -- ^ rot3 + , poniEntrySpline :: (Maybe Text) -- ^ spline file + , poniEntryWavelength :: WaveLength -- ^ wavelength + } + deriving (Show) + +type Poni = [PoniEntry] + +class ToPoni a where + toPoni ∷ a → Text + +instance ToPoni ADetector where + toPoni (ADetector v) = toPyFAI v + +instance ToPoni Double where + toPoni v = pack $ show v + +instance ToPoni Text where + toPoni = id + +commentP :: Parser Text +commentP = "#" *> takeTill isEndOfLine <* endOfLine <?> "commentP" + +headerP :: Parser [Text] +headerP = many1 commentP <?> "headerP" + +doubleP :: Text -> Parser Double +doubleP key = string key *> double <* endOfLine <?> "doubleP" + +lengthP :: Text -> Parser (Length Double) +lengthP key = do + value <-doubleP key + pure $ value *~ meter + +angleP :: Text -> Parser (Angle Double) +angleP key = do + value <-doubleP key + pure $ value *~ radian + +detectorP ∷ ToPyFAI a ⇒ a → Parser a +detectorP d = do + _ ← "Detector: " *> string (toPyFAI d) <* endOfLine + pure d + +aDetectorP ∷ Parser ADetector +aDetectorP = (ADetector <$> detectorP Xpad32) <|> (ADetector <$> detectorP ImXpadS140) + +poniEntryP :: Parser PoniEntry +poniEntryP = PoniEntry + <$> headerP + <*> optional aDetectorP + <*> lengthP "PixelSize1: " + <*> lengthP "PixelSize2: " + <*> lengthP "Distance: " + <*> lengthP "Poni1: " + <*> lengthP "Poni2: " + <*> angleP "Rot1: " + <*> angleP "Rot2: " + <*> angleP "Rot3: " + <*> optional ("SplineFile: " *> takeTill isEndOfLine <* endOfLine) + <*> lengthP "Wavelength: " + <?> "poniEntryP" + +poniP :: Parser Poni +poniP = many poniEntryP + +poniToText :: Poni -> Text +poniToText p = Data.Text.intercalate (Data.Text.pack "\n") (map poniEntryToText p) + +poniEntryToText :: PoniEntry -> Text +poniEntryToText p = intercalate (Data.Text.pack "\n") $ + map (Data.Text.append "#") (poniEntryHeader p) + ++ maybe [] (poniLine "Detector: ") (poniEntryDetector p) + ++ poniLine "PixelSize1: " (poniEntryPixelSize1 p /~ meter) + ++ poniLine "PixelSize2: " (poniEntryPixelSize2 p /~ meter) + ++ poniLine "Distance: " (poniEntryDistance p /~ meter) + ++ poniLine "Poni1: " (poniEntryPoni1 p /~ meter) + ++ poniLine "Poni2: " (poniEntryPoni2 p /~ meter) + ++ poniLine "Rot1: " (poniEntryRot1 p /~ radian) + ++ poniLine "Rot2: " (poniEntryRot2 p /~ radian) + ++ poniLine "Rot3: " (poniEntryRot3 p /~ radian) + ++ maybe [] (poniLine "SplineFile: ") (poniEntrySpline p) + ++ poniLine "Wavelength: " (poniEntryWavelength p /~ meter) + where + poniLine :: ToPoni a ⇒ String → a → [Text] + poniLine key v = [Data.Text.append (Data.Text.pack key) (toPoni v)] + +crossprod :: Vector Double -> Matrix Double +crossprod axis = fromLists [[ 0, -z, y], + [ z, 0, -x], + [-y, x, 0]] + where + x = axis `atIndex` 0 + y = axis `atIndex` 1 + z = axis `atIndex` 2 + +fromAxisAndAngle :: Vector Double -> Angle Double -> Matrix Double +fromAxisAndAngle axis angle = ident 3 Prelude.+ s * q Prelude.+ c * (q <> q) + where + c = scalar (1 - cos (angle /~ one)) + s = scalar (sin (angle /~ one)) + q = crossprod axis + +poniEntryFlip :: PoniEntry -> PoniEntry +poniEntryFlip p = p { poniEntryRot3 = new_rot3 } + where + rot3 = poniEntryRot3 p + new_rot3 = rot3 Numeric.Units.Dimensional.Prelude.+ 180 *~ degree + +poniEntryRotation :: PoniEntry -> Matrix Double -- TODO MyMatrix PyFAIB +poniEntryRotation e = Prelude.foldl (<>) (ident 3) rotations + where + rot1 = poniEntryRot1 e + rot2 = poniEntryRot2 e + rot3 = poniEntryRot3 e + rotations = Prelude.map (uncurry fromAxisAndAngle) + [ (fromList [0, 0, 1], rot3) + , (fromList [0, 1, 0], rot2) + , (fromList [1, 0, 0], rot1)] + +poniEntryTranslation :: PoniEntry -> Vector Double +poniEntryTranslation e = fromList ( [ poniEntryPoni1 e + , poniEntryPoni2 e + , poniEntryDistance e + ] /~~ meter ) + +poniEntryMove :: MyMatrix Double -> MyMatrix Double -> PoniEntry -> PoniEntry +poniEntryMove mym1 mym2 e = e { poniEntryRot1 = new_rot1 + , poniEntryRot2 = new_rot2 + , poniEntryRot3 = new_rot3 + } + where + rot1 = poniEntryRot1 e + rot2 = poniEntryRot2 e + rot3 = poniEntryRot3 e + rotations = Prelude.map (uncurry fromAxisAndAngle) + [ (Data.Vector.Storable.fromList [0, 0, 1], rot3) + , (Data.Vector.Storable.fromList [0, 1, 0], rot2) + , (Data.Vector.Storable.fromList [1, 0, 0], rot1)] + -- M1 . R0 = R1 + r1 = Prelude.foldl (<>) (ident 3) rotations -- pyFAIB + -- M2 . R0 = R2 + -- R2 = M2 . M1.T . R1 + r2 = Prelude.foldl (<>) m2 [tr m1, r1] + (new_rot1, new_rot2, new_rot3) = toEulerians r2 + + (MyMatrix _ m1) = changeBase mym1 PyFAIB + (MyMatrix _ m2) = changeBase mym2 PyFAIB + +poniEntrySet ∷ (Length Double) -- ^ distance + → (Length Double) -- ^ poni1 + → (Length Double) -- ^ poni2 + → (Angle Double) -- ^ rot1 + → (Angle Double) -- ^ rot2 + → (Angle Double) -- ^ rot3 + → PoniEntry + → PoniEntry +poniEntrySet d p1 p2 r1 r2 r3 p = + p { poniEntryDistance = d + , poniEntryPoni1 = p1 + , poniEntryPoni2 = p2 + , poniEntryRot1 = r1 + , poniEntryRot2 = r2 + , poniEntryRot3 = r3 + } + +poniEntryFromList :: PoniEntry -> [Double] -> PoniEntry +poniEntryFromList p [rot1, rot2, rot3, poni1, poni2, d] = + p { poniEntryDistance = d *~ meter + , poniEntryPoni1 = poni1 *~ meter + , poniEntryPoni2 = poni2 *~ meter + , poniEntryRot1 = rot1 *~ radian + , poniEntryRot2 = rot2 *~ radian + , poniEntryRot3 = rot3 *~ radian + } +poniEntryFromList _ _ = error "Can not convert to a PoniEntry" + +poniEntryToList :: PoniEntry -> [Double] +poniEntryToList p = [ poniEntryRot1 p /~ radian + , poniEntryRot2 p /~ radian + , poniEntryRot3 p /~ radian + , poniEntryPoni1 p /~ meter + , poniEntryPoni2 p /~ meter + , poniEntryDistance p /~ meter + ] diff --git a/contrib/haskell/src/Hkl/PyFAI/PoniExt.hs b/contrib/haskell/src/Hkl/PyFAI/PoniExt.hs new file mode 100644 index 0000000..63234f1 --- /dev/null +++ b/contrib/haskell/src/Hkl/PyFAI/PoniExt.hs @@ -0,0 +1,41 @@ +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.PyFAI.PoniExt + ( PoniExt(..) + , flip + , move + , set + ) where + +import Numeric.LinearAlgebra (ident) +import Numeric.Units.Dimensional.Prelude (Angle, Length) + +import Hkl.MyMatrix +import Hkl.PyFAI.Poni + +import Prelude hiding (flip) + +data PoniExt = PoniExt Poni Pose deriving (Show) + +flip :: PoniExt -> PoniExt +flip (PoniExt ps mym1) = PoniExt p mym1 + where + p = map poniEntryFlip ps + +set ∷ PoniExt + → (Length Double) -- ^ distance + → (Length Double) -- ^ poni1 + → (Length Double) -- ^ poni2 + → (Angle Double) -- ^ rot1 + → (Angle Double) -- ^ rot2 + → (Angle Double) -- ^ rot3 + → PoniExt +set (PoniExt ps _) d p1 p2 r1 r2 r3 = PoniExt p pose + where + p = map (poniEntrySet d p1 p2 r1 r2 r3) ps + pose = Pose (MyMatrix HklB (ident 3)) + +move :: PoniExt -> Pose -> PoniExt +move (PoniExt p1 (Pose mym1)) (Pose mym2) = PoniExt p (Pose mym2) + where + p = map (poniEntryMove mym1 mym2) p1 diff --git a/contrib/haskell/src/Hkl/Python.hs b/contrib/haskell/src/Hkl/Python.hs new file mode 100644 index 0000000..7ede8ae --- /dev/null +++ b/contrib/haskell/src/Hkl/Python.hs @@ -0,0 +1,30 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Python + ( PyVal(..) ) + where + +import Data.List (intercalate) + +class PyVal a where + toPyVal ∷ a → String + +instance PyVal a ⇒ PyVal (Maybe a) where + toPyVal (Just v) = toPyVal v + toPyVal Nothing = "None" + +instance PyVal String where + toPyVal s = show s + +instance PyVal [String] where + toPyVal vs = "[" ++ intercalate ",\n" (map toPyVal vs) ++ "]" + +instance PyVal Int where + toPyVal i = show i + +instance PyVal [Int] where + toPyVal is = "[" ++ intercalate ",\n" (map toPyVal is) ++ "]" + +instance PyVal Double where + toPyVal d = show d diff --git a/contrib/haskell/src/Hkl/Script.hs b/contrib/haskell/src/Hkl/Script.hs new file mode 100644 index 0000000..bffc3ec --- /dev/null +++ b/contrib/haskell/src/Hkl/Script.hs @@ -0,0 +1,107 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Script + ( Gnuplot + , Py2 + , Sh + , Script(..) + , run + , scriptRun + , scriptSave ) + where + +import Control.Monad (when) +import Data.Bits ((.|.)) +import Data.Text (Text) +import Data.Text.IO (writeFile) +import System.Directory (createDirectoryIfMissing) +import System.Exit ( ExitCode ( ExitSuccess ) ) +import System.FilePath ( (<.>), takeDirectory) +import System.Posix.Files (accessModes, groupModes, ownerModes, setFileMode) +import System.Posix.Types (FileMode) +import System.Process ( rawSystem ) -- callProcess for futur + +import Paths_hkl (getDataFileName) + +#if MIN_VERSION_directory(1, 3, 0) +import System.Directory (withCurrentDirectory) +#else +import Control.Exception.Base (bracket) +import System.Directory (getCurrentDirectory, setCurrentDirectory) +withCurrentDirectory :: FilePath -- ^ Directory to execute in + -> IO a -- ^ Action to be executed + -> IO a +withCurrentDirectory dir action = + bracket getCurrentDirectory setCurrentDirectory $ \ _ -> do + setCurrentDirectory dir + action +#endif + +type Profile = Bool + +data Gnuplot +data Py2 +data Sh + +data Script a where + Py2Script ∷ (Text, FilePath) → Script Py2 + ScriptGnuplot ∷ (Text, FilePath) → Script Gnuplot + ScriptSh ∷ (Text, FilePath) → Script Sh + +scriptSave' ∷ Text → FilePath → FileMode → IO () +scriptSave' c f m = do + createDirectoryIfMissing True (takeDirectory f) + Data.Text.IO.writeFile f c + setFileMode f m + print $ "--> created : " ++ f + +scriptSave ∷ Script a → IO () +scriptSave (Py2Script (c, f)) = scriptSave' c f (ownerModes .|. groupModes) +scriptSave (ScriptGnuplot (c, f)) = scriptSave' c f accessModes +scriptSave (ScriptSh (c, f)) = scriptSave' c f (ownerModes .|. groupModes) + +scriptRun' ∷ FilePath → String → [String] → Bool → IO ExitCode +scriptRun' f prog args d + | d == True = withCurrentDirectory directory go + | otherwise = go + where + go :: IO ExitCode + go = rawSystem prog args + + directory :: FilePath + directory = takeDirectory f + +scriptRun ∷ Script a → Bool → IO ExitCode +scriptRun (Py2Script (_, p)) d = do + ExitSuccess ← scriptRun' p "python" args d + when p' ( do + gprof2dot ← getDataFileName "data/gprof2dot.py" + ExitSuccess ← rawSystem gprof2dot ["-f", "pstats", stats, "-o", stats <.> "dot"] + ExitSuccess ← rawSystem dot ["-Tsvg", "-o", stats <.> "svg", stats <.> "dot"] + return () + ) + return ExitSuccess + where + -- BEWARE once actived the profiling multiply by two the computing time. + p' ∷ Profile + p' = True + + dot ∷ String + dot = "dot" + + stats ∷ String + stats = p <.> "pstats" + + args :: [String] + args + | p' == True = ["-m" , "cProfile", "-o", stats, p] + | otherwise = [p] +scriptRun (ScriptGnuplot (_, p)) d = scriptRun' p "gnuplot" [p] d +scriptRun (ScriptSh (_, p)) d = scriptRun' p p [] d + +run ∷ Script a → Bool → IO ExitCode +run s b = do + scriptSave s + scriptRun s b diff --git a/contrib/haskell/src/Hkl/Tiff.hs b/contrib/haskell/src/Hkl/Tiff.hs new file mode 100644 index 0000000..a604395 --- /dev/null +++ b/contrib/haskell/src/Hkl/Tiff.hs @@ -0,0 +1,10 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Tiff + ( ToTiff(..) ) where + +import Codec.Picture ( DynamicImage ) + +class ToTiff a where + toTiff ∷ a → IO DynamicImage diff --git a/contrib/haskell/src/Hkl/Types.hs b/contrib/haskell/src/Hkl/Types.hs new file mode 100644 index 0000000..adc56d0 --- /dev/null +++ b/contrib/haskell/src/Hkl/Types.hs @@ -0,0 +1,77 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Types ( AbsDirPath + , Beamline(..) + , beamlineUpper + , Mode(..) + , Engine(..) + , SampleName + , Sample(..) + , Source(..) + , Trajectory + , WaveLength + -- hdf5 + , H5Path + , module X + ) where + +import Data.Char (toUpper) + +import Hkl.Types.Parameter as X +import Hkl.H5 +import Hkl.Lattice +import Numeric.Units.Dimensional.Prelude (Length) + +-- Common + +type AbsDirPath = FilePath +type SampleName = String + +-- | Beamline + +data Beamline = Diffabs | Sixs + +instance Show Beamline where + show Diffabs = "diffabs" + show Sixs = "sixs" + +beamlineUpper ∷ Beamline → String +beamlineUpper b = [toUpper x | x ← show b] + +-- | Engine + +data Mode + = Mode + String -- ^ name + [Parameter] -- ^ parameters of the @Mode@ + deriving (Show) + +data Engine + = Engine + String -- ^ name + [Parameter] -- ^ pseudo axes values of the @Engine@ + Mode -- ^ current Mode + deriving (Show) + +-- | Sample + +data Sample a + = Sample + String -- ^ name of the sample + (Lattice a) -- ^ the lattice of the sample + Parameter -- ^ ux + Parameter -- ^ uy + Parameter -- ^ uz + deriving (Show) + +-- | Source + +type WaveLength = Length Double + +data Source = Source WaveLength deriving (Show) + +-- | Trajectory + +type Trajectory = [[Double]] diff --git a/contrib/haskell/src/Hkl/Types/Parameter.hsc b/contrib/haskell/src/Hkl/Types/Parameter.hsc new file mode 100644 index 0000000..e29ecde --- /dev/null +++ b/contrib/haskell/src/Hkl/Types/Parameter.hsc @@ -0,0 +1,85 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE CPP #-} + +module Hkl.Types.Parameter + ( Parameter(..) + , Range(..) + , copyParameter + , unit + ) where + +import Control.Monad (void) +import Foreign (nullPtr, Ptr, ForeignPtr, newForeignPtr, FunPtr) +import Foreign.Marshal.Alloc (alloca) +import Foreign.C ( CInt ( CInt ) + , CDouble ( CDouble ) + ) +import Foreign.C.String ( CString, peekCString ) +import Foreign.Storable ( Storable + , alignment + , sizeOf + , peek + , poke + ) + +#let alignment t = "%lu", (unsigned long)offsetof(struct {char x__; t (y__); }, y__) + +unit :: CInt +unit = 1 + +-- | Range + +data Range + = Range + Double -- ^ minimum value + Double -- ^ maximum value + deriving (Show) + +-- | Parameter + +data Parameter + = Parameter + String -- ^ name + Double -- ^ value + Range -- ^ range + deriving (Show) + +instance Storable Parameter where + alignment _ = #{alignment int} + sizeOf _ = #{size int} + peek ptr = alloca $ \pmin -> + alloca $ \pmax -> do + cname <- c_hkl_parameter_name_get ptr + name <- peekCString cname + value <- c_hkl_parameter_value_get ptr unit + c_hkl_parameter_min_max_get ptr pmin pmax unit + min_ <- peek pmin + max_ <- peek pmax + return (Parameter name value (Range min_ max_)) + poke ptr (Parameter _name value (Range min_ max_)) = do + void $ c_hkl_parameter_value_set ptr (CDouble value) unit nullPtr + void $ c_hkl_parameter_min_max_set ptr (CDouble min_) (CDouble max_) unit nullPtr + +copyParameter :: Ptr Parameter -> IO (ForeignPtr Parameter) +copyParameter p = newForeignPtr c_hkl_parameter_free =<< c_hkl_parameter_new_copy p + +foreign import ccall unsafe "hkl.h hkl_parameter_name_get" + c_hkl_parameter_name_get:: Ptr Parameter -> IO CString + +foreign import ccall unsafe "hkl.h hkl_parameter_value_get" + c_hkl_parameter_value_get:: Ptr Parameter -> CInt -> IO Double + +foreign import ccall unsafe "hkl.h hkl_parameter_min_max_get" + c_hkl_parameter_min_max_get :: Ptr Parameter -> Ptr Double -> Ptr Double -> CInt -> IO () + +foreign import ccall unsafe "hkl.h &hkl_parameter_free" + c_hkl_parameter_free :: FunPtr (Ptr Parameter -> IO ()) + +foreign import ccall unsafe "hkl.h hkl_parameter_new_copy" + c_hkl_parameter_new_copy:: Ptr Parameter -> IO (Ptr Parameter) + +foreign import ccall unsafe "hkl.h hkl_parameter_value_set" + c_hkl_parameter_value_set:: Ptr Parameter -> CDouble -> CInt -> Ptr () -> IO (CInt) + +foreign import ccall unsafe "hkl.h hkl_parameter_min_max_set" + c_hkl_parameter_min_max_set :: Ptr Parameter -> CDouble -> CDouble -> CInt -> Ptr () -> IO (CInt) diff --git a/contrib/haskell/src/Hkl/Utils.hs b/contrib/haskell/src/Hkl/Utils.hs new file mode 100644 index 0000000..130db10 --- /dev/null +++ b/contrib/haskell/src/Hkl/Utils.hs @@ -0,0 +1,17 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Utils + ( hasContent ) + where + +import Data.Text (Text) +import Data.Text.IO (writeFile) +import System.Directory (createDirectoryIfMissing) +import System.FilePath (takeDirectory) + +hasContent ∷ FilePath → Text → IO () +hasContent f c = do + createDirectoryIfMissing True (takeDirectory f) + Data.Text.IO.writeFile f c + print $ "--> created : " ++ f diff --git a/contrib/haskell/src/Hkl/Xrd.hs b/contrib/haskell/src/Hkl/Xrd.hs new file mode 100644 index 0000000..efc682c --- /dev/null +++ b/contrib/haskell/src/Hkl/Xrd.hs @@ -0,0 +1,6 @@ +module Hkl.Xrd ( module X ) where + +import Hkl.Xrd.Calibration as X +import Hkl.Xrd.OneD as X +import Hkl.Xrd.Mesh as X +import Hkl.Xrd.ZeroD as X diff --git a/contrib/haskell/src/Hkl/Xrd/Calibration.hs b/contrib/haskell/src/Hkl/Xrd/Calibration.hs new file mode 100644 index 0000000..30797cc --- /dev/null +++ b/contrib/haskell/src/Hkl/Xrd/Calibration.hs @@ -0,0 +1,355 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Xrd.Calibration + ( NptExt(..) + , XRDCalibrationEntry(..) + , XRDCalibration(..) + , calibrate + , extractEdf + ) where + +import Control.Applicative ((<$>), (<*>), pure) +import Control.Monad.IO.Class (liftIO) +import Data.ByteString.Char8 (pack) +import Data.List (foldl', intercalate) +import Data.Text (unlines, pack) +import Data.Vector.Storable + ( Vector + , head + , concat + , fromList + , slice + , toList + ) +import Numeric.LinearAlgebra + ( Matrix + , (<>) + , atIndex + , ident + ) +import Numeric.GSL.Minimization + ( MinimizeMethod ( NMSimplex2 ) + , minimizeV + ) +import Numeric.Units.Dimensional.Prelude (meter, radian, nano, (/~), (*~)) +import Pipes.Safe ( MonadSafe + , runSafeT, bracket + ) +import System.Exit ( ExitCode( ExitSuccess ) ) +import System.FilePath.Posix ((</>), takeFileName) +import Text.Printf ( printf ) + +-- import Hkl.C ( Geometry ( Geometry ) +-- , Factory ( K6c ) +-- , geometryDetectorRotationGet +-- ) +-- import Hkl.DataSource ( DataItem ( DataItemH5 ) ) +-- import Hkl.Detector ( Detector ( ZeroD ) +-- , coordinates +-- ) +-- import Hkl.Edf ( ExtractEdf() +-- , extractEdf +-- ) +-- import Hkl.H5 ( Dataset, File, H5 +-- , closeDataset +-- , get_position +-- , openDataset +-- , withH5File +-- ) +-- import Hkl.PyFAI ( Calibrant, Npt +-- , NptEntry ( NptEntry ) +-- , Poni +-- , PoniExt ( PoniExt ) +-- , Pose ( Pose ) +-- , fromAxisAndAngle +-- , nptEntries +-- , nptFromFile +-- , nptWavelength +-- , poniEntryFromList +-- , poniEntryToList +-- , toPyFAICalibArg +-- ) +-- import Hkl.Python ( toPyVal ) +-- import Hkl.MyMatrix ( Basis ( HklB, PyFAIB ) +-- , MyMatrix ( MyMatrix ) +-- , changeBase +-- ) +-- import Hkl.Nxs ( DataFrameH5Path ( XrdOneDH5Path ) +-- , Nxs ( Nxs ) +-- ) +-- import Hkl.Script ( Py2, Sh +-- , Script ( Py2Script, ScriptSh ) +-- , run +-- , scriptSave +-- ) +-- import Hkl.Types ( AbsDirPath, SampleName +-- , Source ( Source ) +-- , WaveLength ) +-- import Hkl.Xrd.OneD ( XrdOneD +-- , getPoseEdf +-- ) + +import Hkl.C +import Hkl.DataSource +import Hkl.Detector +import Hkl.Edf +import Hkl.H5 +import Hkl.PyFAI +import Hkl.Python +import Hkl.MyMatrix +import Hkl.Nxs +import Hkl.Script +import Hkl.Types +import Hkl.Xrd.OneD + +#if !MIN_VERSION_hmatrix(0, 17, 0) +(#>) :: Matrix Double -> Vector Double -> Vector Double +(#>) = (<>) +#else +import Numeric.LinearAlgebra ((#>)) +#endif + +-- | Calibration + +data NptExt a = NptExt { nptExtNpt :: Npt + , nptExtPose :: Pose + , nptExtDetector :: Detector a + } + deriving (Show) + +data XRDCalibrationEntry = XRDCalibrationEntryNxs { xrdCalibrationEntryNxs'Nxs :: Nxs XrdOneD + , xrdCalibrationEntryNxs'Idx :: Int + , xrdCalibrationEntryNxs'NptPath :: FilePath + } + | XRDCalibrationEntryEdf { xrdCalibrationEntryEdf'Edf :: FilePath + , xrdCalibrationEntryEdf'NptPath :: FilePath + } + deriving (Show) + +data XRDCalibration a = XRDCalibration { xrdCalibrationName :: SampleName + , xrdCalibrationOutputDir :: AbsDirPath + , xrdCalibrationDetector ∷ Detector a + , xrdCalibrationCalibrant ∷ Calibrant + , xrdCalibrationEntries :: [XRDCalibrationEntry] + } + deriving (Show) + +withDataItem :: MonadSafe m => File -> DataItem H5 -> (Dataset -> m r) -> m r +withDataItem hid (DataItemH5 name _) = bracket (liftIO acquire') (liftIO . release') + where + acquire' :: IO Dataset + acquire' = openDataset hid (Data.ByteString.Char8.pack name) Nothing + + release' :: Dataset -> IO () + release' = closeDataset + +getPoseNxs :: File -> DataFrameH5Path XrdOneD -> Int -> IO Pose -- TODO move to XRD +getPoseNxs f (XrdOneDH5Path _ g d w) i' = runSafeT $ + withDataItem f g $ \g' -> + withDataItem f d $ \d' -> + withDataItem f w $ \w' -> liftIO $ do + let mu = 0.0 + let komega = 0.0 + let kappa = 0.0 + let kphi = 0.0 + gamma <- get_position g' 0 + delta <- get_position d' i' + wavelength <- get_position w' 0 + let source = Source (Data.Vector.Storable.head wavelength *~ nano meter) + let positions = Data.Vector.Storable.concat [mu, komega, kappa, kphi, gamma, delta] + let geometry = Geometry K6c source positions Nothing + let detector = ZeroD + m <- geometryDetectorRotationGet geometry detector + return $ Pose (MyMatrix HklB m) + + +getWavelength ∷ File → DataFrameH5Path XrdOneD → IO WaveLength +getWavelength f (XrdOneDH5Path _ _ _ w) = runSafeT $ + withDataItem f w $ \w' -> liftIO $ do + wavelength <- get_position w' 0 + return $ Data.Vector.Storable.head wavelength *~ nano meter + +readWavelength :: XRDCalibrationEntry -> IO WaveLength +readWavelength e = + withH5File f $ \h5file -> getWavelength h5file p + where + (Nxs f p) = xrdCalibrationEntryNxs'Nxs e + + +readXRDCalibrationEntry :: Detector a -> XRDCalibrationEntry -> IO (NptExt a) +readXRDCalibrationEntry d e@(XRDCalibrationEntryNxs _ _ _) = + withH5File f $ \h5file -> NptExt + <$> nptFromFile (xrdCalibrationEntryNxs'NptPath e) + <*> getPoseNxs h5file p idx + <*> pure d + where + idx = xrdCalibrationEntryNxs'Idx e + (Nxs f p) = xrdCalibrationEntryNxs'Nxs e +readXRDCalibrationEntry d e@(XRDCalibrationEntryEdf _ _) = + NptExt + <$> nptFromFile (xrdCalibrationEntryEdf'NptPath e) + <*> getPoseEdf (xrdCalibrationEntryEdf'Edf e) + <*> pure d + +-- | Poni Calibration + +-- The minimized function is the quadratic difference of the +-- theoretical tth angle and for each pixel, the computed tth angle. + +-- synonyme types use in order to improve the calibration performance + +type NptEntry' = (Double, [Vector Double]) -- tth, detector pixels coordinates +type Npt' = (Double, [NptEntry']) -- wavelength, [NptEntry'] +type NptExt' a = (Npt', Matrix Double, Detector a) + +class ToGsl a where + toGsl ∷ a → Vector Double + +class FromGsl a where + fromGsl ∷ a → Vector Double → a + +class ToGslFunc a where + toGslFunc ∷ a → [NptExt b] → (Vector Double → Double) + +instance ToGsl PoniExt where + toGsl (PoniExt p _) = fromList $ poniEntryToList (last p) + +instance FromGsl PoniExt where + fromGsl (PoniExt p pose) v = PoniExt poni pose + where + poni ∷ Poni + poni = [poniEntryFromList (last p) (toList v)] + +instance ToGslFunc PoniExt where + toGslFunc _ npts = f (preCalibrate npts) + where + preCalibrate''' ∷ Detector a → NptEntry → NptEntry' + preCalibrate''' detector (NptEntry _ tth _ points) = (tth /~ radian, map (coordinates detector) points) + + preCalibrate'' ∷ Npt → Detector a → Npt' + preCalibrate'' n detector = (nptWavelength n /~ meter, map (preCalibrate''' detector) (nptEntries n)) + + preCalibrate' ∷ NptExt a → NptExt' a + preCalibrate' (NptExt n (Pose m) detector) = (preCalibrate'' n detector, m', detector) + where + (MyMatrix _ m') = changeBase m PyFAIB + + preCalibrate ∷ [NptExt a] → [NptExt' a] + preCalibrate = map preCalibrate' + + f :: [NptExt' a] → Vector Double → Double + f ns params = foldl' (f' rotation translation) 0 ns + where + rot1 = params `atIndex` 0 + rot2 = params `atIndex` 1 + rot3 = params `atIndex` 2 + + rotations = map (uncurry fromAxisAndAngle) + [ (fromList [0, 0, 1], rot3 *~ radian) + , (fromList [0, 1, 0], rot2 *~ radian) + , (fromList [1, 0, 0], rot1 *~ radian)] + + rotation = foldl' (<>) (ident 3) rotations + + translation :: Vector Double + translation = slice 3 3 params + + f' ∷ Matrix Double → Vector Double → Double → NptExt' a → Double + f' rotation translation x ((_wavelength, entries), m, _detector) = + foldl' (f'' translation r) x entries + where + r :: Matrix Double + r = m <> rotation + + f'' ∷ Vector Double → Matrix Double → Double → NptEntry' → Double + {-# INLINE f'' #-} + f'' translation r x (tth, pixels) = foldl' (f''' translation r tth) x pixels + + f''' ∷ Vector Double → Matrix Double → Double → Double → Vector Double → Double + {-# INLINE f''' #-} + f''' translation r tth x pixel = x + dtth * dtth + where + kf = r #> (pixel - translation) + x' = kf `atIndex` 0 + y' = kf `atIndex` 1 + z' = kf `atIndex` 2 + + dtth = tth - atan2 (sqrt (x'*x' + y'*y')) (-z') + +calibrate ∷ XRDCalibration a → PoniExt → IO PoniExt +calibrate (XRDCalibration _ _ d _ es) p = do + npts ← mapM (readXRDCalibrationEntry d) es + let guess = toGsl p + let f = toGslFunc p npts + let box = fromList [0.1, 0.1, 0.1, 0.01, 0.01, 0.01] + let (solution, _p) = minimizeV NMSimplex2 1E-16 3000 box f guess + print _p + return $ fromGsl p solution + +-- | Edf extraction before calibration + +edf ∷ AbsDirPath → FilePath → Int → FilePath +edf o n i = o </> f + where + f = (takeFileName n) ++ printf "_%02d.edf" i + +scriptExtractEdf ∷ AbsDirPath → [XRDCalibrationEntry] → Script Py2 +scriptExtractEdf o es = Py2Script (content, scriptPath) + where + content = Data.Text.unlines $ + map Data.Text.pack [ "#!/bin/env python" + , "" + , "from fabio.edfimage import edfimage" + , "from h5py import File" + , "" + , "NEXUSFILES = " ++ toPyVal nxss + , "IDXS = " ++ toPyVal idxs + , "IMAGEPATHS = " ++ toPyVal (imgs ∷ [String]) + , "OUTPUTS = " ++ toPyVal outputs + , "" + , "for filename, i, p, o in zip(NEXUSFILES, IDXS, IMAGEPATHS, OUTPUTS):" + , " with File(filename, mode='r') as f:" + , " edfimage(f[p][i]).write(o)" + ] + + (nxss, idxs, imgs) = unzip3 [(f, i, img) | (XRDCalibrationEntryNxs (Nxs f (XrdOneDH5Path (DataItemH5 img _) _ _ _)) i _) ← es] + + outputs ∷ [FilePath] + outputs = zipWith (edf o) nxss idxs + + scriptPath ∷ FilePath + scriptPath = o </> "pre-calibration.py" + +scriptPyFAICalib ∷ AbsDirPath → XRDCalibrationEntry → Detector a → Calibrant → WaveLength → Script Sh +scriptPyFAICalib o e d c w = ScriptSh (content, scriptPath) + where + content = Data.Text.unlines $ + map Data.Text.pack [ "#!/usr/bin/env sh" + , "" + , "pyFAI-calib " ++ intercalate " " args + ] + + args = [ toPyFAICalibArg w + , toPyFAICalibArg c + , toPyFAICalibArg d + , toPyFAICalibArg (edf o n i) ] + + (XRDCalibrationEntryNxs (Nxs n _) i _) = e + + scriptPath ∷ FilePath + scriptPath = o </> (takeFileName n) ++ printf "_%02d.sh" i + + +instance ExtractEdf (XRDCalibration a) where + extractEdf (XRDCalibration _ o d c es) = do + let script = scriptExtractEdf o es + ExitSuccess ← run script False + mapM_ go es + return () + where + go e = do + w ← readWavelength e + scriptSave $ scriptPyFAICalib o e d c w diff --git a/contrib/haskell/src/Hkl/Xrd/Mesh.hs b/contrib/haskell/src/Hkl/Xrd/Mesh.hs new file mode 100644 index 0000000..b99387a --- /dev/null +++ b/contrib/haskell/src/Hkl/Xrd/Mesh.hs @@ -0,0 +1,270 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Xrd.Mesh + ( XrdMeshSample(..) + , XrdMesh'(..) + , XrdMeshParams(..) + , XrdMeshSource(..) + , integrateMesh + ) where + +import Control.Concurrent.Async (mapConcurrently) +import Control.Monad (void) +import Control.Monad.Trans.Maybe (MaybeT, runMaybeT) +import Data.Array.Repa (Shape, DIM1, ix1, size) +import Data.Maybe (fromJust) +import Data.Vector.Storable (Vector, any, concat, head, singleton) +import Numeric.Units.Dimensional.Prelude (meter, nano, (/~), (*~)) +import System.Exit ( ExitCode( ExitSuccess ) ) +import System.FilePath ((</>), (<.>), dropExtension, splitDirectories, takeFileName) + +import qualified Data.Text as Text (unlines, pack) + +import Prelude hiding + ( any + , concat + , head + , lookup + , readFile + , unlines + ) +import Pipes ( lift ) + +import Hkl.C +import Hkl.DataSource +import Hkl.Detector +import Hkl.Flat +import Hkl.H5 +import Hkl.PyFAI +import Hkl.Python +import Hkl.MyMatrix +import Hkl.Nxs +import Hkl.Script +import Hkl.Types +import Hkl.Utils +import Hkl.Xrd.OneD + +-- | Types + +data XrdMeshSource = XrdMeshSourceNxs (Nxs XrdMesh) + | XrdMeshSourceNxsFly [Nxs XrdMesh] + deriving (Show) + +data XrdMesh' = XrdMesh DIM1 DIM1 (Maybe Threshold) XrdMeshSource deriving (Show) + +data XrdMeshSample = XrdMeshSample SampleName AbsDirPath [XrdMesh'] deriving (Show) -- ^ nxss + +data XrdMeshParams a = XrdMeshParams PoniExt (Maybe (Flat a)) AIMethod + +data XrdMeshFrame = XrdMeshFrame + WaveLength + Pose + deriving (Show) + +class FrameND t where + rowND :: t -> MaybeT IO XrdMeshFrame + +instance FrameND (DataFrameH5 XrdMesh) where + + rowND (XrdMeshH5 _ _ _ _ _ g d w) = do + let mu = 0.0 + let komega = 0.0 + let kappa = 0.0 + let kphi = 0.0 + gamma <- get_position' g (ix1 0) + delta <- get_position' d (ix1 0) + wavelength <- get_position' w (ix1 0) + let source@(Source w') = Source (head wavelength *~ nano meter) + let positions = concat [mu, komega, kappa, kphi, gamma, delta] + let geometry = Geometry K6c source positions Nothing + let detector = ZeroD + m <- lift $ geometryDetectorRotationGet geometry detector + let pose = Pose (MyMatrix HklB m) + return $ XrdMeshFrame w' pose + where + get_position' :: Shape sh => DataSource a -> sh -> MaybeT IO (Vector Double) + get_position' (DataSourceH5 _ a ) b = lift $ do + v <- get_position_new a b + if any isNaN v then fail "File contains Nan" else return v + get_position' (DataSourceConst v) _ = lift $ return $ singleton v + + rowND (XrdMeshFlyH5 _ _ _ _ _ g d w) = do + let mu = 0.0 + let komega = 0.0 + let kappa = 0.0 + let kphi = 0.0 + gamma <- get_position' g (ix1 0) + delta <- get_position' d (ix1 0) + wavelength <- get_position' w (ix1 0) + let source@(Source w') = Source (head wavelength *~ nano meter) + let positions = concat [mu, komega, kappa, kphi, gamma, delta] + let geometry = Geometry K6c source positions Nothing + let detector = ZeroD + m <- lift $ geometryDetectorRotationGet geometry detector + let pose = Pose (MyMatrix HklB m) + return $ XrdMeshFrame w' pose + where + get_position' :: Shape sh => DataSource a -> sh -> MaybeT IO (Vector Double) + get_position' (DataSourceH5 _ a ) b = lift $ do + v <- get_position_new a b + if any isNaN v then fail "File contains Nan" else return v + get_position' (DataSourceConst v) _ = lift $ return $ singleton v + +integrateMesh ∷ XrdMeshParams a → [XrdMeshSample] → IO () +integrateMesh p ss = void $ mapConcurrently (integrateMesh' p) ss + +integrateMesh' ∷ XrdMeshParams a → XrdMeshSample → IO () +integrateMesh' p (XrdMeshSample _ output nxss) = mapM_ (integrateMesh'' p output) nxss + +getWaveLengthAndPoniExt' ∷ XrdMeshParams a → Nxs XrdMesh → IO (WaveLength, PoniExt) +getWaveLengthAndPoniExt' (XrdMeshParams ref _ _) nxs = + withDataSource nxs $ \h -> do + -- read the first frame and get the poni used for all the integration. + d <- runMaybeT $ rowND h + let (XrdMeshFrame w p) = fromJust d + let poniext = move ref p + return (w, poniext) + +getWaveLengthAndPoniExt ∷ XrdMeshParams a → XrdMeshSource → IO (WaveLength, PoniExt) +getWaveLengthAndPoniExt p (XrdMeshSourceNxs nxs) = getWaveLengthAndPoniExt' p nxs +getWaveLengthAndPoniExt p (XrdMeshSourceNxsFly (nxs:_)) = getWaveLengthAndPoniExt' p nxs +getWaveLengthAndPoniExt _ (XrdMeshSourceNxsFly []) = error "getWaveLengthAndPoniExt" + +getOutputPath' ∷ AbsDirPath → FilePath → (FilePath, FilePath, FilePath) +getOutputPath' o d = (poni, h5, py) + where + poni = o </> d </> d <.> "poni" + h5 = o </> d </> d <.> "h5" + py = o </> d </> d <.> "py" + +getOutputPath ∷ AbsDirPath → XrdMeshSource → (FilePath, FilePath, FilePath) +getOutputPath o (XrdMeshSourceNxs (Nxs f _)) = getOutputPath' o dir + where + dir ∷ FilePath + dir = (dropExtension . takeFileName) f +getOutputPath o (XrdMeshSourceNxsFly (Nxs _ h:_)) = getOutputPath' o dir + where + (XrdMeshFlyH5Path (DataItemH5 i _) _ _ _ _ _) = h + dir:_ = splitDirectories i +getOutputPath _ (XrdMeshSourceNxsFly []) = error "getOutputPath" + + +xrdMeshPy'' ∷ Maybe (Flat a) + → AIMethod -- pyFAI azimuthal integration method + → [FilePath] -- nexus files + → H5Path -- image path + → H5Path -- meshx path + → H5Path -- meshy path + → FilePath -- ponipath + → DIM1 -- bins + → (Maybe Threshold) -- threshold + → WaveLength -- wavelength + → FilePath -- output h5 + → FilePath -- script name + → Script Py2 +xrdMeshPy'' mflat m fs i x y p b mt w o scriptPath = Py2Script (content, scriptPath) + where + content = Text.unlines $ + map Text.pack ["#!/bin/env python" + , "" + , "import itertools" + , "import numpy" + , "from h5py import File" + , "from pyFAI import load" + , "" + , "PONIFILE = " ++ toPyVal p + , "NEXUSFILES = " ++ toPyVal fs + , "MESHX = " ++ toPyVal x + , "MESHY = " ++ toPyVal y + , "IMAGEPATH = " ++ toPyVal i + , "N = " ++ toPyVal (size b) + , "OUTPUT = " ++ toPyVal o + , "WAVELENGTH = " ++ toPyVal (w /~ meter) + , "" + , "# Load the flat" + , "flat = " ++ toPyVal mflat + , "" + , "# Load and prepare the common Azimuthal Integrator" + , "ai = load(PONIFILE)" + , "ai.wavelength = WAVELENGTH" + , "ai._empty = numpy.nan" + , "" + , "# Compute the fix part of the mask" + , "mask = numpy.zeros_like(ai.detector.mask, dtype=bool)" + , "mask[0:50, :] = True" + , "mask[910:960, :] = True" + , "mask[:,0:50] = True" + , "mask[:,510:560] = True" + , "if flat is None:" + , " mask = numpy.logical_or(mask, ai.detector.mask)" + , "" + , dummiesForPy mt + , "" + , "# Compute the size of the output" + , "FS = [File(n, mode='r') for n in NEXUSFILES]" + , "NX = 0" + , "NY = 0" + , "for f in FS:" + , " NX = f[MESHX].shape[1]" + , " NY += f[MESHY].shape[0]" + , "" + , "def gen(fs):" + , " for f in fs:" + , " for i in f[IMAGEPATH]:" + , " yield i" + , "" + , "# Create and fill the ouput file" + , "with File(OUTPUT, mode='w') as o:" + , " dataset = o.create_dataset('map', shape=(NY, NX, N), dtype='float')" + , " lines = gen(FS)" + , " for j, line in enumerate(lines):" + , " for i, img in enumerate(line):" + , " tth, I, sigma = ai.integrate1d(img, N, unit=\"2th_deg\"," + , " error_model=\"poisson\", correctSolidAngle=False," + , " method=\"" ++ show m ++ "\"," + , " mask=mask," + , " dummy=DUMMY, delta_dummy=DELTA_DUMMY," + , " safe=False, flat=flat)" + , " dataset[j, i] = I" + ] + +xrdMeshPy' ∷ XrdMeshParams a + → XrdMeshSource -- data source + → FilePath -- ponipath + → DIM1 -- bins + → (Maybe Threshold) -- threshold + → WaveLength -- wavelength + → FilePath -- output h5 + → FilePath -- script name + → Script Py2 +xrdMeshPy' (XrdMeshParams _ mflat m) (XrdMeshSourceNxs (Nxs f h5path)) p b mt w o scriptPath = + xrdMeshPy'' mflat m [f] i x y p b mt w o scriptPath + where + (XrdMeshH5Path (DataItemH5 i _) (DataItemH5 x _) (DataItemH5 y _) _ _ _) = h5path +xrdMeshPy' (XrdMeshParams _ mflat m) (XrdMeshSourceNxsFly nxss) p b mt w o scriptPath = + xrdMeshPy'' mflat m fs i x y p b mt w o scriptPath + where + fs ∷ [FilePath] + fs = [f | (Nxs f _) ← nxss] + + (Nxs _ h5path):_ = nxss + + (XrdMeshFlyH5Path (DataItemH5 i _) (DataItemH5 x _) (DataItemH5 y _) _ _ _) = h5path + +integrateMesh'' ∷ XrdMeshParams a → AbsDirPath → XrdMesh' → IO () +integrateMesh'' p' output (XrdMesh b _ mt s) = do + -- get the poniext for all the scan + (w, PoniExt p _) <- getWaveLengthAndPoniExt p' s + + -- save this poni at the right place + let (ponipath, h5, py) = getOutputPath output s + ponipath `hasContent` poniToText p + + -- create the python script to do the integration + let script = xrdMeshPy' p' s ponipath b mt w h5 py + ExitSuccess ← run script False + + return () diff --git a/contrib/haskell/src/Hkl/Xrd/OneD.hs b/contrib/haskell/src/Hkl/Xrd/OneD.hs new file mode 100644 index 0000000..e3a2ae5 --- /dev/null +++ b/contrib/haskell/src/Hkl/Xrd/OneD.hs @@ -0,0 +1,667 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Xrd.OneD + ( XrdOneD + , XRDRef(..) + , XrdRefSource(..) + , XRDSample(..) + , Threshold(..) + , XrdNxs(..) + , XrdOneDParams(..) + , XrdSource(..) + , PoniExt(..) + -- reference + , getPoseEdf + , getPoniExtRef + -- integration + , integrate + , substract + -- integrateMulti + , integrateMulti + , substractMulti + -- tools + , dummiesForPy + ) where + +import Control.Applicative ((<$>), (<*>), pure) +import Control.Concurrent.Async (mapConcurrently) +import Control.Monad (forM_, forever, void, when, zipWithM_) +import Control.Monad.Morph (hoist) +import Control.Monad.Trans.Maybe (MaybeT, runMaybeT) +import Control.Monad.Trans.State.Strict (StateT, get, put) +import Data.Array.Repa (DIM1, ix1, size) +import Data.Attoparsec.Text (parseOnly) +import qualified Data.List as List (lookup) +import Data.Maybe (fromJust, fromMaybe, isJust) +import Data.Text (Text) +import qualified Data.Text as Text (unlines, pack, intercalate) +import Data.Text.IO (readFile) +import Data.Vector.Storable (concat, head) +import Numeric.LinearAlgebra (fromList) +import Numeric.Units.Dimensional.Prelude (meter, nano, (/~), (*~)) +import System.Exit ( ExitCode( ExitSuccess ) ) +import System.FilePath ((<.>), (</>), dropExtension, replaceExtension, takeFileName, takeDirectory) +import Text.Printf ( printf ) + +import Pipes + ( Consumer + , Pipe + , lift + , (>->) + , runEffect + , await + , yield + ) +import Pipes.Lift ( evalStateP ) +import Pipes.Prelude ( drain, filter, toListM ) +import Pipes.Safe ( runSafeT ) + +import Hkl.C ( Factory ( K6c ) + , Geometry ( Geometry ) + , geometryDetectorRotationGet + ) +import Hkl.DataSource ( DataItem ( DataItemH5 ) + , DataSource( DataSourceH5 ) + , atIndex' + ) +import Hkl.Detector ( Detector ( ZeroD ) ) +import Hkl.Edf ( Edf ( Edf ) + , edfFromFile + ) +import Hkl.Flat ( Flat ) +import Hkl.H5 ( lenH5Dataspace ) +import Hkl.PyFAI ( AIMethod, Poni + , PoniExt ( PoniExt ) + , PoniPath + , Pose ( Pose ) + , move + , poniP + , poniToText + ) +import Hkl.Python ( PyVal + , toPyVal + ) +import Hkl.MyMatrix ( Basis ( HklB ) + , MyMatrix ( MyMatrix ) + ) +import Hkl.Nxs ( DataFrameH5 ( DataFrameH5 ) + , Nxs ( Nxs ) + , XrdOneD + , DataFrameH5Path ( XrdOneDH5Path ) + , withDataFrameH5 + ) +import Hkl.Script ( Gnuplot, Py2 + , Script ( ScriptGnuplot, Py2Script ) + , run + , scriptSave + ) +import Hkl.Types ( AbsDirPath, SampleName + , Source ( Source ) + ) +import Hkl.Utils ( hasContent ) + +-- | TODO +-- * When we skip the last frame there is problem. + +-- * Let's add a method in order to customize the movement of the poni. + +-- | Types + +data Threshold = Threshold Int deriving (Show) + +instance PyVal Threshold where + toPyVal (Threshold i) = toPyVal i + +data XrdRefSource = XrdRefNxs (Nxs XrdOneD) Int + | XrdRefEdf FilePath PoniPath + deriving (Show) + +data XRDRef = XRDRef SampleName AbsDirPath XrdRefSource + deriving (Show) + +data XrdSource = XrdSourceNxs (Nxs XrdOneD) + | XrdSourceEdf [FilePath] + deriving (Show) + +data XrdNxs + = XrdNxs + DIM1 -- bins + DIM1 -- bins for the multibins + (Maybe Threshold) -- threshold use to remove image Intensity + [Int] -- Index of the frames to skip + XrdSource -- data source + deriving (Show) + +data XRDSample = XRDSample SampleName AbsDirPath [XrdNxs] -- ^ nxss + deriving (Show) + +data XrdOneDParams a = XrdOneDParams PoniExt (Maybe (Flat a)) AIMethod + +data DifTomoFrame sh = + DifTomoFrame { difTomoFrameNxs :: Nxs XrdOneD-- ^ nexus of the current frame + , difTomoFrameIdx :: Int -- ^ index of the current frame + , difTomoFrameEOF :: Bool -- ^ is it the eof of the stream + , difTomoFrameGeometry :: Geometry -- ^ diffractometer geometry + , difTomoFramePoniExt :: PoniExt -- ^ the ref poniext + } deriving (Show) + +class Frame t where + len :: t -> IO (Maybe Int) + row :: t -> Int -> MaybeT IO (DifTomoFrame DIM1) + +instance Frame (DataFrameH5 XrdOneD) where + len (DataFrameH5 _ _ _ (DataSourceH5 _ d) _ _) = lenH5Dataspace d + + row d@(DataFrameH5 nxs' _ g d' w ponigen) idx = do + n <- lift $ len d + let eof = fromJust n - 1 == idx + let mu = 0.0 + let komega = 0.0 + let kappa = 0.0 + let kphi = 0.0 + gamma <- g `atIndex'` (ix1 0) + delta <- d' `atIndex'` (ix1 idx) + wavelength <- w `atIndex'` (ix1 0) + let source = Source (Data.Vector.Storable.head wavelength *~ nano meter) + let positions = Data.Vector.Storable.concat [mu, komega, kappa, kphi, gamma, delta] + -- print positions + let geometry = Geometry K6c source positions Nothing + let detector = ZeroD + m <- lift $ geometryDetectorRotationGet geometry detector + let pose = Pose (MyMatrix HklB m) + poniext <- lift $ ponigen pose idx + return $ DifTomoFrame { difTomoFrameNxs = nxs' + , difTomoFrameIdx = idx + , difTomoFrameEOF = eof + , difTomoFrameGeometry = geometry + , difTomoFramePoniExt = poniext + } + +-- type PipeE e a b m r = EitherT e (Pipe a b m) r + +frames :: (Frame a) => Pipe a (DifTomoFrame DIM1) IO () +frames = do + d <- await + (Just n) <- lift $ len d + forM_ [0..n-1] (\i' -> do + f <- lift $ runMaybeT $ row d i' + when (isJust f) (yield (fromJust f))) + +frames' :: (Frame a) => [Int] -> Pipe a (DifTomoFrame DIM1) IO () +frames' is = do + d <- await + forM_ is (\i' -> do + f <- lift $ runMaybeT $ row d i' + when (isJust f) (yield (fromJust f))) + +skip :: [Int] -> DifTomoFrame sh -> Bool +skip is' (DifTomoFrame _ i _ _ _) = notElem i is' + +-- {-# ANN module "HLint: ignore Use camelCase" #-} + + +-- import Graphics.Rendering.Chart.Easy +-- import Graphics.Rendering.Chart.Backend.Diagrams + +-- plotPonies :: FilePath -> [PoniEntry] -> IO () +-- plotPonies f entries = toFile def f $ do +-- layout_title .= "Ponies" +-- setColors [opaque blue] +-- let values = map extract entries +-- plot (line "am" [values [0,(0.5)..400]]) +-- -- plot (points "am points" (signal [0,7..400])) +-- where +-- extract (PoniEntry _ _ (Length poni1) _ _ _ _ _ _) = poni1 + +-- | Usual methods + +dummiesForPy ∷ Maybe Threshold → String +dummiesForPy mt = unlines [ "# Compute the dummy values for the dynamic mask" + , "DUMMY=" ++ dummy + , "DELTA_DUMMY=" ++ delta_dummy + ] + where + dummy = maybe "None" (\_ → "4294967296") mt -- TODO the default value depends on the number od bits per pixels. + delta_dummy = maybe "None" (\(Threshold t) → show (4294967296 - t)) mt + +getScanDir ∷ AbsDirPath → FilePath → FilePath +getScanDir o f = o </> (dropExtension . takeFileName) f + +pgen :: AbsDirPath -> FilePath -> Int -> FilePath +pgen o f i = o </> scandir </> scandir ++ printf "_%02d.poni" i + where + scandir = (dropExtension . takeFileName) f + +getPoseEdf :: FilePath -> IO Pose +getPoseEdf f = do + edf@(Edf lambda _) <- edfFromFile f + let mnes = map Text.pack ["_mu", "_keta", "_kap", "_kphi", "nu", "del"] + let source = Source lambda + let positions = fromList $ map (extract edf) mnes + let geometry = Geometry K6c source positions Nothing + let detector = ZeroD + m <- geometryDetectorRotationGet geometry detector + return $ Pose (MyMatrix HklB m) + where + extract :: Edf -> Text -> Double + extract (Edf _ ms) key = fromMaybe 0.0 (List.lookup key ms) + +poniFromFile :: FilePath -> IO Poni +poniFromFile filename = do + content <- Data.Text.IO.readFile filename + return $ case parseOnly poniP content of + Left _ -> error $ "Can not parse the " ++ filename ++ " poni file" + Right poni -> poni + +getPoniExtRef :: XRDRef -> IO PoniExt +getPoniExtRef (XRDRef _ output (XrdRefNxs nxs'@(Nxs f _) idx)) = do + poniExtRefs <- runSafeT $ + toListM ( withDataFrameH5 nxs' (gen output f) yield + >-> hoist lift ( frames' [idx])) + return $ difTomoFramePoniExt (Prelude.last poniExtRefs) + where + gen :: FilePath -> FilePath -> Pose -> Int -> IO PoniExt + gen root nxs'' p idx' = PoniExt + <$> poniFromFile (root </> scandir ++ printf "_%02d.poni" idx') + <*> pure p + where + scandir = takeFileName nxs'' +getPoniExtRef (XRDRef _ _ (XrdRefEdf e p)) = PoniExt + <$> poniFromFile p + <*> getPoseEdf e + +integrate ∷ XrdOneDParams a → [XRDSample] → IO () +integrate p ss = void $ mapConcurrently (integrate' p) ss + +integrate' ∷ XrdOneDParams a → XRDSample → IO () +integrate' p (XRDSample _ output nxss) = void $ mapConcurrently (integrate'' p output) nxss + +integrate'' ∷ XrdOneDParams a → AbsDirPath → XrdNxs → IO () +integrate'' p output (XrdNxs b _ mt is (XrdSourceNxs nxs'@(Nxs f _))) = do + print f + runSafeT $ runEffect $ + withDataFrameH5 nxs' (gen p) yield + >-> hoist lift (frames + >-> Pipes.Prelude.filter (skip is) + >-> savePonies (pgen output f) + >-> savePy p b mt + >-> saveGnuplot + >-> drain) + where + gen :: XrdOneDParams a -> Pose -> Int -> IO PoniExt + gen (XrdOneDParams ref' _ _) m _idx = return $ move ref' m + +createPy ∷ XrdOneDParams a → DIM1 → Maybe Threshold → FilePath → DifTomoFrame' sh → (Script Py2, FilePath) +createPy (XrdOneDParams _ mflat m) b mt scriptPath (DifTomoFrame' f poniPath) = (Py2Script (script, scriptPath), output) + where + script = Text.unlines $ + map Text.pack ["#!/bin/env python" + , "" + , "import numpy" + , "from h5py import File" + , "from pyFAI import load" + , "" + , "PONIFILE = " ++ toPyVal poniPath + , "NEXUSFILE = " ++ toPyVal nxs' + , "IMAGEPATH = " ++ toPyVal i' + , "IDX = " ++ toPyVal idx + , "N = " ++ toPyVal (size b) + , "OUTPUT = " ++ toPyVal output + , "WAVELENGTH = " ++ toPyVal (w /~ meter) + , "" + , "# load the flat" + , "flat = " ++ toPyVal mflat + , "" + , dummiesForPy mt + , "" + , "ai = load(PONIFILE)" + , "ai.wavelength = WAVELENGTH" + , "ai._empty = numpy.nan" + , "" + , "with File(NEXUSFILE, mode='r') as f:" + , " img = f[IMAGEPATH][IDX]" + , "" + , " # Compute the mask" + , " mask = numpy.zeros_like(img, dtype=bool)" + , " mask[:,550:] = True" + , " #mask_module[0:50, :] = True" + , " #mask_module[910:960, :] = True" + , " #mask_module[:,0:10] = True" + , " if flat is not None: # this should be removed for pyFAI >= 0.13.1 it is now done by PyFAI" + , " mask = numpy.logical_or(mask, flat == 0.0)" + , "" + , " ai.integrate1d(img, N, filename=OUTPUT, unit=\"2th_deg\", error_model=\"poisson\", correctSolidAngle=False, method=\"" ++ show m ++ "\", mask=mask, flat=flat, dummy=DUMMY, delta_dummy=DELTA_DUMMY)" + ] + (Nxs nxs' (XrdOneDH5Path (DataItemH5 i' _) _ _ _)) = difTomoFrameNxs f + idx = difTomoFrameIdx f + output = poniPath `replaceExtension` "dat" + (Geometry _ (Source w) _ _) = difTomoFrameGeometry f + +-- | Pipes + +data DifTomoFrame' sh = DifTomoFrame' { difTomoFrame'DifTomoFrame :: DifTomoFrame sh + , difTomoFrame'PoniPath :: FilePath + } + +savePonies :: (Int -> FilePath) -> Pipe (DifTomoFrame sh) (DifTomoFrame' sh) IO () +savePonies g = forever $ do + f <- await + let filename = g (difTomoFrameIdx f) + let (PoniExt p _) = difTomoFramePoniExt f + lift $ filename `hasContent` (poniToText p) + yield $ DifTomoFrame' { difTomoFrame'DifTomoFrame = f + , difTomoFrame'PoniPath = filename + } + +data DifTomoFrame'' sh = DifTomoFrame'' { difTomoFrame''DifTomoFrame' :: DifTomoFrame' sh + , difTomoFrame''PySCript :: Script Py2 + , difTomoFrame''DataPath :: FilePath + } + +savePy ∷ XrdOneDParams a → DIM1 → Maybe Threshold → Pipe (DifTomoFrame' sh) (DifTomoFrame'' sh) IO () +savePy p b mt = forever $ do + f@(DifTomoFrame' _difTomoFrame poniPath) <- await + let scriptPath = poniPath `replaceExtension`"py" + let (script, dataPath) = createPy p b mt scriptPath f + ExitSuccess <- lift $ run script True + yield $ DifTomoFrame'' { difTomoFrame''DifTomoFrame' = f + , difTomoFrame''PySCript = script + , difTomoFrame''DataPath = dataPath + } + +data DifTomoFrame''' sh = DifTomoFrame''' { difTomoFrame'''DifTomoFrame'' ∷ DifTomoFrame'' sh + , difTomoFrame'''GnuplotScript ∷ Script Gnuplot + , difTomoFrame'''Curves ∷ [FilePath] + } + +mkGnuplot ∷ [FilePath] → FilePath → Script Gnuplot +mkGnuplot fs o = ScriptGnuplot (content, o) + where + content = Text.unlines $ + ["plot \\"] + ++ [Text.intercalate ",\\\n" [ Text.pack (show f ++ " u 1:2 w l") | f <- fs ]] + ++ ["pause -1"] + +saveGnuplot' :: Pipe (DifTomoFrame'' sh) (DifTomoFrame''' sh) (StateT [FilePath] IO) r +saveGnuplot' = forever $ do + curves <- lift get + f@(DifTomoFrame'' (DifTomoFrame' _ poniPath) _ dataPath) <- await + let curves' = curves ++ [dataPath] + let script = mkGnuplot curves' (takeDirectory poniPath </> "plot.gnuplot") + lift . lift $ scriptSave script + lift $ put $! curves' + yield $ DifTomoFrame''' { difTomoFrame'''DifTomoFrame'' = f + , difTomoFrame'''GnuplotScript = script + , difTomoFrame'''Curves = curves' + } + +saveGnuplot :: Pipe (DifTomoFrame'' sh) (DifTomoFrame''' sh) IO r +saveGnuplot = evalStateP [] saveGnuplot' + +-- substract a sample from another one + +substractPy ∷ [FilePath] → [FilePath] → [FilePath] → FilePath → Script Py2 +substractPy fs1 fs2 os scriptPath = Py2Script (content, scriptPath) + where + content ∷ Text + content = Text.unlines $ + map Text.pack ["#!/bin/env python" + , "" + , "import numpy" + , "" + , "S1 = " ++ toPyVal fs1 + , "S2 = " ++ toPyVal fs2 + , "OUTPUTS = " ++ toPyVal os + , "" + , "def substract(f1, f2, o):" + , " a1 = numpy.genfromtxt(f1)" + , " a2 = numpy.genfromtxt(f2)" + , " res = numpy.copy(a2)" + , " res[:,1] -= a1[:,1]" + , " # TODO deal with the error propagation" + , " numpy.savetxt(output, res)" + , "" + , "for (s1, s2, output) in zip(S1, S2, OUTPUTS):" + , " substract(s1, s2, output)" + ] + +targetP ∷ (Int → FilePath) → Pipe (DifTomoFrame sh) FilePath IO () +targetP g = forever $ do + f ← await + let poniPath = g (difTomoFrameIdx f) + let dataPath = poniPath `replaceExtension` "dat" + yield dataPath + +target' ∷ XrdOneDParams a → AbsDirPath → XrdNxs → IO (FilePath, [FilePath]) +target' p output (XrdNxs _ _ _ is (XrdSourceNxs nxs'@(Nxs f _))) = do + fs ← runSafeT $ toListM $ + withDataFrameH5 nxs' (gen p) yield + >-> hoist lift (frames + >-> Pipes.Prelude.filter (skip is) + >-> targetP (pgen output f) + ) + return (getScanDir output f, fs) + where + gen :: XrdOneDParams a -> Pose -> Int -> IO PoniExt + gen (XrdOneDParams ref' _ _) m _idx = return $ move ref' m + +targets ∷ XrdOneDParams a → XRDSample → IO [(FilePath, [FilePath])] +targets p (XRDSample _ output nxss) = mapConcurrently (target' p output) nxss + +substract' ∷ XrdOneDParams a → XRDSample → XRDSample → IO () +substract' p s1@(XRDSample name _ _) s2 = do + -- compute the output of the s1 sample + -- we take only the first list of the sample + f1s:_ ← targets p s1 + -- compute the output of the s2 sample + f2s ← targets p s2 + -- do the substraction via a python script and add the gnuplot file + _ ← mapConcurrently (go f1s) f2s + return () + where + go ∷ (FilePath, [FilePath]) → (FilePath, [FilePath]) → IO () + go (_, f1) (d, f2) = do + -- compute the substracted output file names take into account + -- that f1 and f2 could have different length + let outputs = [dropExtension f ++ "-" ++ name <.> "dat" | (_, f) ← zip f1 f2] + -- compute the script name + let scriptPath = d </> "substract.py" + let script = substractPy f1 f2 outputs scriptPath + ExitSuccess ← run script False + -- gnuplot + let gnuplotPath = d </> "substract.gnuplot" + scriptSave $ mkGnuplot outputs gnuplotPath + return () + +substract ∷ XrdOneDParams a → XRDSample → [XRDSample] → IO () +substract p s ss = mapM_ (substract' p s) ss + +-- | PyFAI MultiGeometry + +integrateMulti ∷ XrdOneDParams a → [XRDSample] → IO () +integrateMulti p samples = mapM_ (integrateMulti' p) samples + +integrateMulti' ∷ XrdOneDParams a → XRDSample → IO () +integrateMulti' p (XRDSample _ output nxss) = mapM_ (integrateMulti'' p output) nxss + +integrateMulti'' ∷ XrdOneDParams a → AbsDirPath → XrdNxs → IO () +integrateMulti'' p output (XrdNxs _ mb mt is (XrdSourceNxs nxs'@(Nxs f _))) = do + print f + runSafeT $ runEffect $ + withDataFrameH5 nxs' (gen p) yield + >-> hoist lift (frames + >-> Pipes.Prelude.filter (skip is) + >-> savePonies (pgen output f) + >-> saveMultiGeometry p mb mt) + where + gen :: XrdOneDParams a -> Pose -> Int -> IO PoniExt + gen (XrdOneDParams ref' _ _) m _idx = return $ move ref' m + +integrateMulti'' p output (XrdNxs b _ mt _ (XrdSourceEdf fs)) = do + -- generate all the ponies + zipWithM_ (go p) fs ponies + + -- generate the multi.py python script + let scriptPath = output </> "multi.py" + let (script, _) = createMultiPyEdf p b mt fs ponies scriptPath (output </> "multi.dat") + scriptSave script + where + ponies = [output </> (dropExtension . takeFileName) f ++ ".poni" | f <- fs] + + go ∷ XrdOneDParams a → FilePath → FilePath → IO () + go (XrdOneDParams ref _ _) f o = do + m <- getPoseEdf f + let (PoniExt p' _) = move ref m + o `hasContent` (poniToText p') + +createMultiPy ∷ XrdOneDParams a → DIM1 → Maybe Threshold → FilePath → DifTomoFrame' sh → [(Int, FilePath)] → (Script Py2, FilePath) +createMultiPy (XrdOneDParams _ mflat _) b mt scriptPath (DifTomoFrame' f _) idxPonies = (Py2Script (content, scriptPath), output) + where + content = Text.unlines $ + map Text.pack ["#!/bin/env python" + , "" + , "import numpy" + , "from h5py import File" + , "from pyFAI.multi_geometry import MultiGeometry" + , "" + , "NEXUSFILE = " ++ toPyVal nxs' + , "IMAGEPATH = " ++ toPyVal i' + , "BINS = " ++ toPyVal (size b) + , "OUTPUT = " ++ toPyVal output + , "WAVELENGTH = " ++ toPyVal (w /~ meter) + , "THRESHOLD = " ++ toPyVal mt + , "" + , "# load the flat" + , "flat = " ++ toPyVal mflat + , "" + , "# Load all images" + , "PONIES = " ++ toPyVal ponies + , "IDXS = " ++ toPyVal idxs + , "" + , "# Read all the images" + , "imgs = []" + , "with File(NEXUSFILE, mode='r') as f:" + , " for idx in IDXS:" + , " imgs.append(f[IMAGEPATH][idx])" + , "" + , "# Compute the mask" + , "mask = numpy.zeros_like(imgs[0], dtype=bool)" + , "mask[:,550:] = True" + , "if flat is not None: # this should be removed for pyFAI >= 0.13.1 it is now done by PyFAI" + , " mask = numpy.logical_or(mask, flat == 0.0)" + , "lst_mask = []" + , "for img in imgs: # remove all pixels above the threshold" + , " if THRESHOLD is not None:" + , " mask_t = numpy.where(img > THRESHOLD, True, False)" + , " lst_mask.append(numpy.logical_or(mask, mask_t))" + , " else:" + , " lst_mask.append(mask)" + , "" + , "# Integration multi-geometry 1D" + , "mg = MultiGeometry(PONIES, unit=\"2th_deg\", radial_range=(0,80))" + , "p = mg.integrate1d(imgs, BINS, lst_mask=lst_mask, lst_flat=flat)" + , "" + , "# Save the datas" + , "numpy.savetxt(OUTPUT, numpy.array(p).T)" + ] + (Nxs nxs' (XrdOneDH5Path (DataItemH5 i' _) _ _ _)) = difTomoFrameNxs f + output = "multi.dat" + (Geometry _ (Source w) _ _) = difTomoFrameGeometry f + (idxs, ponies) = unzip idxPonies + +createMultiPyEdf ∷ XrdOneDParams a → DIM1 → Maybe Threshold → [FilePath] → [FilePath] → FilePath → FilePath → (Script Py2, FilePath) +createMultiPyEdf (XrdOneDParams _ mflat _) b mt edfs ponies scriptPath output = (Py2Script (content, scriptPath), output) + where + content = Text.unlines $ + map Text.pack ["#!/bin/env python" + , "" + , "import numpy" + , "from fabio import open" + , "from pyFAI.multi_geometry import MultiGeometry" + , "" + , "EDFS = " ++ toPyVal edfs + , "PONIES = " ++ toPyVal ponies + , "BINS = " ++ toPyVal (size b) + , "OUTPUT = " ++ toPyVal output + , "THRESHOLD = " ++ toPyVal mt + , "" + , "# load the flat" + , "flat = " ++ toPyVal mflat + , "" + , "# Read all the images" + , "imgs = [open(edf).data for edf in EDFS]" + , "" + , "# Compute the mask" + , "mask = numpy.zeros_like(imgs[0], dtype=bool)" + , "if THRESHOLD is not None:" + , " for img in imgs:" + , " mask_t = numpy.where(img > THRESHOLD, True, False)" + , " mask = numpy.logical_or(mask, mask_t)" + , "" + , "# Integration multi-geometry 1D" + , "mg = MultiGeometry(PONIES, unit=\"2th_deg\", radial_range=(0,80))" + , "p = mg.integrate1d(imgs, BINS, lst_mask=mask)" + , "" + , "# Save the datas" + , "numpy.savetxt(OUTPUT, numpy.array(p).T)" + ] + +saveMulti' ∷ XrdOneDParams a → DIM1 → Maybe Threshold → Consumer (DifTomoFrame' sh) (StateT [(Int, FilePath)] IO) r +saveMulti' p b mt = forever $ do + idxPonies <- lift get + f'@(DifTomoFrame' f@(DifTomoFrame _ idx _ _ _) poniPath) <- await + let directory = takeDirectory poniPath + let filename = directory </> "multi.py" + let (script, _) = createMultiPy p b mt filename f' idxPonies + ExitSuccess ← lift . lift $ if (difTomoFrameEOF f) then (run script True) else return ExitSuccess + lift $ put $! (idxPonies ++ [(idx, poniPath)]) + +saveMultiGeometry ∷ XrdOneDParams a → DIM1 → Maybe Threshold → Consumer (DifTomoFrame' sh) IO r +saveMultiGeometry p b mt = evalStateP [] (saveMulti' p b mt) + + +-- substract a sample from another one + +targetMulti' ∷ XrdOneDParams a → AbsDirPath → XrdNxs → (FilePath, FilePath) +targetMulti' _ output (XrdNxs _ _ _ _ (XrdSourceNxs (Nxs f _))) = (d, o) + where + d = getScanDir output f + o = d </> "multi.dat" + +targetMulti ∷ XrdOneDParams a → XRDSample → [(FilePath, FilePath)] +targetMulti p (XRDSample _ output nxss) = map (targetMulti' p output) nxss + +substractMulti' ∷ XrdOneDParams a → XRDSample → XRDSample → IO () +substractMulti' p s1@(XRDSample name _ _) s2 = do + -- compute the output of the s1 sample + -- we take only the first list of the sample + let f1s:_ = targetMulti p s1 + -- compute the output of the s2 sample + let f2s = targetMulti p s2 + -- do the substraction via a python script and add the gnuplot file + _ ← mapConcurrently (go f1s) f2s + + return () + where + go ∷ (FilePath, FilePath) → (FilePath, FilePath) → IO () + go (_, f1) (d, f2) = do + -- compute the substracted output file names + let outputs = dropExtension f2 ++ "-" ++ name <.> "dat" + -- compute the script name + let scriptPath = d </> "multi-substract.py" + let script = substractPy [f1] [f2] [outputs] scriptPath + ExitSuccess ← run script False + -- gnuplot + let gnuplotPath = d </> "multi-substract.gnuplot" + scriptSave $ mkGnuplot [outputs] gnuplotPath + return () + +substractMulti ∷ XrdOneDParams a → XRDSample → [XRDSample] → IO () +substractMulti p s ss = mapM_ (substractMulti' p s) ss diff --git a/contrib/haskell/src/Hkl/Xrd/ZeroD.hs b/contrib/haskell/src/Hkl/Xrd/ZeroD.hs new file mode 100644 index 0000000..5de2a42 --- /dev/null +++ b/contrib/haskell/src/Hkl/Xrd/ZeroD.hs @@ -0,0 +1,118 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE UnicodeSyntax #-} + +module Hkl.Xrd.ZeroD + ( XrdZeroDCalibration(..) + , XrdZeroDSample(..) + , XrdZeroDSource(..) + , XrdZeroDParams(..) + ) where + +import Data.List (intercalate) +import Data.Text (unlines, pack) +import Numeric.Units.Dimensional.Prelude (meter, nano, (*~)) +import System.Exit ( ExitCode( ExitSuccess ) ) +import System.FilePath.Posix ((</>), takeFileName) +import Text.Printf ( printf ) + +import Hkl.DataSource ( DataItem ( DataItemH5, DataItemConst ) ) +import Hkl.Detector ( Detector ) +import Hkl.Edf ( ExtractEdf, extractEdf ) +import Hkl.Flat ( Flat ) +import Hkl.PyFAI ( AIMethod,Calibrant, PoniExt, Pose + , toPyFAICalibArg ) +import Hkl.Python ( toPyVal ) +import Hkl.Nxs ( DataFrameH5Path( XrdZeroDH5Path ) + , Nxs ( Nxs ) + , XrdZeroD + ) +import Hkl.Script ( Script ( Py2Script, ScriptSh ) + , Py2, Sh + , run + , scriptSave + ) +import Hkl.Types ( AbsDirPath, SampleName, WaveLength ) + +-- | Types + +data XrdZeroDSource = XrdZeroDSourceNxs (Nxs XrdZeroD) deriving (Show) + +data XrdZeroDSample = XrdZeroDSample SampleName AbsDirPath [XrdZeroDSource] deriving (Show) + +data XrdZeroDCalibration a = XrdZeroDCalibration XrdZeroDSample (Detector a) Calibrant deriving (Show) + +data XrdZeroDParams a = XrdZeroDParams PoniExt (Maybe (Flat a)) AIMethod deriving (Show) + +data XrdZeroDFrame = XrdMeshFrame WaveLength Pose deriving (Show) + +edf ∷ AbsDirPath → FilePath → Int → FilePath +edf o n i = o </> f + where + f = (takeFileName n) ++ printf "_%02d.edf" i + +scriptExtractEdf ∷ AbsDirPath → [XrdZeroDSource] → Script Py2 +scriptExtractEdf o es = Py2Script (content, scriptPath) + where + content = Data.Text.unlines $ + map Data.Text.pack [ "#!/usr/bin/env python" + , "" + , "from fabio.edfimage import edfimage" + , "from h5py import File" + , "" + , "NEXUSFILES = " ++ toPyVal nxss + , "IDXS = " ++ toPyVal idxs + , "IMAGEPATHS = " ++ toPyVal (imgs ∷ [String]) + , "OUTPUTS = " ++ toPyVal outputs + , "" + , "for filename, i, p, o in zip(NEXUSFILES, IDXS, IMAGEPATHS, OUTPUTS):" + , " with File(filename, mode='r') as f:" + , " edfimage(f[p][i]).write(o)" + ] + + idx ∷ Int + idx = 0 + + (nxss, idxs, imgs) = unzip3 [(f, idx, img) + | (XrdZeroDSourceNxs (Nxs f (XrdZeroDH5Path (DataItemH5 img _) _))) ← es] + + outputs ∷ [FilePath] + outputs = zipWith (edf o) nxss idxs + + scriptPath ∷ FilePath + scriptPath = o </> "pre-calibration.py" + +scriptPyFAICalib ∷ AbsDirPath → XrdZeroDSource → Detector a → Calibrant → Script Sh +scriptPyFAICalib o e@(XrdZeroDSourceNxs (Nxs n _)) d c = ScriptSh (content, scriptPath) + where + content = Data.Text.unlines $ + map Data.Text.pack [ "#!/usr/bin/env sh" + , "" + , "pyFAI-calib " ++ intercalate " " args + ] + + args = [ toPyFAICalibArg (readWavelength e) + , toPyFAICalibArg c + , toPyFAICalibArg d + , toPyFAICalibArg (edf o n i) ] + + scriptPath ∷ FilePath + scriptPath = o </> (takeFileName n) ++ printf "_%02d.sh" i + + i ∷ Int + i = 0 + +readWavelength :: XrdZeroDSource -> WaveLength +readWavelength (XrdZeroDSourceNxs (Nxs _ (XrdZeroDH5Path _ (DataItemConst w)))) = w *~ nano meter + +instance ExtractEdf (XrdZeroDCalibration a) where + extractEdf (XrdZeroDCalibration s d c) = do + let script = scriptExtractEdf o es + ExitSuccess ← run script False + mapM_ go es + return () + where + go e = scriptSave $ scriptPyFAICalib o e d c + + (XrdZeroDSample _ o es) = s diff --git a/contrib/haskell/src/Tango/DeviceProxy.hsc b/contrib/haskell/src/Tango/DeviceProxy.hsc new file mode 100644 index 0000000..923ccf0 --- /dev/null +++ b/contrib/haskell/src/Tango/DeviceProxy.hsc @@ -0,0 +1,47 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE EmptyDataDecls #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE RecordWildCards #-} + +module Tango.DeviceProxy ( + deviceproxy + , DeviceProxy ) where + +import Control.Exception + +import Foreign.C +import Foreign.Ptr +import Foreign.Storable +import Foreign.Marshal.Array +import Foreign.Marshal.Alloc + +#let alignment t = "%lu", (unsigned long)offsetof(struct {char x__; t (y__); }, y__) +#include "tango.h" + +data DeviceProxy = DeviceProxy + +foreign import ccall "_ZN5Tango11DeviceProxyC1EPKcPN5CORBA3ORBE" deviceproxy_DeviceProxy :: (Ptr DeviceProxy) -> CString -> Ptr a -> IO () + +class New a where + new :: IO (Ptr a) + +instance Storable DeviceProxy where + sizeOf _ = #{size Tango::DeviceProxy} + +deviceproxy :: String -> IO (Ptr DeviceProxy) +deviceproxy d = do + device <- newCString d + dev <- malloc :: IO (Ptr DeviceProxy) + deviceproxy_DeviceProxy dev device nullPtr + return dev + +main :: IO () +main = do + diffractometer <- catch (deviceproxy "toto") + (\e -> do let err = show (e :: IOException) + hPutStr stderr ("Warning: Couldn't open " ++ f ++ ": " ++ err) + return "") + return () diff --git a/contrib/haskell/src/ghkl.hs b/contrib/haskell/src/ghkl.hs new file mode 100644 index 0000000..0ae7f19 --- /dev/null +++ b/contrib/haskell/src/ghkl.hs @@ -0,0 +1,98 @@ +{-# LANGUAGE UnicodeSyntax #-} + +module Main where + +import Data.Vector.Storable (Vector, fromList) +import Numeric.Units.Dimensional.Prelude (nano, meter, (*~)) +import Pipes +import qualified Pipes.Prelude as P + +import Hkl + +testSirius :: IO () +testSirius = runEffect $ fromToPipe 6 from to + >-> enginesTrajectoryPipe engine + >-> solveTrajPipe geometry detector gaAs + >-> P.tee P.print + >-> P.drain + -- >-> computePipe detector gaAs + -- >-> P.print + where + gaAs :: Sample Cubic + gaAs = Sample "GaAs" (Cubic (0.56533 *~ nano meter)) + (Parameter "ux" (-90.003382) (Range (-180) 180)) + (Parameter "uy" 0.12907 (Range (-180) 180)) + (Parameter "uz" (-159.91372) (Range (-180) 180)) + + geometry ∷ Geometry + geometry = Geometry SoleilSiriusKappa (Source (0.1458637 *~ nano meter)) + (fromList [-0.5193202, 64.7853160, 133.5621380, -80.9690000, -0.0223369, 30.0000299]) + (Just [ Parameter "mu" (-0.5193202) (Range (-180) 180) + , Parameter "komega" 64.7853160 (Range (-180) 180) + , Parameter "kappa" 133.5621380 (Range (-180) 180) + , Parameter "kphi" (-80.9690000) (Range (-180) 180) + , Parameter "delta" (-0.0223369) (Range (-180) 180) + , Parameter "gamma" 30.0000299 (Range (-180) 180)]) + + detector ∷ Detector ZeroD + detector = ZeroD + + engine ∷ Engine + engine = Engine "hkl" [ Parameter "h" 0.0 (Range (-1.0) 1.0) + , Parameter "k" 0.0 (Range (-1.0) 1.0) + , Parameter "l" 2.0 (Range (-1.0) 1.0) + ] + (Mode "bissector_vertical" []) + + from ∷ Vector Double + from = fromList [0, 0, 1] + + to ∷ Vector Double + to = fromList [0, 0, 6] + +test :: IO () +test = do + let sample = Sample "test" (Orthorhombic + (1.05394 *~ nano meter) + (0.25560 *~ nano meter) + (1.49050 *~ nano meter)) + (Parameter "ux" (-89.8821) (Range (-180) 180)) + (Parameter "uy" 0.1733 (Range (-180) 180)) + (Parameter "uz" (-84.0081) (Range (-180) 180)) + + let geometry = Geometry Uhv (Source (0.0672929 *~ nano meter)) + (fromList [0.1794, -160.0013, 21.1381, 0.5194]) + (Just [ Parameter "mu" 0.1794 (Range (-180) 180) + , Parameter "omega" (-160.0013) (Range (-180) 180) + , Parameter "delta" 21.1381 (Range (-180) 180) + , Parameter "gamma" 0.5194 (Range (-180) 180)]) + let detector = ZeroD + + -- compute the pseudo axes values + pseudoAxes <- compute geometry detector sample + print pseudoAxes + + -- solve a pseudo axis problem for the given engine + let engine = Engine "hkl" [ Parameter "h" 4.0 (Range (-1.0) 1.0) + , Parameter "k" 1.0 (Range (-1.0) 1.0) + , Parameter "l" 0.3 (Range (-1.0) 1.0) + ] + (Mode "zaxis" []) + + print =<< solve geometry detector sample engine + + -- let from = fromList [0, 0, 1 :: Double] + -- let to = fromList [0, 1, 1 :: Double] + -- runEffect $ fromToPipe 20 from to + -- >-> P.print + -- -- solve a trajectory with Pipes + -- runEffect $ fromToPipe 10000 from to + -- >-> enginesTrajectoryPipe engine + -- >-> solveTrajPipe factory geometry detector sample + -- >-> P.print + -- -- >-> P.drain + + return () + +main :: IO () +main = testSirius diff --git a/contrib/haskell/src/hkl.hs b/contrib/haskell/src/hkl.hs new file mode 100644 index 0000000..4e46a8f --- /dev/null +++ b/contrib/haskell/src/hkl.hs @@ -0,0 +1,73 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} +{- + Copyright : Copyright (C) 2014-2018 Synchrotron Soleil + License : GPL3+ + + Maintainer : picca@synchrotron-soleil.fr + Stability : Experimental + Portability: GHC only? +-} + +import Numeric.LinearAlgebra (Vector, Matrix, + vecdisp, disps, + dispf) + +import Numeric.Units.Dimensional.Prelude (nano, meter, degree, + (*~), + (*~~), (/~~)) + +import Options.Applicative hiding ((<>)) + +import Hkl.Lattice +import Hkl.Diffractometer + +dispv :: Vector Double -> IO () +dispv = putStr . vecdisp (disps 2) + +disp :: Matrix Double -> IO () +disp = putStr . dispf 3 + +-- command parsing +data Command + = Ca Double Double Double -- ca command + +data Options + = Options Command + +withInfo :: Parser a -> String -> ParserInfo a +withInfo opts desc = info (helper <*> opts) $ progDesc desc + +parseCa :: Parser Command +parseCa = Ca + <$> argument auto (metavar "H") + <*> argument auto (metavar "K") + <*> argument auto (metavar "L") + +parseCommand :: Parser Command +parseCommand = subparser $ + command "ca" (parseCa `withInfo` "compute angles for the given hkl") + +parseOptions :: Parser Options +parseOptions = Options <$> parseCommand + +-- Actual program logic +run :: Options -> IO () +run (Options cmd) = + case cmd of + Ca h k l-> do + print (solution /~~ degree) + dispv (computeHkl e4c solution lattice) + disp path + where + (sol, path) = computeAngles e4c angles lattice mode [h, k, l] + s = [30.0, 0.0, 0.0, 0.0, 10.0, 0.0] + d = [60.0] + angles = (s ++ d) *~~ degree + solution = fromMode mode sol angles + lattice = Cubic (1.54 *~ nano meter) + mode = ModeHklE4CConstantPhi + +main :: IO () +main = run =<< execParser + (parseOptions `withInfo` "Interact with hkl API") diff --git a/contrib/haskell/src/hkl3d.hs b/contrib/haskell/src/hkl3d.hs new file mode 100644 index 0000000..751a617 --- /dev/null +++ b/contrib/haskell/src/hkl3d.hs @@ -0,0 +1,8 @@ +import Hkl.Projects + +{-# ANN module "HLint: ignore Use camelCase" #-} + +main :: IO () +-- main = main_calibration +-- main = main_diffabs +main = main_sixs diff --git a/contrib/haskell/src/xrd.hs b/contrib/haskell/src/xrd.hs new file mode 100644 index 0000000..ced23c6 --- /dev/null +++ b/contrib/haskell/src/xrd.hs @@ -0,0 +1,16 @@ +module Main where + +import Hkl.Projects + +main :: IO () +main = do + -- irdrx + -- martinetto' + -- melle + -- d2am + charlier + -- laure + -- hercules + -- hamon + -- schlegel + -- romeden |