diff --git a/c_src/bayesian.c b/c_src/bayesian.c
new file mode 100644
index 0000000..b3db3f6
--- /dev/null
+++ b/c_src/bayesian.c
@@ -0,0 +1,210 @@
+#ifndef BAYESIAN_H
+#define BAYESIAN_H
+
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+#include "random.h"
+#include "random_forest.h"
+#include "datapoint.h"
+
+#define RANDOM_SAMPLE_PROBABILITY 0.05
+#define NUMBER_OF_TREES 10
+#define EVO_POPULATION 50
+#define MUTATE_PROB 0.1
+
+typedef double (*OptFunc)(const double*);
+
+void get_random(const size_t n, const double* lower, const double* upper, double* x) {
+    for (size_t i=0; i<n; ++i) {
+        x[i] = lower[i] + (upper[i] - lower[i]) * random_f64();
+    }
+}
+
+static inline double pdf(const double x) {
+    return exp(-x*x/2.0) / sqrt(2.0 * M_PI);
+}
+
+static inline double cdf(const double x) {
+    return 0.5 * (1.0 + erf(x / M_SQRT2));
+}
+
+static inline double ei(const double y, const double sdev, const double fmin) {
+    const double delta = fmin - y;
+    const double z = delta / sdev;
+    return delta + sdev * pdf(z) + delta * cdf(z);
+}
+
+void optimize_acquisition_function(RandomForest* forest, const double* lower, const double* upper,
+                     double* out, const size_t dims, const double fmin, DataPoint** points) {
+
+    // This will be our workhorse x
+    DataPoint* p = points[EVO_POPULATION];
+
+    double best = -INFINITY;
+    
+    // Initialize the workspace with random points.
+    for (size_t i=0; i<EVO_POPULATION; ++i) {
+        get_random(dims, lower, upper, points[i]->x);
+        double variance;
+        double y = randomforest_full_predict(forest, points[i]->x, &variance);
+        points[i]->y = ei(y, sqrt(variance), fmin);
+
+        if (points[i]->y > best) {
+            best = points[i]->y;
+            for (size_t j=0; j<dims; ++j) { out[j] = points[i]->x[j]; }
+        }
+    }
+
+    size_t oldest = 0;
+    for (size_t i=0; i<10000; ++i) {
+
+        DataPoint* p1 = points[0];
+        DataPoint* p2 = points[1];
+        if (p2->y > p1->y) {
+            DataPoint* temp = p1;
+            p1 = p2;
+            p2 = temp;
+        }
+
+        for (size_t j=2; j<EVO_POPULATION; ++j) {
+            if (points[j]->y > p1->y) {
+                p2 = p1;
+                p1 = points[j];
+            } else if (points[j]->y > p2->y) {
+                p2 = points[j];
+            }
+        }
+
+        const double sum = p1->y + p2->y;
+        const double p1_prob = p1->y / sum;
+
+        for (size_t j=0; j<dims; ++j) {
+            p->x[j] = (random_f64() < p1_prob) ? p1->x[j] : p2->x[j];
+        }
+
+        for (size_t j=0; j<dims; ++j) {
+            if (random_f64() < MUTATE_PROB) {
+                p->x[j] = lower[j] + (upper[j] - lower[j]) * random_f64();
+            }
+        }
+
+        double variance;
+        double y = randomforest_full_predict(forest, p->x, &variance);
+        p->y = ei(y, sqrt(variance), fmin);
+
+        if (p->y > best) {
+            best = p->y;
+            for (size_t j=0; j<dims; ++j) { out[j] = p->x[j]; }
+        }
+
+        DataPoint* temp = points[oldest];
+        points[oldest] = p;
+        p = temp;
+
+        oldest = (oldest + 1) % EVO_POPULATION;
+
+    }
+
+    // Make sure to update the backing pointer array, since p may have been swapped
+    points[EVO_POPULATION] = p;
+}
+
+void shuffle_points(DataPoint** points, const size_t n) {
+    for (size_t i=0; i<n-1; ++i) {
+        const size_t j = i + random_u64() % (n-i);
+        DataPoint* temp = points[i];
+        points[i] = points[j];
+        points[j] = temp;
+    }
+}
+
+double bayesian_optimization(OptFunc f, const double* lower, const double* upper, double* x, const size_t dims,
+                             const size_t doe, const size_t iter) {
+
+    const size_t dp_size = sizeof(DataPoint) + dims * sizeof(double);
+
+    const size_t total_iter = doe + iter;
+
+    // Alloc backing data for all points (1 per iteration)
+    DataPoint* data = calloc(total_iter, dp_size);
+    
+    // Alloc pointers that can point into backing data points
+    DataPoint** points = calloc(total_iter, sizeof(DataPoint*));
+
+    // Assign each pointer to an allocated data point
+    points[0] = data;
+    for (size_t i=1; i<total_iter; ++i) { points[i] = (void*)points[i-1] + dp_size; }
+
+    // Alloc work space for local search algorithm
+    DataPoint* workspace = calloc(EVO_POPULATION + 1, dp_size);
+
+    // Pointers into workspace
+    DataPoint** workspace_points = calloc(EVO_POPULATION + 1, sizeof(DataPoint*));
+
+    workspace_points[0] = workspace;
+    for (size_t i=1; i<EVO_POPULATION+1; ++i) {
+        workspace_points[i] = (void*)workspace_points[i-1] + dp_size;
+    }
+
+    // This will correspond to the point in the parameter x
+    double best_val = INFINITY;
+
+    // Initialize DoE phase as a latin hypercube
+    for (size_t i=0; i<dims; ++i) {
+        for (size_t j=0; j<doe; ++j) {
+            points[j]->x[i] = lower[i] + (upper[i] - lower[i]) * (j + random_f64()) / doe;
+        }
+        shuffle_points(points, doe);
+    }
+
+    // Evaluate DoE points
+    for (size_t i=0; i<doe; ++i) {
+        const double y = f(points[i]->x);
+        points[i]->y = y;
+
+        if (y < best_val) {
+            best_val = y;
+            for (size_t j=0; j<dims; ++j) { x[j] = points[i]->x[j]; }
+        }
+    }
+
+    double mean_iter = 0.0;
+
+    for (size_t i=doe; i<total_iter; ++i) {
+
+        if (random_f64() < RANDOM_SAMPLE_PROBABILITY) {
+
+            get_random(dims, lower, upper, points[i]->x);
+
+        } else {
+            // Else get a guess from random forest model
+
+            // Fit forest
+            RandomForest forest;
+            randomforest_fit(&forest, NUMBER_OF_TREES, points, dims, i);
+
+            optimize_acquisition_function(&forest, lower, upper, points[i]->x, dims, best_val, workspace_points);
+
+            randomforest_free(&forest);
+        }
+
+        const double y = f(points[i]->x);
+        points[i]->y = y;
+
+        if (y < best_val) {
+            best_val = y;
+            for (size_t j=0; j<dims; ++j) { x[j] = points[i]->x[j]; }
+        }
+    }
+
+    free(workspace_points);
+    free(workspace);
+    free(points);
+    free(data);
+
+    return best_val;
+}
+
+#endif
diff --git a/c_src/datapoint.h b/c_src/datapoint.h
new file mode 100644
index 0000000..75d36fb
--- /dev/null
+++ b/c_src/datapoint.h
@@ -0,0 +1,9 @@
+#ifndef DATASET_H
+#define DATASET_H
+
+typedef struct DataPoint {
+    double y;
+    double x[];
+} DataPoint;
+
+#endif
diff --git a/c_src/random.h b/c_src/random.h
new file mode 100644
index 0000000..264489d
--- /dev/null
+++ b/c_src/random.h
@@ -0,0 +1,121 @@
+#ifndef RANDOM_H
+#define RANDOM_H
+
+#include <stdint.h>
+#include <stdbool.h>
+#include <math.h>
+
+/***
+ * This file contains an implementation of the SplitMix generator,
+ * and the xoshiro256** generator.
+ * They are adapted from:
+ * https://prng.di.unimi.it/splitmix64.c and
+ * https://prng.di.unimi.it/xoshiro256starstar.c
+ * respectively.
+ */
+
+/**
+ * Returns a double-precision floating-point number in the range [0, 1),
+ * using bit manipulation on the given 64-bit integer.
+ */
+static inline double u64_to_f64(const uint64_t x) {
+    union {
+        uint64_t u64;
+        double f64;
+    } u = { .u64 = 0x3ff0000000000000ull | (x >> 12) };
+    return 2.0 - u.f64;
+}
+
+
+typedef uint64_t Splitmix;
+
+static inline void splitmix_init(Splitmix* state, const uint64_t seed) {
+    *state = seed;
+}
+
+uint64_t splitmix_u64(Splitmix* state) {
+    uint64_t z = (*state += 0x9e3779b97f4a7c15);
+    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+    return z ^ (z >> 31);
+}
+
+static inline double splitmix_f64(Splitmix* state) {
+    return u64_to_f64(splitmix_u64(state));
+}
+
+typedef struct Xoshiro256 {
+    uint64_t s[4];
+} Xoshiro256;
+
+void xoshiro256_init(Xoshiro256* state, const uint64_t seed) {
+    Splitmix sm;
+    splitmix_init(&sm, seed);
+
+    state->s[0] = splitmix_u64(&sm);
+    state->s[1] = splitmix_u64(&sm);
+    state->s[2] = splitmix_u64(&sm);
+    state->s[3] = splitmix_u64(&sm);
+}
+
+static inline uint64_t rotl64(const uint64_t x, int k) {
+    return (x << k) | (x >> (64 - k));
+}
+
+uint64_t xoshiro256_u64(Xoshiro256* state) {
+    const uint64_t result = rotl64(state->s[1] * 5, 7) * 9;
+
+    const uint64_t t = state->s[1] << 17;
+
+    state->s[2] ^= state->s[0];
+    state->s[3] ^= state->s[1];
+    state->s[1] ^= state->s[2];
+    state->s[0] ^= state->s[3];
+
+    state->s[2] ^= t;
+
+    state->s[3] = rotl64(state->s[3], 45);
+
+    return result;
+}
+
+static inline double xoshiro256_f64(Xoshiro256* state) {
+    return u64_to_f64(xoshiro256_u64(state));
+}
+
+Xoshiro256 global_xoshiro256_state;
+
+extern inline void random_init(const uint64_t seed) {
+    xoshiro256_init(&global_xoshiro256_state, seed);
+}
+
+static inline uint64_t random_u64(void) {
+    return xoshiro256_u64(&global_xoshiro256_state);
+}
+
+static inline double random_f64(void) {
+    return xoshiro256_f64(&global_xoshiro256_state);
+}
+
+double random_normal(void) {
+    double v, u, q;
+    do {
+        u = random_f64();
+        v = 1.7156*(random_f64() - 0.5);
+        const double x = u - 0.449871;
+        const double y = fabs(v) + 0.386595;
+        q = x*x + y*(0.19600*y-0.25472*x);
+    } while (q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
+    return v / u;
+}
+
+double random_cauchy(void) {
+    double v1, v2;
+    do {
+        v1 = 2.0 * random_f64() - 1.0;
+        v2 = random_f64();
+    } while (v1*v1+v2*v2 >= 1.0 || v2 < 0x1p-63);
+    return v1 / v2;
+}
+
+#endif
diff --git a/c_src/random_forest.h b/c_src/random_forest.h
new file mode 100644
index 0000000..7c59e5c
--- /dev/null
+++ b/c_src/random_forest.h
@@ -0,0 +1,348 @@
+#ifndef RANDOM_FOREST_H
+#define RANDOM_FOREST_H
+
+#include <stdbool.h>
+#include <math.h>
+#include <stdlib.h>
+
+#include "random.h"
+#include "datapoint.h"
+
+#define REGTREE_MIN_SPLIT 5
+
+typedef enum RegressionTreeNodeType {
+    REGTREENODE_TYPE_UNINIT = 0,
+    REGTREENODE_TYPE_SPLIT,
+    REGTREENODE_TYPE_LEAF
+} RegressionTreeNodeType;
+
+typedef struct RegressionTreeNode {
+    union {
+        // A split node is constructed from a collection of {x, y} pairs.
+        // Each x is a vector of attributes.
+        // A split node splits the input based on a split index and split value,
+        // such that:
+        //     if (x[split index] <= split value) goto left subtree
+        //     else goto right subtree
+        struct {
+            size_t split_index;
+            double split_value;
+            struct RegressionTreeNode* left;
+            struct RegressionTreeNode* right;
+        } split;
+
+        // A leaf node is constructed from a collection of {x, y} pairs.
+        // All information about x is encoded in the split-nodes that lead to a leaf,
+        // so here we're only concerned about the y-information.
+        struct {
+            double y; ///< Sum of all y
+            double y2; ///< Sum of all y^2
+            double avg; ///< y/n
+            size_t n; ///< Number of {x, y} pairs in this leaf
+        } leaf;
+    };
+    RegressionTreeNodeType type;
+} RegressionTreeNode;
+
+/**
+ * Calculates and returns the variance in y for a collection of points
+ */
+double variance(DataPoint** points, const size_t n) {
+    double y = 0.0;
+    double y2 = 0.0;
+    for (size_t i=0; i<n; ++i) {
+        y += points[i]->y;
+        y2 += points[i]->y * points[i]->y;
+    }
+    return y2 - y * y / n;
+}
+
+/**
+ * Finds the minimum and maximun y-vale among points, and returns them in min and max respectively.
+ */
+void minmax_y(DataPoint** points, const size_t n, double* min, double* max) {
+    double lo = points[0]->y;
+    double hi = points[0]->y;
+    for (size_t i=1; i<n; ++i) {
+        if (points[i]->y < lo) { lo = points[i]->y; }
+        if (points[i]->y > hi) { hi = points[i]->y; }
+    }
+    *min = lo;
+    *max = hi;
+}
+
+/**
+ * Finds the minimum and maximum value in x[index] among points, and returns them in min/max.
+ */
+void minmax_x(DataPoint** points, const size_t n, const size_t index, double* min, double* max) {
+    double lo = points[0]->x[index];
+    double hi = points[0]->x[index];
+    for (size_t i=1; i<n; ++i) {
+        if (points[i]->x[index] < lo) { lo = points[i]->x[index]; }
+        if (points[i]->x[index] > hi) { hi = points[i]->x[index]; }
+    }
+    *min = lo;
+    *max = hi;
+}
+
+double get_split_variance(DataPoint** points, const size_t n,
+                          const size_t split_index, const double split_value) {
+    double y_left = 0.0;
+    double y2_left = 0.0;
+    size_t n_left = 0;
+
+    double y_right = 0.0;
+    double y2_right = 0.0;
+    size_t n_right = 0;
+
+    for (size_t i=0; i<n; ++i) {
+        if (points[i]->x[split_index] <= split_value) {
+            y_left += points[i]->y;
+            y2_left += points[i]->y * points[i]->y;
+            n_left += 1;
+        } else {
+            y_right += points[i]->y;
+            y2_right += points[i]->y * points[i]->y;
+            n_right += 1;
+        }
+    }
+
+    const double var_left = y2_left - y_left * y_left / n_left;
+    const double var_right = y2_right - y_right * y_right / n_right;
+
+    return var_left + var_right;
+}
+
+bool split(DataPoint** points, const size_t dims, const size_t n, size_t* split_index, double* split_value) {
+
+    //fprintf(stderr, "Trying to split %zu points\n", n);
+
+    if (n < REGTREE_MIN_SPLIT) {
+        //fprintf(stderr, "Too few\n");
+        return false;
+    }
+
+    // If the range of values is small,
+    // there is probably no point in splitting,
+    // so we early-out here.
+    double min_y, max_y;
+    minmax_y(points, n, &min_y, &max_y);
+    if (fabs(min_y - max_y) < 0x1p-63) {
+        //fprintf(stderr, "Small range: [%f, %f]\n", min_y, max_y);
+        return false;
+    }
+
+    const double var_before = variance(points, n);
+    //fprintf(stderr, "Variance before: %f\n", var_before);
+
+    double best_split_var = INFINITY;
+    double best_split_val = 0.0;
+    size_t best_split_index = 0;
+    for (size_t i=0; i<dims; ++i) {
+
+        // Get minimum and maximum x
+        double min_x, max_x;
+        minmax_x(points, n, i, &min_x, &max_x);
+
+        // Compute a random splitting value, and get variance for that split
+        const double split_val = min_x + (max_x - min_x) * random_f64();
+        const double split_var = get_split_variance(points, n, i, split_val);
+
+        if (split_var < best_split_var) {
+            best_split_var = split_var;
+            best_split_val = split_val;
+            best_split_index = i;
+        }
+    }
+    //fprintf(stderr, "Found best split in dimension %zu: val = %f, var = %f\n",
+            //best_split_index, best_split_val, best_split_var);
+
+    // No significant reduction in variance
+    if (var_before - best_split_var < 0x1p-63) {
+        //fprintf(stderr, "No significant reduction: before = %f, after = %f\n", var_before, best_split_var);
+        return false;
+    }
+
+    *split_index = best_split_index;
+    *split_value = best_split_val;
+    return true;
+}
+
+void make_leaf(RegressionTreeNode* node, DataPoint** points, const size_t n) {
+    double y = 0.0;
+    double y2 = 0.0;
+    for (size_t i=0; i<n; ++i) {
+        y += points[i]->y;
+        y2 += points[i]->y * points[i]->y;
+    }
+
+    node->type = REGTREENODE_TYPE_LEAF;
+    node->leaf.y = y;
+    node->leaf.y2 = y2;
+    node->leaf.avg = y / n;
+    node->leaf.n = n;
+}
+
+void regressiontreenode_fit(RegressionTreeNode* node,
+                            DataPoint** points, const size_t dims, const size_t n) {
+
+    //fprintf(stderr, "Fitting %zu points\n", n);
+    size_t split_index;
+    double split_value;
+
+    if (!split(points, dims, n, &split_index, &split_value)) {
+        // If we failed to find a working split, make this node a leaf and return
+        //fprintf(stderr, "No split\n");
+        make_leaf(node, points, n);
+        return;
+    }
+
+    // Here we begin to partition points into
+    // left: x[split_index] <= split_value
+    // right: x[split_index > split_value
+
+    // Find the first 'right' point
+    size_t hi = 0;
+    while (hi < n && points[hi]->x[split_index] <= split_value) { hi += 1; }
+
+    // If no 'right' point was found, then all are left,
+    // and there is no split,
+    // so we make a leaf...
+    if (hi >= n) {
+        //fprintf(stderr, "Bad split\n");
+        make_leaf(node, points, n);
+        return;
+    }
+
+    // Perform actual partitioning
+    for (size_t i=hi+1; i<n; ++i) {
+        if (points[i]->x[split_index] <= split_value) {
+            DataPoint* temp = points[i];
+            points[i] = points[hi];
+            points[hi] = temp;
+            hi += 1;
+        }
+    }
+
+    // Partitioning done
+
+    const size_t n_left = hi;
+    const size_t n_right = n - hi;
+
+    //fprintf(stderr, "Left x[%zu] <= %f: %zu\n", split_index, split_value, n_left);
+    //fprintf(stderr, "Right x[%zu] > %f: %zu\n", split_index, split_value, n_right);
+
+    // Fit left subtree
+    //fprintf(stderr, "Fitting left\n");
+    RegressionTreeNode* left = calloc(1, sizeof(RegressionTreeNode));
+    regressiontreenode_fit(left, points, dims, n_left);
+
+    // Fit right subtree
+    //fprintf(stderr, "Fitting right\n");
+    RegressionTreeNode* right = calloc(1, sizeof(RegressionTreeNode));
+    regressiontreenode_fit(right, points+hi, dims, n_right);
+
+    // Construct this node
+    node->type = REGTREENODE_TYPE_SPLIT;
+    node->split.split_index = split_index;
+    node->split.split_value = split_value;
+    node->split.left = left;
+    node->split.right = right;
+}
+
+void regressiontreenode_free(RegressionTreeNode* node) {
+    if (node->type == REGTREENODE_TYPE_SPLIT) {
+        regressiontreenode_free(node->split.left);
+        free(node->split.left);
+        regressiontreenode_free(node->split.right);
+        free(node->split.right);
+    }
+}
+
+typedef struct RegressionTree {
+    RegressionTreeNode* root;
+} RegressionTree;
+
+static inline void regressiontree_fit(RegressionTree* tree,
+                                      DataPoint** points, const size_t dims, const size_t n) {
+    tree->root = calloc(1, sizeof(RegressionTreeNode));
+    regressiontreenode_fit(tree->root, points, dims, n);
+}
+
+RegressionTreeNode* regressiontree_propagate(RegressionTree* tree, const double* x) {
+    RegressionTreeNode* n = tree->root;
+    while (n->type != REGTREENODE_TYPE_LEAF) {
+        if (x[n->split.split_index] <= n->split.split_value) {
+            n = n->split.left;
+        } else {
+            n = n->split.right;
+        }
+    }
+    return n;
+}
+
+static inline double regressiontree_predict(RegressionTree* tree, const double* x) {
+    RegressionTreeNode* n = regressiontree_propagate(tree, x);
+    return n->leaf.avg;
+}
+
+double regressiontree_fulL_predict(RegressionTree* tree, const double* x, double* variance) {
+    RegressionTreeNode* n = regressiontree_propagate(tree, x);
+    const double y = n->leaf.y / n->leaf.n;
+    const double y2 = n->leaf.y2 / n->leaf.n;
+    *variance = y2 - y * y;
+    return y;
+}
+
+void regressiontree_free(RegressionTree* tree) {
+    regressiontreenode_free(tree->root);
+    free(tree->root);
+}
+
+typedef struct RandomForest {
+    RegressionTree* trees;
+    size_t n_trees;
+} RandomForest;
+
+void randomforest_fit(RandomForest* forest, const size_t n_trees,
+                      DataPoint** points, const size_t dims, const size_t n) {
+    RegressionTree* trees = calloc(n_trees, sizeof(RegressionTree));
+    for (size_t i=0; i<n_trees; ++i) {
+        regressiontree_fit(&trees[i], points, dims, n);
+    }
+    forest->trees = trees;
+    forest->n_trees = n_trees;
+}
+
+double randomforest_predict(RandomForest* forest, const double* x) {
+    double tot = 0.0;
+    for (size_t i=0; i<forest->n_trees; ++i) {
+        tot += regressiontree_predict(&forest->trees[i], x);
+    }
+    return tot / forest->n_trees;
+}
+
+double randomforest_full_predict(RandomForest* forest, const double* x, double* variance) {
+    double y = 0.0;
+    double y2 = 0.0;
+    size_t n = 0;
+    for (size_t i=0; i<forest->n_trees; ++i) {
+        RegressionTreeNode* node = regressiontree_propagate(&forest->trees[i], x);
+        y += node->leaf.y;
+        y2 += node->leaf.y2;
+        n += node->leaf.n;
+    }
+    y /= n;
+    y2 /= n;
+    *variance = y2 - y * y;
+    return y;
+}
+
+void randomforest_free(RandomForest* forest) {
+    for (size_t i=0; i<forest->n_trees; ++i) {
+        regressiontree_free(&forest->trees[i]);
+    }
+    free(forest->trees);
+}
+
+#endif
diff --git a/hypermapper/bayesian_c.py b/hypermapper/bayesian_c.py
new file mode 100644
index 0000000..3c1ff95
--- /dev/null
+++ b/hypermapper/bayesian_c.py
@@ -0,0 +1,49 @@
+from ctypes import *
+import time
+import site
+import os
+
+site_packages = site.getusersitepackages()
+libfile = site_packages + '/' + \
+          [x for x in os.listdir(site_packages) if x.find('hypermapper_c_bayesian.') >= 0][0]
+
+bayesian_c = cdll.LoadLibrary(libfile)
+
+bayesian_c.random_init.restype = None
+bayesian_c.random_init(int(time.time()))
+
+bayesian_c.bayesian_optimization.restype = c_double
+
+OptFuncType = CFUNCTYPE(c_double, POINTER(c_double))
+
+class BayesianCWrapper:
+
+    def __init__(self, blackbox, config):
+        self.f = blackbox
+        self.n_doe = config['design_of_experiment']['number_of_samples']
+        self.iter = config['optimization_iterations']
+        self.y = config['optimization_objectives'][0]
+        self.vars = []
+        self.lower = []
+        self.upper = []
+        for name, prop in config['input_parameters'].items():
+            self.vars.append(name)
+            bounds = prop['values']
+            self.lower.append(bounds[0])
+            self.upper.append(bounds[1])
+        self.dims = len(self.vars)
+
+    def f_wrapper(self, x):
+        x_dict = {self.vars[i]: x[i] for i in range(self.dims)}
+        return self.f(x_dict)
+
+    def run(self):
+        x = (self.dims * c_double)(0.0)
+        y = bayesian_c.bayesian_optimization(OptFuncType(self.f_wrapper),
+                                           (self.dims * c_double)(*self.lower),
+                                           (self.dims * c_double)(*self.upper),
+                                           x,
+                                           self.dims, self.n_doe, self.iter)
+        out = {self.vars[i]: [x[i]] for i in range(self.dims)}
+        out[self.y] = [y]
+        return out
diff --git a/hypermapper/bo.py b/hypermapper/bo.py
index 06b1339..f1a9469 100644
--- a/hypermapper/bo.py
+++ b/hypermapper/bo.py
@@ -26,6 +26,7 @@
         create_output_data_file,
         write_data_array,
     )
