{-# LANGUAGE Strict #-}

module Currycarbon.Calibration.Utils where

import           Currycarbon.Types

import           Data.Maybe            (fromMaybe)
import qualified Data.Vector.Unboxed   as VU
import           Numeric.SpecFunctions (logBeta)

-- | Rescale a CalPDF so that the sum of the densities is approx. 1.0
normalizeCalPDF :: CalPDF -> CalPDF
normalizeCalPDF :: CalPDF -> CalPDF
normalizeCalPDF (CalPDF String
name Vector YearBCAD
cals Vector Float
dens) =
    case forall a. (Unbox a, Num a) => Vector a -> a
VU.sum Vector Float
dens of
      Float
0.0 -> String -> Vector YearBCAD -> Vector Float -> CalPDF
CalPDF String
name Vector YearBCAD
cals Vector Float
dens -- product calibration can yield empty calPDFs
      Float
s   -> String -> Vector YearBCAD -> Vector Float -> CalPDF
CalPDF String
name Vector YearBCAD
cals forall a b. (a -> b) -> a -> b
$ forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map (forall a. Fractional a => a -> a -> a
/Float
s) Vector Float
dens

-- | get the density of a normal distribution at a point x
dnorm :: Float -> Float -> Float -> Float
dnorm :: Float -> Float -> Float -> Float
dnorm Float
mu Float
sigma Float
x =
    let a :: Float
a = forall a. Fractional a => a -> a
recip (forall a. Floating a => a -> a
sqrt (Float
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Float
sigma2))
        b :: Float
b = forall a. Floating a => a -> a
exp (-Float
c2 forall a. Fractional a => a -> a -> a
/ (Float
2 forall a. Num a => a -> a -> a
* Float
sigma2))
        c :: Float
c = Float
x forall a. Num a => a -> a -> a
- Float
mu
        c2 :: Float
c2 = Float
c forall a. Num a => a -> a -> a
* Float
c
        sigma2 :: Float
sigma2 = Float
sigma forall a. Num a => a -> a -> a
* Float
sigma
    in Float
aforall a. Num a => a -> a -> a
*Float
b
    -- alternative implemenation with the statistics package:
    -- import Statistics.Distribution (density)
    -- realToFrac $ density (normalDistr (realToFrac mu) (realToFrac sigma)) (realToFrac x)

-- | get the density of student's-t distribution at a point x
dt :: Double -> Float -> Float
dt :: Double -> Float -> Float
dt Double
dof Float
x =
    let xDouble :: Double
xDouble = forall a b. (Real a, Fractional b) => a -> b
realToFrac Float
x
        logDensityUnscaled :: Double
logDensityUnscaled = forall a. Floating a => a -> a
log (Double
dof forall a. Fractional a => a -> a -> a
/ (Double
dof forall a. Num a => a -> a -> a
+ Double
xDoubleforall a. Num a => a -> a -> a
*Double
xDouble)) forall a. Num a => a -> a -> a
* (Double
0.5 forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
+ Double
dof)) forall a. Num a => a -> a -> a
- Double -> Double -> Double
logBeta Double
0.5 (Double
0.5 forall a. Num a => a -> a -> a
* Double
dof)
    in forall a b. (Real a, Fractional b) => a -> b
realToFrac forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
exp Double
logDensityUnscaled forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Double
dof
    -- alternative implemenation with the statistics package:
    -- import Statistics.Distribution.StudentT (studentT)
    -- realToFrac $ density (studentT (realToFrac dof)) (realToFrac x) -- dof: number of degrees of freedom

isOutsideRangeOfCalCurve :: CalCurveBP -> UncalC14 -> Bool
isOutsideRangeOfCalCurve :: CalCurveBP -> UncalC14 -> Bool
isOutsideRangeOfCalCurve (CalCurveBP Vector Word
_ Vector Word
uncals Vector Word
_) (UncalC14 String
_ Word
age Word
_) =
    Word
