{-# LANGUAGE Strict #-}
module Currycarbon.Calibration.Calibration
(
getRelevantCalCurveSegment
, prepareCalCurveSegment
, makeCalCurveMatrix
, uncalToPDF
, calibrateDate
, calibrateDates
, refineCalDates
, refineCalDate
, CalibrateDatesConf (..)
, defaultCalConf
) where
import Currycarbon.Calibration.Utils
import Currycarbon.Calibration.Bchron
import Currycarbon.Calibration.MatrixMult
import Currycarbon.Types
import Currycarbon.Utils
import Data.List (sort, sortBy, groupBy, elemIndex)
import Data.Maybe (fromJust)
import qualified Data.Vector.Unboxed as VU
data CalibrateDatesConf = CalibrateDatesConf {
CalibrateDatesConf -> CalibrationMethod
_calConfMethod :: CalibrationMethod
, CalibrateDatesConf -> Bool
_calConfAllowOutside :: Bool
, CalibrateDatesConf -> Bool
_calConfInterpolateCalCurve :: Bool
} deriving (Int -> CalibrateDatesConf -> ShowS
[CalibrateDatesConf] -> ShowS
CalibrateDatesConf -> String
(Int -> CalibrateDatesConf -> ShowS)
-> (CalibrateDatesConf -> String)
-> ([CalibrateDatesConf] -> ShowS)
-> Show CalibrateDatesConf
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [CalibrateDatesConf] -> ShowS
$cshowList :: [CalibrateDatesConf] -> ShowS
show :: CalibrateDatesConf -> String
$cshow :: CalibrateDatesConf -> String
showsPrec :: Int -> CalibrateDatesConf -> ShowS
$cshowsPrec :: Int -> CalibrateDatesConf -> ShowS
Show, CalibrateDatesConf -> CalibrateDatesConf -> Bool
(CalibrateDatesConf -> CalibrateDatesConf -> Bool)
-> (CalibrateDatesConf -> CalibrateDatesConf -> Bool)
-> Eq CalibrateDatesConf
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: CalibrateDatesConf -> CalibrateDatesConf -> Bool
$c/= :: CalibrateDatesConf -> CalibrateDatesConf -> Bool
== :: CalibrateDatesConf -> CalibrateDatesConf -> Bool
$c== :: CalibrateDatesConf -> CalibrateDatesConf -> Bool
Eq)
defaultCalConf :: CalibrateDatesConf
defaultCalConf :: CalibrateDatesConf
defaultCalConf = CalibrateDatesConf :: CalibrationMethod -> Bool -> Bool -> CalibrateDatesConf
CalibrateDatesConf {
_calConfMethod :: CalibrationMethod
_calConfMethod = Bchron :: CalibrationDistribution -> CalibrationMethod
Bchron { distribution :: CalibrationDistribution
distribution = Double -> CalibrationDistribution
StudentTDist Double
100 }
, _calConfAllowOutside :: Bool
_calConfAllowOutside = Bool
False
, _calConfInterpolateCalCurve :: Bool
_calConfInterpolateCalCurve = Bool
True
}
calibrateDates :: CalibrateDatesConf
-> CalCurveBP
-> [UncalC14]
-> [Either CurrycarbonException CalPDF]
calibrateDates :: CalibrateDatesConf
-> CalCurveBP -> [UncalC14] -> [Either CurrycarbonException CalPDF]
calibrateDates CalibrateDatesConf
_ CalCurveBP
_ [] = []
calibrateDates (CalibrateDatesConf CalibrationMethod
MatrixMultiplication Bool
allowOutside Bool
interpolate) CalCurveBP
calCurve [UncalC14]
uncalDates =
(UncalC14 -> Either CurrycarbonException CalPDF)
-> [UncalC14] -> [Either CurrycarbonException CalPDF]
forall a b. (a -> b) -> [a] -> [b]
map (Bool
-> Bool
-> CalCurveBP
-> UncalC14
-> Either CurrycarbonException CalPDF
calibrateDateMatrixMult Bool
allowOutside Bool
interpolate CalCurveBP
calCurve) [UncalC14]
uncalDates
calibrateDates (CalibrateDatesConf Bchron{distribution :: CalibrationMethod -> CalibrationDistribution
distribution=CalibrationDistribution
distr} Bool
allowOutside Bool
interpolate) CalCurveBP
calCurve [UncalC14]
uncalDates =
(UncalC14 -> Either CurrycarbonException CalPDF)
-> [UncalC14] -> [Either CurrycarbonException CalPDF]
forall a b. (a -> b) -> [a] -> [b]
map (CalibrationDistribution
-> Bool
-> Bool
-> CalCurveBP
-> UncalC14
-> Either CurrycarbonException CalPDF
calibrateDateBchron CalibrationDistribution
distr Bool
allowOutside Bool
interpolate CalCurveBP
calCurve) [UncalC14]
uncalDates
calibrateDate :: CalibrateDatesConf
-> CalCurveBP
-> UncalC14
-> Either CurrycarbonException CalPDF
calibrateDate :: CalibrateDatesConf
-> CalCurveBP -> UncalC14 -> Either CurrycarbonException CalPDF
calibrateDate (CalibrateDatesConf CalibrationMethod
MatrixMultiplication Bool
allowOutside Bool
interpolate) CalCurveBP
calCurve UncalC14
uncalDate =
Bool
-> Bool
-> CalCurveBP
-> UncalC14
-> Either CurrycarbonException CalPDF
calibrateDateMatrixMult Bool
allowOutside Bool
interpolate CalCurveBP
calCurve UncalC14
uncalDate
calibrateDate (CalibrateDatesConf Bchron{distribution :: CalibrationMethod -> CalibrationDistribution
distribution=CalibrationDistribution
distr} Bool
allowOutside Bool
interpolate) CalCurveBP
calCurve UncalC14
uncalDate =
CalibrationDistribution
-> Bool
-> Bool
-> CalCurveBP
-> UncalC14
-> Either CurrycarbonException CalPDF
calibrateDateBchron CalibrationDistribution
distr Bool
allowOutside Bool
interpolate CalCurveBP
calCurve UncalC14
uncalDate
refineCalDates :: [CalPDF] -> [Maybe CalC14]
refineCalDates :: [CalPDF] -> [Maybe CalC14]
refineCalDates = (CalPDF -> Maybe CalC14) -> [CalPDF] -> [Maybe CalC14]
forall a b. (a -> b) -> [a] -> [b]
map CalPDF -> Maybe CalC14
refineCalDate
refineCalDate :: CalPDF -> Maybe CalC14
refineCalDate :: CalPDF -> Maybe CalC14
refineCalDate (CalPDF String
name Vector Int
cals Vector Float
dens) =
if Vector Float -> Float
forall a. (Unbox a, Num a) => Vector a -> a
VU.sum Vector Float
dens Float -> Float -> Bool
forall a. Eq a => a -> a -> Bool
== Float
0 Bool -> Bool -> Bool
|| Vector Float -> Int
forall a. Unbox a => Vector a -> Int
VU.length ((Float -> Bool) -> Vector Float -> Vector Float
forall a. Unbox a => (a -> Bool) -> Vector a -> Vector a
VU.filter (Float -> Float -> Bool
forall a. Ord a => a -> a -> Bool
>= Float
1.0) Vector Float
dens) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
then Maybe CalC14
forall a. Maybe a
Nothing
else CalC14 -> Maybe CalC14
forall a. a -> Maybe a
Just (CalC14 -> Maybe CalC14) -> CalC14 -> Maybe CalC14
forall a b. (a -> b) -> a -> b
$ CalC14 :: String -> CalRangeSummary -> [HDR] -> [HDR] -> CalC14
CalC14 {
_calC14id :: String
_calC14id = String
name
, _calC14RangeSummary :: CalRangeSummary
_calC14RangeSummary = CalRangeSummary :: Int -> Int -> Int -> Int -> Int -> CalRangeSummary
CalRangeSummary {
_calRangeStartTwoSigma :: Int
_calRangeStartTwoSigma = HDR -> Int
_hdrstart (HDR -> Int) -> HDR -> Int
forall a b. (a -> b) -> a -> b
$ [HDR] -> HDR
forall a. [a] -> a
head [HDR]
hdrs95
, _calRangeStartOneSigma :: Int
_calRangeStartOneSigma = HDR -> Int
_hdrstart (HDR -> Int) -> HDR -> Int
forall a b. (a -> b) -> a -> b
$ [HDR] -> HDR
forall a. [a] -> a
head [HDR]
hdrs68
, _calRangeMedian :: Int
_calRangeMedian = Maybe Int -> Int
forall a. HasCallStack => Maybe a -> a
fromJust (Maybe Int -> Int) -> Maybe Int -> Int
forall a b. (a -> b) -> a -> b
$ Vector Int
cals Vector Int -> Maybe Int -> Maybe Int
forall a. Unbox a => Vector a -> Maybe Int -> Maybe a
`indexVU` Float -> [Float] -> Maybe Int
forall a. Eq a => a -> [a] -> Maybe Int
elemIndex ([Float] -> Float
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Float]
distanceTo05) [Float]
distanceTo05
, _calRangeStopOneSigma :: Int
_calRangeStopOneSigma = HDR -> Int
_hdrstop (HDR -> Int) -> HDR -> Int
forall a b. (a -> b) -> a -> b
$ [HDR] -> HDR
forall a. [a] -> a
last [HDR]
hdrs68
, _calRangeStopTwoSigma :: Int
_calRangeStopTwoSigma = HDR -> Int
_hdrstop (HDR -> Int) -> HDR -> Int
forall a b. (a -> b) -> a -> b
$ [HDR] -> HDR
forall a. [a] -> a
last [HDR]
hdrs95
}
, _calC14HDROneSigma :: [HDR]
_calC14HDROneSigma = [HDR]
hdrs68
, _calC14HDRTwoSigma :: [HDR]
_calC14HDRTwoSigma = [HDR]
hdrs95
}
where
cumsumDensities :: [Float]
cumsumDensities = [(Int, Float)] -> [Float]
cumsumDens (Vector (Int, Float) -> [(Int, Float)]
forall a. Unbox a => Vector a -> [a]
VU.toList (Vector (Int, Float) -> [(Int, Float)])
-> Vector (Int, Float) -> [(Int, Float)]
forall a b. (a -> b) -> a -> b
$ Vector Int -> Vector Float -> Vector (Int, Float)
forall a b.
(Unbox a, Unbox b) =>
Vector a -> Vector b -> Vector (a, b)
VU.zip Vector Int
cals Vector Float
dens)
distanceTo05 :: [Float]
distanceTo05 = (Float -> Float) -> [Float] -> [Float]
forall a b. (a -> b) -> [a] -> [b]
map (\Float
x -> Float -> Float
forall a. Num a => a -> a
abs (Float -> Float) -> Float -> Float
forall a b. (a -> b) -> a -> b
$ (Float
x Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
0.5)) [Float]
cumsumDensities
sortedDensities :: [(Int, Float)]
sortedDensities = ((Int, Float) -> (Int, Float) -> Ordering)
-> [(Int, Float)] -> [(Int, Float)]
forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (((Int, Float) -> (Int, Float) -> Ordering)
-> (Int, Float) -> (Int, Float) -> Ordering
forall a b c. (a -> b -> c) -> b -> a -> c
flip (\ (Int
_, Float
dens1) (Int
_, Float
dens2) -> Float -> Float -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Float
dens1 Float
dens2)) (Vector (Int, Float) -> [(Int, Float)]
forall a. Unbox a => Vector a -> [a]
VU.toList (Vector (Int, Float) -> [(Int, Float)])
-> Vector (Int, Float) -> [(Int, Float)]
forall a b. (a -> b) -> a -> b
$ Vector Int -> Vector Float -> Vector (Int, Float)
forall a b.
(Unbox a, Unbox b) =>
Vector a -> Vector b -> Vector (a, b)
VU.zip Vector Int
cals Vector Float
dens)
cumsumSortedDensities :: [Float]
cumsumSortedDensities = [(Int, Float)] -> [Float]
cumsumDens [(Int, Float)]
sortedDensities
isIn68 :: [Bool]
isIn68 = (Float -> Bool) -> [Float] -> [Bool]
forall a b. (a -> b) -> [a] -> [b]
map (Float -> Float -> Bool
forall a. Ord a => a -> a -> Bool
< Float
0.683) [Float]
cumsumSortedDensities
isIn95 :: [Bool]
isIn95 = (Float -> Bool) -> [Float] -> [Bool]
forall a b. (a -> b) -> [a] -> [b]
map (Float -> Float -> Bool
forall a. Ord a => a -> a -> Bool
< Float
0.954) [Float]
cumsumSortedDensities
contextualizedDensities :: [(Int, Float, Bool, Bool)]
contextualizedDensities = [(Int, Float, Bool, Bool)] -> [(Int, Float, Bool, Bool)]
forall a. Ord a => [a] -> [a]
sort ([(Int, Float, Bool, Bool)] -> [(Int, Float, Bool, Bool)])
-> [(Int, Float, Bool, Bool)] -> [(Int, Float, Bool, Bool)]
forall a b. (a -> b) -> a -> b
$ ((Int, Float) -> Bool -> Bool -> (Int, Float, Bool, Bool))
-> [(Int, Float)] -> [Bool] -> [Bool] -> [(Int, Float, Bool, Bool)]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\(Int
y,Float
d) Bool
in68 Bool
in95 -> (Int
y,Float
d,Bool
in68,Bool
in95)) [(Int, Float)]
sortedDensities [Bool]
isIn68 [Bool]
isIn95
hdrs68 :: [HDR]
hdrs68 = [(Int, Float, Bool, Bool)] -> [HDR]
densities2HDR68 [(Int, Float, Bool, Bool)]
contextualizedDensities
hdrs95 :: [HDR]
hdrs95 = [(Int, Float, Bool, Bool)] -> [HDR]
densities2HDR95 [(Int, Float, Bool, Bool)]
contextualizedDensities
indexVU :: Vector a -> Maybe Int -> Maybe a
indexVU Vector a
_ Maybe Int
Nothing = Maybe a
forall a. Maybe a
Nothing
indexVU Vector a
x (Just Int
i) = Vector a
x Vector a -> Int -> Maybe a
forall a. Unbox a => Vector a -> Int -> Maybe a
VU.!? Int
i
cumsumDens :: [(YearBCAD, Float)] -> [Float]
cumsumDens :: [(Int, Float)] -> [Float]
cumsumDens [(Int, Float)]
x = (Float -> Float -> Float) -> [Float] -> [Float]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 Float -> Float -> Float
forall a. Num a => a -> a -> a
(+) ([Float] -> [Float]) -> [Float] -> [Float]
forall a b. (a -> b) -> a -> b
$ ((Int, Float) -> Float) -> [(Int, Float)] -> [Float]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Float) -> Float
forall a b. (a, b) -> b
snd [(Int, Float)]
x
densities2HDR68 :: [(Int, Float, Bool, Bool)] -> [HDR]
densities2HDR68 :: [(Int, Float, Bool, Bool)] -> [HDR]
densities2HDR68 [(Int, Float, Bool, Bool)]
cDensities =
let highDensityGroups :: [[(Int, Float, Bool, Bool)]]
highDensityGroups = ((Int, Float, Bool, Bool) -> (Int, Float, Bool, Bool) -> Bool)
-> [(Int, Float, Bool, Bool)] -> [[(Int, Float, Bool, Bool)]]
forall a. (a -> a -> Bool) -> [a] -> [[a]]
groupBy (\(Int
_,Float
_,Bool
in681,Bool
_) (Int
_,Float
_,Bool
in682,Bool
_) -> Bool
in681 Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== Bool
in682) [(Int, Float, Bool, Bool)]
cDensities
filteredDensityGroups :: [[(Int, Float, Bool, Bool)]]
filteredDensityGroups = ([(Int, Float, Bool, Bool)] -> Bool)
-> [[(Int, Float, Bool, Bool)]] -> [[(Int, Float, Bool, Bool)]]
forall a. (a -> Bool) -> [a] -> [a]
filter (((Int, Float, Bool, Bool) -> Bool)
-> [(Int, Float, Bool, Bool)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Int, Float, Bool, Bool) -> Bool
getIn68) [[(Int, Float, Bool, Bool)]]
highDensityGroups
in ([(Int, Float, Bool, Bool)] -> HDR)
-> [[(Int, Float, Bool, Bool)]] -> [HDR]
forall a b. (a -> b) -> [a] -> [b]
map (\[(Int, Float, Bool, Bool)]
xs -> let yearRange :: [Int]
yearRange = ((Int, Float, Bool, Bool) -> Int)
-> [(Int, Float, Bool, Bool)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Float, Bool, Bool) -> Int
getYear [(Int, Float, Bool, Bool)]
xs in Int -> Int -> HDR
HDR ([Int] -> Int
forall a. [a] -> a
head [Int]
yearRange) ([Int] -> Int
forall a. [a] -> a
last [Int]
yearRange)) [[(Int, Float, Bool, Bool)]]
filteredDensityGroups
densities2HDR95 :: [(Int, Float, Bool, Bool)] -> [HDR]
densities2HDR95 :: [(Int, Float, Bool, Bool)] -> [HDR]
densities2HDR95 [(Int, Float, Bool, Bool)]
cDensities =
let highDensityGroups :: [[(Int, Float, Bool, Bool)]]
highDensityGroups = ((Int, Float, Bool, Bool) -> (Int, Float, Bool, Bool) -> Bool)
-> [(Int, Float, Bool, Bool)] -> [[(Int, Float, Bool, Bool)]]
forall a. (a -> a -> Bool) -> [a] -> [[a]]
groupBy (\(Int
_,Float
_,Bool
_,Bool
in951) (Int
_,Float
_,Bool
_,Bool
in952) -> Bool
in951 Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== Bool
in952) [(Int, Float, Bool, Bool)]
cDensities
filteredDensityGroups :: [[(Int, Float, Bool, Bool)]]
filteredDensityGroups = ([(Int, Float, Bool, Bool)] -> Bool)
-> [[(Int, Float, Bool, Bool)]] -> [[(Int, Float, Bool, Bool)]]
forall a. (a -> Bool) -> [a] -> [a]
filter (((Int, Float, Bool, Bool) -> Bool)
-> [(Int, Float, Bool, Bool)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Int, Float, Bool, Bool) -> Bool
getIn95) [[(Int, Float, Bool, Bool)]]
highDensityGroups
in ([(Int, Float, Bool, Bool)] -> HDR)
-> [[(Int, Float, Bool, Bool)]] -> [HDR]
forall a b. (a -> b) -> [a] -> [b]
map (\[(Int, Float, Bool, Bool)]
xs -> let yearRange :: [Int]
yearRange = ((Int, Float, Bool, Bool) -> Int)
-> [(Int, Float, Bool, Bool)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Float, Bool, Bool) -> Int
getYear [(Int, Float, Bool, Bool)]
xs in Int -> Int -> HDR
HDR ([Int] -> Int
forall a. [a] -> a
head [Int]
yearRange) ([Int] -> Int
forall a. [a] -> a
last [Int]
yearRange)) [[(Int, Float, Bool, Bool)]]
filteredDensityGroups
getIn68 :: (Int, Float, Bool, Bool) -> Bool
getIn68 :: (Int, Float, Bool, Bool) -> Bool
getIn68 (Int
_,Float
_,Bool
x,Bool
_) = Bool
x
getIn95 :: (Int, Float, Bool, Bool) -> Bool
getIn95 :: (Int, Float, Bool, Bool) -> Bool
getIn95 (Int
_,Float
_,Bool
_,Bool
x) = Bool
x
getYear :: (Int, Float, Bool, Bool) -> Int
getYear :: (Int, Float, Bool, Bool) -> Int
getYear (Int
year,Float
_,Bool
_,Bool
_) = Int
year