Index: branches/symbreg-factors-2650/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/SymbolicDataAnalysisExpressionTreeSimplifier.cs
===================================================================
--- branches/symbreg-factors-2650/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/SymbolicDataAnalysisExpressionTreeSimplifier.cs	(revision 14249)
+++ branches/symbreg-factors-2650/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/SymbolicDataAnalysisExpressionTreeSimplifier.cs	(revision 14251)
@@ -1,3 +1,4 @@
 ﻿#region License Information
+
 /* HeuristicLab
  * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
@@ -18,4 +19,5 @@
  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
  */
+
 #endregion
 
@@ -24,4 +26,5 @@
 using System.Linq;
 using HeuristicLab.Common;
+using HeuristicLab.Core;
 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
 
@@ -66,5 +69,6 @@
 
     // the argumentTrees list contains already expanded trees used as arguments for invocations
-    private ISymbolicExpressionTreeNode MacroExpand(ISymbolicExpressionTreeNode root, ISymbolicExpressionTreeNode node, IList<ISymbolicExpressionTreeNode> argumentTrees) {
+    private ISymbolicExpressionTreeNode MacroExpand(ISymbolicExpressionTreeNode root, ISymbolicExpressionTreeNode node,
+      IList<ISymbolicExpressionTreeNode> argumentTrees) {
       List<ISymbolicExpressionTreeNode> subtrees = new List<ISymbolicExpressionTreeNode>(node.Subtrees);
       while (node.SubtreeCount > 0) node.RemoveSubtree(0);
@@ -98,6 +102,6 @@
     }
 
-
     #region symbol predicates
+
     // arithmetic
     private bool IsDivision(ISymbolicExpressionTreeNode node) {
@@ -120,50 +124,65 @@
       return node.Symbol is Average;
     }
+
     // exponential
     private bool IsLog(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Logarithm;
     }
+
     private bool IsExp(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Exponential;
     }
+
     private bool IsRoot(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Root;
     }
+
     private bool IsSquare(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Square;
     }
+
     private bool IsSquareRoot(ISymbolicExpressionTreeNode node) {
       return node.Symbol is SquareRoot;
     }
+
     private bool IsPower(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Power;
     }
+
     // trigonometric
     private bool IsSine(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Sine;
     }
+
     private bool IsCosine(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Cosine;
     }
+
     private bool IsTangent(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Tangent;
     }
+
     // boolean
     private bool IsIfThenElse(ISymbolicExpressionTreeNode node) {
       return node.Symbol is IfThenElse;
     }
+
     private bool IsAnd(ISymbolicExpressionTreeNode node) {
       return node.Symbol is And;
     }
+
     private bool IsOr(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Or;
     }
+
     private bool IsNot(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Not;
     }
+
     // comparison
     private bool IsGreaterThan(ISymbolicExpressionTreeNode node) {
       return node.Symbol is GreaterThan;
     }
+
     private bool IsLessThan(ISymbolicExpressionTreeNode node) {
       return node.Symbol is LessThan;
@@ -182,7 +201,13 @@
       return node.Symbol is Variable;
     }
+
     private bool IsVariableBase(ISymbolicExpressionTreeNode node) {
       return node is VariableTreeNodeBase;
     }
+
+    private bool IsFactor(ISymbolicExpressionTreeNode node) {
+      return node is FactorVariableTreeNode;
+    }
+
     private bool IsConstant(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Constant;
@@ -193,4 +218,5 @@
       return node.Symbol is TimeLag;
     }
+
     private bool IsIntegral(ISymbolicExpressionTreeNode node) {
       return node.Symbol is Integral;
@@ -256,6 +282,6 @@
     }
 
-
     #region specific simplification routines
+
     private ISymbolicExpressionTreeNode SimplifyAny(ISymbolicExpressionTreeNode original) {
       // can't simplify this function but simplify all subtrees 
@@ -305,5 +331,6 @@
         var remaining = original.Subtrees.Skip(2);
         return
-          MakeProduct(GetSimplifiedTree(first), Invert(remaining.Aggregate(GetSimplifiedTree(second), (a, b) => MakeProduct(a, GetSimplifiedTree(b)))));
+          MakeProduct(GetSimplifiedTree(first),
+            Invert(remaining.Aggregate(GetSimplifiedTree(second), (a, b) => MakeProduct(a, GetSimplifiedTree(b)))));
       }
     }
