{-# 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 CurrycarbonException -> Either CurrycarbonException CalPDF
forall a b. a -> Either a b
Left (CurrycarbonException -> Either CurrycarbonException CalPDF)
-> CurrycarbonException -> Either CurrycarbonException CalPDF
forall a b. (a -> b) -> a -> b
$ String -> CurrycarbonException
CurrycarbonCalibrationRangeException (String -> CurrycarbonException) -> String -> CurrycarbonException
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 = -(YearBP -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBP
age)Float -> Float -> Float
forall a. Num a => a -> a -> a
+Float
1950
            ageSd2 :: YearBP
ageSd2 = YearBP
ageSdYearBP -> YearBP -> YearBP
forall a. Num a => a -> a -> a
*YearBP
ageSd
            ageSd2Float :: Float
ageSd2Float = YearBP -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBP
ageSd2
            musFloat :: Vector Float
musFloat = (YearBCAD -> Float) -> Vector YearBCAD -> Vector Float
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map YearBCAD -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBCAD
mus
            tau1sFloat :: Vector Float
tau1sFloat = (YearBP -> Float) -> Vector YearBP -> Vector Float
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map YearBP -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBP
tau1s
            dens :: Vector Float
dens = case CalibrationDistribution
distr of
                CalibrationDistribution
NormalDist -> 
                    (Float -> Float -> Float)
-> Vector Float -> Vector Float -> Vector Float
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 Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
mu) Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float -> Float
forall a. Floating a => a -> a
sqrt (Float
ageSd2Float Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
tau1 Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
tau1))) Vector Float
musFloat Vector Float
tau1sFloat
                StudentTDist Double
degreesOfFreedom -> 
                    (Float -> Float -> Float)
-> Vector Float -> Vector Float -> Vector Float
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 Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
mu) Float -> Float -> Float
forall a. Fractional a => a -> a -> a
/ Float -> Float
forall a. Floating a => a -> a
sqrt (Float
ageSd2Float Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
tau1 Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
tau1))) Vector Float
musFloat Vector Float
tau1sFloat
        in CalPDF -> Either CurrycarbonException CalPDF
forall a b. b -> Either a b
Right (CalPDF -> Either CurrycarbonException CalPDF)
-> CalPDF -> Either CurrycarbonException CalPDF
forall a b. (a -> b) -> a -> b
$ CalPDF -> CalPDF
trimLowDensityEdgesCalPDF (CalPDF -> CalPDF) -> CalPDF -> CalPDF
forall a b. (a -> b) -> a -> b
$ CalPDF -> CalPDF
normalizeCalPDF (CalPDF -> CalPDF) -> CalPDF -> CalPDF
forall a b. (a -> b) -> a -> b
$ String -> Vector YearBCAD -> Vector Float -> CalPDF
CalPDF String
name Vector YearBCAD
cals Vector Float
dens