Ignore:
Timestamp:
11/09/17 18:03:06 (3 years ago)
Author:
gkronber
Message:

#2789 worked on sbart spline

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs

    r15465 r15468  
    148148      var validTypes = new ItemSet<StringValue>(
    149149        new[] {
    150         "Monotone", "Akima", "Catmull-Rom", "Cubic", "Linear",
    151           "Cubic - Natural (Math.NET)",
    152           "Polynomial (Math.NET)",
    153           "Rational (Math.NET)",
    154           "LogLinear (Math.NET)",
    155           "Common (Math.NET)",
    156           "Smoothing Spline (Mine)",
    157150          "Smoothing Spline (Reinsch)",
    158151          "Smoothing Spline (Reinsch with automatic tolerance determination using LOOCV)",
     
    185178        var type = ((StringValue)Parameters["Type"].ActualValue).Value;
    186179        switch (type) {
    187           case "Monotone":
    188             alglib.spline1dbuildmonotone(x, y, out c);
    189             AddAlglibSplineResult(c, inputVars);
    190             break;
    191           case "Akima":
    192             alglib.spline1dbuildakima(x, y, out c); AddAlglibSplineResult(c, inputVars);
    193             break;
    194             ;
    195           case "Catmull-Rom":
    196             alglib.spline1dbuildcatmullrom(x, y, out c); AddAlglibSplineResult(c, inputVars);
    197             break;
    198 
    199           case "Cubic":
    200             alglib.spline1dbuildcubic(x, y, out c); AddAlglibSplineResult(c, inputVars);
    201             break;
    202 
    203           case "Linear":
    204             alglib.spline1dbuildlinear(x, y, out c); AddAlglibSplineResult(c, inputVars);
    205             break;
    206           case "Cubic - Natural (Math.NET)": {
    207               var spline = MathNet.Numerics.Interpolation.CubicSpline.InterpolateNatural(x, y);
    208               AddMathNetSplineResult(spline, inputVars);
    209               break;
    210             }
    211           case "Common (Math.NET)": {
    212               var spline = MathNet.Numerics.Interpolate.Common(x, y);
    213               AddMathNetSplineResult(spline, inputVars);
    214               break;
    215             }
    216           case "LogLinear (Math.NET)": {
    217               var spline = MathNet.Numerics.Interpolate.LogLinear(x, y);
    218               AddMathNetSplineResult(spline, inputVars);
    219               break;
    220             }
    221           case "Polynomial (Math.NET)": {
    222               var spline = MathNet.Numerics.Interpolate.Polynomial(x, y);
    223               AddMathNetSplineResult(spline, inputVars);
    224               break;
    225             }
    226           case "Rational (Math.NET)": {
    227               var spline = MathNet.Numerics.Interpolate.RationalWithoutPoles(x, y);
    228               AddMathNetSplineResult(spline, inputVars);
    229               break;
    230             }
    231           case "Smoothing Spline (Mine)": {
    232               CalculateSmoothingSpline(x, y, inputVars);
    233               break;
    234             }
    235180          case "Smoothing Spline (Reinsch)": {
    236181              CalculateSmoothingSplineReinsch(x, y, inputVars);
     
    264209          case "Finnbar O'Sullivan (sbart)": {
    265210              SBART.SBART_Report report;
    266               var lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value; // controls extent of smoothing
     211           
    267212              var model =
    268213                SBART.CalculateSBART(x, y,
    269                 Problem.ProblemData.TargetVariable, inputVars, (float)lambda, out report);
     214                Problem.ProblemData.TargetVariable, inputVars, out report);
    270215              var targetVar = Problem.ProblemData.TargetVariable;
    271216              var problemData = (IRegressionProblemData)Problem.ProblemData.Clone();
     
    275220              break;
    276221
    277             }
    278           case "B-Spline Smoothing": {
    279               BSplineSmoothing(x, y, inputVars);
    280               break;
    281222            }
    282223          case "Penalized Regression Spline (alglib)": {
     
    371312
    372313      return new RegressionEnsembleModel(models);
    373     }
    374 
    375     public void BSplineSmoothing(double[] x, double[] y, string[] inputVars) {
    376       Array.Sort(x, y);
    377       CubicSpline2(x, y, inputVars);
    378314    }
    379315
     
    414350    }
    415351
    416     public void CubicSpline2(double[] x, double[] y, string[] inputVars) {
    417       // see Pollock: Smoothing Splines
    418       int n = x.Length;
    419       double lambda = 1.0;
    420       double[] sigma = Enumerable.Repeat(1.0, n).ToArray();
    421       SplineData[] S = new SplineData[x.Length + 1];
    422       for (int i = 0; i < x.Length; i++) {
    423         S[i].x = x[i];
    424         S[i].y = y[i];
    425       }
    426       S[x.Length].x = S[x.Length - 1].x;
    427       S[x.Length].y = S[x.Length - 1].y;
    428 
    429       // var
    430       double[] h = new double[n],
    431         r = new double[n],
    432         f = new double[n],
    433         p = new double[n];
    434       var q = new MyArray<double>(-1, n + 3);
    435       var u = new MyArray<double>(-1, n + 3);
    436       var v = new MyArray<double>(-1, n + 3);
    437       var w = new MyArray<double>(-1, n + 3);
    438       double mu;
    439 
    440       mu = 2 * (1 - lambda) / (3 * lambda);
    441       h[0] = S[1].x - S[0].x;
    442       r[0] = 3 / h[0];
    443       for (int i = 1; i < n - 1; i++) {
    444         h[i] = S[i + 1].x - S[i].x;
    445         r[i] = 3 / h[i];
    446         f[i] = -(r[i - 1] + r[i]);
    447         p[i] = 2 * (S[i + 1].x - S[i - 1].x);
    448         q[i] = 3 * (S[i + 1].y - S[i].y) / h[i] - 3 * (S[i].y - S[i - 1].y) / h[i - 1];
    449       }
    450       for (int i = 1; i < n; i++) {
    451         u[i] = Math.Pow(r[i - 1], 2) * sigma[i - 1] + Math.Pow(f[i], 2) * sigma[i] + Math.Pow(r[i], 2) * sigma[i + 1];    // TODO Sigma 1..n
    452         u[i] = mu * u[i] + p[i];
    453         v[i] = f[i] * r[i] * sigma[i] + r[i] * f[i + 1] * sigma[i + 1];
    454         v[i] = mu * v[i] + h[i];
    455         w[i] = mu * r[i] * r[i + 1] * sigma[i + 1];
    456       }
    457       Quincunx(n, u, v, w, q);
    458       // { Spline P arameters}
    459       S[0].d = S[0].y - mu * r[0] * q[1] * sigma[0];
    460       S[1].d = S[1].y - mu * (f[1] * q[1] + r[1] * q[2]) * sigma[0];
    461       S[0].a = q[1] / (3 * h[0]);
    462       S[0].b = 0;
    463       S[0].c = (S[1].d - S[0].d) / h[0] - q[1] * h[0] / 3;
    464       r[0] = 0;
    465       for (int j = 1; j < n; j++) {
    466         S[j].a = (q[j + 1] - q[j]) / (3 * h[j]);
    467         S[j].b = q[j];
    468         S[j].c = (q[j] + q[j - 1]) * h[j - 1] + S[j - 1].c;
    469         S[j].d = r[j - 1] * q[j - 1] + f[j] * q[j] + r[j] * q[j + 1];
    470         S[j].d = y[j] - mu * S[j].d * sigma[j];
    471       }
    472 
    473       //   { SmoothingSpline}
    474     }
    475 
    476352    public void CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars) {
    477353      double s = x.Length + 1;
     
    541417
    542418
    543     // automatically determines tolerance using loo cv, SLOW!
     419    // automatically determines tolerance using cv, SLOW!
    544420    public static IRegressionModel CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars, string targetVar,
    545421      out double optTolerance, out double cvRMSE) {
     
    665541
    666542      return new ReinschSmoothingSplineModel(a, b, c, d, x, targetVar, inputVars);
    667     }
    668 
    669     private void CalculateSmoothingSpline(double[] x, double[] y, string[] inputVars) {
    670       // see Smoothing and Non-Parametric Regression, Germán Rodríguez, 2001 2.3.1
    671       double[] w = Enumerable.Repeat(1.0, x.Length).ToArray();  // weights necessary for sortAndBin but are ignored below (TODO)
    672       SortAndBin(x, y, w, out x, out y, out w, scaling: false);
    673       int n = x.Length;
    674 
    675       SparseMatrix delta = new SparseMatrix(n - 2, n);
    676       // double[,] delta = new double[n - 2, n];
    677       //double[,] W = new double[n - 2, n - 2];
    678       SparseMatrix W = new SparseMatrix(n - 2, n - 2);
    679       Matrix WInvD = new DenseMatrix(n - 2, n);
    680 
    681       // double[,] W_inv_D = new double[n - 2, n];
    682       // double[,] K = new double[n, n];
    683 
    684       // go over successive knots to determine distances and fill Delta and W
    685       for (int i = 0; i < n - 2; i++) {
    686         double h = x[i + 1] - x[i];
    687         double h_next = x[i + 2] - x[i + 1];
    688         delta[i, i] = 1.0 / h;
    689         delta[i, i + 1] = -1.0 / h - 1.0 / h_next;
    690         delta[i, i + 2] = 1.0 / h_next;
    691         W[i, i] = (h + h_next) / 3;
    692         if (i > 0) {
    693           W[i - 1, i] =
    694            W[i, i - 1] = h / 6.0;
    695         }
    696       }
    697 
    698       // WInvD = W.Cholesky().Solve(delta);
    699       var solvResult = W.TrySolveIterative(delta, WInvD, new MlkBiCgStab());
    700 
    701       // alglib.ablas.rmatrixgemm(n - 2, n, n - 2, 1.0, W, 0, 0, 0, delta, 0, 0, 0, 1.0, W_inv_D, 0, 0); // W^-1(M = n-2, K = n-2) D(K = n-2, N=n)   
    702       // alglib.ablas.rmatrixgemm(n, n, n - 2, 1.0, delta, 0, 0, 1, W_inv_D, 0, 0, 0, 1.0, K, 0, 0); // D(M=n-2, K=n)^T * W^-1D (K=n, N=n-2)
    703 
    704       var K = delta.TransposeThisAndMultiply(WInvD);
    705 
    706       double lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value;
    707 
    708       for (int i = 0; i < n; i++) {
    709         for (int j = 0; j < n; j++) {
    710           K[i, j] *= lambda;
    711           if (i == j) K[i, j] += 1;
    712         }
    713       }
    714 
    715       // solve for y
    716       // double[] s;
    717       // int solverInfo;
    718       // alglib.densesolverreport solverRep;
    719       // alglib.rmatrixsolve(K, n, y, out solverInfo, out solverRep, out s);
    720 
    721       var s = K.Solve(new DenseVector(y)).ToArray();
    722 
    723       Results.Add(new Result("Solution", new RegressionSolution(new SmoothingSplineModel(s, x, Problem.ProblemData.TargetVariable, inputVars),
    724   (IRegressionProblemData)Problem.ProblemData.Clone())));
    725543    }
    726544
     
    766584      w2 = wl.ToArray();
    767585    }
    768 
    769     private void AddAlglibSplineResult(alglib.spline1dinterpolant c, string[] inputVars) {
    770       Results.Add(new Result("Solution", new RegressionSolution(new Spline1dModel(c, Problem.ProblemData.TargetVariable, inputVars),
    771         (IRegressionProblemData)Problem.ProblemData.Clone())));
    772 
    773     }
    774     private void AddMathNetSplineResult(IInterpolation c, string[] inputVars) {
    775       Results.Add(new Result("Solution", new RegressionSolution(new MathNetSplineModel(c, Problem.ProblemData.TargetVariable, inputVars),
    776         (IRegressionProblemData)Problem.ProblemData.Clone())));
    777     }
    778586  }
    779587
     
    905713
    906714  // UNFINISHED
    907   public class SmoothingSplineModel : NamedItem, IRegressionModel {
    908     private double[] s;
    909     private IInterpolation interpolation;
    910 
    911     public string TargetVariable { get; set; }
    912 
    913     public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
    914 
    915     public event EventHandler TargetVariableChanged;
    916 
    917     public SmoothingSplineModel(SmoothingSplineModel orig, Cloner cloner) : base(orig, cloner) {
    918       this.TargetVariable = orig.TargetVariable;
    919       this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
    920       this.s = orig.s; // TODO
    921       this.interpolation = orig.interpolation;
    922     }
    923     public SmoothingSplineModel(double[] s, double[] x, string targetVar, string[] inputs) : base("SplineModel", "SplineModel") {
    924       this.s = s;
    925       this.TargetVariable = targetVar;
    926       this.VariablesUsedForPrediction = inputs;
    927       this.interpolation = MathNet.Numerics.Interpolate.CubicSpline(x, s);
    928     }
    929 
    930     public override IDeepCloneable Clone(Cloner cloner) {
    931       return new SmoothingSplineModel(this, cloner);
    932     }
    933 
    934     public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
    935       return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
    936     }
    937 
    938     public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
    939       foreach (var x in dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows)) {
    940 
    941         yield return interpolation.Interpolate(x);
    942 
    943       }
    944     }
    945   }
    946 
    947   // UNFINISHED
    948715  public class Spline1dModel : NamedItem, IRegressionModel {
    949     private alglib.spline1dinterpolant interpolant;
     716    alglib.spline1dinterpolant interpolant;
    950717
    951718    public string TargetVariable { get; set; }
     
    958725      this.TargetVariable = orig.TargetVariable;
    959726      this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
     727
    960728      this.interpolant = (alglib.spline1dinterpolant)orig.interpolant.make_copy();
    961729    }
    962     public Spline1dModel(alglib.spline1dinterpolant interpolant, string targetVar, string[] inputs) : base("SplineModel", "SplineModel") {
    963       this.interpolant = interpolant;
     730    internal Spline1dModel(
     731      alglib.spline1dinterpolant interpolant,
     732      string targetVar, string[] inputs, double offset = 0, double scale = 1) : base("SplineModel", "SplineModel") {
     733      this.interpolant = (alglib.spline1dinterpolant)interpolant.make_copy();
    964734      this.TargetVariable = targetVar;
    965735      this.VariablesUsedForPrediction = inputs;
     
    972742    public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
    973743      return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
     744    }
     745
     746    public double GetEstimatedValue(double xx) {
     747      return alglib.spline1dcalc(interpolant, xx);
    974748    }
    975749
     
    979753        product = product.Zip(dataset.GetDoubleValues(factor, rows), (pi, fi) => pi * fi).ToArray();
    980754      }
    981 
    982       foreach (var x in product) {
    983         yield return alglib.spline1dcalc(interpolant, x);
    984       }
    985     }
    986   }
    987 
    988 
    989   // UNFINISHED
    990   public class MathNetSplineModel : NamedItem, IRegressionModel {
    991     private IInterpolation interpolant;
    992 
    993     public string TargetVariable { get; set; }
    994 
    995     public IEnumerable<string> VariablesUsedForPrediction { get; private set; }
    996 
    997     public event EventHandler TargetVariableChanged;
    998 
    999     public MathNetSplineModel(MathNetSplineModel orig, Cloner cloner) : base(orig, cloner) {
    1000       this.TargetVariable = orig.TargetVariable;
    1001       this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray();
    1002       this.interpolant = orig.interpolant; // TODO COPY!
    1003     }
    1004     public MathNetSplineModel(IInterpolation interpolant, string targetVar, string[] inputs) : base("SplineModel", "SplineModel") {
    1005       this.interpolant = interpolant;
    1006       this.TargetVariable = targetVar;
    1007       this.VariablesUsedForPrediction = inputs;
    1008     }
    1009 
    1010     public override IDeepCloneable Clone(Cloner cloner) {
    1011       return new MathNetSplineModel(this, cloner);
    1012     }
    1013 
    1014     public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
    1015       return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone());
    1016     }
    1017 
    1018     public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
    1019       foreach (var x in dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows)) {
    1020         yield return interpolant.Interpolate(x);
     755      foreach (var xx in product) {
     756        yield return GetEstimatedValue(xx);
    1021757      }
    1022758    }
Note: See TracChangeset for help on using the changeset viewer.