From 8fc3ba94701964c45d50930f80a74679e9e3b404 Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Wed, 26 Jun 2013 14:25:44 +0200
Subject: [PATCH] Fitting: Idiomatic least squares curve fitting for F#

---
 src/FSharp/FSharp.fsproj                 | 13 +++---
 src/FSharp/LeastSquares.fs               | 55 ++++++++++++++++++++++++
 src/FSharpPortable/FSharpPortable.fsproj |  3 ++
 src/FSharpUnitTests/CurveFittingTests.fs |  9 ++--
 4 files changed, 71 insertions(+), 9 deletions(-)
 create mode 100644 src/FSharp/LeastSquares.fs

diff --git a/src/FSharp/FSharp.fsproj b/src/FSharp/FSharp.fsproj
index 6465831c..24b296ff 100644
--- a/src/FSharp/FSharp.fsproj
+++ b/src/FSharp/FSharp.fsproj
@@ -1,4 +1,4 @@
-<?xml version="1.0" encoding="utf-8"?>
+<?xml version="1.0" encoding="utf-8"?>
 <Project ToolsVersion="4.0" DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
   <PropertyGroup>
     <Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
@@ -14,7 +14,7 @@
     <TargetFrameworkProfile />
     <!-- Conditional Strong Name -->
     <AssemblyOriginatorKeyFile>..\MathNet.Numerics.snk</AssemblyOriginatorKeyFile>
-    <SignAssembly Condition=" '$(SignAssembly)' == '' ">false</SignAssembly>    
+    <SignAssembly Condition=" '$(SignAssembly)' == '' ">false</SignAssembly>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
     <DebugSymbols>true</DebugSymbols>
@@ -34,12 +34,12 @@
     <WarningLevel>3</WarningLevel>
     <DocumentationFile>..\..\out\lib\Net40\MathNet.Numerics.FSharp.xml</DocumentationFile>
     <!-- Conditional Strong Name: NO -->
-    <SignAssembly>false</SignAssembly>    
+    <SignAssembly>false</SignAssembly>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)' == 'Release' ">
     <DocumentationFile>MathNet.Numerics.FSharp.XML</DocumentationFile>
     <!-- Conditional Strong Name: NO -->
-    <SignAssembly>false</SignAssembly>    
+    <SignAssembly>false</SignAssembly>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release-Signed|AnyCPU' ">
     <DebugType>pdbonly</DebugType>
@@ -50,7 +50,7 @@
     <DocumentationFile>..\..\out\lib\Net40\MathNet.Numerics.FSharp.xml</DocumentationFile>
     <OutputPath>..\..\out\lib\Net40\</OutputPath>
     <!-- Conditional Strong Name: YES -->
-    <SignAssembly>true</SignAssembly>    
+    <SignAssembly>true</SignAssembly>
   </PropertyGroup>
   <PropertyGroup>
     <MinimumVisualStudioVersion Condition="'$(MinimumVisualStudioVersion)' == ''">11</MinimumVisualStudioVersion>
@@ -71,6 +71,7 @@
     <Compile Include="BigIntegerExtensions.fs" />
     <Compile Include="BigRational.fsi" />
     <Compile Include="BigRational.fs" />
+    <Compile Include="LeastSquares.fs" />
     <Compile Include="RandomVariable.fs" />
   </ItemGroup>
   <ItemGroup>
@@ -94,4 +95,4 @@
 	<Target Name="AfterBuild">
 	</Target>
 	-->
