Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Encodings.SymbolicExpressionTreeEncoding/3.4/Crossovers/SubtreeCrossover.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Encodings.SymbolicExpressionTreeEncoding/3.4/Crossovers/SubtreeCrossover.cs	(revision 16615)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Encodings.SymbolicExpressionTreeEncoding/3.4/Crossovers/SubtreeCrossover.cs	(revision 16616)
@@ -43,4 +43,5 @@
     private const string MaximumSymbolicExpressionTreeLengthParameterName = "MaximumSymbolicExpressionTreeLength";
     private const string MaximumSymbolicExpressionTreeDepthParameterName = "MaximumSymbolicExpressionTreeDepth";
+    private const string CrossoverProbabilityParameterName = "CrossoverProbability";
 
     #region Parameter Properties
@@ -53,4 +54,7 @@
     public IValueLookupParameter<IntValue> MaximumSymbolicExpressionTreeDepthParameter {
       get { return (IValueLookupParameter<IntValue>)Parameters[MaximumSymbolicExpressionTreeDepthParameterName]; }
+    }
+    public IValueParameter<PercentValue> CrossoverProbabilityParameter {
+      get { return (IValueParameter<PercentValue>)Parameters[CrossoverProbabilityParameterName]; }
     }
     #endregion
@@ -64,4 +68,7 @@
     public IntValue MaximumSymbolicExpressionTreeDepth {
       get { return MaximumSymbolicExpressionTreeDepthParameter.ActualValue; }
+    }
+    public PercentValue CrossoverProbability {
+      get { return CrossoverProbabilityParameter.Value; }
     }
     #endregion
@@ -74,8 +81,16 @@
       Parameters.Add(new ValueLookupParameter<IntValue>(MaximumSymbolicExpressionTreeDepthParameterName, "The maximal depth of the symbolic expression tree (a tree with one node has depth = 0)."));
       Parameters.Add(new ValueLookupParameter<PercentValue>(InternalCrossoverPointProbabilityParameterName, "The probability to select an internal crossover point (instead of a leaf node).", new PercentValue(0.9)));
+      Parameters.Add(new ValueParameter<PercentValue>(CrossoverProbabilityParameterName, "The probability to apply crossover (default 100%).", new PercentValue(1)));
     }
 
     public override IDeepCloneable Clone(Cloner cloner) {
       return new SubtreeCrossover(this, cloner);
+    }
+
+    [StorableHook(HookType.AfterDeserialization)]
+    private void AfterDeserialization() {
+      if (!Parameters.ContainsKey(CrossoverProbabilityParameterName)) {
+        Parameters.Add(new ValueParameter<PercentValue>(CrossoverProbabilityParameterName, "The probability to apply crossover (default 100%).", new PercentValue(1)));
+      }
     }
 
@@ -83,10 +98,14 @@
       ISymbolicExpressionTree parent0, ISymbolicExpressionTree parent1) {
       return Cross(random, parent0, parent1, InternalCrossoverPointProbability.Value,
-        MaximumSymbolicExpressionTreeLength.Value, MaximumSymbolicExpressionTreeDepth.Value);
+        MaximumSymbolicExpressionTreeLength.Value, MaximumSymbolicExpressionTreeDepth.Value, CrossoverProbability.Value);
     }
 
     public static ISymbolicExpressionTree Cross(IRandom random,
       ISymbolicExpressionTree parent0, ISymbolicExpressionTree parent1,
-      double internalCrossoverPointProbability, int maxTreeLength, int maxTreeDepth) {
+      double internalCrossoverPointProbability, int maxTreeLength, int maxTreeDepth, double crossoverProbability = 1.0) {
+      if (random.NextDouble() >= crossoverProbability) {
+        return random.NextDouble() < 0.5 ? parent0 : parent1;
+      }
+
       // select a random crossover point in the first parent 
       CutPoint crossoverPoint0;
Index: /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs
===================================================================
--- /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 16615)
+++ /branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs	(revision 16616)
@@ -245,4 +245,5 @@
       nmse = OptimizeParameters(trees, problemData, targetVars, latentVariables, episodes, maxParameterOptIterations, pretunedParameters, numericIntegrationSteps, odeSolver,
         out double[] optTheta);