age forall a. Ord a => a -> a -> Bool
< forall a. (Unbox a, Ord a) => Vector a -> a
VU.minimum Vector Word
uncals Bool -> Bool -> Bool
|| Word
age forall a. Ord a => a -> a -> Bool
> forall a. (Unbox a, Ord a) => Vector a -> a
VU.maximum Vector Word
uncals

-- | Take an uncalibrated date and a raw calibration curve and return
-- the relevant segment of the calibration curve
getRelevantCalCurveSegment :: UncalC14 -> CalCurveBP -> CalCurveBP
getRelevantCalCurveSegment :: UncalC14 -> CalCurveBP -> CalCurveBP
getRelevantCalCurveSegment (UncalC14 String
_ Word
mean Word
std) (CalCurveBP Vector Word
cals Vector Word
uncals Vector Word
sigmas) =
    let start :: Word
start = Word
meanforall a. Num a => a -> a -> a
+Word
6forall a. Num a => a -> a -> a
*Word
std
        stop :: Word
stop = Word
meanforall a. Num a => a -> a -> a
-Word
6forall a. Num a => a -> a -> a
*Word
std
        startIndex :: YearBCAD
startIndex = forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (forall a. Ord a => a -> a -> Bool
<= Word
start) Vector Word
uncals
        stopIndex :: YearBCAD
stopIndex = (forall a. Unbox a => Vector a -> YearBCAD
VU.length Vector Word
uncals forall a. Num a => a -> a -> a
- YearBCAD
1) forall a. Num a => a -> a -> a
- forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 (forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (forall a. Ord a => a -> a -> Bool
>= Word
stop) forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Vector a -> Vector a
VU.reverse Vector Word
uncals)
        toIndex :: YearBCAD
toIndex = YearBCAD
stopIndex forall a. Num a => a -> a -> a
- YearBCAD
startIndex
    in Vector Word -> Vector Word -> Vector Word -> CalCurveBP
CalCurveBP (forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
startIndex YearBCAD
toIndex Vector Word
cals) (forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
startIndex YearBCAD
toIndex Vector Word
uncals) (forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
startIndex YearBCAD
toIndex Vector Word
sigmas)

-- | Modify a calibration curve (segment) with multiple optional steps,
-- including interpolation and transforming dates to BC/AD format
prepareCalCurveSegment :: Bool -> CalCurveBP -> CalCurveBCAD
prepareCalCurveSegment :: Bool -> CalCurveBP -> CalCurveBCAD
prepareCalCurveSegment Bool
interpolate CalCurveBP
calCurve =
    CalCurveBP -> CalCurveBCAD
makeBCADCalCurve forall a b. (a -> b) -> a -> b
$ if Bool
interpolate then CalCurveBP -> CalCurveBP
interpolateCalCurve CalCurveBP
calCurve else CalCurveBP
calCurve

makeBCADCalCurve :: CalCurveBP -> CalCurveBCAD
makeBCADCalCurve :: CalCurveBP -> CalCurveBCAD
makeBCADCalCurve (CalCurveBP Vector Word
cals Vector Word
uncals Vector Word
sigmas) = Vector YearBCAD -> Vector YearBCAD -> Vector Word -> CalCurveBCAD
CalCurveBCAD (Vector Word -> Vector YearBCAD
vectorBPToBCAD Vector Word
cals) (Vector Word -> Vector YearBCAD
vectorBPToBCAD Vector Word
uncals) Vector Word
sigmas

vectorBPToBCAD :: VU.Vector YearBP -> VU.Vector YearBCAD
vectorBPToBCAD :: Vector Word -> Vector YearBCAD
vectorBPToBCAD = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map Word -> YearBCAD
bp2BCAD

bp2BCAD :: YearBP -> YearBCAD
bp2BCAD :: Word -> YearBCAD
bp2BCAD Word
x = -(forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
x) forall a. Num a => a -> a -> a
+ YearBCAD
1950

