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

Timestamp:

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

Author:

gkronber

Message:

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

File:

• branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/CubicSplineGCV.cs (modified) (2 diffs)

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);
+
+
+    }
+
   }
-
 }