// File Library.fs

module NumericalIntegrator

type XYdata =
    {
        X : float
        Y : float list
    }

let listToXYdata theList =
    {X = List.head theList; Y = List.tail theList}

let step theDiffEq delta (initialData:XYdata) =
    let rec helper newY oldY =
        match oldY with
        | y1::y2::tail   ->  helper
                                ( (y1 + delta * y2) :: newY )
                                (y2 :: tail)
        | [yLast]       ->  (yLast + delta * (theDiffEq initialData)) :: newY
                            |> List.rev
        | []            ->  failwith "error"
    {X = initialData.X + delta; Y = helper [] initialData.Y }


let stepMany (numSteps:uint) theDiffEq delta initialConditions =
    let rec helper history (iterationsLeft:uint) =
        match iterationsLeft with
        | 0u -> history
        | iterationsLeft -> helper (step theDiffEq delta history.Head :: history) (iterationsLeft - 1u)
    helper [initialConditions] numSteps
    |> List.rev


// File Tests.fs

module Tests

open System
open Xunit
open Swensen.Unquote

open NumericalIntegrator

let isWithinToleranceOf expectedFunc tolerance (actual:XYdata list) =
    let helper index (data:XYdata) =
        let expected = expectedFunc data.X
        try
            Assert.InRange(data.Y.Head, expected - tolerance, expected + tolerance)
        with
            | :? Sdk.InRangeException ->
                let message =
                    sprintf
                        "  -----\nAt iteration %d exceeded tolerance of %f\nExpected %f\nActual %f\n     -----"
                        index tolerance expected data.Y.Head
                Assert.True(false, message)

    actual
    |> List.iteri helper
    |> ignore

let computationTime numSteps theDiffEq delta initialConditions =
    let crono = Diagnostics.Stopwatch.StartNew()
    stepMany numSteps theDiffEq delta initialConditions
    |> ignore
    crono.Stop()
    crono.ElapsedMilliseconds

let computationTimeMany howMany numSteps theDiffEq delta initialConditions =
    List.init howMany
        (fun _ -> computationTime numSteps theDiffEq delta initialConditions)

let benchmark theDiffEq =
    computationTimeMany 20 100000u theDiffEq 0.01 {X = 0.0; Y = [1.0]}
    |> List.averageBy float



[<Fact>]
let ``Test single step of y1=1`` () =
    let theDiffEq xy = 1.0
    let delta = 0.1
    let initialConditions = listToXYdata [0.0; 0.0]

    listToXYdata [delta; delta] =! step theDiffEq delta initialConditions

[<Fact>]
let ``Test stepMany with y1=1`` () =
    let theDiffEq xy = 1.0
    let delta = 0.1
    let initialConditions = listToXYdata [0.0; 0.0]
    let steps = 2u

    let expected =
        [ [0.0; 0.0]; [delta; delta]; [delta*2.0; delta*2.0] ]
        |> List.map listToXYdata

    [initialConditions] =! stepMany 0u theDiffEq delta initialConditions
    expected =! stepMany steps theDiffEq delta initialConditions

[<Fact>]
let ``y1 = y`` () =
    let theDiffEq (xy:XYdata) = xy.Y.[0]
    let delta = 5e-4
    let initialConditions = listToXYdata [0.0; 1.0]
    let expectedFunc x = Math.Exp x
    let numSteps = 200u
    let tolerance = 1e-4

    let actualResult = stepMany numSteps theDiffEq delta initialConditions

    actualResult
    |> isWithinToleranceOf expectedFunc tolerance

[<Fact>]
let ``y2 = -y`` () =
    let theDiffEq xy = -xy.Y.[0]
    let delta = 5e-4
    let initialConditions = listToXYdata [0.0; 1.0; 0.0]
    let expectedFunc x = Math.Cos x
    let numSteps = 200u
    let tolerance = 1e-4

    let actualResult = stepMany numSteps theDiffEq delta initialConditions

    actualResult
    |> isWithinToleranceOf expectedFunc tolerance

[<Fact>]
let ``Benchmark`` () =
    let theDiffEq xy = -xy.Y.[0]
    let averageMilliseconds = benchmark theDiffEq
    averageMilliseconds <! 60.0
    printfn "\n -- BENCHMARK --\nTook %f ms\n ---------------" averageMilliseconds