{-# LANGUAGE Strict #-}

module Currycarbon.Calibration.Bchron (calibrateDateBchron) where

import           Currycarbon.Calibration.Utils
import           Currycarbon.Parsers
import           Currycarbon.Types
import           Currycarbon.Utils

import qualified Data.Vector.Unboxed           as VU

-- | Intercept calibration as implemented in the Bchron R package (see 'Bchron')
calibrateDateBchron :: CalibrationDistribution -> Bool -> Bool -> CalCurveBP -> UncalC14 -> Either CurrycarbonException CalPDF
calibrateDateBchron :: CalibrationDistribution
-> Bool
-> Bool
-> CalCurveBP
-> UncalC14
-> Either CurrycarbonException CalPDF
calibrateDateBchron CalibrationDistribution
distr Bool
allowOutside Bool
interpolate CalCurveBP
calCurve uncalC14 :: UncalC14
uncalC14@(UncalC14 String
name YearBP
age YearBP
ageSd) =
    if Bool -> Bool
not Bool
allowOutside Bool -> Bool -> Bool
&& CalCurveBP -> UncalC14 -> Bool
isOutsideRangeOfCalCurve CalCurveBP
calCurve UncalC14
uncalC14
    then forall a b. a -> Either a b
Left forall a b. (a -> b) -> a -> b
$ String -> CurrycarbonException
CurrycarbonCalibrationRangeException forall a b. (a -> b) -> a -> b
$ UncalC14 -> String
renderUncalC14 UncalC14
uncalC14
    else
        let rawCalCurveSegment :: CalCurveBP
rawCalCurveSegment = UncalC14 -> CalCurveBP -> CalCurveBP
getRelevantCalCurveSegment UncalC14
uncalC14 CalCurveBP
calCurve
            CalCurveBCAD Vector YearBCAD
cals Vector YearBCAD
mus Vector YearBP
tau1s = Bool -> CalCurveBP -> CalCurveBCAD
prepareCalCurveSegment Bool
interpolate CalCurveBP
rawCalCurveSegment
            ageFloat :: Float
ageFloat = -(forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBP
age)forall a. Num a => a -> a -> a
+Float
1950
            ageSd2 :: YearBP
ageSd2 = YearBP
ageSdforall a. Num a => a -> a -> a
*YearBP
ageSd
            ageSd2Float :: Float
ageSd2Float = forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBP
ageSd2
            musFloat :: Vector Float
musFloat = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBCAD
mus
            tau1sFloat :: Vector Float
tau1sFloat = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBP
tau1s
            dens :: Vector Float
dens = case CalibrationDistribution
distr of
                CalibrationDistribution
NormalDist ->
                    forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VU.zipWith (\Float
mu Float
tau1 -> Float -> Float -> Float -> Float
dnorm Float
0 Float
1 ((Float
ageFloat forall a. Num a => a -> a -> a
- Float
mu) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Float
ageSd2Float forall a. Num a => a -> a -> a
+ Float
tau1 forall a. Num a => a -> a -> a
* Float
tau1))) Vector Float
musFloat Vector Float
tau1sFloat
                StudentTDist Double
degreesOfFreedom ->
                    forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VU.zipWith (\Float
mu Float
tau1 -> Double -> Float -> Float
dt Double
degreesOfFreedom ((Float
ageFloat forall a. Num a => a -> a -> a
- Float
mu) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Float
ageSd2Float forall a. Num a => a -> a -> a
+ Float
tau1 forall a. Num a => a -> a -> a
* Float
tau1))) Vector Float
musFloat Vector Float
tau1sFloat
        in forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ CalPDF -> CalPDF
trimLowDensityEdgesCalPDF forall a b. (a -> b) -> a -> b
$ CalPDF -> CalPDF
normalizeCalPDF forall a b. (a -> b) -> a -> b
$ String -> Vector YearBCAD -> Vector Float -> CalPDF
CalPDF String
name Vector YearBCAD
cals Vector Float
dens