source: branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/SBART.cs @ 15457

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.Common;
using HeuristicLab.Core;
using HeuristicLab.Problems.DataAnalysis;

namespace HeuristicLab.Algorithms.DataAnalysis.Experimental {
  public static class SBART {
    /*
           # A Cubic B-spline Smoothing routine.

    #
    #          The algorithm minimises:
    #
    #      (1/n) * sum ws(i)**2 * (ys(i)-sz(i))**2 + lambda* int ( sz"(xs) )**2 dxs
    #
    #        lambda is a function of the spar which is assumed to be between
    #        0 and 1


            # Input

    #   n               number of data points
    #  ys(n)      vector of length n containing the observations
    #  ws(n)            vector containing the weights given to each data point
    #  xs(n)            vector containing the ordinates of the observations


    #  nk               number of b-spline coefficients to be estimated
    #                   nk <= n+2
    #  knot(nk+4)       vector of knot points defining the cubic b-spline basis.


    #  spar             penalised likelihood smoothing parameter
    #  ispar            indicator saying if spar is supplied or to be estimated
    #  lspar, uspar     lower and upper values for spar 0.,1. are good values
    #  tol              used in Golden Search routine

    #  isetup           setup indicator

    #  icrit            indicator saying which cross validation score
    #       is to be computed

    #  ld4              the leading dimension of abd (ie ld4=4)
    #  ldnk             the leading dimension of p2ip (not referenced)


                    # Output

    #   coef(nk)       vector of spline coefficients
    #   sz(n)          vector of smoothed z-values
    #   lev(n)         vector of leverages
    #   crit           either ordinary of generalized CV score
    #   ier            error indicator
    #                  ier = 0 ___  everything fine
    #                  ier = 1 ___  spar too small or too big
    #                               problem in cholesky decomposition



             # Working arrays/matrix
    #   xwy         X'Wy
    #   hs0,hs1,hs2,hs3   the diagonals of the X'WX matrix
    #   sg0,sg1,sg2,sg3   the diagonals of the Gram matrix
    #   abd(ld4,nk)       [ X'WX+lambda*SIGMA] in diagonal form
    #   p1ip(ld4,nk)       inner products between columns of L inverse
    #   p2ip(ldnk,nk)      all inner products between columns of L inverse
    #                          L'L = [X'WX+lambdaSIGMA]  NOT REFERENCED

    */

    /*
     * sbart(xs,ys,ws,n,knot,nk,
          coef,sz,lev,
          crit,icrit,spar,ispar,lspar,uspar,tol,
          isetup,
          xwy,
          hs0,hs1,hs2,hs3,
          sg0,sg1,sg2,sg3,
          abd,p1ip,p2ip,ld4,ldnk,ier)

    */



    // To build the fortran library (x64) use:
    // > ifort /dll /Qm64 sbart.f /Fesbart_x64.dll
    // check dumpbin /EXPORTS sbart_x64.dll
    // and   dumpbin /IMPORTS sbart_x64.dll
    [DllImport("sbart_x64.dll", CallingConvention = CallingConvention.Cdecl, EntryPoint = "sbart")]
    public static extern void sbart_x64(
     float[] xs, float[] ys, float[] ws, ref int n, float[] knot, ref int nk,
     float[] coeff, float[] sz, float[] lev,
     ref float crit, ref int icrit, ref float spar, ref int ispar, ref float lspar, ref float uspar, ref float tol,
     ref int isetup,
     float[] xwy,
     float[] hs0, float[] hs1, float[] hs2, float[] hs3,
     float[] sg0, float[] sg1, float[] sg2, float[] sg3,
     float[,] abd, float[,] p1ip, float[,] p2ip, ref int ld4, ref int ldnk, ref int ier);

    [DllImport("sbart_x86.dll", CallingConvention = CallingConvention.Cdecl, EntryPoint = "sbart")]
    public static extern void sbart_x86();


    [DllImport("sbart_x64.dll", CallingConvention = CallingConvention.Cdecl, EntryPoint = "sknotl")]
    public static extern void sknotl_x64(float[] x, ref int n, float[] knot, ref int k);
    [DllImport("sbart_x64.dll", CallingConvention = CallingConvention.Cdecl, EntryPoint = "setreg")]
    public static extern void setreg_x64(float[] x, float[] y, float[] w, ref int n, float[] xw, ref int nx, ref float min, ref float range, float[] knot, ref int nk);

    public class SBART_Report {
      public double smoothingParameter;
      public double gcv;
      public double[] leverage;     
    }

