{-# LANGUAGE Strict #-}
module Currycarbon.Calibration.MatrixMult
( calibrateDateMatrixMult
, makeCalCurveMatrix
, uncalToPDF
) where
import Currycarbon.Calibration.Utils
import Currycarbon.Parsers
import Currycarbon.Types
import Currycarbon.Utils
import qualified Data.Vector as V
import Data.Vector.Generic (convert)
import qualified Data.Vector.Unboxed as VU
calibrateDateMatrixMult :: Bool -> Bool -> CalCurveBP -> UncalC14 -> Either CurrycarbonException CalPDF
calibrateDateMatrixMult :: Bool
-> Bool
-> CalCurveBP
-> UncalC14
-> Either CurrycarbonException CalPDF
calibrateDateMatrixMult Bool
allowOutside Bool
interpolate CalCurveBP
calCurve UncalC14
uncalC14 =
if Bool -> Bool
not Bool
allowOutside Bool -> Bool -> Bool
&& CalCurveBP -> UncalC14 -> Bool
isOutsideRangeOfCalCurve CalCurveBP
calCurve UncalC14
uncalC14
then forall a b. a -> Either a b
Left forall a b. (a -> b) -> a -> b
$ String -> CurrycarbonException
CurrycarbonCalibrationRangeException forall a b. (a -> b) -> a -> b
$ UncalC14 -> String
renderUncalC14 UncalC14
uncalC14
else
let rawCalCurveSegment :: CalCurveBP
rawCalCurveSegment = UncalC14 -> CalCurveBP -> CalCurveBP
getRelevantCalCurveSegment UncalC14
uncalC14 CalCurveBP
calCurve
calCurveSegment :: CalCurveBCAD
calCurveSegment = Bool -> CalCurveBP -> CalCurveBCAD
prepareCalCurveSegment Bool
interpolate CalCurveBP
rawCalCurveSegment
uncalPDF :: UncalPDF
uncalPDF = UncalC14 -> UncalPDF
uncalToPDF UncalC14
uncalC14
calCurveMatrix :: CalCurveMatrix
calCurveMatrix = UncalPDF -> CalCurveBCAD -> CalCurveMatrix
makeCalCurveMatrix UncalPDF
uncalPDF CalCurveBCAD
calCurveSegment
calPDF :: CalPDF
calPDF = UncalPDF -> CalCurveMatrix -> CalPDF
projectUncalOverCalCurve UncalPDF
uncalPDF CalCurveMatrix
calCurveMatrix
in forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ CalPDF -> CalPDF
trimLowDensityEdgesCalPDF forall a b. (a -> b) -> a -> b
$ CalPDF -> CalPDF
normalizeCalPDF CalPDF
calPDF
makeCalCurveMatrix :: UncalPDF -> CalCurveBCAD -> CalCurveMatrix
makeCalCurveMatrix :: UncalPDF -> CalCurveBCAD -> CalCurveMatrix
makeCalCurveMatrix (UncalPDF String
_ Vector YearBP
uncals' Vector Float
_) (CalCurveBCAD Vector YearBCAD
cals Vector YearBCAD
uncals Vector YearBP
sigmas) =
let curveUnCalBCADsFloat :: Vector Float
curveUnCalBCADsFloat = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBCAD
uncals
sigmasFloat :: Vector Float
sigmasFloat = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBP
sigmas
uncalBCADs :: Vector YearBCAD
uncalBCADs = Vector YearBP -> Vector YearBCAD
vectorBPToBCAD Vector YearBP
uncals'
uncalBCADsFloat :: Vector Float
uncalBCADsFloat = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBCAD
uncalBCADs
matrix :: Vector (Vector Float)
matrix = Vector Float
-> Vector Float -> Vector Float -> Vector (Vector Float)
buildMatrix Vector Float
curveUnCalBCADsFloat Vector Float
sigmasFloat Vector Float
uncalBCADsFloat
in Vector YearBCAD
-> Vector YearBCAD -> Vector (Vector Float) -> CalCurveMatrix
CalCurveMatrix Vector YearBCAD
uncalBCADs Vector YearBCAD
cals Vector (Vector Float)
matrix
where
buildMatrix :: VU.Vector Float -> VU.Vector Float -> VU.Vector Float -> V.Vector (VU.Vector Float)
buildMatrix :: Vector Float
-> Vector Float -> Vector Float -> Vector (Vector Float)
buildMatrix Vector Float
curveuncal_ Vector Float
sigmas_ Vector Float
uncal_ =
forall a b. (a -> b) -> Vector a -> Vector b
V.map (\(Float, Float)
x -> forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map ((Float, Float) -> Float -> Float
fillCell (Float, Float)
x) Vector Float
uncal_) forall a b. (a -> b) -> a -> b
$
forall a b. Vector a -> Vector b -> Vector (a, b)
V.zip (forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
convert Vector Float
curveuncal_) (forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
convert Vector Float
sigmas_)
fillCell :: (Float, Float) -> Float -> Float
fillCell :: (Float, Float) -> Float -> Float
fillCell (Float
mean, Float
sigma) Float
matrixPosBP =
if forall a. Num a => a -> a
abs (Float
mean forall a. Num a => a -> a -> a
- Float
matrixPosBP) forall a. Ord a => a -> a -> Bool
< Float
6forall a. Num a => a -> a -> a
*Float
sigma
then Float -> Float -> Float -> Float
dnorm Float
mean Float
sigma Float
matrixPosBP
else Float
0
uncalToPDF :: UncalC14 -> UncalPDF
uncalToPDF :: UncalC14 -> UncalPDF
uncalToPDF (UncalC14 String
name YearBP
mean YearBP
std) =
let meanFloat :: Float
meanFloat = forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBP
mean
stdFloat :: Float
stdFloat = forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBP
std
years :: Vector YearBP
years = forall a. Unbox a => Vector a -> Vector a
VU.reverse forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => [a] -> Vector a
VU.fromList [(YearBP
meanforall a. Num a => a -> a -> a
-YearBP
5forall a. Num a => a -> a -> a
*YearBP
std) .. (YearBP
meanforall a. Num a => a -> a -> a
+YearBP
5forall a. Num a => a -> a -> a
*YearBP
std)]
yearsFloat :: Vector Float
yearsFloat = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBP
years
probabilities :: Vector Float
probabilities = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map (Float -> Float -> Float -> Float
dnorm Float
meanFloat Float
stdFloat) Vector Float
yearsFloat
in String -> Vector YearBP -> Vector Float -> UncalPDF
UncalPDF String
name Vector YearBP
years Vector Float
probabilities
projectUncalOverCalCurve :: UncalPDF -> CalCurveMatrix -> CalPDF
projectUncalOverCalCurve :: UncalPDF -> CalCurveMatrix -> CalPDF
projectUncalOverCalCurve (UncalPDF String
name Vector YearBP
_ Vector Float
dens) (CalCurveMatrix Vector YearBCAD
_ Vector YearBCAD
cals Vector (Vector Float)
matrix) =
String -> Vector YearBCAD -> Vector Float -> CalPDF
CalPDF String
name Vector YearBCAD
cals forall a b. (a -> b) -> a -> b
$ Vector Float -> Vector (Vector Float) -> Vector Float
vectorMatrixMultSum Vector Float
dens Vector (Vector Float)
matrix
where
vectorMatrixMultSum :: VU.Vector Float -> V.Vector (VU.Vector Float) -> VU.Vector Float
vectorMatrixMultSum :: Vector Float -> Vector (Vector Float) -> Vector Float
vectorMatrixMultSum Vector Float
vec Vector (Vector Float)
mat =
forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
convert forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> Vector a -> Vector b
V.map (\Vector Float
x -> forall a. (Unbox a, Num a) => Vector a -> a
VU.sum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VU.zipWith forall a. Num a => a -> a -> a
(*) Vector Float
x Vector Float
vec) Vector (Vector Float)
mat