Index: branches/HeuristicLab.Problems.GrammaticalOptimization-gkr/HeuristicLab.Problems.GrammaticalOptimization.SymbReg/SymbolicRegressionProblem.cs
===================================================================
--- branches/HeuristicLab.Problems.GrammaticalOptimization/HeuristicLab.Problems.GrammaticalOptimization.SymbReg/SymbolicRegressionProblem.cs	(revision 12099)
+++ branches/HeuristicLab.Problems.GrammaticalOptimization-gkr/HeuristicLab.Problems.GrammaticalOptimization.SymbReg/SymbolicRegressionProblem.cs	(revision 12290)
@@ -3,4 +3,6 @@
 using System.Collections.Generic;
 using System.Diagnostics;
+using System.Globalization;
+using System.IO;
 using System.Linq;
 using System.Security;
@@ -69,6 +71,10 @@
     private readonly bool useConstantOpt;
     public string Name { get; private set; }
-
-    public SymbolicRegressionProblem(Random random, string partOfName, bool useConstantOpt = true) {
+    private Random random;
+    private double lambda;
+
+
+    // lambda should be tuned using CV
+    public SymbolicRegressionProblem(Random random, string partOfName, double lambda = 1.0, bool useConstantOpt = true) {
       var instanceProviders = new RegressionInstanceProvider[] 
       {new RegressionRealWorldInstanceProvider(),
@@ -80,11 +86,25 @@
     };
       var instanceProvider = instanceProviders.FirstOrDefault(prov => prov.GetDataDescriptors().Any(dd => dd.Name.Contains(partOfName)));
-      if (instanceProvider == null) throw new ArgumentException("instance not found");
+      IRegressionProblemData problemData = null;
+      if (instanceProvider != null) {
+        var dds = instanceProvider.GetDataDescriptors();
+        problemData = instanceProvider.LoadData(dds.Single(ds => ds.Name.Contains(partOfName)));
+
+      } else if (File.Exists(partOfName)) {
+        // check if it is a file
+        var prov = new RegressionCSVInstanceProvider();
+        problemData = prov.ImportData(partOfName);
+        problemData.TrainingPartition.Start = 0;
+        problemData.TrainingPartition.End = problemData.Dataset.Rows;
+        // no test partition
+        problemData.TestPartition.Start = problemData.Dataset.Rows;
+        problemData.TestPartition.End = problemData.Dataset.Rows;
+      } else {
+        throw new ArgumentException("instance not found");
+      }
 
       this.useConstantOpt = useConstantOpt;
 
-      var dds = instanceProvider.GetDataDescriptors();
-      var problemData = instanceProvider.LoadData(dds.Single(ds => ds.Name.Contains(partOfName)));
-      this.Name = problemData.Name;
+      this.Name = problemData.Name + string.Format("lambda={0:N2}", lambda);
 
       this.N = problemData.TrainingIndices.Count();
@@ -94,12 +114,32 @@
       this.y = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();
 
+      var varEst = new OnlineMeanAndVarianceCalculator();
+
+      var means = new double[d];
+      var stdDevs = new double[d];
+
       int i = 0;
+      foreach (var inputVariable in problemData.AllowedInputVariables) {
+        varEst.Reset();
+        problemData.Dataset.GetDoubleValues(inputVariable).ToList().ForEach(varEst.Add);
+        if (varEst.VarianceErrorState != OnlineCalculatorError.None) throw new ArgumentException();
+        means[i] = varEst.Mean;
+        stdDevs[i] = Math.Sqrt(varEst.Variance);
+        i++;
+      }
+
+      i = 0;
       foreach (var r in problemData.TrainingIndices) {
         int j = 0;
         foreach (var inputVariable in problemData.AllowedInputVariables) {
-          x[i, j++] = problemData.Dataset.GetDoubleValue(inputVariable, r);
+          x[i, j] = (problemData.Dataset.GetDoubleValue(inputVariable, r) - means[j]) / stdDevs[j];
+          j++;
         }
         i++;
       }
+
+      this.random = random;
+      this.lambda = lambda;
+
       // initialize ERC values
       erc = Enumerable.Range(0, 10).Select(_ => Rand.RandNormal(random) * 10).ToArray();
@@ -132,5 +172,4 @@
         return OptimizeConstantsAndEvaluate(sentence);
       else {
-
         Debug.Assert(SimpleEvaluate(sentence) == SimpleEvaluate(extender.CanonicalRepresentation(sentence)));
         return SimpleEvaluate(sentence);
@@ -154,7 +193,7 @@
 
     public IEnumerable<Feature> GetFeatures(string phrase) {
-      // throw new NotImplementedException();
-      phrase = CanonicalRepresentation(phrase);
-      return phrase.Split('+').Distinct().Select(t => new Feature(t, 1.0));
+      throw new NotImplementedException();
+      //phrase = CanonicalRepresentation(phrase);
+      //return phrase.Split('+').Distinct().Select(t => new Feature(t, 1.0));
       // return new Feature[] { new Feature(phrase, 1.0) };
     }
@@ -176,39 +215,72 @@
       if (!constants.Any()) return SimpleEvaluate(sentence);
 
-      AutoDiff.IParametricCompiledTerm compiledFunc = func.Compile(constants, variables); // variate constants leave variables fixed to data
-
-      double[] c = constants.Select(_ => 1.0).ToArray(); // start with ones
-
+      // L2 regularization
+      // not possible with lsfit, would need to change to minlm below
+      // func = TermBuilder.Sum(func, lambda * TermBuilder.Sum(constants.Select(con => con * con)));
+
+      AutoDiff.IParametricCompiledTerm compiledFunc = func.Compile(constants, variables); // variate constants, leave variables fixed to data
+
+      // 10 restarts with random starting points
+      double[] bestStart = null;
+      double bestError = double.PositiveInfinity;
+      int info;
       alglib.lsfitstate state;
       alglib.lsfitreport rep;
-      int info;
-
-
-      int k = c.Length;
-
       alglib.ndimensional_pfunc function_cx_1_func = CreatePFunc(compiledFunc);
       alglib.ndimensional_pgrad function_cx_1_grad = CreatePGrad(compiledFunc);
-
-      const int maxIterations = 10;
-      try {
-        alglib.lsfitcreatefg(x, y, c, n, m, k, false, out state);
-        alglib.lsfitsetcond(state, 0.0, 0.0, maxIterations);
-        //alglib.lsfitsetgradientcheck(state, 0.001);
-        alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null);
-        alglib.lsfitresults(state, out info, out c, out rep);
-      } catch (ArithmeticException) {
-        return 0.0;
-      } catch (alglib.alglibexception) {
-        return 0.0;
-      }
-
-      //info == -7  => constant optimization failed due to wrong gradient
-      if (info == -7) throw new ArgumentException();
+      for (int t = 0; t < 10; t++) {
+        double[] cStart = constants.Select(_ => Rand.RandNormal(random) * 10).ToArray();
+        double[] cEnd;
+        // start with normally distributed (N(0, 10)) weights
+
+
+        int k = cStart.Length;
+
+
+        const int maxIterations = 10;
+        try {
+          alglib.lsfitcreatefg(x, y, cStart, n, m, k, false, out state);
+          alglib.lsfitsetcond(state, 0.0, 0.0, maxIterations);
+          //alglib.lsfitsetgradientcheck(state, 0.001);
+          alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null);
+          alglib.lsfitresults(state, out info, out cEnd, out rep);
+          if (info != -7 && rep.rmserror < bestError) {
+            bestStart = cStart;
+            bestError = rep.rmserror;
+          }
+        } catch (ArithmeticException) {
+          return 0.0;
+        } catch (alglib.alglibexception) {
+          return 0.0;
+        }
+      }
+
+      // 100 iteration steps from the best starting point
       {
-        var rowData = new double[d];
-        return HeuristicLab.Common.Extensions.RSq(y, Enumerable.Range(0, N).Select(i => {
-          for (int j = 0; j < d; j++) rowData[j] = x[i, j];
-          return compiledFunc.Evaluate(c, rowData);
-        }));
+        double[] c = bestStart;
+
+        int k = c.Length;
+
+        const int maxIterations = 100;
+        try {
+          alglib.lsfitcreatefg(x, y, c, n, m, k, false, out state);
+          alglib.lsfitsetcond(state, 0.0, 0.0, maxIterations);
+          //alglib.lsfitsetgradientcheck(state, 0.001);
+          alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null);
+          alglib.lsfitresults(state, out info, out c, out rep);
+        } catch (ArithmeticException) {
+          return 0.0;
+        } catch (alglib.alglibexception) {
+          return 0.0;
+        }
+        //info == -7  => constant optimization failed due to wrong gradient
+        if (info == -7) throw new ArgumentException();
+        {
+          var rowData = new double[d];
+          return HeuristicLab.Common.Extensions.RSq(y, Enumerable.Range(0, N).Select(i => {
+            for (int j = 0; j < d; j++) rowData[j] = x[i, j];
+            return compiledFunc.Evaluate(c, rowData);
+          }));
+        }
       }
     }
