module MultidimensionalOptimizationMethods.Newton_sMethod

open DiffSharp
open FsAlg.Generic

open Types

let run func initial_x (targetError: float) =
    let inverseHessian x = (dsharp.hessian func x).inv ()

    let rec loop x_previous i =
        let gradient = dsharp.grad func x_previous
        let actualError = gradient.toArray1D () |> Vector.ofArray |> Vector.l2norm

        if targetError.CompareTo(actualError) > 0 then
            { ResultTensor = x_previous
              StepsCount = i
              Error = actualError }

        else
            let x =
                x_previous
                - dsharp.matmul ((x_previous |> inverseHessian), (dsharp.grad func x_previous))

            loop x (i + 1)

    loop initial_x 0