@@ -346,4 +373,5 @@
       return MakeNot(GetSimplifiedTree(original.GetSubtree(0)));
     }
+
     private ISymbolicExpressionTreeNode SimplifyOr(ISymbolicExpressionTreeNode original) {
       return original.Subtrees
@@ -351,4 +379,5 @@
         .Aggregate(MakeOr);
     }
+
     private ISymbolicExpressionTreeNode SimplifyAnd(ISymbolicExpressionTreeNode original) {
       return original.Subtrees
@@ -356,28 +385,38 @@
         .Aggregate(MakeAnd);
     }
+
     private ISymbolicExpressionTreeNode SimplifyLessThan(ISymbolicExpressionTreeNode original) {
       return MakeLessThan(GetSimplifiedTree(original.GetSubtree(0)), GetSimplifiedTree(original.GetSubtree(1)));
     }
+
     private ISymbolicExpressionTreeNode SimplifyGreaterThan(ISymbolicExpressionTreeNode original) {
       return MakeGreaterThan(GetSimplifiedTree(original.GetSubtree(0)), GetSimplifiedTree(original.GetSubtree(1)));
     }
+
     private ISymbolicExpressionTreeNode SimplifyIfThenElse(ISymbolicExpressionTreeNode original) {
-      return MakeIfThenElse(GetSimplifiedTree(original.GetSubtree(0)), GetSimplifiedTree(original.GetSubtree(1)), GetSimplifiedTree(original.GetSubtree(2)));
-    }
+      return MakeIfThenElse(GetSimplifiedTree(original.GetSubtree(0)), GetSimplifiedTree(original.GetSubtree(1)),
+        GetSimplifiedTree(original.GetSubtree(2)));
+    }
+
     private ISymbolicExpressionTreeNode SimplifyTangent(ISymbolicExpressionTreeNode original) {
       return MakeTangent(GetSimplifiedTree(original.GetSubtree(0)));
     }
+
     private ISymbolicExpressionTreeNode SimplifyCosine(ISymbolicExpressionTreeNode original) {
       return MakeCosine(GetSimplifiedTree(original.GetSubtree(0)));
     }
+
     private ISymbolicExpressionTreeNode SimplifySine(ISymbolicExpressionTreeNode original) {
       return MakeSine(GetSimplifiedTree(original.GetSubtree(0)));
     }
+
     private ISymbolicExpressionTreeNode SimplifyExp(ISymbolicExpressionTreeNode original) {
       return MakeExp(GetSimplifiedTree(original.GetSubtree(0)));
     }
+
     private ISymbolicExpressionTreeNode SimplifySquare(ISymbolicExpressionTreeNode original) {
       return MakeSquare(GetSimplifiedTree(original.GetSubtree(0)));
     }
+
     private ISymbolicExpressionTreeNode SimplifySquareRoot(ISymbolicExpressionTreeNode original) {
       return MakeSquareRoot(GetSimplifiedTree(original.GetSubtree(0)));
@@ -387,4 +426,5 @@
       return MakeLog(GetSimplifiedTree(original.GetSubtree(0)));
     }
+
     private ISymbolicExpressionTreeNode SimplifyRoot(ISymbolicExpressionTreeNode original) {
       return MakeRoot(GetSimplifiedTree(original.GetSubtree(0)), GetSimplifiedTree(original.GetSubtree(1)));
@@ -394,4 +434,5 @@
       return MakePower(GetSimplifiedTree(original.GetSubtree(0)), GetSimplifiedTree(original.GetSubtree(1)));
     }
+
     private ISymbolicExpressionTreeNode SimplifyTimeLag(ISymbolicExpressionTreeNode original) {
       var laggedTreeNode = original as ILaggedTreeNode;
@@ -403,4 +444,5 @@
       }
     }
+
     private ISymbolicExpressionTreeNode SimplifyIntegral(ISymbolicExpressionTreeNode original) {
       var laggedTreeNode = original as ILaggedTreeNode;
@@ -416,4 +458,5 @@
 
     #region low level tree restructuring
+
     private ISymbolicExpressionTreeNode MakeTimeLag(ISymbolicExpressionTreeNode subtree, int lag) {
       if (lag == 0) return subtree;
@@ -446,5 +489,6 @@
       } else if (!IsBoolean(t)) {
         var gtNode = gtSymbol.CreateTreeNode();
-        gtNode.AddSubtree(t); gtNode.AddSubtree(MakeConstant(0.0));
+        gtNode.AddSubtree(t);
+        gtNode.AddSubtree(MakeConstant(0.0));
         var notNode = notSymbol.CreateTreeNode();
         notNode.AddSubtree(gtNode);
@@ -486,4 +530,5 @@
       }
     }
