diff --git a/eo/contrib/mathsym/GNUmakefile b/eo/contrib/mathsym/GNUmakefile
index c3f635c4..880c72c4 100644
--- a/eo/contrib/mathsym/GNUmakefile
+++ b/eo/contrib/mathsym/GNUmakefile
@@ -1,7 +1,8 @@
 COMPILEFLAGS=-Wno-deprecated -g -Wall #-DINTERVAL_DEBUG
-OPTFLAGS= #-O3 -DNDEBUG 
+OPTFLAGS= -O3 DNDEBUG 
 
 PROFILE_FLAGS=#-pg 
+LDFLAGS=#-a
 
 INCLUDES=-I. -Isym -Ifun -Igen -Ieval -Iregression -I../../src -Ieo_interface  -I..
 
@@ -18,7 +19,7 @@ CXXSOURCES=FunDef.cpp Sym.cpp SymImpl.cpp SymOps.cpp sym_compile.cpp TreeBuilder
 	Dataset.cpp ErrorMeasure.cpp Scaling.cpp TargetInfo.cpp BoundsCheck.cpp util.cpp NodeSelector.cpp\
 	eoSymCrossover.cpp sym_operations.cpp eoSymMutate.cpp
 
-TESTPROGRAMS=test/test_compile test/testeo test/test_simplify test/test_diff
+TESTPROGRAMS=test/test_compile test/testeo test/test_simplify test/test_diff test/test_lambda
 
 OBJS= $(CXXSOURCES:.cpp=.o) c_compile.o
 
@@ -27,19 +28,19 @@ all: tcc/ symreg
 include $(CXXSOURCES:.cpp=.d) symreg.d 
 
 clean:
