diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,7 @@ +- V 0.2.1.0: Added a mechanism to detect terminal encoding and fall back on a simpler CLI plot if it is not UTF-8 +- V 0.2.0.1: Brought sample names back to default CLI output +- V 0.2.0.0: Added sum (and product) calibration and made the necessary changes to various interfaces (including CLI) to make this functionality accessible +- V 0.1.2.0: Added simple summary data (CalRangeSummary with calibrated median age + begin and end of 1- and 2-sigma ranges) to CalC14 and the cli output and plot. The latter got refactored and enhanced in the process. HDRs are now "ordered", so _hdrstart actually stores the older and _hdrstop the younger date - V 0.1.1.0: Complete rewrite of the cli output handling to avoid a memory leak - V 0.1.0.0: Switch to PVP versioning (https://pvp.haskell.org/) - V 0.24.4: Removed big dependencies bytestring and statistics diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ ### CLI app -For stable release versions we automatically prepare binaries that can be downloaded and run. +For stable release versions we automatically prepare statically built binaries that can be downloaded and run directly. You can download them here: [ [Linux 📥](https://github.com/nevrome/currycarbon/releases/latest/download/currycarbon-Linux) | [macOS 📥](https://github.com/nevrome/currycarbon/releases/latest/download/currycarbon-macOS) | [Windows 📥](https://github.com/nevrome/currycarbon/releases/latest/download/currycarbon-Windows.exe) ]. Older release versions are available [here](https://github.com/nevrome/currycarbon/releases). @@ -25,11 +25,34 @@ # make it executable chmod +x currycarbon-Linux # test it -./currycarbon-Linux "Sample1,4000,30" +./currycarbon-Linux "Sample1,4990,30" ``` ``` -Usage: currycarbon [--version] [DATES] [-i|--inputFile ARG] +currycarbon v0.2.1.0 (UTF-8) +Method: Bchron {distribution = StudentTDist {ndf = 100.0}} +Curve: IntCal20 +Calibrating... +DATE: Sample1:4990±30BP +Calibrated: 3936BC >> 3794BC > 3757BC < 3662BC << 3654BC +1-sigma: 3794-3707BC, 3666-3662BC +2-sigma: 3936-3874BC, 3804-3697BC, 3684-3654BC + ▁▁▁ ▁▁▁▁ + ▁▁▒▒▒▁▁▁▁▒▒▒▒▁ + ▒▒▒▒▒▒▒▒▒▒▒▒▒▒ ▁ + ▁▁ ▁▒▒▒▒▒▒▒▒▒▒▒▒▒▒▁ ▁▒ + ▁▁▒▒▁ ▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒ ▁▒▒▁ + ▁▁▁▁▁▒▒▒▒▒▁ ▁▁▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▁▁▁▁▒▒▒▒▁ + ▁▁▁▒▒▒▒▒▒▒▒▒▒▒▁▁▁▁▁▁▁▁▁▁▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▁ + -3950 ┄─────────┬───────────────┬────────────────┬─────────┄ -3640 + > > ^ < < + ──────────────── ─ + ─────────── ────────────────── ────── +Done. +``` + +``` +Usage: currycarbon [--version] [DATE] [-i|--inputFile ARG] [--calibrationCurveFile ARG] [--method ARG] [--allowOutside] [--noInterpolation] [-q|--quiet] [--densityFile ARG] [--hdrFile ARG] [--calCurveSegmentFile ARG] @@ -39,15 +62,21 @@ Available options: -h,--help Show this help text --version Show version - DATES A string with one or multiple uncalibrated dates of + DATE A string with one or multiple uncalibrated dates of the form ",,;..." where is - optional. So for example - "S1,4000,50;3000,25;S3,1000,20". - -i,--inputFile ARG A file with a list of uncalibrated dates. Formated - just as DATES, but with a new line for each input - date. DATES and --inputFile can be combined and you - can provide multiple instances of --inputFile + standard deviation>" where is optional + (e.g. "S1,4000,50"). Multiple dates can be listed + separated by ";" (e.g. "S1,4000,50; 3000,25; + S3,1000,20"). To sum or multiply the post calibration + probability distributions, dates can be combined with + "+" or "*" (e.g. "4000,50 + 4100,100"). These + expressions can be combined arbitrarily. Parentheses + can be added to specify the order of operations (e.g. + "(4000,50 + 4100,100) * 3800,50") + -i,--inputFile ARG A file with a list of calibration expressions. + Formated just as DATE, but with a new line for each + input date. DATE and --inputFile can be combined and + you can provide multiple instances of --inputFile --calibrationCurveFile ARG Path to an calibration curve file in .14c format. The calibration curve will be read and used for @@ -82,9 +111,9 @@ (!) input date in a wide matrix format ``` -### For developers +### For developers who want to edit the code -To install the latest development version you can follow these steps: +To install the latest development version (!) you can follow these steps: 1. Install the Haskell build tool [Stack](https://docs.haskellstack.org/en/stable/README/) 2. Clone the repository diff --git a/currycarbon.cabal b/currycarbon.cabal --- a/currycarbon.cabal +++ b/currycarbon.cabal @@ -1,5 +1,5 @@ name: currycarbon -version: 0.1.1.0 +version: 0.2.1.0 synopsis: A package for simple, fast radiocarbon calibration description: Radiocarbon calibration with the intercept method optimised for fast calibration of many dates. homepage: https://github.com/nevrome/currycarbon @@ -27,6 +27,7 @@ Currycarbon.Calibration.Calibration Currycarbon.CLI.RunCalibrate Currycarbon.Parsers + Currycarbon.SumCalibration Currycarbon.Types Currycarbon.Utils hs-source-dirs: @@ -55,20 +56,4 @@ default-language: Haskell2010 ghc-options: - -threaded -with-rtsopts=-N - -Test-Suite currycarbon-test - type: - exitcode-stdio-1.0 - main-is: - Spec.hs - hs-source-dirs: - test - build-depends: - currycarbon - , base - , doctest >= 0.17 && < 0.18 - default-language: - Haskell2010 - ghc-options: -threaded -with-rtsopts=-N diff --git a/src-executables/Main-currycarbon.hs b/src-executables/Main-currycarbon.hs --- a/src-executables/Main-currycarbon.hs +++ b/src-executables/Main-currycarbon.hs @@ -11,8 +11,7 @@ import Data.Version (showVersion) import qualified Options.Applicative as OP import System.Exit (exitFailure) -import System.IO (hPutStrLn, stderr) - +import System.IO (hPutStrLn, stderr, stdout, hGetEncoding) -- * CLI interface configuration -- -- $cliInterface @@ -25,8 +24,13 @@ -- CLI interface configuration main :: IO () main = do + -- check stdout encoding for the CLI plot + stdOutEncoding <- hGetEncoding stdout + let encoding = maybe "unknown" show stdOutEncoding + hPutStrLn stderr $ "currycarbon v" ++ showVersion version ++ " (" ++ encoding ++ ")" + -- prepare input parsing cmdOpts <- OP.customExecParser p optParserInfo - catch (runCmd cmdOpts) handler + catch (runCmd encoding cmdOpts) handler where p = OP.prefs OP.showHelpOnEmpty handler :: CurrycarbonException -> IO () @@ -34,9 +38,9 @@ hPutStrLn stderr $ renderCurrycarbonException e exitFailure -runCmd :: Options -> IO () -runCmd o = case o of - CmdCalibrate opts -> runCalibrate opts +runCmd :: String -> Options -> IO () +runCmd enc o = case o of + CmdCalibrate opts -> runCalibrate opts {_calibrateStdOutEncoding = enc} optParserInfo :: OP.ParserInfo Options optParserInfo = OP.info (OP.helper <*> versionOption <*> optParser) ( @@ -51,13 +55,14 @@ optParser = CmdCalibrate <$> calibrateOptParser calibrateOptParser :: OP.Parser CalibrateOptions -calibrateOptParser = CalibrateOptions <$> optParseUncalC14String - <*> optParseUncalC14FromFile +calibrateOptParser = CalibrateOptions <$> optParseCalExprString + <*> optParseCalExprFromFile <*> optParseCalCurveFromFile <*> optParseCalibrationMethod <*> optParseAllowOutside <*> optParseDontInterpolateCalCurve <*> optParseQuiet + <*> pure "unknown" <*> optParseDensityFile <*> optParseHDRFile <*> optParseCalCurveSegmentFile @@ -69,22 +74,26 @@ -- -- These functions define and handle the CLI input arguments -optParseUncalC14String :: OP.Parser [UncalC14] -optParseUncalC14String = concat <$> OP.many (OP.argument (OP.eitherReader readUncalC14) ( - OP.metavar "DATES" <> +optParseCalExprString :: OP.Parser [CalExpr] +optParseCalExprString = concat <$> OP.many (OP.argument (OP.eitherReader readCalExpr) ( + OP.metavar "DATE" <> OP.help "A string with one or multiple uncalibrated dates of \ - \the form \",,;...\" \ - \where is optional. \ - \So for example \"S1,4000,50;3000,25;S3,1000,20\"." + \the form \",,\" \ + \where is optional (e.g. \"S1,4000,50\"). \ + \Multiple dates can be listed separated by \";\" (e.g. \"S1,4000,50; 3000,25; S3,1000,20\"). \ + \To sum or multiply the post calibration probability distributions, dates can be combined with \ + \\"+\" or \"*\" (e.g. \"4000,50 + 4100,100\"). \ + \These expressions can be combined arbitrarily. Parentheses can be added to specify the order \ + \of operations (e.g. \"(4000,50 + 4100,100) * 3800,50\")" )) -optParseUncalC14FromFile :: OP.Parser [FilePath] -optParseUncalC14FromFile = OP.many (OP.strOption ( +optParseCalExprFromFile :: OP.Parser [FilePath] +optParseCalExprFromFile = OP.many (OP.strOption ( OP.long "inputFile" <> OP.short 'i' <> - OP.help "A file with a list of uncalibrated dates. \ - \Formated just as DATES, but with a new line for each input date. \ - \DATES and --inputFile can be combined and you can provide multiple instances of --inputFile" + OP.help "A file with a list of calibration expressions. \ + \Formated just as DATE, but with a new line for each input date. \ + \DATE and --inputFile can be combined and you can provide multiple instances of --inputFile" )) optParseCalCurveFromFile :: OP.Parser (Maybe FilePath) diff --git a/src/Currycarbon.hs b/src/Currycarbon.hs --- a/src/Currycarbon.hs +++ b/src/Currycarbon.hs @@ -37,13 +37,23 @@ refineCalDates, CalC14 (..), writeCalC14s, - renderCalDatesPretty + renderCalDatePretty, + + -- * Sum (and product) calibration + -- $sumcal + evalCalExpr, + CalExpr (..), + addPDFs, + multiplyPDFs, + normalizeCalPDF ) where import Currycarbon.Calibration.Calibration import Currycarbon.Parsers +import Currycarbon.SumCalibration import Currycarbon.Types import Currycarbon.CalCurves.Intcal20 +import Currycarbon.Calibration.Utils (normalizeCalPDF) {- $calibration @@ -104,4 +114,17 @@ 'renderCalDatesPretty' finally combines 'UncalC14', 'CalPDF' and 'CalC14' to produce nice command line output summarising the calibration result for a given sample. +-} + +{- $sumcal + +Calculating the sum or product of two calibration curves is a common +application, which currycarbon supports with a custom algebraic data type +'CalExpr'. It encodes a language to describe (very simple) chronological +models, to be evaluated to a single 'CalPDF' with 'evalCalExpr'. + +A more basic interface is available with 'addPDFs' and 'multiplyPDFs', +which allow to combine two 'CalPDF's with the respective operation. +Depending on the application, 'normalizeCalPDF' will come in handy here, +to normalize the output density distributions. -} diff --git a/src/Currycarbon/CLI/RunCalibrate.hs b/src/Currycarbon/CLI/RunCalibrate.hs --- a/src/Currycarbon/CLI/RunCalibrate.hs +++ b/src/Currycarbon/CLI/RunCalibrate.hs @@ -6,24 +6,24 @@ import Currycarbon.CalCurves.Intcal20 import Currycarbon.Calibration.Calibration import Currycarbon.Parsers +import Currycarbon.SumCalibration import Currycarbon.Types import Currycarbon.Utils import Control.Monad (when, unless) -import Data.Either (rights, lefts, isRight) -import Data.Foldable (forM_) -import Data.Maybe (fromJust, isJust) -import System.IO (hPutStrLn, stderr) +import Data.Maybe (fromJust, isJust, fromMaybe) +import System.IO (hPutStrLn, stderr, stdout) -- | A data type to represent the options to the CLI module function runCalibrate data CalibrateOptions = CalibrateOptions { - _calibrateUncalC14 :: [UncalC14] -- ^ Uncalibrated dates that should be calibrated - , _calibrateUncalC14File :: [FilePath] -- ^ List of files with uncalibrated dates to be calibrated + _calibrateExprs :: [CalExpr] -- ^ String listing the uncalibrated dates that should be calibrated + , _calibrateExprFiles :: [FilePath] -- ^ List of files with uncalibrated dates to be calibrated , _calibrateCalCurveFile :: Maybe FilePath -- ^ Path to a .14c file , _calibrateCalibrationMethod :: CalibrationMethod -- ^ Calibration algorithm that should be used , _calibrateAllowOutside :: Bool -- ^ Allow calibration to run outside of the range of the calibration curve , _calibrateDontInterpolateCalCurve :: Bool -- ^ Don't interpolate the calibration curve , _calibrateQuiet :: Bool -- ^ Suppress the printing of calibration results to the command line + , _calibrateStdOutEncoding :: String -- ^ Encoding of the stdout stream (show TextEncoding) , _calibrateDensityFile :: Maybe FilePath -- ^ Path to an output file (see CLI documentation) , _calibrateHDRFile :: Maybe FilePath -- ^ Path to an output file , _calibrateCalCurveSegmentFile :: Maybe FilePath -- ^ Path to an output file @@ -32,69 +32,86 @@ -- | Interface function to trigger calibration from the command line runCalibrate :: CalibrateOptions -> IO () -runCalibrate (CalibrateOptions uncalDates uncalFile calCurveFile method allowOutside noInterpolate quiet densityFile hdrFile calCurveSegmentFile calCurveMatrixFile) = do +runCalibrate (CalibrateOptions exprs exprFiles calCurveFile method allowOutside noInterpolate quiet encoding densityFile hdrFile calCurveSegmentFile calCurveMatrixFile) = do + let ascii = encoding /= "UTF-8" -- compile dates - entitiesFromFile <- mapM readUncalC14FromFile uncalFile - let uncalDatesRenamed = replaceEmptyNames $ uncalDates ++ concat entitiesFromFile - if null uncalDatesRenamed + exprsFromFile <- mapM readCalExprFromFile exprFiles + let exprsRenamed = replaceEmptyNames $ exprs ++ concat exprsFromFile + if null exprsRenamed then hPutStrLn stderr "Nothing to calibrate. See currycarbon -h for help" else do - -- calibration - hPutStrLn stderr "Calibrating..." + -- prep data + hPutStrLn stderr $ "Method: " ++ show method + hPutStrLn stderr $ "Curve: " ++ fromMaybe "IntCal20" calCurveFile calCurve <- maybe (return intcal20) readCalCurveFromFile calCurveFile let calConf = defaultCalConf { _calConfMethod = method , _calConfAllowOutside = allowOutside , _calConfInterpolateCalCurve = not noInterpolate } - let errorOrCalPDFs = calibrateDates calConf calCurve uncalDatesRenamed - handleDates True calCurve $ zip uncalDatesRenamed errorOrCalPDFs - where - handleDates :: Bool -> CalCurveBP -> [(UncalC14, Either CurrycarbonException CalPDF)] -> IO () - handleDates _ _ [] = hPutStrLn stderr "Done." - handleDates True calCurve (x:xs) = case x of - (_, Left ex) -> printEx ex >> handleDates True calCurve xs - (uncal, Right cPDF) -> handleFirstDate calCurve uncal cPDF >> handleDates False calCurve xs - handleDates False calCurve (x:xs) = case x of - (_, Left ex) -> printEx ex >> handleDates False calCurve xs - (uncal, Right cPDF) -> handleOtherDate uncal cPDF >> handleDates False calCurve xs - handleFirstDate :: CalCurveBP -> UncalC14 -> CalPDF -> IO () - handleFirstDate calCurve uncal calPDF = do - -- calcurve segment or calcurve matrix file - if isJust calCurveSegmentFile || isJust calCurveMatrixFile - then do - hPutStrLn stderr $ - "The calCurveSegment file and the calCurveMatrix file only consider the first date, " ++ - renderUncalC14 uncal - let calCurveSegment = prepareCalCurveSegment (not noInterpolate) $ getRelevantCalCurveSegment uncal calCurve - when (isJust calCurveSegmentFile) $ - writeCalCurve (fromJust calCurveSegmentFile) calCurveSegment - when (isJust calCurveMatrixFile) $ - writeCalCurveMatrix (fromJust calCurveMatrixFile) $ - makeCalCurveMatrix (uncalToPDF uncal) calCurveSegment - else do - -- other output - let calC14 = refineCalDate calPDF - unless quiet $ putStrLn $ renderCalDatePretty (uncal, calPDF, calC14) - when (isJust hdrFile) $ writeCalC14 (fromJust hdrFile) calC14 - when (isJust densityFile) $ writeCalPDF (fromJust densityFile) calPDF - handleOtherDate :: UncalC14 -> CalPDF -> IO () - handleOtherDate uncal calPDF = do - let calC14 = refineCalDate calPDF - unless quiet $ putStrLn $ renderCalDatePretty (uncal, calPDF, calC14) - when (isJust hdrFile) $ appendCalC14 (fromJust hdrFile) calC14 - when (isJust densityFile) $ appendCalPDF (fromJust densityFile) calPDF - printEx :: CurrycarbonException -> IO () - printEx ex = hPutStrLn stderr $ renderCurrycarbonException ex + -- run calibration + hPutStrLn stderr "Calibrating..." + let errorOrCalPDFs = map (evalCalExpr calConf calCurve) exprsRenamed + handleDates ascii True calCurve $ zip exprsRenamed errorOrCalPDFs + where + -- the bool manages if a date is the first, calibratable date + handleDates :: Bool -> Bool -> CalCurveBP -> [(CalExpr, Either CurrycarbonException CalPDF)] -> IO () + handleDates _ _ _ [] = hPutStrLn stderr "Done." + handleDates _ascii True calCurve (firstDate:otherDates) = case firstDate of + (_, Left e) -> printE e >> handleDates _ascii True calCurve otherDates + (calExpr, Right cPDF) -> firstOut _ascii calCurve calExpr cPDF >> handleDates _ascii False calCurve otherDates + handleDates _ascii False calCurve (firstDate:otherDates) = case firstDate of + (_, Left e) -> printE e >> handleDates _ascii False calCurve otherDates + (calExpr, Right cPDF) -> otherOut _ascii calExpr cPDF >> handleDates _ascii False calCurve otherDates + firstOut :: Bool -> CalCurveBP -> CalExpr -> CalPDF -> IO () + firstOut _ascii calCurve calExpr@(UnCalDate uncal) calPDF = do + flexOut _ascii calExpr calPDF writeCalPDF writeCalC14 + when (isJust calCurveSegmentFile || isJust calCurveMatrixFile) $ do + hPutStrLn stderr $ + "Warning: The calCurveSegment file and the calCurveMatrix file only consider the first date, " ++ + renderUncalC14 uncal + let calCurveSegment = prepareCalCurveSegment (not noInterpolate) $ getRelevantCalCurveSegment uncal calCurve + when (isJust calCurveSegmentFile) $ + writeCalCurve (fromJust calCurveSegmentFile) calCurveSegment + when (isJust calCurveMatrixFile) $ + writeCalCurveMatrix (fromJust calCurveMatrixFile) $ + makeCalCurveMatrix (uncalToPDF uncal) calCurveSegment + firstOut _ascii _ calExpr calPDF = do + flexOut _ascii calExpr calPDF writeCalPDF writeCalC14 + when (isJust calCurveSegmentFile || isJust calCurveMatrixFile) $ do + hPutStrLn stderr $ "Warning: The calCurveSegment file and the calCurveMatrix file can only be produced for simple dates" + otherOut :: Bool -> CalExpr -> CalPDF -> IO () + otherOut _ascii calExpr calPDF = + flexOut _ascii calExpr calPDF appendCalPDF appendCalC14 + flexOut :: Bool -> CalExpr -> CalPDF -> (FilePath -> CalPDF -> IO ()) -> (FilePath -> CalC14 -> IO ()) -> IO () + flexOut _ascii calExpr calPDF calPDFToFile calC14ToFile = do + case refineCalDate calPDF of + Nothing -> do + unless quiet $ do + hPutStrLn stdout $ renderCalExpr calExpr + hPutStrLn stderr "Warning: Could not calculate meaningful HDRs for this expression. Check --densityFile." + when (isJust hdrFile) $ unless quiet $ hPutStrLn stderr "Nothing written to the HDR file" + when (isJust densityFile) $ calPDFToFile (fromJust densityFile) calPDF + Just calC14 -> do + unless quiet $ hPutStrLn stdout $ renderCalDatePretty _ascii (calExpr, calPDF, calC14) + when (isJust hdrFile) $ calC14ToFile (fromJust hdrFile) calC14 + when (isJust densityFile) $ calPDFToFile (fromJust densityFile) calPDF + printE :: CurrycarbonException -> IO () + printE e = hPutStrLn stderr $ renderCurrycarbonException e -- | Helper function to replace empty input names with a sequence of numbers, -- to get each input date an unique identifier -replaceEmptyNames :: [UncalC14] -> [UncalC14] -replaceEmptyNames xs = - zipWith replaceName xs [1..] +replaceEmptyNames :: [CalExpr] -> [CalExpr] +replaceEmptyNames = zipWith (replaceName . show) ([1..] :: [Integer]) where - replaceName :: UncalC14 -> Int -> UncalC14 - replaceName (UncalC14 name mean std) number = + replaceName :: String -> CalExpr -> CalExpr + replaceName i (UnCalDate (UncalC14 name x y)) = if name == "unknownSampleName" - then UncalC14 (show number) mean std - else UncalC14 name mean std + then UnCalDate $ UncalC14 i x y + else UnCalDate $ UncalC14 name x y + replaceName i (CalDate (CalPDF name x y)) = + if name == "unknownSampleName" + then CalDate $ CalPDF i x y + else CalDate $ CalPDF name x y + replaceName i (SumCal a b) = SumCal (replaceName (i ++ "s") a) (replaceName (i ++ "S") b) + replaceName i (ProductCal a b) = ProductCal (replaceName (i ++ "p") a) (replaceName (i ++ "P") b) diff --git a/src/Currycarbon/Calibration/Bchron.hs b/src/Currycarbon/Calibration/Bchron.hs --- a/src/Currycarbon/Calibration/Bchron.hs +++ b/src/Currycarbon/Calibration/Bchron.hs @@ -2,8 +2,6 @@ module Currycarbon.Calibration.Bchron (calibrateDateBchron) where - -import Currycarbon.CalCurves.Intcal20 -- only for the doctest import Currycarbon.Calibration.Utils import Currycarbon.Parsers import Currycarbon.Types @@ -12,13 +10,6 @@ import qualified Data.Vector.Unboxed as VU -- | Intercept calibration as implemented in the Bchron R package (see 'Bchron') --- --- >>> :{ --- let calResult = calibrateDateBchron NormalDist False True intcal20 (UncalC14 "" 3000 30) --- simplifyResult = \(Right (CalPDF _ age dens)) -> VU.zip (VU.take 3 age) (VU.take 3 dens) --- in simplifyResult calResult --- :} --- [(-1413,1.0703295e-5),(-1412,1.4653518e-5),(-1411,1.799142e-5)] calibrateDateBchron :: CalibrationDistribution -> Bool -> Bool -> CalCurveBP -> UncalC14 -> Either CurrycarbonException CalPDF calibrateDateBchron distr allowOutside interpolate calCurve uncalC14@(UncalC14 name age ageSd) = if not allowOutside && isOutsideRangeOfCalCurve calCurve uncalC14 diff --git a/src/Currycarbon/Calibration/Calibration.hs b/src/Currycarbon/Calibration/Calibration.hs --- a/src/Currycarbon/Calibration/Calibration.hs +++ b/src/Currycarbon/Calibration/Calibration.hs @@ -10,6 +10,7 @@ , prepareCalCurveSegment , makeCalCurveMatrix , uncalToPDF + , calibrateDate , calibrateDates , refineCalDates , refineCalDate @@ -23,7 +24,8 @@ import Currycarbon.Types import Currycarbon.Utils -import Data.List (sort, sortBy, groupBy) +import Data.List (sort, sortBy, groupBy, elemIndex) +import Data.Maybe (fromJust) import qualified Data.Vector.Unboxed as VU -- | A data type to cover the configuration options of the calibrateDates function @@ -43,9 +45,6 @@ -- | A default configuration that should yield almost identical calibration results -- to the [Bchron R package](https://github.com/andrewcparnell/Bchron) --- --- >>> defaultCalConf --- CalibrateDatesConf {_calConfMethod = Bchron {distribution = StudentTDist {ndf = 100.0}}, _calConfAllowOutside = False, _calConfInterpolateCalCurve = True} defaultCalConf :: CalibrateDatesConf defaultCalConf = CalibrateDatesConf { _calConfMethod = Bchron { distribution = StudentTDist 100 } @@ -66,27 +65,55 @@ calibrateDates (CalibrateDatesConf Bchron{distribution=distr} allowOutside interpolate) calCurve uncalDates = map (calibrateDateBchron distr allowOutside interpolate calCurve) uncalDates +-- | Calibrates a date with the provided calibration curve +calibrateDate :: CalibrateDatesConf -- ^ Configuration options to consider + -> CalCurveBP -- ^ A calibration curve + -> UncalC14 -- ^ An uncalibrated radiocarbon date + -> Either CurrycarbonException CalPDF -- ^ The function returns either an exception if the + -- calibration failed for some reason, or a 'CalPDF' +calibrateDate (CalibrateDatesConf MatrixMultiplication allowOutside interpolate) calCurve uncalDate = + calibrateDateMatrixMult allowOutside interpolate calCurve uncalDate +calibrateDate (CalibrateDatesConf Bchron{distribution=distr} allowOutside interpolate) calCurve uncalDate = + calibrateDateBchron distr allowOutside interpolate calCurve uncalDate + -- | Transforms the raw, calibrated probability density table to a meaningful representation of a -- calibrated radiocarbon date --- --- >>> :{ --- let years = VU.fromList [-3000,-2900,-2800,-2700,-2600] --- densities = VU.fromList [0.1,0.2,0.3,0.2,0.1] --- in refineCalDates [CalPDF {_calPDFid = "", _calPDFCals = years, _calPDFDens = densities}] --- :} --- [CalC14 {_calC14id = "", _calC14HDROneSigma = [HDR {_hdrstart = -2800, _hdrstop = -2900}], _calC14HDRTwoSigma = [HDR {_hdrstart = -2600, _hdrstop = -3000}]}] -refineCalDates :: [CalPDF] -> [CalC14] +refineCalDates :: [CalPDF] -> [Maybe CalC14] refineCalDates = map refineCalDate -refineCalDate :: CalPDF -> CalC14 +refineCalDate :: CalPDF -> Maybe CalC14 refineCalDate (CalPDF name cals dens) = - let sortedDensities = sortBy (flip (\ (_, dens1) (_, dens2) -> compare dens1 dens2)) (VU.toList $ VU.zip cals dens) - cumsumDensities = scanl1 (+) $ map snd sortedDensities - isIn68 = map (< 0.683) cumsumDensities - isIn95 = map (< 0.954) cumsumDensities - contextualizedDensities = reverse $ sort $ zipWith3 (\(y,d) in68 in95 -> (y,d,in68,in95)) sortedDensities isIn68 isIn95 - in CalC14 name (densities2HDR68 contextualizedDensities) (densities2HDR95 contextualizedDensities) + if VU.sum dens == 0 || VU.length (VU.filter (>= 1.0) dens) == 1 -- don't calculate CalC14, if it's not meaningful + then Nothing + else Just $ CalC14 { + _calC14id = name + , _calC14RangeSummary = CalRangeSummary { + _calRangeStartTwoSigma = _hdrstart $ head hdrs95 + , _calRangeStartOneSigma = _hdrstart $ head hdrs68 + , _calRangeMedian = fromJust $ cals `indexVU` elemIndex (minimum distanceTo05) distanceTo05 + , _calRangeStopOneSigma = _hdrstop $ last hdrs68 + , _calRangeStopTwoSigma = _hdrstop $ last hdrs95 + } + , _calC14HDROneSigma = hdrs68 + , _calC14HDRTwoSigma = hdrs95 + } where + -- simple density cumsum for median age + cumsumDensities = cumsumDens (VU.toList $ VU.zip cals dens) + distanceTo05 = map (\x -> abs $ (x - 0.5)) cumsumDensities + -- sorted density cumsum for hdrs + sortedDensities = sortBy (flip (\ (_, dens1) (_, dens2) -> compare dens1 dens2)) (VU.toList $ VU.zip cals dens) + cumsumSortedDensities = cumsumDens sortedDensities + isIn68 = map (< 0.683) cumsumSortedDensities + isIn95 = map (< 0.954) cumsumSortedDensities + contextualizedDensities = sort $ zipWith3 (\(y,d) in68 in95 -> (y,d,in68,in95)) sortedDensities isIn68 isIn95 + hdrs68 = densities2HDR68 contextualizedDensities + hdrs95 = densities2HDR95 contextualizedDensities + -- helper functions + indexVU _ Nothing = Nothing + indexVU x (Just i) = x VU.!? i + cumsumDens :: [(YearBCAD, Float)] -> [Float] + cumsumDens x = scanl1 (+) $ map snd x densities2HDR68 :: [(Int, Float, Bool, Bool)] -> [HDR] densities2HDR68 cDensities = let highDensityGroups = groupBy (\(_,_,in681,_) (_,_,in682,_) -> in681 == in682) cDensities diff --git a/src/Currycarbon/Calibration/Utils.hs b/src/Currycarbon/Calibration/Utils.hs --- a/src/Currycarbon/Calibration/Utils.hs +++ b/src/Currycarbon/Calibration/Utils.hs @@ -8,10 +8,14 @@ import Data.Maybe (fromMaybe) import Numeric.SpecFunctions (logBeta) +-- | Rescale a CalPDF so that the sum of the densities is approx. 1.0 +normalizeCalPDF :: CalPDF -> CalPDF +normalizeCalPDF (CalPDF name cals dens) = + case VU.sum dens of + 0.0 -> CalPDF name cals dens -- product calibration can yield empty calPDFs + s -> CalPDF name cals $ VU.map (/s) dens + -- | get the density of a normal distribution at a point x --- --- >>> dnorm 1.0 1.0 1.0 --- 0.3989423 dnorm :: Float -> Float -> Float -> Float dnorm mu sigma x = let a = recip (sqrt (2 * pi * sigma2)) @@ -25,9 +29,6 @@ -- realToFrac $ density (normalDistr (realToFrac mu) (realToFrac sigma)) (realToFrac x) -- | get the density of student's-t distribution at a point x --- --- >>> dt 1.0 1.0 --- 0.15915494 dt :: Double -> Float -> Float dt dof x = let xDouble = realToFrac x @@ -95,12 +96,6 @@ in (xPred, y1 + xPredRel * yDiffPerxDiff) uncurry3 :: (a -> b -> c -> d) -> ((a, b, c) -> d) uncurry3 f ~(a,b,c) = f a b c - -normalizeCalPDF :: CalPDF -> CalPDF -normalizeCalPDF (CalPDF name cals dens) = - let sumDens = VU.sum dens - normalizedDens = VU.map (/ sumDens) dens - in CalPDF name cals normalizedDens trimLowDensityEdgesCalPDF :: CalPDF -> CalPDF trimLowDensityEdgesCalPDF (CalPDF name cals dens) = diff --git a/src/Currycarbon/Parsers.hs b/src/Currycarbon/Parsers.hs --- a/src/Currycarbon/Parsers.hs +++ b/src/Currycarbon/Parsers.hs @@ -19,34 +19,6 @@ -- This module contains a number of functions to manage data input and -- output plumbing for different datatypes --- | Combine 'UncalC14', 'CalPDF' and 'CalC14' to render pretty command line output --- like this: --- --- @ --- Sample: 1 ~\> [5000±30BP] --- 1-sigma: 3894-3880BC, 3797-3709BC --- 2-sigma: 3941-3864BC, 3810-3700BC, 3680-3655BC --- *** --- *'''** ***** --- *** *'''''''***''''''* --- ******'''* '''''''''''''''' ** --- *''''''''''** **''''''''''''''''* **''* --- ***''''''''''''''********'''''''''''''''''''''***''''''** --- -3960 \<~~~~~~~~~|~~~~~~~~~~~~~~~~|~~~~~~~~~~~~~~~|~~~~~~~~~~\> -3640 --- @ --- -renderCalDatesPretty :: [(UncalC14, CalPDF, CalC14)] -> String -renderCalDatesPretty xs = - intercalate "\n" $ map renderCalDatePretty xs - -renderCalDatePretty :: (UncalC14, CalPDF, CalC14) -> String -renderCalDatePretty (uncalC14, calPDF, calC14) = - intercalate "\n" [ - renderUncalC14 uncalC14 - , renderCalC14 calC14 - , renderCLIPlotCalPDF 5 50 calPDF - ] - -- CalibrationMethod readCalibrationMethod :: String -> Either String CalibrationMethod readCalibrationMethod s = @@ -72,6 +44,86 @@ _ <- P.string "MatrixMult" return MatrixMultiplication +-- | Combine 'CalExpr', 'CalPDF' and 'CalC14' to render pretty command line output +-- like this: +-- +-- @ +-- DATE: (5000±30BP + 5100±100BP) +-- Calibrated: 4150BC \>\> 3941BC \> 3814BC \< 3660BC \<\< 3651BC +-- 1-sigma: 3941-3864BC, 3810-3707BC, 3667-3660BC +-- 2-sigma: 4150-4148BC, 4048-3651BC +-- ▁ +-- ▒▁ ▁▁ +-- ▁▁▁ ▁▒▒▁▒▒ +-- ▁▁▒▒▒ ▒▒▒▒▒▒ +-- ▁▁▒▒▒▒▒▁▁▁▁▒▒▒▒▒▒▁ ▁ +-- ▁▁▁▁▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▁▒▁ +-- ▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▁▁▁▁▁▁▁▁ +-- -4330 ┄──┬─────┬─────┬─────┬──────┬─────┬─────┬─────┬─────┄ -3530 +-- \> \> \^ \< +-- ────── ──────── ── +-- ─ ────────────────────────── +-- @ +-- +renderCalDatePretty :: + Bool -- ^ Should the CLI plot be restricted to (boring) ASCII symbols? + -> (CalExpr, CalPDF, CalC14) + -> String +renderCalDatePretty ascii (calExpr, calPDF, calC14) = + "DATE: " ++ intercalate "\n" [ + renderCalExpr calExpr + , renderCalC14 calC14 + , renderCLIPlotCalPDF ascii 6 50 calPDF calC14 + ] + +renderCalExpr :: CalExpr -> String +renderCalExpr (UnCalDate a) = renderUncalC14 a +renderCalExpr (CalDate (CalPDF name _ _)) = name +renderCalExpr (SumCal a b) = "(" ++ renderCalExpr a ++ " + " ++ renderCalExpr b ++ ")" +renderCalExpr (ProductCal a b) = "(" ++ renderCalExpr a ++ " * " ++ renderCalExpr b ++ ")" + +-- https://gist.github.com/abhin4v/017a36477204a1d57745 +spaceChar :: Char -> P.Parser Char +spaceChar c = P.between P.spaces P.spaces (P.char c) +--spaceChar = P.char + +add :: P.Parser CalExpr +add = SumCal <$> term <*> (spaceChar '+' *> expr) + +mul :: P.Parser CalExpr +mul = ProductCal <$> factor <*> (spaceChar '*' *> term) + +parens :: P.Parser CalExpr +parens = P.between (spaceChar '(') (spaceChar ')') expr + +factor :: P.Parser CalExpr +factor = parens P.<|> (UnCalDate <$> parseUncalC14) + +term :: P.Parser CalExpr +term = P.try mul P.<|> factor + +expr :: P.Parser CalExpr +expr = P.try add P.<|> term -- <* P.eof + +readCalExpr :: String -> Either String [CalExpr] +readCalExpr s = + case P.runParser parseCalExprSepBySemicolon () "" s of + Left err -> Left $ renderCurrycarbonException $ CurrycarbonCLIParsingException $ show err + Right x -> Right x + where + parseCalExprSepBySemicolon :: P.Parser [CalExpr] + parseCalExprSepBySemicolon = P.sepBy expr (P.char ';' <* P.spaces) <* P.eof + +readCalExprFromFile :: FilePath -> IO [CalExpr] +readCalExprFromFile uncalFile = do + s <- readFile uncalFile + case P.runParser parseCalExprSepByNewline () "" s of + Left err -> throwIO $ CurrycarbonCLIParsingException $ show err + Right x -> return x + where + parseCalExprSepByNewline :: P.Parser [CalExpr] + parseCalExprSepByNewline = P.endBy expr (P.newline <* P.spaces) <* P.eof + -- CalC14 -- | Write 'CalC14's to the file system. The output file is a long .csv file with the following structure: -- @@ -106,7 +158,7 @@ appendFile path $ "\n" ++ renderCalC14ForFile calC14 renderCalC14ForFile :: CalC14 -> String -renderCalC14ForFile (CalC14 name hdrs68 hdrs95) = +renderCalC14ForFile (CalC14 name _ hdrs68 hdrs95) = intercalate "\n" $ map renderRow $ zip3 (repeat name) (repeat "1") (renderHDRsForFile hdrs68) ++ @@ -121,22 +173,38 @@ ++ intercalate "\n" (map renderCalC14 xs) renderCalC14 :: CalC14 -> String -renderCalC14 (CalC14 _ hdrs68 hdrs95) = - "1-sigma: " ++ renderHDRs (reverse hdrs68) ++ "\n" - ++ "2-sigma: " ++ renderHDRs (reverse hdrs95) +renderCalC14 (CalC14 _ rangeSummary hdrs68 hdrs95) = + "Calibrated: " ++ renderCalRangeSummary rangeSummary ++ "\n" + ++ "1-sigma: " ++ renderHDRs hdrs68 ++ "\n" + ++ "2-sigma: " ++ renderHDRs hdrs95 +renderCalRangeSummary :: CalRangeSummary -> String +renderCalRangeSummary s = + renderYearBCAD (_calRangeStartTwoSigma s) ++ " >> " + ++ renderYearBCAD (_calRangeStartOneSigma s) ++ " > " + ++ renderYearBCAD (_calRangeMedian s) ++ " < " + ++ renderYearBCAD (_calRangeStopOneSigma s) ++ " << " + ++ renderYearBCAD (_calRangeStopTwoSigma s) + +-- BCAD +renderYearBCAD :: YearBCAD -> String +renderYearBCAD x + | x < 0 = show (-x) ++ "BC" + | x >= 0 = show x ++ "AD" + | otherwise = error $ "This should never happen: " ++ show x + -- HDR renderHDRsForFile :: [HDR] -> [(String, String)] renderHDRsForFile = map renderHDRForFile renderHDRForFile :: HDR -> (String, String) -renderHDRForFile (HDR start stop) = (show stop, show start) +renderHDRForFile (HDR start stop) = (show start, show stop) renderHDRs :: [HDR] -> String renderHDRs xs = intercalate ", " (map renderHDR xs) renderHDR :: HDR -> String -renderHDR (HDR stop start) +renderHDR (HDR start stop) | start < 0 && stop <= 0 = show (-start) ++ "-" ++ show (-stop) ++ "BC" | start < 0 && stop > 0 = show (-start) ++ "BC-" ++ show stop ++ "AD" | start >= 0 && stop >= 0 = show start ++ "-" ++ show stop ++ "AD" @@ -198,19 +266,27 @@ where makeRow (x,y) = show name ++ "," ++ show x ++ "," ++ show y ++ "\n" -renderCLIPlotCalPDF :: Int -> Int -> CalPDF -> String -renderCLIPlotCalPDF rows cols (CalPDF _ cals dens) = - let binWidth = quot (VU.length dens) cols - -- last bin will often be shorter, which renders the whole plot - -- slightly incorrect for the last column - binDens = meanBinDens (fromIntegral rows) binWidth dens - plotRows = map (replicate 8 ' ' ++) $ map (\x -> map (getSymbol x) binDens) $ reverse [0..rows] - xAxis = constructXAxis (VU.head cals) (VU.last cals) (length binDens) binWidth +data PlotSymbol = HistFill | HistTop | AxisEnd | AxisLine | AxisTick | HDRLine + +renderCLIPlotCalPDF :: Bool -> Int -> Int -> CalPDF -> CalC14 -> String +renderCLIPlotCalPDF ascii rows cols (CalPDF _ cals dens) c14 = + let startYear = VU.head cals + stopYear = VU.last cals + yearsPerCol = case quot (VU.length cals) cols of + 0 -> 1 -- relevant for very short PDFs + 1 -> 2 + q -> q + -- last bin will often be shorter, which renders the whole plot slightly incorrect for the last column + meanDensPerCol = calculateMeanDens yearsPerCol dens + effectiveCols = length meanDensPerCol + plotRows = map (replicate 8 ' ' ++) $ map (\x -> map (getHistSymbol x) meanDensPerCol) $ reverse [0..rows] + xAxis = constructXAxis startYear stopYear effectiveCols yearsPerCol in intercalate "\n" plotRows ++ "\n" ++ xAxis where - meanBinDens :: Float -> Int -> VU.Vector Float -> [Int] - meanBinDens scaling binWidth dens_ = - let meanDens = map (\x -> sum x / fromIntegral (length x)) $ splitEvery binWidth $ VU.toList dens_ + calculateMeanDens :: Int -> VU.Vector Float -> [Int] + calculateMeanDens yearsPerCol dens_ = + let scaling = fromIntegral rows + meanDens = map (\x -> sum x / fromIntegral (length x)) $ splitEvery yearsPerCol $ VU.toList dens_ maxDens = maximum meanDens in map (\x -> round $ (x / maxDens) * scaling) meanDens splitEvery :: Int -> [a] -> [[a]] -- https://stackoverflow.com/a/8681226/3216883 @@ -219,30 +295,73 @@ where (first,rest) = splitAt n list padString :: Int -> String -> String padString l x = replicate (l - length x) ' ' ++ x - getSymbol :: Int -> Int -> Char - getSymbol x y - | x == y = '*' - | x < y = '\'' + getSymbol :: Bool -> PlotSymbol -> Char + getSymbol True HistFill = '*' + getSymbol False HistFill = '▒' + getSymbol True HistTop = '_' + getSymbol False HistTop = '▁' + getSymbol True AxisEnd = '+' + getSymbol False AxisEnd = '┄' + getSymbol True AxisLine = '-' + getSymbol False AxisLine = '─' + getSymbol True AxisTick = '|' + getSymbol False AxisTick = '┬' + getSymbol True HDRLine = '-' + getSymbol False HDRLine = '─' + getHistSymbol :: Int -> Int -> Char + getHistSymbol x y + | x == y = getSymbol ascii HistTop + | x < y = getSymbol ascii HistFill | otherwise = ' ' constructXAxis :: Int -> Int -> Int -> Int -> String - constructXAxis start stop l binWidth = - let startS = padString 6 (show $ roundTo10 start) - stopS = show (roundTo10 stop) - tickFreq = if abs (start - stop) < 1500 then 100 else 1000 - axis = zipWith (axisSymbol binWidth tickFreq) [0 .. (l - 1)] [1 .. l] - in startS ++ " <" ++ axis ++ "> " ++ stopS - where - axisSymbol axisL tickFreq a b = if hasTick tickFreq (start + axisL * a + 1) (start + axisL * b) then '|' else '~' - hasTick tickFreq a b = any (\x -> rem (abs x) tickFreq == 0) [a..b] - roundTo10 :: Int -> Int - roundTo10 x = - let (dec,rest) = quotRem (abs x) 10 - roundedDec = if rest >= 5 then dec + 1 else dec - in roundedDec * 10 * signum x + constructXAxis startYear stopYear effCols yearsPerCol = + let startS = padString 6 (show $ roundTo10 startYear) + stopS = show (roundTo10 stopYear) + tickFreq = if abs (startYear - stopYear) < 1500 then 100 else 1000 + colStartYears = map (\a -> startYear + yearsPerCol * a) [0..(effCols - 1)] + colStopYears = map (\b -> startYear + yearsPerCol * b - 1) [1..effCols] + axis = zipWith (getAxisSymbol tickFreq) colStartYears colStopYears + simpleRange = zipWith (getRangeSymbol (_calC14RangeSummary c14)) colStartYears colStopYears + hdrOne = zipWith (getHDRSymbol (_calC14HDROneSigma c14)) colStartYears colStopYears + hdrTwo = zipWith (getHDRSymbol (_calC14HDRTwoSigma c14)) colStartYears colStopYears + in startS ++ (" " ++ [getSymbol ascii AxisEnd]) ++ axis ++ ([getSymbol ascii AxisEnd] ++ " ") ++ stopS ++ "\n" ++ + replicate 8 ' ' ++ simpleRange ++ "\n" ++ + replicate 8 ' ' ++ hdrOne ++ "\n" ++ + replicate 8 ' ' ++ hdrTwo + where + roundTo10 :: Int -> Int + roundTo10 x = + let (dec,rest) = quotRem (abs x) 10 + roundedDec = if rest >= 5 then dec + 1 else dec + in roundedDec * 10 * signum x + getAxisSymbol :: Int -> Int -> Int -> Char + getAxisSymbol tickFreq colStartYear colStopYear + | any (\x -> rem x tickFreq == 0) [colStartYear..colStopYear] = getSymbol ascii AxisTick + | otherwise = getSymbol ascii AxisLine + getRangeSymbol :: CalRangeSummary -> Int -> Int -> Char + getRangeSymbol range colStartYear colStopYear + | colStartYear <= _calRangeMedian range && colStopYear >= _calRangeMedian range = '^' + | colStartYear <= _calRangeStartOneSigma range && colStopYear >= _calRangeStartOneSigma range = '>' + | colStartYear <= _calRangeStopOneSigma range && colStopYear >= _calRangeStopOneSigma range = '<' + | colStartYear <= _calRangeStartTwoSigma range && colStopYear >= _calRangeStartTwoSigma range = '>' + | colStartYear <= _calRangeStopTwoSigma range && colStopYear >= _calRangeStopTwoSigma range = '<' + | otherwise = ' ' + getHDRSymbol :: [HDR] -> Int -> Int -> Char + getHDRSymbol hdr colStartYear colStopYear + | any (doesOverlap colStartYear colStopYear) hdr = getSymbol ascii HDRLine + | otherwise = ' ' + where + doesOverlap :: Int -> Int -> HDR -> Bool + doesOverlap a b h = + let ha = _hdrstart h; hb = _hdrstop h + in (a >= ha && a <= hb) || (b >= ha && b <= hb) || (a <= ha && b >= hb) -- UncalC14 +renderUncalC14WithoutName :: UncalC14 -> String +renderUncalC14WithoutName (UncalC14 _ bp sigma) = show bp ++ "±" ++ show sigma ++ "BP" + renderUncalC14 :: UncalC14 -> String -renderUncalC14 (UncalC14 name bp sigma) = "Sample: " ++ name ++ " ~> [" ++ show bp ++ "±" ++ show sigma ++ "BP]" +renderUncalC14 (UncalC14 name bp sigma) = name ++ ":" ++ show bp ++ "±" ++ show sigma ++ "BP" -- | Read uncalibrated radiocarbon dates from a file. The file should feature one radiocarbon date -- per line in the form "\,\,\", where diff --git a/src/Currycarbon/SumCalibration.hs b/src/Currycarbon/SumCalibration.hs new file mode 100644 --- /dev/null +++ b/src/Currycarbon/SumCalibration.hs @@ -0,0 +1,71 @@ +{-# LANGUAGE FlexibleInstances #-} + +module Currycarbon.SumCalibration where + +import Currycarbon.Types +import Currycarbon.Utils +import Currycarbon.Calibration.Calibration +import Currycarbon.Calibration.Utils + +import Data.Foldable (foldl') +import qualified Data.Vector.Unboxed as VU +import Data.List (sortBy, groupBy) +import Data.Ord (comparing) + +-- | Evaluate a dating expression by calibrating the individual dates and forming the respective +-- sums and products of post-calibration density distributions +evalCalExpr :: CalibrateDatesConf -> CalCurveBP -> CalExpr -> Either CurrycarbonException CalPDF +evalCalExpr conf curve calExpr = mapEither id normalizeCalPDF $ evalE calExpr + where + evalE :: CalExpr -> Either CurrycarbonException CalPDF + evalE (UnCalDate a) = calibrateDate conf curve a + evalE (CalDate a) = Right a + evalE (SumCal a b) = eitherCombinePDFs (+) 0 (evalE a) (evalE b) + evalE (ProductCal a b) = mapEither id normalizeCalPDF $ eitherCombinePDFs (*) 1 + (mapEither id normalizeCalPDF $ evalE a) (mapEither id normalizeCalPDF $ evalE b) -- product needs extra normalization + -- https://hackage.haskell.org/package/either-5.0.2/docs/Data-Either-Combinators.html + mapEither :: (a -> c) -> (b -> d) -> Either a b -> Either c d + mapEither f _ (Left x) = Left (f x) + mapEither _ f (Right x) = Right (f x) + +eitherCombinePDFs :: + (Float -> Float -> Float) -> Float -> + Either CurrycarbonException CalPDF -> + Either CurrycarbonException CalPDF -> + Either CurrycarbonException CalPDF +eitherCombinePDFs _ _ (Left e) _ = Left e +eitherCombinePDFs _ _ _ (Left e) = Left e +eitherCombinePDFs f initVal (Right a) (Right b) = Right $ combinePDFs f initVal a b + +-- | Add two probabilty densities +addPDFs :: CalPDF -> CalPDF -> CalPDF +addPDFs = combinePDFs (+) 0 + +-- | Multiply two probabilty densities +multiplyPDFs :: CalPDF -> CalPDF -> CalPDF +multiplyPDFs = combinePDFs (*) 1 + +-- Combine probability densities +combinePDFs :: (Float -> Float -> Float) -> Float -> CalPDF -> CalPDF -> CalPDF +combinePDFs f initVal (CalPDF name1 cals1 dens1) (CalPDF name2 cals2 dens2) = + let minC1 = VU.minimum cals1 + minC2 = VU.minimum cals2 + maxC1 = VU.maximum cals1 + maxC2 = VU.maximum cals2 + emptyC1 = getMiss minC1 maxC1 minC2 maxC2 + emptyC2 = getMiss minC2 maxC2 minC1 maxC1 + c1 = VU.toList (VU.zip cals1 dens1) ++ zip emptyC1 (repeat (0 :: Float)) + c2 = VU.toList (VU.zip cals2 dens2) ++ zip emptyC2 (repeat (0 :: Float)) + pdfSorted = sortBy (comparing fst) (c1 ++ c2) + pdfGrouped = groupBy (\a b -> fst a == fst b) pdfSorted + pdfRes = map foldYearGroup pdfGrouped + in CalPDF (name1 ++ ":" ++ name2) (VU.fromList $ map fst pdfRes) (VU.fromList $ map snd pdfRes) + where + getMiss :: YearBCAD -> YearBCAD -> YearBCAD -> YearBCAD -> [YearBCAD] + getMiss a1 a2 b1 b2 + | a1 < b1 && a2 > b2 = [a1..b1] ++ [b2..a2] + | a1 < b1 && a2 <= b2 = [a1..b1] + | a1 >= b1 && a2 > b2 = [b2..a2] + | otherwise = [] + foldYearGroup :: [(YearBCAD, Float)] -> (YearBCAD, Float) + foldYearGroup oneYear = (fst $ head oneYear, foldl' f initVal $ map snd oneYear) diff --git a/src/Currycarbon/Types.hs b/src/Currycarbon/Types.hs --- a/src/Currycarbon/Types.hs +++ b/src/Currycarbon/Types.hs @@ -108,17 +108,43 @@ , _calPDFCals :: VU.Vector YearBCAD -- | Probability densities for each year in '_calPDFCals' , _calPDFDens :: VU.Vector Float - } deriving Show + } deriving (Show, Eq) +-- | A data type to represent an expression for sum- or product calibration +data CalExpr = + UnCalDate UncalC14 + | CalDate CalPDF + | SumCal CalExpr CalExpr + | ProductCal CalExpr CalExpr + deriving Show +-- http://www.cse.chalmers.se/edu/year/2018/course/TDA452/lectures/RecursiveDataTypes.html + -- | A data type to represent a human readable summary of a calibrated radiocarbon date data CalC14 = CalC14 { -- | Identifier, e.g. a lab number _calC14id :: String + -- | Summary of the range of the calibrated date + , _calC14RangeSummary :: CalRangeSummary -- | One-sigma high density regions , _calC14HDROneSigma :: [HDR] -- | Two-sigma high density regions , _calC14HDRTwoSigma :: [HDR] } deriving Show + +-- | A data type to store a simple range summary of the calibrated date, including +-- the median age +data CalRangeSummary = CalRangeSummary { + -- | Start of the two-sigma age range + _calRangeStartTwoSigma :: YearBCAD + -- | Start of the one-sigma age range + , _calRangeStartOneSigma :: YearBCAD + -- | Median age + , _calRangeMedian :: YearBCAD + -- | End of the one-sigma age range + , _calRangeStopOneSigma :: YearBCAD + -- | End of the two-sigma age range + , _calRangeStopTwoSigma :: YearBCAD +} deriving Show -- | A data type to represent a high density region of a probability distribution. -- A high density region is here defined as an age range, within which the respective diff --git a/test/Spec.hs b/test/Spec.hs deleted file mode 100644 --- a/test/Spec.hs +++ /dev/null @@ -1,11 +0,0 @@ -import Test.DocTest - -import Currycarbon.CalCurves.Intcal20 - -main :: IO () -main = doctest [ - "-isrc" - , "src/Currycarbon/Calibration/Bchron.hs" - , "src/Currycarbon/Calibration/Calibration.hs" - , "src/Currycarbon/Calibration/Utils.hs" - ]