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