-</Project>
+</Project>
\ No newline at end of file
diff --git a/src/FSharp/LeastSquares.fs b/src/FSharp/LeastSquares.fs
new file mode 100644
index 00000000..71f92609
--- /dev/null
+++ b/src/FSharp/LeastSquares.fs
@@ -0,0 +1,55 @@
+// <copyright file="LeastSquares.fs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics
+
+open System
+open MathNet.Numerics.LinearAlgebra.Double
+open MathNet.Numerics.LinearAlgebra.Generic.Factorization
+
+[<CompilationRepresentation(CompilationRepresentationFlags.ModuleSuffix)>]
+module Fit =
+
+    let private tofs (f:Func<_,_>) = fun a -> f.Invoke(a)
+
+    let line x y = let p = LeastSquares.FitToLine(x,y) in (p.[0],p.[1])
+    let linef x y = LeastSquares.FitToLineFunc(x,y) |> tofs
+
+    let polynomial order x y = LeastSquares.FitToPolynomial(x,y,order)
+    let polynomialf order x y = LeastSquares.FitToPolynomialFunc(x,y,order) |> tofs
+
+    let linear functions (x:float[]) (y:float[]) =
+        functions
+        |> List.map (fun f -> List.init (Array.length x) (fun i -> f x.[i]))
+        |> DenseMatrix.ofColumnsList (Array.length x) (List.length functions)
+        |> fun m -> m.QR(QRMethod.Thin).Solve(DenseVector(y)).ToArray()
+    let linearf functions x y =
+        let parameters = linear functions x y |> List.ofArray
+        in fun z -> functions |> List.zip parameters |> List.fold (fun s (p,f) -> s+p*(f z)) 0.0
diff --git a/src/FSharpPortable/FSharpPortable.fsproj b/src/FSharpPortable/FSharpPortable.fsproj
index 6ec5ad09..af6ac149 100644
--- a/src/FSharpPortable/FSharpPortable.fsproj
+++ b/src/FSharpPortable/FSharpPortable.fsproj
@@ -77,6 +77,9 @@
     <Compile Include="..\FSharp\BigRational.fs">
       <Link>BigRational.fs</Link>
     </Compile>
+    <Compile Include="..\FSharp\LeastSquares.fs">
+      <Link>LeastSquares.fs</Link>
+    </Compile>
     <Compile Include="..\FSharp\RandomVariable.fs">
       <Link>RandomVariable.fs</Link>
     </Compile>
diff --git a/src/FSharpUnitTests/CurveFittingTests.fs b/src/FSharpUnitTests/CurveFittingTests.fs
index 3e5c2a6e..44aec84a 100644
--- a/src/FSharpUnitTests/CurveFittingTests.fs
+++ b/src/FSharpUnitTests/CurveFittingTests.fs
@@ -15,8 +15,10 @@ module CurveFittingTests =
         let x = Array.append [| 1.0 .. 2.0 .. 10.0 |] [| -1.0 .. -1.0 .. -5.0  |]
         let y = x |> Array.map f
 
-        LeastSquares.FitToLine(x,y)
-        |> should (equalWithin 1.0e-12) [| 4.0; -1.5 |]
+        // LeastSquares.FitToLine(x,y)
+        let a, b = Fit.line x y
+        a |> should (equalWithin 1.0e-12) 4.0
+        b |> should (equalWithin 1.0e-12) -1.5
 
         let fres = LeastSquares.FitToLineFunc(x,y) |> tofs
         in x |> Array.iter (fun x -> fres x |> should (equalWithin 1.0e-12) (f x))
@@ -30,7 +32,8 @@ module CurveFittingTests =
         let x = [| 1.0 .. 6.0 |]
         let y = [| 4.986; 2.347; 2.061; -2.995; -2.352; -5.782 |]
 
-        LeastSquares.FitToLinearCombination(x, y, (fun z -> 1.0), (fun z -> Math.Sin(z)), (fun z -> Math.Cos(z)))
+        // LeastSquares.FitToLinearCombination(x, y, (fun z -> 1.0), (fun z -> Math.Sin(z)), (fun z -> Math.Cos(z)))
+        (x,y) ||> Fit.linear [(fun _ -> 1.0); (Math.Sin); (Math.Cos) ]
         |> should (equalWithin 1.0e-4) [| -0.287476; 4.02159; -1.46962 |]
 
         let fres = LeastSquares.FitToLinearCombinationFunc(x, y, (fun z -> 1.0), (fun z -> Math.Sin(z)), (fun z -> Math.Cos(z))) |> tofs