    public static IRegressionModel CalculateSBART(double[] x, double[] y,
      string targetVariable, string[] inputVars, float smoothingParameter,
      out SBART_Report report) {
      var w = Enumerable.Repeat(1.0, x.Length).ToArray();

      int n = x.Length;
      float[] leverage = new float[n];
      float[] y_smoothed = new float[n];
      int ic = n - 1;
      int ier = -99;
      int nk = n;
      float[] knots = new float[nk + 6];

      float crit = -99.0f;
      int icrit = 0; // 0..don't calc CV,  1 .. GCV, 2 CV

      // float smoothingParameter = -99.0f;
      int smoothingParameterIndicator = 1; 
      float lowerSmoothingParameter = 0f;
      float upperSmoothingParameter = 1.0f;
      float tol = 0.01f;
      int isetup = 0; // not setup?

      float min = -99.0f;
      float range = -99.0f;

      if (Environment.Is64BitProcess) {
        float[] xw = new float[n];
        int nx = -99;
        float[] xs = x.Select(xi => (float)xi).ToArray();
        float[] ys = y.Select(yi => (float)yi).ToArray();
        float[] ws = w.Select(wi => (float)wi).ToArray();
        SBART.setreg_x64(xs, ys, ws,
          ref n, xw, ref nx, ref min, ref range, knots, ref nk);
        if(nx < 4) {
          report = new SBART_Report();
          report.leverage = new double[0];
          return new ConstantModel(ys.Take(nx).Average(), targetVariable);
        }

        float[] coeff = new float[nk];

        // working arrays for sbart
        float[] xwy = new float[nk];
        float[] hs0 = new float[nk];
        float[] hs1 = new float[nk];
        float[] hs2 = new float[nk];
        float[] hs3 = new float[nk];
        float[] sg0 = new float[nk];
        float[] sg1 = new float[nk];
        float[] sg2 = new float[nk];
        float[] sg3 = new float[nk];
        int ld4 = 4;
        float[,] adb = new float[nk, ld4];

        float[,] p1ip = new float[nk, ld4];
        int ldnk = nk + 4;
        float[,] p2ip = new float[nk, ldnk];

        SBART.sbart_x64(xs, ys, ws, ref nx, knots, ref nk, coeff, y_smoothed, leverage, ref crit, ref icrit, 
          ref smoothingParameter, ref smoothingParameterIndicator, ref lowerSmoothingParameter, ref upperSmoothingParameter,
          ref tol, ref isetup, 
          xwy, hs0, hs1, hs2, hs3, sg0, sg1, sg2, sg3, adb, p1ip, p2ip, ref ld4, ref ldnk, ref ier);

        if (ier > 0) throw new ArgumentException(ier.ToString());

        report = new SBART_Report();
        report.gcv = crit;
        report.smoothingParameter = smoothingParameter;
        report.leverage = leverage.Select(li => (double)li).ToArray();

        return new BartRegressionModel(xs.Take(nx).ToArray(), y_smoothed.Take(nx).ToArray(), targetVariable, inputVars, min, range);

      } else {
        throw new NotSupportedException();
      }

    }

    public class BartRegressionModel : NamedItem, IRegressionModel {
      private float[] x;
      private float[] y;
      private double min;
      private double range;
      public string TargetVariable { get; set; }

      private string[] variablesUsedForPrediction;
      public IEnumerable<string> VariablesUsedForPrediction {
        get {
          return variablesUsedForPrediction;
        }
      }

      public BartRegressionModel(BartRegressionModel orig, Cloner cloner) {
        this.x = orig.x;
        this.y = orig.y;
        this.min = orig.min;
        this.range = orig.range;
        this.TargetVariable = orig.TargetVariable;
        this.variablesUsedForPrediction = orig.variablesUsedForPrediction;
      }
      public BartRegressionModel(float[] x, float[] y_smoothed, string targetVariable, string[] inputVars, double min, double range) {
        this.variablesUsedForPrediction = inputVars;
        this.TargetVariable = targetVariable;
        this.x = x;
        this.y = y_smoothed;
        this.range = range;
        this.min = min;
      }

      public event EventHandler TargetVariableChanged;

      public IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
        return new RegressionSolution(this, (IRegressionProblemData) problemData.Clone());
      }

      public double GetEstimatedValue(double xx) {
        xx = (xx - min) / range;
        // piecewise linear approximation
        int n = x.Length;
        if (xx <= x[0]) {
          double h = xx - x[0];
          return h * (y[1] - y[0]) / (x[1] - x[0]) + y[0];
        } else if (xx >= x[n-1]) {
          double h = xx - x[n-1];
          return h * (y[n-1] - y[n-2]) / (x[n-1] - x[n-2]) + y[n-1];
        } else {
          // binary search
          int lower = 0;
          int upper = n-1;
          while (true) {
            if (upper < lower) throw new InvalidProgramException();
            int i = lower + (upper - lower) / 2;
            if (x[i] <= xx && xx < x[i + 1]) {
              double h = xx - x[i];
              double k = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
              return h * k + y[i];
            } else if (xx < x[i]) {
              upper = i - 1;
            } else {
              lower = i + 1;
            }
          }
        }
        return 0.0;
      }

      public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
        foreach(var x in dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows)) {
          yield return GetEstimatedValue(x);
        }
      }

      public override IDeepCloneable Clone(Cloner cloner) {
        return new BartRegressionModel(this, cloner);
      }
    }

  }
}