Changeset 15449 for branches/MathNetNumerics-Exploration-2789

Timestamp:

11/03/17 20:28:37 (7 years ago)

Author:

gkronber

Message:

#2789 created x64 build of CUBGCV algorithm, fixed some bugs

Location:

branches/MathNetNumerics-Exploration-2789

Files:

• HeuristicLab.Algorithms.DataAnalysis.Experimental/642.f (modified) (1 diff)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/642_x64.dll (added)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/CubicSplineGCV.cs (modified) (2 diffs)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs (modified) (4 diffs)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj (modified) (1 diff)
• HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs (modified) (6 diffs)
• Main/Main.csproj (modified) (2 diffs)
• Main/Program.cs (modified) (1 diff)

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/642.f

@@ -130 +130 @@
 C-----------------------------------------------------------------------
 C
-       SUBROUTINE CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER) 
+      SUBROUTINE CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER) 
      . BIND(C, NAME='cubgcv')
       USE ISO_C_BINDING
+      !DEC$ ATTRIBUTES DLLEXPORT::CUBGCV
 C
 C---SPECIFICATIONS FOR ARGUMENTS---
-      INTEGER N,IC,JOB,IER
+      INTEGER(KIND=4) N,IC,JOB,IER
       DOUBLE PRECISION X(N),F(N),DF(N),Y(N),C(IC,3),SE(N),VAR,
      .                 WK(0:N+1,7)

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/CubicSplineGCV.cs

@@ -1 +1 @@
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Linq;
 using System.Runtime.InteropServices;
 using System.Text;
 using System.Threading.Tasks;
+using HeuristicLab.Problems.DataAnalysis;
 
 namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
@@ -105 +107 @@
     //                           IER = 133, JOB IS NOT 0 OR 1.
 
