currycarbon 0.2.1.2 → 0.3.0.0

raw patch · 16 files changed

+1213/−351 lines, 16 filesdep +MonadRandomdep +hspecdep +processdep ~basedep ~filepathPVP ok

version bump matches the API change (PVP)

Dependencies added: MonadRandom, hspec, process, random

Dependency ranges changed: base, filepath

API changes (from Hackage documentation)

Files

• CHANGELOG.md +66/−49
• README.md +97/−56
• currycarbon.cabal +24/−3
• src-executables/Main-currycarbon.hs +127/−37
• src/Currycarbon.hs +19/−1
• src/Currycarbon/CLI/RunCalibrate.hs +128/−55
• src/Currycarbon/Calibration/Calibration.hs +59/−8
• src/Currycarbon/Calibration/Utils.hs +4/−1
• src/Currycarbon/ParserHelpers.hs +168/−0
• src/Currycarbon/Parsers.hs +267/−127
• src/Currycarbon/SumCalibration.hs +31/−7
• src/Currycarbon/Types.hs +26/−2
• src/Currycarbon/Utils.hs +17/−5
• test/GoldenSpec.hs +62/−0
• test/ParserSpec.hs +117/−0
• test/Spec.hs +1/−0

@@ -1,49 +1,66 @@-- V 0.2.1.2: Maintenance: Switched to a newer compiler/resolver version, lifted some dependency restrictions, ran stylish-haskell on the entire codebase, updated the github actions, deprecated the haddock documentation for the dev version on GitHub-- V 0.2.1.1: Lifted some restrictions regarding the upper version bounds of dependencies-- 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-- V 0.24.3: Multiple changes in .cabal to make cabal check happy-- V 0.24.2: Found and fixed another severe bug in renderCalCurve-- V 0.24.1: Fixed a serious bug in renderCalCurveMatrix-- V 0.24.0: Introduced more precise data types to distinguish years BP and years BC/AD-- V 0.23.1: Small changes to the instances of some general types-- V 0.23.0: Renamed multiple functions to make the naming of operations for parsing, reading, from-file reading, rendering and writing consistent across data types-- V 0.22.0: Changed the interface of the important calibrateDates function with a new config data type CalibrateDatesConf-- V 0.21.3: Refactored the calibration curve interpolation-- V 0.21.2: Introduced doctest and added some tiny examples/tests to try it out-- V 0.21.1: Split up the calibration module for better readability-- V 0.21.0: Added a neat CLI density plot for calibrated dates-- V 0.20.2: Some performance improvements for the calibration of large numbers of dates-- V 0.20.1: Better (parsing) error handling-- V 0.20.0: Added an option --allowOutside to allow for calibrations to run outside the range of the calibration curve-- V 0.19.0: Added functionality to filter out dates outside of the range of the calibration curve and report an error in this case-- V 0.18.0: Implemented calibration with a StudentT distribution to mimic Bchron and established that as the new default. Reimplemented the --method option of the CLI tool to reflect that change-- V 0.17.0: Changed argument order in CalCurve data type to adjust to the order in .14C files-- V 0.16.0: Refactoring in the library to simplify and clarify the interface-- V 0.15.0: Added another calibration algorithm (following the implementation by Andrew Parnell in Bchron) and a method switch for the CLI-- V 0.14.0: Introduced strictness, which brought a significant increase in performance. See the discussion here: https://old.reddit.com/r/haskell/comments/picjy6/how_could_i_improve_the_performance_of_my/-- V 0.13.0: Major rewrite with the vector library - includes multiple bugfixes, but is surprisingly slow-- V 0.12.0: Renamed some core functions-- V 0.11.0: Made calibration curve interpolation optional and turned it off by default-- V 0.10.0: Simplified CLI interface by dropping the "calibrate" subcommand (currycarbon is sufficient now) and by repurposing -q from --quickOut to --quiet-- V 0.9.0: Made --hdrFile output a lot more machine-readable-- V 0.8.0: Added option --calibrationCurveFile to calibrate with different calibration curves-- V 0.7.2: More documentation, small changes in code layout and renamed CLI module that provides runCalibrate-- V 0.7.1: Added type documentation with haddock and replaced the existing types with record types-- V 0.7.0: Changed the date input interface once more-- V 0.6.0: Changed the date input interface, because parenthesis can be part of valid lab numbers-- V 0.5.2: Fixed parallel evalutation (deepseq forced memory-intensive, non-lazy behaviour)-- V 0.5.1: Added github release action (copied from poseidon-hs)-- V 0.5.0: Added file input for dates to calibrate-- V 0.4.0: Made output calibrated dates negative numbers for BC and positive for AD - and adjusted HDR printing accordingly-- V 0.3.2: Some optimisation-- V 0.3.1: Added automatic filling of unknown sample names-- V 0.3.0: Simplified interface-- V 0.2.1: Removed ascii plot functionality-- V 0.2.0: Added parallel processing for the main calibration operation-- V 0.1.0: First basically working version+- V 0.3.0.0: Major update with multiple breaking changes and new features:+	- Added a new mechanism to draw random age samples from a CalPDF (`sampleAgesFromCalPDF :: AgeSamplingConf -> CalPDF -> RandomAgeSample`). This is available from the command line with the options `samplesFile`, `--seed`, and `-n`/`--nrSamples`.+	- Added a new concept to the `CalExpr` data type: Age ranges with uniform probability for each year in the range (`TimeWindowBP` and `TimeWindowBCAD`).+	- Reworked the encoding and evaluation mechanism for calibration expressions:+		- Introduced the `NamedCalExpr` as a wrapper around `CalExpr` with an identifier, and then adjusted the ID generation for `CalPDF`s to prioritize this identifier.+		- Reworked the CLI DSL to support a standardized configuration language syntax implemented in a new module `ParserHelpers.hs`. This introduces a set of flexible functions (`calExpr()`, `uncalC14()`, `rangeBP()`, `rangeBCAD()`, `sum()` and `product()`) which generally complement the previously available syntax and operators. The old syntax is mostly preserved as syntactic sugar for the new, more standardized syntax. Unfortunately this is not entirely seamless: The change breaks some expressions that were valid before (e.g. `"3000,30 + 3020,50"`). They now require additional parentheses to pass (so e.g. `"(3000,30) + (3020,50)"`).+		- Added some unit tests to cover the increasingly complex DSL.+	- Changed the output files from .csv to .tsv and to a more meaningful and consistent set of column names.+	- Slightly adjusted the rendering of the pretty, human-focussed command line output.+	- Updated and improved the command line documentation.+	- Renamed some CLI arguments:+		- `--calibrationCurveFile` -> `--calCurveFile`+		- `--calCurveSegmentFile` -> `--calCurveSegFile`+		- `--calCurveMatrixFile` -> `--calCurveMatFile`+	- Changed the CLI behaviour with  `--calCurveSegFile` and `--calCurveMatFile`: currycarbon now fails with these options if the first sample is not a single, uncalibrated radiocarbon date (so `uncalC14()`).+	- Added a simple golden test system with some basic calls to the `currycarbon` CLI tool.+	- Switched to a new GHC version (v9.4.7) and stackage resolver version (lts-21.17).+- V 0.2.1.2: Maintenance: Switched to a newer compiler/resolver version, lifted some dependency restrictions, ran stylish-haskell on the entire codebase, updated the github actions, deprecated the haddock documentation for the dev version on GitHub.+- V 0.2.1.1: Lifted some restrictions regarding the upper version bounds of dependencies.+- 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. `HDR`s 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.+- V 0.24.3: Multiple changes in .cabal to make cabal check happy.+- V 0.24.2: Found and fixed another severe bug in `renderCalCurve`.+- V 0.24.1: Fixed a serious bug in `renderCalCurveMatrix`.+- V 0.24.0: Introduced more precise data types to distinguish years BP and years BC/AD.+- V 0.23.1: Small changes to the instances of some general types.+- V 0.23.0: Renamed multiple functions to make the naming of operations for parsing, reading, from-file reading, rendering and writing consistent across data types.+- V 0.22.0: Changed the interface of the important `calibrateDates` function with a new config data type `CalibrateDatesConf`.+- V 0.21.3: Refactored the calibration curve interpolation.+- V 0.21.2: Introduced doctest and added some tiny examples/tests to try it out.+- V 0.21.1: Split up the calibration module for better readability.+- V 0.21.0: Added a neat CLI density plot for calibrated dates.+- V 0.20.2: Some performance improvements for the calibration of large numbers of dates.+- V 0.20.1: Better (parsing) error handling.+- V 0.20.0: Added an option `--allowOutside` to allow for calibrations to run outside the range of the calibration curve.+- V 0.19.0: Added functionality to filter out dates outside of the range of the calibration curve and report an error in this case.+- V 0.18.0: Implemented calibration with a StudentT distribution to mimic Bchron and established that as the new default. Reimplemented the `--method` option of the CLI tool to reflect that change.+- V 0.17.0: Changed argument order in the `CalCurve` data type to adjust to the order in `.14C` files.+- V 0.16.0: Refactoring in the library to simplify and clarify the interface.+- V 0.15.0: Added another calibration algorithm (following the implementation by Andrew Parnell in Bchron) and a method switch for the CLI.+- V 0.14.0: Introduced strictness, which brought a significant increase in performance. See the discussion here: https://old.reddit.com/r/haskell/comments/picjy6/how_could_i_improve_the_performance_of_my/.+- V 0.13.0: Major rewrite with the vector library - includes multiple bugfixes, but is surprisingly slow.+- V 0.12.0: Renamed some core functions.+- V 0.11.0: Made calibration curve interpolation optional and turned it off by default.+- V 0.10.0: Simplified CLI interface by dropping the `calibrate` subcommand (`currycarbon` is sufficient now) and by repurposing `-q` from `--quickOut` to `--quiet`.+- V 0.9.0: Made `--hdrFile` output a lot more machine-readable.+- V 0.8.0: Added option `--calibrationCurveFile` to calibrate with different calibration curves.+- V 0.7.2: More documentation, small changes in code layout and renamed CLI module that provides `runCalibrate`.+- V 0.7.1: Added type documentation with haddock and replaced the existing types with record types.+- V 0.7.0: Changed the date input interface once more.+- V 0.6.0: Changed the date input interface, because parenthesis can be part of valid lab numbers.+- V 0.5.2: Fixed parallel evaluation (deepseq forced memory-intensive, non-lazy behaviour).+- V 0.5.1: Added github release action (copied from poseidon-hs).+- V 0.5.0: Added file input for dates to calibrate.+- V 0.4.0: Made output calibrated dates negative numbers for BC and positive for AD - and adjusted HDR printing accordingly.+- V 0.3.2: Some optimisation.+- V 0.3.1: Added automatic filling of unknown sample names.+- V 0.3.0: Simplified interface.+- V 0.2.1: Removed ascii plot functionality.+- V 0.2.0: Added parallel processing for the main calibration operation.+- V 0.1.0: First basically working version.

