Changeset 15468 for branches/MathNetNumericsExploration2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs
 Timestamp:
 11/09/17 18:03:06 (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/MathNetNumericsExploration2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs
r15465 r15468 148 148 var validTypes = new ItemSet<StringValue>( 149 149 new[] { 150 "Monotone", "Akima", "CatmullRom", "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)",157 150 "Smoothing Spline (Reinsch)", 158 151 "Smoothing Spline (Reinsch with automatic tolerance determination using LOOCV)", … … 185 178 var type = ((StringValue)Parameters["Type"].ActualValue).Value; 186 179 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 "CatmullRom":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 }235 180 case "Smoothing Spline (Reinsch)": { 236 181 CalculateSmoothingSplineReinsch(x, y, inputVars); … … 264 209 case "Finnbar O'Sullivan (sbart)": { 265 210 SBART.SBART_Report report; 266 var lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value; // controls extent of smoothing211 267 212 var model = 268 213 SBART.CalculateSBART(x, y, 269 Problem.ProblemData.TargetVariable, inputVars, (float)lambda,out report);214 Problem.ProblemData.TargetVariable, inputVars, out report); 270 215 var targetVar = Problem.ProblemData.TargetVariable; 271 216 var problemData = (IRegressionProblemData)Problem.ProblemData.Clone(); … … 275 220 break; 276 221 277 }278 case "BSpline Smoothing": {279 BSplineSmoothing(x, y, inputVars);280 break;281 222 } 282 223 case "Penalized Regression Spline (alglib)": { … … 371 312 372 313 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);378 314 } 379 315 … … 414 350 } 415 351 416 public void CubicSpline2(double[] x, double[] y, string[] inputVars) {417 // see Pollock: Smoothing Splines418 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 // var430 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..n452 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 476 352 public void CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars) { 477 353 double s = x.Length + 1; … … 541 417 542 418 543 // automatically determines tolerance using loocv, SLOW!419 // automatically determines tolerance using cv, SLOW! 544 420 public static IRegressionModel CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars, string targetVar, 545 421 out double optTolerance, out double cvRMSE) { … … 665 541 666 542 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 NonParametric Regression, Germán Rodríguez, 2001 2.3.1671 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 W685 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 = n2, K = n2) D(K = n2, 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=n2, K=n)^T * W^1D (K=n, N=n2)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 y716 // 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())));725 543 } 726 544 … … 766 584 w2 = wl.ToArray(); 767 585 } 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 }778 586 } 779 587 … … 905 713 906 714 // 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; // TODO921 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 // UNFINISHED948 715 public class Spline1dModel : NamedItem, IRegressionModel { 949 privatealglib.spline1dinterpolant interpolant;716 alglib.spline1dinterpolant interpolant; 950 717 951 718 public string TargetVariable { get; set; } … … 958 725 this.TargetVariable = orig.TargetVariable; 959 726 this.VariablesUsedForPrediction = orig.VariablesUsedForPrediction.ToArray(); 727 960 728 this.interpolant = (alglib.spline1dinterpolant)orig.interpolant.make_copy(); 961 729 } 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(); 964 734 this.TargetVariable = targetVar; 965 735 this.VariablesUsedForPrediction = inputs; … … 972 742 public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) { 973 743 return new RegressionSolution(this, (IRegressionProblemData)problemData.Clone()); 744 } 745 746 public double GetEstimatedValue(double xx) { 747 return alglib.spline1dcalc(interpolant, xx); 974 748 } 975 749 … … 979 753 product = product.Zip(dataset.GetDoubleValues(factor, rows), (pi, fi) => pi * fi).ToArray(); 980 754 } 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); 1021 757 } 1022 758 }
Note: See TracChangeset
for help on using the changeset viewer.