-    // TODO: Build x64 version of the library using mingw
-    // TODO: statically link all required libraries (libgcc_s_dw2-1.dll is still required)
-    // To build the fortran library use:
-    // > gfortran -static -Lc:/mingw/bin/libgcc_s_dw2-1.dll -m32 642.f -shared -o 642_x86.dll
-    // this means that libgcc_s_dw2-1 is statically linked, i386 code is produced and a shared library is produced
+    // To build the fortran library (x64) use:
+    // > gfortran -m64 642.f -shared -o 642_x64.dll  or
+    // > ifort /dll /Qm64 642.f /Fe642_x64.dll
+    // check dumpbin /EXPORTS 642_x64.dll
+    // and   dumpbin /IMPORTS 642_x64.dll
     [DllImport("642_x64.dll", CallingConvention = CallingConvention.Cdecl, EntryPoint = "cubgcv")]
     public static extern void cubgcv_x64(
-      double[] x, double[] f, double[] df, ref int n, double[] y, double[,] c, ref int ic, ref double var, ref int job, double[] se, double[,] wk);
-
-    // CUBGCV(X, F, DF, N, Y, C, IC, VAR, JOB, SE, WK, IER)
+      double[] x,
+      double[] f,
+      double[] df,
+      ref int n,
+      [Out] double[] y,
+      [Out] double[,] c,
+      ref int ic,
+      ref double var,
+      ref int job,
+      [Out] double[] se,
+      double[,] wk,
+      ref int ier);
+
+    // To build the fortran library (x86) use:
+    // > gfortran -static -Lc:/mingw/bin/libgcc_s_dw2-1.dll -m32 642.f -shared -o 642_x86.dll
     [DllImport("642_x86.dll", CallingConvention = CallingConvention.Cdecl, EntryPoint = "cubgcv")]
     public static extern void cubgcv_x86(
-      double[] x, 
-      double[] f, 
-      double[] df, 
-      ref int n, 
-      double[] y, 
-      double[,] c, 
-      ref int ic, 
-      ref double var, 
-      ref int job, 
-      double[] se, 
+      double[] x,
+      double[] f,
+      double[] df,
+      ref Int32 n,
+      double[] y,
+      double[,] c,
+      ref Int32 ic,
+      ref double var,
+      ref Int32 job,
+      double[] se,
       double[,] wk,
-      ref int ier);
+      ref Int32 ier);
+
+
+    public class CubGcvReport {
+      //
+      //   
+      //
+      //
+      //   IF WK(1)=0 (RHO=0) AN INTERPOLATING NATURAL CUBIC
+      //   SPLINE HAS BEEN CALCULATED.
+      //   IF WK(1)=1 (RHO=INFINITE) A LEAST SQUARES
+      //   REGRESSION LINE HAS BEEN CALCULATED.
+      //
+      public double smoothingParameter;
+      //   WK(2) IS AN ESTIMATE OF THE NUMBER OF DEGREES OF
+      //   FREEDOM OF THE RESIDUAL WHICH REDUCES TO THE
+      //   USUAL VALUE OF N-2 WHEN A LEAST SQUARES REGRESSION
+      //   LINE IS CALCULATED.
+      public double estimatedRSSDegreesOfFreedom;
+      //   WK(3),WK(4),WK(5) ARE CALCULATED WITH THE DF(I)
+      //   SCALED TO HAVE MEAN SQUARE VALUE 1.  THE
+      //   UNSCALED VALUES OF WK(3),WK(4),WK(5) MAY BE
+      //   CALCULATED BY DIVIDING BY WK(7).
+      public double generalizedCrossValidation;
+      public double meanSquareResiudal;
+      public double estimatedTrueMeanSquaredErrorAtDataPoints;
+
+      //   WK(6) COINCIDES WITH THE OUTPUT VALUE OF VAR IF
+      //   VAR IS NEGATIVE ON INPUT.  IT IS CALCULATED WITH
+      //   THE UNSCALED VALUES OF THE DF(I) TO FACILITATE
+      //   COMPARISONS WITH A PRIORI VARIANCE ESTIMATES.
+      public double estimatedErrorVariance;
+      public double meanSquareOfDf;
+      public double[] se;
+    }
+    public static ReinschSmoothingSplineModel CalculateCubicSpline(double[] x, double[] y,
+      string targetVariable, string[] inputVars,
+      out CubGcvReport report) {
+      var w = Enumerable.Repeat(1.0, x.Length).ToArray();
+      double[] x2, y2, w2;
+      Splines.SortAndBin(x, y, w, out x2, out y2, out w2);
+      // TODO: use weights correctly
+
+      int n = x2.Length;
+      double[] df = Enumerable.Repeat(1.0, n).ToArray();   // set each df if actual standard deviation of data points is not known
+      double[,] c = new double[3, n - 1];
+      double var = -99.0;
+      double[] se = new double[n]; // standard errors in points
+      double[,] wk = new double[7, (n + 2)]; // work array;
+      double[] y_smoothed = new double[n];
+      int job = 1; // calc estimates of standard errors in points
+      int ic = n - 1;
+      int ier = -99;
+      if (Environment.Is64BitProcess) {
+        CubicSplineGCV.cubgcv_x64(x2, y2, df, ref n, y_smoothed,
+          c, ref ic, ref var, ref job, se, wk, ref ier);
+      } else {
+        CubicSplineGCV.cubgcv_x86(x2, y2, df, ref n, y_smoothed,
+          c, ref ic, ref var, ref job, se, wk, ref ier);
+
+
+      }
+
+      /*
+       *                  IER    - ERROR PARAMETER. (OUTPUT)
+       *                           TERMINAL ERROR
+       *                             IER = 129, IC IS LESS THAN N-1.
+       *                             IER = 130, N IS LESS THAN 3.
+       *                             IER = 131, INPUT ABSCISSAE ARE NOT
+       *                               ORDERED SO THAT X(I).LT.X(I+1).
+       *                             IER = 132, DF(I) IS NOT POSITIVE FOR SOME I.
+       *                             IER = 133, JOB IS NOT 0 OR 1.
+       *
+       */
+
+      if (ier == 131) {
+        throw new ArgumentException("x is not ordered");
+      } else if (ier != 0) throw new ArgumentException("Error in CUBGCV " + ier);
+
+      //                           OF THE SPLINE APPROXIMATION AT T IS
+      //                           S(T)=((C(I,3)*D+C(I,2))*D+C(I,1))*D+Y(I)
+      //                           WHERE X(I).LE.T.LT.X(I+1) AND
+      //                           D = T-X(I).
+
+      //                           WK(1) = SMOOTHING PARAMETER (= RHO/(RHO + 1))
+      //                           WK(2) = ESTIMATE OF THE NUMBER OF DEGREES OF
+      //                                   FREEDOM OF THE RESIDUAL SUM OF SQUARES
+      //                           WK(3) = GENERALIZED CROSS VALIDATION
+      //                           WK(4) = MEAN SQUARE RESIDUAL
+      //                           WK(5) = ESTIMATE OF THE TRUE MEAN SQUARE ERROR
+      //                                   AT THE DATA POINTS
+      //                           WK(6) = ESTIMATE OF THE ERROR VARIANCE
+      //                           WK(7) = MEAN SQUARE VALUE OF THE DF(I)
+      report = new CubGcvReport();
+      report.smoothingParameter = wk[0, 0];
+      report.estimatedRSSDegreesOfFreedom = wk[0, 1];
+      report.generalizedCrossValidation = wk[0, 2];
+      report.meanSquareResiudal = wk[0, 3];
+      report.estimatedTrueMeanSquaredErrorAtDataPoints = wk[0, 4];
+      report.estimatedErrorVariance = wk[0, 5];
+      report.meanSquareOfDf = wk[0, 6];
+      report.se = se;
+      Debug.Assert(var == report.estimatedErrorVariance);
+
+
+
+      double[] dArr = new double[n];
+      double[] cArr = new double[n];
+      double[] bArr = new double[n];
+
+      for (int i = 0; i < n - 1; i++) {
+        bArr[i] = c[0, i];
+        cArr[i] = c[1, i];
+        dArr[i] = c[2, i];
+      }
+
+      // extrapolate linearly for xx > x[n]
+      bArr[n - 1] = bArr[n - 2];
+      // extrapolate linearly for xx < x[1]
+      dArr[0] = 0;
+
+      return new ReinschSmoothingSplineModel(
+        new MyArray<double>(1, y_smoothed),
+        new MyArray<double>(1, bArr),
+        new MyArray<double>(1, cArr),
+        new MyArray<double>(1, dArr),
+        new MyArray<double>(1, x2),
+        targetVariable, inputVars);
+
+
+    }
+
   }