@@ -1,4 +1,5 @@ [![GitHub Workflow Status](https://github.com/nevrome/currycarbon/actions/workflows/normalCheck.yml/badge.svg)](https://github.com/nevrome/currycarbon/actions/workflows/normalCheck.yml)+[![Coverage Status](https://img.shields.io/codecov/c/github/nevrome/currycarbon/master.svg)](https://codecov.io/github/nevrome/currycarbon?branch=master) [![GitHub release (latest by date including pre-releases)](https://img.shields.io/github/v/release/nevrome/currycarbon?include_prereleases)  ![GitHub all releases](https://img.shields.io/github/downloads/nevrome/currycarbon/total)](https://github.com/nevrome/currycarbon/releases) @@ -28,86 +29,126 @@ ```  ```-currycarbon v0.2.1.0 (UTF-8)+currycarbon v0.3.0.0 (UTF-8) Method: Bchron {distribution = StudentTDist {ndf = 100.0}} Curve: IntCal20 Calibrating...-DATE: Sample1:4990±30BP+CalEXPR: [1] 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] -                   [--calCurveMatrixFile ARG]+Usage: currycarbon [--version] [CalEXPRs] [-i|--inputFile FILE]+                   [--calCurveFile FILE] [--method DSL] [--allowOutside]+                   [--noInterpolation] [-q|--quiet] [--densityFile FILE]+                   [--hdrFile FILE]+                   [[--seed INT] (-n|--nrSamples INT) --samplesFile FILE]+                   [--calCurveSegFile FILE] [--calCurveMatFile FILE]+   Intercept calibration of radiocarbon dates  Available options:   -h,--help                Show this help text   --version                Show version-  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-                           (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+  CalEXPRs                 ---+                           A string to specify "calibration expressions", so+                           small chronological models for individual events.+                           These can include uncalibrated radiocarbon ages,+                           uniform age ranges and operations to combine the+                           resulting age probability distribution as sums or+                           products.+                           The expression language includes the following+                           functions:++                           - calExpr(id = STRING, expr = EXPR)+                           - uncalC14(id = STRING, yearBP = INT, sigma = INT)+                           - rangeBP(id = STRING, start = INT, stop = INT)+                           - rangeBCAD(id = STRING, start = INT, stop = INT)+                           - sum(a = EXPR, b = EXPR)+                           - product(a = EXPR, b = EXPR)++                           The order of arguments is fixed, but the argument+                           names '<arg> =' can be left out. The 'id' arguments+                           are optional. Some functions can be shortened with+                           syntactic sugar:++                           - calExpr(STRING, EXPR) -> id: EXPR+                           - uncalC14(STRING, INT, INT) -> STRING,INT,INT+                           - sum(EXPR, EXPR) -> EXPR + EXPR+                           - product(EXPR, EXPR) -> EXPR * EXPR++                           Parentheses '()' can be used to specify the+                           evaluation order within an expression. Multiple+                           expressions can be chained, separated by ';'.++                           Examples:+                           1. Calibrate a single radiocarbon date with a mean+                           age BP and a one sigma standard deviation:+                           "3000,30" or "uncalC14(yearBP = 3000, sigma = 30)"+                           2. Calibrate two radiocarbon dates and sum them:+                           "(3000,30) + (3100,40)" or+                           "sum(uncalC14(3000,30), uncalC14(3100,40))"+                           3. Compile a complex, named expression:+                           "Ex3: ((3000,30) + (3100,40)) * rangeBP(3200,3000)"+                           ---+  -i,--inputFile FILE      A file with a list of calibration expressions.+                           Formatted just as CalEXPRs, but with a new line for+                           each input expression. CalEXPRs and --inputFile can+                           be combined and you can provide multiple instances of+                           --inputFile. Note that syntactic sugar allows to read+                           simple radiocarbon dates from a headless .csv file+                           with one sample per row: <sample name>,<mean age+                           BP>,<one sigma standard deviation>.+  --calCurveFile FILE      Path to an calibration curve file in '.14c' format.+                           The calibration curve will be read and used for                            calibration. If no file is provided, currycarbon will-                           use the intcal20 curve.-  --method ARG             The calibration algorithm that should be used:-                           "<Method>,<Distribution>,<NumberOfDegreesOfFreedom>".+                           use the 'intcal20' curve.+  --method DSL             The calibration algorithm that should be used:+                           '<Method>,<Distribution>,<NumberOfDegreesOfFreedom>'.                            The default setting is equivalent to                            "Bchron,StudentT,100" which copies the algorithm-                           implemented in the Bchron R package. Alternatively we-                           implemented "MatrixMult", which comes without further-                           arguments. For the Bchron algorithm with a normal-                           distribution ("Bchron,Normal") the degrees of freedom-                           argument is not relevant+                           implemented in the Bchron R package. For the Bchron+                           algorithm with a normal distribution+                           ("Bchron,Normal") the degrees of freedom argument is+                           not relevant+                           Alternatively we implemented "MatrixMult", which+                           comes without further arguments.   --allowOutside           Allow calibrations to run outside the range of the-                           calibration curve-  --noInterpolation        Don't interpolate the calibration curve+                           calibration curve.+  --noInterpolation        Do not interpolate the calibration curve.   -q,--quiet               Suppress the printing of calibration results to the-                           command line-  --densityFile ARG        Path to an output file which stores output densities-                           per sample and calender year-  --hdrFile ARG            Path to an output file which stores the high-                           probability density regions for each sample-  --calCurveSegmentFile ARG-                           Path to an output file which stores the relevant,+                           command line.+  --densityFile FILE       Path to an output file to store output densities per+                           CalEXPR and calender year.+  --hdrFile FILE           Path to an output file to store the high probability+                           density regions for each CalEXPR.+  --seed INT               Seed for the random number generator for age+                           sampling. The default causes currycarbon to fall back+                           to a random seed. (default: Nothing)+  -n,--nrSamples INT       Number of age samples to draw per CalEXPR.+  --samplesFile FILE       Path to an output file to store age samples for each+                           CalEXPR.+  --calCurveSegFile FILE   Path to an output file to store the relevant,                            interpolated calibration curve segment for the first-                           (!) input date in a long format. This option as well-                           as --calCurveMatrixFile are mostly meant for-                           debugging-  --calCurveMatrixFile ARG Path to an output file which stores the relevant,+                           (!) input date. This option as well as+                           --calCurveMatFile are meant for debugging.+  --calCurveMatFile FILE   Path to an output file which stores the relevant,                            interpolated calibration curve segment for the first-                           (!) input date in a wide matrix format+                           (!) input date in a wide matrix format. ```  ### For developers who want to edit the code

@@ -1,5 +1,5 @@ name:                currycarbon-version:             0.2.1.2+version:             0.3.0.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@@ -26,6 +26,7 @@       Currycarbon.Calibration.Bchron       Currycarbon.Calibration.Calibration       Currycarbon.CLI.RunCalibrate+      Currycarbon.ParserHelpers       Currycarbon.Parsers       Currycarbon.SumCalibration       Currycarbon.Types@@ -38,6 +39,8 @@       , parsec >= 3.1 && < 3.2       , vector >= 0.12 && < 0.14       , math-functions >= 0.3 && < 0.4+      , MonadRandom >= 0.6 && < 1+      , random > 1.2 && < 1.3     default-language:       Haskell2010 @@ -48,12 +51,30 @@       src-executables     build-depends:               currycarbon-      , base+      , base >= 4.14 && < 5       , optparse-applicative >= 0.16 && < 0.19-      , filepath+      , filepath >= 1.4 && < 1.5     other-modules:       Paths_currycarbon     default-language:       Haskell2010     ghc-options:        -threaded -with-rtsopts=-N -optP-Wno-nonportable-include-path++test-suite spec+    hs-source-dirs:+      test+    main-is:+      Spec.hs+    type:+      exitcode-stdio-1.0+    build-depends:+        base >= 4.9 && < 5+      , currycarbon+      , hspec >= 2 && < 3+      , process >= 1.6 && < 1.7+    other-modules:+      ParserSpec,+      GoldenSpec+    default-language:+      Haskell2010

@@ -10,6 +10,7 @@ import           Control.Exception            (catch) import           Data.Version                 (showVersion) import qualified Options.Applicative          as OP+import qualified Options.Applicative.Help     as OH import           System.Exit                  (exitFailure) import           System.IO                    (hGetEncoding, hPutStrLn, stderr,                                                stdout)@@ -56,8 +57,8 @@ optParser = CmdCalibrate <$> calibrateOptParser  calibrateOptParser :: OP.Parser CalibrateOptions-calibrateOptParser = CalibrateOptions <$> optParseCalExprString-                                      <*> optParseCalExprFromFile+calibrateOptParser = CalibrateOptions <$> optParseNamedCalExprString+                                      <*> optParseNamedCalExprFromFile                                       <*> optParseCalCurveFromFile                                       <*> optParseCalibrationMethod                                       <*> optParseAllowOutside@@ -66,6 +67,7 @@                                       <*> pure "unknown"                                       <*> optParseDensityFile                                       <*> optParseHDRFile+                                      <*> optParseAgeSamplingSettings                                       <*> optParseCalCurveSegmentFile                                       <*> optParseCalCurveMatrixFile @@ -75,97 +77,185 @@ -- -- These functions define and handle the CLI input arguments -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 (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\")"+optParseNamedCalExprString :: OP.Parser [NamedCalExpr]+optParseNamedCalExprString = concat <$> OP.many (OP.argument (OP.eitherReader readNamedCalExprs) (+    OP.metavar "CalEXPRs" <>+    OP.helpDoc ( Just (+        "---"+        <> OH.hardline+        <> s2d "A string to specify \"calibration expressions\", so small chronological \+           \models for individual events. These can include uncalibrated radiocarbon ages, \+           \uniform age ranges and operations to combine the resulting age probability \+           \distribution as sums or products."+        <> OH.hardline <>+            s2d "The expression language includes the following functions:"+        <> OH.hardline+        <> OH.hardline <> "- calExpr(id = STRING, expr = EXPR)"+        <> OH.hardline <> "- uncalC14(id = STRING, yearBP = INT, sigma = INT)"+        <> OH.hardline <> "- rangeBP(id = STRING, start = INT, stop = INT)"+        <> OH.hardline <> "- rangeBCAD(id = STRING, start = INT, stop = INT)"+        <> OH.hardline <> "- sum(a = EXPR, b = EXPR)"+        <> OH.hardline <> "- product(a = EXPR, b = EXPR)"+        <> OH.hardline+        <> OH.hardline+        <> s2d "The order of arguments is fixed, but the argument names '<arg> =' \+           \can be left out. The 'id' arguments are optional. \+           \Some functions can be shortened with syntactic sugar:"+        <> OH.hardline+        <> OH.hardline <> "- calExpr(STRING, EXPR) -> id: EXPR"+        <> OH.hardline <> "- uncalC14(STRING, INT, INT) -> STRING,INT,INT"+        <> OH.hardline <> "- sum(EXPR, EXPR) -> EXPR + EXPR"+        <> OH.hardline <> "- product(EXPR, EXPR) -> EXPR * EXPR"+        <> OH.hardline+        <> OH.hardline+        <> s2d "Parentheses '()' can be used to specify the evaluation order within \+           \an expression. Multiple expressions can be chained, separated by ';'."+        <> OH.hardline+        <> OH.hardline+        <> "Examples:"+        <> OH.hardline <> s2d "1. Calibrate a single radiocarbon date with a mean age BP \+                          \and a one sigma standard deviation:"+        <> OH.hardline <> "\"3000,30\" or \"uncalC14(yearBP = 3000, sigma = 30)\""+        <> OH.hardline <> s2d "2. Calibrate two radiocarbon dates and sum them:"+        <> OH.hardline <> "\"(3000,30) + (3100,40)\" or"+        <> OH.hardline <> "\"sum(uncalC14(3000,30), uncalC14(3100,40))\""+        <> OH.hardline <> s2d "3. Compile a complex, named expression:"+        <> OH.hardline <> "\"Ex3: ((3000,30) + (3100,40)) * rangeBP(3200,3000)\""+        <> OH.hardline+        <> "---"     ))+    )) -optParseCalExprFromFile :: OP.Parser [FilePath]-optParseCalExprFromFile = OP.many (OP.strOption (+s2d :: String -> OH.Doc+s2d str = OH.fillSep $ map OH.pretty $ words str++optParseNamedCalExprFromFile :: OP.Parser [FilePath]+optParseNamedCalExprFromFile = OP.many (OP.strOption (     OP.long "inputFile" <>     OP.short 'i' <>+    OP.metavar "FILE" <>     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"+            \Formatted just as CalEXPRs, but with a new line for each input expression. \+            \CalEXPRs and --inputFile can be combined and you can provide multiple \+            \instances of --inputFile. \+            \Note that syntactic sugar allows to read simple radiocarbon dates from \+            \a headless .csv file with one sample per row: \+            \<sample name>,<mean age BP>,<one sigma standard deviation>."     ))  optParseCalCurveFromFile :: OP.Parser (Maybe FilePath) optParseCalCurveFromFile = OP.option (Just <$> OP.str) (-    OP.long "calibrationCurveFile" <>-    OP.help "Path to an calibration curve file in .14c format. \+    OP.long "calCurveFile" <>+    OP.metavar "FILE" <>+    OP.help "Path to an calibration curve file in '.14c' format. \             \The calibration curve will be read and used for calibration. \-            \If no file is provided, currycarbon will use the intcal20 curve." <>+            \If no file is provided, currycarbon will use the 'intcal20' curve." <>     OP.value Nothing     )  optParseCalibrationMethod :: OP.Parser CalibrationMethod optParseCalibrationMethod = OP.option (OP.eitherReader readCalibrationMethod) (     OP.long "method" <>-    OP.help "The calibration algorithm that should be used: \-            \\"<Method>,<Distribution>,<NumberOfDegreesOfFreedom>\". \-            \The default setting is equivalent to \"Bchron,StudentT,100\" \+    OP.metavar "DSL" <>+    OP.helpDoc ( Just (+            s2d "The calibration algorithm that should be used: \+            \'<Method>,<Distribution>,<NumberOfDegreesOfFreedom>'. "+        <> OH.hardline <>+            s2d "The default setting is equivalent to \"Bchron,StudentT,100\" \             \which copies the algorithm implemented in the Bchron R package. \-            \Alternatively we implemented  \"MatrixMult\", which comes without further arguments. \             \For the Bchron algorithm with a normal distribution (\"Bchron,Normal\") \-            \the degrees of freedom argument is not relevant" <>+            \the degrees of freedom argument is not relevant"+        <> OH.hardline <>+            s2d "Alternatively we implemented  \"MatrixMult\", which comes without further \+            \arguments."+    )) <>     OP.value (Bchron $ StudentTDist 100)     )  optParseAllowOutside :: OP.Parser (Bool) optParseAllowOutside = OP.switch (     OP.long "allowOutside" <>-    OP.help "Allow calibrations to run outside the range of the calibration curve"+    OP.help "Allow calibrations to run outside the range of the calibration curve."     )  optParseDontInterpolateCalCurve :: OP.Parser (Bool) optParseDontInterpolateCalCurve = OP.switch (     OP.long "noInterpolation" <>-    OP.help "Don't interpolate the calibration curve"+    OP.help "Do not interpolate the calibration curve."     )  optParseQuiet :: OP.Parser (Bool) optParseQuiet = OP.switch (     OP.long "quiet" <>     OP.short 'q' <>-    OP.help "Suppress the printing of calibration results to the command line"+    OP.help "Suppress the printing of calibration results to the command line."     )  optParseDensityFile :: OP.Parser (Maybe FilePath) optParseDensityFile = OP.option (Just <$> OP.str) (     OP.long "densityFile" <>-    OP.help "Path to an output file which stores output densities per sample and calender year" <>+    OP.metavar "FILE" <>+    OP.help "Path to an output file to store output densities per CalEXPR and calender \+            \year." <>     OP.value Nothing     )  optParseHDRFile :: OP.Parser (Maybe FilePath) optParseHDRFile = OP.option (Just <$> OP.str) (     OP.long "hdrFile" <>-    OP.help "Path to an output file which stores the high probability density regions for each \-            \sample" <>+    OP.metavar "FILE" <>+    OP.help "Path to an output file to store the high probability density regions for each \+            \CalEXPR." <>     OP.value Nothing     ) +optParseAgeSamplingSettings :: OP.Parser (Maybe (Maybe Word, Word, FilePath))+optParseAgeSamplingSettings =+    OP.optional $ (,,) <$>+            optParseAgeSamplingConfSeed+        <*> optParseAgeSamplingConfNrOfSamples+        <*> optParseAgeSamplingFile++optParseAgeSamplingConfSeed :: OP.Parser (Maybe Word)+optParseAgeSamplingConfSeed = OP.option (Just <$> OP.auto) (+       OP.long  "seed"+    <> OP.metavar "INT"+    <> OP.help  "Seed for the random number generator for age sampling. \+                \The default causes currycarbon to fall back to a random seed."+    <> OP.value Nothing+    <> OP.showDefault+    )++optParseAgeSamplingConfNrOfSamples :: OP.Parser Word+optParseAgeSamplingConfNrOfSamples = OP.option OP.auto (+       OP.short 'n'+    <> OP.long "nrSamples"+    <> OP.metavar "INT"+    <> OP.help "Number of age samples to draw per CalEXPR."+    )++optParseAgeSamplingFile :: OP.Parser FilePath+optParseAgeSamplingFile = OP.strOption (+    OP.long "samplesFile" <>+    OP.metavar "FILE" <>+    OP.help "Path to an output file to store age samples for each CalEXPR."+    )+ optParseCalCurveSegmentFile :: OP.Parser (Maybe FilePath) optParseCalCurveSegmentFile = OP.option (Just <$> OP.str) (-    OP.long "calCurveSegmentFile" <>-    OP.help "Path to an output file which stores the relevant, interpolated calibration curve \-            \segment for the first (!) input date in a long format. \-            \This option as well as --calCurveMatrixFile are mostly meant for debugging" <>+    OP.long "calCurveSegFile" <>+    OP.metavar "FILE" <>+    OP.help "Path to an output file to store the relevant, interpolated calibration curve \+            \segment for the first (!) input date. \+            \This option as well as --calCurveMatFile are meant for debugging." <>     OP.value Nothing     )  optParseCalCurveMatrixFile :: OP.Parser (Maybe FilePath) optParseCalCurveMatrixFile = OP.option (Just <$> OP.str) (-    OP.long "calCurveMatrixFile" <>+    OP.long "calCurveMatFile" <>+    OP.metavar "FILE" <>     OP.help "Path to an output file which stores the relevant, interpolated calibration curve \-            \segment for the first (!) input date in a wide matrix format" <>+            \segment for the first (!) input date in a wide matrix format." <>     OP.value Nothing     )

@@ -45,7 +45,13 @@     CalExpr (..),     addPDFs,     multiplyPDFs,-    normalizeCalPDF+    normalizeCalPDF,++    -- * Drawing random samples from CalPDFs+    -- $randsamp+    AgeSamplingConf (..),+    sampleAgesFromCalPDF,+    RandomAgeSample (..)     ) where  import           Currycarbon.CalCurves.Intcal20@@ -127,4 +133,16 @@ 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.+-}++{- $randsamp++Another common requirement for archaeological data analysis is temporal resampling,+where random age samples are drawn from 'CalPDF's according to the probability+density distribution.++currycarbon supports this with 'sampleAgesFromCalPDF', which takes a configuration+data type 'AgeSamplingConf' including a random number generator and the number of+requested age samples, and an arbitrary 'CalPDF'. It returns an object of type+'RandomAgeSample' with a vector of sampled 'YearBCAD's. -}

@@ -1,5 +1,3 @@-{-# LANGUAGE BangPatterns #-}- module Currycarbon.CLI.RunCalibrate     (CalibrateOptions (..), runCalibrate) where @@ -10,14 +8,16 @@ import           Currycarbon.Types import           Currycarbon.Utils +import           Control.Exception                   (throwIO) import           Control.Monad                       (unless, when) import           Data.Maybe                          (fromJust, fromMaybe,                                                       isJust)-import           System.IO                           (hPutStrLn, stderr, stdout)+import           System.IO                           (hPutStrLn, stderr)+import qualified System.Random                       as R  -- | A data type to represent the options to the CLI module function runCalibrate data CalibrateOptions = CalibrateOptions {-        _calibrateExprs                   :: [CalExpr] -- ^ String listing the uncalibrated dates that should be calibrated+        _calibrateExprs                   :: [NamedCalExpr] -- ^ 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@@ -27,17 +27,27 @@       , _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+      , _calibrateAgeSampling             :: Maybe (Maybe Word, Word, FilePath) -- ^ Settings for the age sampling       , _calibrateCalCurveSegmentFile     :: Maybe FilePath -- ^ Path to an output file       , _calibrateCalCurveMatrixFile      :: Maybe FilePath -- ^ Path to an output file     }  -- | Interface function to trigger calibration from the command line runCalibrate :: CalibrateOptions -> IO ()-runCalibrate (CalibrateOptions exprs exprFiles calCurveFile method allowOutside noInterpolate quiet encoding densityFile hdrFile calCurveSegmentFile calCurveMatrixFile) = do+runCalibrate (+        CalibrateOptions+            exprs exprFiles+            calCurveFile method allowOutside noInterpolate+            quiet encoding+            densityFile hdrFile+            ageSampling+            calCurveSegmentFile calCurveMatrixFile+        ) = do     let ascii = encoding /= "UTF-8"     -- compile dates-    exprsFromFile <- mapM readCalExprFromFile exprFiles-    let exprsRenamed = replaceEmptyNames $ exprs ++ concat exprsFromFile+    exprsFromFile <- mapM readNamedCalExprsFromFile exprFiles+    let exprsCombined = exprs ++ concat exprsFromFile+        exprsRenamed = replaceEmptyNames exprsCombined     if null exprsRenamed     then hPutStrLn stderr "Nothing to calibrate. See currycarbon -h for help"     else do@@ -50,68 +60,131 @@             , _calConfAllowOutside = allowOutside             , _calConfInterpolateCalCurve = not noInterpolate             }+        -- handle the special debug cases+        when (isJust calCurveSegmentFile || isJust calCurveMatrixFile) $ do+            case exprsRenamed of+                [NamedCalExpr _ (UnCalDate uncal)] -> do+                    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+                _ -> do+                    throwIO $ CurrycarbonCLIException+                        "--calCurveSegFile and --calCurveMatFile only work with \+                        \a single uncalibrated radiocarbon date."         -- run calibration         hPutStrLn stderr "Calibrating..."-        let errorOrCalPDFs = map (evalCalExpr calConf calCurve) exprsRenamed-        handleDates ascii True calCurve $ zip exprsRenamed errorOrCalPDFs+        let errorOrCalPDFs = map (evalNamedCalExpr calConf calCurve) exprsRenamed+        -- prepare random number generator for age sampling+        maybeRNG <- case ageSampling of+            Nothing -> pure Nothing+            Just (maybeSeed, _, _) -> case maybeSeed of+                Nothing   -> Just <$> R.initStdGen+                Just seed -> return $ Just $ R.mkStdGen (fromIntegral seed)+        -- prepare and write the output per expression+        handleExprs ascii True calCurve maybeRNG $ 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++        -- loop over first and subsequent expressions+        handleExprs ::+               Bool -- encoding+            -> Bool -- is this expression the first in the list of expressions?+            -> CalCurveBP+            -> Maybe R.StdGen -- rng for the age sampling seeds+            -> [(NamedCalExpr, Either CurrycarbonException CalPDF)]+            -> IO ()+        handleExprs _ _ _ _ [] = hPutStrLn stderr "Done."+        -- first expression+        handleExprs _ascii True calCurve maybeRNG (firstDate:otherDates) =+            case firstDate of+                (_, Left e) -> do+                    printE e+                    handleExprs _ascii True calCurve maybeRNG otherDates+                (namedCalExpr, Right cPDF) -> do+                    let (sampleSeed, newRNG) = drawSeed maybeRNG+                    flexOut _ascii namedCalExpr cPDF sampleSeed writeCalPDF writeCalC14 writeRandomAgeSample+                    handleExprs _ascii False calCurve newRNG otherDates+        -- subsequent expression+        handleExprs _ascii False calCurve maybeRNG (nextDate:otherDates) =+            case nextDate of+                (_, Left e) -> do+                    printE e+                    handleExprs _ascii False calCurve maybeRNG otherDates+                (namedCalExpr, Right cPDF) -> do+                    let (sampleSeed, newRNG) = drawSeed maybeRNG+                    flexOut _ascii namedCalExpr cPDF sampleSeed appendCalPDF appendCalC14 appendRandomAgeSample+                    handleExprs _ascii False calCurve newRNG otherDates+         printE :: CurrycarbonException -> IO ()         printE e = hPutStrLn stderr $ renderCurrycarbonException e +        drawSeed :: Maybe R.StdGen -> (Maybe Int, Maybe R.StdGen)+        drawSeed maybeRNG = (\x -> (fromIntegral . fst <$> x, snd <$> x)) (R.genWord32 <$> maybeRNG)++        -- flexible expression handler+        flexOut ::+               Bool+            -> NamedCalExpr+            -> CalPDF+            -> Maybe Int+            -> (FilePath -> CalPDF -> IO ())+            -> (FilePath -> CalC14 -> IO ())+            -> (FilePath -> RandomAgeSample -> IO ())+            -> IO ()+        flexOut _ascii namedCalExpr calPDF maybeSeed calPDFToFile calC14ToFile randomAgeSampleToFile = do+            case refineCalDate calPDF of+                Left e -> do+                    unless quiet $ do+                        putStrLn ("CalEXPR: " ++ renderNamedCalExpr namedCalExpr)+                        printE e+                    when (isJust hdrFile) $ unless quiet $+                        hPutStrLn stderr "<!> Error: Can not create --hdrFile"+                Right calC14 -> do+                    unless quiet $ do+                        putStrLn (renderCalDatePretty _ascii (namedCalExpr, calPDF, calC14))+                    when (isJust hdrFile) $+                        calC14ToFile (fromJust hdrFile) calC14+            when (isJust ageSampling && isJust maybeSeed) $ do+                let (_, nrOfSamples, path) = fromJust ageSampling+                    rng = R.mkStdGen (fromJust maybeSeed)+                    conf = AgeSamplingConf rng nrOfSamples+                case sampleAgesFromCalPDF conf calPDF of+                    Left e -> do+                        unless quiet $ do+                            printE e+                            hPutStrLn stderr "<!> Error: Can not create --samplesFile"+                    Right res -> randomAgeSampleToFile path res+            when (isJust densityFile) $+                calPDFToFile (fromJust densityFile) calPDF+ -- | Helper function to replace empty input names with a sequence of numbers, -- to get each input date an unique identifier-replaceEmptyNames :: [CalExpr] -> [CalExpr]-replaceEmptyNames = zipWith (replaceName . show) ([1..] :: [Integer])+replaceEmptyNames :: [NamedCalExpr] -> [NamedCalExpr]+replaceEmptyNames = zipWith (modifyNamedExpr . show) ([1..] :: [Integer])     where+        modifyNamedExpr :: String -> NamedCalExpr -> NamedCalExpr+        modifyNamedExpr i nexpr =+            if _exprID nexpr == ""+            then nexpr { _exprID = i, _expr = replaceName i (_expr nexpr) }+            else nexpr {              _expr = replaceName i (_expr nexpr) }         replaceName :: String -> CalExpr -> CalExpr         replaceName i (UnCalDate (UncalC14 name x y)) =-            if name == "unknownSampleName"+            if name == ""             then UnCalDate $ UncalC14 i x y             else UnCalDate $ UncalC14 name x y+        replaceName i (WindowBP (TimeWindowBP name start stop)) =+            if name == ""+            then WindowBP $ TimeWindowBP i start stop+            else WindowBP $ TimeWindowBP name start stop+        replaceName i (WindowBCAD (TimeWindowBCAD name start stop)) =+            if name == ""+            then WindowBCAD $ TimeWindowBCAD i start stop+            else WindowBCAD $ TimeWindowBCAD name start stop         replaceName i (CalDate (CalPDF name x y)) =-            if name == "unknownSampleName"+            if name == ""             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)

@@ -16,6 +16,8 @@       , refineCalDate       , CalibrateDatesConf (..)       , defaultCalConf+      , AgeSamplingConf (..)+      , sampleAgesFromCalPDF     ) where  import           Currycarbon.Calibration.Bchron@@ -24,10 +26,12 @@ import           Currycarbon.Types import           Currycarbon.Utils +import qualified Control.Monad.Random               as CMR import           Data.List                          (elemIndex, groupBy, sort,                                                      sortBy) import           Data.Maybe                         (fromJust) import qualified Data.Vector.Unboxed                as VU+import qualified System.Random                      as R  -- | A data type to cover the configuration options of the calibrateDates function data CalibrateDatesConf = CalibrateDatesConf {@@ -79,29 +83,49 @@  -- | Transforms the raw, calibrated probability density table to a meaningful representation of a -- calibrated radiocarbon date-refineCalDates :: [CalPDF] -> [Maybe CalC14]+refineCalDates :: [CalPDF] -> [Either CurrycarbonException CalC14] refineCalDates = map refineCalDate -refineCalDate :: CalPDF -> Maybe CalC14-refineCalDate (CalPDF name cals dens) =-    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 {+refineCalDate :: CalPDF -> Either CurrycarbonException CalC14+refineCalDate calPDF@(CalPDF name cals dens)+    -- don't calculate CalC14, if it's not meaningful+    | isInvalidCalPDF calPDF =+        Left $ CurrycarbonInvalidCalPDFException "refinement"+    -- for simple uniform age ranges+    | VU.length (VU.uniq dens) == 1 =+        let start = VU.head cals+            stop  = VU.last cals+        in Right $ CalC14 {           _calC14id           = name         , _calC14RangeSummary = CalRangeSummary {+              _calRangeStartTwoSigma = start+            , _calRangeStartOneSigma = start+            , _calRangeMedian        = median+            , _calRangeStopOneSigma  = stop+            , _calRangeStopTwoSigma  = stop+            }+        , _calC14HDROneSigma  = [HDR start stop]+        , _calC14HDRTwoSigma  = [HDR start stop]+        }+    -- for normal post-calibration probability distributions+    | otherwise =+        Right $ CalC14 {+          _calC14id           = name+        , _calC14RangeSummary = CalRangeSummary {               _calRangeStartTwoSigma = _hdrstart $ head hdrs95             , _calRangeStartOneSigma = _hdrstart $ head hdrs68-            , _calRangeMedian        = fromJust $ cals `indexVU` elemIndex (minimum distanceTo05) distanceTo05+            , _calRangeMedian        = median             , _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+        median = fromJust $ cals `indexVU` elemIndex (minimum distanceTo05) distanceTo05         -- sorted density cumsum for hdrs         sortedDensities = sortBy (flip (\ (_, dens1) (_, dens2) -> compare dens1 dens2)) (VU.toList $ VU.zip cals dens)         cumsumSortedDensities = cumsumDens sortedDensities@@ -131,3 +155,30 @@         getIn95 (_,_,_,x) = x         getYear :: (Int, Float, Bool, Bool) -> Int         getYear (year,_,_,_) = year++-- age sampling++-- | A data type to define the settings for age sampling+data AgeSamplingConf = AgeSamplingConf {+    -- | Random number generator+      _assRNG             :: R.StdGen+    -- | Number of samples that should be drawn per sample+    , _assNumberOfSamples :: Word+    } deriving (Show, Eq)++-- | Draw random samples from a probability density table+sampleAgesFromCalPDF :: AgeSamplingConf -> CalPDF -> Either CurrycarbonException RandomAgeSample+sampleAgesFromCalPDF (AgeSamplingConf rng n) calPDF@(CalPDF calPDFid cals dens) =+    let weightedList = zip (VU.toList cals) (map toRational $ VU.toList dens)+        infSamplesList = sampleWeightedList rng weightedList+        samples = take (fromIntegral n) infSamplesList+    in if isInvalidCalPDF calPDF+       then Left $ CurrycarbonInvalidCalPDFException "random age sampling"+       else Right $ RandomAgeSample calPDFid (VU.fromList samples)+    where+        sampleWeightedList :: CMR.RandomGen g => g -> [(a, Rational)] -> [a]+        sampleWeightedList gen weights = CMR.evalRand m gen+            where m = sequence . repeat . CMR.fromList $ weights++isInvalidCalPDF :: CalPDF -> Bool+isInvalidCalPDF (CalPDF _ _ dens) = VU.sum dens == 0 || VU.any (>= 1.0) dens

@@ -63,7 +63,10 @@ makeBCADCalCurve (CalCurveBP cals uncals sigmas) = CalCurveBCAD (vectorBPToBCAD cals) (vectorBPToBCAD uncals) sigmas  vectorBPToBCAD :: VU.Vector YearBP -> VU.Vector YearBCAD-vectorBPToBCAD = VU.map (\x -> -(fromIntegral x) + 1950)+vectorBPToBCAD = VU.map bp2BCAD++bp2BCAD :: YearBP -> YearBCAD+bp2BCAD x = -(fromIntegral x) + 1950  interpolateCalCurve :: CalCurveBP -> CalCurveBP interpolateCalCurve (CalCurveBP cals uncals sigmas) =

@@ -0,0 +1,168 @@+module Currycarbon.ParserHelpers where++import qualified Text.Parsec        as P+import qualified Text.Parsec.Error  as P+import qualified Text.Parsec.String as P++-- * High level building blocks++parseRecordType :: String -> P.Parser a -> P.Parser a+parseRecordType typeName parser = do+    _ <- P.string typeName+    parseInParens parser++parseNamedVector :: P.Parser a -> P.Parser b -> P.Parser [(a,b)]+parseNamedVector parseKey parseValue =+    parseVector $ parseKeyValuePair parseKey parseValue++parseVector :: P.Parser a -> P.Parser [a]+parseVector parser = do+    _ <- P.char 'c'+    parseInParens (P.sepBy parser consumeCommaSep)++parseArgumentWithDefault :: String -> P.Parser b -> b -> P.Parser b+parseArgumentWithDefault argumentName parseValue defaultValue =+    P.option defaultValue (parseArgument argumentName parseValue)++parseArgumentOptional :: String -> P.Parser b -> P.Parser (Maybe b)+parseArgumentOptional argumentName parseValue =+    P.optionMaybe $ P.try (parseArgument argumentName parseValue)++parseArgument :: String -> P.Parser b -> P.Parser b+parseArgument argumentName parseValue = do+    res <- parseArgumentWithoutComma argumentName parseValue+    P.optional consumeCommaSep+    return res++parseNamedArgumentOptional :: String -> P.Parser b -> P.Parser (Maybe b)+parseNamedArgumentOptional argumentName parseValue =+    P.optionMaybe $ P.try (parseNamedArgument argumentName parseValue)++-- * Low level blocks++parseArgumentWithoutComma :: String -> P.Parser b -> P.Parser b+parseArgumentWithoutComma argumentName parseValue =+    P.try (parseNamedArgument argumentName parseValue) P.<|> parseUnnamedArgument parseValue++parseNamedArgument :: String -> P.Parser b -> P.Parser b+parseNamedArgument argumentName parseValue = do+    (_,b) <- parseKeyValuePair (P.string argumentName) parseValue+    return b++parseUnnamedArgument :: P.Parser b -> P.Parser b+parseUnnamedArgument parseValue = parseValue++parseKeyValuePair :: P.Parser a -> P.Parser b -> P.Parser (a,b)+parseKeyValuePair parseKey parseValue = do+    key <- parseKey+    consumeEqualSep+    value <- parseValue+    return (key, value)++parseInParens :: P.Parser b -> P.Parser b+parseInParens parser = do+    _ <- P.char '('+    _ <- P.spaces+    res <- parser+    _ <- P.spaces+    _ <- P.char ')'+    return res++consumeEqualSep :: P.Parser ()+consumeEqualSep = do+    _ <- P.spaces *> P.char '=' <* P.spaces+    return ()+consumeCommaSep :: P.Parser ()+consumeCommaSep = do+    _ <- P.spaces *> P.char ',' <* P.spaces+    return ()++parseCharInSpace :: Char -> P.Parser Char+parseCharInSpace c = P.between P.spaces P.spaces (P.char c)++parseAnyString :: P.Parser String+parseAnyString =+    P.try inDoubleQuotes P.<|> P.try inSingleQuotes P.<|> inNoQuotes+    where+        inDoubleQuotes = P.between (P.char '"') (P.char '"') (P.many P.anyChar)+        inSingleQuotes = P.between (P.char '\'') (P.char '\'') (P.many P.anyChar)+        inNoQuotes = P.many (P.noneOf ",):")++-- * Sequence parsers++parseDoubleSequence :: P.Parser [Double]+parseDoubleSequence = do+    start <- parseDouble+    _ <- P.oneOf ":"+    stop <- parseDouble+    _ <- P.oneOf ":"+    by <- parsePositiveFloatNumber+    return [start,(start+by)..stop]++-- * Number parsers++parseDouble :: P.Parser Double+parseDouble = do+    P.try parseNegativeFloatNumber P.<|> parsePositiveFloatNumber++parseNegativeFloatNumber :: P.Parser Double+parseNegativeFloatNumber = do+    _ <- P.oneOf "-"+    i <- parsePositiveFloatNumber+    return (-i)++parseFraction :: P.Parser Double+parseFraction = do+    num <- parsePositiveFloatNumber+    if num > 1+    then fail "must be between zero and one"+    else return num++parsePositiveFloatNumber :: P.Parser Double+parsePositiveFloatNumber = do+    num <- parseNumber+    optionalMore <- P.option "" $ (:) <$> P.char '.' <*> parseNumber+    return $ read $ num ++ optionalMore++parseIntegerSequence :: P.Parser [Int]+parseIntegerSequence = do+    start <- parseInt+    _ <- P.oneOf ":"+    stop <- parseInt+    _ <- P.oneOf ":"+    by <- fromIntegral <$> parsePositiveInt+    return [start,(start+by)..stop]++parseInt :: P.Parser Int+parseInt = do+    P.try parseNegativeInt P.<|> parsePositiveInt++parseNegativeInt :: P.Parser Int+parseNegativeInt = do+    _ <- P.oneOf "-"+    i <- parsePositiveInt+    return (-i)++parsePositiveInt :: P.Parser Int+parsePositiveInt = fromIntegral <$> parseWord++-- https://hackage.haskell.org/package/base-4.19.0.0/docs/Data-Word.html+parseWord :: P.Parser Word+parseWord = do+    read <$> parseNumber++parsePositiveDouble :: P.Parser Double+parsePositiveDouble = do+    read <$> parseNumber++parseNumber :: P.Parser [Char]+parseNumber = P.many1 P.digit++-- * Error helpers++showParsecErr :: P.ParseError -> String+showParsecErr err =+    P.showErrorMessages+        "or" "unknown parse error"+        "expecting" "unexpected" "end of input"+        (P.errorMessages err)

@@ -2,15 +2,16 @@  module Currycarbon.Parsers where +import           Currycarbon.ParserHelpers import           Currycarbon.Types import           Currycarbon.Utils -import           Control.Exception   (throwIO)-import           Data.List           (intercalate, transpose)-import qualified Data.Vector         as V-import qualified Data.Vector.Unboxed as VU-import qualified Text.Parsec         as P-import qualified Text.Parsec.String  as P+import           Control.Exception         (throwIO)+import           Data.List                 (intercalate, transpose)+import qualified Data.Vector               as V+import qualified Data.Vector.Unboxed       as VU+import qualified Text.Parsec               as P+import qualified Text.Parsec.String        as P  -- * Parsing, rendering and writing functions --@@ -19,12 +20,13 @@ -- This module contains a number of functions to manage data input and -- output plumbing for different datatypes --- CalibrationMethod+-- read the calibration method+ readCalibrationMethod :: String -> Either String CalibrationMethod readCalibrationMethod s =     case P.runParser parseCalibrationMethod () "" s of-        Left err -> Left $ renderCurrycarbonException $ CurrycarbonCLIParsingException $ show err-        Right x -> Right x+        Left err -> Left $ showParsecErr err+        Right x  -> Right x  parseCalibrationMethod :: P.Parser CalibrationMethod parseCalibrationMethod = do@@ -35,7 +37,7 @@             P.try studentT P.<|> normal         studentT = do             _ <- P.string "StudentT,"-            dof <- read <$> P.many1 P.digit+            dof <- parsePositiveDouble             return (Bchron $ StudentTDist dof)         normal = do             _ <- P.string "Normal"@@ -44,11 +46,13 @@             _ <- P.string "MatrixMult"             return MatrixMultiplication +-- pretty printing+ -- | Combine 'CalExpr', 'CalPDF' and 'CalC14' to render pretty command line output -- like this: -- -- @--- DATE: (5000±30BP + 5100±100BP)+-- CalEXPR: [Ex2] (S1:5000±30BP + S2:5100±100BP) -- Calibrated: 4150BC \>\> 3941BC \> 3814BC \< 3660BC \<\< 3651BC -- 1-sigma: 3941-3864BC, 3810-3707BC, 3667-3660BC -- 2-sigma: 4150-4148BC, 4048-3651BC@@ -67,90 +71,230 @@ -- renderCalDatePretty ::        Bool -- ^ Should the CLI plot be restricted to (boring) ASCII symbols?-    -> (CalExpr, CalPDF, CalC14)+    -> (NamedCalExpr, CalPDF, CalC14)     -> String renderCalDatePretty ascii (calExpr, calPDF, calC14) =-    "DATE: " ++ intercalate "\n" [-          renderCalExpr calExpr+    "CalEXPR: " ++ intercalate "\n" [+          renderNamedCalExpr calExpr         , renderCalC14 calC14         , renderCLIPlotCalPDF ascii 6 50 calPDF calC14         ] +-- write and read calibration expressions++renderNamedCalExpr :: NamedCalExpr -> String+renderNamedCalExpr (NamedCalExpr exprID calExpr) = renderExprID exprID ++ " " ++ renderCalExpr calExpr++renderExprID :: String -> String+renderExprID s = "[" ++ s ++ "]"+ renderCalExpr :: CalExpr -> String renderCalExpr (UnCalDate a)               = renderUncalC14 a+renderCalExpr (WindowBP a)                = renderTimeWindowBP a+renderCalExpr (WindowBCAD a)              = renderTimeWindowBCAD a renderCalExpr (CalDate (CalPDF name _ _)) = name renderCalExpr (SumCal a b)                = "(" ++ renderCalExpr a ++ " + " ++ renderCalExpr b ++ ")" renderCalExpr (ProductCal a b)            = "(" ++ renderCalExpr a ++ " * " ++ renderCalExpr b ++ ")" +renderTimeWindowBP :: TimeWindowBP -> String+renderTimeWindowBP (TimeWindowBP name start stop) =+    name ++ ":" ++ renderYearBP start ++ "-" ++ renderYearBP stop++renderTimeWindowBCAD :: TimeWindowBCAD -> String+renderTimeWindowBCAD (TimeWindowBCAD name start stop) =+    name ++ ":" ++ renderYearBCAD start ++ "-" ++ renderYearBCAD stop++parseTimeWindowBP :: P.Parser TimeWindowBP+parseTimeWindowBP = parseRecordType "rangeBP" $ P.try long P.<|> short+    where+        long = do+            name  <- parseArgument "id" parseAnyString+            start <- parseArgument "start" parseWord+            stop  <- parseArgument "stop" parseWord+            construct name start stop+        short = do+            start <- parseArgument "start" parseWord+            stop  <- parseArgument "stop" parseWord+            construct "" start stop+        construct name start stop = do+            if start >= stop+            then return (TimeWindowBP name start stop)+            else fail "the BP stop date can not be larger than the start date"++parseTimeWindowBCAD :: P.Parser TimeWindowBCAD+parseTimeWindowBCAD = parseRecordType "rangeBCAD" $ P.try long P.<|> short+    where+        long = do+            name  <- parseArgument "id" parseAnyString+            start <- parseArgument "start" parseInt+            stop  <- parseArgument "stop" parseInt+            construct name start stop+        short = do+            start <- parseArgument "start" parseInt+            stop  <- parseArgument "stop" parseInt+            construct "" start stop+        construct name start stop = do+            if start <= stop+            then return (TimeWindowBCAD name start stop)+            else fail "the BC/AD stop date can not be smaller than the start date"+ -- 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+addFun :: P.Parser CalExpr+addFun = parseRecordType "sum" $ do+    a <- parseArgument "a" term+    b <- parseArgument "b" expr+    return $ SumCal a b -add :: P.Parser CalExpr-add = SumCal <$> term <*> (spaceChar '+' *> expr)+addOperator :: P.Parser CalExpr+addOperator = SumCal <$> term <*> (parseCharInSpace '+' *> expr) -mul :: P.Parser CalExpr-mul = ProductCal <$> factor <*> (spaceChar '*' *> term)+mulFun :: P.Parser CalExpr+mulFun = parseRecordType "product" $ do+    a <- parseArgument "a" factor+    b <- parseArgument "b" term+    return $ ProductCal a b +mulOperator :: P.Parser CalExpr+mulOperator = ProductCal <$> factor <*> (parseCharInSpace '*' *> term)+ parens :: P.Parser CalExpr-parens = P.between (spaceChar '(') (spaceChar ')') expr+parens = P.between (parseCharInSpace '(') (parseCharInSpace ')') expr  factor :: P.Parser CalExpr-factor = parens P.<|> (UnCalDate <$> parseUncalC14)+factor =      P.try parens+        P.<|> P.try addFun+        P.<|> P.try mulFun+        P.<|> P.try (WindowBP <$> parseTimeWindowBP)+        P.<|> P.try (WindowBCAD <$> parseTimeWindowBCAD)+        P.<|> (UnCalDate <$> parseUncalC14)  term :: P.Parser CalExpr-term = P.try mul P.<|> factor+term = P.try mulOperator P.<|> factor  expr :: P.Parser CalExpr-expr = P.try add P.<|> term -- <* P.eof+expr = P.try addOperator P.<|> term -- <* P.eof -readCalExpr :: String -> Either String [CalExpr]-readCalExpr s =+namedExpr :: P.Parser NamedCalExpr+namedExpr = P.try nameBeforeColon P.<|> P.try record P.<|> onlyExpr+    where+        nameBeforeColon = do+            name <- parseAnyString+            _ <- P.char ':'+            _ <- P.spaces+            ex <- expr+            return (NamedCalExpr name ex)+        record = parseRecordType "calExpr" $ P.try long P.<|> short+        long = do+            name <- parseArgument "id" parseAnyString+            ex   <- parseArgument "expr" expr+            return (NamedCalExpr name ex)+        short = do+            ex   <- parseArgument "expr" expr+            return (NamedCalExpr "" ex)+        onlyExpr = NamedCalExpr "" <$> expr++readNamedCalExprs :: String -> Either String [NamedCalExpr]+readNamedCalExprs s =     case P.runParser parseCalExprSepBySemicolon () "" s of-        Left err -> Left $ renderCurrycarbonException $ CurrycarbonCLIParsingException $ show err-        Right x -> Right x+        Left err -> Left $ showParsecErr err+        Right x  -> Right x         where-        parseCalExprSepBySemicolon :: P.Parser [CalExpr]-        parseCalExprSepBySemicolon = P.sepBy expr (P.char ';' <* P.spaces) <* P.eof+        parseCalExprSepBySemicolon :: P.Parser [NamedCalExpr]+        parseCalExprSepBySemicolon = P.sepBy namedExpr (P.char ';' <* P.spaces) <* P.eof -readCalExprFromFile :: FilePath -> IO [CalExpr]-readCalExprFromFile uncalFile = do+readOneNamedCalExpr :: String -> Either String NamedCalExpr+readOneNamedCalExpr s =+    case P.runParser namedExpr () "" s of+        Left err -> Left $ showParsecErr err+        Right x  -> Right x++readNamedCalExprsFromFile :: FilePath -> IO [NamedCalExpr]+readNamedCalExprsFromFile uncalFile = do+    ss <- lines <$> readFile uncalFile+    mapM readOneLine ss+    where+        readOneLine :: String -> IO NamedCalExpr+        readOneLine s = case readOneNamedCalExpr s of+            Left err -> throwIO $ CurrycarbonCLIParsingException $ err ++ "\nin \"" ++ s ++ "\""+            Right x  -> return x++-- UncalC14+renderUncalC14WithoutName :: UncalC14 -> String+renderUncalC14WithoutName (UncalC14 _ bp sigma) = show bp ++ "±" ++ show sigma ++ "BP"++renderUncalC14 :: UncalC14 -> String+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+-- \<sample name\> is optional. A valid file could look like this:+--+-- @+-- Sample1,5000,30+-- 6000,50+-- Sample3,4000,25+-- @+--+readUncalC14FromFile :: FilePath -> IO [UncalC14]+readUncalC14FromFile uncalFile = do     s <- readFile uncalFile-    case P.runParser parseCalExprSepByNewline () "" s of-        Left err -> throwIO $ CurrycarbonCLIParsingException $ show err+    case P.runParser uncalC14SepByNewline () "" s of+        Left err -> throwIO $ CurrycarbonCLIParsingException $ showParsecErr err         Right x  -> return x     where-        parseCalExprSepByNewline :: P.Parser [CalExpr]-        parseCalExprSepByNewline = P.endBy expr (P.newline <* P.spaces) <* P.eof+        uncalC14SepByNewline :: P.Parser [UncalC14]+        uncalC14SepByNewline = P.endBy parseUncalC14 (P.newline <* P.spaces) <* P.eof +readUncalC14 :: String -> Either String [UncalC14]+readUncalC14 s =+    case P.runParser uncalC14SepBySemicolon () "" s of+        Left err -> Left $ showParsecErr err+        Right x  -> Right x+    where+        uncalC14SepBySemicolon :: P.Parser [UncalC14]+        uncalC14SepBySemicolon = P.sepBy parseUncalC14 (P.char ';' <* P.spaces) <* P.eof++parseUncalC14 :: P.Parser UncalC14+parseUncalC14 = P.try record P.<|> P.try long P.<|> short+    where+        record = parseRecordType "uncalC14" $ P.try long P.<|> short+        long = do+            name  <- parseArgument "id" parseAnyString+            age   <- parseArgument "yearBP" parseWord+            sigma <- parseArgument "sigma" parseWord+            return (UncalC14 name age sigma)+        short = do+            age   <- parseArgument "yearBP" parseWord+            sigma <- parseArgument "sigma" parseWord+            return (UncalC14 "" age sigma)+ -- CalC14--- | Write 'CalC14's to the file system. The output file is a long .csv file with the following structure:+-- | Write 'CalC14's to the file system. The output file is a long .tsv file with the following structure: -- -- @--- sample,hdrSigma,hdrStartBCAD,hdrStopBCAD--- Sample1,1,-3797,-3709--- Sample1,1,-3894,-3880--- Sample1,2,-3680,-3655--- Sample1,2,-3810,-3700--- Sample1,2,-3941,-3864--- Sample2,1,-1142,-1130--- Sample2,1,-1173,-1161--- Sample2,1,-1293,-1194--- Sample2,1,-1368,-1356--- Sample2,2,-1061,-1059--- Sample2,2,-1323,-1112--- Sample2,2,-1393,-1334+-- id  hdrSigmaLevel  hdrStartYearBCAD  hdrStopYearBCAD+-- Sample1  1	-3797	-3709+-- Sample1  1	-3894	-3880+-- Sample1  2	-3680	-3655+-- Sample1  2	-3810	-3700+-- Sample1  2	-3941	-3864+-- Sample2  1	-1142	-1130+-- Sample2	1	-1173	-1161+-- Sample2	1	-1293	-1194+-- Sample2	1	-1368	-1356+-- Sample2	2	-1061	-1059+-- Sample2	2	-1323	-1112+-- Sample2	2	-1393	-1334 -- @ -- writeCalC14s :: FilePath -> [CalC14] -> IO () writeCalC14s path calC14s = writeFile path $-    "sample,hdrSigma,hdrStartBCAD,hdrStopBCAD\n"+    "id\thdrSigmaLevel\thdrStartYearBCAD\thdrStopYearBCAD\n"     ++ intercalate "\n" (map renderCalC14ForFile calC14s)  writeCalC14 :: FilePath -> CalC14 -> IO () writeCalC14 path calC14 = writeFile path $-    "sample,hdrSigma,hdrStartBCAD,hdrStopBCAD\n"+    "id\thdrSigmaLevel\thdrStartYearBCAD\thdrStopYearBCAD\n"     ++ renderCalC14ForFile calC14  appendCalC14 :: FilePath -> CalC14 -> IO ()@@ -165,7 +309,7 @@         zip3 (repeat name) (repeat "2") (renderHDRsForFile hdrs95)     where         renderRow :: (String, String, (String, String)) -> String-        renderRow (a, b, (c, d)) = intercalate "," [a,b,c,d]+        renderRow (a, b, (c, d)) = intercalate "\t" [a,b,c,d]  renderCalC14s :: [CalC14] -> String renderCalC14s xs =@@ -186,6 +330,11 @@     ++ renderYearBCAD (_calRangeStopOneSigma s) ++ " << "     ++ renderYearBCAD (_calRangeStopTwoSigma s) +-- BP+renderYearBP :: YearBP -> String+renderYearBP x =+    show x ++ "BP" -- ++ " (" ++ (renderYearBCAD $ bp2BCAD x) ++ ")"+ -- BCAD renderYearBCAD :: YearBCAD -> String renderYearBCAD x@@ -217,40 +366,40 @@  renderCalCurveMatrix :: CalCurveMatrix -> String renderCalCurveMatrix (CalCurveMatrix uncals cals curveDensities) =-    let header = "," ++ intercalate "," (map show $ VU.toList cals) ++ "\n"+    let header = "\t" ++ intercalate "\t" (map show $ VU.toList cals) ++ "\n"         body = zipWith makeRow (VU.toList uncals) (transpose $ V.toList (V.map VU.toList curveDensities))     in header ++ intercalate "\n" body     where-      makeRow uncal dens = show uncal ++ "," ++ intercalate "," (map show dens)+      makeRow uncal dens = show uncal ++ "\t" ++ intercalate "\t" (map show dens)  -- CalPDF--- | Write 'CalPDF's to the file system. The output file is a long .csv file with the following structure:+-- | Write 'CalPDF's to the file system. The output file is a long .tsv file with the following structure: -- -- @--- sample,calBCAD,density+-- id  yearBCAD  density -- ...--- Sample1,-1391,2.8917924e-4--- Sample1,-1390,3.3285577e-4--- Sample1,-1389,3.5674628e-4--- Sample1,-1388,3.750703e-4+-- Sample1  -1391   2.8917924e-4+-- Sample1  -1390   3.3285577e-4+-- Sample1  -1389   3.5674628e-4+-- Sample1  -1388   3.750703e-4 -- ...--- Sample2,-3678,1.8128564e-3--- Sample2,-3677,1.9512239e-3--- Sample2,-3676,2.0227064e-3--- Sample2,-3675,2.095691e-3+-- Sample2  -3678   1.8128564e-3+-- Sample2  -3677   1.9512239e-3+-- Sample2  -3676   2.0227064e-3+-- Sample2  -3675   2.095691e-3 -- ... -- @ -- writeCalPDFs :: FilePath -> [CalPDF] -> IO () writeCalPDFs path calPDFs =     writeFile path $-        "sample,calBCAD,density\n"+        "id\tyearBCAD\tdensity\n"         ++ renderCalPDFs calPDFs  writeCalPDF :: FilePath -> CalPDF -> IO () writeCalPDF path calPDF =     writeFile path $-        "sample,calBCAD,density\n"+        "id\tyearBCAD\tdensity\n"         ++ renderCalPDF calPDF  appendCalPDF :: FilePath -> CalPDF -> IO ()@@ -264,8 +413,9 @@ renderCalPDF (CalPDF name cals dens) =     concatMap makeRow $ VU.toList $ VU.zip cals dens     where-      makeRow (x,y) = show name ++ "," ++ show x ++ "," ++ show y ++ "\n"+      makeRow (x,y) = name ++ "\t" ++ show x ++ "\t" ++ show y ++ "\n" +-- cli plot data PlotSymbol = HistFill | HistTop | AxisEnd | AxisLine | AxisTick | HDRLine  renderCLIPlotCalPDF :: Bool -> Int -> Int -> CalPDF -> CalC14 -> String@@ -356,59 +506,6 @@                             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) = 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--- \<sample name\> is optional. A valid file could look like this:------ @--- Sample1,5000,30--- 6000,50--- Sample3,4000,25--- @----readUncalC14FromFile :: FilePath -> IO [UncalC14]-readUncalC14FromFile uncalFile = do-    s <- readFile uncalFile-    case P.runParser uncalC14SepByNewline () "" s of-        Left err -> throwIO $ CurrycarbonCLIParsingException $ show err-        Right x  -> return x-    where-        uncalC14SepByNewline :: P.Parser [UncalC14]-        uncalC14SepByNewline = P.endBy parseUncalC14 (P.newline <* P.spaces) <* P.eof--readUncalC14 :: String -> Either String [UncalC14]-readUncalC14 s =-    case P.runParser uncalC14SepBySemicolon () "" s of-        Left err -> Left $ renderCurrycarbonException $ CurrycarbonCLIParsingException $ show err-        Right x -> Right x-    where-        uncalC14SepBySemicolon :: P.Parser [UncalC14]-        uncalC14SepBySemicolon = P.sepBy parseUncalC14 (P.char ';' <* P.spaces) <* P.eof--parseUncalC14 :: P.Parser UncalC14-parseUncalC14 = do-    P.try long P.<|> short-    where-        long = do-            name <- P.many (P.noneOf ",")-            _ <- P.oneOf ","-            mean <- read <$> P.many1 P.digit-            _ <- P.oneOf ","-            std <- read <$> P.many1 P.digit-            return (UncalC14 name mean std)-        short = do-            mean <- read <$> P.many1 P.digit-            _ <- P.oneOf ","-            std <- read <$> P.many1 P.digit-            return (UncalC14 "unknownSampleName" mean std)- -- CalCurve writeCalCurve :: FilePath -> CalCurveBCAD -> IO () writeCalCurve path calCurve =@@ -416,11 +513,11 @@  renderCalCurve :: CalCurveBCAD -> String renderCalCurve (CalCurveBCAD cals uncals sigmas) =-    let header = "calBCAD,uncalBCAD,Sigma\n"+    let header = "calYearBCAD\tuncalYearBCAD\tsigma\n"         body = map makeRow $ VU.toList $ VU.zip3 cals uncals sigmas     in header ++ intercalate "\n" body     where-      makeRow (x,y,z) = show x ++ "," ++ show y ++ "," ++ show z+      makeRow (x,y,z) = show x ++ "\t" ++ show y ++ "\t" ++ show z  -- | Read a calibration curve file. The file must adhere to the current version of the -- .c14 file format (e.g. [here](http://intcal.org/curves/intcal20.14c)). Look@@ -446,11 +543,11 @@  parseCalCurveLine :: P.Parser (YearBP, YearBP, YearRange) parseCalCurveLine = do-  calBP <- read <$> P.many1 P.digit+  calBP <- parseWord   _ <- P.oneOf ","-  bp <- read <$> P.many1 P.digit+  bp <- parseWord   _ <- P.oneOf ","-  sigma <- read <$> P.many1 P.digit+  sigma <- parseWord   return (calBP, bp, sigma)  comments :: P.Parser String@@ -458,3 +555,46 @@     _ <- P.string "#"     _ <- P.manyTill P.anyChar P.newline     return ""++-- RandomAgeSamples+-- | Write 'RandomAgeSamples's to the file system. The output file is a long .tsv file with the following structure:+--+-- @+-- id  yearBCAD+-- ...+-- Sample1  -1221+-- Sample1  -1211+-- Sample1  -1230+-- Sample1  -1225+-- ...+-- Sample2  -3763+-- Sample2  -3788+-- Sample2  -3767+-- Sample2  -3774+-- ...+-- @+--+writeRandomAgeSamples :: FilePath -> [RandomAgeSample] -> IO ()+writeRandomAgeSamples path calPDFs =+    writeFile path $+        "id\tyearBCAD\n"+        ++ renderRandomAgeSamples calPDFs++writeRandomAgeSample :: FilePath -> RandomAgeSample -> IO ()+writeRandomAgeSample path calPDF =+    writeFile path $+        "id\tyearBCAD\n"+        ++ renderRandomAgeSample calPDF++appendRandomAgeSample :: FilePath -> RandomAgeSample -> IO ()+appendRandomAgeSample path calPDF =+    appendFile path $ renderRandomAgeSample calPDF++renderRandomAgeSamples :: [RandomAgeSample] -> String+renderRandomAgeSamples = concatMap renderRandomAgeSample++renderRandomAgeSample :: RandomAgeSample -> String+renderRandomAgeSample (RandomAgeSample name samples) =+    concatMap makeRow $ VU.toList samples+    where+      makeRow x = name ++ "\t" ++ show x ++ "\n"

@@ -12,17 +12,30 @@ import           Data.Ord                            (comparing) import qualified Data.Vector.Unboxed                 as VU +evalNamedCalExpr :: CalibrateDatesConf -> CalCurveBP -> NamedCalExpr -> Either CurrycarbonException CalPDF+evalNamedCalExpr conf curve (NamedCalExpr exprID expr) =+    case evalCalExpr conf curve expr of+        Left err     -> Left err+        Right calPDF -> Right calPDF { _calPDFid = exprID }+ -- | 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+evalCalExpr conf curve calExpr = norm $ evalE calExpr     where         evalE :: CalExpr -> Either CurrycarbonException CalPDF+        -- these are already normalized by their constructors         evalE (UnCalDate a)    = calibrateDate conf curve a-        evalE (CalDate a)      = Right a+        evalE (WindowBP a)     = Right $ windowBP2CalPDF a+        evalE (WindowBCAD a)   = Right $ windowBCAD2CalPDF a+        -- this can theoretically be non-normalized input+        evalE (CalDate a)      = norm $ Right a+        -- sums must not be normalized         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+        -- products must be normalized (and their input, in case it's a sum)+        evalE (ProductCal a b) = norm $ eitherCombinePDFs (*) 1 (norm $ evalE a) (norm $ evalE b)+        norm :: Either CurrycarbonException CalPDF -> Either CurrycarbonException CalPDF+        norm = mapEither id normalizeCalPDF         -- 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)@@ -37,11 +50,11 @@ eitherCombinePDFs _ _ _ (Left e) = Left e eitherCombinePDFs f initVal (Right a) (Right b) = Right $ combinePDFs f initVal a b --- | Add two probabilty densities+-- | Add two probability densities addPDFs :: CalPDF -> CalPDF -> CalPDF addPDFs = combinePDFs (+) 0 --- | Multiply two probabilty densities+-- | Multiply two probability densities multiplyPDFs :: CalPDF -> CalPDF -> CalPDF multiplyPDFs = combinePDFs (*) 1 @@ -59,7 +72,7 @@             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)+        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@@ -69,3 +82,14 @@                 | otherwise = []             foldYearGroup :: [(YearBCAD, Float)] -> (YearBCAD, Float)             foldYearGroup oneYear = (fst $ head oneYear, foldl' f initVal $ map snd oneYear)++-- | Create pseudo-CalPDF from RangeBCAD+windowBCAD2CalPDF :: TimeWindowBCAD -> CalPDF+windowBCAD2CalPDF (TimeWindowBCAD name start stop) =+    let years = VU.fromList $ [start..stop]+        dens = VU.replicate (VU.length years) 1+    in normalizeCalPDF $ CalPDF name years dens++windowBP2CalPDF :: TimeWindowBP -> CalPDF+windowBP2CalPDF (TimeWindowBP name start stop) =+    windowBCAD2CalPDF (TimeWindowBCAD name (bp2BCAD start) (bp2BCAD stop))

@@ -110,15 +110,31 @@     , _calPDFDens :: VU.Vector Float     } deriving (Show, Eq) +-- | A data type for named calibration expressions+data NamedCalExpr = NamedCalExpr {+    -- | Expression identifier+      _exprID :: String+    -- | Expression+    , _expr   :: CalExpr+    } deriving (Show, Eq)+ -- | A data type to represent an expression for sum- or product calibration data CalExpr =       UnCalDate UncalC14+    | WindowBP TimeWindowBP+    | WindowBCAD TimeWindowBCAD     | CalDate CalPDF     | SumCal CalExpr CalExpr     | ProductCal CalExpr CalExpr-    deriving Show+    deriving (Show, Eq) -- http://www.cse.chalmers.se/edu/year/2018/course/TDA452/lectures/RecursiveDataTypes.html +data TimeWindowBP = TimeWindowBP String YearBP YearBP+    deriving (Show, Eq)++data TimeWindowBCAD = TimeWindowBCAD String YearBCAD YearBCAD+    deriving (Show, Eq)+ -- | A data type to represent a human readable summary of a calibrated radiocarbon date data CalC14 = CalC14 {     -- | Identifier, e.g. a lab number@@ -148,7 +164,7 @@  -- | 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--- cummulative probability (e.g. of an calibrated radiocarbon date density curve)+-- cumulative probability (e.g. of an calibrated radiocarbon date density curve) -- is above a certain threshold data HDR = HDR {     -- | Start of the high density region in years calBCAD@@ -156,3 +172,11 @@     -- | End of the high density region in years calBCAD     , _hdrstop  :: YearBCAD     } deriving (Show, Eq)++-- | A data type to store random samples drawn from a calPDF+data RandomAgeSample = RandomAgeSample {+    -- | Identifier+      _rasId      :: String+    -- | Random samples+    , _rasSamples :: VU.Vector YearBCAD+    } deriving Show

@@ -7,9 +7,15 @@  -- | Different exceptions for currycarbon data CurrycarbonException =-      CurrycarbonCLIParsingException String -- ^ An exception to describe an issue in the currycarbon CLI input parsing-    | CurrycarbonCalibrationRangeException String -- ^ An exection to describe the case that a-                                                  -- date is not in the range of the supplied calibration curve+    -- | An exception to describe an issue in the currycarbon CLI input parsing+      CurrycarbonCLIParsingException String+    -- | An exception to describe the case that a date is not in the range of+    -- the supplied calibration curve+    | CurrycarbonCalibrationRangeException String+    -- | An exception for CalPDFs that are unsuitable for certain purposes+    | CurrycarbonInvalidCalPDFException String+    -- | An exception for any issues with the CLI+    | CurrycarbonCLIException String     deriving (Show)  instance Exception CurrycarbonException@@ -17,5 +23,11 @@ renderCurrycarbonException :: CurrycarbonException -> String renderCurrycarbonException (CurrycarbonCLIParsingException s) =     "<!> Error: Input can not be parsed\n" ++ s-renderCurrycarbonException (CurrycarbonCalibrationRangeException s) =-    s ++ " <!> Error: Date outside of calibration range"+renderCurrycarbonException (CurrycarbonCalibrationRangeException i) =+    "<!> Error: Date outside of calibration range. Date ID: " ++ i+renderCurrycarbonException (CurrycarbonInvalidCalPDFException o) =+    "<!> Error: Invalid CalPDF for " ++ o +++    ", either because all densities are 0 or one density is > 1"+renderCurrycarbonException (CurrycarbonCLIException s) =+    "<!> Error: " ++ s+

@@ -0,0 +1,62 @@+module GoldenSpec (spec) where++import           Control.Applicative+import           Control.Monad+import           System.IO+import           System.Process+import           Test.Hspec          (Spec, describe, it, shouldBe)++spec :: Spec+spec = goldenTest++goldenTest :: Spec+goldenTest =++  describe "currycarbon cli test" $ do++    let stdout_stderr_tests = [+              "single_radiocarbon_date"+            , "cli_error_behaviour"+            ]++    runStdoutStderrTests stdout_stderr_tests++    let file_output_tests = [+              "density_file"+            , "hdr_file"+            , "samples_file"+            , "cal_curve_seg_file"+            ]++    runFileOutputTests file_output_tests++runStdoutStderrTests :: [String] -> Spec+runStdoutStderrTests tests = do+    forM_ tests $ \test -> do+        it (test ++ " should yield the correct stdout stderr output") $ do+            let cp = (shell ("bash " ++ test ++ ".sh")) {+                  cwd = Just "test/golden",+                  std_out = CreatePipe,+                  std_err = CreatePipe+                }+            (_, Just out, Just err, _) <- createProcess cp+            hSetBuffering out NoBuffering+            hSetBuffering err NoBuffering+            outActually <- liftA2 (++) (hGetContents err) (hGetContents out)+            outExpected <- readFile $ "test/golden/expected_data/" ++ test ++ ".out"+            outActually `shouldBe` outExpected++runFileOutputTests :: [String] -> Spec+runFileOutputTests tests = do+    forM_ tests $ \(test) -> do+        it (test ++ " should produce the correct output file") $ do+            let cp = (shell ("bash " ++ test ++ ".sh")) {+                  cwd = Just "test/golden",+                  std_out = CreatePipe,+                  std_err = CreatePipe+                }+            (_, _, _, exitCode) <- createProcess cp+            _ <- waitForProcess exitCode+            outActually <- readFile ("test/golden/actual_data/" ++ test ++ ".tsv")+            outExpected <- readFile ("test/golden/expected_data/" ++ test ++ ".tsv")+            outActually `shouldBe` outExpected

@@ -0,0 +1,117 @@+module ParserSpec (spec) where++import           Currycarbon.Parsers+import           Currycarbon.Types++import           Test.Hspec          (Spec, describe, it, shouldBe)++spec :: Spec+spec = do+  testReadNamedExpression++uncalC14N :: String -> CalExpr+uncalC14N s = UnCalDate (UncalC14 s 3000 30)+windowBPN :: String -> CalExpr+windowBPN s = WindowBP (TimeWindowBP s 3000 2000)+windowBCADN :: String -> CalExpr+windowBCADN s = WindowBCAD (TimeWindowBCAD s (-1050) (-50))+uncalC14 :: CalExpr+uncalC14 = uncalC14N ""+windowBP :: CalExpr+windowBP = windowBPN ""+windowBCAD :: CalExpr+windowBCAD = windowBCADN ""++testReadNamedExpression :: Spec+testReadNamedExpression =+    describe "Currycarbon.Parsers.readOneNamedCalExpr" $ do+    it "should read uncalibrated C14 dates correctly" $ do+        readOneNamedCalExpr "3000,30"+            `shouldBe`+            Right (NamedCalExpr "" uncalC14)+        readOneNamedCalExpr "uncalC14(3000,30)"+            `shouldBe`+            Right (NamedCalExpr "" uncalC14)+        readOneNamedCalExpr "uncalC14(test,3000,30)"+            `shouldBe`+            Right (NamedCalExpr "" (uncalC14N "test"))+    it "should read named function arguments correctly" $ do+        readOneNamedCalExpr "uncalC14(id = test, yearBP = 3000, sigma = 30)"+            `shouldBe`+            Right (NamedCalExpr "" (uncalC14N "test"))+    it "should read partially named function arguments correctly" $ do+        readOneNamedCalExpr "uncalC14(3000,sigma=30)"+            `shouldBe`+            Right (NamedCalExpr "" uncalC14)+    it "should read time windows correctly" $ do+        readOneNamedCalExpr "rangeBP(3000,2000)"+            `shouldBe`+            Right (NamedCalExpr "" windowBP)+        readOneNamedCalExpr "rangeBCAD(-1050,-50)"+            `shouldBe`+            Right (NamedCalExpr "" windowBCAD)+        readOneNamedCalExpr "rangeBP(test,3000,2000)"+            `shouldBe`+            Right (NamedCalExpr "" (windowBPN "test"))+        readOneNamedCalExpr "rangeBCAD(test,-1050,-50)"+            `shouldBe`+            Right (NamedCalExpr "" (windowBCADN "test"))+    it "should read sums with + operator correctly " $ do+        readOneNamedCalExpr "uncalC14(3000,30) + rangeBP(3000,2000)"+            `shouldBe`+            Right (NamedCalExpr "" $ SumCal uncalC14 windowBP)+        readOneNamedCalExpr "uncalC14(3000,30) + rangeBP(3000,2000) + rangeBCAD(-1050,-50)"+            `shouldBe`+            Right (NamedCalExpr "" $ SumCal uncalC14 (SumCal windowBP windowBCAD))+        readOneNamedCalExpr "uncalC14(3000,30) + rangeBP(3000,2000) + rangeBCAD(-1050,-50) + uncalC14(3000,30)"+            `shouldBe`+            Right (NamedCalExpr "" $ SumCal uncalC14 (SumCal windowBP (SumCal windowBCAD uncalC14)))+    it "should read sums with sum() function and + operator correctly " $ do+        readOneNamedCalExpr "sum(uncalC14(3000,30), rangeBP(3000,2000))"+            `shouldBe`+            Right (NamedCalExpr "" $ SumCal uncalC14 windowBP)+        readOneNamedCalExpr "sum(uncalC14(3000,30), rangeBP(3000,2000)) + rangeBCAD(-1050,-50)"+            `shouldBe`+            Right (NamedCalExpr "" $ SumCal (SumCal uncalC14 windowBP) windowBCAD)+        readOneNamedCalExpr "uncalC14(3000,30) + sum(rangeBP(3000,2000), rangeBCAD(-1050,-50)) + uncalC14(3000,30)"+            `shouldBe`+            Right (NamedCalExpr "" $ SumCal uncalC14 (SumCal (SumCal windowBP windowBCAD) uncalC14))+    it "should read products with * operator correctly " $ do+        readOneNamedCalExpr "uncalC14(3000,30) * rangeBP(3000,2000)"+            `shouldBe`+            Right (NamedCalExpr "" $ ProductCal uncalC14 windowBP)+        readOneNamedCalExpr "uncalC14(3000,30) * rangeBP(3000,2000) * rangeBCAD(-1050,-50)"+            `shouldBe`+            Right (NamedCalExpr "" $ ProductCal uncalC14 (ProductCal windowBP windowBCAD))+        readOneNamedCalExpr "uncalC14(3000,30) * rangeBP(3000,2000) * rangeBCAD(-1050,-50) * uncalC14(3000,30)"+            `shouldBe`+            Right (NamedCalExpr "" $ ProductCal uncalC14 (ProductCal windowBP (ProductCal windowBCAD uncalC14)))+    it "should read products with product() function and * operator correctly " $ do+        readOneNamedCalExpr "product(uncalC14(3000,30), rangeBP(3000,2000))"+            `shouldBe`+            Right (NamedCalExpr "" $ ProductCal uncalC14 windowBP)+        readOneNamedCalExpr "product(uncalC14(3000,30), rangeBP(3000,2000)) * rangeBCAD(-1050,-50)"+            `shouldBe`+            Right (NamedCalExpr "" $ ProductCal (ProductCal uncalC14 windowBP) windowBCAD)+        readOneNamedCalExpr "uncalC14(3000,30) * product(rangeBP(3000,2000), rangeBCAD(-1050,-50)) * uncalC14(3000,30)"+            `shouldBe`+            Right (NamedCalExpr "" $ ProductCal uncalC14 (ProductCal (ProductCal windowBP windowBCAD) uncalC14))+    it "should understand parenthesis correctly" $ do+        readOneNamedCalExpr "(uncalC14(3000,30) + rangeBP(3000,2000)) * rangeBCAD(-1050,-50)"+            `shouldBe`+            Right (NamedCalExpr "" $ ProductCal (SumCal uncalC14 windowBP) windowBCAD)+    it "should read unnamed and named calibration expressions correctly" $ do+        readOneNamedCalExpr "test: 3000,30"+            `shouldBe`+            Right (NamedCalExpr "test" uncalC14)+        readOneNamedCalExpr "calExpr(test,3000,30)"+            `shouldBe`+            Right (NamedCalExpr "test" uncalC14)+        readOneNamedCalExpr "calExpr(3000,30)"+            `shouldBe`+            Right (NamedCalExpr "" uncalC14)+    it "should be able to handle complex, nested queries" $ do+        readOneNamedCalExpr "calExpr(id = test, sum(uncalC14(3000,30), product(rangeBP(3000,2000), rangeBCAD(-1050,-50))) * uncalC14(3000,30))"+            `shouldBe`+            Right (NamedCalExpr "test" $ ProductCal (SumCal uncalC14 (ProductCal windowBP windowBCAD)) uncalC14)+

@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}