Repository: emlearn/emlearn (language: python)

Issue:
Added support for Regression Models
I added support for regression models for a university project I worked for.
 
I made a copy of the file trees.py called regressor_trees.py that has been modified to support these models.


Discussion / hints:
Hi @SaverioMonaco - thank you for this! Regression has been requested by several others in #31 and #41 .

To integrated it into the main version we would need to avoid copying all of the trees.py module. Instead the existing module should be modified to support both regression and classification, with minimal divergence in codepaths between the two cases.
We also need to have a few tests, to ensure that the implementation is correct. See tests/test_trees.py for examples.

Would you be interested in making these changes, so we could integrate this?


I would like to add the support for regression properly, however it might take some time due to my exams.
Thats great! Best of luck on your exams :)
Hi @SaverioMonaco - are you still interested in working on this?

Fix patch:
diff --git a/bindings/eml_trees.cpp b/bindings/eml_trees.cpp
index d1fd614..5ad4ba2 100644
--- a/bindings/eml_trees.cpp
+++ b/bindings/eml_trees.cpp
@@ -69,6 +69,32 @@ class TreesClassifier {
         return classes;
     }
 
+    py::array_t<int32_t>
+    regress(py::array_t<float, py::array::c_style | py::array::forcecast> in) {
+        if (in.ndim() != 2) {
+            throw std::runtime_error("predict input must have dimensions 2");
+        }
+
+        const int64_t n_samples = in.shape()[0];
+        const int32_t n_features = in.shape()[1];
+
+        auto outputs = py::array_t<float>(n_samples);
+        //auto s = in.unchecked();
+        auto r = outputs.mutable_unchecked<1>();
+        float out[1];
+        for (int i=0; i<n_samples; i++) {
+            const float *v = in.data(i);
+            const EmlError err = eml_trees_regress(&forest, v, n_features, out, 1);
+            if (err != EmlOk) {
+                const std::string msg = eml_error_str(err);
+                throw std::runtime_error(msg);
+            }
+            r(i) = out[0];
+        }
+
+        return outputs;
+    }
+
 };
 
 PYBIND11_MODULE(eml_trees, m) {
@@ -76,6 +102,7 @@ PYBIND11_MODULE(eml_trees, m) {
 
     py::class_<TreesClassifier>(m, "Classifier")
         .def(py::init<std::vector<float>, std::vector<int32_t>>())
+        .def("regress", &TreesClassifier::regress)
         .def("predict", &TreesClassifier::predict);
 }
 
diff --git a/emlearn/bayes.py b/emlearn/bayes.py
index 685750b..4761cd8 100644
--- a/emlearn/bayes.py
+++ b/emlearn/bayes.py
@@ -78,10 +78,19 @@ def __init__(self, estimator, method):
         # FIXME: use var,mean numpy arrays directly
         n_classes, n_features = estimator.theta_.shape
         model = numpy.ndarray(shape=(n_classes, n_features, 3), dtype=float)
+
+        if hasattr(estimator, 'var_'):
+            # 1.0 and later
+            variance = estimator.var_
+        elif hasattr(estimator, 'sigma_'):
+            # pre 1.0
+            variance = estimator.sigma_
+
+
         for class_n in range(0, n_classes):
             for feature_n in range(0, n_features):
                 mean = estimator.theta_[class_n,feature_n]
-                std = numpy.sqrt(estimator.sigma_[class_n,feature_n])
+                std = numpy.sqrt(variance[class_n,feature_n])
                 std_log2 = numpy.log2(std)
                 model[class_n,feature_n] = (mean, std, std_log2)
         self.model = model
diff --git a/emlearn/common.py b/emlearn/common.py
index d5c6288..07bee79 100644
--- a/emlearn/common.py
+++ b/emlearn/common.py
@@ -60,15 +60,24 @@ def build_classifier(cmodel, name, temp_dir, include_dir, func=None, test_functi
 
     return bin_path
 
-def run_classifier(bin_path, data, out_dtype='int'):
+def run_classifier(bin_path, data, out_dtype='int', float_precision=8):
+
+    # Serialize input data as CSV
+    def serialize_value(v):
+        return '{:.{prec}f}'.format(v, prec=float_precision)
+
     lines = []
     for row in data:
-        lines.append(",".join(str(v) for v in row))
+        lines.append(",".join(serialize_value(v) for v in row))
     stdin = '\n'.join(lines)
 
+    assert len(lines) == len(data), (len(lines), data.shape)
+
+    # Run as subprocess
     args = [ bin_path ]
     out = subprocess.check_output(args, input=stdin, encoding='utf8', universal_newlines=True)
 
+    # Parse output
     outputs = []
     for line in out.split('\n'):
         if line:
@@ -81,7 +90,7 @@ def run_classifier(bin_path, data, out_dtype='int'):
                 out_ = out_dtype(out_)
             outputs.append(out_)
 
-    assert len(outputs) == len(data)
+    assert len(outputs) == len(data), (len(outputs), len(data), out)
 
     return outputs
 
@@ -96,6 +105,10 @@ def __init__(self, cmodel, name, call=None, include_dir=None, temp_dir='tmp', te
     def predict(self, X):
         return run_classifier(self.bin_path, X, out_dtype=self._out_dtype)
 
+    def regress(self, X):
+        return self.predict(X)
+
+
 
 def compile_executable(code_file,
                     out_dir,
diff --git a/emlearn/convert.py b/emlearn/convert.py
index 1fcd3c9..38b40e2 100644
--- a/emlearn/convert.py
+++ b/emlearn/convert.py
@@ -12,7 +12,7 @@ def convert(estimator, kind=None, method='pymodule', dtype='float'):
         kind = type(estimator).__name__
 
     # Use name instead of instance to avoid hard dependency on the libraries
-    if kind in ['RandomForestClassifier', 'ExtraTreesClassifier', 'DecisionTreeClassifier']:
+    if kind in set(trees.SUPPORTED_ESTIMATORS):
         return trees.Wrapper(estimator, method, dtype=dtype)
     elif kind in ['EllipticEnvelope']:
         return distance.Wrapper(estimator, method, dtype=dtype)
diff --git a/emlearn/eml_trees.h b/emlearn/eml_trees.h
index b8719f6..c53c017 100644
--- a/emlearn/eml_trees.h
+++ b/emlearn/eml_trees.h
@@ -3,6 +3,8 @@
 #define EMTREES_H
 
 #include <stdint.h>
+#include <math.h>
+#include <eml_common.h>
 
 typedef struct _EmlTreesNode {
     int8_t feature;
@@ -90,4 +92,40 @@ eml_trees_predict(const EmlTrees *forest, const float *features, int8_t features
     return most_voted_class;
 }
 
+
+EmlError
+eml_trees_regress(const EmlTrees *forest,
+        const float *features, int8_t features_length,
+        float *out, int8_t out_length)
+{
+    if (out_length < 1) {
+        return EmlSizeMismatch;
+    }
+
+    float sum = 0; 
+
+    for (int32_t i=0; i<forest->n_trees; i++) {
+        const float val = eml_trees_predict_tree(forest, forest->tree_roots[i], features, features_length);
+        sum += val;
+    }
+
+    out[0] = sum / forest->n_trees;
+
+    return EmlOk;
+}
+
+float
+eml_trees_regress1(const EmlTrees *forest,
+        const float *features, int8_t features_length)
+{
+    float out[1];
+    EmlError err = eml_trees_regress(forest,
+        features, features_length,
+        out, 1);
+    if (err != EmlOk) {    
+        return nanf("");
+    }
+    return out[0];
+}
+
 #endif // EMTREES_H
diff --git a/emlearn/trees.py b/emlearn/trees.py
index 14dbadf..07b62cc 100644
--- a/emlearn/trees.py
+++ b/emlearn/trees.py
@@ -7,12 +7,20 @@
 
 from . import common
 
+SUPPORTED_ESTIMATORS=[
+    'RandomForestClassifier',
+    'ExtraTreesClassifier',
+    'DecisionTreeClassifier',
+    'RandomForestRegressor',
+    'ExtraTreesRegressor',
+    'DecisionTreeRegressor',
+]
 
 # Tree representation as 2d array
 # feature, value, left_child, right_child
 # Leaf node: -1, class, -1, -1
 # 	fields = { 'feature': 0, 'value': 1, 'left': 2, 'right': 3 }
-def flatten_tree(tree):
+def flatten_tree(tree, leaf='argmax'):
     flat = []
 
     assert tree.node_count == len(tree.value)
@@ -20,8 +28,12 @@ def flatten_tree(tree):
 
     for left, right, feature, th, value in zip(tree.children_left, tree.children_right, tree.feature, tree.threshold, tree.value):
         if left == -1 and right == -1:
-            cls = numpy.argmax(value[0])
-            n = [ -1, cls, -1, -1 ] # leaf
+            if leaf == 'argmax':
+                val = numpy.argmax(value[0])
+            elif leaf == 'value':
+                val = value[0][0]
+
+            n = [ -1, val, -1, -1 ] # leaf
         else:
             n = [ feature, th, left, right ]
 
@@ -30,13 +42,13 @@ def flatten_tree(tree):
     return flat
 
 
-def flatten_forest(trees):
+def flatten_forest(trees, leaf='argmax'):
     tree_roots = []
     tree_offset = 0
     forest_nodes = []
 
     for tree in trees: 
-        flat = flatten_tree(tree)
+        flat = flatten_tree(tree, leaf=leaf)
 
         # Offset the nodes in tree, so they can be stored in one array 
         root = 0 + tree_offset
@@ -198,7 +210,7 @@ def node(n):
 
     return out
 
-def generate_c_inlined(forest, name, dtype='float'):
+def generate_c_inlined(forest, name, dtype='float', classifier=True):
     nodes, roots = forest
 
     def is_leaf(n):
@@ -207,9 +219,13 @@ def class_value(n):
       assert is_leaf(n)
       return n[1]
 
-    class_values = set(map(class_value, filter(is_leaf, nodes)))
-    assert min(class_values) == 0
-    n_classes = max(class_values)+1
+    if classifier:
+        class_values = set(map(class_value, filter(is_leaf, nodes)))
+        assert min(class_values) == 0
+        n_classes = max(class_values)+1
+    else:
+        n_classes = numpy.shape(forest[1])[0]
+
     tree_names = [ name + '_tree_{}'.format(i) for i,_ in enumerate(roots) ]
 
     ctype = dtype
@@ -236,7 +252,7 @@ def c_node(nid, depth):
             return c_leaf(n, depth+1)
         return c_internal(n, depth+1)
 
-    def tree_func(name, root):
+    def tree_func(name, root, return_type='int32_t'):
         return """static inline int32_t {function_name}(const {ctype} *features, int32_t features_length) {{
         {code}
         }}
@@ -244,14 +260,31 @@ def tree_func(name, root):
             'function_name': name,
             'code': c_node(root, 0),
             'ctype': ctype,
+            'return_type': return_type 
         })
 
-    def tree_vote(name):
+    def tree_vote_classifier(name):
         return '_class = {}(features, features_length); votes[_class] += 1;'.format(name)
 
-    tree_votes = [ tree_vote(n) for n in tree_names ]
+    def tree_vote_regressor(name):
+        return 'avg += {}(features, features_length); '.format(name)
+
+    forest_regressor_func = """float {function_name}(const {ctype} *features, int32_t features_length) {{
 
-    forest_func = """int32_t {function_name}(const {ctype} *features, int32_t features_length) {{
+        float avg = 0;
+
+        {tree_predictions}
+        
+        return avg/{n_classes};
+    }}
+    """.format(**{
+      'function_name': name,
+      'n_classes': n_classes,
+      'tree_predictions': '\n    '.join([ tree_vote_regressor(n) for n in tree_names ]),
+      'ctype': ctype,
+    })
+
+    forest_classifier_func = """int32_t {function_name}(const {ctype} *features, int32_t features_length) {{
 
         int32_t votes[{n_classes}] = {{0,}};
         int32_t _class = -1;
@@ -272,15 +305,22 @@ def tree_vote(name):
     """.format(**{
       'function_name': name,
       'n_classes': n_classes,
-      'tree_predictions': '\n    '.join(tree_votes),
+      'tree_predictions': '\n    '.join([ tree_vote_classifier(n) for n in tree_names ]),
       'ctype': ctype,
     })
+
+    return_type = 'int32_t'
+    forest_func = forest_classifier_func
     
-    tree_funcs = [tree_func(n, r) for n,r in zip(tree_names, roots)]
+    if not classifier:
+        return_type = 'float'
+        forest_func = forest_regressor_func
+
+    tree_funcs = [tree_func(n, r, return_type=return_type) for n,r in zip(tree_names, roots)]
 
     return '\n\n'.join(tree_funcs + [forest_func])
 
-def generate_c_forest(forest, name='myclassifier', dtype='float'):
+def generate_c_forest(forest, name='myclassifier', dtype='float', classifier=True):
     nodes, roots = forest
 
     nodes_name = name+'_nodes'
@@ -305,7 +345,7 @@ def generate_c_forest(forest, name='myclassifier', dtype='float'):
     #include <eml_trees.h>
     """
 
-    inline = generate_c_inlined(forest, name+'_predict', dtype=dtype)
+    inline = generate_c_inlined(forest, name+'_predict', dtype=dtype, classifier=classifier)
 
     return '\n\n'.join([head, nodes_c, tree_roots, forest_struct, inline]) 
 
@@ -317,12 +357,21 @@ def __init__(self, estimator, classifier, dtype='float'):
 
         self.dtype = dtype
 
+        kind = type(estimator).__name__
+        leaf = 'argmax'
+        self.is_classifier = True
+        out_dtype = "int"
+        if 'Regressor' in kind:
+            leaf = 'value'
+            self.is_classifier = False
+            out_dtype = "float"
+
         if hasattr(estimator, 'estimators_'):
             trees = [ e.tree_ for e in estimator.estimators_]
         else:
             trees = [ estimator.tree_ ]
 
-        self.forest_ = flatten_forest(trees)
+        self.forest_ = flatten_forest(trees, leaf=leaf)
         self.forest_ = remove_duplicate_leaves(self.forest_)
 
         if classifier == 'pymodule':
@@ -339,19 +388,27 @@ def __init__(self, estimator, classifier, dtype='float'):
 
         elif classifier == 'loadable':
             name = 'mytree'
-            func = 'eml_trees_predict(&{}, values, length)'.format(name)
+
+            if self.is_classifier:
+                func = 'eml_trees_predict(&{}, values, length)'.format(name)
+            else:
+                func = 'eml_trees_regress1(&{}, values, length)'.format(name)
             code = self.save(name=name)
-            self.classifier_ = common.CompiledClassifier(code, name=name, call=func)
+            self.classifier_ = common.CompiledClassifier(code, name=name, call=func, out_dtype=out_dtype)
         elif classifier == 'inline':
             name = 'myinlinetree'
             func = '{}_predict(values, length)'.format(name)
             code = self.save(name=name)
-            self.classifier_ = common.CompiledClassifier(code, name=name, call=func)
+            self.classifier_ = common.CompiledClassifier(code, name=name, call=func, out_dtype=out_dtype)
         else:
             raise ValueError("Unsupported classifier method '{}'".format(classifier))
 
     def predict(self, X):
-        predictions = self.classifier_.predict(X)
+        if self.is_classifier:
+            predictions = self.classifier_.predict(X)
+        else:
+            predictions = self.classifier_.regress(X)            
+
         return predictions
 
     def save(self, name=None, file=None):
@@ -361,7 +418,7 @@ def save(self, name=None, file=None):
             else:
                 name = os.path.splitext(os.path.basename(file))[0]
 
-        code = generate_c_forest(self.forest_, name, dtype=self.dtype)
+        code = generate_c_forest(self.forest_, name, dtype=self.dtype, classifier=self.is_classifier)
         if file:
             with open(file, 'w') as f:
                 f.write(code)


Test patch:
diff --git a/emlearn/eml_test.h b/emlearn/eml_test.h
index 3bee2ed..dc65062 100644
--- a/emlearn/eml_test.h
+++ b/emlearn/eml_test.h
@@ -9,7 +9,7 @@
 #include "eml_common.h"
 
 #ifndef EML_CSV_BUFFER_LENGTH
-#define EML_CSV_BUFFER_LENGTH 1024
+#define EML_CSV_BUFFER_LENGTH 2*1024
 #endif
 #ifndef EML_CSV_VALUES_LENGTH
 #define EML_CSV_VALUES_LENGTH 256
diff --git a/test/test_trees.py b/test/test_trees.py
index 64b63ee..0cb53a4 100644
--- a/test/test_trees.py
+++ b/test/test_trees.py
@@ -4,8 +4,9 @@
 import numpy.testing
 
 from sklearn import datasets
-from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
-from sklearn.tree import DecisionTreeClassifier
+from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
+from sklearn.ensemble import ExtraTreesClassifier, ExtraTreesRegressor
+from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor
 
 import emlearn
 import pytest
@@ -13,23 +14,35 @@
 random = numpy.random.randint(0, 1000)
 print('random_state={}'.format(random))
 
-MODELS = {
+CLASSIFICATION_MODELS = {
     'RFC': RandomForestClassifier(n_estimators=10, random_state=random),
     'ETC': ExtraTreesClassifier(n_estimators=10, random_state=random),
     'DTC': DecisionTreeClassifier(random_state=random),
 }
-DATASETS = {
+
+REGRESSION_MODELS = {
+    'RFR': RandomForestRegressor(n_estimators=10, random_state=random),
+    'ERR': ExtraTreesRegressor(n_estimators=10, random_state=random),
+    'DTR': DecisionTreeRegressor(random_state=random),
+}
+
+CLASSIFICATION_DATASETS = {
     'binary': datasets.make_classification(n_classes=2, n_samples=100, random_state=random),
     '5way': datasets.make_classification(n_classes=5, n_informative=5, n_samples=100, random_state=random),
 }
+
+REGRESSION_DATASETS = {
+    '1out': datasets.make_regression(n_targets=1, n_samples=100, random_state=random),
+}
+
 METHODS = ['pymodule', 'loadable', 'inline']
 
-@pytest.mark.parametrize("data", DATASETS.keys())
-@pytest.mark.parametrize("model", MODELS.keys())
+@pytest.mark.parametrize("data", CLASSIFICATION_DATASETS.keys())
+@pytest.mark.parametrize("model", CLASSIFICATION_MODELS.keys())
 @pytest.mark.parametrize("method", METHODS)
-def test_trees_sklearn_predict(data, model, method):
-    X, y = DATASETS[data]
-    estimator = MODELS[model]
+def test_trees_sklearn_classifier_predict(data, model, method):
+    X, y = CLASSIFICATION_DATASETS[data]
+    estimator = CLASSIFICATION_MODELS[model]
 
     estimator.fit(X, y)
     cmodel = emlearn.convert(estimator, method=method)
@@ -39,6 +52,21 @@ def test_trees_sklearn_predict(data, model, method):
 
     numpy.testing.assert_equal(pred_c, pred_original)
 
+@pytest.mark.parametrize("data", REGRESSION_DATASETS.keys())
+@pytest.mark.parametrize("model", REGRESSION_MODELS.keys())
+@pytest.mark.parametrize("method", METHODS)
+def test_trees_sklearn_regressor_predict(data, model, method):
+    X, y = REGRESSION_DATASETS[data]
+    estimator = REGRESSION_MODELS[model]
+
+    estimator.fit(X, y)
+    cmodel = emlearn.convert(estimator, method=method)
+
+    pred_original = estimator.predict(X[:5])
+    pred_c = cmodel.predict(X[:5])
+
+    numpy.testing.assert_allclose(pred_c, pred_original, rtol=1e-3, atol=2)
+
 
 def test_deduplicate_single_tree():
     nodes = [