+      // var optTheta = pretunedParameters;
 
       if (double.IsNaN(nmse) ||
@@ -297,15 +298,22 @@
             alglib.minlmreport report;
             var p = new double[initialTheta[treeIdx].Length];
-            var lowerBounds = Enumerable.Repeat(-100.0, p.Length).ToArray();
-            var upperBounds = Enumerable.Repeat(100.0, p.Length).ToArray();
+            var lowerBounds = Enumerable.Repeat(-1000.0, p.Length).ToArray();
+            var upperBounds = Enumerable.Repeat(1000.0, p.Length).ToArray();
             Array.Copy(initialTheta[treeIdx], p, p.Length);
             alglib.minlmcreatevj(targetValuesDiff.Count, p, out state);
             alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
             alglib.minlmsetbc(state, lowerBounds, upperBounds);
-            // alglib.minlmsetgradientcheck(state, 1.0e-3);
+#if DEBUG
+            //alglib.minlmsetgradientcheck(state, 1.0e-7);
+#endif
             alglib.minlmoptimize(state, EvaluateObjectiveVector, EvaluateObjectiveVectorAndJacobian, null, myState);
 
             alglib.minlmresults(state, out optTheta, out report);
-            if (report.terminationtype < 0) { optTheta = initialTheta[treeIdx]; }
+            if (report.terminationtype < 0) {
+#if DEBUG
+              if (report.terminationtype == -7) throw new InvalidProgramException("gradient calculation fail!");
+#endif
+              optTheta = initialTheta[treeIdx];
+            }
           } catch (alglib.alglibexception) {
             optTheta = initialTheta[treeIdx];
@@ -342,6 +350,6 @@
 
       if (initialTheta.Length > 0) {
-        var lowerBounds = Enumerable.Repeat(-100.0, initialTheta.Length).ToArray();
-        var upperBounds = Enumerable.Repeat(100.0, initialTheta.Length).ToArray();
+        var lowerBounds = Enumerable.Repeat(-1000.0, initialTheta.Length).ToArray();
+        var upperBounds = Enumerable.Repeat(1000.0, initialTheta.Length).ToArray();
         try {
           alglib.minlmstate state;
@@ -350,5 +358,7 @@
           alglib.minlmsetbc(state, lowerBounds, upperBounds);
           alglib.minlmsetcond(state, 0.0, 0.0, 0.0, maxParameterOptIterations);
-          // alglib.minlmsetgradientcheck(state, 1.0e-3);
+#if DEBUG          
+          //alglib.minlmsetgradientcheck(state, 1.0e-7);
+#endif
           alglib.minlmoptimize(state, IntegrateAndEvaluateObjectiveVector, IntegrateAndEvaluateObjectiveVectorAndJacobian, null, myState);
 
@@ -356,5 +366,7 @@
 
           if (report.terminationtype < 0) {
-            // there was a problem: reset theta and evaluate for inital values
+#if DEBUG
+            if (report.terminationtype == -7) throw new InvalidProgramException("gradient calculation fail!");
+#endif            // there was a problem: reset theta and evaluate for inital values
             optTheta = initialTheta;
           }
@@ -397,5 +409,5 @@
           var pred = InterpretRec(tree.Root.GetSubtree(0).GetSubtree(0), nodeValueLookup);
           var y = optimizationData.targetValues[treeIdx][trainIdx];
-          fi[outputIdx++] = y - pred;
+          fi[outputIdx++] = (y - pred) * optimizationData.inverseStandardDeviation[treeIdx];
         }
       }
@@ -602,8 +614,8 @@
             var y = targetValues[colIdx][r];
 
-            var res = (y - y_pred_f);
+            var res = (y - y_pred_f) * optimizationData.inverseStandardDeviation[colIdx];
             var ressq = res * res;
-            f_i += ressq * optimizationData.inverseStandardDeviation[colIdx];
-            g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res * optimizationData.inverseStandardDeviation[colIdx]));
+            f_i += ressq;
+            g_i.Add(y_pred[colIdx].Item2.Scale(-2.0 * res));
           }
           seRow.Values.Add(f_i);
@@ -858,22 +870,22 @@
 
         var flag = CVODES.CVodeInit(cvode_mem, f, 0.0, y);
-        Debug.Assert(CVODES.CV_SUCCESS == flag);
+        Assert(CVODES.CV_SUCCESS == flag);
 
         double relTol = 1.0e-2;
         double absTol = 1.0;
         flag = CVODES.CVodeSStolerances(cvode_mem, relTol, absTol);  // TODO: probably need to adjust absTol per variable
-        Debug.Assert(CVODES.CV_SUCCESS == flag);
+        Assert(CVODES.CV_SUCCESS == flag);
 
         A = CVODES.SUNDenseMatrix(numberOfEquations, numberOfEquations);