interpolateCalCurve :: CalCurveBP -> CalCurveBP
interpolateCalCurve :: CalCurveBP -> CalCurveBP
interpolateCalCurve (CalCurveBP Vector Word
cals Vector Word
uncals Vector Word
sigmas) =
    let obs :: Vector (Word, Word, Word)
obs = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
Vector a -> Vector b -> Vector c -> Vector (a, b, c)
VU.zip3 Vector Word
cals Vector Word
uncals Vector Word
sigmas
        timeWindows :: Vector ((Word, Word, Word), (Word, Word, Word))
timeWindows = Vector (Word, Word, Word)
-> Vector ((Word, Word, Word), (Word, Word, Word))
getTimeWindows Vector (Word, Word, Word)
obs
        obsFilled :: Vector (Word, Word, Word)
obsFilled = forall a b.
(Unbox a, Unbox b) =>
(a -> Vector b) -> Vector a -> Vector b
VU.concatMap ((Word, Word, Word), (Word, Word, Word))
-> Vector (Word, Word, Word)
fillTimeWindows Vector ((Word, Word, Word), (Word, Word, Word))
timeWindows
    in forall a b c d. (a -> b -> c -> d) -> (a, b, c) -> d
uncurry3 Vector Word -> Vector Word -> Vector Word -> CalCurveBP
CalCurveBP forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
Vector (a, b, c) -> (Vector a, Vector b, Vector c)
VU.unzip3 Vector (Word, Word, Word)
obsFilled
    where
        getTimeWindows :: VU.Vector (YearBP,YearBP,YearRange) -> VU.Vector ((YearBP,YearBP,YearRange),(YearBP,YearBP,YearRange))
        getTimeWindows :: Vector (Word, Word, Word)
-> Vector ((Word, Word, Word), (Word, Word, Word))
getTimeWindows Vector (Word, Word, Word)
xs = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VU.zipWith (,) (forall a. Unbox a => Vector a -> Vector a
VU.init Vector (Word, Word, Word)
xs) (forall a. Unbox a => Vector a -> Vector a
VU.tail Vector (Word, Word, Word)
xs)
        fillTimeWindows :: ((YearBP,YearBP,YearRange),(YearBP,YearBP,YearRange)) -> VU.Vector (YearBP,YearBP,YearRange)
        fillTimeWindows :: ((Word, Word, Word), (Word, Word, Word))
-> Vector (Word, Word, Word)
fillTimeWindows ((Word
calbp1,Word
bp1,Word
sigma1),(Word
calbp2,Word
bp2,Word
sigma2)) =
            if Word
calbp1 forall a. Eq a => a -> a -> Bool
== Word
calbp2 Bool -> Bool -> Bool
|| Word
calbp1forall a. Num a => a -> a -> a
+Word
1 forall a. Eq a => a -> a -> Bool
== Word
calbp2 Bool -> Bool -> Bool
|| Word
calbp1forall a. Num a => a -> a -> a
-Word
1 forall a. Eq a => a -> a -> Bool
== Word
calbp2
            then forall a. Unbox a => a -> Vector a
VU.singleton (Word
calbp1,Word
bp1,Word
sigma1)
            else
                let newCals :: Vector Word
newCals = forall a. Unbox a => [a] -> Vector a
VU.fromList [Word
calbp1,Word
calbp1forall a. Num a => a -> a -> a
-Word
1..Word
calbp2forall a. Num a => a -> a -> a
+Word
1] -- range definition like this to trigger counting down
                    newBPs :: Vector Word
newBPs = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map (forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Word, Word) -> (Word, Word) -> Word -> (Word, Word)
getInBetweenPointsInt (Word
calbp1,Word
bp1) (Word
calbp2,Word
bp2)) Vector Word
newCals
                    newSigmas :: Vector Word
newSigmas = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map (forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Word, Word) -> (Word, Word) -> Word -> (Word, Word)
getInBetweenPointsInt (Word
calbp1,Word
sigma1) (Word
calbp2,Word
sigma2)) Vector Word
newCals
                in forall a b c.
