source: branches/HeuristicLab.Algorithms.DataAnalysis.Glmnet/3.4/ElasticNetLinearRegression.cs @ 13927

using System;
using System.Linq;
using System.Runtime.InteropServices;
using HeuristicLab.Algorithms.DataAnalysis;
using HeuristicLab.Analysis;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;

namespace HeuristicLab.LibGlmNet {
  [Item("Elastic-net Linear Regression (LR)", "Linear regression with elastic-net regularization (wrapper for glmnet)")]
  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 110)]
  [StorableClass]
  public sealed class ElasticNetLinearRegression : FixedDataAnalysisAlgorithm<IRegressionProblem> {
    private const string PenalityParameterName = "Penality";
    #region parameters
    public IFixedValueParameter<DoubleValue> PenalityParameter {
      get { return (IFixedValueParameter<DoubleValue>)Parameters[PenalityParameterName]; }
    }
    #endregion
    #region properties
    public double Penality {
      get { return PenalityParameter.Value.Value; }
      set { PenalityParameter.Value.Value = value; }
    }
    #endregion

    [StorableConstructor]
    private ElasticNetLinearRegression(bool deserializing) : base(deserializing) { }
    private ElasticNetLinearRegression(ElasticNetLinearRegression original, Cloner cloner)
      : base(original, cloner) {
    }
    public ElasticNetLinearRegression() : base() {
      Problem = new RegressionProblem();
      Parameters.Add(new FixedValueParameter<DoubleValue>(PenalityParameterName, "Penalty factor for balancing between ridge (0.0) and lasso (1.0) regression", new DoubleValue(0.5)));
    }

    [StorableHook(HookType.AfterDeserialization)]
    private void AfterDeserialization() { }

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

    protected override void Run() {
      double[] lambda;
      double[] rsq;
      double[,] coeff;
      double[] intercept;
      RunElasticNetLinearRegression(Problem.ProblemData, Penality, out lambda, out rsq, out coeff, out intercept);

      var coeffTable = new DataTable("Coefficient Paths", "The paths of coefficient values over different lambda values");
      var nLambdas = lambda.Length;
      var nCoeff = coeff.GetLength(1);
      var dataRows = new DataRow[nCoeff];
      for (int i = 0; i < nCoeff; i++) {
        var coeffId = string.Format("c_{0}", i);
        dataRows[i] = new DataRow(coeffId, coeffId, Enumerable.Range(0, nLambdas).Select(r => coeff[r, i]).ToArray());
        coeffTable.Rows.Add(dataRows[i]);
      }

      Results.Add(new Result(coeffTable.Name, coeffTable.Description, coeffTable));

      var rsqTable = new DataTable("R-Squared", "Path of R² values over different lambda values");
      rsqTable.Rows.Add(new DataRow("R-Squared", "Path of R² values over different lambda values", rsq));
      Results.Add(new Result(rsqTable.Name, rsqTable.Description, coeffTable));
    }

    //     private static void TestTower() {
    //       var prov = new HeuristicLab.Problems.Instances.DataAnalysis.RegressionRealWorldInstanceProvider();
    //       var towerDesc = prov.GetDataDescriptors().Where(dd => dd.Name.Contains("Tower")).Single();
    //       try {
    //         var problemData = prov.LoadData(towerDesc);
    // 
    //         double[] lambda;
    //         double[] rsq;
    //         double[,] coeff;
    //         double[] intercept;
    //         RunElasticNetLinearRegression(problemData, 0.5, out lambda, out rsq, out coeff, out intercept);
    // 
    //         for (int i = 0; i < rsq.Length; i++) {
    //           Console.WriteLine(rsq[i]);
    //         }
    // 
    //         double rsq0;
    //         CreateElasticNetLinearRegressionSolution(problemData, 0.5, lambda.Skip(20).First(), out rsq0);
    // 
    //         CreateElasticNetLinearRegressionSolution(problemData, 0.5, lambda.Skip(20).Take(10).ToArray(), out rsq);
    // 
    // 
    //       }
    //       catch (Exception e) {
    //       }
    //     }

    public static double[] CreateElasticNetLinearRegressionSolution(IRegressionProblemData problemData, double penalty, double lambda,
            out double rsq,
            double coeffLowerBound = double.NegativeInfinity, double coeffUpperBound = double.PositiveInfinity) {
      double[] rsqs;
      // run for exactly one lambda
      var coeffs = CreateElasticNetLinearRegressionSolution(problemData, penalty, new double[] { lambda }, out rsqs, coeffLowerBound, coeffUpperBound);
      rsq = rsqs[0];
      return coeffs[0];
    }
    public static double[][] CreateElasticNetLinearRegressionSolution(IRegressionProblemData problemData, double penalty, double[] lambda,
            out double[] rsq,
            double coeffLowerBound = double.NegativeInfinity, double coeffUpperBound = double.PositiveInfinity,
            int maxVars = -1) {
      // run for multiple user-supplied lambdas
      double[,] coeff;
      double[] intercept;
      RunElasticNetLinearRegression(problemData, penalty, lambda.Length, 1.0, lambda, out lambda, out rsq, out coeff, out intercept, coeffLowerBound, coeffUpperBound, maxVars);

      int nRows = intercept.Length;
      int nCols = coeff.GetLength(1) + 1;
      double[][] sols = new double[nRows][];
      for (int solIdx = 0; solIdx < nRows; solIdx++) {
        sols[solIdx] = new double[nCols];
        for (int cIdx = 0; cIdx < nCols - 1; cIdx++) {
          sols[solIdx][cIdx] = coeff[solIdx, cIdx];
        }
        sols[solIdx][nCols - 1] = intercept[solIdx];
      }
      return sols;
    }

    public static void RunElasticNetLinearRegression(IRegressionProblemData problemData, double penalty,
      out double[] lambda, out double[] rsq, out double[,] coeff, out double[] intercept,
      double coeffLowerBound = double.NegativeInfinity, double coeffUpperBound = double.PositiveInfinity,
      int maxVars = -1
      ) {
      double[] userLambda = new double[0];
      // automatically determine lambda values (maximum 100 different lambda values)
      RunElasticNetLinearRegression(problemData, penalty, 100, 0.0, userLambda, out lambda, out rsq, out coeff, out intercept, coeffLowerBound, coeffUpperBound, maxVars);
    }

    /// <summary>
    /// Elastic net with squared-error-loss for dense predictor matrix, runs the full path of all lambdas
    /// </summary>
    /// <param name="problemData">Predictor target matrix x and target vector y</param>
    /// <param name="penalty">Penalty for balance between ridge (0.0) and lasso (1.0) regression</param>
    /// <param name="nlam">Maximum number of lambda values (default 100)</param>
    /// <param name="flmin">User control of lambda values (&lt;1.0 => minimum lambda = flmin * (largest lambda value), >= 1.0 => use supplied lambda values</param>
    /// <param name="ulam">User supplied lambda values</param>
    /// <param name="lambda">Output lambda values</param>
    /// <param name="rsq">Vector of R² values for each set of coefficients along the path</param>
    /// <param name="coeff">Vector of coefficient vectors for each solution along the path</param>
    /// <param name="intercept">Vector of intercepts for each solution along the path</param>
    /// <param name="coeffLowerBound">Optional lower bound for all coefficients</param>
    /// <param name="coeffUpperBound">Optional upper bound for all coefficients</param>
    /// <param name="maxVars">Maximum allowed number of variables in each solution along the path (-1 => all variables are allowed)</param>
    private static void RunElasticNetLinearRegression(IRegressionProblemData problemData, double penalty,
  int nlam, double flmin, double[] ulam, out double[] lambda, out double[] rsq, out double[,] coeff, out double[] intercept,
  double coeffLowerBound = double.NegativeInfinity, double coeffUpperBound = double.PositiveInfinity,
  int maxVars = -1
  ) {
      if (penalty < 0.0 || penalty > 1.0) throw new ArgumentException("0 <= penalty <= 1", "penalty");

      double[,] x;
      double[] y;
      int numObs;
      int numVars;
      PrepareData(problemData, out x, out y, out numObs, out numVars);

      int ka = 1; // => covariance updating algorithm
      double parm = penalty;
      double[] w = Enumerable.Repeat(1.0, numObs).ToArray(); // all observations have the same weight
      int[] jd = new int[1]; // do not force to use any of the variables
      double[] vp = Enumerable.Repeat(1.0, numVars).ToArray(); // all predictor variables are unpenalized
      double[,] cl = new double[numVars, 2]; // use the same bounds for all coefficients
      for (int i = 0; i < numVars; i++) {
        cl[i, 0] = coeffLowerBound;
        cl[i, 1] = coeffUpperBound;
      }

      int ne = maxVars > 0 ? maxVars : numVars;
      int nx = numVars;
      double thr = 1.0e-5; // default value as recommended in glmnet
      int isd = 0; //  => regression on original predictor variables
      int intr = 1;  // => do include intercept in model
      int maxit = 100000; // default value as recommended in glmnet
      // outputs
      int lmu = -1;
      double[,] ca;
      int[] ia;
      int[] nin;
      int nlp = -99;
      int jerr = -99;

      elnet(ka, parm, numObs, numVars, x, y, w, jd, vp, cl, ne, nx, nlam, flmin, ulam, thr, isd, intr, maxit, out lmu, out intercept, out ca, out ia, out nin, out rsq, out lambda, out nlp, out jerr);

      coeff = new double[lmu, numVars];
      for (int solIdx = 0; solIdx < lmu; solIdx++) {
        // uncompress coefficients of solution
        int selectedNin = nin[solIdx];
        double[] coefficients;
        double[] selectedCa = new double[nx];
        for (int i = 0; i < nx; i++) selectedCa[i] = ca[solIdx, i];

        uncomp(numVars, selectedCa, ia, selectedNin, out coefficients);
        for (int i = 0; i < coefficients.Length; i++) {
          coeff[solIdx, i] = coefficients[i];
        }
      }
    }

    private static void PrepareData(IRegressionProblemData problemData, out double[,] x, out double[] y, out int numObs, out int numVars) {
      numVars = problemData.AllowedInputVariables.Count();
      numObs = problemData.TrainingIndices.Count();

      x = new double[numVars, numObs];
      y = new double[numObs];
      var ds = problemData.Dataset;
      var targetVar = problemData.TargetVariable;
      int rIdx = 0;
      foreach (var row in problemData.TrainingIndices) {
        int cIdx = 0;
        foreach (var var in problemData.AllowedInputVariables) {
          x[cIdx, rIdx] = ds.GetDoubleValue(var, row);
          cIdx++;
        }
        y[rIdx] = ds.GetDoubleValue(targetVar, row);
        rIdx++;
      }
    }


    #region dllimport
    /// <summary>Wrapper for elnet procedure in glmnet library</summary>
    /// (see: https://cran.r-project.org/web/packages/glmnet/index.html)
    /// 
    ///  ka = algorithm flag
    ///       ka=1 => covariance updating algorithm
    ///       ka=2 => naive algorithm
    ///  parm = penalty member index(0&lt;= parm &lt;= 1)
    ///         = 0.0 => ridge
    ///  = 1.0 => lasso
    ///    no = number of observations
    ///    ni = number of predictor variables
    ///  y(no) = response vector(overwritten)
    ///  w(no)= observation weights(overwritten)
    ///  jd(jd(1)+1) = predictor variable deletion flag
    ///  jd(1) = 0  => use all variables
    ///       jd(1) != 0 => do not use variables jd(2)...jd(jd(1)+1)
    ///  vp(ni) = relative penalties for each predictor variable
    ///       vp(j) = 0 => jth variable unpenalized
    ///    cl(2, ni) = interval constraints on coefficient values(overwritten)
    ///  cl(1, j) = lower bound for jth coefficient value(&lt;= 0.0)
    ///  cl(2, j) = upper bound for jth coefficient value(>= 0.0)
    ///  ne = maximum number of variables allowed to enter largest model
    /// (stopping criterion)
    ///  nx = maximum number of variables allowed to enter all modesl
    ///  along path(memory allocation, nx > ne).
    ///  nlam = (maximum)number of lamda values
    ///    flmin = user control of lamda values(>=0)
    ///  flmin&lt; 1.0 => minimum lamda = flmin * (largest lamda value)
    ///  flmin >= 1.0 => use supplied lamda values(see below)
    ///  ulam(nlam) = user supplied lamda values(ignored if flmin&lt; 1.0)
    ///  thr = convergence threshold for each lamda solution.
    ///  iterations stop when the maximum reduction in the criterion value
    ///       as a result of each parameter update over a single pass
    ///       is less than thr times the null criterion value.
    /// (suggested value, thr= 1.0e-5)
    ///  isd = predictor variable standarization flag:
    ///  isd = 0 => regression on original predictor variables
    ///       isd = 1 => regression on standardized predictor variables
    ///       Note: output solutions always reference original
    ///             variables locations and scales.
    ///    intr = intercept flag
    ///       intr = 0 / 1 => don't/do include intercept in model
    ///  maxit = maximum allowed number of passes over the data for all lambda
    ///  values (suggested values, maxit = 100000)
    /// 
    ///  output:
    /// 
    ///    lmu = actual number of lamda values(solutions)
    ///  a0(lmu) = intercept values for each solution
    ///  ca(nx, lmu) = compressed coefficient values for each solution
    ///  ia(nx) = pointers to compressed coefficients
    ///  nin(lmu) = number of compressed coefficients for each solution
    ///  rsq(lmu) = R**2 values for each solution
    ///  alm(lmu) = lamda values corresponding to each solution
    ///  nlp = actual number of passes over the data for all lamda values
    ///    jerr = error flag:
    ///  jerr = 0 => no error
    ///  jerr > 0 => fatal error - no output returned
    ///          jerr&lt; 7777 => memory allocation error
    ///          jerr = 7777 => all used predictors have zero variance
    ///          jerr = 10000 => maxval(vp) &lt;= 0.0
    ///  jerr&lt; 0 => non fatal error - partial output:
    ///  Solutions for larger lamdas (1:(k-1)) returned.
    ///  jerr = -k => convergence for kth lamda value not reached
    ///             after maxit(see above) iterations.
    ///  jerr = -10000 - k => number of non zero coefficients along path
    ///             exceeds nx(see above) at kth lamda value.
    /// 
    private static void elnet(
      int ka,
      double parm,
      int no,
      int ni,
      double[,] x,
      double[] y,
      double[] w,
      int[] jd,
      double[] vp,
      double[,] cl,
      int ne,
      int nx,
      int nlam,
      double flmin,
      double[] ulam,
      double thr,
      int isd,
      int intr,
      int maxit,
      // outputs
      out int lmu,
      out double[] a0,
      out double[,] ca,
      out int[] ia,
      out int[] nin,
      out double[] rsq,
      out double[] alm,
      out int nlp,
      out int jerr
      ) {
      // initialize output values and allocate arrays big enough
      a0 = new double[nlam];
      ca = new double[nlam, nx];
      ia = new int[nx];
      nin = new int[nlam];
      rsq = new double[nlam];
      alm = new double[nlam];
      nlp = -1;
      jerr = -1;
      lmu = -1;

      // load correct version of native dll based on process (x86/x64)
      if (Environment.Is64BitProcess) {
        elnet_x64(ref ka, ref parm, ref no, ref ni, x, y, w, jd, vp, cl, ref ne, ref ni, ref nlam, ref flmin, ulam, ref thr, ref isd, ref intr, ref maxit, ref lmu, a0, ca, ia, nin, rsq, alm, ref nlp, ref jerr);
      } else {
        elnet_x86(ref ka, ref parm, ref no, ref ni, x, y, w, jd, vp, cl, ref ne, ref ni, ref nlam, ref flmin, ulam, ref thr, ref isd, ref intr, ref maxit, ref lmu, a0, ca, ia, nin, rsq, alm, ref nlp, ref jerr);
      }
      //  jerr = error flag:
      //  jerr = 0 => no error
      //  jerr > 0 => fatal error -no output returned
      //  jerr < 7777 => memory allocation error
      //          jerr = 7777 => all used predictors have zero variance
      //  jerr = 10000 => maxval(vp) <= 0.0
      //  jerr < 0 => non fatal error - partial output:
      //      c Solutions for larger lamdas (1:(k - 1)) returned.
      //  jerr = -k => convergence for kth lamda value not reached
      //             after maxit(see above) iterations.
      //          jerr = -10000 - k => number of non zero coefficients along path
      //             exceeds nx(see above) at kth lamda value.
      if (jerr != 0) {
        if (jerr > 0 && jerr < 7777) throw new InvalidOperationException("glmnet: memory allocation error");
        else if (jerr == 7777) throw new InvalidOperationException("glmnet: all used predictors have zero variance");
        else if (jerr == 10000) throw new InvalidOperationException("glmnet: maxval(vp) <= 0.0");
        else if (jerr < 0 && jerr > -1000) throw new InvalidOperationException(string.Format("glmnet: convergence for {0}th lamda value not reached after maxit iterations ", -jerr));
        else if (jerr <= -10000) throw new InvalidOperationException(string.Format("glmnet: number of non zero coefficients along path exceeds number of maximally allowed variables (nx) at {0}th lamda value", -jerr - 10000));
        else throw new InvalidOperationException(string.Format("glmnet: error {0}", jerr));
      }


      // resize arrays to the capacity that is acutally necessary for the results
      Array.Resize(ref a0, lmu);
      Array.Resize(ref nin, lmu);
      Array.Resize(ref rsq, lmu);
      Array.Resize(ref alm, lmu);
    }

    [DllImport("glmnet-x86.dll", EntryPoint = "elnet_", CallingConvention = CallingConvention.Cdecl)]
    private static extern void elnet_x86(
      ref int ka,
      ref double parm,
      ref int no,
      ref int ni,
      double[,] x,
      double[] y,
      double[] w,
      int[] jd,
      double[] vp,
      double[,] cl,
      ref int ne,
      ref int nx,
      ref int nlam,
      ref double flmin,
      double[] ulam,
      ref double thr,
      ref int isd,
      ref int intr,
      ref int maxit,
      // outputs:
      ref int lmu,
      [Out] double[] a0,
      [Out] double[,] ca,
      [Out] int[] ia,
      [Out] int[] nin,
      [Out] double[] rsq,
      [Out] double[] alm,
      ref int nlp,
      ref int jerr
      );
    [DllImport("glmnet-x64.dll", EntryPoint = "elnet_", CallingConvention = CallingConvention.Cdecl)]
    private static extern void elnet_x64(
      ref int ka,
      ref double parm,
      ref int no,
      ref int ni,
      double[,] x,
      double[] y,
      double[] w,
      int[] jd,
      double[] vp,
      double[,] cl,
      ref int ne,
      ref int nx,
      ref int nlam,
      ref double flmin,
      double[] ulam,
      ref double thr,
      ref int isd,
      ref int intr,
      ref int maxit,
      // outputs:
      ref int lmu,
      [Out] double[] a0,
      [Out] double[,] ca,
      [Out] int[] ia,
      [Out] int[] nin,
      [Out] double[] rsq,
      [Out] double[] alm,
      ref int nlp,
      ref int jerr
      );


    /// <summary>Wrapper for uncompress coefficient vector for particular solution in glmnet</summary>
    /// (see: https://cran.r-project.org/web/packages/glmnet/index.html)
    ///
    /// call uncomp(ni, ca, ia, nin, a)
    ///
    /// input:
    ///
    ///    ni = total number of predictor variables
    ///    ca(nx) = compressed coefficient values for the solution
    /// ia(nx) = pointers to compressed coefficients
    /// nin = number of compressed coefficients for the solution
    ///
    /// output:
    ///
    ///    a(ni) =  uncompressed coefficient vector
    ///             referencing original variables
    ///
    private static void uncomp(int numVars, double[] ca, int[] ia, int nin, out double[] a) {
      a = new double[numVars];
      // load correct version of native dll based on process (x86/x64)
      if (Environment.Is64BitProcess) {
        uncomp_x64(ref numVars, ca, ia, ref nin, a);
      } else {
        uncomp_x86(ref numVars, ca, ia, ref nin, a);
      }
    }

    [DllImport("glmnet-x86.dll", EntryPoint = "uncomp_", CallingConvention = CallingConvention.Cdecl)]
    private static extern void uncomp_x86(ref int numVars, double[] ca, int[] ia, ref int nin, double[] a);
    [DllImport("glmnet-x64.dll", EntryPoint = "uncomp_", CallingConvention = CallingConvention.Cdecl)]
    private static extern void uncomp_x64(ref int numVars, double[] ca, int[] ia, ref int nin, double[] a);

    #endregion
  }
}