diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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 ",," where is optional - (e.g. "S1,4000,50"). Multiple dates can be listed - separated by ";" (e.g. "S1,4000,50; 3000,25; - S3,1000,20"). To sum or multiply the post calibration - probability distributions, dates can be combined with - "+" or "*" (e.g. "4000,50 + 4100,100"). These - expressions can be combined arbitrarily. Parentheses - can be added to specify the order of operations (e.g. - "(4000,50 + 4100,100) * 3800,50") - -i,--inputFile ARG A file with a list of calibration expressions. - Formated just as DATE, but with a new line for each - input date. DATE and --inputFile can be combined and - you can provide multiple instances of --inputFile - --calibrationCurveFile ARG - Path to an calibration curve file in .14c format. The - calibration curve will be read and used for + 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 ' =' 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: ,,. + --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: - ",,". + use the 'intcal20' curve. + --method DSL The calibration algorithm that should be used: + ',,'. 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 diff --git a/currycarbon.cabal b/currycarbon.cabal --- a/currycarbon.cabal +++ b/currycarbon.cabal @@ -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 diff --git a/src-executables/Main-currycarbon.hs b/src-executables/Main-currycarbon.hs --- a/src-executables/Main-currycarbon.hs +++ b/src-executables/Main-currycarbon.hs @@ -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 \",,\" \ - \where is optional (e.g. \"S1,4000,50\"). \ - \Multiple dates can be listed separated by \";\" (e.g. \"S1,4000,50; 3000,25; S3,1000,20\"). \ - \To sum or multiply the post calibration probability distributions, dates can be combined with \ - \\"+\" or \"*\" (e.g. \"4000,50 + 4100,100\"). \ - \These expressions can be combined arbitrarily. Parentheses can be added to specify the order \ - \of operations (e.g. \"(4000,50 + 4100,100) * 3800,50\")" +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 ' =' \ + \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: \ + \,,." )) 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: \ - \\",,\". \ - \The default setting is equivalent to \"Bchron,StudentT,100\" \ + OP.metavar "DSL" <> + OP.helpDoc ( Just ( + s2d "The calibration algorithm that should be used: \ + \',,'. " + <> 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 ) diff --git a/src/Currycarbon.hs b/src/Currycarbon.hs --- a/src/Currycarbon.hs +++ b/src/Currycarbon.hs @@ -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. -} diff --git a/src/Currycarbon/CLI/RunCalibrate.hs b/src/Currycarbon/CLI/RunCalibrate.hs --- a/src/Currycarbon/CLI/RunCalibrate.hs +++ b/src/Currycarbon/CLI/RunCalibrate.hs @@ -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) diff --git a/src/Currycarbon/Calibration/Calibration.hs b/src/Currycarbon/Calibration/Calibration.hs --- a/src/Currycarbon/Calibration/Calibration.hs +++ b/src/Currycarbon/Calibration/Calibration.hs @@ -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 diff --git a/src/Currycarbon/Calibration/Utils.hs b/src/Currycarbon/Calibration/Utils.hs --- a/src/Currycarbon/Calibration/Utils.hs +++ b/src/Currycarbon/Calibration/Utils.hs @@ -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) = diff --git a/src/Currycarbon/ParserHelpers.hs b/src/Currycarbon/ParserHelpers.hs new file mode 100644 --- /dev/null +++ b/src/Currycarbon/ParserHelpers.hs @@ -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) diff --git a/src/Currycarbon/Parsers.hs b/src/Currycarbon/Parsers.hs --- a/src/Currycarbon/Parsers.hs +++ b/src/Currycarbon/Parsers.hs @@ -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 "\,\,\", where +-- \ 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 "\,\,\", where --- \ 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" diff --git a/src/Currycarbon/SumCalibration.hs b/src/Currycarbon/SumCalibration.hs --- a/src/Currycarbon/SumCalibration.hs +++ b/src/Currycarbon/SumCalibration.hs @@ -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)) diff --git a/src/Currycarbon/Types.hs b/src/Currycarbon/Types.hs --- a/src/Currycarbon/Types.hs +++ b/src/Currycarbon/Types.hs @@ -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 diff --git a/src/Currycarbon/Utils.hs b/src/Currycarbon/Utils.hs --- a/src/Currycarbon/Utils.hs +++ b/src/Currycarbon/Utils.hs @@ -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 + diff --git a/test/GoldenSpec.hs b/test/GoldenSpec.hs new file mode 100644 --- /dev/null +++ b/test/GoldenSpec.hs @@ -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 diff --git a/test/ParserSpec.hs b/test/ParserSpec.hs new file mode 100644 --- /dev/null +++ b/test/ParserSpec.hs @@ -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) + diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,1 @@ +{-# OPTIONS_GHC -F -pgmF hspec-discover #-}