Repository: emlearn/emlearn (language: python)

Issue:
Mixture model does not support export to c file
Thank you for your awesome repo. I tried to apply Gaussian Mixture from Sklearn to emlearn and I have tried to convert emlearn to C code file. But I got the problem. It's raised "raise NotImplementedError("TODO implement save()")". Pls help me to solve it


Discussion / hints:
Hi @cnxdeveloper, thank you for reporting. I can confirm this issue. It seems that the save() function is not hooked up to the model code generation. Will try to look into it over the next days.

Fix patch:
diff --git a/emlearn/mixture.py b/emlearn/mixture.py
index 0f2e1e7..daad75c 100644
--- a/emlearn/mixture.py
+++ b/emlearn/mixture.py
@@ -1,148 +1,28 @@
 
-import os.path
-import os
-
-import numpy
-
-from . import common, cgen
-
-# Ref 
 """
-References
-
-https://github.com/scikit-learn/scikit-learn/blob/95119c13af77c76e150b753485c662b7c52a41a2/sklearn/mixture/_gaussian_mixture.py#L380
-
-https://github.com/scikit-learn/scikit-learn/blob/95119c13af77c76e150b753485c662b7c52a41a2/sklearn/mixture/_base.py
-
-_estimate_log_gaussian_prob
-
-log probability is used
-implementation depends on covariance type. Can be known at fit time
-
-Looks like log_det can be known at fit time
-one constant per component
+Gaussian Mixture models
 
-Stores
-means_ means of each component
-weights_ weights of each component
+Implements the GMM models from scikit-learn,
+and reuses some code for compile-time processing.
 
-precisions_cholesky_
-which is a Cholesky decomposition of the precision, the inverse of the covariance matrix
-
-for full. components, features, features
-for spherical. components
-for diag. components, features
-for tied, features, features 
-
-BaysianGaussianMixture
+https://github.com/scikit-learn/scikit-learn/blob/95119c13af77c76e150b753485c662b7c52a41a2/sklearn/mixture/_gaussian_mixture.py#L380
 https://github.com/scikit-learn/scikit-learn/blob/95119c13af77c76e150b753485c662b7c52a41a2/sklearn/mixture/_bayesian_mixture.py
+https://github.com/scikit-learn/scikit-learn/blob/95119c13af77c76e150b753485c662b7c52a41a2/sklearn/mixture/_base.py
 
-Seem to reuse _estimate_log_gaussian_prob from GMM
-Has an additional term log_lambda
-does not seem to depend on X?
-
+For other GMM implementations in C, see vlfeat
 
 https://www.vlfeat.org/api/gmm.html
-implements only diagonal covariance matrix
-
 https://github.com/vlfeat/vlfeat/blob/master/vl/gmm.c#L712
-
+Note: implements only diagonal covariance matrix
 """
 
