currycarbon 0.1.1.0 → 0.2.1.0

raw patch · 13 files changed

+509/−224 lines, 13 filesdep −doctestdep ~basePVP ok

version bump matches the API change (PVP)

Dependencies removed: doctest

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• CHANGELOG.md +4/−0
• README.md +42/−13
• currycarbon.cabal +2/−17
• src-executables/Main-currycarbon.hs +28/−19
• src/Currycarbon.hs +24/−1
• src/Currycarbon/CLI/RunCalibrate.hs +75/−58
• src/Currycarbon/Calibration/Bchron.hs +0/−9
• src/Currycarbon/Calibration/Calibration.hs +46/−19
• src/Currycarbon/Calibration/Utils.hs +7/−12
• src/Currycarbon/Parsers.hs +183/−64
• src/Currycarbon/SumCalibration.hs +71/−0
• src/Currycarbon/Types.hs +27/−1
• test/Spec.hs +0/−11

@@ -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

@@ -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 "<sample name>,<mean age BP>,<one sigma-                           standard deviation>;..." where <sample name> 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 <sample name> 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

@@ -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

@@ -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 \"<sample name>,<mean age BP>,<one sigma standard deviation>;...\" \-            \where <sample name> is optional. \-            \So for example \"S1,4000,50;3000,25;S3,1000,20\"."+            \the form \"<sample name>,<mean age BP>,<one sigma standard deviation>\" \+            \where <sample name> 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)

@@ -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. -}

@@ -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)

@@ -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

@@ -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

@@ -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) =

@@ -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 "\<sample name\>,\<mean age BP\>,\<one sigma standard deviation\>", where 

@@ -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)

@@ -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 

@@ -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"-    ]