binoculars-ng/src/Hkl/Binoculars/Projections/Hkl.hs

   1 {-# LANGUAGE DataKinds             #-}
   2 {-# LANGUAGE DeriveAnyClass        #-}
   3 {-# LANGUAGE DeriveGeneric         #-}
   4 {-# LANGUAGE FlexibleContexts      #-}
   5 {-# LANGUAGE FlexibleInstances     #-}
   6 {-# LANGUAGE MultiParamTypeClasses #-}
   7 {-# LANGUAGE OverloadedStrings     #-}
   8 {-# LANGUAGE ScopedTypeVariables   #-}
   9 {-# LANGUAGE StandaloneDeriving    #-}
  10 {-# LANGUAGE TypeFamilies          #-}
  11
  12 {-# OPTIONS_GHC -fno-warn-orphans #-}
  13
  14 {-
  15     Copyright  : Copyright (C) 2014-2024 Synchrotron SOLEIL
  16                                          L'Orme des Merisiers Saint-Aubin
  17                                          BP 48 91192 GIF-sur-YVETTE CEDEX
  18     License    : GPL3+
  19
  20     Maintainer : Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
  21     Stability  : Experimental
  22     Portability: GHC only (not tested)
  23 -}
  24
  25 module Hkl.Binoculars.Projections.Hkl
  26     ( Config(..)
  27     , DataFrameHkl'(..)
  28     , newHkl
  29     , processHkl
  30     , updateHkl
  31     ) where
  32
  33
  34 import           Control.Concurrent.Async           (mapConcurrently)
  35 import           Control.Monad.Catch                (MonadThrow)
  36 import           Control.Monad.IO.Class             (MonadIO (liftIO))
  37 import           Control.Monad.Logger               (MonadLogger, logDebugN,
  38                                                      logInfoN)
  39 import           Control.Monad.Reader               (MonadReader, ask, forM_,
  40                                                      forever)
  41 import           Data.Aeson                         (FromJSON, ToJSON,
  42                                                      eitherDecode', encode)
  43 import           Data.ByteString.Lazy               (fromStrict, toStrict)
  44 import           Data.Functor.Identity              (Identity)
  45 import           Data.HashMap.Strict                (fromList)
  46 import           Data.Ini                           (Ini (..))
  47 import           Data.Ini.Config.Bidir              (FieldValue (..))
  48 import           Data.Text                          (pack, unpack)
  49 import           Data.Text.Encoding                 (decodeUtf8, encodeUtf8)
  50 import           Data.Text.IO                       (putStr)
  51 import           Data.Vector.Storable.Mutable       (unsafeWith)
  52 import           Foreign.C.Types                    (CDouble (..))
  53 import           Foreign.ForeignPtr                 (withForeignPtr)
  54 import           GHC.Generics                       (Generic)
  55 import           Path                               (Abs, Dir, Path)
  56 import           Pipes                              (await, each, runEffect,
  57                                                      (>->))
  58 import           Pipes.Prelude                      (filter, map, tee, toListM)
  59 import           Pipes.Safe                         (runSafeP, runSafeT)
  60 import           Text.Printf                        (printf)
  61
  62 import           Hkl.Binoculars.Common
  63 import           Hkl.Binoculars.Config
  64 import           Hkl.Binoculars.Config.Common
  65 import           Hkl.Binoculars.Config.Sample
  66 import           Hkl.Binoculars.Pipes
  67 import           Hkl.Binoculars.Projections
  68 import           Hkl.Binoculars.Projections.QCustom
  69 import           Hkl.C.Binoculars
  70 import           Hkl.DataSource
  71 import           Hkl.Detector
  72 import           Hkl.H5
  73 import           Hkl.HKD
  74 import           Hkl.Image
  75 import           Hkl.Orphan                         ()
  76 import           Hkl.Pipes
  77 import           Hkl.Repa
  78 import           Hkl.Sample
  79 import           Hkl.Utils
  80
  81 ----------------
  82 -- DataPath's --
  83 ----------------
  84
  85 data DataFrameHkl' f
  86   = DataFrameHkl
  87     { dataFrameHkl'DataFrameQCustom :: HKD f DataFrameQCustom
  88     , dataFrameHkl'Sample           :: HKD f Sample
  89     }
  90   deriving (Generic)
  91
  92 deriving instance Show (DataFrameHkl' DataSourcePath)
  93 instance FromJSON (DataFrameHkl' DataSourcePath)
  94 instance ToJSON (DataFrameHkl' DataSourcePath)
  95
  96 defaultDataSourcePath'DataFrameHkl :: DataFrameHkl' DataSourcePath
  97 defaultDataSourcePath'DataFrameHkl = DataFrameHkl
  98                                      default'DataSourcePath'DataFrameQCustom
  99                                      default'DataSourcePath'Sample
 100
 101 instance HasFieldComment (DataFrameHkl' DataSourcePath) where
 102   fieldComment _ = [ "`datapath` internal value used to find the data in the data file."
 103                    , ""
 104                    , "This value is for expert only."
 105                    , ""
 106                    , "default value: <not set>"
 107                    ]
 108
 109 instance HasFieldValue (DataFrameHkl' DataSourcePath) where
 110   fieldvalue = FieldValue
 111                { fvParse = eitherDecode' . fromStrict . encodeUtf8
 112                , fvEmit = decodeUtf8 . toStrict . encode
 113                }
 114
 115 ------------
 116 -- Config --
 117 ------------
 118
 119 data instance Config 'HklProjection
 120   = BinocularsConfig'Hkl
 121     { binocularsConfig'Hkl'Common                 :: BinocularsConfig'Common
 122     , binocularsConfig'Hkl'Sample                 :: BinocularsConfig'Sample
 123     , binocularsConfig'Hkl'ProjectionType         :: ProjectionType
 124     , binocularsConfig'Hkl'ProjectionResolution   :: Resolutions DIM3
 125     , binocularsConfig'Hkl'ProjectionLimits       :: Maybe (RLimits DIM3)
 126     , binocularsConfig'Hkl'DataPath               :: DataFrameHkl' DataSourcePath
 127     } deriving (Generic)
 128
 129 newtype instance Args 'HklProjection = Args'HklProjection (Maybe ConfigRange)
 130
 131 default'BinocularsConfig'Hkl :: Config 'HklProjection
 132 default'BinocularsConfig'Hkl
 133   = BinocularsConfig'Hkl
 134     { binocularsConfig'Hkl'Common = default'BinocularsConfig'Common
 135     , binocularsConfig'Hkl'Sample = default'BinocularsConfig'Sample
 136     , binocularsConfig'Hkl'ProjectionType = HklProjection
 137     , binocularsConfig'Hkl'ProjectionResolution = Resolutions3 0.01 0.01 0.01
 138     , binocularsConfig'Hkl'ProjectionLimits  = Nothing
 139     , binocularsConfig'Hkl'DataPath = defaultDataSourcePath'DataFrameHkl
 140     }
 141
 142 overload'DataSourcePath'DataFrameHkl :: BinocularsConfig'Common
 143                                      -> BinocularsConfig'Sample
 144                                      -> DataFrameHkl' DataSourcePath
 145                                      -> DataFrameHkl' DataSourcePath
 146 overload'DataSourcePath'DataFrameHkl common sample (DataFrameHkl qCustomPath samplePath)
 147   = DataFrameHkl newQCustomPath newSamplePath
 148   where
 149     newQCustomPath = overload'DataSourcePath'DataFrameQCustom common Nothing qCustomPath
 150     newSamplePath = overload'DataSourcePath'Sample sample samplePath
 151
 152 instance HasIniConfig 'HklProjection where
 153
 154   getConfig content@(ConfigContent cfg) (Args'HklProjection mr) capabilities = do
 155
 156     common <- parse'BinocularsConfig'Common cfg mr capabilities
 157     sample <- parse'BinocularsConfig'Sample cfg
 158
 159     BinocularsConfig'Hkl
 160       <$> pure common
 161       <*> pure sample
 162       <*> parseFDef cfg "projection" "type" (binocularsConfig'Hkl'ProjectionType default'BinocularsConfig'Hkl)
 163       <*> parseFDef cfg "projection" "resolution" (binocularsConfig'Hkl'ProjectionResolution default'BinocularsConfig'Hkl)
 164       <*> parseMb cfg "projection" "limits"
 165       <*> (pure $ eitherF (const $ guess'DataSourcePath'DataFrameHkl common sample content) (parse' cfg "input" "datapath")
 166            (\md -> case md of
 167                     Nothing -> guess'DataSourcePath'DataFrameHkl common sample content
 168                     Just d  ->  overload'DataSourcePath'DataFrameHkl common sample d))
 169
 170
 171 instance ToIni (Config 'HklProjection) where
 172   toIni c = toIni (binocularsConfig'Hkl'Common c)
 173             `mergeIni`
 174             toIni (binocularsConfig'Hkl'Sample c)
 175             `mergeIni`
 176             Ini { iniSections = fromList [ ("input", elemFDef' "datapath" binocularsConfig'Hkl'DataPath c default'BinocularsConfig'Hkl)
 177                                          , ("projection", elemFDef' "type" binocularsConfig'Hkl'ProjectionType c default'BinocularsConfig'Hkl
 178                                                           <> elemFDef' "resolution" binocularsConfig'Hkl'ProjectionResolution c default'BinocularsConfig'Hkl
 179                                                           <> elemFMbDef' "limits" binocularsConfig'Hkl'ProjectionLimits c default'BinocularsConfig'Hkl
 180                                            )
 181                                          ]
 182                 , iniGlobals = []
 183                 }
 184
 185
 186 ----------------
 187 -- Projection --
 188 ----------------
 189
 190 {-# INLINE spaceHkl #-}
 191 spaceHkl :: Detector b DIM2 -> Array F DIM3 Double -> Resolutions DIM3 -> Maybe Mask -> Maybe (RLimits DIM3) -> Bool -> Space DIM3 -> DataFrameHkl' Identity -> IO (DataFrameSpace DIM3)
 192 spaceHkl det pixels rs mmask' mlimits doPolarizationCorrection space@(Space fSpace) (DataFrameHkl (DataFrameQCustom att g img _) samplePath) = do
 193   withNPixels det $ \nPixels ->
 194     withForeignPtr g $ \geometry ->
 195     withSample samplePath $ \sample ->
 196     withForeignPtr (toForeignPtr pixels) $ \pix ->
 197     withResolutions rs $ \nr r ->
 198     withMaybeMask mmask' $ \ mask'' ->
 199     withPixelsDims pixels $ \ndim dims ->
 200     withMaybeLimits mlimits rs $ \nlimits limits ->
 201     withForeignPtr fSpace $ \pSpace -> do
 202     case img of
 203       (ImageInt32 arr) -> unsafeWith arr $ \i -> do
 204         {-# SCC "hkl_binoculars_space_hkl_int32_t" #-} c'hkl_binoculars_space_hkl_int32_t pSpace geometry sample i nPixels (CDouble . unAttenuation $ att) pix (toEnum ndim) dims r (toEnum nr) mask'' limits (toEnum nlimits) (toEnum . fromEnum $ doPolarizationCorrection)
 205       (ImageWord16 arr) -> unsafeWith arr $ \i -> do
 206         {-# SCC "hkl_binoculars_space_hkl_uint16_t" #-} c'hkl_binoculars_space_hkl_uint16_t pSpace geometry sample i nPixels (CDouble . unAttenuation $ att) pix (toEnum ndim) dims r (toEnum nr) mask'' limits (toEnum nlimits) (toEnum . fromEnum $ doPolarizationCorrection)
 207       (ImageWord32 arr) -> unsafeWith arr $ \i -> do
 208         {-# SCC "hkl_binoculars_space_hkl_uint32_t" #-} c'hkl_binoculars_space_hkl_uint32_t pSpace geometry sample i nPixels (CDouble . unAttenuation $ att) pix (toEnum ndim) dims r (toEnum nr) mask'' limits (toEnum nlimits) (toEnum . fromEnum $ doPolarizationCorrection)
 209     return (DataFrameSpace img space att)
 210
 211 ----------
 212 -- Pipe --
 213 ----------
 214
 215 processHklP :: (MonadIO m, MonadLogger m, MonadReader (Config 'HklProjection) m, MonadThrow m)
 216             => m ()
 217 processHklP = do
 218   conf :: Config 'HklProjection <- ask
 219
 220   -- directly from the common config
 221   let common = binocularsConfig'Hkl'Common conf
 222
 223   let overwrite = binocularsConfig'Common'Overwrite common
 224   let det = binocularsConfig'Common'Detector common
 225   let (NCores cap) =  binocularsConfig'Common'NCores common
 226   let destination = binocularsConfig'Common'Destination common
 227   let centralPixel' = binocularsConfig'Common'Centralpixel common
 228   let (Meter sampleDetectorDistance) = binocularsConfig'Common'Sdd common
 229   let (Degree detrot) = binocularsConfig'Common'Detrot common
 230   let mImageSumMax = binocularsConfig'Common'ImageSumMax common
 231   let inputRange = binocularsConfig'Common'InputRange common
 232   let nexusDir = binocularsConfig'Common'Nexusdir common
 233   let tmpl = binocularsConfig'Common'Tmpl common
 234   let maskMatrix = binocularsConfig'Common'Maskmatrix common
 235   let mSkipFirstPoints = binocularsConfig'Common'SkipFirstPoints common
 236   let mSkipLastPoints = binocularsConfig'Common'SkipLastPoints common
 237   let doPolarizationCorrection = binocularsConfig'Common'PolarizationCorrection common
 238
 239   -- directly from the specific config
 240   let mlimits = binocularsConfig'Hkl'ProjectionLimits conf
 241   let res = binocularsConfig'Hkl'ProjectionResolution conf
 242   let datapaths = binocularsConfig'Hkl'DataPath conf
 243   let projectionType = binocularsConfig'Hkl'ProjectionType conf
 244
 245
 246   -- built from the config
 247   output' <- liftIO $ destination' projectionType Nothing inputRange mlimits destination overwrite
 248   filenames <- InputFn'List <$> files nexusDir (Just inputRange) tmpl
 249   mask' <- getMask maskMatrix det
 250   pixels <- liftIO $ getPixelsCoordinates det centralPixel' sampleDetectorDistance detrot NoNormalisation
 251
 252   let fns = concatMap (replicate 1) (toList filenames)
 253   chunks <- liftIO $ runSafeT $ toListM $ each fns >-> chunkP mSkipFirstPoints mSkipLastPoints datapaths
 254   let ntot = sum (Prelude.map clength chunks)
 255   let jobs = chunk (quot ntot cap) chunks
 256
 257   -- log parameters
 258
 259   logDebugNSH filenames
 260   logDebugNSH datapaths
 261   logDebugNSH chunks
 262   logDebugN "start gessing final cube size"
 263
 264   -- guess the final cube dimensions (To optimize, do not create the cube, just extract the shape)
 265
 266   guessed <- liftIO $ withCubeAccumulator EmptyCube $ \c ->
 267     runSafeT $ runEffect $
 268     each chunks
 269     >-> Pipes.Prelude.map (\(Chunk fn f t) -> (fn, [f, quot (f + t) 4, quot (f + t) 4 * 2, quot (f + t) 4 * 3, t]))
 270     >-> framesP datapaths
 271     >-> project det 3 (spaceHkl det pixels res mask' mlimits doPolarizationCorrection)
 272     >-> accumulateP c
 273
 274   logDebugN "stop gessing final cube size"
 275
 276   logInfoN (pack $ printf "let's do an Hkl projection of %d %s image(s) on %d core(s)" ntot (show det) cap)
 277
 278   liftIO $ withProgressBar ntot $ \pb -> do
 279     r' <- mapConcurrently (\job -> withCubeAccumulator guessed $ \c ->
 280                              runEffect $ runSafeP $
 281                              each job
 282                              >-> Pipes.Prelude.map (\(Chunk fn f t) -> (fn, [f..t]))
 283                              -- >-> tee Pipes.Prelude.print
 284                              >-> framesP datapaths
 285                              >-> Pipes.Prelude.filter (\(DataFrameHkl (DataFrameQCustom _ _ img _) _) -> filterSumImage mImageSumMax img)
 286                              >-> project det 3 (spaceHkl det pixels res mask' mlimits doPolarizationCorrection)
 287                              >-> tee (accumulateP c)
 288                              >-> progress pb
 289                          ) jobs
 290     saveCube output' (unpack . serializeConfig $ conf) r'
 291
 292 -- FramesHklP
 293
 294 instance ChunkP (DataFrameHkl' DataSourcePath) where
 295   chunkP sf sl (DataFrameHkl p _) = chunkP sf sl p
 296
 297 instance FramesP (DataFrameHkl' DataSourcePath) (DataFrameHkl' Identity) where
 298   framesP (DataFrameHkl qcustom sample) = skipMalformed $ forever $ do
 299     (fp, js) <- await
 300     withFileP (openFile' fp) $ \f ->
 301       withDataSourceP f qcustom $ \qcustomAcq ->
 302       withDataSourceP f sample $ \sampleAcq ->
 303       forM_ js (\j -> tryYield ( DataFrameHkl
 304                                 <$> extract1DStreamValue qcustomAcq j
 305                                 <*> extract0DStreamValue sampleAcq
 306                               ))
 307
 308 ------------
 309 -- Inputs --
 310 ------------
 311
 312 guess'DataSourcePath'DataFrameHkl :: BinocularsConfig'Common
 313                                   -> BinocularsConfig'Sample
 314                                   -> ConfigContent
 315                                   -> DataFrameHkl' DataSourcePath
 316 guess'DataSourcePath'DataFrameHkl common sample content
 317   = DataFrameHkl
 318     (guess'DataSourcePath'DataFrameQCustom common Nothing content)
 319     (guess'DataSourcePath'Sample common sample)
 320
 321 ---------
 322 -- Cmd --
 323 ---------
 324
 325 processHkl :: (MonadLogger m, MonadThrow m, MonadIO m) => Maybe FilePath -> Maybe ConfigRange -> m ()
 326 processHkl mf mr = cmd processHklP mf (Args'HklProjection mr)
 327
 328 newHkl :: (MonadIO m, MonadLogger m, MonadThrow m)
 329        => Path Abs Dir -> m ()
 330 newHkl cwd = do
 331   let conf = default'BinocularsConfig'Hkl
 332              { binocularsConfig'Hkl'Common = default'BinocularsConfig'Common
 333                                              { binocularsConfig'Common'Nexusdir = Just cwd }
 334              }
 335   liftIO $ Data.Text.IO.putStr $ serializeConfig conf
 336
 337
 338 updateHkl :: (MonadIO m, MonadLogger m, MonadThrow m)
 339           => Maybe FilePath -> Maybe ConfigRange -> m ()
 340 updateHkl mf mr = cmd (pure ()) mf (Args'HklProjection mr)