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

Last change on this file since 13961 was 13961, checked in by gkronber, 5 years ago

#745: minor change to compile with current trunk

File size: 26.7 KB
Line 
1using System;
2using System.Linq;
3using System.Runtime.InteropServices;
4using HeuristicLab.Algorithms.DataAnalysis;
5using HeuristicLab.Analysis;
6using HeuristicLab.Common;
7using HeuristicLab.Core;
8using HeuristicLab.Data;
9using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
10using HeuristicLab.Optimization;
11using HeuristicLab.Parameters;
12using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
13using HeuristicLab.Problems.DataAnalysis;
14using HeuristicLab.Problems.DataAnalysis.Symbolic;
15using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
16
17namespace HeuristicLab.LibGlmNet {
18  [Item("Elastic-net Linear Regression (LR)", "Linear regression with elastic-net regularization (wrapper for glmnet)")]
19  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 110)]
20  [StorableClass]
21  public sealed class ElasticNetLinearRegression : FixedDataAnalysisAlgorithm<IRegressionProblem> {
22    private const string PenalityParameterName = "Penality";
23    private const string LogLambdaParameterName = "Log10(Lambda)";
24    #region parameters
25    public IFixedValueParameter<DoubleValue> PenalityParameter {
26      get { return (IFixedValueParameter<DoubleValue>)Parameters[PenalityParameterName]; }
27    }
28    public IValueParameter<DoubleValue> LogLambdaParameter {
29      get { return (IValueParameter<DoubleValue>)Parameters[LogLambdaParameterName]; }
30    }
31    #endregion
32    #region properties
33    public double Penality {
34      get { return PenalityParameter.Value.Value; }
35      set { PenalityParameter.Value.Value = value; }
36    }
37    public DoubleValue LogLambda {
38      get { return LogLambdaParameter.Value; }
39      set { LogLambdaParameter.Value = value; }
40    }
41    #endregion
42
43    [StorableConstructor]
44    private ElasticNetLinearRegression(bool deserializing) : base(deserializing) { }
45    private ElasticNetLinearRegression(ElasticNetLinearRegression original, Cloner cloner)
46      : base(original, cloner) {
47    }
48    public ElasticNetLinearRegression()
49      : base() {
50      Problem = new RegressionProblem();
51      Parameters.Add(new FixedValueParameter<DoubleValue>(PenalityParameterName, "Penalty factor for balancing between ridge (0.0) and lasso (1.0) regression", new DoubleValue(0.5)));
52      Parameters.Add(new OptionalValueParameter<DoubleValue>(LogLambdaParameterName, "Optional: the value of lambda for which to calculate an elastic-net solution. lambda == null => calculate the whole path of all lambdas"));
53    }
54
55    [StorableHook(HookType.AfterDeserialization)]
56    private void AfterDeserialization() { }
57
58    public override IDeepCloneable Clone(Cloner cloner) {
59      return new ElasticNetLinearRegression(this, cloner);
60    }
61
62    protected override void Run() {
63      if (LogLambda == null) {
64        CreateSolutionPath();
65      } else {
66        CreateSolution(LogLambda.Value);
67      }
68    }
69
70    private void CreateSolution(double logLambda) {
71      double trainRsq;
72      double testRsq;
73      var coeff = CreateElasticNetLinearRegressionSolution(Problem.ProblemData, Penality, Math.Pow(10, logLambda), out trainRsq, out testRsq);
74      Results.Add(new Result("R² (train)", new DoubleValue(trainRsq)));
75      Results.Add(new Result("R² (test)", new DoubleValue(testRsq)));
76
77      // copied from LR => TODO: reuse code (but skip coefficients = 0.0)
78      ISymbolicExpressionTree tree = new SymbolicExpressionTree(new ProgramRootSymbol().CreateTreeNode());
79      ISymbolicExpressionTreeNode startNode = new StartSymbol().CreateTreeNode();
80      tree.Root.AddSubtree(startNode);
81      ISymbolicExpressionTreeNode addition = new Addition().CreateTreeNode();
82      startNode.AddSubtree(addition);
83
84      int col = 0;
85      foreach (string column in Problem.ProblemData.AllowedInputVariables) {
86        if (!coeff[col].IsAlmost(0.0)) {
87          VariableTreeNode vNode = (VariableTreeNode)new HeuristicLab.Problems.DataAnalysis.Symbolic.Variable().CreateTreeNode();
88          vNode.VariableName = column;
89          vNode.Weight = coeff[col];
90          addition.AddSubtree(vNode);
91        }
92        col++;
93      }
94
95      if (!coeff[coeff.Length - 1].IsAlmost(0.0)) {
96        ConstantTreeNode cNode = (ConstantTreeNode)new Constant().CreateTreeNode();
97        cNode.Value = coeff[coeff.Length - 1];
98        addition.AddSubtree(cNode);
99      }
100
101      SymbolicRegressionSolution solution = new SymbolicRegressionSolution(
102        new SymbolicRegressionModel(Problem.ProblemData.TargetVariable, tree, new SymbolicDataAnalysisExpressionTreeInterpreter()),
103        (IRegressionProblemData)Problem.ProblemData.Clone());
104      solution.Model.Name = "Elastic-net Linear Regression Model";
105      solution.Name = "Elastic-net Linear Regression Solution";
106
107      Results.Add(new Result(solution.Name, solution.Description, solution));
108    }
109
110    private void CreateSolutionPath() {
111      double[] lambda;
112      double[] trainRsq;
113      double[] testRsq;
114      double[,] coeff;
115      double[] intercept;
116      RunElasticNetLinearRegression(Problem.ProblemData, Penality, out lambda, out trainRsq, out testRsq, out coeff, out intercept);
117
118      var coeffTable = new DataTable("Coefficient Paths", "The paths of standarized coefficient values over different lambda values");
119      var nLambdas = lambda.Length;
120      var nCoeff = coeff.GetLength(1);
121      var dataRows = new DataRow[nCoeff];
122      var allowedVars = Problem.ProblemData.AllowedInputVariables.ToArray();
123      for (int i = 0; i < nCoeff; i++) {
124        var coeffId = allowedVars[i];
125        double sigma = Problem.ProblemData.Dataset.GetDoubleValues(coeffId).StandardDeviation();
126        var path = Enumerable.Range(0, nLambdas).Select(r => coeff[r, i] * sigma).ToArray();
127        dataRows[i] = new DataRow(coeffId, coeffId, path);
128        coeffTable.Rows.Add(dataRows[i]);
129      }
130
131      Results.Add(new Result(coeffTable.Name, coeffTable.Description, coeffTable));
132
133      var rsqPlot = new ScatterPlot("R-Squared", "Path of R² values over different lambda values");
134      rsqPlot.VisualProperties.YAxisMaximumAuto = false;
135      rsqPlot.VisualProperties.YAxisMinimumAuto = false;
136      rsqPlot.VisualProperties.XAxisMaximumAuto = false;
137      rsqPlot.VisualProperties.XAxisMinimumAuto = false;
138
139      rsqPlot.VisualProperties.YAxisMinimumFixedValue = 0;
140      rsqPlot.VisualProperties.YAxisMaximumFixedValue = 1.0;
141      rsqPlot.VisualProperties.XAxisTitle = "Log10(Lambda)";
142      rsqPlot.VisualProperties.YAxisTitle = "R²";
143      rsqPlot.Rows.Add(new ScatterPlotDataRow("R² (train)", "Path of R² values over different lambda values", lambda.Zip(trainRsq, (l, r) => new Point2D<double>(Math.Log10(l), r))));
144      rsqPlot.Rows.Add(new ScatterPlotDataRow("R² (test)", "Path of R² values over different lambda values", lambda.Zip(testRsq, (l, r) => new Point2D<double>(Math.Log10(l), r))));
145      if (lambda.Length > 2) {
146        rsqPlot.VisualProperties.XAxisMinimumFixedValue = Math.Floor(Math.Log10(lambda.Last()));
147        rsqPlot.VisualProperties.XAxisMaximumFixedValue = Math.Ceiling(Math.Log10(lambda.Skip(1).First()));
148      }
149      rsqPlot.Rows["R² (train)"].VisualProperties.PointSize = 5;
150      rsqPlot.Rows["R² (test)"].VisualProperties.PointSize = 5;
151
152      Results.Add(new Result(rsqPlot.Name, rsqPlot.Description, rsqPlot));
153    }
154
155    public static double[] CreateElasticNetLinearRegressionSolution(IRegressionProblemData problemData, double penalty, double lambda,
156            out double trainRsq, out double testRsq,
157            double coeffLowerBound = double.NegativeInfinity, double coeffUpperBound = double.PositiveInfinity) {
158      double[] trainRsqs;
159      double[] testRsqs;
160      // run for exactly one lambda
161      var coeffs = CreateElasticNetLinearRegressionSolution(problemData, penalty, new double[] { lambda }, out trainRsqs, out testRsqs, coeffLowerBound, coeffUpperBound);
162      trainRsq = trainRsqs[0];
163      testRsq = testRsqs[0];
164      return coeffs[0];
165    }
166    public static double[][] CreateElasticNetLinearRegressionSolution(IRegressionProblemData problemData, double penalty, double[] lambda,
167            out double[] trainRsq, out double[] testRsq,
168            double coeffLowerBound = double.NegativeInfinity, double coeffUpperBound = double.PositiveInfinity,
169            int maxVars = -1) {
170      // run for multiple user-supplied lambdas
171      double[,] coeff;
172      double[] intercept;
173      RunElasticNetLinearRegression(problemData, penalty, lambda.Length, 1.0, lambda, out lambda, out trainRsq, out testRsq, out coeff, out intercept, coeffLowerBound, coeffUpperBound, maxVars);
174
175      int nRows = intercept.Length;
176      int nCols = coeff.GetLength(1) + 1;
177      double[][] sols = new double[nRows][];
178      for (int solIdx = 0; solIdx < nRows; solIdx++) {
179        sols[solIdx] = new double[nCols];
180        for (int cIdx = 0; cIdx < nCols - 1; cIdx++) {
181          sols[solIdx][cIdx] = coeff[solIdx, cIdx];
182        }
183        sols[solIdx][nCols - 1] = intercept[solIdx];
184      }
185      return sols;
186    }
187
188    public static void RunElasticNetLinearRegression(IRegressionProblemData problemData, double penalty,
189      out double[] lambda, out double[] trainRsq, out double[] testRsq, out double[,] coeff, out double[] intercept,
190      double coeffLowerBound = double.NegativeInfinity, double coeffUpperBound = double.PositiveInfinity,
191      int maxVars = -1
192      ) {
193      double[] userLambda = new double[0];
194      // automatically determine lambda values (maximum 100 different lambda values)
195      RunElasticNetLinearRegression(problemData, penalty, 100, 0.0, userLambda, out lambda, out trainRsq, out testRsq, out coeff, out intercept, coeffLowerBound, coeffUpperBound, maxVars);
196    }
197
198    /// <summary>
199    /// Elastic net with squared-error-loss for dense predictor matrix, runs the full path of all lambdas
200    /// </summary>
201    /// <param name="problemData">Predictor target matrix x and target vector y</param>
202    /// <param name="penalty">Penalty for balance between ridge (0.0) and lasso (1.0) regression</param>
203    /// <param name="nlam">Maximum number of lambda values (default 100)</param>
204    /// <param name="flmin">User control of lambda values (&lt;1.0 => minimum lambda = flmin * (largest lambda value), >= 1.0 => use supplied lambda values</param>
205    /// <param name="ulam">User supplied lambda values</param>
206    /// <param name="lambda">Output lambda values</param>
207    /// <param name="trainRsq">Vector of R² values on the training set for each set of coefficients along the path</param>
208    /// <param name="testRsq">Vector of R² values on the test set for each set of coefficients along the path</param>
209    /// <param name="coeff">Vector of coefficient vectors for each solution along the path</param>
210    /// <param name="intercept">Vector of intercepts for each solution along the path</param>
211    /// <param name="coeffLowerBound">Optional lower bound for all coefficients</param>
212    /// <param name="coeffUpperBound">Optional upper bound for all coefficients</param>
213    /// <param name="maxVars">Maximum allowed number of variables in each solution along the path (-1 => all variables are allowed)</param>
214    private static void RunElasticNetLinearRegression(IRegressionProblemData problemData, double penalty,
215  int nlam, double flmin, double[] ulam, out double[] lambda, out double[] trainRsq, out double[] testRsq, out double[,] coeff, out double[] intercept,
216  double coeffLowerBound = double.NegativeInfinity, double coeffUpperBound = double.PositiveInfinity,
217  int maxVars = -1
218  ) {
219      if (penalty < 0.0 || penalty > 1.0) throw new ArgumentException("0 <= penalty <= 1", "penalty");
220
221      double[,] trainX;
222      double[,] testX;
223
224      double[] trainY;
225      double[] testY;
226      int numTrainObs, numTestObs;
227      int numVars;
228      PrepareData(problemData, out trainX, out trainY, out numTrainObs, out testX, out testY, out numTestObs, out numVars);
229
230      int ka = 1; // => covariance updating algorithm
231      double parm = penalty;
232      double[] w = Enumerable.Repeat(1.0, numTrainObs).ToArray(); // all observations have the same weight
233      int[] jd = new int[1]; // do not force to use any of the variables
234      double[] vp = Enumerable.Repeat(1.0, numVars).ToArray(); // all predictor variables are unpenalized
235      double[,] cl = new double[numVars, 2]; // use the same bounds for all coefficients
236      for (int i = 0; i < numVars; i++) {
237        cl[i, 0] = coeffLowerBound;
238        cl[i, 1] = coeffUpperBound;
239      }
240
241      int ne = maxVars > 0 ? maxVars : numVars;
242      int nx = numVars;
243      double thr = 1.0e-5; // default value as recommended in glmnet
244      int isd = 1; //  => regression on standardized predictor variables
245      int intr = 1;  // => do include intercept in model
246      int maxit = 100000; // default value as recommended in glmnet
247      // outputs
248      int lmu = -1;
249      double[,] ca;
250      int[] ia;
251      int[] nin;
252      int nlp = -99;
253      int jerr = -99;
254
255      elnet(ka, parm, numTrainObs, numVars, trainX, trainY, w, jd, vp, cl, ne, nx, nlam, flmin, ulam, thr, isd, intr, maxit, out lmu, out intercept, out ca, out ia, out nin, out trainRsq, out lambda, out nlp, out jerr);
256
257      testRsq = new double[lmu];
258      coeff = new double[lmu, numVars];
259      for (int solIdx = 0; solIdx < lmu; solIdx++) {
260        // uncompress coefficients of solution
261        int selectedNin = nin[solIdx];
262        double[] coefficients;
263        double[] selectedCa = new double[nx];
264        for (int i = 0; i < nx; i++) {
265          selectedCa[i] = ca[solIdx, i];
266        }
267
268        // apply to test set to calculate test R² values for each lambda step
269        double[] fn;
270        modval(intercept[solIdx], selectedCa, ia, selectedNin, numTestObs, testX, out fn);
271        OnlineCalculatorError error;
272        var r = OnlinePearsonsRCalculator.Calculate(testY, fn, out error);
273        if (error != OnlineCalculatorError.None) r = 0;
274        testRsq[solIdx] = r * r;
275
276        // uncompress coefficients
277        uncomp(numVars, selectedCa, ia, selectedNin, out coefficients);
278        for (int i = 0; i < coefficients.Length; i++) {
279          coeff[solIdx, i] = coefficients[i];
280        }
281      }
282    }
283
284    private static void PrepareData(IRegressionProblemData problemData, out double[,] trainX, out double[] trainY, out int numTrainObs,
285      out double[,] testX, out double[] testY, out int numTestObs, out int numVars) {
286      numVars = problemData.AllowedInputVariables.Count();
287      numTrainObs = problemData.TrainingIndices.Count();
288      numTestObs = problemData.TestIndices.Count();
289
290      trainX = new double[numVars, numTrainObs];
291      trainY = new double[numTrainObs];
292      testX = new double[numVars, numTestObs];
293      testY = new double[numTestObs];
294      var ds = problemData.Dataset;
295      var targetVar = problemData.TargetVariable;
296      // train
297      int rIdx = 0;
298      foreach (var row in problemData.TrainingIndices) {
299        int cIdx = 0;
300        foreach (var var in problemData.AllowedInputVariables) {
301          trainX[cIdx, rIdx] = ds.GetDoubleValue(var, row);
302          cIdx++;
303        }
304        trainY[rIdx] = ds.GetDoubleValue(targetVar, row);
305        rIdx++;
306      }
307      // test
308      rIdx = 0;
309      foreach (var row in problemData.TestIndices) {
310        int cIdx = 0;
311        foreach (var var in problemData.AllowedInputVariables) {
312          testX[cIdx, rIdx] = ds.GetDoubleValue(var, row);
313          cIdx++;
314        }
315        testY[rIdx] = ds.GetDoubleValue(targetVar, row);
316        rIdx++;
317      }
318    }
319
320
321    #region dllimport
322    /// <summary>Wrapper for elnet procedure in glmnet library</summary>
323    /// (see: https://cran.r-project.org/web/packages/glmnet/index.html)
324    ///
325    ///  ka = algorithm flag
326    ///       ka=1 => covariance updating algorithm
327    ///       ka=2 => naive algorithm
328    ///  parm = penalty member index(0&lt;= parm &lt;= 1)
329    ///         = 0.0 => ridge
330    ///  = 1.0 => lasso
331    ///    no = number of observations
332    ///    ni = number of predictor variables
333    ///  y(no) = response vector(overwritten)
334    ///  w(no)= observation weights(overwritten)
335    ///  jd(jd(1)+1) = predictor variable deletion flag
336    ///  jd(1) = 0  => use all variables
337    ///       jd(1) != 0 => do not use variables jd(2)...jd(jd(1)+1)
338    ///  vp(ni) = relative penalties for each predictor variable
339    ///       vp(j) = 0 => jth variable unpenalized
340    ///    cl(2, ni) = interval constraints on coefficient values(overwritten)
341    ///  cl(1, j) = lower bound for jth coefficient value(&lt;= 0.0)
342    ///  cl(2, j) = upper bound for jth coefficient value(>= 0.0)
343    ///  ne = maximum number of variables allowed to enter largest model
344    /// (stopping criterion)
345    ///  nx = maximum number of variables allowed to enter all modesl
346    ///  along path(memory allocation, nx > ne).
347    ///  nlam = (maximum)number of lamda values
348    ///    flmin = user control of lamda values(>=0)
349    ///  flmin&lt; 1.0 => minimum lamda = flmin * (largest lamda value)
350    ///  flmin >= 1.0 => use supplied lamda values(see below)
351    ///  ulam(nlam) = user supplied lamda values(ignored if flmin&lt; 1.0)
352    ///  thr = convergence threshold for each lamda solution.
353    ///  iterations stop when the maximum reduction in the criterion value
354    ///       as a result of each parameter update over a single pass
355    ///       is less than thr times the null criterion value.
356    /// (suggested value, thr= 1.0e-5)
357    ///  isd = predictor variable standarization flag:
358    ///  isd = 0 => regression on original predictor variables
359    ///       isd = 1 => regression on standardized predictor variables
360    ///       Note: output solutions always reference original
361    ///             variables locations and scales.
362    ///    intr = intercept flag
363    ///       intr = 0 / 1 => don't/do include intercept in model
364    ///  maxit = maximum allowed number of passes over the data for all lambda
365    ///  values (suggested values, maxit = 100000)
366    ///
367    ///  output:
368    ///
369    ///    lmu = actual number of lamda values(solutions)
370    ///  a0(lmu) = intercept values for each solution
371    ///  ca(nx, lmu) = compressed coefficient values for each solution
372    ///  ia(nx) = pointers to compressed coefficients
373    ///  nin(lmu) = number of compressed coefficients for each solution
374    ///  rsq(lmu) = R**2 values for each solution
375    ///  alm(lmu) = lamda values corresponding to each solution
376    ///  nlp = actual number of passes over the data for all lamda values
377    ///    jerr = error flag:
378    ///  jerr = 0 => no error
379    ///  jerr > 0 => fatal error - no output returned
380    ///          jerr&lt; 7777 => memory allocation error
381    ///          jerr = 7777 => all used predictors have zero variance
382    ///          jerr = 10000 => maxval(vp) &lt;= 0.0
383    ///  jerr&lt; 0 => non fatal error - partial output:
384    ///  Solutions for larger lamdas (1:(k-1)) returned.
385    ///  jerr = -k => convergence for kth lamda value not reached
386    ///             after maxit(see above) iterations.
387    ///  jerr = -10000 - k => number of non zero coefficients along path
388    ///             exceeds nx(see above) at kth lamda value.
389    ///
390    private static void elnet(
391      int ka,
392      double parm,
393      int no,
394      int ni,
395      double[,] x,
396      double[] y,
397      double[] w,
398      int[] jd,
399      double[] vp,
400      double[,] cl,
401      int ne,
402      int nx,
403      int nlam,
404      double flmin,
405      double[] ulam,
406      double thr,
407      int isd,
408      int intr,
409      int maxit,
410      // outputs
411      out int lmu,
412      out double[] a0,
413      out double[,] ca,
414      out int[] ia,
415      out int[] nin,
416      out double[] rsq,
417      out double[] alm,
418      out int nlp,
419      out int jerr
420      ) {
421      // initialize output values and allocate arrays big enough
422      a0 = new double[nlam];
423      ca = new double[nlam, nx];
424      ia = new int[nx];
425      nin = new int[nlam];
426      rsq = new double[nlam];
427      alm = new double[nlam];
428      nlp = -1;
429      jerr = -1;
430      lmu = -1;
431
432      // load correct version of native dll based on process (x86/x64)
433      if (Environment.Is64BitProcess) {
434        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);
435      } else {
436        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);
437      }
438      //  jerr = error flag:
439      //  jerr = 0 => no error
440      //  jerr > 0 => fatal error -no output returned
441      //  jerr < 7777 => memory allocation error
442      //          jerr = 7777 => all used predictors have zero variance
443      //  jerr = 10000 => maxval(vp) <= 0.0
444      //  jerr < 0 => non fatal error - partial output:
445      //      c Solutions for larger lamdas (1:(k - 1)) returned.
446      //  jerr = -k => convergence for kth lamda value not reached
447      //             after maxit(see above) iterations.
448      //          jerr = -10000 - k => number of non zero coefficients along path
449      //             exceeds nx(see above) at kth lamda value.
450      if (jerr != 0) {
451        if (jerr > 0 && jerr < 7777) throw new InvalidOperationException("glmnet: memory allocation error");
452        else if (jerr == 7777) throw new InvalidOperationException("glmnet: all used predictors have zero variance");
453        else if (jerr == 10000) throw new InvalidOperationException("glmnet: maxval(vp) <= 0.0");
454        else if (jerr < 0 && jerr > -1000) throw new InvalidOperationException(string.Format("glmnet: convergence for {0}th lamda value not reached after maxit iterations ", -jerr));
455        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));
456        else throw new InvalidOperationException(string.Format("glmnet: error {0}", jerr));
457      }
458
459
460      // resize arrays to the capacity that is acutally necessary for the results
461      Array.Resize(ref a0, lmu);
462      Array.Resize(ref nin, lmu);
463      Array.Resize(ref rsq, lmu);
464      Array.Resize(ref alm, lmu);
465    }
466
467    [DllImport("glmnet-x86.dll", EntryPoint = "elnet_", CallingConvention = CallingConvention.Cdecl)]
468    private static extern void elnet_x86(
469      ref int ka,
470      ref double parm,
471      ref int no,
472      ref int ni,
473      double[,] x,
474      double[] y,
475      double[] w,
476      int[] jd,
477      double[] vp,
478      double[,] cl,
479      ref int ne,
480      ref int nx,
481      ref int nlam,
482      ref double flmin,
483      double[] ulam,
484      ref double thr,
485      ref int isd,
486      ref int intr,
487      ref int maxit,
488      // outputs:
489      ref int lmu,
490      [Out] double[] a0,
491      [Out] double[,] ca,
492      [Out] int[] ia,
493      [Out] int[] nin,
494      [Out] double[] rsq,
495      [Out] double[] alm,
496      ref int nlp,
497      ref int jerr
498      );
499    [DllImport("glmnet-x64.dll", EntryPoint = "elnet_", CallingConvention = CallingConvention.Cdecl)]
500    private static extern void elnet_x64(
501      ref int ka,
502      ref double parm,
503      ref int no,
504      ref int ni,
505      double[,] x,
506      double[] y,
507      double[] w,
508      int[] jd,
509      double[] vp,
510      double[,] cl,
511      ref int ne,
512      ref int nx,
513      ref int nlam,
514      ref double flmin,
515      double[] ulam,
516      ref double thr,
517      ref int isd,
518      ref int intr,
519      ref int maxit,
520      // outputs:
521      ref int lmu,
522      [Out] double[] a0,
523      [Out] double[,] ca,
524      [Out] int[] ia,
525      [Out] int[] nin,
526      [Out] double[] rsq,
527      [Out] double[] alm,
528      ref int nlp,
529      ref int jerr
530      );
531
532
533    /// <summary>Wrapper for uncompress coefficient vector for particular solution in glmnet</summary>
534    /// (see: https://cran.r-project.org/web/packages/glmnet/index.html)
535    ///
536    /// call uncomp(ni, ca, ia, nin, a)
537    ///
538    /// input:
539    ///
540    ///    ni = total number of predictor variables
541    ///    ca(nx) = compressed coefficient values for the solution
542    /// ia(nx) = pointers to compressed coefficients
543    /// nin = number of compressed coefficients for the solution
544    ///
545    /// output:
546    ///
547    ///    a(ni) =  uncompressed coefficient vector
548    ///             referencing original variables
549    ///
550    private static void uncomp(int numVars, double[] ca, int[] ia, int nin, out double[] a) {
551      a = new double[numVars];
552      // load correct version of native dll based on process (x86/x64)
553      if (Environment.Is64BitProcess) {
554        uncomp_x64(ref numVars, ca, ia, ref nin, a);
555      } else {
556        uncomp_x86(ref numVars, ca, ia, ref nin, a);
557      }
558    }
559
560    [DllImport("glmnet-x86.dll", EntryPoint = "uncomp_", CallingConvention = CallingConvention.Cdecl)]
561    private static extern void uncomp_x86(ref int numVars, double[] ca, int[] ia, ref int nin, double[] a);
562    [DllImport("glmnet-x64.dll", EntryPoint = "uncomp_", CallingConvention = CallingConvention.Cdecl)]
563    private static extern void uncomp_x64(ref int numVars, double[] ca, int[] ia, ref int nin, double[] a);
564
565    private static void modval(double a0, double[] ca, int[] ia, int nin, int numObs, double[,] x, out double[] fn) {
566      fn = new double[numObs];
567      if (Environment.Is64BitProcess) {
568        modval_x64(ref a0, ca, ia, ref nin, ref numObs, x, fn);
569      } else {
570        modval_x86(ref a0, ca, ia, ref nin, ref numObs, x, fn);
571      }
572    }
573    // evaluate linear model from compressed coefficients and
574    // uncompressed predictor matrix:
575    //
576    // call modval(a0, ca, ia, nin, n, x, f);
577    //   c
578    //   c input:
579    //
580    //    a0 = intercept
581    //    ca(nx) = compressed coefficient values for a solution
582    // ia(nx) = pointers to compressed coefficients
583    // nin = number of compressed coefficients for solution
584    //    n = number of predictor vectors(observations)
585    // x(n, ni) = full(uncompressed) predictor matrix
586    //
587    // output:
588    //
589    //    f(n) = model predictions
590    [DllImport("glmnet-x86.dll", EntryPoint = "modval_", CallingConvention = CallingConvention.Cdecl)]
591    private static extern void modval_x86(ref double a0, double[] ca, int[] ia, ref int nin, ref int numObs, [Out] double[,] x, double[] fn);
592    [DllImport("glmnet-x64.dll", EntryPoint = "modval_", CallingConvention = CallingConvention.Cdecl)]
593    private static extern void modval_x64(ref double a0, double[] ca, int[] ia, ref int nin, ref int numObs, [Out] double[,] x, double[] fn);
594
595    #endregion
596  }
597}
Note: See TracBrowser for help on using the repository browser.