+
     private ISymbolicExpressionTreeNode MakeAnd(ISymbolicExpressionTreeNode a, ISymbolicExpressionTreeNode b) {
       if (IsConstant(a) && IsConstant(b)) {
@@ -515,5 +560,7 @@
       }
     }
-    private ISymbolicExpressionTreeNode MakeLessThan(ISymbolicExpressionTreeNode leftSide, ISymbolicExpressionTreeNode rightSide) {
+
+    private ISymbolicExpressionTreeNode MakeLessThan(ISymbolicExpressionTreeNode leftSide,
+      ISymbolicExpressionTreeNode rightSide) {
       if (IsConstant(leftSide) && IsConstant(rightSide)) {
         var lsConst = leftSide as ConstantTreeNode;
@@ -528,5 +575,7 @@
       }
     }
-    private ISymbolicExpressionTreeNode MakeGreaterThan(ISymbolicExpressionTreeNode leftSide, ISymbolicExpressionTreeNode rightSide) {
+
+    private ISymbolicExpressionTreeNode MakeGreaterThan(ISymbolicExpressionTreeNode leftSide,
+      ISymbolicExpressionTreeNode rightSide) {
       if (IsConstant(leftSide) && IsConstant(rightSide)) {
         var lsConst = leftSide as ConstantTreeNode;
@@ -541,5 +590,7 @@
       }
     }
-    private ISymbolicExpressionTreeNode MakeIfThenElse(ISymbolicExpressionTreeNode condition, ISymbolicExpressionTreeNode trueBranch, ISymbolicExpressionTreeNode falseBranch) {
+
+    private ISymbolicExpressionTreeNode MakeIfThenElse(ISymbolicExpressionTreeNode condition,
+      ISymbolicExpressionTreeNode trueBranch, ISymbolicExpressionTreeNode falseBranch) {
       if (IsConstant(condition)) {
         var constT = condition as ConstantTreeNode;
@@ -552,5 +603,6 @@
         } else {
           var gtNode = gtSymbol.CreateTreeNode();
-          gtNode.AddSubtree(condition); gtNode.AddSubtree(MakeConstant(0.0));
+          gtNode.AddSubtree(condition);
+          gtNode.AddSubtree(MakeConstant(0.0));
           ifNode.AddSubtree(gtNode);
         }
@@ -571,4 +623,5 @@
       }
     }
+
     private ISymbolicExpressionTreeNode MakeTangent(ISymbolicExpressionTreeNode node) {
       if (IsConstant(node)) {
@@ -581,4 +634,5 @@
       }
     }
+
     private ISymbolicExpressionTreeNode MakeCosine(ISymbolicExpressionTreeNode node) {
       if (IsConstant(node)) {
@@ -591,8 +645,13 @@
       }
     }
+
     private ISymbolicExpressionTreeNode MakeExp(ISymbolicExpressionTreeNode node) {
       if (IsConstant(node)) {
         var constT = node as ConstantTreeNode;
         return MakeConstant(Math.Exp(constT.Value));
+      } else if (IsFactor(node)) {
+        var factNode = node as FactorVariableTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++) factNode.Weights[i] = Math.Exp(factNode.Weights[i]);
+        return factNode;
       } else if (IsLog(node)) {
         return node.GetSubtree(0);
@@ -612,4 +671,8 @@
         var constT = node as ConstantTreeNode;
         return MakeConstant(constT.Value * constT.Value);
+      } else if (IsFactor(node)) {
+        var factNode = node as FactorVariableTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++) factNode.Weights[i] = factNode.Weights[i] * factNode.Weights[i];
+        return factNode;
       } else if (IsSquareRoot(node)) {
         return node.GetSubtree(0);
@@ -620,8 +683,13 @@
       }
     }
