diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,11 @@ +- V 0.6.0.0: + - New input method: Read .tsv files in a Poseidon .janno file-inspired format. + - Simplified CLI input options: Either read calibration expressions directly from the command line, or read them from one expression file (now called `--inputExprFile`), or read them from one .tsv file (`--inputTSVFile`). Combining the input methods or reading from multiple files at once is not possible any more. + - Added a lower limit for the size of the calibration curve segment in `getRelevantCalCurveSegment`. + - Added a mechanism to collect and show calibration errors at the end of a sequence of dates. + - Separated command line output with visual separators of the form `--- 1 ---`, `--- 2 ---`, etc. + - Fixed the calibration curve command line plot: It was always set to an interpolated `intcal20`, no matter the context. Now it properly considers the input settings. + - Turned off command line plots for the `--noInterpolation` case; the plotting code assumes year-wise resolution and is not trivial to adjust. - V 0.5.0.0: Cleaned up parser error behaviour: - Parsing errors are now more helpful with better positional information, using a modified `showParsecErr` and a new `showParsecErrOneLine`. - Breaking change: The type of `readCalCurve` changed. It is now safe, returning `Either String CalCurveBP`. The old version is preserved in `readCalCurveUnsafe`. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -32,10 +32,11 @@ ``` ``` -currycarbon v0.4.0.2 (UTF-8) +currycarbon v0.6.0.0 (UTF-8) Method: Bchron {distribution = StudentTDist {ndf = 100.0}} Curve: IntCal20 Calibrating... +--- 1 --- CalEXPR: [1] Sample1:4990±30BP Calibrated: 3936BC >> 3794BC > 3757BC < 3662BC << 3654BC 1-sigma: 3794-3707BC, 3666-3662BC @@ -63,10 +64,14 @@ ──────────────── ─ ─────────── ────────────────── ────── Done. +--- +No breaking issues. ``` ``` -Usage: currycarbon [--version] [CalEXPRs] [-i|--inputFile FILE] +Usage: currycarbon [--version] + [[CalEXPRs] | (-i|--inputExprFile FILE) | + (-t|--inputTSVFile FILE) [--combStrat Sum | Product]] [--calCurve IntCal20 | SHCal20 | Marine20 | FILE] [--method DSL] [--allowOutside] [--noInterpolation] [--noTrimCalCurve] [--noTrimOutCalPDF] [-q|--quiet] @@ -120,14 +125,34 @@ 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. + -i,--inputExprFile 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: ,,. + each input expression. Note that syntactic sugar + allows to read simple radiocarbon dates from a + headless .csv file with one sample per row: ,,. + -t,--inputTSVFile FILE A tab-separated file with dates to calibrate. Can + read Poseidon .janno files, which inspire this input + format. The rows are transformed to calibration + expressions upon reading, following a priority order + based on the provided columns: + - (or ): + Mandatory identifier, used in calExpr(id = ). + - + : + Gets transformed to uncalC14(id = , yearBP = + , sigma = ). Each of these columns can be a + list column with multiple entries separated by ;. + Multiple uncalibrated dates are then combined with + the strategy set in --combStrat when they are + evaluated. + - : + Gets transformed to rangeBCAD(start = , stop = + ). + --combStrat Sum | Product + Strategy that should be applied to combine multiple + C14 ages (uncalC14) from one row when reading from + --inputTSVFile. (default: Product) --calCurve IntCal20 | SHCal20 | Marine20 | FILE Either one of the included calibration curves, or a file path to an calibration curve file in '.14c' @@ -204,11 +229,16 @@ `stack` allows to upload a release candidate with ``` -stack haddock --haddock-for-hackage -stack upload . --test-tarball --candidate --documentation --no-save-hackage-creds +stack upload . --test-tarball --candidate --no-save-hackage-creds ``` -using my Hackage credentials. +using my Hackage credentials. It can then be published at https://hackage.haskell.org/package/currycarbon/candidates + +The building of the haddock documentation can be tested with + +``` +stack haddock --haddock-for-hackage +``` #### Preparing a new stable release diff --git a/currycarbon.cabal b/currycarbon.cabal --- a/currycarbon.cabal +++ b/currycarbon.cabal @@ -1,5 +1,5 @@ name: currycarbon -version: 0.5.0.0 +version: 0.6.0.0 synopsis: A package for simple, fast radiocarbon calibration description: Radiocarbon calibration library and command line tool. homepage: https://github.com/nevrome/currycarbon @@ -30,6 +30,7 @@ Currycarbon.ParserHelpers Currycarbon.Parsers Currycarbon.SumCalibration + Currycarbon.TSV Currycarbon.Types Currycarbon.Utils hs-source-dirs: @@ -41,8 +42,12 @@ , vector >= 0.12 && < 0.14 , math-functions >= 0.3 && < 0.4 , MonadRandom >= 0.6 && < 1 - , random > 1.2 && < 1.4 - , file-embed > 0.0.11.1 && < 0.1.0.0 + , random >= 1.2 && < 1.4 + , file-embed >= 0.0.11.1 && < 0.1 + , cassava >= 0.5 && < 0.6 + , bytestring >= 0.10.8 && < 0.13 + , unordered-containers >= 0.2.7.1 && < 0.3 + , text >= 1.2.2.1 && < 2.2 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 @@ -1,9 +1,11 @@ {-# LANGUAGE OverloadedStrings #-} import Currycarbon.CalCurves -import Currycarbon.CLI.RunCalibrate (CalibrateOptions (..), +import Currycarbon.CLI.RunCalibrate (CalibrateInput (..), + CalibrateOptions (..), runCalibrate) import Currycarbon.Parsers +import Currycarbon.TSV import Currycarbon.Types import Currycarbon.Utils import Paths_currycarbon (version) @@ -57,8 +59,7 @@ optParser = CmdCalibrate <$> calibrateOptParser calibrateOptParser :: OP.Parser CalibrateOptions -calibrateOptParser = CalibrateOptions <$> optParseNamedCalExprString - <*> optParseNamedCalExprFromFile +calibrateOptParser = CalibrateOptions <$> optParseCalibrateInput <*> optParseCalCurveSelection <*> optParseCalibrationMethod <*> optParseAllowOutside @@ -79,6 +80,11 @@ -- $inputParsing -- -- These functions define and handle the CLI input arguments +optParseCalibrateInput :: OP.Parser CalibrateInput +optParseCalibrateInput = + (CalibrateExprs <$> optParseNamedCalExprString) + OP.<|> (CalibrateExprFile <$> optParseNamedCalExprFromFile) + OP.<|> (CalibrateTSVFile <$> optParseInputTSVFile <*> optParseCombStrat) optParseNamedCalExprString :: OP.Parser [NamedCalExpr] optParseNamedCalExprString = concat <$> OP.many (OP.argument (OP.eitherReader readNamedCalExprs) ( @@ -132,19 +138,53 @@ 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" <> +optParseNamedCalExprFromFile :: OP.Parser FilePath +optParseNamedCalExprFromFile = OP.strOption ( + OP.long "inputExprFile" <> OP.short 'i' <> OP.metavar "FILE" <> OP.help "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: \ \,,." - )) + ) + +optParseInputTSVFile :: OP.Parser FilePath +optParseInputTSVFile = OP.strOption ( + OP.long "inputTSVFile" <> + OP.short 't' <> + OP.metavar "FILE" <> + OP.helpDoc ( Just ( + s2d "A tab-separated file with dates to calibrate. Can read Poseidon .janno files, \ + \which inspire this input format. The rows are transformed to calibration expressions \ + \upon reading, following a priority order based on the provided columns:" + <> OH.hardline <> + s2d "- (or ):" + <> OH.hardline <> + s2d "Mandatory identifier, used in calExpr(id = )." + <> OH.hardline <> + s2d "- :" + <> OH.hardline <> + s2d "Gets transformed to uncalC14(id = , yearBP = , sigma = ). \ + \Each of these columns can be a list column with multiple entries separated by ;. \ + \Multiple uncalibrated dates are then combined with the strategy set in --combStrat \ + \when they are evaluated." + <> OH.hardline <> + s2d "- :" + <> OH.hardline <> + s2d "Gets transformed to rangeBCAD(start = , stop = )." + ))) + +optParseCombStrat :: OP.Parser CombinationStrategy +optParseCombStrat = OP.option (OP.eitherReader readCombinationStrategy) ( + OP.long "combStrat" <> + OP.metavar "Sum | Product" <> + OP.help "Strategy that should be applied to combine multiple C14 ages (uncalC14) \ + \from one row when reading from --inputTSVFile." <> + OP.value CombProduct <> + OP.showDefault + ) optParseCalCurveSelection :: OP.Parser CalCurveSelection optParseCalCurveSelection = OP.option (OP.eitherReader readCalCurveSelection) ( diff --git a/src/Currycarbon.hs b/src/Currycarbon.hs --- a/src/Currycarbon.hs +++ b/src/Currycarbon.hs @@ -5,17 +5,17 @@ module Currycarbon ( - -- * Calibration + -- * Basic calibration -- $calibration calibrateDates, -- ** Configuration CalibrationMethod (..), CalibrateDatesConf (..), defaultCalConf, - -- ** Input + -- ** Basic input UncalC14 (..), readUncalC14FromFile, - -- ** Output + -- ** Basic output CalPDF (..), writeCalPDFs, @@ -42,20 +42,30 @@ writeCalC14HDRs, renderCalDatePretty, + -- * Drawing random samples from CalPDFs + -- $randsamp + AgeSamplingConf (..), + sampleAgesFromCalPDF, + RandomAgeSample (..), + writeRandomAgeSamples, + -- * Sum (and product) calibration -- $sumcal evalCalExpr, + evalNamedCalExpr, CalExpr (..), addPDFs, multiplyPDFs, normalizeCalPDF, trimLowDensityEdgesCalPDF, - -- * Drawing random samples from CalPDFs - -- $randsamp - AgeSamplingConf (..), - sampleAgesFromCalPDF, - RandomAgeSample (..) + -- * Reading input from .tsv files + -- $tsv + TSV (..), + TSVRow (..), + readTSV, + tsv2NamedCalExprs, + CombinationStrategy (..) ) where import Currycarbon.CalCurves @@ -64,6 +74,7 @@ trimLowDensityEdgesCalPDF) import Currycarbon.Parsers import Currycarbon.SumCalibration +import Currycarbon.TSV import Currycarbon.Types {- $calibration @@ -128,6 +139,18 @@ result for a given sample. -} +{- $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. +-} + {- $sumcal Calculating the sum or product of two calibration curves is a common @@ -142,14 +165,12 @@ -} -{- $randsamp +{- $tsv -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 can read calibration expressions from a tab-separated (.tsv) +file format inspired by the [Poseidon](https://www.poseidon-adna.org) .janno file. +`tsv2NamedCalExprs` turns the 'TSV' input to 'NamedCalExpr's, prioritizing +C14 age information over "contextual" start and stop ranges, to hook it into +the calibration pipeline. -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,24 +1,27 @@ +{-# LANGUAGE LambdaCase #-} + module Currycarbon.CLI.RunCalibrate - (CalibrateOptions (..), runCalibrate) where + (CalibrateOptions (..), CalibrateInput (..), runCalibrate) where import Currycarbon.CalCurves import Currycarbon.Calibration.Calibration import Currycarbon.Calibration.Utils import Currycarbon.Parsers import Currycarbon.SumCalibration +import Currycarbon.TSV import Currycarbon.Types import Currycarbon.Utils import Control.Exception (throwIO) import Control.Monad (unless, when) import Data.Maybe (fromJust, isJust) +import System.Exit (exitFailure, exitSuccess) 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 :: [NamedCalExpr] -- ^ String listing the uncalibrated dates that should be calibrated - , _calibrateExprFiles :: [FilePath] -- ^ List of files with uncalibrated dates to be calibrated + _calibrateInput :: CalibrateInput , _calibrateCalCurveFile :: CalCurveSelection -- ^ Either a preloaded calibration curve or a path to a .14c file , _calibrateCalibrationMethod :: CalibrationMethod -- ^ Calibration algorithm that should be used , _calibrateAllowOutside :: Bool -- ^ Allow calibration to run outside of the range of the calibration curve @@ -35,11 +38,18 @@ , _calibrateCalCurveMatrixFile :: Maybe FilePath -- ^ Path to an output file } +-- | A data type with input options for the calibration +data CalibrateInput = + CalibrateExprs [NamedCalExpr] -- ^ String listing the uncalibrated dates in CalExpr format that should be calibrated + | CalibrateExprFile FilePath -- ^ A file with uncalibrated dates in CalExpr to be calibrated + | CalibrateTSVFile FilePath CombinationStrategy -- ^ A file with uncalibrated dates in .tsv format to be calibrated and + -- a strategy to combine multiple C14 ages + -- | Interface function to trigger calibration from the command line runCalibrate :: CalibrateOptions -> IO () runCalibrate ( CalibrateOptions - exprs exprFiles + input calCurveSelection method allowOutside noInterpolate noTrimCalCurve noTrimOutCalPDF quiet encoding basicFile densityFile hdrFile @@ -48,9 +58,17 @@ ) = do let ascii = encoding /= "UTF-8" -- compile dates - exprsFromFile <- mapM readNamedCalExprsFromFile exprFiles - let exprsCombined = exprs ++ concat exprsFromFile - exprsRenamed = replaceEmptyNames exprsCombined + (errOrExprs,_) <- case input of + CalibrateExprs exprs -> do + return (map Right exprs, Nothing) + CalibrateExprFile path -> do + exprs <- readNamedCalExprsFromFile path + return (map Right exprs, Nothing) + CalibrateTSVFile path combStrat -> do + tsv <- readTSV path + let exprs = tsv2NamedCalExprs combStrat tsv + return (exprs, Just tsv) + let exprsRenamed = replaceEmptyNames errOrExprs if null exprsRenamed then hPutStrLn stderr "Nothing to calibrate. See currycarbon -h for help" else do @@ -69,7 +87,8 @@ -- handle the special debug cases when (isJust calCurveSegmentFile || isJust calCurveMatrixFile) $ do case exprsRenamed of - [NamedCalExpr _ (UnCalDate uncal)] -> do + [Left e] -> throwIO e + [Right (NamedCalExpr _ (UnCalDate uncal))] -> do let calCurveSegment = prepareCalCurveSegment (not noInterpolate) $ getRelevantCalCurveSegment uncal calCurve when (isJust calCurveSegmentFile) $ @@ -83,7 +102,7 @@ \a single uncalibrated radiocarbon date." -- run calibration hPutStrLn stderr "Calibrating..." - let errorOrCalPDFs = map (evalNamedCalExpr method calConf calCurve) exprsRenamed + let errorOrCalPDFs = fmap (either Left (evalNamedCalExpr method calConf calCurve)) exprsRenamed -- trim output calRes = if not noTrimOutCalPDF then map (mapEither id trimLowDensityEdgesCalPDF) errorOrCalPDFs @@ -95,41 +114,67 @@ 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 calRes - where + handleExprs ascii True maybeRNG calCurve $ zip3 [1..] exprsRenamed calRes + -- final conclusion + let calibrationErrors = + [ (i, e) + | (i, Left e) <- zip [1..] calRes + ] + refinementErrors = + [ (i, e) + | (i, Right calPDF) <- zip [1..] calRes + , Left e <- [validateCalPDF calPDF] + ] + errors = calibrationErrors ++ refinementErrors + if null errors + then do + hPutStrLn stderr "---" + hPutStrLn stderr "No breaking issues." + exitSuccess + else do + hPutStrLn stderr "---" + hPutStrLn stderr "Failed expressions:" + mapM_ (uncurry printE) errors + exitFailure + + where -- 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)] + -> CalCurveBP + -> [(Integer, Either CurrycarbonException NamedCalExpr, Either CurrycarbonException CalPDF)] -> IO () handleExprs _ _ _ _ [] = hPutStrLn stderr "Done." -- first expression - handleExprs _ascii True calCurve maybeRNG (firstDate:otherDates) = + handleExprs _ascii True maybeRNG calCurve (firstDate:otherDates) = case firstDate of - (_, Left e) -> do - printE e - handleExprs _ascii True calCurve maybeRNG otherDates - (namedCalExpr, Right cPDF) -> do + (i, Right ex, Right cPDF) -> do + unless quiet $ hPutStrLn stderr $ "--- " ++ show i ++ " ---" let (sampleSeed, newRNG) = drawSeed maybeRNG - flexOut _ascii namedCalExpr cPDF sampleSeed writeCalPDF writeCalC14CalRangeSummary writeCalC14HDR writeRandomAgeSample - handleExprs _ascii False calCurve newRNG otherDates + flexOut _ascii calCurve i ex cPDF sampleSeed writeCalPDF writeCalC14CalRangeSummary writeCalC14HDR writeRandomAgeSample + handleExprs _ascii False newRNG calCurve otherDates + (i, _, Left e) -> do + printE i e + handleExprs _ascii True maybeRNG calCurve otherDates + _ -> error "can not happen" -- subsequent expression - handleExprs _ascii False calCurve maybeRNG (nextDate:otherDates) = + handleExprs _ascii False maybeRNG calCurve (nextDate:otherDates) = case nextDate of - (_, Left e) -> do - printE e - handleExprs _ascii False calCurve maybeRNG otherDates - (namedCalExpr, Right cPDF) -> do + (i, Right ex, Right cPDF) -> do + unless quiet $ hPutStrLn stderr $ "--- " ++ show i ++ " ---" let (sampleSeed, newRNG) = drawSeed maybeRNG - flexOut _ascii namedCalExpr cPDF sampleSeed appendCalPDF appendCalC14CalRangeSummary appendCalC14HDR appendRandomAgeSample - handleExprs _ascii False calCurve newRNG otherDates + flexOut _ascii calCurve i ex cPDF sampleSeed appendCalPDF appendCalC14CalRangeSummary appendCalC14HDR appendRandomAgeSample + handleExprs _ascii False newRNG calCurve otherDates + (i, _, Left e) -> do + printE i e + handleExprs _ascii False maybeRNG calCurve otherDates + _ -> error "can not happen" - printE :: CurrycarbonException -> IO () - printE e = hPutStrLn stderr $ renderCurrycarbonException e + printE :: Integer -> CurrycarbonException -> IO () + printE i e = hPutStrLn stderr $ show i ++ ": " ++ renderCurrycarbonException e drawSeed :: Maybe R.StdGen -> (Maybe Int, Maybe R.StdGen) drawSeed maybeRNG = (\x -> (fromIntegral . fst <$> x, snd <$> x)) (R.genWord32 <$> maybeRNG) @@ -137,6 +182,8 @@ -- flexible expression handler flexOut :: Bool + -> CalCurveBP + -> Integer -> NamedCalExpr -> CalPDF -> Maybe Int @@ -145,19 +192,19 @@ -> (FilePath -> CalC14 -> IO ()) -> (FilePath -> RandomAgeSample -> IO ()) -> IO () - flexOut _ascii namedCalExpr calPDF maybeSeed calPDFToFile calC14CalRangeSummaryToFile calC14HDRToFile randomAgeSampleToFile = do + flexOut _ascii calCurve i namedCalExpr calPDF maybeSeed calPDFToFile calC14CalRangeSummaryToFile calC14HDRToFile randomAgeSampleToFile = do case refineCalDate calPDF of Left e -> do unless quiet $ do putStrLn ("CalEXPR: " ++ renderNamedCalExpr namedCalExpr) - printE e + printE i e when (isJust basicFile) $ unless quiet $ hPutStrLn stderr " Error: Can not create --basicFile" when (isJust hdrFile) $ unless quiet $ hPutStrLn stderr " Error: Can not create --hdrFile" Right calC14 -> do unless quiet $ do - putStrLn (renderCalDatePretty _ascii (namedCalExpr, calPDF, calC14)) + putStrLn (renderCalDatePretty _ascii calCurve (not noInterpolate) (namedCalExpr, calPDF, calC14)) when (isJust basicFile) $ calC14CalRangeSummaryToFile (fromJust basicFile) calC14 when (isJust hdrFile) $ @@ -169,7 +216,7 @@ case sampleAgesFromCalPDF conf calPDF of Left e -> do unless quiet $ do - printE e + printE i e hPutStrLn stderr " Error: Can not create --samplesFile" Right res -> randomAgeSampleToFile path res when (isJust densityFile) $ @@ -177,9 +224,12 @@ -- | Helper function to replace empty input names with a sequence of numbers, -- to get each input date an unique identifier -replaceEmptyNames :: [NamedCalExpr] -> [NamedCalExpr] -replaceEmptyNames = zipWith (modifyNamedExpr . show) ([1..] :: [Integer]) +replaceEmptyNames :: [Either CurrycarbonException NamedCalExpr] -> [Either CurrycarbonException NamedCalExpr] +replaceEmptyNames = zipWith modify [1..] where + modify :: Integer -> Either CurrycarbonException NamedCalExpr -> Either CurrycarbonException NamedCalExpr + modify _ (Left e) = Left e + modify i (Right x) = Right $ modifyNamedExpr (show i) x modifyNamedExpr :: String -> NamedCalExpr -> NamedCalExpr modifyNamedExpr i nexpr = if _exprID nexpr == "" diff --git a/src/Currycarbon/CalCurves.hs b/src/Currycarbon/CalCurves.hs --- a/src/Currycarbon/CalCurves.hs +++ b/src/Currycarbon/CalCurves.hs @@ -109,3 +109,37 @@ marine20 = readCalCurveUnsafe marine20String marine20String :: String marine20String = $(FE.makeRelativeToProject "data/marine20.14c" >>= FE.embedStringFile) + +-- | Calibration curves usually do not have yearly resolution. +-- This function completes the curves with a most basic linear interpolation. +interpolateCalCurve :: CalCurveBP -> CalCurveBP +interpolateCalCurve (CalCurveBP cals uncals sigmas) = + let obs = VU.zip3 cals uncals sigmas + timeWindows = getTimeWindows obs + obsFilled = VU.concatMap fillTimeWindows timeWindows + in uncurry3 CalCurveBP $ VU.unzip3 obsFilled + where + getTimeWindows :: VU.Vector (YearBP,YearBP,YearRange) -> VU.Vector ((YearBP,YearBP,YearRange),(YearBP,YearBP,YearRange)) + getTimeWindows xs = VU.zipWith (,) (VU.init xs) (VU.tail xs) + fillTimeWindows :: ((YearBP,YearBP,YearRange),(YearBP,YearBP,YearRange)) -> VU.Vector (YearBP,YearBP,YearRange) + fillTimeWindows ((calbp1,bp1,sigma1),(calbp2,bp2,sigma2)) = + if calbp1 == calbp2 || calbp1+1 == calbp2 || calbp1-1 == calbp2 + then VU.singleton (calbp1,bp1,sigma1) + else + let newCals = VU.fromList [calbp1,calbp1-1..calbp2+1] -- range definition like this to trigger counting down + newBPs = VU.map (snd . getInBetweenPointsInt (calbp1,bp1) (calbp2,bp2)) newCals + newSigmas = VU.map (snd . getInBetweenPointsInt (calbp1,sigma1) (calbp2,sigma2)) newCals + in VU.zip3 newCals newBPs newSigmas + getInBetweenPointsInt :: (Word, Word) -> (Word, Word) -> Word -> (Word, Word) + getInBetweenPointsInt (x1,y1) (x2,y2) xPred = + let (_,yPred) = getInBetweenPoints (fromIntegral x1,fromIntegral y1) (fromIntegral x2,fromIntegral y2) $ fromIntegral xPred + in (xPred, round yPred) + getInBetweenPoints :: (Double, Double) -> (Double, Double) -> Double -> (Double, Double) + getInBetweenPoints (x1,y1) (x2,y2) xPred = + let yDiff = y2 - y1 + xDiff = abs $ x1 - x2 + yDiffPerxDiff = yDiff/xDiff + xPredRel = x1 - xPred + in (xPred, y1 + xPredRel * yDiffPerxDiff) + uncurry3 :: (a -> b -> c -> d) -> ((a, b, c) -> d) + uncurry3 f ~(a,b,c) = f a b c 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 @@ -14,6 +14,7 @@ , calibrateDates , refineCalDates , refineCalDate + , validateCalPDF , CalibrateDatesConf (..) , defaultCalConf , AgeSamplingConf (..) @@ -56,21 +57,28 @@ -- calibration failed for some reason, or a 'CalPDF' calibrateDate method config calCurve uncalDate = head $ calibrateDates method config calCurve [uncalDate] +validateCalPDF :: CalPDF -> Either CurrycarbonException () +validateCalPDF calPDF + | isInvalidCalPDF calPDF = Left $ CurrycarbonInvalidCalPDFException "refinement" + | otherwise = Right () + -- | Transforms the raw, calibrated probability density table to a meaningful representation of a -- calibrated radiocarbon date refineCalDates :: [CalPDF] -> [Either CurrycarbonException CalC14] refineCalDates = map refineCalDate 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" +refineCalDate calPDF = do + validateCalPDF calPDF + pure (refineValidCalDate calPDF) + +refineValidCalDate :: CalPDF -> CalC14 +refineValidCalDate (CalPDF name cals dens) -- for simple uniform age ranges | VU.length (VU.uniq dens) == 1 = let start = VU.head cals stop = VU.last cals - in Right $ CalC14 { + in CalC14 { _calC14id = name , _calC14RangeSummary = CalRangeSummary { _calRangeStartTwoSigma = start @@ -84,7 +92,7 @@ } -- for normal post-calibration probability distributions | otherwise = - Right $ CalC14 { + CalC14 { _calC14id = name , _calC14RangeSummary = CalRangeSummary { _calRangeStartTwoSigma = _hdrstart $ head hdrs95 @@ -156,4 +164,4 @@ where m = sequence . repeat . CMR.fromList $ weights isInvalidCalPDF :: CalPDF -> Bool -isInvalidCalPDF (CalPDF _ _ dens) = VU.sum dens == 0 || VU.any (>= 1.0) dens +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 @@ -2,6 +2,7 @@ module Currycarbon.Calibration.Utils where +import Currycarbon.CalCurves (interpolateCalCurve) import Currycarbon.Types import Data.Maybe (fromMaybe) @@ -51,8 +52,9 @@ -- the relevant segment of the calibration curve getRelevantCalCurveSegment :: UncalC14 -> CalCurveBP -> CalCurveBP getRelevantCalCurveSegment (UncalC14 _ mean std) (CalCurveBP cals uncals sigmas) = - let start = mean+6*std - stop = mean-6*std + let std' = max std 10 + start = mean+6*std' + stop = mean-6*std' startIndex = fromMaybe 0 $ VU.findIndex (<= start) uncals stopIndex = (VU.length uncals - 1) - fromMaybe 0 (VU.findIndex (>= stop) $ VU.reverse uncals) toIndex = stopIndex - startIndex @@ -86,38 +88,6 @@ bcad2BP :: YearBCAD -> YearBP bcad2BP y = 1950 - fromIntegral y - -interpolateCalCurve :: CalCurveBP -> CalCurveBP -interpolateCalCurve (CalCurveBP cals uncals sigmas) = - let obs = VU.zip3 cals uncals sigmas - timeWindows = getTimeWindows obs - obsFilled = VU.concatMap fillTimeWindows timeWindows - in uncurry3 CalCurveBP $ VU.unzip3 obsFilled - where - getTimeWindows :: VU.Vector (YearBP,YearBP,YearRange) -> VU.Vector ((YearBP,YearBP,YearRange),(YearBP,YearBP,YearRange)) - getTimeWindows xs = VU.zipWith (,) (VU.init xs) (VU.tail xs) - fillTimeWindows :: ((YearBP,YearBP,YearRange),(YearBP,YearBP,YearRange)) -> VU.Vector (YearBP,YearBP,YearRange) - fillTimeWindows ((calbp1,bp1,sigma1),(calbp2,bp2,sigma2)) = - if calbp1 == calbp2 || calbp1+1 == calbp2 || calbp1-1 == calbp2 - then VU.singleton (calbp1,bp1,sigma1) - else - let newCals = VU.fromList [calbp1,calbp1-1..calbp2+1] -- range definition like this to trigger counting down - newBPs = VU.map (snd . getInBetweenPointsInt (calbp1,bp1) (calbp2,bp2)) newCals - newSigmas = VU.map (snd . getInBetweenPointsInt (calbp1,sigma1) (calbp2,sigma2)) newCals - in VU.zip3 newCals newBPs newSigmas - getInBetweenPointsInt :: (Word, Word) -> (Word, Word) -> Word -> (Word, Word) - getInBetweenPointsInt (x1,y1) (x2,y2) xPred = - let (_,yPred) = getInBetweenPoints (fromIntegral x1,fromIntegral y1) (fromIntegral x2,fromIntegral y2) $ fromIntegral xPred - in (xPred, round yPred) - getInBetweenPoints :: (Double, Double) -> (Double, Double) -> Double -> (Double, Double) - getInBetweenPoints (x1,y1) (x2,y2) xPred = - let yDiff = y2 - y1 - xDiff = abs $ x1 - x2 - yDiffPerxDiff = yDiff/xDiff - xPredRel = x1 - xPred - in (xPred, y1 + xPredRel * yDiffPerxDiff) - uncurry3 :: (a -> b -> c -> d) -> ((a, b, c) -> d) - uncurry3 f ~(a,b,c) = f a b c -- | Subset the calibration curve to the non-zero density range. The threshold is set to 0.00001. trimLowDensityEdgesCalPDF :: CalPDF -> CalPDF diff --git a/src/Currycarbon/Parsers.hs b/src/Currycarbon/Parsers.hs --- a/src/Currycarbon/Parsers.hs +++ b/src/Currycarbon/Parsers.hs @@ -2,7 +2,7 @@ module Currycarbon.Parsers where -import Currycarbon.CalCurves (intcal20) +import Currycarbon.CalCurves (interpolateCalCurve) import Currycarbon.Calibration.Utils import Currycarbon.ParserHelpers import Currycarbon.Types @@ -75,16 +75,25 @@ -- renderCalDatePretty :: Bool -- ^ Should the CLI plot be restricted to (boring) ASCII symbols? + -> CalCurveBP + -> Bool -> (NamedCalExpr, CalPDF, CalC14) -> String -renderCalDatePretty ascii (calExpr, calPDF, calC14) = - "CalEXPR: " ++ intercalate "\n" [ +renderCalDatePretty ascii calCurve interpolate (calExpr, calPDF, calC14) = + "CalEXPR: " ++ intercalate "" [ renderNamedCalExpr calExpr + , "\n" , renderCalC14 calC14 - , "" - , renderCLIPlotCalCurve ascii 8 50 calPDF calExpr - , renderCLIPlotCalPDF ascii 6 50 calPDF calC14 - ] + , "\n" + ] ++ + if interpolate + then intercalate "" [ + "\n" + , renderCLIPlotCalCurve ascii 8 50 calCurve calPDF calExpr + , "\n" + , renderCLIPlotCalPDF ascii 6 50 calPDF calC14 + ] + else [] -- write and read calibration expressions @@ -117,15 +126,11 @@ name <- parseArgument "id" parseAnyString start <- parseArgument "start" parseWord stop <- parseArgument "stop" parseWord - construct name start stop + eitherToFail $ makeTimeWindowBP 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" + eitherToFail $ makeTimeWindowBP "" start stop parseTimeWindowBCAD :: P.Parser TimeWindowBCAD parseTimeWindowBCAD = parseRecordType "rangeBCAD" $ P.try long P.<|> short @@ -134,15 +139,11 @@ name <- parseArgument "id" parseAnyString start <- parseArgument "start" parseInt stop <- parseArgument "stop" parseInt - construct name start stop + eitherToFail $ makeTimeWindowBCAD 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" + eitherToFail $ makeTimeWindowBCAD "" start stop -- https://gist.github.com/abhin4v/017a36477204a1d57745 addFun :: P.Parser CalExpr @@ -500,14 +501,16 @@ roundedDec = if rest >= 5 then dec + 1 else dec in roundedDec * 10 * signum x -renderCLIPlotCalCurve :: Bool -> Int -> Int -> CalPDF -> NamedCalExpr -> String +renderCLIPlotCalCurve :: Bool -> Int -> Int -> CalCurveBP -> CalPDF -> NamedCalExpr -> String renderCLIPlotCalCurve - ascii rows cols (CalPDF _ cals _) + ascii rows cols + calCurve + (CalPDF _ cals _) (NamedCalExpr _ (UnCalDate (UncalC14 _ yearBP sigma))) = let startYear = VU.head cals stopYear = VU.last cals -- prepare calcurve - calcurvePrep = makeBCADCalCurve $ interpolateCalCurve intcal20 + calcurvePrep = makeBCADCalCurve $ interpolateCalCurve calCurve calCurveSegment = punchOutCalCurveBCAD startYear stopYear calcurvePrep calCurveUncals = VU.map fromIntegral $ _calCurveBCADUnCals calCurveSegment calCurveUncalStart = bcad2BP $ round $ VU.head calCurveUncals @@ -553,7 +556,7 @@ | otherwise = ' ' makeTick :: (Integral n) => n -> String makeTick n = padString 6 (show $ roundTo10 $ fromIntegral n) ++ " " ++ getSymbol ascii YAxisTick : " " -renderCLIPlotCalCurve _ _ _ _ _ = "" +renderCLIPlotCalCurve _ _ _ _ _ _ = "" renderCLIPlotCalPDF :: Bool -> Int -> Int -> CalPDF -> CalC14 -> String renderCLIPlotCalPDF ascii rows cols (CalPDF _ cals dens) c14 = diff --git a/src/Currycarbon/TSV.hs b/src/Currycarbon/TSV.hs new file mode 100644 --- /dev/null +++ b/src/Currycarbon/TSV.hs @@ -0,0 +1,171 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Currycarbon.TSV where + +import Currycarbon.ParserHelpers +import Currycarbon.Types +import Currycarbon.Utils + +import Control.Applicative (empty) +import Control.Exception (throwIO) +import qualified Data.ByteString.Char8 as B8 +import qualified Data.ByteString.Lazy as BL +import Data.Char (ord) +import qualified Data.Csv as Csv +import qualified Data.HashMap.Strict as HM +import Data.Maybe (catMaybes) +import qualified Data.Vector as V +import qualified Text.Parsec as P +import qualified Text.Parsec.String as P + +-- | Turn 'TSV' to 'NamedCalExpr's for further processing (see 'evalNamedCalExpr') +tsv2NamedCalExprs :: CombinationStrategy -> TSV -> [Either CurrycarbonException NamedCalExpr] +tsv2NamedCalExprs combStrat (TSV _ _ rows) = V.toList $ V.map (tsvRow2NamedCalExprs combStrat) rows + +tsvRow2NamedCalExprs :: CombinationStrategy -> TSVRow -> Either CurrycarbonException NamedCalExpr +-- C14 age with n dates and a labcode +tsvRow2NamedCalExprs combStrat + (TSVRow i (Just (ListColumn lcs)) (Just (ListColumn bps)) (Just (ListColumn errs)) _ _) | + (length lcs == length bps) && (length bps == length errs) = + Right $ NamedCalExpr i $ foldC14 combStrat $ zip3 lcs bps errs +-- C14 age with n dates and no labcode +tsvRow2NamedCalExprs combStrat + (TSVRow i _ (Just (ListColumn bps)) (Just (ListColumn errs)) _ _) | + length bps == length errs = + Right $ NamedCalExpr i $ foldC14 combStrat $ zip3 (repeat "") bps errs +-- contextual age +tsvRow2NamedCalExprs _ (TSVRow i _ _ _ (Just start) (Just stop)) = + NamedCalExpr i . WindowBCAD <$> makeTimeWindowBCAD "" start stop +-- error case if nothing fits +tsvRow2NamedCalExprs _ (TSVRow i _ _ _ _ _) = + Left $ CurrycarbonTSV2CalExprException i + +-- | A data type for the strategy to combine multiple C14 dates. +data CombinationStrategy = + CombSum -- ^ Form a normalised sum of the 'CalPDF's. + | CombProduct -- ^ Form a normalised product of the 'CalPDF's. + +instance Show CombinationStrategy where + show CombSum = "Sum" + show CombProduct = "Product" + +readCombinationStrategy :: String -> Either String CombinationStrategy +readCombinationStrategy s = + case P.runParser parseCombinationStrategy () s s of + Left err -> Left $ showParsecErrOneLine err + Right x -> Right x + +parseCombinationStrategy :: P.Parser CombinationStrategy +parseCombinationStrategy = do + x <- P.many P.anyChar + case x of + "Sum" -> pure CombSum + "Product" -> pure CombProduct + _ -> fail "must be either Sum or Product" + +foldC14 :: CombinationStrategy -> [(String, Word, Word)] -> CalExpr +foldC14 CombSum xs = foldl1 SumCal $ map (\(lc,bp,err) -> UnCalDate $ UncalC14 lc bp err) xs +foldC14 CombProduct xs = foldl1 ProductCal $ map (\(lc,bp,err) -> UnCalDate $ UncalC14 lc bp err) xs + +-- | Read calibration expressions from a .tsv file. The file should have the following columns: +-- +-- @ +-- Date_ID Date_C14_Labnr Date_C14_Uncal_BP Date_C14_Uncal_BP_Err Date_BC_AD_Start Date_BC_AD_Stop +-- Sample1 -200 100 +-- Sample2 TEST-1 1000 30 +-- Sample3 TEST-2;TEST-3 3000;3200 30;50 +-- @ +-- +-- Subsets like + + or +-- + + are also sufficient. +-- The structure is inspired by the Poseidon .janno file. +-- +readTSV :: FilePath -> IO TSV +readTSV path = do + bs <- BL.readFile path + case Csv.decodeByNameWith decodingOptions bs of + Left s -> throwIO $ CurrycarbonTSVParsingException s + Right (header, rows) -> return (TSV path header rows) + +decodingOptions :: Csv.DecodeOptions +decodingOptions = Csv.defaultDecodeOptions { Csv.decDelimiter = fromIntegral (ord '\t') } + +-- | A data type for tab-separated (.tsv) files +data TSV = TSV { + _tsvFile :: FilePath + , _tsvHeader :: Csv.Header + , _tsvRows :: V.Vector TSVRow + } + +-- | A data type to represent the relevant fields of one row in a .tsv file +data TSVRow = TSVRow { + -- | Identifier of the expression + _tsvRowID :: String + -- | Lab codes of radiocarbon dates + , _tsvRowDateC14Labnr :: Maybe (ListColumn String) + -- | C14 ages in years BP + , _tsvRowDateC14UncalBP :: Maybe (ListColumn Word) + -- | C14 standard deviations (one sigma in years) + , _tsvRowDateC14UncalBPErr :: Maybe (ListColumn Word) + -- | Start of time window in years BC + , _tsvRowDateBCADStart :: Maybe Int + -- | End of time window in years BC + , _tsvRowDateBCADStop :: Maybe Int + -- for any other columns + -- , _tsvRowAllColumns :: Csv.NamedRecord + } + deriving Show + +instance Csv.FromNamedRecord TSVRow where + parseNamedRecord m = do + i <- filterLookupMulti m ["Date_ID", "Poseidon_ID"] + labnr <- filterLookupOptional m "Date_C14_Labnr" + uncalBP <- filterLookupOptional m "Date_C14_Uncal_BP" + uncalBPErr <- filterLookupOptional m "Date_C14_Uncal_BP_Err" + start <- filterLookupOptional m "Date_BC_AD_Start" + stop <- filterLookupOptional m "Date_BC_AD_Stop" + pure $ TSVRow { + _tsvRowID = i + , _tsvRowDateC14Labnr = labnr + , _tsvRowDateC14UncalBP = uncalBP + , _tsvRowDateC14UncalBPErr = uncalBPErr + , _tsvRowDateBCADStart = start + , _tsvRowDateBCADStop = stop + -- for any other columns + -- , _tsvRowAllColumns = m + } + +-- lookup column by name +filterLookup :: Csv.FromField a => Csv.NamedRecord -> B8.ByteString -> Csv.Parser a +filterLookup m name = maybe empty Csv.parseField $ cleanInput $ HM.lookup name m + +-- lookup optional column by name +filterLookupOptional :: Csv.FromField a => Csv.NamedRecord -> B8.ByteString -> Csv.Parser (Maybe a) +filterLookupOptional m name = maybe (pure Nothing) (\bs -> Just <$> Csv.parseField bs) $ + cleanInput $ HM.lookup name m + +-- lookup column by multiple different names and keep the first match +filterLookupMulti :: Csv.FromField a => Csv.NamedRecord -> [B8.ByteString] -> Csv.Parser a +filterLookupMulti m names = + maybe empty Csv.parseField $ cleanInput $ lookupMulti names + where + lookupMulti :: [B8.ByteString] -> Maybe B8.ByteString + lookupMulti ns = + let vals = map (`HM.lookup` m) ns + in case catMaybes vals of + [] -> Nothing + (x:_) -> Just x + +cleanInput :: Maybe B8.ByteString -> Maybe B8.ByteString +cleanInput Nothing = Nothing +cleanInput (Just rawInputBS) = transNA rawInputBS + where + transNA ".." = Nothing + transNA "" = Nothing + transNA "n/a" = Nothing + transNA x = Just x + +newtype ListColumn a = ListColumn {getListColumn :: [a]} + deriving (Eq, Ord, Show) +instance (Csv.FromField a) => Csv.FromField (ListColumn a) where + parseField x = fmap ListColumn . mapM Csv.parseField $ B8.splitWith (==';') x diff --git a/src/Currycarbon/Types.hs b/src/Currycarbon/Types.hs --- a/src/Currycarbon/Types.hs +++ b/src/Currycarbon/Types.hs @@ -2,6 +2,8 @@ module Currycarbon.Types where +import Currycarbon.Utils + import qualified Data.Vector as V import qualified Data.Vector.Unboxed as VU @@ -164,8 +166,22 @@ data TimeWindowBP = TimeWindowBP String YearBP YearBP deriving (Show, Eq) +makeTimeWindowBP :: String -> YearBP -> YearBP -> Either CurrycarbonException TimeWindowBP +makeTimeWindowBP n start stop + | start < stop = Left $ CurrycarbonMakeCalExprException "the BP start age can not be younger than stop age" + | otherwise = Right $ TimeWindowBP n start stop + data TimeWindowBCAD = TimeWindowBCAD String YearBCAD YearBCAD deriving (Show, Eq) + +makeTimeWindowBCAD :: String -> YearBCAD -> YearBCAD -> Either CurrycarbonException TimeWindowBCAD +makeTimeWindowBCAD n start stop + | start > stop = Left $ CurrycarbonMakeCalExprException "the BC/AD start age can not be younger than the stop age" + | otherwise = Right $ TimeWindowBCAD n start stop + +eitherToFail :: MonadFail m => Either CurrycarbonException b -> m b +eitherToFail (Left e) = fail $ renderCurrycarbonException e +eitherToFail (Right x) = return x -- | A data type to represent a human readable summary of a calibrated radiocarbon date data CalC14 = CalC14 { diff --git a/src/Currycarbon/Utils.hs b/src/Currycarbon/Utils.hs --- a/src/Currycarbon/Utils.hs +++ b/src/Currycarbon/Utils.hs @@ -16,6 +16,12 @@ | CurrycarbonInvalidCalPDFException String -- | An exception for any issues with the CLI | CurrycarbonCLIException String + -- | An exception for issues with the .tsv file parsing + | CurrycarbonTSVParsingException String + -- | An exception for incomplete rows in .tsv files + | CurrycarbonTSV2CalExprException String + -- | An exception for issues when making CalEXPRs + | CurrycarbonMakeCalExprException String deriving (Show) instance Exception CurrycarbonException @@ -30,4 +36,9 @@ ", either because all densities are 0 or one density is > 1" renderCurrycarbonException (CurrycarbonCLIException s) = " Error: " ++ s - +renderCurrycarbonException (CurrycarbonTSVParsingException s) = + " Error: " ++ s +renderCurrycarbonException (CurrycarbonTSV2CalExprException i) = + " Error: Incomplete information in .tsv entry " ++ show i +renderCurrycarbonException (CurrycarbonMakeCalExprException s) = + " Error: Issue when making CalExpr: " ++ s diff --git a/test/GoldenSpec.hs b/test/GoldenSpec.hs --- a/test/GoldenSpec.hs +++ b/test/GoldenSpec.hs @@ -25,6 +25,7 @@ , "hdr_file" , "samples_file" , "cal_curve_seg_file" + , "tsv_file_reading" ] runFileOutputTests file_output_tests