Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/HEAL.VarPro.csproj
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/HEAL.VarPro.csproj	(revision 18003)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/HEAL.VarPro.csproj	(revision 18003)
@@ -0,0 +1,70 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <PropertyGroup>
+    <TargetFramework>netstandard2.1</TargetFramework>
+    <SignAssembly>false</SignAssembly>
+    <AssemblyOriginatorKeyFile>HeuristicLab.snk</AssemblyOriginatorKeyFile>
+    <DelaySign>false</DelaySign>
+  </PropertyGroup>
+
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|AnyCPU'">
+    <AllowUnsafeBlocks>true</AllowUnsafeBlocks>
+  </PropertyGroup>
+
+  <ItemGroup>
+    <Compile Remove="MKLTest.cs" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x64\libgcc_s_seh-1.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x64\libgfortran-3.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x64\libopenblas.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x64\libquadmath-0.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x64\MathNet.Numerics.OpenBLAS.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x86\libgcc_s_sjlj-1.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x86\libgfortran-3.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x86\libopenblas.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x86\libquadmath-0.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Content Remove="C:\Users\gabri\.nuget\packages\mathnet.numerics.openblas.win\0.2.0\build\x86\MathNet.Numerics.OpenBLAS.dll" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Reference Include="alglibnet2">
+      <HintPath>alglibnet2.dll</HintPath>
+    </Reference>
+  </ItemGroup>
+
+  <ItemGroup>
+    <None Update="alglibnet2.dll">
+      <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
+    </None>
+  </ItemGroup>
+
+</Project>
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/VariableProjection.cs
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/VariableProjection.cs	(revision 18003)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/VariableProjection.cs	(revision 18003)
@@ -0,0 +1,363 @@
+﻿using System;
+using System.Linq;
+using System.Threading;
+
+namespace HEAL.VarPro {
+  public class VariableProjection {
+    public delegate void FeatureFunc(double[] alpha, ref double[,] phi);
+    public delegate void JacobianFunc(double[] alpha, ref double[,] Jac, ref int[,] ind);
+    public class Report {
+      public int iter;
+      public double residNorm;
+      public double residNormSqr;
+      public double lambda;
+      public double w;
+      public double lineSearchStep;
+      public double[] alpha;
+      public double[] coeff;
+      public double gradNorm;
+    }
+
+    private VariableProjection(FeatureFunc Phi, JacobianFunc Jac, int nAlpha, double[] y, Action<Report, CancellationToken> iterationCallback = null) {
+      this.y = y;
+      this.Phi = Phi;
+      this.Jac = Jac;
+      this.iterationCallback = iterationCallback;
+      this.r = new double[y.Length]; // residuals
+      this.m = y.Length;
+      this.nAlpha = nAlpha;
+      J = new double[m, nAlpha];
+      Jac2 = new double[m, nAlpha];
+      grad = new double[nAlpha];
+      U = new double[m, m];
+    }
+
+    public bool WGCV { get; set; } = true;
+
+    private readonly double[] y;
+    private readonly FeatureFunc Phi;
+    private readonly JacobianFunc Jac;
+
+
+    private Action<Report, CancellationToken> iterationCallback;
+    private double eps;
+    private int maxIters;
+    private int nAlpha;
+    private int m;
+
+    // buffers
+    private double[] r, s, alpha, coeff, grad;
+    private double[,] F, dF, J, U, VT, Jac2;
+    private int[,] ind;
+
+    public static void Fit(FeatureFunc Phi, JacobianFunc Jac, double[] y, double[] alpha, out double[] coeff, out Report report,
+      int maxIters = 20, Action<Report, CancellationToken> iterationCallback = null, bool useWGCV = true, double eps = 1e-16) {
+      // Step 1: Choose theta_N_0 \in R^k, a small positive number eps and 
+      //         the maximum iteration N. Set k = 0;
+
+      var varpro = new VariableProjection(Phi, Jac, alpha.Length, y);
+      varpro.WGCV = useWGCV;
+      varpro.iterationCallback = iterationCallback;
+      varpro.eps = eps;
+      varpro.maxIters = maxIters;
+      varpro.alpha = alpha; // don't copy alpha. Fit() updates alpha
+
+      varpro.FitInternal(out report);
+      coeff = varpro.coeff;
+    }
+    private void FitInternal(out Report report) {
+      int k = 0;
+      int rank;
+      var w = 1.0;
+      var lambda = 0.0;
+      C(y, alpha, lambda, w, ref coeff, ref r, ref U, ref s, ref VT, out rank, out var nextLambda, out var nextW, out var resNormSqr, out var coeffNormSqr);
+      lambda = nextLambda;
+      w = nextW;
+      var cancellationTokenSource = new CancellationTokenSource();
+      report = null;
+      while (k <= maxIters) {
+        CalculateJacobian(Jac, alpha, coeff, U, s, VT, rank, r, ref J);
+
+        // Step 4: Calculate the gradient of the objective function (17)
+        //         with resprect to the nonlinear parameters
+        //         ∇C = - J^T r
+        alglib.rmatrixgemv(nAlpha, m, -1.0, J, 0, 0, 1, r, 0, 0.0, ref grad, 0);
+        // Step 5: if ||∇C||^2 < eps, terminate the algorithm
+        var gradNorm = 0.0;
+        for (int i = 0; i < grad.Length; i++) gradNorm += grad[i] * grad[i];
+
+        report = new Report() {
+          iter = k,
+          residNorm = Math.Sqrt(resNormSqr),
+          lambda = lambda,
+          residNormSqr = resNormSqr,
+          w = w,
+          alpha = (double[])alpha.Clone(),
+          coeff = (double[])coeff.Clone(),
+          gradNorm = gradNorm
+        };
+
+        if (gradNorm < eps) break;
+
+        // Step 6: Calculate the search direction d using (20), then
+        //         update the nonlinear parameters theta_N by the line
+        //         search procedure and (21)
+        var d = SolveLS(J, r);
+
+
+        #region line search
+        // Backtracking line search (Boyd, Convex optimization, Chapter 9)
+        // given descent direction d, for f a x \in dom(f), a \in (0, 0.5), b \in (0, 1)
+
+        // line search parameters a and b
+        // recommended settings from Boyd, Chapter 9
+        var a = 0.1;
+        var b = 0.5;
+
+        var t = 1.0;
+
+        var predictedChange = 0.0;
+        for (int i = 0; i < grad.Length; i++) predictedChange += grad[i] * d[i];
+        if (predictedChange >= 0) throw new InvalidProgramException("Descent direction does not improve objective. Check your Jacobian calculation");
+
+        var newAlpha = (double[])alpha.Clone(); // TODO allocate once
+        for (int i = 0; i < newAlpha.Length; i++) newAlpha[i] -= t * d[i];
+
+        C(y, newAlpha, lambda, w, ref coeff, ref r, ref U, ref s, ref VT, out rank, out nextLambda, out nextW, out var newResNormSqr, out var newCoeffNormSqr);
+        var fx = resNormSqr + nextLambda * coeffNormSqr;
+        var fx_new = newResNormSqr + nextLambda * newCoeffNormSqr;
+
+        while (t > 0 && fx_new > fx + a * t * predictedChange) {
+          t = b * t;
+          for (int i = 0; i < newAlpha.Length; i++) newAlpha[i] = alpha[i] - t * d[i];
+          C(y, newAlpha, lambda, w, ref coeff, ref r, ref U, ref s, ref VT, out rank, out nextLambda, out nextW, out newResNormSqr, out newCoeffNormSqr);
+          fx = resNormSqr + nextLambda * coeffNormSqr;
+          fx_new = newResNormSqr + nextLambda * newCoeffNormSqr;
+        }
+        report.lineSearchStep = t;
+
+        lambda = nextLambda;
+        w = nextW;
+        Array.Copy(newAlpha, alpha, alpha.Length);
+        resNormSqr = newResNormSqr;
+        coeffNormSqr = newCoeffNormSqr;
+        #endregion
+
+        // Step 7: If k > N, terminate the algorithm; else k = k + 1,
+        //         turn to step 2.
+        k++;
+
+        
+        iterationCallback?.Invoke(report, cancellationTokenSource.Token);
+
+        if (cancellationTokenSource.IsCancellationRequested) break;
+
+      }
+    }
+
+
+    // solves the regularized problem for the linear coefficients
+    // the error and penalty terms of the objective funcion C are returned separately
+    // additionally produces the updated WGCV parameters lambda and w
+    private void C(double[] y, double[] alpha,
+      double lambda, double w,
+      ref double[] coeff,
+      ref double[] r,
+      ref double[,] U,
+      ref double[] s,
+      ref double[,] VT,
+      out int rank,
+      out double nextLambda, out double nextW,
+      out double resNormSqr,
+      out double coeffNormSqr) {
+
+      Phi(alpha, ref F);
+      var n = F.GetLength(1);
+
+      var s_full = new double[n];
+      alglib.svd.rmatrixsvd(F, m, n, uneeded: 2, vtneeded: 2, additionalmemory: 2, w: ref s, u: ref U, vt: ref VT, null);
+
+      var tol = m * 2.2204460492503131E-16;  // the difference between 1.0 and the next larger double value
+      var s0 = s[0];
+      rank = s.Count(si => si > tol * s0);
+
+      // left-apply U to y
+      var uy = new double[m];
+      alglib.rmatrixgemv(m, m, alpha: 1.0, a: U, ia: 0, ja: 0, opa: 1,
+        x: y, ix: 0,
+        beta: 0.0, y: ref uy, iy: 0);
+
+      if (WGCV) {
+        // Step 2: Calculate the weight parameter w_k and
+        //         regularization parameter lambda_k using (24) and (23)
+        nextLambda = CalculateLambda(w, lambda, s, rank, uy);
+        nextW = CalculateW(nextLambda, s, rank, uy);
+      } else {
+        nextLambda = lambda;
+        nextW = w;
+      }
+
+      // Step 3: Compute the linear parameters theta_L using (16).
+      //         Obtain the residual vector r and its approximated
+      //         Jacobian Matrix J using J_GP, J_Kau, or J_R.
+
+      // eqn 16. 
+      // (using SVD)
+      // Golub & Van Loan, Matrix Computations, 4th Edition, Section 6.1.4 Ridge Regression, page 307.
+      if (coeff == null) coeff = new double[n];
+      double[] temp = new double[rank];
+      for (int i = 0; i < rank; i++) temp[i] = s[i] * uy[i] / (s[i] * s[i] + nextLambda);
+      alglib.rmatrixgemv(n, rank, 1.0, VT, 0, 0, 1, temp, 0, 0.0, ref coeff, 0);
+
+      // calculate residuals
+      r = (double[])y.Clone();
+      alglib.rmatrixgemv(m, n, -1.0, F, 0, 0, 0, coeff, 0, 1.0, ref r, 0); // r = y - F coeff
+      
+      resNormSqr = 0.0;
+      for (int i = 0; i < r.Length; i++) resNormSqr += r[i] * r[i];
+
+      coeffNormSqr = 0.0;
+      for (int i = 0; i < coeff.Length; i++) coeffNormSqr += coeff[i] * coeff[i];
+    }
+
+
+    // set full=false to drop the second term from the Jacobian (Kaufmann)
+    private void CalculateJacobian(JacobianFunc Jac, double[] alpha, double[] coeff, double[,] U, double[] s, double[,] VT, int rank, double[] resid, ref double[,] J, bool full = true) {
+      // code is based on VarPro.m by O'Leary and Rust
+      // https://www.cs.umd.edu/~oleary/software/varpro/
+
+      int numL = coeff.Length;
+      int numN = alpha.Length;
+
+      Jac(alpha, ref dF, ref ind);
+      var p = dF.GetLength(1);
+      double[] dF_r = new double[p]; // TODO allocate once
+      alglib.rmatrixgemv(p, m, 1.0, dF, 0, 0, 1, resid, 0, 0.0, ref dF_r, 0);
+      var T2 = new double[numL, numN]; // TODO allocate once
+
+      Array.Clear(J, 0, J.Length);
+      for (int j = 0; j < numN; j++) {                        // for each nonlinear parameter 
+        for (int col = 0; col < ind.GetLength(1); col++) {
+          if (ind[1, col] == j) {                             // columns of dPhi relevant to alpha(j)
+            var termIdx = ind[0, col];                        // relevant element of coeff
+            for (int k = 0; k < m; k++) {
+              J[k, j] += coeff[termIdx] * dF[k, col];
+            }
+            T2[termIdx, j] = dF_r[col];
+          }
+        }
+      }
+
+
+      double[,] tmp = new double[m - rank, numN]; 
+      alglib.rmatrixgemm(m - rank, numN, m, 1.0, U, 0, rank, 1, J, 0, 0, 0, 0.0, ref tmp, 0, 0);
+      alglib.rmatrixgemm(m, numN, m - rank, -1.0, U, 0, rank, 0, tmp, 0, 0, 0, 0.0, ref J, 0, 0);
+
+
+      double[,] tempT2 = new double[rank, numN];
+      alglib.rmatrixgemm(rank, numN, numL, 1.0, VT, 0, 0, 0, T2, 0, 0, 0, 0.0, ref tempT2, 0, 0);
+
+      for (int i = 0; i < rank; i++)
+        for (int j = 0; j < numN; j++)
+          tempT2[i, j] *= 1.0 / s[i];
+
+      alglib.rmatrixgemm(m, numN, rank, -1.0, U, 0, 0, 0, tempT2, 0, 0, 0, 0.0, ref Jac2, 0, 0);
+
+      if (full)
+        for (int i = 0; i < m; i++) {
+          for (int j = 0; j < numN; j++) {
+            J[i, j] += Jac2[i, j];
+          }
+        }
+    }
+
+    private static double[] SolveLS(double[,] A, double[] y) {
+      var m = A.GetLength(0);
+      var n = A.GetLength(1);
+      alglib.rmatrixsolvels(A, m, n, y, 0.0, out var info, out var report, out var x);
+      if (info < 0) throw new InvalidProgramException("Problem when solving linear system (LS)");
+      return x;
+    }
+
+    // Since we do not know lambda_opt we use the regularization parameter lambda_k obtained
+    // in the last iteration instead of lambda_opt. The initial value w_0 = 1.
+    // Solve Eqn (24). See Appendix
+    // length of s is limited to the rank
+    private static double CalculateW(double lambda_opt, double[] s, int rank, double[] uy) {
+      var lO = lambda_opt;
+      var m = uy.Length;
+      double h1(double lambda) {
+        var t1 = 0.0;
+        for (int i = 0; i < rank; i++) {
+          t1 += lambda * lambda * uy[i] * uy[i] / (s[i] * s[i] + lambda);
+        }
+        var t2 = 0.0;
+        for (int i = rank; i < m; i++) {
+          t2 += uy[i] * uy[i];
+        }
+        return t1 + t2;
+      }
+
+      // Eq 36
+      var s1 = 0.0; // numerator sum
+      var s2 = 0.0; // first term in denominator
+      var s3 = 0.0; // second term in denominator
+
+      var f3 = 0.0;
+      for (int i = 0; i < rank; i++) {
+        var div = s[i] * s[i] + lO;
+        s2 += s[i] * s[i] / (div * div);
+        f3 += s[i] * s[i] / div; // Tikhonov filter factor 
+      }
+
+      for (int i = 0; i < rank; i++) {
+        var suy = s[i] * uy[i];
+        var div = s[i] * s[i] + lO;
+        s1 += lO * suy * suy / (div * div * div);
+        s3 += s[i] * s[i] * lO * uy[i] * uy[i] / (div * div * div);
+      }
+
+      return m * s1 / (s2 * h1(lO) + f3 * s3);
+    }
+
+    // length of s is limited to rank
+    private static double CalculateLambda(double w, double lambda, double[] s, int rank, double[] uy) {
+      var m = uy.Length;
+      // eqn (23)
+      void G(double[] _x, ref double func, object _) {
+        var _lambda = _x[0];
+        var t1 = 0.0; // first term in numerator
+        var t2 = 0.0;
+        var dt1 = 0.0; // first term in denominator
+        var dt2 = 0.0;
+        double div;
+        for (int i = 0; i < rank; i++) {
+          div = 1.0 / (s[i] * s[i] + _lambda);
+          var t1_i = _lambda * uy[i] * div;
+          t1 += t1_i * t1_i;
+
+          dt1 += (1 - w) * s[i] * s[i] * div; // (1-w) varphi_i 
+          dt2 += _lambda * div;
+        }
+        for (int i = rank; i < m; i++) {
+          t2 += uy[i] * uy[i];
+        }
+
+        div = 1.0 / (dt1 + dt2 + m - rank);
+        func = rank * (t1 + t2) * div * div;
+      }
+
+      // When w is fixed the optimal lambda can be obtained by solving (23) easily
+      // using standard minimization algorithms.
+      double[] x = new double[] { lambda }; // initial value
+      alglib.minlbfgscreatef(1, x, 1e-6, out var state);
+      // alglib.minlbfgsoptguardgradient(state, 1e-5); // for debugging
+      alglib.minlbfgsoptimize(state, func: G, rep: null, obj: null);
+      // alglib.minlbfgsoptguardresults(state, out var optGuardReport);
+      // if (optGuardReport.badgradsuspected || optGuardReport.nonc0suspected || optGuardReport.nonc1suspected) throw new InvalidProgramException();
+      alglib.minlbfgsresults(state, out var lambda_opt, out var report);
+      if (report.terminationtype < 0) throw new InvalidProgramException();
+      return lambda_opt[0];
+    }
+  }
+}
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/alglib-license-gpl3.txt
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/alglib-license-gpl3.txt	(revision 18003)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HEAL.VarPro/alglib-license-gpl3.txt	(revision 18003)
@@ -0,0 +1,674 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+  The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works.  By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.  We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors.  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+  To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights.  Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received.  You must make sure that they, too, receive
+or can get the source code.  And you must show them these terms so they
+know their rights.
+
+  Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+  For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software.  For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+  Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so.  This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software.  The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable.  Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products.  If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+  Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary.  To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                       TERMS AND CONDITIONS
+
+  0. Definitions.
+
+  "This License" refers to version 3 of the GNU General Public License.
+
+  "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+  "The Program" refers to any copyrightable work licensed under this
+License.  Each licensee is addressed as "you".  "Licensees" and
+"recipients" may be individuals or organizations.
+
+  To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy.  The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+  A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+  To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy.  Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+  To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies.  Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+  An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License.  If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+  1. Source Code.
+
+  The "source code" for a work means the preferred form of the work
+for making modifications to it.  "Object code" means any non-source
+form of a work.
+
+  A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+  The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form.  A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+  The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities.  However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work.  For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+  The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+  The Corresponding Source for a work in source code form is that
+same work.
+
+  2. Basic Permissions.
+
+  All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met.  This License explicitly affirms your unlimited
+permission to run the unmodified Program.  The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work.  This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+  You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force.  You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright.  Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+  Conveying under any other circumstances is permitted solely under
+the conditions stated below.  Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+  No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+  When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+  4. Conveying Verbatim Copies.
+
+  You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+  You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+  5. Conveying Modified Source Versions.
+
+  You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+    a) The work must carry prominent notices stating that you modified
+    it, and giving a relevant date.
+
+    b) The work must carry prominent notices stating that it is
+    released under this License and any conditions added under section
+    7.  This requirement modifies the requirement in section 4 to
+    "keep intact all notices".
+
+    c) You must license the entire work, as a whole, under this
+    License to anyone who comes into possession of a copy.  This
+    License will therefore apply, along with any applicable section 7
+    additional terms, to the whole of the work, and all its parts,
+    regardless of how they are packaged.  This License gives no
+    permission to license the work in any other way, but it does not
+    invalidate such permission if you have separately received it.
+
+    d) If the work has interactive user interfaces, each must display
+    Appropriate Legal Notices; however, if the Program has interactive
+    interfaces that do not display Appropriate Legal Notices, your
+    work need not make them do so.
+
+  A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit.  Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+  6. Conveying Non-Source Forms.
+
+  You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+    a) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by the
+    Corresponding Source fixed on a durable physical medium
+    customarily used for software interchange.
+
+    b) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by a
+    written offer, valid for at least three years and valid for as
+    long as you offer spare parts or customer support for that product
+    model, to give anyone who possesses the object code either (1) a
+    copy of the Corresponding Source for all the software in the
+    product that is covered by this License, on a durable physical
+    medium customarily used for software interchange, for a price no
+    more than your reasonable cost of physically performing this
+    conveying of source, or (2) access to copy the
+    Corresponding Source from a network server at no charge.
+
+    c) Convey individual copies of the object code with a copy of the
+    written offer to provide the Corresponding Source.  This
+    alternative is allowed only occasionally and noncommercially, and
+    only if you received the object code with such an offer, in accord
+    with subsection 6b.
+
+    d) Convey the object code by offering access from a designated
+    place (gratis or for a charge), and offer equivalent access to the
+    Corresponding Source in the same way through the same place at no
+    further charge.  You need not require recipients to copy the
+    Corresponding Source along with the object code.  If the place to
+    copy the object code is a network server, the Corresponding Source
+    may be on a different server (operated by you or a third party)
+    that supports equivalent copying facilities, provided you maintain
+    clear directions next to the object code saying where to find the
+    Corresponding Source.  Regardless of what server hosts the
+    Corresponding Source, you remain obligated to ensure that it is
+    available for as long as needed to satisfy these requirements.
+
+    e) Convey the object code using peer-to-peer transmission, provided
+    you inform other peers where the object code and Corresponding
+    Source of the work are being offered to the general public at no
+    charge under subsection 6d.
+
+  A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+  A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling.  In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage.  For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product.  A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+  "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source.  The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+  If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information.  But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+  The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed.  Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+  Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+  7. Additional Terms.
+
+  "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law.  If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+  When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it.  (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.)  You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+  Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+    a) Disclaiming warranty or limiting liability differently from the
+    terms of sections 15 and 16 of this License; or
+
+    b) Requiring preservation of specified reasonable legal notices or
+    author attributions in that material or in the Appropriate Legal
+    Notices displayed by works containing it; or
+
+    c) Prohibiting misrepresentation of the origin of that material, or
+    requiring that modified versions of such material be marked in
+    reasonable ways as different from the original version; or
+
+    d) Limiting the use for publicity purposes of names of licensors or
+    authors of the material; or
+
+    e) Declining to grant rights under trademark law for use of some
+    trade names, trademarks, or service marks; or
+
+    f) Requiring indemnification of licensors and authors of that
+    material by anyone who conveys the material (or modified versions of
+    it) with contractual assumptions of liability to the recipient, for
+    any liability that these contractual assumptions directly impose on
+    those licensors and authors.
+
+  All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10.  If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term.  If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+  If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+  Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+  8. Termination.
+
+  You may not propagate or modify a covered work except as expressly
+provided under this License.  Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+  However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+  Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+  Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License.  If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+  9. Acceptance Not Required for Having Copies.
+
+  You are not required to accept this License in order to receive or
+run a copy of the Program.  Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance.  However,
+nothing other than this License grants you permission to propagate or
+modify any covered work.  These actions infringe copyright if you do
+not accept this License.  Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+  10. Automatic Licensing of Downstream Recipients.
+
+  Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License.  You are not responsible
+for enforcing compliance by third parties with this License.
+
+  An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations.  If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+  You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License.  For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+  11. Patents.
+
+  A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based.  The
+work thus licensed is called the contributor's "contributor version".
+
+  A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version.  For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+  Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+  In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement).  To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+  If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients.  "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+  If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+  A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License.  You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+  Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+  12. No Surrender of Others' Freedom.
+
+  If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all.  For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+  13. Use with the GNU Affero General Public License.
+
+  Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work.  The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+  14. Revised Versions of this License.
+
+  The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+  Each version is given a distinguishing version number.  If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation.  If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+  If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+  Later license versions may give you additional or different
+permissions.  However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+  15. Disclaimer of Warranty.
+
+  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. Limitation of Liability.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+  17. Interpretation of Sections 15 and 16.
+
+  If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+  If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+    <program>  Copyright (C) <year>  <name of author>
+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+  You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<http://www.gnu.org/licenses/>.
+
+  The GNU General Public License does not permit incorporating your program
+into proprietary programs.  If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.  But first, please read
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HeuristicLab.Problems.VarProMRGP-3.3.csproj
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HeuristicLab.Problems.VarProMRGP-3.3.csproj	(revision 17998)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/HeuristicLab.Problems.VarProMRGP-3.3.csproj	(revision 18003)
@@ -74,4 +74,8 @@
   </PropertyGroup>
   <ItemGroup>