-	rm *.o *.d $(TESTPROGRAMS) $(SYMLIB) symreg || true
+	rm *.o *.d $(TESTPROGRAMS) $(SYMLIB) symreg test/*.o || true
 
 distclean: clean
 	rm -rf tcc
 
 symreg: libsym.a symreg.o $(EXTLIBS)
-	$(CXX) -o symreg symreg.o libsym.a $(LIBS) $(PROFILE_FLAGS) 
+	$(CXX) -o symreg symreg.o libsym.a $(LIBS) $(PROFILE_FLAGS) ${LDFLAGS}
 
 libsym.a: $(OBJS) 
 	rm libsym.a; ar cq $(SYMLIB) $(OBJS) 
 
 check: $(TESTPROGRAMS)
-	test/test_compile && test/testeo && test/test_simplify && test/test_diff && echo "all tests succeeded"
+	test/test_compile && test/testeo && test/test_simplify && test/test_diff && test/test_lambda && echo "all tests succeeded"
 
 test/test_compile: test/test_compile.o ${SYMLIB}
 	$(CXX) -o test/test_compile test/test_compile.o $(SYMLIB) ${LIBS}
@@ -53,6 +54,9 @@ test/test_simplify: test/test_simplify.o $(SYMLIB)
 test/test_diff: test/test_diff.o $(SYMLIB)
 	$(CXX) -o test/test_diff test/test_diff.o  $(SYMLIB) ${LIBS}
 
+test/test_lambda: test/test_lambda.o $(SYMLIB)
+	$(CXX) -o test/test_lambda test/test_lambda.o  $(SYMLIB) ${LIBS}
+
 # eo
 ../../src/libeo.a:
 	make -C ../../src
diff --git a/eo/contrib/mathsym/eo_interface/eoSymCrossover.cpp b/eo/contrib/mathsym/eo_interface/eoSymCrossover.cpp
index 63d46c20..412ac78e 100644
--- a/eo/contrib/mathsym/eo_interface/eoSymCrossover.cpp
+++ b/eo/contrib/mathsym/eo_interface/eoSymCrossover.cpp
@@ -23,11 +23,11 @@
 #include <utils/eoRNG.h>
 
 bool subtree_quad(Sym& a, Sym& b, NodeSelector& select) {
-    unsigned i = select.select_node(a);
-    unsigned j = select.select_node(b);
+    NodeSelector::NodeSelection sel_a = select.select_node(a);
+    NodeSelector::NodeSelection sel_b = select.select_node(b);
 
-    Sym aprime = insert_subtree(a, i, get_subtree(b, j));
-    Sym bprime = insert_subtree(b, j, get_subtree(a, i));
+    Sym aprime = insert_subtree(a, sel_a.idx(), sel_b.subtree() );
+    Sym bprime = insert_subtree(b, sel_b.idx(), sel_a.subtree() );
 
     a = aprime;
     b = bprime;
@@ -35,16 +35,18 @@ bool subtree_quad(Sym& a, Sym& b, NodeSelector& select) {
 }
 
 bool subtree_bin(Sym& a, const Sym& b, NodeSelector& select) {
-    unsigned i = select.select_node(a);
-    unsigned j = select.select_node(b);
+    NodeSelector::NodeSelection sel_a = select.select_node(a);
+    NodeSelector::NodeSelection sel_b = select.select_node(b);
 
-    a = insert_subtree(a, i, get_subtree(b,j));
+    a = insert_subtree(a, sel_a.idx(), sel_b.subtree());
 
     return true;
 }
 
 Sym homologous_binimpl(Sym a, Sym b) {
-
+    
+    if(a == b) { return a; } // no point 
+    
     bool use_a = rng.random(2);
 
     token_t head = (use_a? a : b).token();
diff --git a/eo/contrib/mathsym/eo_interface/eoSymMutate.h b/eo/contrib/mathsym/eo_interface/eoSymMutate.h
index 33623a12..aabe9141 100644
--- a/eo/contrib/mathsym/eo_interface/eoSymMutate.h
+++ b/eo/contrib/mathsym/eo_interface/eoSymMutate.h
@@ -36,7 +36,7 @@ class eoSymSubtreeMutate : public eoMonOp<EoType> {
 	
 	
 	bool operator()(EoType& tomutate) {
-	    unsigned xover_point = node_selector.select_node(tomutate);
+	    unsigned xover_point = node_selector.select_node(tomutate).idx();
 	    // create subtree
 	    Sym newtree = subtree_builder.build_tree(6, true); // TODO, parameterize
 	    static_cast<Sym&>(tomutate) = insert_subtree(tomutate, xover_point, newtree);
diff --git a/eo/contrib/mathsym/fun/FunDef.cpp b/eo/contrib/mathsym/fun/FunDef.cpp
index 1cbffd09..7c11010b 100644
--- a/eo/contrib/mathsym/fun/FunDef.cpp
+++ b/eo/contrib/mathsym/fun/FunDef.cpp
@@ -127,13 +127,12 @@ typedef hash_map<double, token_t, HashDouble> DoubleSet;
 #endif
 
 static DoubleSet doubleSet; // for quick checking if a constant already exists
-static vector<bool> is_constant_flag;
 static vector<double> token_value;
 
 static std::vector<token_t> free_list;
 
-void delete_val(token_t token) { // clean up the information about this constant
-    if (is_constant_flag[token]) {
+void delete_val(token_t token) { // clean up the information about this value
+    if (is_constant(token)) {
 	double value = token_value[token];
 	
 	delete language[token];
@@ -142,9 +141,26 @@ void delete_val(token_t token) { // clean up the information about this constant
 	doubleSet.erase(value);
 	free_list.push_back(token);
     }
+
+    if (is_lambda(token)) {
+	delete language[token];
+	language[token] = 0;
+	free_list.push_back(token);
+    }
 }
+
+
 FunDef* make_const(double value);
+
+void extend_free_list() {
+    unsigned sz = language.size();
+    language.resize(sz + sz+1); // double
     
+    for (unsigned i = sz; i < language.size(); ++i) {
+       free_list.push_back(i); 
+    }
+}
+
 Sym SymConst(double value) { 
    
     Sym::set_extra_dtor(delete_val);
@@ -157,14 +173,8 @@ Sym SymConst(double value) {
        
     
     if (free_list.empty()) { // make space for tokens;
-	unsigned sz = language.size();
-	language.resize(sz + sz+1); // double
-	is_constant_flag.resize(language.size(), false);
+	extend_free_list();
 	token_value.resize(language.size(), 0.0);
-	
-	for (unsigned i = sz; i < language.size(); ++i) {
-	   free_list.push_back(i); 
-	}
     }
 
     token_t token = free_list.back();
@@ -174,17 +184,11 @@ Sym SymConst(double value) {
     
     language[token] = make_const(value);
     doubleSet[value] = token;
-    is_constant_flag[token] = true;
     token_value[token] = value;
     
     return Sym(token);
 }
 
-bool is_constant(token_t token) {
-    if (token >= is_constant_flag.size()) return false;
-    return is_constant_flag[token];
-}
-
 /* LanguageTable depends on this one, XXX move somewhere safe.*/
 #include <utils/eoRNG.h>
 extern Sym default_const() { return SymConst(rng.normal()); }