-        Debug.Assert(A != IntPtr.Zero);
+        Assert(A != IntPtr.Zero);
 
         linearSolver = CVODES.SUNDenseLinearSolver(y, A);
-        Debug.Assert(linearSolver != IntPtr.Zero);
+        Assert(linearSolver != IntPtr.Zero);
 
         flag = CVODES.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
-        Debug.Assert(CVODES.CV_SUCCESS == flag);
+        Assert(CVODES.CV_SUCCESS == flag);
 
         flag = CVODES.CVDlsSetJacFn(cvode_mem, jac);
-        Debug.Assert(CVODES.CV_SUCCESS == flag);
+        Assert(CVODES.CV_SUCCESS == flag);
 
         yS0 = CVODES.N_VCloneVectorArray_Serial(ns, y); // clone the output vector for each parameter
@@ -889,8 +901,8 @@
 
         flag = CVODES.CVodeSensInit(cvode_mem, ns, CVODES.CV_SIMULTANEOUS, sensF, yS0);
-        Debug.Assert(CVODES.CV_SUCCESS == flag);
+        Assert(CVODES.CV_SUCCESS == flag);
 
         flag = CVODES.CVodeSensEEtolerances(cvode_mem);
-        Debug.Assert(CVODES.CV_SUCCESS == flag);
+        Assert(CVODES.CV_SUCCESS == flag);
 
         // make one forward integration step
@@ -898,9 +910,9 @@
         flag = CVODES.CVode(cvode_mem, t, y, ref tout, CVODES.CV_NORMAL);
         if (flag == CVODES.CV_SUCCESS) {
-          Debug.Assert(t == tout);
+          Assert(t == tout);
 
           // get sensitivities
           flag = CVODES.CVodeGetSens(cvode_mem, ref tout, yS0);
-          Debug.Assert(CVODES.CV_SUCCESS == flag);
+          Assert(CVODES.CV_SUCCESS == flag);
 
           // update variableValues based on integration results
@@ -1149,15 +1161,15 @@
         return nodeValues.NodeValue(node);
       } else if (node.Symbol is Addition) {
-        Debug.Assert(node.SubtreeCount == 2);
+        Assert(node.SubtreeCount == 2);
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         var g = InterpretRec(node.GetSubtree(1), nodeValues);
         return f + g;
       } else if (node.Symbol is Multiplication) {
-        Debug.Assert(node.SubtreeCount == 2);
+        Assert(node.SubtreeCount == 2);
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         var g = InterpretRec(node.GetSubtree(1), nodeValues);
         return f * g;
       } else if (node.Symbol is Subtraction) {
-        Debug.Assert(node.SubtreeCount <= 2);
+        Assert(node.SubtreeCount <= 2);
         if (node.SubtreeCount == 1) {
           var f = InterpretRec(node.GetSubtree(0), nodeValues);
@@ -1170,5 +1182,5 @@
         }
       } else if (node.Symbol is Division) {
-        Debug.Assert(node.SubtreeCount <= 2);
+        Assert(node.SubtreeCount <= 2);
 
         if (node.SubtreeCount == 1) {
@@ -1192,29 +1204,35 @@
         }
       } else if (node.Symbol is Sine) {
-        Debug.Assert(node.SubtreeCount == 1);
+        Assert(node.SubtreeCount == 1);
 
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         return Math.Sin(f);
       } else if (node.Symbol is Cosine) {
-        Debug.Assert(node.SubtreeCount == 1);
+        Assert(node.SubtreeCount == 1);
 
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         return Math.Cos(f);
       } else if (node.Symbol is Square) {
-        Debug.Assert(node.SubtreeCount == 1);
+        Assert(node.SubtreeCount == 1);
 
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         return f * f;
       } else if (node.Symbol is Exponential) {
-        Debug.Assert(node.SubtreeCount == 1);
+        Assert(node.SubtreeCount == 1);
 
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         return Math.Exp(f);
       } else if (node.Symbol is Logarithm) {
-        Debug.Assert(node.SubtreeCount == 1);
+        Assert(node.SubtreeCount == 1);
 
         var f = InterpretRec(node.GetSubtree(0), nodeValues);
         return Math.Log(f);
       } else throw new NotSupportedException("unsupported symbol");
+    }
+
+    private static void Assert(bool cond) {
+#if DEBUG
+      if (!cond) throw new InvalidOperationException("Assertion failed");
+#endif
     }
 
