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

-- | Intercept calibration implemented with matrix multiplication (see 'MatrixMultiplication')
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

-- | Construct a matrix representation of a calibration curve for a given date
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

-- | Transform an uncalibrated date to an uncalibrated
-- probability density table
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