Add SpecialFunctions folder

FSharp.Stats.sln

@@ -30,6 +30,7 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "content", "content", "{8E6D
 	ProjectSection(SolutionItems) = preProject
 		docs\content\index.fsx = docs\content\index.fsx
 		docs\content\Matrix_Vector.fsx = docs\content\Matrix_Vector.fsx
+		docs\content\Optimization.fsx = docs\content\Optimization.fsx
 		docs\content\tutorial.fsx = docs\content\tutorial.fsx
 	EndProjectSection
 EndProject

docs/content/Matrix_Vector.fsx

@@ -37,6 +37,14 @@ let rv =
     rowvec [|2.0; 20.0; 1.|]
 
 
+let f x = x**2.
+
+"2.;3.;".Split(';') 
+|> Array.filter (fun x -> x <> "") 
+|> Array.map float 
+|> Array.map f
+
+
 
 A+B
 A-B
@@ -61,3 +69,65 @@ Matrix.dot A B
 // 
 // http://blog.ploeh.dk/2011/04/27/Providerisnotapattern/
 // http://blog.ploeh.dk/2011/04/27/Providerisnotapattern/
+
+
+open System.Runtime.InteropServices
+
+// [<DllImport("C:/Program Files (x86)/Sho 2.1/bin/bin64/mkl.dll",EntryPoint="dgemm_")>]
+[<DllImport("D:\\libopenblas.dll",EntryPoint="dgemm_")>]
+extern void dgemm_
+  ( char *transa, char *transb, 
+    int *m, int *n, int *k, 
+    double *alpha, double *a, int *lda, 
+    double *b, int *ldb, double *beta, 
+    double *c, int *ldc );
+
+
+
+
+#nowarn "51" 
+
+open Microsoft.FSharp.NativeInterop
+
+
+let matmul_blas (a:float[,]) (b:float[,]) = 
+    // Get dimensions of the input matrices
+    let m = Array2D.length1 a
+    let k = Array2D.length2 a
+    let n = Array2D.length2 b
+
+    // Allocate array for the result
+    let c = Array2D.create n m 0.0
+
+    // Declare arguments for the call
+    let mutable arg_transa = 't'
+    let mutable arg_transb = 't'
+    let mutable arg_m = m
+    let mutable arg_n = n
+    let mutable arg_k = k
+    let mutable arg_alpha = 1.0
+    let mutable arg_ldk = k
+    let mutable arg_ldn = n
+    let mutable arg_beta = 1.0
+    let mutable arg_ldm = m
+
+    // Temporarily pin the arrays
+    use arg_a = PinnedArray2.of_array2D(a)
+    use arg_b = PinnedArray2.of_array2D(b)
+    use arg_c = PinnedArray2.of_array2D(c)
+
+    // Invoke the native routine
+    dgemm_( &&arg_transa, &&arg_transb,
+            &&arg_m, &&arg_n, &&arg_k,
+            &&arg_alpha, arg_a.Ptr, &&arg_ldk,
+            arg_b.Ptr, &&arg_ldn, &&arg_beta,
+            arg_c.Ptr, &&arg_ldm )
+
+    // Transpose the result to get m*n matrix 
+    Array2D.init m n (fun i j -> c.[j,i])
+
+
+
+
+matmul_blas (A.ToArray2D()) (B.ToArray2D())
+

docs/content/Optimization.fsx

@@ -0,0 +1,102 @@
+(*** hide ***)
+// This block of code is omitted in the generated HTML documentation. Use 
+// it to define helpers that you do not want to show in the documentation.
+//#I "../../bin"
+#r "../../packages/build/FSharp.Plotly/lib/net40/Fsharp.Plotly.dll"
+open FSharp.Plotly
+(**
+dfdf
+*)
+#r "D:/Source/FSharp.Stats/bin/FSharp.Stats.dll"
+//#r "FSharp.Stats.dll"
+//open FSharp
+open Microsoft.FSharp.Math
+
+
+open FSharp.Stats.Optimization
+
+
+// http://fsharpnews.blogspot.de/2011/01/gradient-descent.html
+
+let rosenbrock (xs: vector) =
+    let x, y = xs.[0], xs.[1]
+    pown (1.0 - x) 2 + 100.0 * pown (y - pown x 2) 2
+
+
+// The minimum at (1, 1) may be found quickly and easily using the functions defined above as follows:
+
+let xs =
+    vector[0.0; 0.0]
+    |> GradientDescent.minimize rosenbrock (GradientDescent.grad rosenbrock)
+
+
+
+
+
+
+
+
+let x = vector [|-2. ..0.1.. 2.|]
+let y = vector [|-2. ..0.1.. 2.|]
+
+let rosen (x,y) =
+    pown (1.0 - x) 2 + 100.0 * pown (y - pown x 2) 2
+
+
+let z = 
+    Array.init y.Length (fun i -> 
+        Array.init x.Length (fun j -> rosen (x.[j], y.[i]) )
+                    )
+
+
+
+(*** define-output:rosenContour ***)
+Chart.Surface(z,x,y)
+|> Chart.withSize(600.,600.)
+(*** include-it:rosenContour ***)
+|> Chart.Show
+
+
+
+
+// we define a small number that we be used to calculate numerical approximations to derivatives:
+let eps = System.Double.Epsilon ** (1.0 / 3.0)
+
+// The following function repeatedly applies the given function to the given initial value until the result stops changing:
+let rec fixedPoint f x =
+    let f_x = f x
+    if f_x = x then x else fixedPoint f f_x
+
+
+// The numerical approximation to the grad of a scalar field is built up from partial derivatives in each direction:
+let partialD f_xs f (xs : vector) i xi =
+    xs.[i] <- xi + eps
+    try (f xs - f_xs) / eps finally
+    xs.[i] <- xi
+
+
+// The following function performs a single iteration of gradient descent by scaling the step size lambda by either 'a' or 'b' if the result increases or decreases the function being minimized, respectively:
+let descend a b f (f': _ -> vector) (lambda, xs, f_xs) =
+    let xs_2 = xs - lambda * f' xs
+    let f_xs_2 = f xs_2
+    if f_xs_2 >= f_xs then
+        a * lambda, xs, f_xs
+    else
+        b * lambda, xs_2, f_xs_2
+
+
+/// radient descent algorithm to minimize a given function and derivative
+let minimize f f' xs =
+    //let _, xs, _ = fixedPoint (descend 0.5 1.1 f f') (eps, xs, f xs)
+    //xs
+    fixedPoint (descend 0.5 1.1 f f') (eps, xs, f xs)
+
+/// Computes a numerical approximation to the derivative of a function
+let grad f xs =
+    Vector.mapi (partialD (f xs) f xs) xs
+
+
+let xs' =
+    vector[0.0; 0.0]
+    |> minimize rosenbrock (grad rosenbrock)
+