Start of I-Vector implementation (wip)

Added tests with expected results.

bob/learn/em/__init__.py

@@ -2,6 +2,7 @@ import bob.extension
 
 from .factor_analysis import ISVMachine, JFAMachine
 from .gmm import GMMMachine, GMMStats
+from .ivector import IVectorMachine
 from .kmeans import KMeansMachine
 from .linear_scoring import linear_scoring  # noqa: F401
 from .wccn import WCCN
@@ -30,6 +31,13 @@ def __appropriate__(*args):
 
 
 __appropriate__(
-    KMeansMachine, GMMMachine, GMMStats, WCCN, Whitening, ISVMachine, JFAMachine
+    KMeansMachine,
+    GMMMachine,
+    GMMStats,
+    IVectorMachine,
+    WCCN,
+    Whitening,
+    ISVMachine,
+    JFAMachine,
 )
 __all__ = [_ for _ in dir() if not _.startswith("_")]

bob/learn/em/ivector.py

+#!/usr/bin/env python
+# @author: Yannick Dayer <yannick.dayer@idiap.ch>
+# @date: Fri 06 May 2022 14:18:25 UTC+02
+
+import numpy as np
+
+from sklearn.base import BaseEstimator
+
+from bob.learn.em import GMMMachine, GMMStats
+
+
+class IVectorStats:
+    """This class is a container for the statistics of an IVectorMachine."""
+
+    def __init__(self) -> None:
+        self.nij_sigma_wij2 = {}
+        self.fnorm_sigma_wij = {}
+        self.snorm = np.zeros(shape=(self.dim_c * self.dim_d,), dtype=float)
+        self.n = np.zeros(shape=(self.dim_c,), dtype=float)
+
+
+class IVectorMachine(BaseEstimator):
+    """Trains and projects data using I-Vector."""
+
+    def __init__(
+        self,
+        ubm: GMMMachine,
+        dim_t: int = 2,
+        convergence_threshold: float = 1e-5,
+        max_iterations: int = 25,
+        **kwargs
+    ) -> None:
+        """Initializes the IVectorMachine object.
+
+        Parameters
+        ----------
+        ubm:
+            The Universal Background Model.
+        dim_t:
+            The dimension of the i-vector.
+        """
+
+        super().__init__(**kwargs)
+        self.ubm = ubm
+        self.dim_t = dim_t
+        self.convergence_threshold = convergence_threshold
+        self.max_iterations = max_iterations
+        # TODO: add params
+        # self.compute_likelihood = compute_likelihood
+        # self.sigma_update = sigma_update
+        # self.variance_floor = variance_floor
+
+    def initialize(self, features: np.ndarray) -> None:
+        """Initializes the I-Vector parameters at fit time."""
+        # TODO implement
+        # self.
+        pass
+
+    def e_step(self, data: np.ndarray) -> GMMStats:
+        """Computes the expectation step of the e-m algorithm."""
+        # n_samples = len(data)
+        # self.m_acc_Nij_Sigma_wij2  = {}
+        # self.m_acc_Fnorm_Sigma_wij = {}
+        # self.m_acc_Snorm = numpy.zeros(shape=(self.m_dim_c*self.m_dim_d,), dtype=numpy.float64)
+        # self.m_N = numpy.zeros(shape=(self.m_dim_c,), dtype=numpy.float64)
+
+        # for c in range(self.m_dim_c):
+        #   self.m_acc_Nij_Sigma_wij2[c]  = numpy.zeros(shape=(self.m_dim_t,self.m_dim_t), dtype=numpy.float64)
+        #   self.m_acc_Fnorm_Sigma_wij[c] = numpy.zeros(shape=(self.m_dim_d,self.m_dim_t), dtype=numpy.float64)
+
+        # for n in range(n_samples):
+        #   Nij = data[n].n
+        #   Fij = data[n].sum_px
+        #   Sij = data[n].sum_pxx
+
+        #   # Estimate latent variables
+        #   TtSigmaInv_Fnorm = machine.__compute_TtSigmaInvFnorm__(data[n])
+        #   I_TtSigmaInvNT = machine.__compute_Id_TtSigmaInvT__(data[n])
+
+        #   Fnorm = numpy.zeros(shape=(self.m_dim_c*self.m_dim_d,), dtype=numpy.float64)
+        #   Snorm = numpy.zeros(shape=(self.m_dim_c*self.m_dim_d,), dtype=numpy.float64)
+
+        #   # Compute normalized statistics
+        #   for c in range(self.m_dim_c):
+        #     start            = c*self.m_dim_d
+        #     end              = (c+1)*self.m_dim_d
+
+        #     Fc               = Fij[c,:]
+        #     Sc               = Sij[c,:]
+        #     mc               = self.m_meansupervector[start:end]
+
+        #     Fc_mc            = Fc * mc
+        #     Nc_mc_mcT        = Nij[c] * mc * mc
+
+        #     Fnorm[start:end] = Fc - Nij[c] * mc
+        #     Snorm[start:end] = Sc - (2 * Fc_mc) + Nc_mc_mcT
+
+        #   # Latent variables
+        #   I_TtSigmaInvNT_inv = numpy.linalg.inv(I_TtSigmaInvNT)
+        #   E_w_ij             = numpy.dot(I_TtSigmaInvNT_inv, TtSigmaInv_Fnorm)
+        #   E_w_ij2            = I_TtSigmaInvNT_inv + numpy.outer(E_w_ij, E_w_ij)
+
+        #   # Do the accumulation for each component
+        #   self.m_acc_Snorm   = self.m_acc_Snorm + Snorm    # (dim_c*dim_d)
+        #   for c in range(self.m_dim_c):
+        #     start            = c*self.m_dim_d
+        #     end              = (c+1)*self.m_dim_d
+        #     current_Fnorm    = Fnorm[start:end]            # (dim_d)
+        #     self.m_acc_Nij_Sigma_wij2[c]  = self.m_acc_Nij_Sigma_wij2[c] + Nij[c] * E_w_ij2                    # (dim_t, dim_t)
+        #     self.m_acc_Fnorm_Sigma_wij[c] = self.m_acc_Fnorm_Sigma_wij[c] + numpy.outer(current_Fnorm, E_w_ij) # (dim_d, dim_t)
+        #     self.m_N[c]                   = self.m_N[c] + Nij[c]
+
+        # TODO convert
+        pass
+
+    def m_step(self, stats: GMMStats) -> None:
+        """Updates the Machine with the maximization step of the e-m algorithm."""
+        # A = self.m_acc_Nij_Sigma_wij2
+
+        # T = numpy.zeros(shape=(self.m_dim_c*self.m_dim_d,self.m_dim_t), dtype=numpy.float64)
+        # Told = machine.t
+        # if self.m_sigma_update:
+        #   sigma = numpy.zeros(shape=self.m_acc_Snorm.shape, dtype=numpy.float64)
+        # for c in range(self.m_dim_c):
+        #   start = c*self.m_dim_d;
+        #   end   = (c+1)*self.m_dim_d;
+        #   # T update
+        #   A     = self.m_acc_Nij_Sigma_wij2[c].transpose()
+        #   B     = self.m_acc_Fnorm_Sigma_wij[c].transpose()
+        #   if numpy.array_equal(A, numpy.zeros(A.shape)):
+        #     X = numpy.zeros(shape=(self.m_dim_t,self.m_dim_d), dtype=numpy.float64)
+        #   else:
+        #     X = numpy.linalg.solve(A, B)
+        #   T[start:end,:] = X.transpose()
+        #   # Sigma update
+        #   if self.m_sigma_update:
+        #     Told_c           = Told[start:end,:].transpose()
+        #     # warning: Use of the new T estimate! (revert second next line if you don't want that)
+        #     Fnorm_Ewij_Tt    = numpy.diag(numpy.dot(self.m_acc_Fnorm_Sigma_wij[c], X))
+        #     #Fnorm_Ewij_Tt = numpy.diag(numpy.dot(self.m_acc_Fnorm_Sigma_wij[c], Told_c))
+        #     sigma[start:end] = (self.m_acc_Snorm[start:end] - Fnorm_Ewij_Tt) / self.m_N[c]
+        # TODO convert
+        pass
+
+    def fit(self, data: np.ndarray) -> "IVectorMachine":
+        """Trains the IVectorMachine."""
+        # TODO implement
+        return self

