docs: def compute_G_from_H( H, mdimG = None, mode_inv = "svd" ):
      """
      Parameters
      ----------
      H : the inverse R. Metric
      if mode_inv == 'svd':
         also returns Hvv, Hsvals, Gsvals the (transposed) eigenvectors of
         H and the singular values of H and G
      if mdimG < H.shape[2]:
         G.shape = [ n_samples, mdimG, mdimG ] with n_samples = H.shape[0]
  
      Notes
      ------
      currently Hvv, Hsvals are n_dim = H.shape[2], and Gsvals, G are mdimG
      (This contradicts the documentation of riemann_metric which states
      that riemann_metric and h_dual_metric have the same dimensions)
  
      See the notes in RiemannMetric
      """
      n_samples = H.shape[0]
      n_dim = H.shape[2]
      if mode_inv is 'svd':
          Huu, Hsvals, Hvv = np.linalg.svd( H )
          if mdimG is None:
              Gsvals = 1./Hsvals
              G = np.zeros((n_samples, n_dim, n_dim))
              for i in np.arange(n_samples):
                  G[i,:,:] = np.dot(Huu[i,:,:], np.dot( np.diag(Gsvals[i,:]), Hvv[i,:,:]))
          elif mdimG < n_dim:
              Gsvals[:,:mdimG] = 1./Hsvals[:,:mdimG]
              Gsvals[:,mdimG:] = 0.
              # this can be redone with np.einsum() but it's barbaric
              G = np.zeros((n_samples, mdimG, mdimG))
              for i in np.arange(n_samples):
                  G[i,:,:mdimG] = np.dot(Huu[i,:,mdimG], np.dot( np.diag(Gsvals[i,:mdimG]), Hvv[i,:,:mdimG]))
          else:
              raise ValueError('mdimG must be <= H.shape[1]')
          return G, Hvv, Hsvals, Gsvals
      else:
          riemann_metric = np.linalg.inv(h_dual_metric)
          return riemann_metric, None, None, None