+
     private ISymbolicExpressionTreeNode MakeSquareRoot(ISymbolicExpressionTreeNode node) {
       if (IsConstant(node)) {
         var constT = node as ConstantTreeNode;
         return MakeConstant(Math.Sqrt(constT.Value));
+      } else if (IsFactor(node)) {
+        var factNode = node as FactorVariableTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++) factNode.Weights[i] = Math.Sqrt(factNode.Weights[i]);
+        return factNode;
       } else if (IsSquare(node)) {
         return node.GetSubtree(0);
@@ -637,4 +705,8 @@
         var constT = node as ConstantTreeNode;
         return MakeConstant(Math.Log(constT.Value));
+      } else if (IsFactor(node)) {
+        var factNode = node as FactorVariableTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++) factNode.Weights[i] = Math.Log(factNode.Weights[i]);
+        return factNode;
       } else if (IsExp(node)) {
         return node.GetSubtree(0);
@@ -647,4 +719,5 @@
       }
     }
+
     private ISymbolicExpressionTreeNode MakeRoot(ISymbolicExpressionTreeNode a, ISymbolicExpressionTreeNode b) {
       if (IsConstant(a) && IsConstant(b)) {
@@ -652,4 +725,22 @@
         var constB = b as ConstantTreeNode;
         return MakeConstant(Math.Pow(constA.Value, 1.0 / Math.Round(constB.Value)));
+      } else if (IsFactor(a) && IsConstant(b)) {
+        var factNode = a as FactorVariableTreeNode;
+        var constNode = b as ConstantTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++)
+          factNode.Weights[i] = Math.Pow(factNode.Weights[i], 1.0 / Math.Round(constNode.Value));
+        return factNode;
+      } else if (IsConstant(a) && IsFactor(b)) {
+        var constNode = a as ConstantTreeNode;
+        var factNode = b as FactorVariableTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++)
+          factNode.Weights[i] = Math.Pow(constNode.Value, 1.0 / Math.Round(factNode.Weights[i]));
+        return factNode;
+      } else if (IsFactor(a) && IsFactor(b) && AreSameVariable(a, b)) {
+        var node0 = a as FactorVariableTreeNode;
+        var node1 = b as FactorVariableTreeNode;
+        for (int i = 0; i < node0.Weights.Length; i++)
+          node0.Weights[i] = Math.Pow(node0.Weights[i], 1.0 / Math.Round(node1.Weights[i]));
+        return node0;
       } else if (IsConstant(b)) {
         var constB = b as ConstantTreeNode;
@@ -679,4 +770,5 @@
       }
     }