@@ -1231,4 +1249,6 @@
         dz = Vector.CreateNew(nodeValues.NodeGradient(node)); // original gradient vectors are never changed by evaluation
       } else if (node.Symbol is Addition) {
+
+        Assert(node.SubtreeCount == 2);
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
@@ -1236,4 +1256,6 @@
         dz = df.Add(dg);
       } else if (node.Symbol is Multiplication) {
+
+        Assert(node.SubtreeCount == 2);
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         InterpretRec(node.GetSubtree(1), nodeValues, out g, out dg);
@@ -1242,4 +1264,6 @@
 
       } else if (node.Symbol is Subtraction) {
+
+        Assert(node.SubtreeCount <= 2);
         if (node.SubtreeCount == 1) {
           InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
@@ -1254,4 +1278,6 @@
 
       } else if (node.Symbol is Division) {
+
+        Assert(node.SubtreeCount <= 2);
         if (node.SubtreeCount == 1) {
           InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
@@ -1281,4 +1307,6 @@
 
       } else if (node.Symbol is Sine) {
+
+        Assert(node.SubtreeCount == 1);
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         z = Math.Sin(f);
@@ -1286,16 +1314,24 @@
 
       } else if (node.Symbol is Cosine) {
+
+        Assert(node.SubtreeCount == 1);
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         z = Math.Cos(f);
         dz = df.Scale(-Math.Sin(f));
       } else if (node.Symbol is Square) {
+
+        Assert(node.SubtreeCount == 1);
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         z = f * f;
         dz = df.Scale(2.0 * f);
       } else if (node.Symbol is Exponential) {
+
+        Assert(node.SubtreeCount == 1);
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         z = Math.Exp(f);
         dz = df.Scale(Math.Exp(f));
       } else if (node.Symbol is Logarithm) {
+
+        Assert(node.SubtreeCount == 1);
         InterpretRec(node.GetSubtree(0), nodeValues, out f, out df);
         z = Math.Log(f);
@@ -1443,4 +1479,11 @@
         // make sure our multi-manipulator is the only manipulator
         e.Operators = new IOperator[] { multiManipulator }.Concat(filteredOperators);
+
+        // set the crossover probability to reduce likelihood that multiple trees are crossed at the same time
+        var subtreeCrossovers = e.Operators.OfType<SubtreeCrossover>();
+        foreach (var xover in subtreeCrossovers) {
+          xover.CrossoverProbability.Value = 0.3;
+        }
+
         encoding = encoding.Add(e); // only limit by length
       }
@@ -1451,4 +1494,11 @@
         // make sure our multi-manipulator is the only manipulator
         e.Operators = new IOperator[] { multiManipulator }.Concat(filteredOperators);
+
+        // set the crossover probability to reduce likelihood that multiple trees are crossed at the same time
+        var subtreeCrossovers = e.Operators.OfType<SubtreeCrossover>();
+        foreach (var xover in subtreeCrossovers) {
+          xover.CrossoverProbability.Value = 0.3;
+        }
+
         encoding = encoding.Add(e);
       }
@@ -1589,19 +1639,19 @@
         this.targetValues = targetValues;
         this.variables = inputVariables;
-        // if (targetValues != null) {
-        //   this.inverseStandardDeviation = new double[targetValues.Length];
-        //   for (int i = 0; i < targetValues.Length; i++) {
-        //     // calculate variance for each episode separately and calc the average 
-        //     var epStartIdx = 0;
-        //     var stdevs = new List<double>();
-        //     foreach (var ep in episodes) {
-        //       var epValues = targetValues[i].Skip(epStartIdx).Take(ep.Size);
-        //       stdevs.Add(epValues.StandardDeviation());
-        //       epStartIdx += ep.Size;
-        //     }
-        //     inverseStandardDeviation[i] = 1.0 / stdevs.Average();
-        //   }
-        // } else
-        this.inverseStandardDeviation = Enumerable.Repeat(1.0, trees.Length).ToArray();
+        if (targetValues != null) {
+          this.inverseStandardDeviation = new double[targetValues.Length];
+          for (int i = 0; i < targetValues.Length; i++) {
+            // calculate variance for each episode separately and calc the average 
+            var epStartIdx = 0;
+            var stdevs = new List<double>();
+            foreach (var ep in episodes) {
+              var epValues = targetValues[i].Skip(epStartIdx).Take(ep.Size);
+              stdevs.Add(epValues.StandardDeviation());
+              epStartIdx += ep.Size;
+            }
+            inverseStandardDeviation[i] = 1.0 / stdevs.Average();
+          }
+        } else
+          this.inverseStandardDeviation = Enumerable.Repeat(1.0, trees.Length).ToArray();
         this.episodes = episodes;
         this.numericIntegrationSteps = numericIntegrationSteps;
