Merge branch 'master' of https://github.com/HugoGranstrom/numericalnim

Author: HugoGranstrom

.gitignore

@@ -2,3 +2,6 @@ nimcache/
 *.exe
 bin/
 .vscode/
+*.code-workspace
+*.html
+!*.*

README.md

@@ -264,12 +264,58 @@ proc adaptiveGauss*[T](f: proc(x: float, optional: seq[T]): T, xStart, xEnd: flo
 ```
 If you don't understand what the "T" stands for, you can replace it with "float" in your head and read up on "Generics" in Nim.
 
+# Optimization
+## Optimization methods
+## 1 dimensional function optimization
+So far only a few methods have been implemented:
+
+### One Dimensional optimization methods
+- `steepest_descent` - Standard method for local minimum finding over a 2D plane
+- `conjugate_gradient` - iterative implementation of solving Ax = b
+- `newtons` - Newton-Raphson implementation for 1-dimensional functions
+
+## Usage
+Using default parameters the methods need 3 things: the function f(t, y) = y'(t, y), the initial values and the timepoints that you want the solution y(t) at.
+
+## Quick Tutorial
+
+Say we have some differentiable function and we would like to find one of its roots
+
+f = $\frac{1}{3}$x$^{3}$ - 2x$^{2}$ + 3x
+
+$\frac{df}{dx}$ = x$^{2}$ - 4x + 3
+
+
+
+If we translate this to code we get:
+
+```nim
+import math
+import numericalnim
+
+proc f(x:float64): float64 = (1.0 / 3.0) * x ^ 3 - 2 * x ^ 2 + 3 * x
+proc df(x:float64): float64 = x ^ 2 - 4 * x + 3
+```
+
+now given a starting point (and optional precision) we can estimate a nearby root
+We know for this function our actual root is 0
+
+```nim
+import numericalnim
+var start = 0.5
+result = newtons(f, df, start)
+echo result
+
+-1.210640218782444e-23
+```
+Pretty close!
 
 # Utils
 I have included a few handy tools in `numericalnim/utils`.
 ## Vector
 Hurray! Yet another vector library! This was mostly done for my own use but I figured it could come in handy if one wanted to just throw something together. It's main purpose is to enable you to solve systems of ODEs using your own types (for example arbitrary precision numbers). The `Vector` type is just a glorified seq with operator overload. No vectorization (unless the compiler does it automatically) sadly. Maybe can get OpenMP to work (or you maybe you, dear reader, can fix it :wink).
 The following operators and procs are supported:
+
 - `+` : Addition between `Vector`s and floats.
 - `-` : Addition between `Vector`s and floats.
 - `+=`, `-=` : inplace addition and subtraction.
@@ -281,10 +327,12 @@ The following operators and procs are supported:
 - `.*=`, `./=` : inplace elementwise multiplication and division between `Vector`s. (not nested `Vector`s)
 - `-` : negation (-Vector).
 - `dot` : Same as `*` between `Vector`s. It is recursive so it will not be a matrix dot product if nested `Vector`s are used.
-- `abs` : The magnitude of the `Vector`. It is recursive so it may not be one of the usual norms.
+- `abs` : The magnitude of the `Vector`. It is recursive so it may not be one of the usual norms. Equivalent to norm(`Vector`, 2)
 - `[]` : Use `v[i]` to get the i:th element of the `Vector`.
 - `==` : Compares two `Vector`s to see if they are equal.
 - `@` : Unpacks the Vector to (nested) seqs. Works with 1, 2 and 3 dimensional Vectors. 
+- `^` : Element-wise exponentiation, works with natural and floating point powers, returns a new Vector object
+- `norm` : General vector norm function
 
 A `Vector` is created using the `newVector` proc and is passed an `openArray` of the elements:
 ```nim
@@ -340,3 +388,4 @@ If you want to use Arraymancer with NumericalNim, the basics should work but I h
 - Add more integration methods (Gaussian Quadrature on it's way. Done).
 - Make the existing code more efficient and robust.
 - Add parallelization of some kind to speed it up. `Vector` would probably benefit from it.
+- More optimization methods!

src/numericalnim.nim

@@ -4,3 +4,5 @@ import numericalnim/integrate
 export integrate
 import numericalnim/ode
 export ode
+import numericalnim/optimize
+export optimize

src/numericalnim/optimize.nim

@@ -0,0 +1,108 @@
+import strformat
+import arraymancer
+import sequtils
+import math
+
+proc steepest_descent*(deriv: proc(x: float64): float64, start: float64, gamma: float64 = 0.01, precision: float64 = 1e-5, max_iters: Natural = 1000):float64 {.inline.} =
+    ## Gradient descent optimization algorithm for finding local minimums of a function with derivative 'deriv'
+    ## 
+    ## Assuming that a multivariable function F is defined and differentiable near a minimum, F(x) decreases fastest
+    ## when going in the direction negative to the gradient of F(a), similar to how water might traverse down a hill 
+    ## following the path of least resistance.
+    ## can benefit from preconditioning if the condition number of the coefficient matrix is ill-conditioned
+    ## Input:
+    ##   - deriv: derivative of a multivariable function F
+    ##   - start: starting point near F's minimum
+    ##   - gamma: step size multiplier, used to control the step size between iterations
+    ##   - precision: numerical precision
+    ##   - max_iters: maximum iterations
+    ##
+    ## Returns:
+    ##   - float64.
+    var
+        current = 0.0
+        x = start
+
+    for i in 0 .. max_iters:
+        # calculate the next direction to propogate
+        current = x
+        x = current - gamma * deriv(current)
+        
+        # If we haven't moved much since the last iteration, break
+        if abs(x - current) <= precision:
+            break
+
+        if i == max_iters:
+            raise newException(ArithmeticError, "Maximum iterations for Steepest descent method exceeded")
+
+    return x
+
+proc conjugate_gradient*[T](A, b, x_0: Tensor[T], tolerance: float64): Tensor[T] =
+    ## Conjugate Gradient method.
+    ## Given a Symmetric and Positive-Definite matrix A, solve the linear system Ax = b
+    ## Symmetric Matrix: Square matrix that is equal to its transpose, transpose(A) == A
+    ## Positive Definite: Square matrix such that transpose(x)Ax > 0 for all x in R^n
+    ## 
+    ## Input:
+    ##   - A: NxN square matrix
+    ##   - b: vector on the right side of Ax=b
+    ##   - x_0: Initial guess vector
+    ##
+    ## Returns:
+    ##   - Tensor.
+    
+    var r = b - (A * x_0)
+    var p = r
+    var rsold = (r.transpose() * r)[0,0] # multiplication returns a Tensor, so need the first element
+    
+    result = x_0
+    
+    var
+        Ap = A
+        alpha = 0.0
+        rsnew = 0.0
+        Ap_p = 0.0
+
+    for i in 1 .. b.shape[0]:
+        Ap = A * p
+        Ap_p = (p.transpose() * Ap)[0,0] 
+        alpha = rsold / Ap_p
+        result = result + alpha * p
+        r = r - alpha * Ap
+        rsnew = (r.transpose() * r)[0,0]
+        if sqrt(rsnew) < tolerance:
+            break
+        p = r + (rsnew / rsold) * p
+        rsold = rsnew
+    
+
+proc newtons*(f: proc(x: float64): float64, deriv: proc(x: float64): float64, start: float64, precision: float64 = 1e-5, max_iters: Natural = 1000): float64 {.raises: [ArithmeticError].} =
+    ## Newton-Raphson implementation for 1-dimensional functions
+
+    ## Given a single variable function f and it's derivative, calcuate an approximation to f(x) = 0
+    ## Input:
+    ##   - f: "Well behaved" function of a single variable with a known root
+    ##   - deriv: derivative of f with respect to x
+    ##   - start: starting x 
+    ##   - precision: numerical precision
+    ##   - max_iters: maxmimum number of iterations
+    ##
+    ## Returns:
+    ##   - float64.
+    var 
+        x_iter = start
+        i = 0
+        current_f = f(start)
+
+    while abs(current_f) >= precision and i <= max_iters:
+        current_f = f(x_iter)
+        x_iter = x_iter - (current_f / deriv(x_iter))
+        i += 1
+        if i == max_iters:
+            raise newException(ArithmeticError, "Maximum iterations for Newtons method exceeded")
+
+    return x_iter - (current_f / deriv(x_iter))
+
+
+
+

src/numericalnim/utils.nim

@@ -192,6 +192,20 @@ proc abs*[T](v1: Vector[T]): float {.inline.} =
     result = sqrt(result)
 proc mean_squared_error*[T](v1: Vector[T]): float {.inline.} = abs(v1) / v1.len.toFloat
 
+proc `^`*[float](v: Vector[float], power: Natural): Vector[float] {.inline.} =
+    ## Returns a Vector object after raising each element to a power (Natural number powers)
+    var newComponents = newSeq[float](v.len)
+    for i in 0 .. v.components.high:
+        newComponents[i] = v[i] ^ power
+    result = newVector(newComponents)
+
+proc `^`*[float](v: Vector[float], power: float): Vector[float] {.inline.} =
+    ## Returns a Vector object after raising each element to a power (float powers)
+    var newComponents = newSeq[float](v.len)
+    for i in 0 .. v.components.high:
+        newComponents[i] = pow(v[i], power)
+    result = newVector(newComponents)
+ 
 
 proc clone*[T](x: T): T {.inline.} = x
 proc mean_squared_error*[T](y_true, y: T): float {.inline.} = abs(y_true - y)
@@ -206,6 +220,7 @@ proc hermiteSpline*[T](x, x1, x2: float, y1, y2, dy1, dy2: T): T {.inline.}=
     let h11 = t ^ 3 - t ^ 2
     result = h00 * y1 + h10 * (x2 - x1) * dy1 + h01 * y2 + h11 * (x2 - x1) * dy2
 
+
 proc hermiteInterpolate*[T](x: openArray[float], t: openArray[float],
                             y, dy: openArray[T]): seq[T] {.inline.} =
     # loop over each interval and check if x is in there, if x is sorted
@@ -237,6 +252,8 @@ proc hermiteInterpolate*[T](x: openArray[float], t: openArray[float],
 
 
 
+
+
 proc sortDataset*[T](X: openArray[float], Y: openArray[T]): seq[(float, T)] {.inline.} =
     if X.len != Y.len:
         raise newException(ValueError, "X and Y must have the same length")
@@ -252,6 +269,7 @@ proc isClose*[T](y1, y2: T, tol: float = 1e-3): bool {.inline.} =
     else:
         return false
 
+
 proc arange*(x1, x2, dx: float, includeStart = true, includeEnd = false): seq[float] {.inline.} =
     let dx = abs(dx) * sgn(x2 - x1).toFloat
     if dx == 0.0:
@@ -275,6 +293,26 @@ proc linspace*(x1, x2: float, N: int): seq[float] {.inline.} =
         result.add(x1 + dx * i.toFloat)
     result.add(x2)
 
+
+proc norm*(v1: Vector, p: int): float64 =
+    ## Calculate various norms of our Vector class
+    
+    # we have to make a case for p = 0 to avoid division by zero, may as well flesh them all out
+    case p:
+        of 0:
+            # max(v1) Infinity norm
+            result = float64(max(@v1))
+        of 1:
+            # sum(v1) Taxicab norm
+            result = float64(sum(@v1))
+        of 2:
+            # sqrt(sum([v ^ 2 for v in v1])) Euclidean norm
+            result = sqrt(sum(@(v1 ^ 2)))
+        else:
+            # pow(sum([v ^ p for v in v1]), 1.0/p) P norm
+            result = pow(sum(@(v1 ^ p)), 1.0 / float64(p))
+
+
 template timeit*(s: untyped, n = 100, msg = ""): untyped =
     var tTotal = 0.0
     for i in 1 .. n:

tests/test_optimize.nim

@@ -0,0 +1,53 @@
+import unittest, math, sequtils, arraymancer
+import numericalnim
+
+test "steepest_descent func":
+    proc df(x: float): float = 4 * x^3 - 9.0 * x^2
+    let start = 6.0
+    let gamma = 0.01
+    let precision = 0.00001
+    let max_iters = 10000
+    let correct = 2.24996
+    let value = steepest_descent(df, start, gamma, precision, max_iters)
+    check isClose(value, correct, tol = 1e-5)
+
+test "steepest_descent func starting at zero":
+    proc df(x: float): float = 4 * x^3 - 9.0 * x^2 + 4
+    let start = 0.0
+    let correct = -0.59301
+    let value = steepest_descent(df, start)
+    check isClose(value, correct, tol = 1e-5)
+
+test "conjugate_gradient func":
+    var A = toSeq([4.0, 1.0, 1.0, 3.0]).toTensor.reshape(2,2).astype(float64)
+    var x = toSeq([2.0, 1.0]).toTensor.reshape(2,1)
+    var b = toSeq([1.0,2.0]).toTensor.reshape(2,1)
+    let tol = 0.001
+    let correct = toSeq([0.090909, 0.636363]).toTensor.reshape(2,1).astype(float64)
+
+    let value = conjugate_gradient(A, b, x, tol)
+    check isClose(value, correct, tol = 1e-5)
+
+test "Newtons 1 dimension func":
+    proc f(x:float64): float64 = (1.0 / 3.0) * x ^ 3 - 2 * x ^ 2 + 3 * x
+    proc df(x:float64): float64 = x ^ 2 - 4 * x + 3
+    let x = 0.5
+    let correct = 0.0
+    let value = newtons(f, df, x, 0.000001, 1000)
+    check isClose(value, correct, tol=1e-5)
+
+test "Newtons 1 dimension func default args":
+    proc f(x:float64): float64 = (1.0 / 3.0) * x ^ 3 - 2 * x ^ 2 + 3 * x
+    proc df(x:float64): float64 = x ^ 2 - 4 * x + 3
+    let x = 0.5
+    let correct = 0.0
+    let value = newtons(f, df, x)
+    check isClose(value, correct, tol=1e-5)
+
+test "Newtons unable to find a root":
+    proc bad_f(x:float64): float64 = pow(E, x) + 1
+    proc bad_df(x:float64): float64 = pow(E, x)
+    expect(ArithmeticError):
+        discard newtons(bad_f, bad_df, 0, 0.000001, 1000)
+
+