+
     private ISymbolicExpressionTreeNode MakePower(ISymbolicExpressionTreeNode a, ISymbolicExpressionTreeNode b) {
       if (IsConstant(a) && IsConstant(b)) {
@@ -684,4 +776,22 @@
         var constB = b as ConstantTreeNode;
         return MakeConstant(Math.Pow(constA.Value, Math.Round(constB.Value)));
+      } else if (IsFactor(a) && IsConstant(b)) {
+        var factNode = a as FactorVariableTreeNode;
+        var constNode = b as ConstantTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++)
+          factNode.Weights[i] = Math.Pow(factNode.Weights[i], Math.Round(constNode.Value));
+        return factNode;
+      } else if (IsConstant(a) && IsFactor(b)) {
+        var constNode = a as ConstantTreeNode;
+        var factNode = b as FactorVariableTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++)
+          factNode.Weights[i] = Math.Pow(constNode.Value, Math.Round(factNode.Weights[i]));
+        return factNode;
+      } else if (IsFactor(a) && IsFactor(b) && AreSameVariable(a, b)) {
+        var node0 = a as FactorVariableTreeNode;
+        var node1 = b as FactorVariableTreeNode;
+        for (int i = 0; i < node0.Weights.Length; i++)
+          node0.Weights[i] = Math.Pow(node0.Weights[i], Math.Round(node1.Weights[i]));
+        return node0;
       } else if (IsConstant(b)) {
         var constB = b as ConstantTreeNode;
@@ -718,5 +828,5 @@
         // fold constants
         return MakeConstant(((ConstantTreeNode)a).Value / ((ConstantTreeNode)b).Value);
-      } if (IsConstant(a) && !((ConstantTreeNode)a).Value.IsAlmost(1.0)) {
+      } else if ((IsConstant(a) && !((ConstantTreeNode)a).Value.IsAlmost(1.0))) {
         return MakeFraction(MakeConstant(1.0), MakeProduct(b, Invert(a)));
       } else if (IsVariableBase(a) && IsConstant(b)) {
@@ -725,4 +835,16 @@
         ((VariableTreeNodeBase)a).Weight /= constB;
         return a;
+      } else if (IsFactor(a) && IsConstant(b)) {
+        var factNode = a as FactorVariableTreeNode;
+        var constNode = b as ConstantTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++) factNode.Weights[i] = factNode.Weights[i] / constNode.Value;
+        return factNode;
+      } else if (IsFactor(a) && IsFactor(b) && AreSameVariable(a, b)) {
+        var node0 = a as FactorVariableTreeNode;
+        var node1 = b as FactorVariableTreeNode;
+        for (int i = 0; i < node0.Weights.Length; i++) node0.Weights[i] = node0.Weights[i] / node1.Weights[i];
+        return node0;
+      } else if (IsFactor(a)) {
+        return MakeFraction(MakeConstant(1.0), MakeProduct(b, Invert(a)));
       } else if (IsVariableBase(a) && IsVariableBase(b) && AreSameVariable(a, b)) {
         // cancel variables
@@ -733,6 +855,6 @@
         return a.Subtrees
           .Select(x => GetSimplifiedTree(x))
-         .Select(x => MakeFraction(x, b))
-         .Aggregate((c, d) => MakeSum(c, d));
+          .Select(x => MakeFraction(x, b))
+          .Aggregate((c, d) => MakeSum(c, d));
       } else if (IsMultiplication(a) && IsConstant(b)) {
         return MakeProduct(a, Invert(b));
@@ -769,4 +891,14 @@
         // x + 0 => x
         return a;
+      } else if (IsFactor(a) && IsConstant(b)) {
+        var factNode = a as FactorVariableTreeNode;
+        var constNode = b as ConstantTreeNode;
+        for (int i = 0; i < factNode.Weights.Length; i++) factNode.Weights[i] = factNode.Weights[i] + constNode.Value;
+        return factNode;
+      } else if (IsFactor(a) && IsFactor(b) && AreSameVariable(a, b)) {
+        var node0 = a as FactorVariableTreeNode;
+        var node1 = b as FactorVariableTreeNode;
+        for (int i = 0; i < node0.Weights.Length; i++) node0.Weights[i] = node0.Weights[i] + node1.Weights[i];
+        return node0;
       } else if (IsAddition(a) && IsAddition(b)) {
         // merge additions
@@ -835,19 +967,42 @@
       var subtrees = new List<ISymbolicExpressionTreeNode>(sum.Subtrees);
       while (sum.Subtrees.Any()) sum.RemoveSubtree(0);
-      var groupedVarNodes = from node in subtrees.OfType<VariableTreeNodeBase>()
+      var groupedVarNodes = from node in subtrees.OfType<IVariableTreeNode>()
+                            where node.SubtreeCount == 0
+                            // only consider terminal nodes
                             let lag = (node is LaggedVariableTreeNode) ? ((LaggedVariableTreeNode)node).Lag : 0
-                            let cat = (node is BinaryFactorVariableTreeNode) ? ((BinaryFactorVariableTreeNode)node).VariableValue : string.Empty
-                            group node by node.VariableName + cat + lag into g
+                            let cat =
+                              (node is BinaryFactorVariableTreeNode) ? ((BinaryFactorVariableTreeNode)node).VariableValue : string.Empty
+                            group node by node.VariableName + cat + lag
+        into g
                             select g;
-      var unchangedSubtrees = subtrees.Where(t => !(t is VariableTreeNodeBase));
+      var constant = (from node in subtrees.OfType<ConstantTreeNode>()
+                      select node.Value).DefaultIfEmpty(0.0).Sum();
+      var unchangedSubtrees = subtrees.Where(t => t.SubtreeCount > 0 || !(t is IVariableTreeNode) && !(t is ConstantTreeNode));
 
       foreach (var variableNodeGroup in groupedVarNodes) {
-        var weightSum = variableNodeGroup.Select(t => t.Weight).Sum();
-        var representative = variableNodeGroup.First();
-        representative.Weight = weightSum;
-        sum.AddSubtree(representative);
+        var firstNode = variableNodeGroup.First();
+        if (firstNode is VariableTreeNodeBase) {
+          var representative = firstNode as VariableTreeNodeBase;
+          var weightSum = variableNodeGroup.Cast<VariableTreeNodeBase>().Select(t => t.Weight).Sum();
+          representative.Weight = weightSum;
+          sum.AddSubtree(representative);
+        } else if (firstNode is FactorVariableTreeNode) {
+          var representative = firstNode as FactorVariableTreeNode;
+          foreach (var node in variableNodeGroup.Skip(1).Cast<FactorVariableTreeNode>()) {
+            for (int j = 0; j < representative.Weights.Length; j++) {
+              representative.Weights[j] += node.Weights[j];
+            }
+          }
+          for (int j = 0; j < representative.Weights.Length; j++) {
+            representative.Weights[j] += constant;
+          }
+          sum.AddSubtree(representative);
+        }
       }
       foreach (var unchangedSubtree in unchangedSubtrees)
         sum.AddSubtree(unchangedSubtree);
+      if (!constant.IsAlmost(0.0)) {
+        sum.AddSubtree(MakeConstant(constant));
+      }
     }
 
