code: def fda(X, y, p=2, reg=1e-16):
      
  
      mx = np.mean(X)
      X -= mx.reshape((1, -1))
  
      # data split between classes
      d = X.shape[1]
      xc = split_classes(X, y)
      nc = len(xc)
  
      p = min(nc - 1, p)
  
      Cw = 0
      for x in xc:
          Cw += np.cov(x, rowvar=False)
      Cw /= nc
  
      mxc = np.zeros((d, nc))
  
      for i in range(nc):
          mxc[:, i] = np.mean(xc[i])
  
      mx0 = np.mean(mxc, 1)
      Cb = 0
      for i in range(nc):
          Cb += (mxc[:, i] - mx0).reshape((-1, 1)) * \
              (mxc[:, i] - mx0).reshape((1, -1))
  
      w, V = linalg.eig(Cb, Cw + reg * np.eye(d))
  
      idx = np.argsort(w.real)
  
      Popt = V[:, idx[-p:]]
  
      def proj(X):
          return (X - mx.reshape((1, -1))).dot(Popt)
  
      return Popt, proj
docs: Fisher Discriminant Analysis
  
  
      Parameters
      ----------
      X : numpy.ndarray (n,d)
          Training samples
      y : np.ndarray (n,)
          labels for training samples
      p : int, optional
          size of dimensionnality reduction
      reg : float, optional
          Regularization term >0 (ridge regularization)
  
  
      Returns
      -------
      P : (d x p) ndarray
          Optimal transportation matrix for the given parameters
      proj : fun
          projection function including mean centering