-from sklearn.mixture._gaussian_mixture import _compute_log_det_cholesky
-from sklearn.utils.extmath import row_norms
-np = numpy
-
-
-def _estimate_log_gaussian_prob(X, means, precisions_chol, covariance_type):
-    """Estimate the log Gaussian probability.
-    Parameters
-    ----------
-    X : array-like of shape (n_samples, n_features)
-    means : array-like of shape (n_components, n_features)
-    precisions_chol : array-like
-        Cholesky decompositions of the precision matrices.
-        'full' : shape of (n_components, n_features, n_features)
-        'tied' : shape of (n_features, n_features)
-        'diag' : shape of (n_components, n_features)
-        'spherical' : shape of (n_components,)
-    covariance_type : {'full', 'tied', 'diag', 'spherical'}
-    Returns
-    -------
-    log_prob : array, shape (n_samples, n_components)
-    """
-    n_samples, n_features = X.shape
-    n_components, _ = means.shape
-    # det(precision_chol) is half of det(precision)
-    log_det = _compute_log_det_cholesky(
-        precisions_chol, covariance_type, n_features)
-
-    if covariance_type == 'full':
-        log_prob = np.empty((n_samples, n_components))
-
-        #print('mm', n_samples, n_features, n_components)
-
-        for i, x in enumerate(X):
-            #print('x', i, x)
-            for k, (mu, prec_chol) in enumerate(zip(means, precisions_chol)):
-
-                pp = 0.0
-   
-                for f in range(x.shape[0]):
-                    dot_m = 0.0
-                    dot_x = 0.0
-
-                    for p in range(prec_chol.shape[0]):
-                        dot_m += (mu[p] * prec_chol[p,f])
-                        dot_x += (x[p] * prec_chol[p,f])
-
-                    y = (dot_x - dot_m)
-                    pp += ( y * y ) 
-
-                #print('k', k, '\n', mu, '\n', prec_chol)
-                dot_x = np.dot(x, prec_chol)
-                dot_m = np.dot(mu, prec_chol)
-                y = dot_x - dot_m
-
-                #print('dot_x', dot_x)
-                #print('dot_m', dot_m)
-                #print('y', y)
-
-                p = np.sum(np.square(y), axis=0) # sum over features
-                #assert p == pp, (p, pp)
-                #print("log_prob", i, k, p)
-
-                log_prob[i, k] = p
-
-    elif covariance_type == 'tied':
-        log_prob = np.empty((n_samples, n_components))
-        for k, mu in enumerate(means):
-            y = np.dot(X, precisions_chol) - np.dot(mu, precisions_chol)
-            log_prob[:, k] = np.sum(np.square(y), axis=1)
-
-    elif covariance_type == 'diag':
-        precisions = precisions_chol ** 2
-        log_prob = (np.sum((means ** 2 * precisions), 1) -
-                    2. * np.dot(X, (means * precisions).T) +
-                    np.dot(X ** 2, precisions.T))
-
-    elif covariance_type == 'spherical':
-        precisions = precisions_chol ** 2
-        log_prob = (np.sum(means ** 2, 1) * precisions -
-                    2 * np.dot(X, means.T * precisions) +
-                    np.outer(row_norms(X, squared=True), precisions))
-
-
-    s = -.5 * (n_features * np.log(2 * np.pi) + log_prob) + log_det
-    #print('s', s, 'log_det\n', log_det)
-
-    return s
 
+import os.path
+import tempfile
 
-def init_model():
-    
+import numpy
 
-    pass
+from . import common, cgen
 
 
 def get_covariance_type(s):
@@ -247,7 +127,7 @@ def predict_one(x):
     y = [ predict_one(x) for x in numpy.array(X) ]
     return numpy.array(y)
 
-def build_executable_proba(wrapper, name='gmm'):
+def build_executable_proba(wrapper, out_dir, name='gmm'):
     n_components, n_features = wrapper._means.shape
 
     model_code = generate_code(wrapper, name=name)
@@ -293,7 +173,6 @@ def build_executable_proba(wrapper, name='gmm'):
     
 
     # Compile the xor.c example program
-    out_dir = './examples'
     src_path = os.path.join(out_dir, 'gmm.c')
 
     with open(src_path, 'w') as f:
@@ -304,7 +183,7 @@ def build_executable_proba(wrapper, name='gmm'):
 
     return bin_path
 
-def build_executable_score(wrapper, name='gmm'):
+def build_executable_score(wrapper, out_dir, name='gmm'):
     n_components, n_features = wrapper._means.shape
 
     model_code = generate_code(wrapper, name=name)