bob/learn/em/test/test_ivector.py

+#!/usr/bin/env python
+# @author: Yannick Dayer <yannick.dayer@idiap.ch>
+# @date: Fri 06 May 2022 12:59:21 UTC+02
+
+import numpy as np
+
+from bob.learn.em import GMMMachine, GMMStats, IVectorMachine
+
+
+def test_ivector_machine_projection():
+    # Preset IVector Machine parameters
+    t = np.array(
+        [[1, 2], [4, 1], [0, 3], [5, 8], [7, 10], [11, 1]], dtype=float
+    )
+    sigma = np.array([1, 2, 1, 3, 2, 4], dtype=float)
+
+    # Create the UBM and set its values manually
+    ubm = GMMMachine(n_gaussians=2)
+    ubm.weights = np.array([0.4, 0.6], dtype=float)
+    ubm.means = np.array([[1, 7, 4], [4, 5, 3]], dtype=float)
+    ubm.variances = np.array([[0.5, 1.0, 1.5], [1.0, 1.5, 2.0]], dtype=float)
+
+    # Manually create a feature projected on the UBM
+    gmm_projection = GMMStats(ubm.n_gaussians, ubm.means.shape[-1])
+    gmm_projection.t = 1
+    gmm_projection.n = np.array([0.4, 0.6], dtype=float)
+    gmm_projection.sum_px = np.array([[1, 2, 3], [2, 4, 3]], dtype=float)
+    gmm_projection.sum_pxx = np.array([[10, 20, 30], [40, 50, 60]], dtype=float)
+
+    machine = IVectorMachine(ubm=ubm, dim_t=2)
+    machine.t = t
+    machine.sigma = sigma
+
+    # Reference from C++ implementation
+    ivector_projection_ref = np.array([-0.04213415, 0.21463343])
+    ivector_projection = machine.project(gmm_projection)
+    np.testing.assert_almost_equal(
+        ivector_projection_ref, ivector_projection, decimal=5
+    )
+
+
+def test_ivector_machine_training():
+    assert False, "Not implemented"