@@ -252,11 +256,12 @@ class Const : public FunDef {
 
 	string name() const { return "parameter"; }
 };
+
 } // namespace 
     
 void get_constants(Sym sym, vector<double>& ret) {
     token_t token = sym.token();
-    if (is_constant_flag[token]) {
+    if (is_constant(token)) {
 	double val = static_cast<const Const*>(language[token])->value;
 	ret.push_back(val);
     }
@@ -280,7 +285,7 @@ vector<double> get_constants(Sym sym) {
 Sym set_constants(Sym sym, vector<double>::const_iterator& it) {
     
     token_t token = sym.token();
-    if (is_constant_flag[token]) {
+    if (is_constant(token)) {
 	return SymConst(*it++);
     }
 
@@ -303,7 +308,7 @@ vector<const FunDef*> get_defined_functions() {
     for (unsigned i = 0; i < language.size(); ++i) {
 	res.push_back(language[i]);
 
-	if (dynamic_cast<const Const*>(language[i]) != 0 || (dynamic_cast<const Var*>(language[i]) != 0) ) {
+	if (is_constant(i) || is_variable(i)) {
 	    res.back() = 0; // erase
 	}
     }
@@ -314,6 +319,103 @@ vector<const FunDef*> get_defined_functions() {
 FunDef* make_var(int idx) { return new Var(idx); }
 FunDef* make_const(double value) { return new Const(value); }
 
+bool is_constant(token_t token) {
+    const Const* cnst = dynamic_cast<const Const*>( language[token] );
+    return cnst != 0;
+}
+
+bool is_variable(token_t token) {
+    const Var* var = dynamic_cast<const Var*>( language[token] );
+    return var != 0;
+}
+
+
+namespace {
+class Lambda : public FunDef {
+    public:
+	Sym expression;
+	int arity;
+	
+	Lambda(Sym expr, int arity_) : expression(expr), arity(arity_) {}
+    
+	double eval(const vector<double>& vals, const vector<double>& _) const {
+	    return ::eval(expression, vals); 
+	}
+	
+	string c_print(const vector<string>& args, const vector<string>& str) const {
+	    return ::c_print(expression, args);
+	}
+	
+	Interval eval(const vector<Interval>& args, const vector<Interval>& inputs) const { 
+	    return ::eval(expression, args);
+	}
+
+	unsigned min_arity() const { return arity; }
+	string name() const { return "F"; }
+	
+};
+    Sym normalize(Sym sym, SymVec& args) {
+	// check if it's a variable
+	token_t token = sym.token();
+	const Var* var = dynamic_cast< const Var*>(language[token]);
+    
+	if (var != 0) {
+	    for (unsigned i = 0; i < args.size(); ++i) {
+		if (sym == args[i]) {
+		    return SymVar(i); // replace with reference to arg
+		}
+	    }
+
+	    // not replaced, add it
+	    args.push_back(sym);
+	    return SymVar(args.size()-1);
+	    
+	}
+	
+	SymVec a = sym.args();
+	for (unsigned i = 0; i < a.size(); ++i) {
+	    a[i] = normalize(a[i], args);
+	}
+	
+	return Sym(token, a);
+    }
+}
+
+
+bool is_lambda(token_t token) {
+    const Lambda* lambda = dynamic_cast<const Lambda*>( language[token]); 
+    return lambda != 0;
+}
+
+
+Sym SymLambda(Sym expression) {
+    vector<Sym> args;
+    Sym expr = normalize(expression, args);
+    
+    // check if expression is already present as a lambda expression
+    for (unsigned i = 0; i < language.size(); ++i) {
+	const Lambda* lambda = dynamic_cast<const Lambda*>(language[i]);
+	if (lambda != 0 && lambda->expression == expr) {
+	    return Sym(i, args);
+	}
+    }
+    // else add it
+    Lambda* lambda = new Lambda(expr, args.size());
+    
+    
+    // insert in language table
+    if (free_list.empty()) {
+	extend_free_list();
+    }
+    
+    token_t lambda_token = free_list.back();
+    free_list.pop_back();
+    language[lambda_token] = lambda;
+    
+    
+    return Sym(lambda_token, args);
+}
+
 
 namespace {
 
@@ -550,7 +652,7 @@ void write_raw(ostream& os, const Sym& sym) {
     token_t token = sym.token();
     const SymVec& args = sym.args();
     
-    if (is_constant_flag.size() > token && is_constant_flag[token]) {
+    if (is_constant(token)) {
 	os << "c" << language[token]->c_print(vector<string>(), vector<string>()); 
     } else {
 
diff --git a/eo/contrib/mathsym/fun/FunDef.h b/eo/contrib/mathsym/fun/FunDef.h
index 6f517602..bbb7f39e 100644
--- a/eo/contrib/mathsym/fun/FunDef.h
+++ b/eo/contrib/mathsym/fun/FunDef.h
@@ -52,7 +52,10 @@ class FunDef {
     
 };
 
+/** Gets out all function that are defined (excluding constants and variables) */
 extern std::vector<const FunDef*> get_defined_functions();
+
+/** Gets a specific function (including vars and constants) out */
 extern const FunDef& get_element(token_t token);
 
 /** Single case evaluation */
@@ -70,8 +73,24 @@ extern std::string c_print(const Sym& sym, const std::vector<std::string>& var_n
 /** Pretty printer streamer */
 inline std::ostream& operator<<(std::ostream& os, Sym sym) { return os << c_print(sym); }
 
+/* Support for Ephemeral Random Constants (ERC) */
+
 /** Create constant with this value, memory is managed. If reference count drops to zero value is deleted.  */
 extern Sym SymConst(double value);
+/** Create variable */
+extern Sym SymVar(unsigned idx);
+
+/** Create 'lambda expression; 
+ * This is a neutral operation. It will replace
+ * all variables in the expression by arguments,
+ * wrap the expression in a Lambda function
+ * and returns a tree applying the lambda function
+ * to the original variable.
+ *
+ * A call like SymLambda( SymVar(1) + SymVar(1) * 3.1) will result in 
+ * a Lambda function (a0 + a0 * 3.1) with one argument: SymVar(1)*/
+
+extern Sym SymLambda(Sym expression);
 
 /** Get out the values for all constants in the expression */
 std::vector<double> get_constants(Sym sym);
@@ -80,13 +99,15 @@ std::vector<double> get_constants(Sym sym);
  * The argument isn't touched, it will return a new sym with the constants set. */
 Sym set_constants(Sym sym, const std::vector<double>& constants);
 
-/** check if a token is constant (for instance Sym sym; is_constant(sym.token()); ) */
+/** check if a token is a constant */
 extern bool is_constant(token_t token);
+/** check if a token is a variable */
+extern bool is_variable(token_t token);
+/** check if a token is a user/automatically defined function */
+extern bool is_lambda(token_t token);
 
-/** Create variable */
-extern Sym SymVar(unsigned idx);
 
-/** simplifies a sym (sym_operations.cpp) */
+/** simplifies a sym (sym_operations.cpp) Currently only simplifies constants */
 extern Sym simplify(Sym sym);
 
 /** differentiates a sym to a token (sym_operations.cpp)
@@ -114,6 +135,7 @@ enum {
     sqr_token, sqrt_token
 };
 
+/* Defition of function overloads: for example, this define the function 'Sym sin(Sym)' */
 
 #define HEADERFUNC(name) inline Sym name(Sym arg) { return Sym(name##_token, arg); }
 
@@ -139,9 +161,10 @@ HEADERFUNC(log);
 HEADERFUNC(sqr);
 HEADERFUNC(sqrt);
 
+/* Get the prototype functions out, this is needed for compilation */
 extern std::string get_prototypes();
 
-// reading and writing in internal format
+// reading and writing in internal format, no parser for symbolic functions implemented yet
 extern std::string write_raw(const Sym& sym);
 extern void write_raw(std::ostream& os, const Sym& sym);
 extern Sym read_raw(std::string str);
diff --git a/eo/contrib/mathsym/gen/NodeSelector.cpp b/eo/contrib/mathsym/gen/NodeSelector.cpp
index db41448f..aefa4dd6 100644
--- a/eo/contrib/mathsym/gen/NodeSelector.cpp
+++ b/eo/contrib/mathsym/gen/NodeSelector.cpp
@@ -20,21 +20,31 @@
 
 #include <utils/eoRNG.h>
 
-unsigned RandomNodeSelector::select_node(Sym sym) const {
-    return rng.random(sym.size());
+// If subtree is not set (by randomnodeselector for instance, find it now
+Sym NodeSelector::NodeSelection::subtree() {
+    if (subtree_.empty()) {
+	subtree_ = get_subtree(root_, subtree_index_);
+    }
+    return subtree_;
 }
 
-unsigned BiasedNodeSelector::select_node(Sym sym) const {
+NodeSelector::NodeSelection RandomNodeSelector::select_node(Sym sym) const {
+    unsigned idx = rng.random(sym.size());
+    return  NodeSelection(sym, idx, Sym() ); // empty subtree, find it when needed
+}
+
+NodeSelector::NodeSelection BiasedNodeSelector::select_node(Sym sym) const {
     
     unsigned p = rng.random(sym.size());
+    Sym res;
     for (unsigned i = 0; i < nRounds; ++i) {
-	Sym res = get_subtree(sym, p);
+	res = get_subtree(sym, p);
 	
 	if (res.args().size() > 0) break;
 	
 	p = rng.random(sym.size());
     }
     
-    return p;
+    return NodeSelection(sym, p, res);
 }
 
diff --git a/eo/contrib/mathsym/gen/NodeSelector.h b/eo/contrib/mathsym/gen/NodeSelector.h
index 5a93575b..e30abd9f 100644
--- a/eo/contrib/mathsym/gen/NodeSelector.h
+++ b/eo/contrib/mathsym/gen/NodeSelector.h
@@ -18,21 +18,37 @@
 #ifndef NODESELECTOR_H
 #define NODESELECTOR_H
 
-class Sym;
+#include <Sym.h>
 
 /** Base class for selecting nodes */
 class NodeSelector {
     public:
+	
+    class NodeSelection {
+	Sym root_;
+	unsigned subtree_index_;
+	Sym subtree_;
+
+	public :
+	    NodeSelection(Sym r, unsigned idx, Sym s) 
+		: root_(r), subtree_index_(idx), subtree_(s) {}
+
+	    Sym root() const { return root_; }
+	    unsigned idx()  const { return subtree_index_; }
+	    Sym subtree(); 
+	
+    };
+
     virtual ~NodeSelector() {}
 	
-    virtual unsigned select_node(Sym sym) const = 0;  
+    virtual NodeSelection select_node(Sym sym) const = 0;  
 };
 
 
 /** Select nodes uniformly */
 class RandomNodeSelector : public NodeSelector {
     public:
-    unsigned select_node(Sym sym) const;
+    NodeSelection select_node(Sym sym) const;
 };
 
 /** A node selector that does a specified number of rounds ignoring terminals */
@@ -43,7 +59,7 @@ class BiasedNodeSelector : public NodeSelector {
     BiasedNodeSelector() : nRounds(3) {} // 3: for binary trees 87.5% chance of selecting an internal node
     BiasedNodeSelector(unsigned n) : nRounds(n) {}
     
-    unsigned select_node(Sym sym) const;
+    NodeSelection select_node(Sym sym) const;
 };
 
 #endif
diff --git a/eo/contrib/mathsym/symreg.cpp b/eo/contrib/mathsym/symreg.cpp
index 406546ea..5254f386 100644
--- a/eo/contrib/mathsym/symreg.cpp
+++ b/eo/contrib/mathsym/symreg.cpp
@@ -100,7 +100,7 @@ class DagSizeStat : public eoStat<EoType, unsigned> {
 };
 
 int main(int argc, char* argv[]) {
-    
+   
     eoParser parser(argc, argv);
   
     /* Language */
diff --git a/eo/contrib/mathsym/test/test_lambda.cpp b/eo/contrib/mathsym/test/test_lambda.cpp
new file mode 100644
index 00000000..dc8e1ac7
--- /dev/null
+++ b/eo/contrib/mathsym/test/test_lambda.cpp
@@ -0,0 +1,34 @@
+#include <FunDef.h>
+
+using namespace std;
+
+int main() {
+
+    Sym x = SymVar(0);
+    Sym y = SymVar(1);
+
+    Sym f = y + x*x;
+
+    Sym l = SymLambda(f);
+    
+    SymVec args = l.args();
+    args[0] = x;
+    args[1] = y;
+    l = Sym(l.token(), args);
+    
+    vector<double> v(3);
+    v[0] = 2.0;
+    v[1] = 3.0;
+    v[2] = 4.0;
+
+    double v1 = eval(f,v);
+    double v2 = eval(l,v);
+
+    cout << v1 << ' ' << v2 << endl;
+    cout << f << endl;
+    cout << l << endl;
+
+    if (v1 != 7.0) return 1;
+    if (v2 != 11.0) return 1;
+}
+