{-# 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.Unboxed as VU
import qualified Data.Vector as V
import Data.Vector.Generic (convert)

-- | 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 CurrycarbonException -> Either CurrycarbonException CalPDF
forall a b. a -> Either a b
Left (CurrycarbonException -> Either CurrycarbonException CalPDF)
-> CurrycarbonException -> Either CurrycarbonException CalPDF
forall a b. (a -> b) -> a -> b
$ String -> CurrycarbonException
CurrycarbonCalibrationRangeException (String -> CurrycarbonException) -> String -> CurrycarbonException
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 CalPDF -> Either CurrycarbonException CalPDF
forall a b. b -> Either a b
Right (CalPDF -> Either CurrycarbonException CalPDF)
-> CalPDF -> Either CurrycarbonException CalPDF
forall a b. (a -> b) -> a -> b
$ CalPDF -> CalPDF
trimLowDensityEdgesCalPDF (CalPDF -> CalPDF) -> CalPDF -> CalPDF
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 = (YearBCAD -> Float) -> Vector YearBCAD -> Vector Float
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map YearBCAD -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBCAD
uncals
        sigmasFloat :: Vector Float
sigmasFloat = (YearBP -> Float) -> Vector YearBP -> Vector Float
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map YearBP -> Float
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 = (YearBCAD -> Float) -> Vector YearBCAD -> Vector Float
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map YearBCAD -> Float
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_ =
          ((Float, Float) -> Vector Float)
-> Vector (Float, Float) -> Vector (Vector Float)
forall a b. (a -> b) -> Vector a -> Vector b
V.map (\(Float, Float)
x -> (Float -> Float) -> Vector Float -> Vector Float
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_) (Vector (Float, Float) -> Vector (Vector Float))
-> Vector (Float, Float) -> Vector (Vector Float)
forall a b. (a -> b) -> a -> b
$ 
            Vector Float -> Vector Float -> Vector (Float, Float)
forall a b. Vector a -> Vector b -> Vector (a, b)
V.zip (Vector Float -> Vector Float
forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
convert Vector Float
curveuncal_) (Vector Float -> Vector Float
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 Float -> Float
forall a. Num a => a -> a
abs (Float
mean Float -> Float -> Float
forall a. Num a => a -> a -> a
- Float
matrixPosBP) Float -> Float -> Bool
forall a. Ord a => a -> a -> Bool
< Float
6Float -> Float -> Float
forall 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 = YearBP -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBP
mean
        stdFloat :: Float
stdFloat = YearBP -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral YearBP
std
        years :: Vector YearBP
years = Vector YearBP -> Vector YearBP
forall a. Unbox a => Vector a -> Vector a
VU.reverse (Vector YearBP -> Vector YearBP) -> Vector YearBP -> Vector YearBP
forall a b. (a -> b) -> a -> b
$ [YearBP] -> Vector YearBP
forall a. Unbox a => [a] -> Vector a
VU.fromList [(YearBP
meanYearBP -> YearBP -> YearBP
forall a. Num a => a -> a -> a
-YearBP
5YearBP -> YearBP -> YearBP
forall a. Num a => a -> a -> a
*YearBP
std) .. (YearBP
meanYearBP -> YearBP -> YearBP
forall a. Num a => a -> a -> a
+YearBP
5YearBP -> YearBP -> YearBP
forall a. Num a => a -> a -> a
*YearBP
std)]
        yearsFloat :: Vector Float
yearsFloat = (YearBP -> Float) -> Vector YearBP -> Vector Float
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
VU.map YearBP -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral Vector YearBP
years
        probabilities :: Vector Float
probabilities = (Float -> Float) -> Vector Float -> Vector Float
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 (Vector Float -> CalPDF) -> Vector Float -> CalPDF
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 =
          Vector Float -> Vector Float
forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
convert (Vector Float -> Vector Float) -> Vector Float -> Vector Float
forall a b. (a -> b) -> a -> b
$ (Vector Float -> Float) -> Vector (Vector Float) -> Vector Float
forall a b. (a -> b) -> Vector a -> Vector b
V.map (\Vector Float
x -> Vector Float -> Float
forall a. (Unbox a, Num a) => Vector a -> a
VU.sum (Vector Float -> Float) -> Vector Float -> Float
forall a b. (a -> b) -> a -> b
$ (Float -> Float -> Float)
-> Vector Float -> Vector Float -> Vector Float
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
VU.zipWith Float -> Float -> Float
forall a. Num a => a -> a -> a
(*) Vector Float
x Vector Float
vec) Vector (Vector Float)
mat