initial experimentation with genetic programming

9 files changed, 921 insertions, 5 deletions

diff --git a/eval/CMakeLists.txt b/eval/CMakeLists.txt
index d107ebfed40..f063faf19c3 100644
--- a/eval/CMakeLists.txt
+++ b/eval/CMakeLists.txt
@@ -23,11 +23,12 @@ vespa_define_module(
     src/tests/eval/tensor_spec
     src/tests/eval/value_cache
     src/tests/eval/value_type
+    src/tests/gp/ponder_nov2017
     src/tests/tensor/dense_dot_product_function
-    src/tests/tensor/dense_xw_product_function
     src/tests/tensor/dense_tensor_address_combiner
     src/tests/tensor/dense_tensor_builder
     src/tests/tensor/dense_tensor_function_compiler
+    src/tests/tensor/dense_xw_product_function
     src/tests/tensor/sparse_tensor_builder
     src/tests/tensor/tensor_address
     src/tests/tensor/tensor_conformance
@@ -39,11 +40,12 @@ vespa_define_module(
     LIBS
     src/vespa/eval
     src/vespa/eval/eval
+    src/vespa/eval/eval/llvm
     src/vespa/eval/eval/test
     src/vespa/eval/eval/value_cache
-    src/vespa/eval/eval/llvm
+    src/vespa/eval/gp
     src/vespa/eval/tensor
-    src/vespa/eval/tensor/sparse
     src/vespa/eval/tensor/dense
     src/vespa/eval/tensor/serialization
+    src/vespa/eval/tensor/sparse
 )
diff --git a/eval/src/tests/gp/ponder_nov2017/.gitignore b/eval/src/tests/gp/ponder_nov2017/.gitignore
new file mode 100644
index 00000000000..913fb857d36
--- /dev/null
+++ b/eval/src/tests/gp/ponder_nov2017/.gitignore
@@ -0,0 +1 @@
+/eval_gp_ponder_nov2017_app
diff --git a/eval/src/tests/gp/ponder_nov2017/CMakeLists.txt b/eval/src/tests/gp/ponder_nov2017/CMakeLists.txt
new file mode 100644
index 00000000000..7c915b46de7
--- /dev/null
+++ b/eval/src/tests/gp/ponder_nov2017/CMakeLists.txt
@@ -0,0 +1,7 @@
+# Copyright 2017 Yahoo Holdings. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+vespa_add_executable(eval_gp_ponder_nov2017_app
+    SOURCES
+    ponder_nov2017.cpp
+    DEPENDS
+    vespaeval
+)
diff --git a/eval/src/tests/gp/ponder_nov2017/ponder_nov2017.cpp b/eval/src/tests/gp/ponder_nov2017/ponder_nov2017.cpp
new file mode 100644
index 00000000000..c1af2f2b74c
--- /dev/null
+++ b/eval/src/tests/gp/ponder_nov2017/ponder_nov2017.cpp
@@ -0,0 +1,285 @@
+// Copyright 2017 Yahoo Holdings. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+
+#include <vespa/vespalib/testkit/test_kit.h>
+#include <vespa/vespalib/util/signalhandler.h>
+#include <vespa/eval/gp/gp.h>
+#include <limits.h>
+#include <algorithm>
+
+using namespace vespalib;
+using namespace vespalib::gp;
+
+// Inspired by the great and sometimes frustrating puzzles posed to us
+// by IBM; what about automatically evolving a solution instead of
+// figuring it out on our own. Turns out GP is no free lunch, but
+// rather a strange and interesting adventure all of its own...
+
+// problem: https://www.research.ibm.com/haifa/ponderthis/challenges/November2017.html
+// solution: https://www.research.ibm.com/haifa/ponderthis/solutions/November2017.html
+
+// illegal div/mod will result in 0
+bool div_ok(int a, int b) {
+    if ((a == INT_MIN) && (b == -1)) {
+        return false;
+    }
+    return (b != 0);
+}
+int my_add(int a, int b) { return a + b; }
+int my_sub(int a, int b) { return a - b; }
+int my_mul(int a, int b) { return a * b; }
+int my_div(int a, int b) { return div_ok(a, b) ? (a / b) : 0; }
+int my_mod(int a, int b) { return div_ok(a, b) ? (a % b) : 0; }
+int my_pow(int a, int b) { return pow(a,b); }
+int my_and(int a, int b) { return a & b; }
+int my_or(int a, int b)  { return a | b; }
+int my_xor(int a, int b) { return a ^ b; }
+
+struct Dist {
+    std::vector<int> slots;
+    static size_t need_slots(size_t num_outputs) {
+        size_t result = 6; // z
+        if (num_outputs > 1) {
+            result *= 2; // y
+            if (num_outputs > 2) {
+                result *= 2; // x
+                ASSERT_EQUAL(num_outputs, 3u);
+            }
+        }
+        return result;
+    }
+    Dist(size_t num_outputs) : slots(need_slots(num_outputs), 0) {}
+    void sample(int z) {
+        int post_z = (size_t(z) % 6);
+        ASSERT_GREATER_EQUAL(post_z, 0);
+        ASSERT_LESS(post_z, 6);
+        int slot = post_z;
+        ASSERT_LESS(size_t(slot), slots.size());
+        ++slots[slot];
+    }
+    void sample(int z, int y) {
+        int post_y = (y & 1);
+        int post_z = (size_t(z) % 6);
+        ASSERT_GREATER_EQUAL(post_y, 0);
+        ASSERT_GREATER_EQUAL(post_z, 0);
+        ASSERT_LESS(post_y, 2);
+        ASSERT_LESS(post_z, 6);
+        int slot = (post_z<<1) | (post_y);
+        ASSERT_LESS(size_t(slot), slots.size());
+        ++slots[slot];
+    }
+    void sample(int z, int y, int x) {
+        int post_x = (x & 1);
+        int post_y = (y & 1);
+        int post_z = (size_t(z) % 6);
+        ASSERT_GREATER_EQUAL(post_x, 0);
+        ASSERT_GREATER_EQUAL(post_y, 0);
+        ASSERT_GREATER_EQUAL(post_z, 0);
+        ASSERT_LESS(post_x, 2);
+        ASSERT_LESS(post_y, 2);
+        ASSERT_LESS(post_z, 6);
+        int slot = (post_z<<2) | (post_y<<1) | (post_x);
+        ASSERT_LESS(size_t(slot), slots.size());
+        ++slots[slot];
+    }
+    size_t error() const {
+        size_t err = 0;
+        int expect = (216 / slots.size());
+        ASSERT_EQUAL(216 % slots.size(), 0u);
+        for (int cnt: slots) {
+            err += (std::max(cnt, expect) - std::min(cnt, expect));
+        }
+        return err;
+    }
+};
+
+Feedback find_weakness(const MultiFunction &fun) {
+    size_t num_outputs = fun.num_outputs();
+    std::vector<Dist> state(fun.num_alternatives(), Dist(num_outputs));
+    for (int d1 = 1; d1 <= 6; ++d1) {
+        for (int d2 = 1; d2 <= 6; ++d2) {
+            for (int d3 = 1; d3 <= 6; ++d3) {
+                Input input({d1, d2, d3});
+                std::sort(input.begin(), input.end());
+                if (fun.num_inputs() == 6) {
+                    // add const values for hand-crafted case
+                    input.push_back(2);
+                    input.push_back(1502);
+                    input.push_back(70677);
+                }
+                Result result = fun.execute(input);
+                ASSERT_EQUAL(result.size(), state.size());
+                for (size_t i = 0; i < result.size(); ++i) {
+                    const Output &output = result[i];
+                    switch(output.size()) {
+                    case 1:
+                        state[i].sample(output[0]); // z
+                        break;
+                    case 2:
+                        state[i].sample(output[0], output[1]); // z,y
+                        break;
+                    default:
+                        ASSERT_EQUAL(output.size(), 3u);
+                        state[i].sample(output[0], output[1], output[2]); // z,y,x
+                    }
+                }
+            }
+        }
+    }
+    Feedback feedback;
+    for (const Dist &dist: state) {
+        feedback.push_back(dist.error());
+    }
+    return feedback;
+}
+
+OpRepo my_repo() {
+    return OpRepo(find_weakness)
+        .add("add", my_add)  // 1
+        .add("sub", my_sub)  // 2
+        .add("mul", my_mul)  // 3
+        .add("div", my_div)  // 4
+        .add("mod", my_mod)  // 5
+        .add("pow", my_pow)  // 6
+        .add("and", my_and)  // 7
+        .add("or",  my_or)   // 8
+        .add("xor", my_xor); // 9
+}
+
+// Featured solution (Bert Dobbelaere):
+//
+// d=2**(((c-a)*(c+a))/2)
+//     x=(1502/d)%2
+//     y=(70677/d)%2
+//     z=(a+b+c)%6+1
+
+const size_t add_id = 1;
+const size_t sub_id = 2;
+const size_t mul_id = 3;
+const size_t div_id = 4;
+const size_t pow_id = 6;
+
+using Ref = Program::Ref;
+using Op = Program::Op;
+
+TEST("evaluating hand-crafted solution") {
+    // constants are modeled as inputs
+    Program prog(my_repo(), 6, 3, 2, 0);
+    auto a = Ref::in(0);                   // a
+    auto b = Ref::in(1);                   // b
+    auto c = Ref::in(2);                   // c
+    auto k1 = Ref::in(3);                  // 2
+    auto k2 = Ref::in(4);                  // 1502
+    auto k3 = Ref::in(5);                  // 70677
+    auto _1 = prog.add_op(sub_id, c, a);   // _1 = c-a
+    auto _2 = prog.add_op(add_id, c, a);   // _2 = c+a
+    auto _3 = prog.add_op(mul_id, _1, _2); // _3 = (c-a)*(c+a)
+    // (zero-cost forwarding, for testing)
+    _1 = prog.add_forward(_1);
+    _2 = prog.add_forward(_2);
+    _3 = prog.add_forward(_3);
+    auto _4 = prog.add_op(div_id, _3, k1); // _4 = ((c-a)*(c+a))/2
+    auto d = prog.add_op(pow_id, k1, _4);  // d = 2**(((c-a)*(c+a))/2)
+    auto _5 = prog.add_op(add_id, a, b);   // _5 = a+b
+    // --- alt 0 (dummy outputs, for testing)
+    prog.add_forward(_1);
+    prog.add_forward(_2);
+    prog.add_forward(_3);
+    // --- alt 1 (correct output)
+    auto z = prog.add_op(add_id, _5, c);   // z = a+b+c
+    auto y = prog.add_op(div_id, k3, d);   // y = 70677/d
+    auto x = prog.add_op(div_id, k2, d);   // x = 1502/d
+    // '%2' (for x and y) and '%6+1' (for z) done outside program
+    //--- verify sub-expressions
+    EXPECT_EQUAL(prog.as_string(a), "i0");
+    EXPECT_EQUAL(prog.as_string(k2), "i4");
+    EXPECT_EQUAL(prog.as_string(d), "pow(i3,div(mul(sub(i2,i0),add(i2,i0)),i3))");
+    EXPECT_EQUAL(prog.as_string(x), "div(i4,pow(i3,div(mul(sub(i2,i0),add(i2,i0)),i3)))");
+    EXPECT_EQUAL(prog.as_string(y), "div(i5,pow(i3,div(mul(sub(i2,i0),add(i2,i0)),i3)))");
+    EXPECT_EQUAL(prog.as_string(z), "add(add(i0,i1),i2)");
+    //--- verify (expression) sizes
+    EXPECT_EQUAL(prog.size_of(a), 1u);
+    EXPECT_EQUAL(prog.size_of(k2), 1u);
+    EXPECT_EQUAL(prog.size_of(d), 11u);
+    EXPECT_EQUAL(prog.size_of(x), 13u);
+    EXPECT_EQUAL(prog.size_of(y), 13u);
+    EXPECT_EQUAL(prog.size_of(z), 5u);
+    //--- verify costs
+    EXPECT_EQUAL(prog.get_cost(0), 3u);
+    EXPECT_EQUAL(prog.get_cost(1), 9u);
+    //--- evaluate
+    Random dummy;
+    prog.handle_feedback(dummy, find_weakness(prog));
+    EXPECT_EQUAL(prog.stats().weakness, 0.0);
+    EXPECT_EQUAL(prog.stats().cost, 9u);
+    EXPECT_EQUAL(prog.stats().alt, 1u);
+}
+
+void maybe_newline(bool &partial_line) {
+    if (partial_line) {
+        fprintf(stderr, "\n");
+        partial_line = false;
+    }
+}
+
+Program try_evolve(const Params &params, size_t max_idle, const Program *program = nullptr) {
+    Population population(params, my_repo(), Random().make_seed());
+    if (program != nullptr) {
+        population.init(*program);
+    }
+    bool partial_line = false;
+    size_t ticks = 0;
+    size_t sample_tick = ticks;
+    Program::Stats best_sample = population._programs[0].stats();
+    while (!SignalHandler::INT.check() &&
+           ((best_sample.weakness > 0) ||
+            ((ticks - sample_tick) < max_idle)))
+    {
+        ++ticks;
+        population.tick();
+        if ((ticks % 500) == 0) {
+            maybe_newline(partial_line);
+            population.print_stats();
+        } else if ((ticks % 10) == 0) {
+            fprintf(stderr, ".");
+            partial_line = true;
+        }
+        Program::Stats sample = population._programs[0].stats();
+        best_sample.born = sample.born;
+        if (sample < best_sample) {
+            best_sample = sample;
+            sample_tick = ticks;
+        }
+    }
+    if (SignalHandler::INT.check()) {
+        fprintf(stderr, "<INT>\n");
+        SignalHandler::INT.clear();
+    }
+    maybe_newline(partial_line);
+    Program::Stats best = population._programs[0].stats();
+    fprintf(stderr, "best stats after %zu ticks: (weakness=%g,cost=%zu)\n",
+            ticks, best.weakness, best.cost);
+    return population._programs[0];
+}
+
+// best stats: (weakness=0,cost=9)
+// x(size=21): mod(add(div(add(i2,i0),i0),and(mod(mul(i1,add(i1,add(i2,i0))),add(i2,i0)),i2)),i2)
+// y(size=13): sub(mod(mul(i1,add(i1,add(i2,i0))),add(i2,i0)),i2)
+// z(size=5): add(i1,add(i2,i0))
+
+TEST("trying to evolve a solution automatically") {
+    fprintf(stderr, "training f(a,b,c) -> (z)...\n");
+    Program best_z = try_evolve(Params(3, 1, 8, 8, 8), 10 * 1000);
+    fprintf(stderr, "training f(a,b,c) -> (z,y)...\n");
+    Program best_zy = try_evolve(Params(3, 2, 8, 8, 8), 100 * 1000, &best_z);
+    fprintf(stderr, "training f(a,b,c) -> (z,y,x)...\n");
+    Program best = try_evolve(Params(3, 3, 8, 8, 8), 1000 * 1000 * 1000, &best_zy);
+    auto refs = best.get_refs(best.stats().alt);
+    fprintf(stderr, "x(size=%zu): %s\n", best.size_of(refs[2]), best.as_string(refs[2]).c_str());
+    fprintf(stderr, "y(size=%zu): %s\n", best.size_of(refs[1]), best.as_string(refs[1]).c_str());
+    fprintf(stderr, "z(size=%zu): %s\n", best.size_of(refs[0]), best.as_string(refs[0]).c_str());
+}
+
+TEST_MAIN() {
+    SignalHandler::INT.hook();
+    TEST_RUN_ALL();
+}
diff --git a/eval/src/vespa/eval/CMakeLists.txt b/eval/src/vespa/eval/CMakeLists.txt
index 8a2cf16bf87..94a46015b7d 100644
--- a/eval/src/vespa/eval/CMakeLists.txt
+++ b/eval/src/vespa/eval/CMakeLists.txt
@@ -2,13 +2,14 @@
 vespa_add_library(vespaeval
     SOURCES
     $<TARGET_OBJECTS:eval_eval>
+    $<TARGET_OBJECTS:eval_eval_llvm>
     $<TARGET_OBJECTS:eval_eval_test>
     $<TARGET_OBJECTS:eval_eval_value_cache>
-    $<TARGET_OBJECTS:eval_eval_llvm>
+    $<TARGET_OBJECTS:eval_gp>
     $<TARGET_OBJECTS:eval_tensor>
-    $<TARGET_OBJECTS:eval_tensor_sparse>
     $<TARGET_OBJECTS:eval_tensor_dense>
     $<TARGET_OBJECTS:eval_tensor_serialization>
+    $<TARGET_OBJECTS:eval_tensor_sparse>
     INSTALL lib64
     DEPENDS
     LLVM-${VESPA_LLVM_VERSION}
diff --git a/eval/src/vespa/eval/gp/CMakeLists.txt b/eval/src/vespa/eval/gp/CMakeLists.txt
new file mode 100644
index 00000000000..ea61ca3dd27
--- /dev/null
+++ b/eval/src/vespa/eval/gp/CMakeLists.txt
@@ -0,0 +1,5 @@
+# Copyright 2017 Yahoo Holdings. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+vespa_add_library(eval_gp OBJECT
+    SOURCES
+    gp.cpp
+)
diff --git a/eval/src/vespa/eval/gp/README b/eval/src/vespa/eval/gp/README
new file mode 100644
index 00000000000..97bce8d0935
--- /dev/null
+++ b/eval/src/vespa/eval/gp/README
@@ -0,0 +1,53 @@
+Experimental code related to genetic programming.
+
+The goal is to evolve a function satisfying predefined constraints as
+closely as possible with minimal computational effort. We are willing
+to spend more CPU finding a function taking less CPU to execute.
+
+Possible solutions are represented by simulated individuals that
+evolve over time. Each individual has a weakness value indicating how
+well it solves the problem. A weakness of 0 indicates that the
+individual solves the problem perfectly. Worse solutions have higher
+weakness. Each individual also has information about the cost of its
+solution as well as the age of the individual itself.
+
+Simulated individuals are grouped together in a population. The
+triplet (weakness,cost,age) enables us to impose a strict weak
+ordering of individuals within a population. This ordering is used to
+skew the selection of which individuals to further evolve in order to
+find the best possible solution. We prefer young solutions with low
+cost and little weakness.
+
+Each simulated individual may represent multiple possible solutions
+for the function we are trying to evolve. The fitness evaluation of an
+individual will evaluate the weakness of each possible solution
+represented by the individual and give direct feedback to the
+individual by presenting it with the weakness value for each
+alternative solution. This enables the individual to update the
+weakness and cost values used to order it within the population. It
+also enables some self-reflection about the less awesome solutions it
+contains.
+
+Having multiple possible solutions represented by a single individual
+enables even further sorting criteria to be imposed within a
+population based on the potential quality of unexpressed
+solutions. This concept is inverted and called 'waste'.
+
+When training a solution with multiple outputs, it might make sense to
+start training a single output and add more outputs as we go. To focus
+the training on the newly added outputs and how they can best re-use
+the calculations needed by the already trained outputs the code used
+to produce the other outputs is frozen during training. To avoid the
+quality of training to depend on the order which outputs are trained,
+there needs to be more advanced ways to freeze/thaw operations and
+outputs. Maybe we need to do circular freeze/thaw to simulate
+multi-dimensional gradient decent. Maybe we need to utilize multiple
+parallel generations with different free outputs. Maybe freeze/thaw
+itself is a mutation.
+
+Another theme that should be investigated is the cross-over to
+reinforcement learning and intelligent design. Can we perform error
+push-back through a genetically learned program? Can we efficiently
+express neural networks in the selected programming model? Can we
+perform smart mutation with greater than random chance of improving
+the solution?
diff --git a/eval/src/vespa/eval/gp/gp.cpp b/eval/src/vespa/eval/gp/gp.cpp
new file mode 100644
index 00000000000..5ef9ea8f004
--- /dev/null
+++ b/eval/src/vespa/eval/gp/gp.cpp
@@ -0,0 +1,331 @@
+// Copyright 2017 Yahoo Holdings. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+
+#include "gp.h"
+#include <algorithm>
+#include <vespa/vespalib/util/stringfmt.h>
+#include <map>
+
+namespace vespalib::gp {
+
+namespace {
+
+Value get(const Input &input, const std::vector<Value> &values, Program::Ref ref) {
+    return ref.is_input() ? input[ref.in_idx()] : values[ref.op_idx()];
+}
+
+size_t get(const std::vector<size_t> &sizes, Program::Ref ref) {
+    return ref.is_input() ? 1 : sizes[ref.op_idx()];
+}
+
+Program::Ref map(const std::map<Program::Ref,Program::Ref> &ref_map, Program::Ref ref) {
+    if (ref.is_input()) {
+        return ref;
+    }
+    auto pos = ref_map.find(ref);
+    assert(pos != ref_map.end());
+    return pos->second;
+}
+
+} // namespace vespalib::gp::<unnamed>
+
+void
+Program::assert_valid(Ref ref, size_t limit) const
+{
+    assert(ref.is_input() != ref.is_operation());
+    if (ref.is_input()) {
+        assert(ref.in_idx() < _in_cnt);
+    }
+    if (ref.is_operation()) {
+        assert(ref.op_idx() < limit);
+    }
+}
+
+Program::Ref
+Program::add_op(size_t code, Ref lhs, Ref rhs)
+{
+    size_t op_idx = _program.size();
+    assert(code <= _repo.max_op());
+    assert_valid(lhs, op_idx);
+    assert_valid(rhs, op_idx);
+    _program.emplace_back(code, lhs, rhs);
+    return Ref::op(op_idx);
+}
+
+Program::Ref
+Program::add_forward(Ref ref)
+{
+    return add_op(0, ref, Ref::nop());
+}
+
+void
+Program::init(const Program &src)
+{
+    assert(src._out_cnt < _out_cnt);
+    std::map<Ref,Ref> ref_map;
+    auto used = src.get_used_ops(src.stats().alt);
+    for (size_t i = 0; i < used.size(); ++i) {
+        if (used[i]) {
+            const Op &op = src._program[i];
+            if (op.code == 0) { // forward
+                auto res = ref_map.emplace(Ref::op(i), map(ref_map, op.lhs));
+                assert(res.second);
+            } else {
+                auto res = ref_map.emplace(Ref::op(i), Ref::op(_program.size()));
+                assert(res.second);
+                _program.emplace_back(op.code,
+                                      map(ref_map, op.lhs),
+                                      map(ref_map, op.rhs));
+            }
+        }
+    }
+    _frozen = _program.size();
+    for (Ref ref: src.get_refs(src.stats().alt)) {
+        _bound.push_back(map(ref_map, ref));
+    }
+}
+
+void
+Program::grow(Random &rnd, size_t op_cnt)
+{
+    for (size_t i = 0; i < op_cnt; ++i) {
+        size_t op_idx = _program.size();
+        add_op(rnd_op(rnd),
+               rnd_ref(rnd, op_idx),
+               rnd_ref(rnd, op_idx));
+    }
+    size_t prefix = _program.size();
+    size_t suffix = _alt_cnt * get_alt_size();
+    for (size_t i = 0; i < suffix; ++i) {
+        add_op(rnd_op(rnd),
+               rnd_ref(rnd, prefix),
+               rnd_ref(rnd, prefix));
+    }
+}
+
+void
+Program::mutate(Random &rnd, size_t mut_idx)
+{
+    size_t prefix = get_alt_offset(0);
+    Op &op = _program[mut_idx];
+    size_t sel = rnd.get(0,2);
+    if (sel == 0) {
+        op.code = rnd_op(rnd);
+    } else if (sel == 1) {
+        op.lhs = rnd_ref(rnd, std::min(mut_idx, prefix));
+    } else {
+        assert(sel == 2);
+        op.rhs = rnd_ref(rnd, std::min(mut_idx, prefix));
+    }
+}
+
+void
+Program::mutate(Random &rnd)
+{
+    assert(_frozen < _program.size());
+    mutate(rnd, rnd.get(_frozen, _program.size() - 1));
+}
+
+std::vector<Program::Ref>
+Program::get_refs(size_t alt) const
+{
+    std::vector<Ref> refs;
+    refs.reserve(_out_cnt);
+    refs = _bound;
+    size_t offset = get_alt_offset(alt);
+    while (refs.size() < _out_cnt) {
+        refs.push_back(Ref::op(offset++));
+    }
+    return refs;
+}
+
+std::vector<bool>
+Program::get_used_ops(size_t alt) const
+{
+    std::vector<bool> used(_program.size(), false);
+    std::vector<Ref> todo = get_refs(alt);
+    while (!todo.empty()) {
+        Ref ref = todo.back();
+        todo.pop_back();
+        if (ref.is_operation() && !used[ref.op_idx()]) {
+            const Op &op = _program[ref.op_idx()];
+            todo.push_back(op.lhs);
+            if (op.code > 0) {
+                todo.push_back(op.rhs);
+            }
+            used[ref.op_idx()] = true;
+        }
+    }
+    return used;
+}
+
+size_t
+Program::get_cost(size_t alt) const
+{
+    size_t cost = 0;
+    auto used = get_used_ops(alt);
+    for (size_t i = 0; i < used.size(); ++i) {
+        if (used[i]) {
+            cost += _repo.cost_of(_program[i].code);
+        }
+    }
+    return cost;
+}
+
+size_t
+Program::size_of(Ref ref) const
+{
+    assert_valid(ref, _program.size());
+    if (ref.is_input()) {
+        return 1;
+    }
+    std::vector<size_t> sizes;
+    for (size_t i = 0; i <= ref.op_idx(); ++i) {
+        const Op &op = _program[i];
+        if (op.code == 0) {
+            sizes.push_back(get(sizes, op.lhs)); // forward
+        } else {
+            sizes.push_back(1 + get(sizes, op.lhs) + get(sizes, op.rhs));
+        }
+    }
+    return sizes.back();
+}
+
+vespalib::string
+Program::as_string(Ref ref) const
+{
+    assert_valid(ref, _program.size());
+    size_t expr_size = size_of(ref);
+    if (expr_size > 9000) {
+        // its over 9000!
+        return vespalib::make_string("expr(%zu nodes)", expr_size);
+    } else if (ref.is_input()) {
+        return vespalib::make_string("i%zu", ref.in_idx());
+    } else {
+        const Op &my_op = _program[ref.op_idx()];
+        if (my_op.code == 0) {
+            return as_string(my_op.lhs); // forward
+        } else {
+            return vespalib::make_string("%s(%s,%s)", _repo.name_of(my_op.code).c_str(),
+                                         as_string(my_op.lhs).c_str(), as_string(my_op.rhs).c_str());
+        }
+    }
+}
+
+Result
+Program::execute(const Input &input) const
+{
+    Result result;
+    std::vector<Value> values;
+    size_t prefix = get_alt_offset(0);
+    values.reserve(prefix);
+    size_t idx = 0;
+    for (; idx < prefix; ++idx) {
+        const Op &op = _program[idx];
+        values.push_back(_repo.perform(op.code,
+                                       get(input, values, op.lhs),
+                                       get(input, values, op.rhs)));
+    }
+    for (size_t i = 0; i < _alt_cnt; ++i) {
+        std::vector<Value> out;
+        out.reserve(_out_cnt);
+        for (Ref ref: _bound) {
+            out.push_back(get(input, values, ref));
+        }
+        while (out.size() < _out_cnt) {
+            const Op &op = _program[idx++];
+            out.push_back(_repo.perform(op.code,
+                                        get(input, values, op.lhs),
+                                        get(input, values, op.rhs)));
+        }
+        result.push_back(out);
+    }
+    assert(idx == _program.size());
+    return result;
+}
+
+void
+Program::handle_feedback(Random &rnd, const Feedback &feedback)
+{
+    assert(feedback.size() == _alt_cnt);
+    std::vector<Stats> my_stats;
+    my_stats.reserve(_alt_cnt);
+    for (size_t i = 0; i < _alt_cnt; ++i) {
+        my_stats.emplace_back(feedback[i], get_cost(i), _stats.born, i);
+    }
+    std::sort(my_stats.begin(), my_stats.end());
+    _stats = my_stats[0];
+    for (size_t i = 1; i < my_stats.size(); ++i) {
+        if ((i + 1) == my_stats.size()) { // worst
+            size_t len = get_alt_size();
+            size_t src = get_alt_offset(my_stats.front().alt);
+            size_t dst = get_alt_offset(my_stats.back().alt);
+            for (size_t j = 0; j < len; ++j) {
+                _program[dst + j] = _program[src + j];
+            }
+            mutate(rnd, rnd.get(dst, dst + len - 1));
+        } else { // not best, not worst; mediocre
+            double my_waste = (my_stats[i].weakness + 1) *
+                              (my_stats[i].cost + 1);
+            _waste = std::min(my_waste, (i == 1) ? my_waste : _waste);
+        }
+    }
+}
+
+void
+Population::grow(size_t cnt)
+{
+    while (_programs.size() < cnt) {
+        _programs.emplace_back(_repo, _params.in_cnt, _params.out_cnt, _params.alt_cnt, _gen);
+        _programs.back().grow(_rnd, _params.op_cnt);
+        _repo.find_weakness(_rnd, _programs.back());
+    }
+    std::sort(_programs.begin(), _programs.end());
+}
+
+void
+Population::print_stats() const
+{
+    const Program::Stats &best = _programs.front().stats();
+    const Program::Stats &worst = _programs.back().stats();
+    fprintf(stderr, "[%zu] best(weakness=%g,cost=%zu,age=%zu), "
+            "worst(weakness=%g,cost=%zu,age=%zu)\n", _gen,
+            best.weakness, best.cost, (_gen - best.born),
+            worst.weakness, worst.cost, (_gen - worst.born));
+}
+
+Program
+Population::mutate(const Program &a)
+{
+    Program new_prog = a;
+    do {
+        new_prog.mutate(_rnd);
+    } while(_rnd.get(0,99) < 80);
+    new_prog.reborn(_gen);
+    return new_prog;
+}
+
+void
+Population::init(const Program &program)
+{
+    _programs.clear();
+    _programs.emplace_back(_repo, _params.in_cnt, _params.out_cnt, _params.alt_cnt, _gen);
+    _programs.back().init(program);
+    _programs.back().grow(_rnd, _params.op_cnt);
+    _repo.find_weakness(_rnd, _programs.back());
+}
+
+void
+Population::tick()
+{
+    ++_gen;
+    while (_programs.size() > 1) {
+        _programs.pop_back();
+    }
+    while (_programs.size() < _params.pop_cnt) {
+        _programs.push_back(mutate(_programs[0]));
+        _repo.find_weakness(_rnd, _programs.back());
+    }
+    std::sort(_programs.begin(), _programs.end());
+}
+
+} // namespace vespalib::gp
diff --git a/eval/src/vespa/eval/gp/gp.h b/eval/src/vespa/eval/gp/gp.h
new file mode 100644
index 00000000000..88d187a40b6
--- /dev/null
+++ b/eval/src/vespa/eval/gp/gp.h
@@ -0,0 +1,231 @@
+// Copyright 2017 Yahoo Holdings. Licensed under the terms of the Apache 2.0 license. See LICENSE in the project root.
+
+#pragma once
+
+#include <vespa/vespalib/stllike/string.h>
+#include <random>
+#include <chrono>
+#include <assert.h>
+
+namespace vespalib::gp {
+
+using Value = int; // all input/output/intermediate values have this type (should be template parameter)
+using Weakness = double;
+
+// high level training parameters
+struct Params {
+    size_t in_cnt;  // # function inputs
+    size_t out_cnt; // # function outputs
+    size_t op_cnt;  // # internal operations per individual
+    size_t alt_cnt; // # alternative outputs
+    size_t pop_cnt; // # individuals in population
+    Params(size_t in_cnt_in, size_t out_cnt_in,
+           size_t op_cnt_in, size_t alt_cnt_in,
+           size_t pop_cnt_in)
+        : in_cnt(in_cnt_in), out_cnt(out_cnt_in),
+          op_cnt(op_cnt_in), alt_cnt(alt_cnt_in),
+          pop_cnt(pop_cnt_in) {}
+};
+
+using Input = std::vector<Value>;       // input values
+using Output = std::vector<Value>;      // output values
+using Result = std::vector<Output>;     // alternative output values
+using Feedback = std::vector<Weakness>; // weakness per result alternative
+
+// simple random generator
+struct Random {
+    std::mt19937 gen;
+    Random(int seed) : gen(seed) {}
+    Random() : Random(std::chrono::system_clock::now().time_since_epoch().count()) {}
+    int get(int min, int max) {
+        std::uniform_int_distribution<int> dist(min, max);
+        return dist(gen);
+    }
+    int make_seed() {
+        return get(std::numeric_limits<int>::lowest(),
+                   std::numeric_limits<int>::max());
+    }
+};
+
+// Multiple alternatives for a function taking multiple inputs
+// producing multiple outputs
+struct MultiFunction {
+    virtual size_t num_inputs() const = 0;
+    virtual size_t num_outputs() const = 0;
+    virtual size_t num_alternatives() const = 0;
+    virtual Result execute(const Input &input) const = 0;
+    virtual ~MultiFunction() {}
+};
+
+// simulated individual representing a multifunction
+struct Sim : public MultiFunction {
+    virtual void handle_feedback(Random &rnd, const Feedback &feedback) = 0;
+};
+
+// available operations
+struct OpRepo {
+    using feedback_fun = Feedback (*)(const MultiFunction &multi_fun);
+    using value_op2 = Value (*)(Value lhs, Value rhs);
+    static Value forward_op(Value lhs, Value) { return lhs; }
+    struct Entry {
+        vespalib::string name;
+        value_op2 fun;
+        size_t cost;
+        Entry(const vespalib::string &name_in, value_op2 fun_in, size_t cost_in)
+            : name(name_in), fun(fun_in), cost(cost_in) {}
+    };
+    feedback_fun _find_weakness;
+    std::vector<Entry> _list;
+    OpRepo(feedback_fun find_weakness_in)
+        : _find_weakness(find_weakness_in), _list()
+    {
+        _list.emplace_back("forward", forward_op, 0);
+    }
+    OpRepo &add(const vespalib::string &name, value_op2 fun) {
+        _list.emplace_back(name, fun, 1);
+        return *this;
+    }
+    const vespalib::string &name_of(size_t op) const { return _list[op].name; }
+    size_t cost_of(size_t op) const { return _list[op].cost; }
+    size_t max_op() const { return (_list.size() - 1); }
+    void find_weakness(Random &rnd, Sim &sim) const {
+        sim.handle_feedback(rnd, _find_weakness(sim));
+    }
+    Value perform(size_t op, Value lhs, Value rhs) const {
+        return _list[op].fun(lhs, rhs);
+    }
+};
+
+// specific simulated individual implementation
+struct Program : public Sim {
+    class Ref {
+    private:
+        int idx; // negative: input, zero/positive: operation_result
+        Ref(int idx_in) : idx(idx_in) {}
+    public:
+        bool is_input() const { return (idx < 0); }
+        bool is_operation() const { return (idx >= 0); }
+        size_t in_idx() const { return -(idx + 1); }
+        size_t op_idx() const { return idx; }
+        bool operator<(const Ref &rhs) const { return (idx < rhs.idx); }
+        static Ref in(size_t idx_in) { return Ref(-int(idx_in + 1)); }
+        static Ref op(size_t idx_in) { return Ref(idx_in); }
+        static Ref nop() { return in(0); }
+        static Ref rnd(Random &rnd_in, size_t in_cnt, size_t op_cnt) {
+            return Ref(rnd_in.get(-in_cnt, op_cnt - 1));
+        }
+    };
+    struct Op {
+        size_t code;
+        Ref    lhs;
+        Ref    rhs;
+        Op(size_t code_in, Ref lhs_in, Ref rhs_in)
+            : code(code_in), lhs(lhs_in), rhs(rhs_in) {}
+    };
+    struct Stats {
+        Weakness weakness;
+        size_t   cost;
+        size_t   born;
+        size_t   alt;
+        Stats(size_t gen) : weakness(0.0), cost(0), born(gen), alt(0) {}
+        Stats(Weakness weakness_in, size_t cost_in, size_t born_in, size_t alt_in)
+            : weakness(weakness_in), cost(cost_in), born(born_in), alt(alt_in) {}
+        bool operator<(const Stats &rhs) const {
+            if (weakness != rhs.weakness) {
+                return (weakness < rhs.weakness);
+            }
+            if (cost != rhs.cost) {
+                return (cost < rhs.cost);
+            }
+            return (born > rhs.born); // younger is better
+        }
+    };
+
+    OpRepo           _repo;
+    Stats            _stats;
+    double           _waste;
+    size_t           _in_cnt;
+    size_t           _out_cnt;
+    size_t           _alt_cnt;
+    std::vector<Op>  _program;
+    size_t           _frozen;
+    std::vector<Ref> _bound;
+
+    size_t get_alt_size() const { return (_out_cnt - _bound.size()); }
+    size_t get_alt_offset(size_t alt) const {
+        assert(alt < _alt_cnt);
+        size_t r_offset = (_alt_cnt - alt) * get_alt_size();
+        assert(_program.size() >= r_offset);
+        return (_program.size() - r_offset);
+    }
+
+    void assert_valid(Ref ref, size_t limit) const;
+
+    size_t rnd_op(Random &rnd) { return rnd.get(0, _repo.max_op()); }
+    Ref rnd_ref(Random &rnd, size_t limit) { return Ref::rnd(rnd, _in_cnt, limit); }
+
+    Program(Program &&) = default;
+    Program &operator=(Program &&) = default;
+    Program(const Program &) = default;
+    Program &operator=(const Program &) = delete;
+    ~Program() {}
+
+    Program(const OpRepo &repo, size_t in_cnt, size_t out_cnt, size_t alt_cnt, size_t gen)
+        : _repo(repo), _stats(gen), _waste(0.0),
+          _in_cnt(in_cnt), _out_cnt(out_cnt), _alt_cnt(alt_cnt),
+          _program(), _frozen(0), _bound() {}
+    Ref add_op(size_t code, Ref lhs, Ref rhs);
+    Ref add_forward(Ref ref);
+
+    void init(const Program &src);
+    void grow(Random &rnd, size_t op_cnt);
+    void mutate(Random &rnd, size_t mut_idx);
+    void mutate(Random &rnd);
+    void reborn(size_t gen) { _stats.born = gen; }
+    const Stats &stats() const { return _stats; }
+    bool operator<(const Program &rhs) const {
+        return (stats() < rhs.stats() ||
+                ((_waste < rhs._waste) && !(rhs.stats() < stats())));
+    }
+
+    std::vector<Ref> get_refs(size_t alt) const;
+    std::vector<bool> get_used_ops(size_t alt) const;
+    size_t get_cost(size_t alt) const;
+
+    size_t size_of(Ref ref) const;
+    vespalib::string as_string(Ref ref) const;
+
+    // implementation of the Sim interface
+    size_t num_inputs() const override { return _in_cnt; }
+    size_t num_outputs() const override { return _out_cnt; }
+    size_t num_alternatives() const override { return _alt_cnt; }
+    Result execute(const Input &input) const override;
+    void handle_feedback(Random &rnd, const Feedback &feedback) override;
+};
+
+struct Population
+{
+    Random _rnd;
+    size_t _gen;
+    Params _params;
+    OpRepo _repo;
+    std::vector<Program> _programs;
+
+    void grow(size_t cnt);
+    void print_stats() const;
+
+    Population(const Params &params, const OpRepo &repo, int seed)
+        : _rnd(seed),
+          _gen(0),
+          _params(params),
+          _repo(repo),
+          _programs()
+    {
+        grow(1);
+    }
+    Program mutate(const Program &a);
+    void init(const Program &program);
+    void tick();
+};
+
+} // namespace vespalib::gp