-
 }

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs

@@ -107 +107 @@
       var avgY = y.Average();
       var inputVars = Problem.ProblemData.AllowedInputVariables.ToArray();
-      var nTerms = inputVars.Length; // LR
+      var nTerms = 0; // inputVars.Length; // LR
       for (int i = 1; i <= maxInteractions; i++) {
         nTerms += inputVars.Combinations(i).Count();
@@ -136 +136 @@
       while (iters++ < maxIters) {
         int j = 0;
-        foreach (var inputVar in inputVars) {
-          var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
-          rss[j] = res.Variance();
-          terms[j] = inputVar;
-          f[j] = RegressLR(problemData, inputVar, res);
-          j++;
-        }
+        //foreach (var inputVar in inputVars) {
+        //  var res = CalculateResiduals(problemData, f, j, avgY, problemData.TrainingIndices);
+        //  rss[j] = res.Variance();
+        //  terms[j] = inputVar;
+        //  f[j] = RegressLR(problemData, inputVar, res);
+        //  j++;
+        //}
 
         for (int interaction = 1; interaction <= maxInteractions; interaction++) {
@@ -202 +202 @@
     }
 
-    private IRegressionModel RegressSpline(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
-      if (problemData.Dataset.VariableHasType<double>(inputVar)) {
-        // Umständlich!
-        return Splines.CalculatePenalizedRegressionSpline(
-          problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray(),
-          (double[])target.Clone(), lambda,
-          problemData.TargetVariable, new string[] { inputVar }
-          );
-      } else return new ConstantModel(target.Average(), problemData.TargetVariable);
-    }
+    // private IRegressionModel RegressSpline(IRegressionProblemData problemData, string inputVar, double[] target, double lambda) {
+    //   if (problemData.Dataset.VariableHasType<double>(inputVar)) {
+    //     // Umständlich!
+    //     return Splines.CalculatePenalizedRegressionSpline(
+    //       problemData.Dataset.GetDoubleValues(inputVar, problemData.TrainingIndices).ToArray(),
+    //       (double[])target.Clone(), lambda,
+    //       problemData.TargetVariable, new string[] { inputVar }
+    //       );
+    //   } else return new ConstantModel(target.Average(), problemData.TargetVariable);
+    // }
     private IRegressionModel RegressSpline(IRegressionProblemData problemData, string[] inputVars, double[] target, double lambda) {
       if (inputVars.All(problemData.Dataset.VariableHasType<double>)) {
@@ -218 +218 @@
           product = product.Zip(problemData.Dataset.GetDoubleValues(inputVars[i], problemData.TrainingIndices), (pi, vi) => pi * vi).ToArray();
         }
-        double optTolerance, looRMSE;
-        return Splines.CalculateSmoothingSplineReinsch(
+        CubicSplineGCV.CubGcvReport report;
+        return CubicSplineGCV.CalculateCubicSpline(
           product,
-          (double[])target.Clone(), inputVars,
-          problemData.TargetVariable, 
-          out optTolerance, out looRMSE
+          (double[])target.Clone(), 
+          problemData.TargetVariable, inputVars, out report
           );
       } else return new ConstantModel(target.Average(), problemData.TargetVariable);

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/HeuristicLab.Algorithms.DataAnalysis.Experimental.csproj

@@ -169 +169 @@
   </ItemGroup>
   <ItemGroup>
+    <Content Include="642_x64.dll">
+      <CopyToOutputDirectory>Always</CopyToOutputDirectory>
+    </Content>
     <Content Include="642_x86.dll">
       <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>

branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs

@@ -68 +68 @@
           "Smoothing Spline (Reinsch with automatic tolerance determination using LOOCV)",
           "B-Spline Smoothing",
-          "Penalized Regression Spline (alglib)"
+          "Penalized Regression Spline (alglib)",
+          "CUBGCV"
       }.Select(s => new StringValue(s)));
 
-      Parameters.Add(new ConstrainedValueParameter<StringValue>("Type", "The type of spline (as supported by alglib)", validTypes, validTypes.First()));
+      Parameters.Add(new ConstrainedValueParameter<StringValue>("Type", "The type of spline (as supported by alglib)", validTypes, validTypes.Last()));
       Parameters.Add(new ValueParameter<DoubleValue>("Lambda", "Regularization parameter for smoothing splines (0..+inf)", new DoubleValue(100)));
       Parameters.Add(new ValueParameter<DoubleValue>("StdDev (noise)", "Known error in y values. Used only be Reinsch Smoothing Splines", new DoubleValue(0.1)));
@@ -151 +152 @@
               Results.Add(new Result("Optimal tolerance", new DoubleValue(optTol)));
               Results.Add(new Result("RMSE (LOO-CV)", new DoubleValue(looRMSE)));
+              break;
+            }
+          case "CUBGCV": {
+              CubicSplineGCV.CubGcvReport report;
+              var model = 
+                CubicSplineGCV.CalculateCubicSpline(x, y, 
+                Problem.ProblemData.TargetVariable, inputVars, out report);
+              var targetVar = Problem.ProblemData.TargetVariable;
+              var problemData = (IRegressionProblemData)Problem.ProblemData.Clone();
+              Results.Add(new Result("Solution", model.CreateRegressionSolution(problemData)));
+              Results.Add(new Result("GCV", new DoubleValue(report.generalizedCrossValidation)));
+              Results.Add(new Result("Estimated error variance", new DoubleValue(report.estimatedErrorVariance)));
+              Results.Add(new Result("Estimated RSS degrees of freedom", new DoubleValue(report.estimatedRSSDegreesOfFreedom)));
+              Results.Add(new Result("Estimated treu mean squared error at data points", new DoubleValue(report.estimatedTrueMeanSquaredErrorAtDataPoints)));
+              Results.Add(new Result("Mean square of DF", new DoubleValue(report.meanSquareOfDf)));
+              Results.Add(new Result("Mean square residual", new DoubleValue(report.meanSquareResiudal)));
+              Results.Add(new Result("Smoothing parameter (optimized for GCV)", new DoubleValue(report.smoothingParameter)));
               break;
             }
