{-# 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)

-- 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 :: forall a c b d. (a -> c) -> (b -> d) -> Either a b -> Either c d
mapEither a -> c
f b -> d
_ (Left a
x)  = c -> Either c d
forall a b. a -> Either a b
Left (a -> c
f a
x)
mapEither a -> c
_ b -> d
f (Right b
x) = d -> Either c d
forall a b. b -> Either a b
Right (b -> d
f b
x)

-- | 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 Double
dens) =
    case Vector Double -> Double
forall a. (Unbox a, Num a) => Vector a -> a
VU.sum Vector Double
dens of
      Double
0.0 -> String -> Vector YearBCAD -> Vector Double -> CalPDF
CalPDF String
name Vector YearBCAD
cals Vector Double
dens -- product calibration can yield empty calPDFs
      Double
s   -> String -> Vector YearBCAD -> Vector Double -> CalPDF
CalPDF String
name Vector YearBCAD
cals (Vector Double -> CalPDF) -> Vector Double -> CalPDF
forall a b. (a -> b) -> a -> b
$ (Double -> Double) -> Vector Double -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
s) Vector Double
dens

-- | get the density of a normal distribution at a point x
dnorm :: Double -> Double -> Double -> Double
dnorm :: Double -> Double -> Double -> Double
dnorm Double
mu Double
sigma Double
x =
    let a :: Double
a = Double -> Double
forall a. Fractional a => a -> a
recip (Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sigma2))
        b :: Double
b = Double -> Double
forall a. Floating a => a -> a
exp (-Double
c2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sigma2))
        c :: Double
c = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mu
        c2 :: Double
c2 = Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c
        sigma2 :: Double
sigma2 = Double
sigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sigma
    in Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
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 -> Double -> Double
dt :: Double -> Double -> Double
dt Double
dof Double
x =
    let xDouble :: Double
xDouble = Double -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
x
        logDensityUnscaled :: Double
logDensityUnscaled = Double -> Double
forall a. Floating a => a -> a
log (Double
dof Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
dof Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
xDoubleDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dof)) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double
logBeta Double
0.5 (Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dof)
    in Double -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
exp Double
logDensityUnscaled Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
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 Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
< Vector Word -> Word
forall a. (Unbox a, Ord a) => Vector a -> a
VU.minimum Vector Word
uncals Bool -> Bool -> Bool
|| Word
age Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
> Vector Word -> Word
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
meanWord -> Word -> Word
forall a. Num a => a -> a -> a
+Word
6Word -> Word -> Word
forall a. Num a => a -> a -> a
*Word
std
        stop :: Word
stop = Word
meanWord -> Word -> Word
forall a. Num a => a -> a -> a
-Word
6Word -> Word -> Word
forall a. Num a => a -> a -> a
*Word
std
        startIndex :: YearBCAD
startIndex = YearBCAD -> Maybe YearBCAD -> YearBCAD
forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 (Maybe YearBCAD -> YearBCAD) -> Maybe YearBCAD -> YearBCAD
forall a b. (a -> b) -> a -> b
$ (Word -> Bool) -> Vector Word -> Maybe YearBCAD
forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
<= Word
start) Vector Word
uncals
        stopIndex :: YearBCAD
stopIndex = (Vector Word -> YearBCAD
forall a. Unbox a => Vector a -> YearBCAD
VU.length Vector Word
uncals YearBCAD -> YearBCAD -> YearBCAD
forall a. Num a => a -> a -> a
- YearBCAD
1) YearBCAD -> YearBCAD -> YearBCAD
forall a. Num a => a -> a -> a
- YearBCAD -> Maybe YearBCAD -> YearBCAD
forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 ((Word -> Bool) -> Vector Word -> Maybe YearBCAD
forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
>= Word
stop) (Vector Word -> Maybe YearBCAD) -> Vector Word -> Maybe YearBCAD
forall a b. (a -> b) -> a -> b
$ Vector Word -> Vector Word
forall a. Unbox a => Vector a -> Vector a
VU.reverse Vector Word
uncals)
        toIndex :: YearBCAD
