code: def apply_M(self, ax, ay):
          
          jac = numpy.array(
              [[self.dx.dot(ax), self.dy.dot(ax)], [self.dx.dot(ay), self.dy.dot(ay)]]
          )
  
          # jacs and J are of shape (2, 2, k). M must be of the same shape and
          # contain the result of the k 2x2 dot products. Perhaps there's a
          # dot() for this.
          M = numpy.einsum("ijl,jkl->ikl", jac, self.J)
          # M = numpy.array([
          #     [
          #         jac[0][0]*self.J[0][0] + jac[0][1]*self.J[1][0],
          #         jac[0][0]*self.J[0][1] + jac[0][1]*self.J[1][1],
          #     ],
          #     [
          #         jac[1][0]*self.J[0][0] + jac[1][1]*self.J[1][0],
          #         jac[1][0]*self.J[0][1] + jac[1][1]*self.J[1][1],
          #     ],
          #     ])
  
          # One could use
          #
          #     M = numpy.moveaxis(M, -1, 0)
          #     _, sigma, _ = numpy.linalg.svd(M)
          #
          # but computing the singular values explicitly via
          # <https://scicomp.stackexchange.com/a/14103/3980> is faster and more
          # explicit.
          a = (M[0, 0] + M[1, 1]) / 2
          b = (M[0, 0] - M[1, 1]) / 2
          c = (M[1, 0] + M[0, 1]) / 2
          d = (M[1, 0] - M[0, 1]) / 2
  
          return a, b, c, d
docs: Linear operator that converts ax, ay to abcd.