class_docstr: Implements of the Sasaki metric on the tangent bundle TM of a Riem. manifold M.
  
  The Sasaki metric is characterized by the following three properties:
  
  * the canonical projection of TM becomes a Riemannian submersion,
  * parallel vector fields along curves are orthogonal to their fibres, and
  * its restriction to any tangent space is Euclidean.
  
  Geodesic computations are realized via a discrete formulation of the
  geodesic equation on TM that involve geodesics, parallel translation, and
  the curvature tensor on the base manifold M (see [1]_ for details).
  However, as the implemented energy in the discrete-geodesics-optimization
  as well as the approximations of its gradient slightly differ from those
  proposed in [1]_, we also refer to [2]_ for additional details.
  
  Parameters
  ----------
  space : Manifold
      Tangent bundle.
  n_jobs: int
      Number of jobs for parallel computing.
      Optional, default: 1.
  n_steps : int
      Number of discrete time steps.
      Optional, default: 3.
  
  References
  ----------
  .. [1] Muralidharan, P., & Fletcher, P. T. "Sasaki metrics for analysis of
      longitudinal data on manifolds", IEEE CVPR 2012, pp. 1027-1034
      https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4270017/
  .. [2] Nava-Yazdani, E., Hanik, M., Ambellan, F., & von Tycowicz, C.
      "On Gradient Formulas in an Algorithm for the Logarithm of the Sasaki
      Metric", Technical Report Zuse-Institut Berlin, 2022
      https://nbn-resolving.org/urn/resolver.pl?urn:nbn:de:0297-zib-87174
class_docstr_tok: ["Implements","of","the","Sasaki","metric","on","the","tangent","bundle","TM","of","a","Riem",".","manifold","M.","The","Sasaki","metric","is","characterized","by","the","following","three","properties",":","*","the","canonical","projection","of","TM","becomes","a","Riemannian","submersion",",","*","parallel","vector","fields","along","curves","are","orthogonal","to","their","fibres",",","and","*","its","restriction","to","any","tangent","space","is","Euclidean",".","Geodesic","computations","are","realized","via","a","discrete","formulation","of","the","geodesic","equation","on","TM","that","involve","geodesics",",","parallel","translation",",","and","the","curvature","tensor","on","the","base","manifold","M","(","see","[","1","]","_","for","details",")",".","However",",","as","the","implemented","energy","in","the","discrete-geodesics-optimization","as","well","as","the","approximations","of","its","gradient","slightly","differ","from","those","proposed","in","[","1","]","_",",","we","also","refer","to","[","2","]","_","for","additional","details",".","Parameters","--","--","--","--","--","space",":","Manifold","Tangent","bundle",".","n_jobs",":","int","Number","of","jobs","for","parallel","computing",".","Optional",",","default",":","1.","n_steps",":","int","Number","of","discrete","time","steps",".","Optional",",","default",":","3",".","References","--","--","--","--","--","..","[","1","]","Muralidharan",",","P.",",","&","Fletcher",",","P.","T.","``","Sasaki","metrics","for","analysis","of","longitudinal","data","on","manifolds","''",",","IEEE","CVPR","2012",",","pp",".","1027-1034","https",":","\/\/www.ncbi.nlm.nih.gov\/pmc\/articles\/PMC4270017\/","..","[","2","]","Nava-Yazdani",",","E.",",","Hanik",",","M.",",","Ambellan",",","F.",",","&","von","Tycowicz",",","C.","``","On","Gradient","Formulas","in","an","Algorithm","for","the","Logarithm","of","the","Sasaki","Metric","''",",","Technical","Report","Zuse-Institut","Berlin",",","2022","https",":","\/\/nbn-resolving.org\/urn\/resolver.pl","?","urn",":","nbn",":","de:0297-zib-87174"]
docstr: Compute Sakai geodesic employing a variational time discretization.
  
  Parameters
  ----------
  end_points : array-like, shape=[..., 2, M.shape]
      Points in the tangent bundle TM of manifold M.
  initial_points : array-like, shape=[..., 2, M.shape]
      Points in the tangent bundle TM of manifold M.
  
  Returns
  -------
  geodesic : array-like, shape=[..., n_steps + 1, 2, M.shape]
      Discrete geodesics of form x(s)=(p(s), u(s)) in Sasaki metric
      connecting initial_point = x(0) and end_point = x(1).