(Unbox a, Unbox b, Unbox c) =>
Vector a -> Vector b -> Vector c -> Vector (a, b, c)
VU.zip3 Vector Word
newCals Vector Word
newBPs Vector Word
newSigmas
        getInBetweenPointsInt :: (Word, Word) -> (Word, Word) -> Word -> (Word, Word)
        getInBetweenPointsInt :: (Word, Word) -> (Word, Word) -> Word -> (Word, Word)
getInBetweenPointsInt (Word
x1,Word
y1) (Word
x2,Word
y2) Word
xPred =
            let (Float
_,Float
yPred) = (Float, Float) -> (Float, Float) -> Float -> (Float, Float)
getInBetweenPoints (forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
x1,forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
y1) (forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
x2,forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
y2) forall a b. (a -> b) -> a -> b
$ forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
xPred
            in (Word
xPred, forall a b. (RealFrac a, Integral b) => a -> b
round Float
yPred)
        getInBetweenPoints :: (Float, Float) -> (Float, Float) -> Float -> (Float, Float)
        getInBetweenPoints :: (Float, Float) -> (Float, Float) -> Float -> (Float, Float)
getInBetweenPoints (Float
x1,Float
y1) (Float
x2,Float
y2) Float
xPred =
            let yDiff :: Float
yDiff = Float
y2 forall a. Num a => a -> a -> a
- Float
y1
                xDiff :: Float
xDiff = forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ Float
x1 forall a. Num a => a -> a -> a
- Float
x2
                yDiffPerxDiff :: Float
yDiffPerxDiff = Float
yDiffforall a. Fractional a => a -> a -> a
/Float
xDiff
                xPredRel :: Float
xPredRel = Float
x1 forall a. Num a => a -> a -> a
- Float
xPred
            in (Float
xPred, Float
y1 forall a. Num a => a -> a -> a
+ Float
xPredRel forall a. Num a => a -> a -> a
* Float
yDiffPerxDiff)
        uncurry3 :: (a -> b -> c -> d) -> ((a, b, c) -> d)
        uncurry3 :: forall a b c d. (a -> b -> c -> d) -> (a, b, c) -> d
uncurry3 a -> b -> c -> d
f ~(a
a,b
b,c
c) = a -> b -> c -> d
f a
a b
b c
c

trimLowDensityEdgesCalPDF :: CalPDF -> CalPDF
trimLowDensityEdgesCalPDF :: CalPDF -> CalPDF
trimLowDensityEdgesCalPDF (CalPDF String
name Vector YearBCAD
cals Vector Float
dens) =
    let firstAboveThreshold :: YearBCAD
firstAboveThreshold = forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 (forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (forall a. Ord a => a -> a -> Bool
> Float
0.00001) Vector Float
dens)
        lastAboveThreshold :: YearBCAD
lastAboveThreshold = forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 (forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (forall a. Ord a => a -> a -> Bool
> Float
0.00001) forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Vector a -> Vector a
VU.reverse Vector Float
dens)
        untilLastAboveThreshold :: YearBCAD
untilLastAboveThreshold = forall a. Unbox a => Vector a -> YearBCAD
VU.length Vector Float
dens forall a. Num a => a -> a -> a
- YearBCAD
firstAboveThreshold forall a. Num a => a -> a -> a
- YearBCAD
lastAboveThreshold
        calsSlice :: Vector YearBCAD
calsSlice = forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
firstAboveThreshold YearBCAD
untilLastAboveThreshold Vector YearBCAD
cals
        densSlice :: Vector Float
densSlice = forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
firstAboveThreshold YearBCAD
untilLastAboveThreshold Vector Float
dens
    in String -> Vector YearBCAD -> Vector Float -> CalPDF
CalPDF String
name Vector YearBCAD
calsSlice Vector Float
densSlice