{-# LANGUAGE Strict #-}
module Currycarbon.Calibration.Calibration
(
getRelevantCalCurveSegment
, prepareCalCurveSegment
, makeCalCurveMatrix
, uncalToPDF
, calibrateDate
, calibrateDates
, refineCalDates
, refineCalDate
, CalibrateDatesConf (..)
, defaultCalConf
, AgeSamplingConf (..)
, sampleAgesFromCalPDF
) where
import Currycarbon.Calibration.Bchron
import Currycarbon.Calibration.MatrixMult
import Currycarbon.Calibration.Utils
import Currycarbon.Types
import Currycarbon.Utils
import qualified Control.Monad.Random as CMR
import Data.List (elemIndex, groupBy, sort,
sortBy)
import Data.Maybe (fromJust)
import qualified Data.Vector.Unboxed as VU
import qualified System.Random as R
data CalibrateDatesConf = CalibrateDatesConf {
CalibrateDatesConf -> CalibrationMethod
_calConfMethod :: CalibrationMethod
, CalibrateDatesConf -> Bool
_calConfAllowOutside :: Bool
, CalibrateDatesConf -> Bool
_calConfInterpolateCalCurve :: Bool
} deriving (YearBCAD -> CalibrateDatesConf -> ShowS
[CalibrateDatesConf] -> ShowS
CalibrateDatesConf -> String
forall a.
(YearBCAD -> a -> ShowS)
-> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [CalibrateDatesConf] -> ShowS
$cshowList :: [CalibrateDatesConf] -> ShowS
show :: CalibrateDatesConf -> String
$cshow :: CalibrateDatesConf -> String
showsPrec :: YearBCAD -> CalibrateDatesConf -> ShowS
$cshowsPrec :: YearBCAD -> CalibrateDatesConf -> ShowS
Show, CalibrateDatesConf -> CalibrateDatesConf -> Bool
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 {
_calConfMethod :: CalibrationMethod
_calConfMethod = 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 =
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 =
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] -> [Either CurrycarbonException CalC14]
refineCalDates :: [CalPDF] -> [Either CurrycarbonException CalC14]
refineCalDates = forall a b. (a -> b) -> [a] -> [b]
map CalPDF -> Either CurrycarbonException CalC14
refineCalDate
refineCalDate :: CalPDF -> Either CurrycarbonException CalC14
refineCalDate :: CalPDF -> Either CurrycarbonException CalC14
refineCalDate calPDF :: CalPDF
calPDF@(CalPDF String
name Vector YearBCAD
cals Vector Float
dens)
| CalPDF -> Bool
isInvalidCalPDF CalPDF
calPDF =
forall a b. a -> Either a b
Left forall a b. (a -> b) -> a -> b
$ String -> CurrycarbonException
CurrycarbonInvalidCalPDFException String
"refinement"
| forall a. Unbox a => Vector a -> YearBCAD
VU.length (forall a. (Unbox a, Eq a) => Vector a -> Vector a
VU.uniq Vector Float
dens) forall a. Eq a => a -> a -> Bool
== YearBCAD
1 =
let start :: YearBCAD
start = forall a. Unbox a => Vector a -> a
VU.head Vector YearBCAD
cals
stop :: YearBCAD
stop = forall a. Unbox a => Vector a -> a
VU.last Vector YearBCAD
cals
in forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ CalC14 {
_calC14id :: String
_calC14id = String
name
, _calC14RangeSummary :: CalRangeSummary
_calC14RangeSummary = CalRangeSummary {
_calRangeStartTwoSigma :: YearBCAD
_calRangeStartTwoSigma = YearBCAD
start
, _calRangeStartOneSigma :: YearBCAD
_calRangeStartOneSigma = YearBCAD
start
, _calRangeMedian :: YearBCAD
_calRangeMedian = YearBCAD
median
, _calRangeStopOneSigma :: YearBCAD
_calRangeStopOneSigma = YearBCAD
stop
, _calRangeStopTwoSigma :: YearBCAD
_calRangeStopTwoSigma = YearBCAD
stop
}
, _calC14HDROneSigma :: [HDR]
_calC14HDROneSigma = [YearBCAD -> YearBCAD -> HDR
HDR YearBCAD
start YearBCAD
stop]
, _calC14HDRTwoSigma :: [HDR]
_calC14HDRTwoSigma = [YearBCAD -> YearBCAD -> HDR
HDR YearBCAD
start YearBCAD
stop]
}
| Bool
otherwise =
forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ CalC14 {
_calC14id :: String
_calC14id = String
name
, _calC14RangeSummary :: CalRangeSummary
_calC14RangeSummary = CalRangeSummary {
_calRangeStartTwoSigma :: YearBCAD
_calRangeStartTwoSigma = HDR -> YearBCAD
_hdrstart forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
head [HDR]
hdrs95
, _calRangeStartOneSigma :: YearBCAD
_calRangeStartOneSigma = HDR -> YearBCAD
_hdrstart forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
head [HDR]
hdrs68
, _calRangeMedian :: YearBCAD
_calRangeMedian = YearBCAD
median
, _calRangeStopOneSigma :: YearBCAD
_calRangeStopOneSigma = HDR -> YearBCAD
_hdrstop forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
last [HDR]
hdrs68
, _calRangeStopTwoSigma :: YearBCAD
_calRangeStopTwoSigma = HDR -> YearBCAD
_hdrstop forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
last [HDR]
hdrs95
}
, _calC14HDROneSigma :: [HDR]
_calC14HDROneSigma = [HDR]
hdrs68
, _calC14HDRTwoSigma :: [HDR]
_calC14HDRTwoSigma = [HDR]
hdrs95
}
where
cumsumDensities :: [Float]
cumsumDensities = [(YearBCAD, Float)] -> [Float]
cumsumDens (forall a. Unbox a => Vector a -> [a]
VU.toList forall a b. (a -> b) -> a -> b
$ forall a b.
(Unbox a, Unbox b) =>
Vector a -> Vector b -> Vector (a, b)
VU.zip Vector YearBCAD
cals Vector Float
dens)
distanceTo05 :: [Float]
distanceTo05 = forall a b. (a -> b) -> [a] -> [b]
map (\Float
x -> forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ (Float
x forall a. Num a => a -> a -> a
- Float
0.5)) [Float]
cumsumDensities
median :: YearBCAD
median = forall a. HasCallStack => Maybe a -> a
fromJust forall a b. (a -> b) -> a -> b
$ Vector YearBCAD
cals forall {a}. Unbox a => Vector a -> Maybe YearBCAD -> Maybe a
`indexVU` forall a. Eq a => a -> [a] -> Maybe YearBCAD
elemIndex (forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Float]
distanceTo05) [Float]
distanceTo05
sortedDensities :: [(YearBCAD, Float)]
sortedDensities = forall a. (a -> a -> Ordering) -> [a] -> [a]
sortBy (forall a b c. (a -> b -> c) -> b -> a -> c
flip (\ (YearBCAD
_, Float
dens1) (YearBCAD
_, Float
dens2) -> forall a. Ord a => a -> a -> Ordering
compare Float
dens1 Float
dens2)) (forall a. Unbox a => Vector a -> [a]
VU.toList forall a b. (a -> b) -> a -> b
$ forall a b.
(Unbox a, Unbox b) =>
Vector a -> Vector b -> Vector (a, b)
VU.zip Vector YearBCAD
cals Vector Float
dens)
cumsumSortedDensities :: [Float]
cumsumSortedDensities = [(YearBCAD, Float)] -> [Float]
cumsumDens [(YearBCAD, Float)]
sortedDensities
isIn68 :: [Bool]
isIn68 = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Ord a => a -> a -> Bool
< Float
0.683) [Float]
cumsumSortedDensities
isIn95 :: [Bool]
isIn95 = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Ord a => a -> a -> Bool
< Float
0.954) [Float]
cumsumSortedDensities
contextualizedDensities :: [(YearBCAD, Float, Bool, Bool)]
contextualizedDensities = forall a. Ord a => [a] -> [a]
sort forall a b. (a -> b) -> a -> b
$ forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\(YearBCAD
y,Float
d) Bool
in68 Bool
in95 -> (YearBCAD
y,Float
d,Bool
in68,Bool
in95)) [(YearBCAD, Float)]
sortedDensities [Bool]
isIn68 [Bool]
isIn95
hdrs68 :: [HDR]
hdrs68 = [(YearBCAD, Float, Bool, Bool)] -> [HDR]
densities2HDR68 [(YearBCAD, Float, Bool, Bool)]
contextualizedDensities
hdrs95 :: [HDR]
hdrs95 = [(YearBCAD, Float, Bool, Bool)] -> [HDR]
densities2HDR95 [(YearBCAD, Float, Bool, Bool)]
contextualizedDensities
indexVU :: Vector a -> Maybe YearBCAD -> Maybe a
indexVU Vector a
_ Maybe YearBCAD
Nothing = forall a. Maybe a
Nothing
indexVU Vector a
x (Just YearBCAD
i) = Vector a
x forall a. Unbox a => Vector a -> YearBCAD -> Maybe a
VU.!? YearBCAD
i
cumsumDens :: [(YearBCAD, Float)] -> [Float]
cumsumDens :: [(YearBCAD, Float)] -> [Float]
cumsumDens [(YearBCAD, Float)]
x = forall a. (a -> a -> a) -> [a] -> [a]
scanl1 forall a. Num a => a -> a -> a
(+) forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd [(YearBCAD, Float)]
x
densities2HDR68 :: [(Int, Float, Bool, Bool)] -> [HDR]
densities2HDR68 :: [(YearBCAD, Float, Bool, Bool)] -> [HDR]
densities2HDR68 [(YearBCAD, Float, Bool, Bool)]
cDensities =
let highDensityGroups :: [[(YearBCAD, Float, Bool, Bool)]]
highDensityGroups = forall a. (a -> a -> Bool) -> [a] -> [[a]]
groupBy (\(YearBCAD
_,Float
_,Bool
in681,Bool
_) (YearBCAD
_,Float
_,Bool
in682,Bool
_) -> Bool
in681 forall a. Eq a => a -> a -> Bool
== Bool
in682) [(YearBCAD, Float, Bool, Bool)]
cDensities
filteredDensityGroups :: [[(YearBCAD, Float, Bool, Bool)]]
filteredDensityGroups = forall a. (a -> Bool) -> [a] -> [a]
filter (forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (YearBCAD, Float, Bool, Bool) -> Bool
getIn68) [[(YearBCAD, Float, Bool, Bool)]]
highDensityGroups
in forall a b. (a -> b) -> [a] -> [b]
map (\[(YearBCAD, Float, Bool, Bool)]
xs -> let yearRange :: [YearBCAD]
yearRange = forall a b. (a -> b) -> [a] -> [b]
map (YearBCAD, Float, Bool, Bool) -> YearBCAD
getYear [(YearBCAD, Float, Bool, Bool)]
xs in YearBCAD -> YearBCAD -> HDR
HDR (forall a. [a] -> a
head [YearBCAD]
yearRange) (forall a. [a] -> a
last [YearBCAD]
yearRange)) [[(YearBCAD, Float, Bool, Bool)]]
filteredDensityGroups
densities2HDR95 :: [(Int, Float, Bool, Bool)] -> [HDR]
densities2HDR95 :: [(YearBCAD, Float, Bool, Bool)] -> [HDR]
densities2HDR95 [(YearBCAD, Float, Bool, Bool)]
cDensities =
let highDensityGroups :: [[(YearBCAD, Float, Bool, Bool)]]
highDensityGroups = forall a. (a -> a -> Bool) -> [a] -> [[a]]
groupBy (\(YearBCAD
_,Float
_,Bool
_,Bool
in951) (YearBCAD
_,Float
_,Bool
_,Bool
in952) -> Bool
in951 forall a. Eq a => a -> a -> Bool
== Bool
in952) [(YearBCAD, Float, Bool, Bool)]
cDensities
filteredDensityGroups :: [[(YearBCAD, Float, Bool, Bool)]]
filteredDensityGroups = forall a. (a -> Bool) -> [a] -> [a]
filter (forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (YearBCAD, Float, Bool, Bool) -> Bool
getIn95) [[(YearBCAD, Float, Bool, Bool)]]
highDensityGroups
in forall a b. (a -> b) -> [a] -> [b]
map (\[(YearBCAD, Float, Bool, Bool)]
xs -> let yearRange :: [YearBCAD]
yearRange = forall a b. (a -> b) -> [a] -> [b]
map (YearBCAD, Float, Bool, Bool) -> YearBCAD
getYear [(YearBCAD, Float, Bool, Bool)]
xs in YearBCAD -> YearBCAD -> HDR
HDR (forall a. [a] -> a
head [YearBCAD]
yearRange) (forall a. [a] -> a
last [YearBCAD]
yearRange)) [[(YearBCAD, Float, Bool, Bool)]]
filteredDensityGroups
getIn68 :: (Int, Float, Bool, Bool) -> Bool
getIn68 :: (YearBCAD, Float, Bool, Bool) -> Bool
getIn68 (YearBCAD
_,Float
_,Bool
x,Bool
_) = Bool
x
getIn95 :: (Int, Float, Bool, Bool) -> Bool
getIn95 :: (YearBCAD, Float, Bool, Bool) -> Bool
getIn95 (YearBCAD
_,Float
_,Bool
_,Bool
x) = Bool
x
getYear :: (Int, Float, Bool, Bool) -> Int
getYear :: (YearBCAD, Float, Bool, Bool) -> YearBCAD
getYear (YearBCAD
year,Float
_,Bool
_,Bool
_) = YearBCAD
year
data AgeSamplingConf = AgeSamplingConf {
AgeSamplingConf -> StdGen
_assRNG :: R.StdGen
, AgeSamplingConf -> Word
_assNumberOfSamples :: Word
} deriving (YearBCAD -> AgeSamplingConf -> ShowS
[AgeSamplingConf] -> ShowS
AgeSamplingConf -> String
forall a.
(YearBCAD -> a -> ShowS)
-> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [AgeSamplingConf] -> ShowS
$cshowList :: [AgeSamplingConf] -> ShowS
show :: AgeSamplingConf -> String
$cshow :: AgeSamplingConf -> String
showsPrec :: YearBCAD -> AgeSamplingConf -> ShowS
$cshowsPrec :: YearBCAD -> AgeSamplingConf -> ShowS
Show, AgeSamplingConf -> AgeSamplingConf -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: AgeSamplingConf -> AgeSamplingConf -> Bool
$c/= :: AgeSamplingConf -> AgeSamplingConf -> Bool
== :: AgeSamplingConf -> AgeSamplingConf -> Bool
$c== :: AgeSamplingConf -> AgeSamplingConf -> Bool
Eq)
sampleAgesFromCalPDF :: AgeSamplingConf -> CalPDF -> Either CurrycarbonException RandomAgeSample
sampleAgesFromCalPDF :: AgeSamplingConf
-> CalPDF -> Either CurrycarbonException RandomAgeSample
sampleAgesFromCalPDF (AgeSamplingConf StdGen
rng Word
n) calPDF :: CalPDF
calPDF@(CalPDF String
calPDFid Vector YearBCAD
cals Vector Float
dens) =
let weightedList :: [(YearBCAD, Rational)]
weightedList = forall a b. [a] -> [b] -> [(a, b)]
zip (forall a. Unbox a => Vector a -> [a]
VU.toList Vector YearBCAD
cals) (forall a b. (a -> b) -> [a] -> [b]
map forall a. Real a => a -> Rational
toRational forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Vector a -> [a]
VU.toList Vector Float
dens)
infSamplesList :: [YearBCAD]
infSamplesList = forall g a. RandomGen g => g -> [(a, Rational)] -> [a]
sampleWeightedList StdGen
rng [(YearBCAD, Rational)]
weightedList
samples :: [YearBCAD]
samples = forall a. YearBCAD -> [a] -> [a]
take (forall a b. (Integral a, Num b) => a -> b
fromIntegral Word
n) [YearBCAD]
infSamplesList
in if CalPDF -> Bool
isInvalidCalPDF CalPDF
calPDF
then forall a b. a -> Either a b
Left forall a b. (a -> b) -> a -> b
$ String -> CurrycarbonException
CurrycarbonInvalidCalPDFException String
"random age sampling"
else forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ String -> Vector YearBCAD -> RandomAgeSample
RandomAgeSample String
calPDFid (forall a. Unbox a => [a] -> Vector a
VU.fromList [YearBCAD]
samples)
where
sampleWeightedList :: CMR.RandomGen g => g -> [(a, Rational)] -> [a]
sampleWeightedList :: forall g a. RandomGen g => g -> [(a, Rational)] -> [a]
sampleWeightedList g
gen [(a, Rational)]
weights = forall g a. Rand g a -> g -> a
CMR.evalRand RandT g Identity [a]
m g
gen
where m :: RandT g Identity [a]
m = forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. a -> [a]
repeat forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall (m :: * -> *) a. MonadRandom m => [(a, Rational)] -> m a
CMR.fromList forall a b. (a -> b) -> a -> b
$ [(a, Rational)]
weights
isInvalidCalPDF :: CalPDF -> Bool
isInvalidCalPDF :: CalPDF -> Bool
isInvalidCalPDF (CalPDF String
_ Vector YearBCAD
_ Vector Float
dens) = forall a. (Unbox a, Num a) => Vector a -> a
VU.sum Vector Float
dens forall a. Eq a => a -> a -> Bool
== Float
0 Bool -> Bool -> Bool
|| forall a. Unbox a => (a -> Bool) -> Vector a -> Bool
VU.any (forall a. Ord a => a -> a -> Bool
>= Float
1.0) Vector Float
dens