+    <Reference Include="ALGLIB-3.17.0, Version=3.17.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\bin\ALGLIB-3.17.0.dll</HintPath>
+    </Reference>
     <Reference Include="HeuristicLab.Problems.DataAnalysis.Symbolic.NativeInterpreter-0.1, Version=0.0.0.1, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
       <SpecificVersion>False</SpecificVersion>
@@ -89,4 +93,5 @@
   </ItemGroup>
   <ItemGroup>
+    <Compile Include="HEAL.VarPro\VariableProjection.cs" />
     <Compile Include="VarProGrammar.cs" />
     <Compile Include="Plugin.cs" />
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Plugin.cs
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Plugin.cs	(revision 17998)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Plugin.cs	(revision 18003)
@@ -23,5 +23,5 @@
 
 namespace HeuristicLab.Problems.VarProMRGP {
-  [Plugin("HeuristicLab.Problems.VarProMRGP", "3.3.16.0")]
+  [Plugin("HeuristicLab.Problems.VarProMRGP", "3.3.16.17998")]
   [PluginFile("HeuristicLab.Problems.VarProMRGP-3.3.dll", PluginFileType.Assembly)]
   [PluginDependency("HeuristicLab.Algorithms.DataAnalysis", "3.4")]
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Problem.cs
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Problem.cs	(revision 17998)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Problem.cs	(revision 18003)
@@ -38,5 +38,4 @@
 using System.Runtime.InteropServices;
 using HeuristicLab.Random;
-using HeuristicLab.Analysis.Statistics;
 
 namespace HeuristicLab.Problems.VarProMRGP {
@@ -91,9 +90,8 @@
 
     #region not cloned or stored
-    private int[] featureIdx = null; // code[featureIdx[i]] is the last operation of feature i
-    [ExcludeFromObjectGraphTraversal]
-    private NativeInterpreter.NativeInstruction[] code = null;
-    private Dictionary<string, GCHandle> cachedData = null;
-    private LinearTransformation[] transformations;
+    ISymbolicExpressionTree[] features;
+    private List<TreeToAutoDiffTermConverter.ParametricFunctionGradient> featCode; // AutoDiff delegates for the features
+    private List<double[]> featParam; // parameters for the features
+    private List<double[][]> featVariables;
     #endregion
 
@@ -108,4 +106,6 @@
     public Problem() {
       var g = new VarProGrammar();
+
+      // TODO optionally: scale dataset
 
       Parameters.Add(new ValueParameter<IRegressionProblemData>("ProblemData", "", new RegressionProblemData()));
@@ -140,19 +140,17 @@
     // ProblemData
     //   |
-    //   +------------------+
-    //   |                  |
-    // cachedData         Grammar             MaxSize           MaxDepth          MaxInteractions
-    // cachedDataSet        |                    |                 |                 |
-    //                      +--------------------+-----------------+-----------------+
-    //                      |
-    //                     Features
-    //                      Code
-    //                      |
-    //                     Encoding (Length)
-    //                      |
-    //                      +--------------------+
-    //                      |                    |
-    //                   BestKnownSolution      Operators (Parameter)
-    //                   BestKnownQuality
+    // Grammar             MaxSize           MaxDepth          MaxInteractions
+    //   |                    |                 |                 |
+    //   +--------------------+-----------------+-----------------+
+    //   |
+    //  Features
+    //   Code
+    //   |
+    //  Encoding (Length)
+    //   |
+    //   +--------------------+
+    //   |                    |
+    // BestKnownSolution      Operators (Parameter)
+    // BestKnownQuality
 
     private void RegisterEventHandlers() {
@@ -187,5 +185,4 @@
 
     private void RegressionProblemData_Changed(object sender, EventArgs e) {
-      UpdateDataCache();
       Grammar.ConfigureVariableSymbols(RegressionProblemData);
     }
@@ -193,5 +190,4 @@
     private void RegressionProblemDataParameter_ValueChanged(object sender, EventArgs e) {
       RegressionProblemData.Changed += RegressionProblemData_Changed;
-      UpdateDataCache();
       Grammar.ConfigureVariableSymbols(RegressionProblemData);
     }
@@ -209,55 +205,7 @@
     }
 
-    private void UpdateDataCache() {
-      if (cachedData != null) {
-        foreach (var gch in cachedData.Values) {
-          gch.Free();
-        }
-        cachedData = null;
-      }
-
-      var dataset = RegressionProblemData.Dataset;
-
-      // free handles to old data
-      if (cachedData != null) {
-        foreach (var gch in cachedData.Values) {
-          gch.Free();
-        }
-        cachedData = null;
-      }
-
-      // cache new data
-      cachedData = new Dictionary<string, GCHandle>();
-      transformations = new LinearTransformation[dataset.DoubleVariables.Count()];
-      int varIdx = 0;
-      foreach (var v in dataset.DoubleVariables) {
-        var values = dataset.GetDoubleValues(v).ToArray();
-        if (v == RegressionProblemData.TargetVariable) {
-          // do not scale target
-          var linTransform = new LinearTransformation(dataset.DoubleVariables);
-          linTransform.Addend = 0;
-          linTransform.Multiplier = 1.0;
-          transformations[varIdx++] = linTransform;
-        } else {
-          // Scale to 0 .. 1
-          var max = values.Max();
-          var min = values.Min();
-          var range = max - min;
-          var linTransform = new LinearTransformation(dataset.DoubleVariables);
-          transformations[varIdx++] = linTransform;
-          linTransform.Addend = -1 * min / range - 0.0;
-          linTransform.Multiplier = 1.0 / range;
-          for (int i = 0; i < values.Length; i++) {
-            values[i] = values[i] * linTransform.Multiplier + linTransform.Addend;
-          }
-        }
-        var gch = GCHandle.Alloc(values, GCHandleType.Pinned);
-        cachedData[v] = gch;
-      }
-    }
-
     private void UpdateFeaturesAndCode() {
-      var features = GenerateFeaturesSystematic(10000, new MersenneTwister(31415), Grammar, MaxSize, MaxDepth, maxVariables: 3);
-      GenerateCode(features, RegressionProblemData.Dataset);
+      features = GenerateFeaturesSystematic(10000, new MersenneTwister(31415), Grammar, MaxSize, MaxDepth, maxVariables: 3);
+      GenerateCode(features, RegressionProblemData);
       var formatter = new InfixExpressionFormatter();
       Features = new ItemArray<StringValue>(features.Select(fi => new StringValue(formatter.Format(fi, System.Globalization.NumberFormatInfo.InvariantInfo, formatString: "0.0")).AsReadOnly())).AsReadOnly();
@@ -268,23 +216,8 @@
 
     public override double Evaluate(Individual individual, IRandom random) {
-      if (cachedData == null || code == null) {
-        // after loading from file or cloning the data cache and code must be restored
-        UpdateDataCache();
+      if (featCode == null) {
         UpdateFeaturesAndCode();
       }
       var b = individual.BinaryVector(Encoding.Name);
-      var fittingOptions = new NativeInterpreter.SolverOptions();
-      fittingOptions.Iterations = 10;
-      fittingOptions.Algorithm = NativeInterpreter.Algorithm.Krogh;
-      fittingOptions.LinearSolver = NativeInterpreter.CeresTypes.LinearSolverType.DENSE_QR;
-      fittingOptions.Minimizer = NativeInterpreter.CeresTypes.MinimizerType.TRUST_REGION;
-      fittingOptions.TrustRegionStrategy = NativeInterpreter.CeresTypes.TrustRegionStrategyType.LEVENBERG_MARQUARDT;
-
-      var evalOptions = new NativeInterpreter.SolverOptions();
-      evalOptions.Iterations = 0;
-      evalOptions.Algorithm = NativeInterpreter.Algorithm.Krogh;
-      evalOptions.LinearSolver = NativeInterpreter.CeresTypes.LinearSolverType.DENSE_QR;
-      evalOptions.Minimizer = NativeInterpreter.CeresTypes.MinimizerType.TRUST_REGION;
-      evalOptions.TrustRegionStrategy = NativeInterpreter.CeresTypes.TrustRegionStrategyType.LEVENBERG_MARQUARDT;
 
       var rows = RegressionProblemData.TrainingIndices.ToArray();
@@ -295,107 +228,96 @@
       for (int i = 0; i < b.Length; i++) {
         if (b[i] == true) {
-          termIndexList.Add(featureIdx[i]);
-        }
-      }
-
-      var coefficients = new double[termIndexList.Count + 1]; // coefficients for all terms + offset 
-                                                              // var avgCoefficients = new double[Features.Length];
-
-      // 5-fold CV
-      var nFolds = 5;
-      var nFoldRows = nRows / nFolds;
-      var fittingRows = new int[(nFolds - 1) * nFoldRows];
-      var fittingResult = new double[fittingRows.Length];
-
-      var testRows = new int[nFoldRows];
-      var testResult = new double[nFoldRows];
-      var cvMSE = new List<double>();
-      var trainMSE = new List<double>();
-      for (int testFold = 0; testFold < nFolds; testFold++) {
-        // TODO: THREAD SAFETY FOR PARALLEL EVALUATION
-        // CODE IS MANIPULATED
-        lock (code) {
-          var oldParameterValues = ExtractParameters(code, termIndexList);
-
-          #region fit to training folds
-          // copy rows from the four training folds
-          var nextFoldIdx = 0;
-          for (int fold = 0; fold < nFolds; fold++) {
-            if (fold == testFold) {
-              Array.Copy(allRows, fold * nFoldRows, testRows, 0, nFoldRows);
-            } else {
-              Array.Copy(allRows, fold * nFoldRows, fittingRows, nextFoldIdx, nFoldRows);
-              nextFoldIdx += nFoldRows;
+          termIndexList.Add(i);
+        }
+      }
+
+      var oldParameterValues = ExtractParameters(termIndexList);
+      var alpha = (double[])oldParameterValues.Clone();
+
+      var target = RegressionProblemData.TargetVariableTrainingValues.ToArray();
+
+      // local function for feature evaluation
+      void Phi(double[] a, ref double[,] phi) {
+        if (phi == null) {
+          phi = new double[nRows, termIndexList.Count + 1]; // + offset term
+          // last term is constant = 1
+          for (int i = 0; i < nRows; i++)
+            phi[i, termIndexList.Count] = 1.0;
+        }
+        var offset = 0;
+        // for each term
+        for (int i = 0; i < termIndexList.Count; i++) {
+          var termIdx = termIndexList[i];
+          var numFeatParam = this.featParam[termIdx].Length;
+          var variableValues = new double[featVariables[termIdx].Length];
+          var featParam = new double[numFeatParam];
+          Array.Copy(a, offset, featParam, 0, featParam.Length);
+          // for each row
+          for (int j = 0; j < nRows; j++) {
+            // copy row values
+            for (int k = 0; k < variableValues.Length; k++) {
+              variableValues[k] = featVariables[termIdx][k][j]; // featVariables is column-order
             }
-          }
-
-          var fittingTarget = RegressionProblemData.Dataset.GetDoubleValues(RegressionProblemData.TargetVariable, fittingRows).ToArray();
-
-          Array.Clear(fittingResult, 0, fittingResult.Length);
-          Array.Clear(coefficients, 0, coefficients.Length);
-
-
-          NativeInterpreter.NativeWrapper.GetValuesVarPro(code, code.Length, termIndexList.ToArray(), termIndexList.Count, fittingRows, fittingRows.Length,
-            coefficients, fittingOptions, fittingResult, fittingTarget, out var optSummary);
-          #endregion
-
-          if (optSummary.InitialCost < 0 || optSummary.FinalCost < 0) throw new InvalidProgramException();
-          trainMSE.Add(optSummary.FinalCost * 2.0 / fittingRows.Length); // ceres cost is 0.5 * sum of squared residuals
-                                                                         // average of all coefficients
-                                                                         // TODO return a statistic for the relevance of the term instead of the coefficient
-                                                                         // for (int k = 0; k < termIndexList.Count; k++) {
-                                                                         //   avgCoefficients[Array.IndexOf(featureIdx, termIndexList[k])] += coefficients[k] / nFolds; // TODO perf
-                                                                         // }
-
-          #region calculate output on test fold to determine CV MSE
-
-          // evaluate all terms for test rows
-          var sumOutput = new double[testResult.Length];
-          for (int r = 0; r < sumOutput.Length; r++) {
-            sumOutput[r] = coefficients[coefficients.Length - 1]; // offset
-          }
-
-          for (int k = 0; k < termIndexList.Count; k++) {
-            var termIdx = termIndexList[k];
-            // copy relevant part of the code
-            var termCode = new NativeInterpreter.NativeInstruction[code[termIdx].length];
-            Array.Copy(code, termIdx - termCode.Length + 1, termCode, 0, termCode.Length);
-            NativeInterpreter.NativeWrapper.GetValues(termCode, termCode.Length, testRows, testRows.Length, evalOptions, testResult, null, out var evalSummary);
-
-            for (int r = 0; r < sumOutput.Length; r++) {
-              sumOutput[r] += coefficients[k] * testResult[r];
+            var tup = featCode[termIdx].Invoke(featParam, variableValues); // TODO for phi we do not actually need g
+            phi[j, i] = tup.Item2;
+          }
+          offset += numFeatParam;
+        }
+      }
+
+      // local function for Jacobian evaluation
+      void Jac(double[] a, ref double[,] J, ref int[,] ind) {
+        if (J == null) {
+          J = new double[nRows, featParam.Sum(fp => fp.Length)]; // all parameters
+          ind = new int[2, featParam.Sum(fp => fp.Length)];
+        }
+        var offset = 0;
+        // for each term
+        for (int i = 0; i < termIndexList.Count; i++) {
+          var termIdx = termIndexList[i];
+          var numFeatParam = this.featParam[termIdx].Length;
+          var variableValues = new double[featVariables[termIdx].Length];
+          var featParam = new double[numFeatParam];
+          Array.Copy(a, offset, featParam, 0, featParam.Length);
+
+          // for each parameter
+          for (int k = 0; k < featParam.Length; k++) {
+            ind[0, offset + k] = i; // column idx in phi
+            ind[1, offset + k] = offset + k; // parameter idx (no parameter is used twice)
+          }
+
+          // for each row
+          for (int j = 0; j < nRows; j++) {
+            // copy row values
+            for (int k = 0; k < variableValues.Length; k++) {
+              variableValues[k] = featVariables[termIdx][k][j]; // featVariables is column-order
             }
-          }
-
-          var evalTarget = RegressionProblemData.Dataset.GetDoubleValues(RegressionProblemData.TargetVariable, testRows);
-          if (sumOutput.Any(d => double.IsNaN(d))) cvMSE.Add(evalTarget.VariancePop()); // variance of target is MSE of constant model
-          else cvMSE.Add(sumOutput.Zip(evalTarget, (ri, ti) => { var res = ti - ri; return res * res; }).Average());
-
-          #endregion
-
-          #region prepare for next varpro call
-          var newParameterValues = ExtractParameters(code, termIndexList);
-          // keep the updated parameter values only when the coefficient is significantly <> 0
-          for (int i = 0; i < termIndexList.Count; i++) {
-            // TODO would actually need to take variance of term into account
-            if (Math.Abs(coefficients[i]) > 1e-5) {
-              oldParameterValues[i] = newParameterValues[i];
+            var tup = featCode[termIdx].Invoke(featParam, variableValues);
+            // for each parameter
+            for (int k = 0; k < featParam.Length; k++) {
+              J[j, offset + k] = tup.Item1[k];
             }
           }
-          UpdateParameter(code, termIndexList, oldParameterValues);
-        }
-        #endregion
-      }
-
-
-      individual["Train MSE (avg)"] = new DoubleValue(trainMSE.Average());
-      individual["Parameter"] = new DoubleArray(ExtractParameters(code, termIndexList).SelectMany(arr => arr).ToArray()); // store the parameter which worked for this individual for solution creation
-                                                                                                                          // individual["Coefficients"] = new DoubleArray(avgCoefficients); // for allele frequency analysis
-                                                                                                                          // return new double[] { avgCost, coefficients.Count(ci => Math.Abs(ci) > 1e-9) };
-                                                                                                                          // average of averages 
-      individual["CV MSE (avg)"] = new DoubleValue(cvMSE.Average());
-      individual["CV MSE (stdev)"] = new DoubleValue(cvMSE.StandardDeviation());
-      return cvMSE.Average() + cvMSE.StandardDeviation();
-      // return new double[] { cvMSE.Average() + cvMSE.StandardDeviation(), termIndexList.Count };
+          offset += numFeatParam;
+        }
+      }
+
+      try {
+        HEAL.VarPro.VariableProjection.Fit(Phi, Jac, target, alpha, out var coeff, out var report);
+
+
+        if (report.residNorm < 0) throw new InvalidProgramException();
+        UpdateParameter(termIndexList, alpha);
+
+
+        individual["Parameter"] = new DoubleArray(alpha); // store the parameter which worked for this individual for solution creation
+        individual["Coeff"] = new DoubleArray(coeff);
+
+        return report.residNormSqr / nRows;
+      } catch (Exception _) {
+        individual["Parameter"] = new DoubleArray(alpha); // store the parameter which worked for this individual for solution creation
+        individual["Coeff"] = new DoubleArray(termIndexList.Count + 1);
+        return double.MaxValue;
+      }
     }
 
@@ -415,105 +337,48 @@
       if (double.IsNaN(curBestQuality) || bestQ < curBestQuality) {
         var bestVector = bestIndividual.BinaryVector(Encoding.Name);
-
-        var options = new NativeInterpreter.SolverOptions();
-        options.Iterations = 10; // optimize on full dataset
-        options.Algorithm = NativeInterpreter.Algorithm.Krogh;
-        options.LinearSolver = NativeInterpreter.CeresTypes.LinearSolverType.DENSE_QR;
-        options.Minimizer = NativeInterpreter.CeresTypes.MinimizerType.TRUST_REGION;
-        options.TrustRegionStrategy = NativeInterpreter.CeresTypes.TrustRegionStrategyType.LEVENBERG_MARQUARDT;
-
-        var rows = RegressionProblemData.TrainingIndices.ToArray();
-        var target = RegressionProblemData.TargetVariableTrainingValues.ToArray();
-
-        var rowsArray = rows.ToArray();
-        var nRows = rowsArray.Length;
-        var result = new double[nRows];
+        var bestParams = ((DoubleArray)bestIndividual["Parameter"]).ToArray();
+        var bestCoeff = ((DoubleArray)bestIndividual["Coeff"]).ToArray();
+        // var rows = RegressionProblemData.TrainingIndices.ToArray();
+        // var target = RegressionProblemData.TargetVariableTrainingValues.ToArray();
+        // 
+        // var rowsArray = rows.ToArray();
+        // var nRows = rowsArray.Length;
+        // var result = new double[nRows];
         var termIndexList = new List<int>();
-        var predictorNames = new List<string>();
+        // var predictorNames = new List<string>();
         for (int i = 0; i < bestVector.Length; i++) {
           if (bestVector[i] == true) {
-            termIndexList.Add(featureIdx[i]);
-            predictorNames.Add(Features[i].Value);
-          }
-        }
-        var coefficients = new double[termIndexList.Count + 1]; // coefficients for all terms + offset 
-                                                                // TODO: thread-safety
-        lock (code) {
-          UpdateParameter(code, termIndexList, ((DoubleArray)bestIndividual["Parameter"]).ToArray());
-
-          NativeInterpreter.NativeWrapper.GetValuesVarPro(code, code.Length, termIndexList.ToArray(), termIndexList.Count, rows, nRows, coefficients, options, result, target, out var optSummary);
-        }
-
-        #region linear model statistics
-        var xy = new double[nRows, termIndexList.Count + 1];
-        for(int j=0;j<termIndexList.Count;j++) {
-          var t = termIndexList[j];
-          var termCode = new NativeInterpreter.NativeInstruction[code[t].length];
-          Array.Copy(code, t - code[t].length + 1, termCode, 0, code[t].length);
-          options.Iterations = 0;
-          NativeInterpreter.NativeWrapper.GetValues(termCode, termCode.Length, rows, nRows, options, result, null, out _);
-          for(int i=0;i<nRows;i++) { xy[i, j] = result[i]; } // copy to matrix
-        }
-        // copy target to matrix
-        for (int i = 0; i < nRows; i++) { xy[i, termIndexList.Count] = target[i]; } // copy to matrix
-        predictorNames.Add("<const>"); // last coefficient is offset
-        // var stats = Statistics.CalculateLinearModelStatistics(xy, coefficients);
-        // results.AddOrUpdateResult("Statistics", stats.AsResultCollection(predictorNames));
-        #endregion
-
-        results.AddOrUpdateResult("Solution", CreateRegressionSolution(code, termIndexList, coefficients, RegressionProblemData));
+            termIndexList.Add(i);
+          }
+        }
+
+        results.AddOrUpdateResult("Solution", CreateRegressionSolution(termIndexList.ToArray(), bestParams, bestCoeff, RegressionProblemData));
         results.AddOrUpdateResult("BestQuality", new DoubleValue(bestQ));
-        results.AddOrUpdateResult("CV MSE (avg)", bestIndividual["CV MSE (avg)"]);
-        results.AddOrUpdateResult("CV MSE (stdev)", bestIndividual["CV MSE (stdev)"]);
       }
     }
 
     #region retrieval / update of non-linear parameters
-    private double[][] ExtractParameters(NativeInterpreter.NativeInstruction[] code, List<int> termIndexList) {
-      var p = new double[termIndexList.Count][];
+    private double[] ExtractParameters(List<int> termIndexList) {
+      var p = new List<double>();
       for (int i = 0; i < termIndexList.Count; i++) {
-        var termIdx = termIndexList[i];
-        var pList = new List<double>();
-        var start = termIdx - code[termIdx].length + 1;
-        for (int codeIdx = start; codeIdx <= termIdx; codeIdx++) {
-          if (code[codeIdx].optimize)
-            pList.Add(code[codeIdx].value);
-        }
-        p[i] = pList.ToArray();
-      }
-      return p;
-    }
-
-    // parameters are given as array for each term
-    private void UpdateParameter(NativeInterpreter.NativeInstruction[] code, List<int> termIndexList, double[][] p) {
+        p.AddRange(featParam[termIndexList[i]]);
+      }
+      return p.ToArray();
+    }
+
+
+    // parameters are given as a flat array
+    private void UpdateParameter(List<int> termIndexList, double[] p) {
+      var offset = 0;
       for (int i = 0; i < termIndexList.Count; i++) {
-        if (p[i].Length == 0) continue; // nothing to update
-        var termIdx = termIndexList[i];
-        var pIdx = 0;
-        var start = termIdx - code[termIdx].length + 1;
-        for (int codeIdx = start; codeIdx <= termIdx; codeIdx++) {
-          if (code[codeIdx].optimize)
-            code[codeIdx].value = p[i][pIdx++];
-        }
-      }
-    }
-
-    // parameters are given as a flat array
-    private void UpdateParameter(NativeInterpreter.NativeInstruction[] code, List<int> termIndexList, double[] p) {
-      var pIdx = 0;
-      for (int i = 0; i < termIndexList.Count; i++) {
-        var termIdx = termIndexList[i];
-        var start = termIdx - code[termIdx].length + 1;
-        for (int codeIdx = start; codeIdx <= termIdx; codeIdx++) {
-          if (code[codeIdx].optimize)
-            code[codeIdx].value = p[pIdx++];
-        }
+        var numFeatParam = featParam[termIndexList[i]].Length;
+        Array.Copy(p, offset, featParam[termIndexList[i]], 0, numFeatParam);
+        offset += numFeatParam;
       }
     }
     #endregion
 
-
-
     #region feature generation
+    /*
     private static ISymbolicExpressionTree[] GenerateFeatures(int n, IRandom random, ISymbolicDataAnalysisGrammar grammar, int maxSize, int maxDepth) {
       var features = new ISymbolicExpressionTree[n];
@@ -531,4 +396,5 @@
       return features;
     }
+    */
     private static ISymbolicExpressionTree[] GenerateFeaturesSystematic(int n, IRandom random, ISymbolicDataAnalysisGrammar grammar, int maxSize, int maxDepth, int maxVariables) {
       var hashes = new HashSet<ulong>();
@@ -607,22 +473,15 @@
 
 
-    private void GenerateCode(ISymbolicExpressionTree[] features, IDataset dataset) {
-      byte mapSupportedSymbols(ISymbolicExpressionTreeNode node) {
-        var opCode = OpCodes.MapSymbolToOpCode(node);
-        if (supportedOpCodes.Contains(opCode)) return opCode;
-        else throw new NotSupportedException($"VarPro does not support {node.Symbol.Name}");
-      };
-
-      var i = 0;
-      featureIdx = new int[features.Length];
-
-      var code = new List<NativeInterpreter.NativeInstruction>(capacity: (int)(features.Sum(fi => fi.Length - 2) * 1.2));
+    private void GenerateCode(ISymbolicExpressionTree[] features, IRegressionProblemData problemData) {
+      this.featCode = new List<TreeToAutoDiffTermConverter.ParametricFunctionGradient>();
+      this.featParam = new List<double[]>();
+      this.featVariables = new List<double[][]>();
       foreach (var f in features) {
-        var featureCode = Compile(f, mapSupportedSymbols);
-
-        code.AddRange(featureCode);
-        featureIdx[i++] = code.Count - 1;
-      }
-      this.code = code.ToArray();
+        var featureCode = Compile(f, problemData, out var initialParamValues, out var variableValues);
+
+        featCode.Add(featureCode);
+        featParam.Add(initialParamValues);
+        featVariables.Add(variableValues);
+      }
     }
 
@@ -650,75 +509,44 @@
       (byte)OpCode.AnalyticQuotient
     };
-    private NativeInterpreter.NativeInstruction[] Compile(ISymbolicExpressionTree tree, Func<ISymbolicExpressionTreeNode, byte> opCodeMapper) {
-      var root = tree.Root.GetSubtree(0).GetSubtree(0);
-      var code = new List<NativeInterpreter.NativeInstruction>();
-      foreach (var n in root.IterateNodesPrefix()) {
-        var instr = new NativeInterpreter.NativeInstruction { narg = (ushort)n.SubtreeCount, opcode = opCodeMapper(n), length = 1 }; // length is updated in a second pass below
-        if (n is VariableTreeNode variable) {
-          instr.value = variable.Weight;
-          instr.optimize = false;
-          instr.data = cachedData[variable.VariableName].AddrOfPinnedObject();
-        } else if (n is ConstantTreeNode constant) {
-          instr.value = constant.Value;
-          if (n.Symbol.Name != "<1.0>") // HACK TODO this depends on the name given in the grammar!
-            instr.optimize = true; // the VarPro grammar is specifically designed to make sure we have all necessary and only necessary non-linear parameters
-        }
-        if (n.Symbol is Logarithm) {
-          // for log(f(x)) generate code log ( sqrt(f(x)²)) to ensure argument to log is positive
-          code.Add(instr);
-          code.Add(new NativeInterpreter.NativeInstruction { narg = 1, opcode = (byte)OpCode.SquareRoot, length = 1 }); // length is updated in a second pass below
-          code.Add(new NativeInterpreter.NativeInstruction { narg = 1, opcode = (byte)OpCode.Square, length = 1 });
-        } else {
-          // all other operations are added verbatim
-          code.Add(instr);
-        }
-      }
-
-      code.Reverse();
-      var codeArr = code.ToArray();
-
-      // second pass to calculate lengths
-      for (int i = 0; i < codeArr.Length; i++) {
-        var c = i - 1;
-        for (int j = 0; j < codeArr[i].narg; ++j) {
-          codeArr[i].length += codeArr[c].length;
-          c -= codeArr[c].length;
-        }
-      }
-
-      return codeArr;
+    private TreeToAutoDiffTermConverter.ParametricFunctionGradient Compile(ISymbolicExpressionTree tree, IRegressionProblemData problemData,
+      out double[] initialParameterValues, out double[][] variableValues) {
+      TreeToAutoDiffTermConverter.TryConvertToAutoDiff(tree, makeVariableWeightsVariable: false, addLinearScalingTerms: false,
+        out var parameters, out initialParameterValues, out var func, out var func_grad);
+      variableValues = new double[parameters.Count][];
+      for (int i = 0; i < parameters.Count; i++) {
+        variableValues[i] = problemData.Dataset.GetDoubleValues(parameters[i].variableName, problemData.TrainingIndices).ToArray(); // TODO: we could reuse the arrays
+      }
+      return func_grad;
     }
     #endregion
 
     #region solution creation
-    private IRegressionSolution CreateRegressionSolution(NativeInterpreter.NativeInstruction[] code, IList<int> termIndexList, double[] coefficients, IRegressionProblemData problemData) {
-      // parse back solution from code (required because we need to used optimized parameter values)
-      var addSy = symbols[(byte)OpCode.Add];
-      var mulSy = symbols[(byte)OpCode.Mul];
-      var sum = addSy.CreateTreeNode();
-      for (int i = 0; i < termIndexList.Count; i++) {
-        if (Math.Abs(coefficients[i]) < 1e-8) continue;
-        var termIdx = termIndexList[i];
-        var prod = mulSy.CreateTreeNode();
-        var constNode = (ConstantTreeNode)constSy.CreateTreeNode();
-        constNode.Value = coefficients[i];
-        var term = CreateTree(code, termIdx);
-        prod.AddSubtree(constNode);
-        prod.AddSubtree(term);
-        sum.AddSubtree(prod);
-      }
-      {
-        var constNode = (ConstantTreeNode)constSy.CreateTreeNode();
-        constNode.Value = coefficients.Last();
-        sum.AddSubtree(constNode);
-      }
+    private IRegressionSolution CreateRegressionSolution(int[] featIdx, double[] parameters, double[] coefficients, IRegressionProblemData problemData) {
       var root = (new ProgramRootSymbol()).CreateTreeNode();
       var start = (new StartSymbol()).CreateTreeNode();
+      var add = (new Addition()).CreateTreeNode();
       root.AddSubtree(start);
-      start.AddSubtree(sum);
+      start.AddSubtree(add);
+      var offset = 0;
+      for (int i = 0; i < featIdx.Length; i++) {
+        var term = (ISymbolicExpressionTreeNode)features[featIdx[i]].Root.GetSubtree(0).GetSubtree(0).Clone();
+
+        var termParameters = new double[featParam[featIdx[i]].Length];
+        Array.Copy(parameters, offset, termParameters, 0, termParameters.Length);
+        ReplaceParameters(term, termParameters);
+        offset += termParameters.Length;
+
+        var mul = (new Multiplication()).CreateTreeNode();
+        mul.AddSubtree(term);
+        mul.AddSubtree(CreateConstant(coefficients[i]));
+        add.AddSubtree(mul);
+      }
+      // last coeff is offset
+      add.AddSubtree(CreateConstant(coefficients[coefficients.Length - 1]));
+
       var tree = new SymbolicExpressionTree(root);
       var ds = problemData.Dataset;
-      var scaledDataset = new Dataset(ds.DoubleVariables, ds.ToArray(ds.DoubleVariables, transformations, Enumerable.Range(0, ds.Rows)));
-      var scaledProblemData = new RegressionProblemData(scaledDataset, problemData.AllowedInputVariables, problemData.TargetVariable, transformations);
+      var scaledDataset = new Dataset(ds.DoubleVariables, ds.ToArray(ds.DoubleVariables, Enumerable.Range(0, ds.Rows)));
+      var scaledProblemData = new RegressionProblemData(scaledDataset, problemData.AllowedInputVariables, problemData.TargetVariable);
       scaledProblemData.TrainingPartition.Start = problemData.TrainingPartition.Start;
       scaledProblemData.TrainingPartition.End = problemData.TrainingPartition.End;
@@ -727,4 +555,21 @@
       return new SymbolicRegressionSolution(
         new SymbolicRegressionModel(problemData.TargetVariable, tree, new SymbolicDataAnalysisExpressionTreeNativeInterpreter()), scaledProblemData);
+    }
+
+    private void ReplaceParameters(ISymbolicExpressionTreeNode term, double[] termParameters) {
+      // Autodiff converter extracts parameters using a pre-order tree traversal.
+      // Therefore, we must use a pre-order tree traversal here as well.
+      // Only ConstantTreeNode values are optimized.
+      var paramIdx = 0;
+      foreach (var node in term.IterateNodesPrefix().OfType<ConstantTreeNode>()) {
+        node.Value = termParameters[paramIdx++];
+      }
+      if (paramIdx != termParameters.Length) throw new InvalidProgramException();
+    }
+
+    private ISymbolicExpressionTreeNode CreateConstant(double coeff) {
+      var constNode = (ConstantTreeNode)(new Constant()).CreateTreeNode();
+      constNode.Value = coeff;
+      return constNode;
     }
 
@@ -753,51 +598,4 @@
 
 
-    private ISymbolicExpressionTreeNode CreateTree(NativeInterpreter.NativeInstruction[] code, int i) {
-      switch (code[i].opcode) {
-        case (byte)OpCode.Add:
-        case (byte)OpCode.Sub:
-        case (byte)OpCode.Mul:
-        case (byte)OpCode.Div:
-        case (byte)OpCode.Exp:
-        case (byte)OpCode.Log:
-        case (byte)OpCode.Sin:
-        case (byte)OpCode.Cos:
-        case (byte)OpCode.Tan:
-        case (byte)OpCode.Tanh:
-        case (byte)OpCode.Square:
-        case (byte)OpCode.SquareRoot:
-        case (byte)OpCode.Cube:
-        case (byte)OpCode.CubeRoot:
-        case (byte)OpCode.Absolute:
-        case (byte)OpCode.AnalyticQuotient: {
-            var node = symbols[code[i].opcode].CreateTreeNode();
-            var c = i - 1;
-            for (int childIdx = 0; childIdx < code[i].narg; childIdx++) {
-              node.AddSubtree(CreateTree(code, c));
-              c = c - code[c].length;
-            }
-            return node;
-          }
-        case (byte)OpCode.Constant: {
-            var node = (ConstantTreeNode)constSy.CreateTreeNode();
-            node.Value = code[i].value;
-            return node;
-          }
-        case (byte)OpCode.Variable: {
-            var node = (VariableTreeNode)varSy.CreateTreeNode();
-            node.Weight = code[i].value;
-            // TODO perf
-            node.VariableName = string.Empty;
-            foreach (var tup in this.cachedData) {
-              if (tup.Value.AddrOfPinnedObject() == code[i].data) {
-                node.VariableName = tup.Key;
-                break;
-              }
-            }
-            return node;
-          }
-        default: throw new NotSupportedException("unknown opcode");
-      }
-    }
     #endregion
 
@@ -809,13 +607,13 @@
       Operators.Add(new AlleleFrequencyAnalyzer());
 
-      var cvMSEAnalyzer = new BestAverageWorstQualityAnalyzer();
-      cvMSEAnalyzer.Name = "CVMSE Analzer";
-      ParameterizeAnalyzer(cvMSEAnalyzer, "CV MSE (avg)");
-      Operators.Add(cvMSEAnalyzer);
-      
-      var trainingMSEAnalyzer = new BestAverageWorstQualityAnalyzer();
-      trainingMSEAnalyzer.Name = "Training MSE Analzer";
-      ParameterizeAnalyzer(trainingMSEAnalyzer, "Train MSE (avg)");
-      Operators.Add(trainingMSEAnalyzer);
+      // var cvMSEAnalyzer = new BestAverageWorstQualityAnalyzer();
+      // cvMSEAnalyzer.Name = "CVMSE Analzer";
+      // ParameterizeAnalyzer(cvMSEAnalyzer, "CV MSE (avg)");
+      // Operators.Add(cvMSEAnalyzer);
+      // 
+      // var trainingMSEAnalyzer = new BestAverageWorstQualityAnalyzer();
+      // trainingMSEAnalyzer.Name = "Training MSE Analzer";
+      // ParameterizeAnalyzer(trainingMSEAnalyzer, "Train MSE (avg)");
+      // Operators.Add(trainingMSEAnalyzer);
 
       ParameterizeOperators();
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Properties/AssemblyInfo.cs
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Properties/AssemblyInfo.cs	(revision 17998)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/Properties/AssemblyInfo.cs	(revision 18003)
@@ -53,3 +53,3 @@
 // [assembly: AssemblyVersion("1.0.*")]
 [assembly: AssemblyVersion("3.3.0.0")]
-[assembly: AssemblyFileVersion("3.3.16.0")]
+[assembly: AssemblyFileVersion("3.3.16.17998")]
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/VarProGrammar.cs
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/VarProGrammar.cs	(revision 17998)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Problems.VarProMRGP/3.3/VarProGrammar.cs	(revision 18003)
@@ -51,4 +51,5 @@
       Symbol aqSy = new AnalyticQuotient();
       Symbol logSy = new Logarithm();
+      Symbol absSy = new Absolute();
       Symbol expSy = new Exponential();
 
@@ -69,5 +70,5 @@
       varSy.WeightSigma = 0.0;
 
-      var allSymbols = new List<Symbol>() { sumForLogSy, mulSy, simpleMulSy, aqSy, logSy, expSy, constSy, constOneSy, varSy };
+      var allSymbols = new List<Symbol>() { sumForLogSy, mulSy, simpleMulSy, aqSy, logSy, absSy, expSy, constSy, constOneSy, varSy };
 
       foreach (var symb in allSymbols)
@@ -79,5 +80,5 @@
 
 
-      foreach (var funSymb in new[] { logSy, expSy }) {
+      foreach (var funSymb in new[] { logSy, expSy, absSy }) {
         SetSubtreeCount(funSymb, 1, 1);
       }
@@ -99,6 +100,9 @@
       }
 
+      // log(abs(sum))
       // allowed for log:
-      AddAllowedChildSymbol(logSy, sumForLogSy);
+      AddAllowedChildSymbol(logSy, absSy);
+      // allowed for abs:
+      AddAllowedChildSymbol(absSy, sumForLogSy);
 
       // allowed for exp:
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/VarProMRGPTest.cs
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/VarProMRGPTest.cs	(revision 18003)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Tests/HeuristicLab.Algorithms.DataAnalysis-3.4/VarProMRGPTest.cs	(revision 18003)
@@ -0,0 +1,53 @@
+﻿#region License Information
+/* HeuristicLab
+ * Copyright (C) Heuristic and Evolutionary Algorithms Laboratory (HEAL)
+ *
+ * This file is part of HeuristicLab.
+ *
+ * HeuristicLab is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * HeuristicLab is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
+ */
+#endregion
+
+
+using System.Linq;
+using Microsoft.VisualStudio.TestTools.UnitTesting;
+
+namespace HeuristicLab.Tests {
+  [TestClass]
+  public class VarProMRGPTest {
+
+    [TestMethod]
+    [TestCategory("Algorithms.DataAnalysis")]
+    [TestProperty("Time", "medium")]
+    public void VarProMRGPPoly10() {
+      var prov = new HeuristicLab.Problems.Instances.DataAnalysis.VariousInstanceProvider();
+      var poly10 = prov.GetDataDescriptors().First(dd => dd.Name.Contains("Poly"));
+      var problemData = prov.LoadData(poly10);
+      var prob = new HeuristicLab.Problems.VarProMRGP.Problem();
+      prob.RegressionProblemData = problemData;
+      var alg = new HeuristicLab.Algorithms.GeneticAlgorithm.GeneticAlgorithm();
+      alg.Problem = prob;
+
+      alg.MaximumGenerations.Value = 2;
+      alg.PopulationSize.Value = 10;
+      alg.MutationProbability.Value = 0.1;
+      alg.SetSeedRandomly.Value = false;
+      alg.Seed.Value = 1234;
+
+      alg.Start();
+
+
+    }
+  }
+}
Index: branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Tests/HeuristicLab.Tests.csproj
===================================================================
--- branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Tests/HeuristicLab.Tests.csproj	(revision 17998)
+++ branches/3127-MRGP-VarPro-Exploration/HeuristicLab.Tests/HeuristicLab.Tests.csproj	(revision 18003)
@@ -410,4 +410,8 @@
       <HintPath>..\bin\HeuristicLab.Problems.TravelingSalesman-3.3.dll</HintPath>
       <Private>False</Private>
+    </Reference>
+    <Reference Include="HeuristicLab.Problems.VarProMRGP-3.3, Version=3.3.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\bin\HeuristicLab.Problems.VarProMRGP-3.3.dll</HintPath>
     </Reference>
     <Reference Include="HeuristicLab.Problems.VehicleRouting-3.4">
@@ -499,4 +503,5 @@
     <Compile Include="HeuristicLab.Algorithms.DataAnalysis-3.4\CrossValidationTest.cs" />
     <Compile Include="HeuristicLab.Algorithms.DataAnalysis-3.4\NcaAlgorithmTest.cs" />
+    <Compile Include="HeuristicLab.Algorithms.DataAnalysis-3.4\VarProMRGPTest.cs" />
     <Compile Include="HeuristicLab.Algorithms.DataAnalysis-3.4\SupportVectorMachineTest.cs" />
     <Compile Include="HeuristicLab.Algorithms.DataAnalysis-3.4\GaussianProcessModelTest.cs" />