@@ -346,7 +225,6 @@ def build_executable_score(wrapper, name='gmm'):
     """
 
     # Compile
-    out_dir = './examples'
     src_path = os.path.join(out_dir, 'gmm.c')
 
     with open(src_path, 'w') as f:
@@ -358,6 +236,12 @@ def build_executable_score(wrapper, name='gmm'):
     return bin_path
 
 def convert_to_full(means, precisions_chol, covariance_type):
+    """
+    Convert the different covariance structures into "full" covariance matrix form.
+
+    This means only a single inference implementation is needed.
+    Likely less optimal in terms of run-time, but much simpler.
+    """
 
     n_components, n_features = means.shape
 
@@ -414,11 +298,11 @@ def get_log_weights(estimator):
     # https://github.com/scikit-learn/scikit-learn/blob/7e1e6d09bcc2eaeba98f7e737aac2ac782f0e5f1/sklearn/mixture/_bayesian_mixture.py#L768
     is_bayesian = hasattr(estimator, 'degrees_of_freedom_') 
     if is_bayesian:
-        # Remove `n_features * np.log(self.degrees_of_freedom_)` because the precision matrix is normalized
-        log_weights = log_weights - (0.5 * n_features * np.log(estimator.degrees_of_freedom_))
+        # Remove `n_features * numpy.log(self.degrees_of_freedom_)` because the precision matrix is normalized
+        log_weights = log_weights - (0.5 * n_features * numpy.log(estimator.degrees_of_freedom_))
 
-        log_lambda = n_features * np.log(2.0) + np.sum(
-            digamma(0.5* (estimator.degrees_of_freedom_ - np.arange(0, n_features)[:, np.newaxis])),
+        log_lambda = n_features * numpy.log(2.0) + numpy.sum(
+            digamma(0.5* (estimator.degrees_of_freedom_ - numpy.arange(0, n_features)[:, numpy.newaxis])),
             0,
         )
 
@@ -456,15 +340,11 @@ def __init__(self, estimator, classifier, dtype='float', verbose=0):
 
 
     def predict_proba(self, X):
-        #from sklearn.mixture._gaussian_mixture import _estimate_log_gaussian_prob
 
-        py_predictions = _estimate_log_gaussian_prob(X, self._means,
-                                    self._precisions_col, self._covariance_type)
+        with tempfile.TemporaryDirectory() as out_dir:
 
-        py_predictions += self._log_weights
-
-        bin_path = build_executable_proba(self)
-        c_predictions = predict(bin_path, X, verbose=self.verbose)
+            bin_path = build_executable_proba(self, out_dir=out_dir)
+            c_predictions = predict(bin_path, X, verbose=self.verbose)
 
         # XXX: note, this is actually log probabilities
         return c_predictions
@@ -478,8 +358,11 @@ def score_samples(self, X):
         from scipy.special import logsumexp
 
         if True:
-            bin_path = build_executable_score(self)
-            predictions = predict(bin_path, X, verbose=self.verbose)
+            with tempfile.TemporaryDirectory() as out_dir:
+
+                bin_path = build_executable_score(self, out_dir=out_dir)
+                predictions = predict(bin_path, X, verbose=self.verbose)
+
             score = predictions[:,0]
 
         else:
@@ -496,11 +379,10 @@ def save(self, name=None, file=None):
             else:
                 name = os.path.splitext(os.path.basename(file))[0]
 
-        code = "" # Implement
+        code = generate_code(self, name=name)
         if file:
             with open(file, 'w') as f:
                 f.write(code)
 
-        raise NotImplementedError("TODO implement save()")
         return code
 


Test patch:
diff --git a/test/test_mixture.py b/test/test_mixture.py
index fbf4c33..f57dba5 100644
--- a/test/test_mixture.py
+++ b/test/test_mixture.py
@@ -20,6 +20,7 @@
 
 import pytest
 
+here = os.path.dirname(__file__)
 random = numpy.random.randint(0, 1000)
 random = 1
 print('random_state={}'.format(random))
@@ -46,7 +47,6 @@
     #'loadable',
 ]
 
-
 @pytest.mark.parametrize("data", DATASETS.keys())
 @pytest.mark.parametrize("model", MODELS.keys())
 @pytest.mark.parametrize("method", METHODS)
@@ -70,3 +70,10 @@ def test_gaussian_mixture_equals_sklearn(data, model, method):
   
     numpy.testing.assert_equal(pred_c, pred_original)
 
+    # check that code can be generated
+    out_dir = os.path.join(here, 'out', 'gaussian_mixture', f'{data}_{model}_{method}.c')
+    save_path = os.path.join(out_dir, f'{data}_{model}_{method}.c')
+    if not os.path.exists(out_dir):
+        os.makedirs(out_dir)
+    cmodel.save(file=save_path, name='my_test_model')
+