@@ -489 +507 @@
       }
 
+      // extrapolate for xx > x[n2]
+      b[b.Length] = b[b.Length - 1];
+      d[1] = 0;
+
       return new ReinschSmoothingSplineModel(a, b, c, d, x, targetVar, inputVars);
     }
@@ -550 +572 @@
     }
 
-    private static void SortAndBin(double[] x, double[] y, double[] w, out double[] x2, out double[] y2, out double[] w2, bool scaling = false) {
+    public static void SortAndBin(double[] x, double[] y, double[] w, out double[] x2, out double[] y2, out double[] w2, bool scaling = false) {
       var sortedIdx = Enumerable.Range(0, x.Length).ToArray();
       // sort by x
@@ -679 +701 @@
       this.scale = scale;
       this.offset = offset;
-
-      // extrapolate for xx > x[n2]
-      b[b.Length] = b[b.Length - 1];
-      d[1] = 0;
-      // d[b.Length] = -d[d.Length - 1];
     }
 
@@ -695 +712 @@
 
     public double GetEstimatedValue(double xx) {
-      int n = x.Length;
+      int n = a.Length;
       if (xx <= x[1]) {
         double h = xx - x[1];

branches/MathNetNumerics-Exploration-2789/Main/Main.csproj

@@ -23 +23 @@
     <ErrorReport>prompt</ErrorReport>
     <WarningLevel>4</WarningLevel>
+    <Prefer32Bit>false</Prefer32Bit>
   </PropertyGroup>
   <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
@@ -40 +41 @@
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Common-3.3.dll</HintPath>
     </Reference>
+    <Reference Include="HeuristicLab.Core-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec" />
     <Reference Include="HeuristicLab.Problems.DataAnalysis-3.4, Version=3.4.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
       <SpecificVersion>False</SpecificVersion>

branches/MathNetNumerics-Exploration-2789/Main/Program.cs

@@ -13 +13 @@
     static void Main(string[] args) {
       var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.1, includeEnd: true).ToList();
-      var ys = xs.Select(xi => 1.0 / Math.Sqrt(2 *Math.PI) * Math.Exp(-0.5 * xi*xi)).ToArray(); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+      var ys = xs.Select(xi => 1.0 / Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi)).ToArray(); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
 
+      int n = xs.Count();
       alglib.hqrndstate state;
       alglib.hqrndseed(1234, 5678, out state);
       var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.1).ToList();
 
-      int n = xs.Count;
-      double[] ys_smoothed = new double[n];
-      double[] df = Enumerable.Repeat(1.0, n).ToArray();   // set each df if actual standard deviation of data points is not known
-      double[,] c = new double[n - 1, 3];
-      double var = -99.0;
-      int ic = n - 1;
-      int job = 1; // calc estimates of standard errors in points
-      double[] se = new double[n]; // standard errors in points
-      double[,] wk = new double[7, (n + 2)]; // work array;
-      int ier = -99;
-      if (Environment.Is64BitProcess) {
-        // CubicSplineGCV.cubgcv_x64(xs.ToArray(), f, df, ref n, ys_noise.ToArray(), c, ref ic, ref var, ref job, se, wk);
-        Console.WriteLine("x64 version not supported");
-      } else {
-        CubicSplineGCV.cubgcv_x86(xs.ToArray(), ys_noise.ToArray(), df, ref n, ys_smoothed,
-          c, ref ic, ref var, ref job, se, wk, ref ier);
+      CubicSplineGCV.CubGcvReport report;
+      var model = CubicSplineGCV.CalculateCubicSpline(
+        xs.ToArray(), ys_noise.ToArray(), "y", new string[] { "x" }, out report);
 
+      Console.WriteLine("Smoothing Parameter (= RHO/(RHO + 1) {0}", report.smoothingParameter);
+      Console.WriteLine("Estimated DOF of RSS {0}", report.estimatedRSSDegreesOfFreedom);
+      Console.WriteLine("GCV {0}", report.generalizedCrossValidation);
+      Console.WriteLine("Mean squared residual {0}", report.meanSquareResiudal);
+      Console.WriteLine("Estimate of true MSE at data points {0}", report.estimatedTrueMeanSquaredErrorAtDataPoints);
+      Console.WriteLine("Estimate of error variance {0}", report.estimatedErrorVariance);
+      Console.WriteLine("Mean square value of DF(I) {0}", report.meanSquareOfDf);
 
-        //                           WK(1) = SMOOTHING PARAMETER (= RHO/(RHO + 1))
-        //                           WK(2) = ESTIMATE OF THE NUMBER OF DEGREES OF
-        //                                   FREEDOM OF THE RESIDUAL SUM OF SQUARES
-        //                           WK(3) = GENERALIZED CROSS VALIDATION
-        //                           WK(4) = MEAN SQUARE RESIDUAL
-        //                           WK(5) = ESTIMATE OF THE TRUE MEAN SQUARE ERROR
-        //                                   AT THE DATA POINTS
-        //                           WK(6) = ESTIMATE OF THE ERROR VARIANCE
-        //                           WK(7) = MEAN SQUARE VALUE OF THE DF(I)
+      OnlineCalculatorError error;
+      var ys_smoothed = xs.Select(xi => model.GetEstimatedValue(xi)).ToArray();
+      var mse = OnlineMeanSquaredErrorCalculator.Calculate(ys, ys_smoothed, out error);
+      Console.WriteLine("MSE(ys, ys_smooth) = {0}", mse);
+      mse = OnlineMeanSquaredErrorCalculator.Calculate(ys, ys_noise, out error);
+      Console.WriteLine("MSE(ys, ys_noise) = {0}", mse);
 
-        Console.WriteLine("Smoothing Parameter (= RHO/(RHO + 1) {0}", wk[0, 0]);
-        Console.WriteLine("Estimated DOF of RSS {0}", wk[0, 1]);
-        Console.WriteLine("GCV {0}", wk[0, 2]);
-        Console.WriteLine("Mean squared residual {0}", wk[0, 3]);
-        Console.WriteLine("Estimate of true MSE at data points {0}", wk[0, 4]);
-        Console.WriteLine("Estimate of error variance {0}", wk[0, 5]);
-        Console.WriteLine("Mean square value of DF(I) {0}", wk[0, 6]);
+      Thread.CurrentThread.CurrentCulture = CultureInfo.InvariantCulture;
 
-        OnlineCalculatorError error;
-        var mse = OnlineMeanSquaredErrorCalculator.Calculate(ys, ys_smoothed, out error);
-        Console.WriteLine("MSE(ys, ys_smooth) = {0}", mse);
-        mse = OnlineMeanSquaredErrorCalculator.Calculate(ys, ys_noise, out error);
-        Console.WriteLine("MSE(ys, ys_noise) = {0}", mse);
-
-        Thread.CurrentThread.CurrentCulture = CultureInfo.InvariantCulture;
-
-        for (int i=0;i<n;i++) {
-          Console.WriteLine("{0}\t{1}\t{2}\t{3}\t{4}\t{5}", 
-            xs[i], ys[i], ys_smoothed[i], ys_noise[i], ys_smoothed[i] + 1.96 * se[i], ys_smoothed[i] - 1.96 * se[i]);
-        }
+      for (int i = 0; i < n; i++) {
+        Console.WriteLine("{0}\t{1}\t{2}\t{3}\t{4}\t{5}",
+          xs[i], ys[i], ys_smoothed[i], ys_noise[i], ys_smoothed[i] + 1.96 * report.se[i], ys_smoothed[i] - 1.96 * report.se[i]);
       }
     }