@@ -861,4 +1016,14 @@
         // a * $ => $ * a
         return MakeProduct(b, a);
+      } else if (IsFactor(a) && IsFactor(b) && AreSameVariable(a, b)) {
+        var node0 = a as FactorVariableTreeNode;
+        var node1 = b as FactorVariableTreeNode;
+        for (int i = 0; i < node0.Weights.Length; i++) node0.Weights[i] = node0.Weights[i] * node1.Weights[i];
+        return node0;
+      } else if (IsFactor(a) && IsConstant(b)) {
+        var node0 = a as FactorVariableTreeNode;
+        var node1 = b as ConstantTreeNode;
+        for (int i = 0; i < node0.Weights.Length; i++) node0.Weights[i] = node0.Weights[i] * node1.Value;
+        return node0;
       } else if (IsConstant(b) && ((ConstantTreeNode)b).Value.IsAlmost(1.0)) {
         // $ * 1.0 => $
@@ -868,5 +1033,7 @@
         ((VariableTreeNodeBase)a).Weight *= ((ConstantTreeNode)b).Value;
         return a;
-      } else if (IsConstant(b) && IsAddition(a)) {
+      } else if (IsConstant(b) && IsAddition(a) ||
+          IsFactor(a) && IsAddition(b) ||
+          IsAddition(b) && IsFactor(a)) {
         // multiply constants into additions
         return a.Subtrees.Select(x => MakeProduct(x, b)).Aggregate((c, d) => MakeSum(c, d));
@@ -902,8 +1069,9 @@
       }
     }
+
     #endregion
 
-
     #region helper functions
+
     private bool ContainsVariableCondition(ISymbolicExpressionTreeNode node) {
       if (node.Symbol is VariableCondition) return true;
@@ -940,5 +1108,5 @@
       if (aLaggedVar != null && bLaggedVar != null) {
         return aLaggedVar.VariableName == bLaggedVar.VariableName &&
-          aLaggedVar.Lag == bLaggedVar.Lag;
+               aLaggedVar.Lag == bLaggedVar.Lag;
       }
       var aVar = a as VariableTreeNode;
@@ -947,11 +1115,15 @@
         return aVar.VariableName == bVar.VariableName;
       }
-      var aFactor = a as BinaryFactorVariableTreeNode;
-      var bFactor = b as BinaryFactorVariableTreeNode;
+      var aFactor = a as FactorVariableTreeNode;
+      var bFactor = b as FactorVariableTreeNode;
       if (aFactor != null && bFactor != null) {
-        return aFactor.VariableName == bFactor.VariableName &&
-          aFactor.VariableValue == bFactor.VariableValue;
-      }
-
+        return aFactor.VariableName == bFactor.VariableName;
+      }
+      var aBinFactor = a as BinaryFactorVariableTreeNode;
+      var bBinFactor = b as BinaryFactorVariableTreeNode;
+      if (aBinFactor != null && bBinFactor != null) {
+        return aBinFactor.VariableName == bBinFactor.VariableName &&
+               aBinFactor.VariableValue == bBinFactor.VariableValue;
+      }
       return false;
     }
@@ -961,31 +1133,50 @@
       var subtrees = new List<ISymbolicExpressionTreeNode>(prod.Subtrees);
       while (prod.Subtrees.Any()) prod.RemoveSubtree(0);