+    from hypermapper import bayesian_c
 
 except ImportError:
     if os.getenv("HYPERMAPPER_HOME"):  # noqa
@@ -74,6 +75,7 @@
         create_output_data_file,
         write_data_array,
     )
+    from hypermapper import bayesian_c
 
 
 def main(config, black_box_function=None, profiling=None):
@@ -130,6 +132,10 @@ def main(config, black_box_function=None, profiling=None):
     normalize_objectives = False
     debug = False
 
+    if config["use_c_backend"]:
+        c_wrapper = bayesian_c.BayesianCWrapper(black_box_function, config)
+        return c_wrapper.run()
+
     if "feasible_output" in config:
         feasible_output = config["feasible_output"]
         feasible_output_name = feasible_output["name"]
diff --git a/hypermapper/optimizer.py b/hypermapper/optimizer.py
index 05781a0..de141d5 100644
--- a/hypermapper/optimizer.py
+++ b/hypermapper/optimizer.py
@@ -120,6 +120,7 @@ def optimize(parameters_file, black_box_function=None):
         (optimization_method == "random_scalarizations")
         or (optimization_method == "bayesian_optimization")
         or (optimization_method == "prior_guided_optimization")
+        or (optimization_method == "bayesian_optimization_c_backend")
     ):
         data_array = bo.main(
             config, black_box_function=black_box_function, profiling=profiling
diff --git a/hypermapper/schema.json b/hypermapper/schema.json
index b99f5de..05e443e 100644
--- a/hypermapper/schema.json
+++ b/hypermapper/schema.json
@@ -240,6 +240,11 @@
         "enum": ["bayesian_optimization", "random_scalarizations", "local_search", "prior_guided_optimization", "evolutionary_optimization"],
         "default": "bayesian_optimization"
       },
+      "use_c_backend":{
+          "type": "boolean",
+          "description": "Use C backend. The C backend is optimized for speed, but more restrictive in terms of functionality. Currently only available for bayesian_optimization",
+          "default": false
+      },
       "local_search_starting_points":{
         "type": "integer",
         "description": "number of starting points for the multi-start local search used to optimize the acquisition functions.",
diff --git a/setup.py b/setup.py
index 91e6be9..51146b9 100644
--- a/setup.py
+++ b/setup.py
@@ -1,4 +1,4 @@
-from setuptools import setup
+from setuptools import setup, Extension
 
 with open("README.md", "r") as fh:
     long_description = fh.read()
@@ -45,4 +45,5 @@
             "hm-plot-optimization-results=hypermapper._cli:_plot_optimization_results_cli",
         ],
     },
+    ext_modules=[Extension('hypermapper_c_bayesian', ['c_src/bayesian.c'])],
 )