toIndex = YearBCAD
stopIndex YearBCAD -> YearBCAD -> YearBCAD
forall a. Num a => a -> a -> a
- YearBCAD
startIndex
    in Vector Word -> Vector Word -> Vector Word -> CalCurveBP
CalCurveBP (YearBCAD -> YearBCAD -> Vector Word -> Vector Word
forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
startIndex YearBCAD
toIndex Vector Word
cals) (YearBCAD -> YearBCAD -> Vector Word -> Vector Word
forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
startIndex YearBCAD
toIndex Vector Word
uncals) (YearBCAD -> YearBCAD -> Vector Word -> Vector Word
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 (CalCurveBP -> CalCurveBCAD) -> CalCurveBP -> CalCurveBCAD
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

punchOutCalCurveBCAD :: Int -> Int -> CalCurveBCAD -> CalCurveBCAD
punchOutCalCurveBCAD :: YearBCAD -> YearBCAD -> CalCurveBCAD -> CalCurveBCAD
punchOutCalCurveBCAD YearBCAD
start YearBCAD
stop (CalCurveBCAD Vector YearBCAD
cals Vector YearBCAD
uncals Vector Word
sigmas) =
    let startIndex :: YearBCAD
startIndex = YearBCAD -> Maybe YearBCAD -> YearBCAD
forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 (Maybe YearBCAD -> YearBCAD) -> Maybe YearBCAD -> YearBCAD
forall a b. (a -> b) -> a -> b
$ (YearBCAD -> Bool) -> Vector YearBCAD -> Maybe YearBCAD
forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (YearBCAD -> YearBCAD -> Bool
forall a. Ord a => a -> a -> Bool
>= YearBCAD
start) Vector YearBCAD
cals
        stopIndex :: YearBCAD
stopIndex = Vector YearBCAD -> YearBCAD
forall a. Unbox a => Vector a -> YearBCAD
VU.length Vector YearBCAD
cals YearBCAD -> YearBCAD -> YearBCAD
forall a. Num a => a -> a -> a
- YearBCAD -> Maybe YearBCAD -> YearBCAD
forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 ((YearBCAD -> Bool) -> Vector YearBCAD -> Maybe YearBCAD
forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (YearBCAD -> YearBCAD -> Bool
forall a. Ord a => a -> a -> Bool
<= YearBCAD
stop) (Vector YearBCAD -> Maybe YearBCAD)
-> Vector YearBCAD -> Maybe YearBCAD
forall a b. (a -> b) -> a -> b
$ Vector YearBCAD -> Vector YearBCAD
forall a. Unbox a => Vector a -> Vector a
VU.reverse Vector YearBCAD
cals)
        toIndex :: YearBCAD
toIndex = YearBCAD
stopIndex YearBCAD -> YearBCAD -> YearBCAD
forall a. Num a => a -> a -> a
- YearBCAD
startIndex
    --in error $ show $ (start, stop, VU.slice startIndex toIndex cals)
    in Vector YearBCAD -> Vector YearBCAD -> Vector Word -> CalCurveBCAD
CalCurveBCAD
       (YearBCAD -> YearBCAD -> Vector YearBCAD -> Vector YearBCAD
forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
startIndex YearBCAD
toIndex Vector YearBCAD
cals)
       (YearBCAD -> YearBCAD -> Vector YearBCAD -> Vector YearBCAD
forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
startIndex YearBCAD
toIndex Vector YearBCAD
uncals)
       (YearBCAD -> YearBCAD -> Vector Word -> Vector Word
forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
startIndex YearBCAD
toIndex Vector Word
sigmas)

vectorBPToBCAD :: VU.Vector YearBP -> VU.Vector YearBCAD
vectorBPToBCAD :: Vector Word -> Vector YearBCAD
vectorBPToBCAD = (Word -> YearBCAD) -> Vector Word -> Vector YearBCAD
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 = -(Word -> YearBCAD
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
x) YearBCAD -> YearBCAD -> YearBCAD
forall a. Num a => a -> a -> a
+ YearBCAD
1950

bcad2BP :: YearBCAD -> YearBP
bcad2BP :: YearBCAD -> Word
bcad2BP YearBCAD
y = Word
1950 Word -> Word -> Word
forall a. Num a => a -> a -> a
- YearBCAD -> Word
forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBCAD
y

interpolateCalCurve :: CalCurveBP -> CalCurveBP
interpolateCalCurve :: CalCurveBP -> CalCurveBP
interpolateCalCurve (CalCurveBP Vector Word
cals Vector Word
uncals Vector Word
sigmas) =
    let obs :: Vector (Word, Word, Word)
obs = Vector Word
-> Vector Word -> Vector Word -> Vector (Word, Word, Word)
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 = (((Word, Word, Word), (Word, Word, Word))
 -> Vector (Word, Word, Word))
-> Vector ((Word, Word, Word), (Word, Word, Word))
-> Vector (Word, Word, Word)
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 (Vector Word -> Vector Word -> Vector Word -> CalCurveBP)
-> (Vector Word, Vector Word, Vector Word) -> CalCurveBP
forall a b c d. (a -> b -> c -> d) -> (a, b, c) -> d
uncurry3 Vector Word -> Vector Word -> Vector Word -> CalCurveBP
CalCurveBP ((Vector Word, Vector Word, Vector Word) -> CalCurveBP)
-> (Vector Word, Vector Word, Vector Word) -> CalCurveBP
forall a b. (a -> b) -> a -> b
$ Vector (Word, Word, Word)
-> (Vector Word, Vector Word, Vector Word)
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 = ((Word, Word, Word)
 -> (Word, Word, Word) -> ((Word, Word, Word), (Word, Word, Word)))
-> Vector (Word, Word, Word)
-> Vector (Word, Word, Word)
-> Vector ((Word, Word, Word), (Word, Word, Word))
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VU.zipWith (,) (Vector (Word, Word, Word) -> Vector (Word, Word, Word)
forall a. Unbox a => Vector a -> Vector a
VU.init Vector (Word, Word, Word)
xs) (Vector (Word, Word, Word) -> Vector (Word, Word, Word)
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 Word -> Word -> Bool
forall a. Eq a => a -> a -> Bool
== Word
calbp2 Bool -> Bool -> Bool
|| Word
calbp1Word -> Word -> Word
forall a. Num a => a -> a -> a
+Word
1 Word -> Word -> Bool
forall a. Eq a => a -> a -> Bool
== Word
calbp2 Bool -> Bool -> Bool
|| Word
calbp1Word -> Word -> Word
forall a. Num a => a -> a -> a
-Word
1 Word -> Word -> Bool
forall a. Eq a => a -> a -> Bool
== Word
calbp2
            then (Word, Word, Word) -> Vector (Word, Word, Word)
forall a. Unbox a => a -> Vector a
VU.singleton (Word
calbp1,Word
bp1,Word
sigma1)
            else
                let newCals :: Vector Word
newCals = [Word] -> Vector Word
forall a. Unbox a => [a] -> Vector a
VU.fromList [Word
calbp1,Word
calbp1Word -> Word -> Word
forall a. Num a => a -> a -> a
-Word
1..Word
calbp2Word -> Word -> Word
forall a. Num a => a -> a -> a
+Word
1] -- range definition like this to trigger counting down
                    newBPs :: Vector Word
newBPs = (Word -> Word) -> Vector Word -> Vector Word
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map ((Word, Word) -> Word
forall a b. (a, b) -> b
snd ((Word, Word) -> Word) -> (Word -> (Word, Word)) -> Word -> Word
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 = (Word -> Word) -> Vector Word -> Vector Word
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map ((Word, Word) -> Word
forall a b. (a, b) -> b
snd ((Word, Word) -> Word) -> (Word -> (Word, Word)) -> Word -> Word
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 Vector Word
-> Vector Word -> Vector Word -> Vector (Word, Word, Word)
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 (Double
_,Double
yPred) = (Double, Double) -> (Double, Double) -> Double -> (Double, Double)
getInBetweenPoints (Word -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
x1,Word -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
y1) (Word -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
x2,Word -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
y2) (Double -> (Double, Double)) -> Double -> (Double, Double)
forall a b. (a -> b) -> a -> b
$ Word -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
xPred
            in (Word
xPred, Double -> Word
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round Double
yPred)
        getInBetweenPoints :: (Double, Double) -> (Double, Double) -> Double -> (Double, Double)
        getInBetweenPoints :: (Double, Double) -> (Double, Double) -> Double -> (Double, Double)
getInBetweenPoints (Double
x1,Double
y1) (Double
x2,Double
y2) Double
xPred =
            let yDiff :: Double
yDiff = Double
y2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y1
                xDiff :: Double
xDiff = Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
x1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x2
                yDiffPerxDiff :: Double
yDiffPerxDiff = Double
yDiffDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
xDiff
                xPredRel :: Double
xPredRel = Double
x1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xPred
            in (Double
xPred, Double
y1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
xPredRel Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
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

-- | Subset the calibration curve to the non-zero density range. The threshold is set to 0.00001.
trimLowDensityEdgesCalPDF :: CalPDF -> CalPDF
trimLowDensityEdgesCalPDF :: CalPDF -> CalPDF
trimLowDensityEdgesCalPDF (CalPDF String
name Vector YearBCAD
cals Vector Double
dens) =
    let firstAboveThreshold :: YearBCAD
firstAboveThreshold = YearBCAD -> Maybe YearBCAD -> YearBCAD
forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 ((Double -> Bool) -> Vector Double -> Maybe YearBCAD
forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.00001) Vector Double
dens)
        lastAboveThreshold :: YearBCAD
lastAboveThreshold = YearBCAD -> Maybe YearBCAD -> YearBCAD
forall a. a -> Maybe a -> a
fromMaybe YearBCAD
0 ((Double -> Bool) -> Vector Double -> Maybe YearBCAD
forall a. Unbox a => (a -> Bool) -> Vector a -> Maybe YearBCAD
VU.findIndex (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.00001) (Vector Double -> Maybe YearBCAD)
-> Vector Double -> Maybe YearBCAD
forall a b. (a -> b) -> a -> b
$ Vector Double -> Vector Double
forall a. Unbox a => Vector a -> Vector a
VU.reverse Vector Double
dens)
        untilLastAboveThreshold :: YearBCAD
untilLastAboveThreshold = Vector Double -> YearBCAD
forall a. Unbox a => Vector a -> YearBCAD
VU.length Vector Double
dens YearBCAD -> YearBCAD -> YearBCAD
forall a. Num a => a -> a -> a
- YearBCAD
firstAboveThreshold YearBCAD -> YearBCAD -> YearBCAD
forall a. Num a => a -> a -> a
- YearBCAD
lastAboveThreshold
        calsSlice :: Vector YearBCAD
calsSlice = YearBCAD -> YearBCAD -> Vector YearBCAD -> Vector YearBCAD
forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
firstAboveThreshold YearBCAD
untilLastAboveThreshold Vector YearBCAD
cals
        densSlice :: Vector Double
densSlice = YearBCAD -> YearBCAD -> Vector Double -> Vector Double
forall a. Unbox a => YearBCAD -> YearBCAD -> Vector a -> Vector a
VU.slice YearBCAD
firstAboveThreshold YearBCAD
untilLastAboveThreshold Vector Double
dens
    in String -> Vector YearBCAD -> Vector Double -> CalPDF
CalPDF String
name Vector YearBCAD
calsSlice Vector Double
densSlice