docs: def triangle_address(fx, pt):
      '''
      triangle_address(FX, P) yields an address coordinate (t,r) for the point P in the triangle
      defined by the (3 x d)-sized coordinate matrix FX, in which each row of the matrix is the
      d-dimensional vector representing the respective triangle vertx for triangle [A,B,C]. The
      resulting coordinates (t,r) (0 <= t <= 1, 0 <= r <= 1) address the point P such that, if t gives
      the fraction of the angle from vector AB to vector AC that is made by the angle between vectors
      AB and AP, and r gives the fraction ||AP||/||AR|| where R is the point of intersection between
      lines AP and BC. If P is a (d x n)-sized matrix of points, then a (2 x n) matrix of addresses
      is returned.
      '''
      fx = np.asarray(fx)
      pt = np.asarray(pt)
      # The triangle vectors...
      ab = fx[1] - fx[0]
      ac = fx[2] - fx[0]
      bc = fx[2] - fx[1]
      ap = np.asarray([pt_i - a_i for (pt_i, a_i) in zip(pt, fx[0])])
      # get the unnormalized distance...
      r = np.sqrt((ap ** 2).sum(0))
      # now we can find the angle...
      unit = 1 - r.astype(bool)
      t0 = vector_angle(ab, ac)
      t = vector_angle(ap + [ab_i * unit for ab_i in ab], ab)
      sint = np.sin(t)
      sindt = np.sin(t0 - t)
      # finding r0 is tricker--we use this fancy formula based on the law of sines
      q0 = np.sqrt((bc ** 2).sum(0))          # B->C distance
      beta = vector_angle(-ab, bc)            # Angle at B
      sinGamma = np.sin(math.pi - beta - t0)
      sinBeta  = np.sin(beta)
      r0 = q0 * sinBeta * sinGamma / (sinBeta * sindt + sinGamma * sint)
      return np.asarray([t/t0, r/r0])