-      var groupedVarNodes = from node in subtrees.OfType<VariableTreeNodeBase>()
+      var groupedVarNodes = from node in subtrees.OfType<IVariableTreeNode>()
+                            where node.SubtreeCount == 0
                             let lag = (node is LaggedVariableTreeNode) ? ((LaggedVariableTreeNode)node).Lag : 0
-                            group node by node.VariableName + lag into g
+                            let cat =
+                              (node is BinaryFactorVariableTreeNode) ? ((BinaryFactorVariableTreeNode)node).VariableValue : string.Empty
+                            group node by node.VariableName + cat + lag
+        into g
                             orderby g.Count()
                             select g;
       var constantProduct = (from node in subtrees.OfType<VariableTreeNodeBase>()
                              select node.Weight)
-                            .Concat(from node in subtrees.OfType<ConstantTreeNode>()
-                                    select node.Value)
-                            .DefaultIfEmpty(1.0)
-                            .Aggregate((c1, c2) => c1 * c2);
+        .Concat(from node in subtrees.OfType<ConstantTreeNode>()
+                select node.Value)
+        .DefaultIfEmpty(1.0)
+        .Aggregate((c1, c2) => c1 * c2);
 
       var unchangedSubtrees = from tree in subtrees
-                              where !(tree is VariableTreeNodeBase)
-                              where !(tree is ConstantTreeNode)
+                              where tree.SubtreeCount > 0 || !(tree is IVariableTreeNode) && !(tree is ConstantTreeNode)
                               select tree;
 
       foreach (var variableNodeGroup in groupedVarNodes) {
-        var representative = variableNodeGroup.First();
-        representative.Weight = 1.0;
-        if (variableNodeGroup.Count() > 1) {
-          var poly = mulSymbol.CreateTreeNode();
-          for (int p = 0; p < variableNodeGroup.Count(); p++) {
-            poly.AddSubtree((ISymbolicExpressionTreeNode)representative.Clone());
+        var firstNode = variableNodeGroup.First();
+        if (firstNode is VariableTreeNodeBase) {
+          var representative = (VariableTreeNodeBase)firstNode;
+          representative.Weight = 1.0;
+          if (variableNodeGroup.Count() > 1) {
+            var poly = mulSymbol.CreateTreeNode();
+            for (int p = 0; p < variableNodeGroup.Count(); p++) {
+              poly.AddSubtree((ISymbolicExpressionTreeNode)representative.Clone());
+            }
+            prod.AddSubtree(poly);
+          } else {
+            prod.AddSubtree(representative);
           }
-          prod.AddSubtree(poly);
-        } else {
+        } else if (firstNode is FactorVariableTreeNode) {
+          var representative = (FactorVariableTreeNode)firstNode;
+          foreach (var node in variableNodeGroup.Skip(1).Cast<FactorVariableTreeNode>()) {
+            for (int j = 0; j < representative.Weights.Length; j++) {
+              representative.Weights[j] *= node.Weights[j];
+            }
+          }
+          for (int j = 0; j < representative.Weights.Length; j++) {
+            representative.Weights[j] *= constantProduct;
+          }
+          constantProduct = 1.0;
+          // if the product already contains a factor it is not necessary to multiply a constant below
           prod.AddSubtree(representative);
         }
@@ -1013,4 +1204,7 @@
         var variableTree = (VariableTreeNodeBase)x;
         variableTree.Weight *= -1.0;
+      } else if (IsFactor(x)) {
+        var factorNode = (FactorVariableTreeNode)x;
+        for (int i = 0; i < factorNode.Weights.Length; i++) factorNode.Weights[i] *= -1;
       } else if (IsAddition(x)) {
         // (x0 + x1 + .. + xn) * -1 => (-x0 + -x1 + .. + -xn)        
@@ -1041,4 +1235,8 @@
       if (IsConstant(x)) {
         return MakeConstant(1.0 / ((ConstantTreeNode)x).Value);
+      } else if (IsFactor(x)) {
+        var factorNode = (FactorVariableTreeNode)x;
+        for (int i = 0; i < factorNode.Weights.Length; i++) factorNode.Weights[i] = 1.0 / factorNode.Weights[i];
+        return factorNode;
       } else if (IsDivision(x)) {
         return MakeFraction(x.GetSubtree(1), x.GetSubtree(0));
@@ -1061,4 +1259,5 @@
       return tree;
     }
+
     #endregion
   }