docstr_tok: ["Compute","Sakai","geodesic","employing","a","variational","time","discretization",".","Parameters","--","--","--","--","--","end_points",":","array-like",",","shape=","[","...",",","2",",","M.shape","]","Points","in","the","tangent","bundle","TM","of","manifold","M.","initial_points",":","array-like",",","shape=","[","...",",","2",",","M.shape","]","Points","in","the","tangent","bundle","TM","of","manifold","M.","Returns","--","--","--","-","geodesic",":","array-like",",","shape=","[","...",",","n_steps","+","1",",","2",",","M.shape","]","Discrete","geodesics","of","form","x","(","s",")","=","(","p","(","s",")",",","u","(","s",")",")","in","Sasaki","metric","connecting","initial_point","=","x","(","0",")","and","end_point","=","x","(","1",")","."]
code: def geodesic_discrete(self, initial_point, end_point):
      """Compute Sakai geodesic employing a variational time discretization.
  
      Parameters
      ----------
      end_points : array-like, shape=[..., 2, M.shape]
          Points in the tangent bundle TM of manifold M.
      initial_points : array-like, shape=[..., 2, M.shape]
          Points in the tangent bundle TM of manifold M.
  
      Returns
      -------
      geodesic : array-like, shape=[..., n_steps + 1, 2, M.shape]
          Discrete geodesics of form x(s)=(p(s), u(s)) in Sasaki metric
          connecting initial_point = x(0) and end_point = x(1).
      """
      metric = self._space.space.metric
      par_trans = metric.parallel_transport
  
      def _grad(pu, i_pt, e_pt):
          """Gradient of discrete geodesic energy."""
          pu = gs.vstack(
              [
                  gs.expand_dims(i_pt, axis=0),
                  pu,
                  gs.expand_dims(e_pt, axis=0),
              ]
          )
          p, u = self._space._unstack(pu)
  
          p1, p2, p3 = p[:-2], p[1:-1], p[2:]
          u1, u2, u3 = u[:-2], u[1:-1], u[2:]
  
          eps = 1 / self.n_steps
  
          v2 = metric.log(p3, p2) / eps
          w2 = (par_trans(u3, p3, end_point=p2) - u2) / eps
  
          gp = (metric.log(p3, p2) + metric.log(p1, p2)) / (
              2 * eps**2
          ) - metric.curvature(u2, w2, v2, p2)
  
          gu = (
              par_trans(u3, p3, end_point=p2)
              - 2 * u2
              + par_trans(u1, p1, end_point=p2)
          ) / eps**2
  
          return -self._space._stack(gp, gu) * eps
  
      def _geodesic_discrete_single(initial_point, end_point):
          """Calculate the discrete geodesic."""
          ijk = "ijk"[: self._space.space.point_ndim]
  
          def _scalarmul(scalar, point):
              return gs.einsum(f"...,...{ijk}->...{ijk}", scalar, point)
  
          p0, u0 = initial_point[0], initial_point[1]
          pL, uL = end_point[0], end_point[1]
  
          v = metric.log(pL, p0)
          s = gs.linspace(0.0, 1.0, self.n_steps + 1)[1:-1]
  
          p_ini = metric.exp(gs.einsum(f"p, {ijk}->p{ijk}", s, v), p0)
          u_ini = _scalarmul(
              (1.0 - s), par_trans(u0, p0, end_point=p_ini)
          ) + _scalarmul(s, par_trans(uL, pL, end_point=p_ini))
          pu_ini = self._space._stack(p_ini, u_ini)
  
          x = self._gradient_descent.minimize(
              pu_ini, initial_point, end_point, _grad, self.exp
          )
  
          return gs.vstack(
              [
                  gs.expand_dims(initial_point, axis=0),
                  x,
                  gs.expand_dims(end_point, axis=0),
              ]
          )
  
      is_batch = check_is_batch(
          self._space.point_ndim, initial_point, end_point
      )
      if not is_batch:
          return _geodesic_discrete_single(initial_point, end_point)
  
      initial_point, end_point = gs.broadcast_arrays(
          initial_point, end_point
      )
  
      with Parallel(
          n_jobs=min(self.n_jobs, len(end_point)), verbose=0
      ) as parallel:
          rslt = parallel(
              delayed(_geodesic_discrete_single)(i_pt, e_pt)
              for i_pt, e_pt in zip(initial_point, end_point)
          )
      return gs.array(rslt)
  
imports: from joblib import Parallel, delayed
  import geomstats.backend
  from geomstats.geometry.base import Manifold
  from geomstats.geometry.riemannian_metric import RiemannianMetric
  from geomstats.vectorization import check_is_batch
label_desc: Use image node_id 4 for calling the SasakiMetric obj's underlying member method code with example usage: obj.geodesic_discrete(initial_point, end_point) and returns: gs, unknown, gs, _geodesic_discrete_single, gs