nodes: 0, module; 1, function_definition; 2, function_name:squad; 3, parameters; 4, block; 5, identifier:R_in; 6, identifier:t_in; 7, identifier:t_out; 8, expression_statement; 9, if_statement; 10, comment:# This list contains an index for each `t_out` such that; 11, comment:# t_in[i-1] <= t_out < t_in[i]; 12, comment:# Note that `side='right'` is much faster in my tests; 13, comment:# i_in_for_out = t_in.searchsorted(t_out, side='left'); 14, comment:# np.clip(i_in_for_out, 0, len(t_in) - 1, out=i_in_for_out); 15, expression_statement; 16, comment:# Now, for each index `i` in `i_in`, we need to compute the; 17, comment:# interpolation "coefficients" (`A_i`, `B_ip1`).; 18, comment:#; 19, comment:# I previously tested an explicit version of the loops below,; 20, comment:# comparing `stride_tricks.as_strided` with explicit; 21, comment:# implementation via `roll` (as seen here).  I found that the; 22, comment:# `roll` was significantly more efficient for simple calculations,; 23, comment:# though the difference is probably totally washed out here.  In; 24, comment:# any case, it might be useful to test again.; 25, comment:#; 26, expression_statement; 27, expression_statement; 28, comment:# Correct the first and last A time steps, and last two B time steps.  We extend R_in with the following wrap-around; 29, comment:# values:; 30, comment:# R_in[0-1] = R_in[0]*(~R_in[1])*R_in[0]; 31, comment:# R_in[n+0] = R_in[-1] * (~R_in[-2]) * R_in[-1]; 32, comment:# R_in[n+1] = R_in[0] * (~R_in[-1]) * R_in[0]; 33, comment:#           = R_in[-1] * (~R_in[-2]) * R_in[-1] * (~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1]; 34, comment:#           = R_in[-1] * (~R_in[-2]) * R_in[-1] * (~R_in[-2]) * R_in[-1]; 35, comment:# A[i] = R_in[i] * np.exp((- np.log((~R_in[i]) * R_in[i+1]); 36, comment:#                          + np.log((~R_in[i-1]) * R_in[i]) * ((t_in[i+1] - t_in[i]) / (t_in[i] - t_in[i-1])); 37, comment:#                          ) * 0.25); 38, comment:# A[0] = R_in[0] * np.exp((- np.log((~R_in[0]) * R_in[1]) + np.log((~R_in[0])*R_in[1]*(~R_in[0])) * R_in[0]) * 0.25); 39, comment:#      = R_in[0]; 40, expression_statement; 41, comment:# A[-1] = R_in[-1] * np.exp((- np.log((~R_in[-1]) * R_in[n+0]); 42, comment:#                          + np.log((~R_in[-2]) * R_in[-1]) * ((t_in[n+0] - t_in[-1]) / (t_in[-1] - t_in[-2])); 43, comment:#                          ) * 0.25); 44, comment:#       = R_in[-1] * np.exp((- np.log((~R_in[-1]) * R_in[n+0]) + np.log((~R_in[-2]) * R_in[-1])) * 0.25); 45, comment:#       = R_in[-1] * np.exp((- np.log((~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1]); 46, comment:#                           + np.log((~R_in[-2]) * R_in[-1])) * 0.25); 47, comment:#       = R_in[-1] * np.exp((- np.log((~R_in[-2]) * R_in[-1]) + np.log((~R_in[-2]) * R_in[-1])) * 0.25); 48, comment:#       = R_in[-1]; 49, expression_statement; 50, comment:# B[i] = R_in[i+1] * np.exp((np.log((~R_in[i+1]) * R_in[i+2]) * ((t_in[i+1] - t_in[i]) / (t_in[i+2] - t_in[i+1])); 51, comment:#                            - np.log((~R_in[i]) * R_in[i+1])) * -0.25); 52, comment:# B[-2] = R_in[-1] * np.exp((np.log((~R_in[-1]) * R_in[0]) * ((t_in[-1] - t_in[-2]) / (t_in[0] - t_in[-1])); 53, comment:#                            - np.log((~R_in[-2]) * R_in[-1])) * -0.25); 54, comment:#       = R_in[-1] * np.exp((np.log((~R_in[-1]) * R_in[0]) - np.log((~R_in[-2]) * R_in[-1])) * -0.25); 55, comment:#       = R_in[-1] * np.exp((np.log((~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1]); 56, comment:#                            - np.log((~R_in[-2]) * R_in[-1])) * -0.25); 57, comment:#       = R_in[-1] * np.exp((np.log((~R_in[-2]) * R_in[-1]) - np.log((~R_in[-2]) * R_in[-1])) * -0.25); 58, comment:#       = R_in[-1]; 59, expression_statement; 60, comment:# B[-1] = R_in[0]; 61, comment:# B[-1] = R_in[0] * np.exp((np.log((~R_in[0]) * R_in[1]) - np.log((~R_in[-1]) * R_in[0])) * -0.25); 62, comment:#       = R_in[-1] * (~R_in[-2]) * R_in[-1]; 63, comment:#         * np.exp((np.log((~(R_in[-1] * (~R_in[-2]) * R_in[-1])) * R_in[-1] * (~R_in[-2]) * R_in[-1] * (~R_in[-2]) * R_in[-1]); 64, comment:#                  - np.log((~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1])) * -0.25); 65, comment:#       = R_in[-1] * (~R_in[-2]) * R_in[-1]; 66, comment:#         * np.exp((np.log(((~R_in[-1]) * R_in[-2] * (~R_in[-1])) * R_in[-1] * (~R_in[-2]) * R_in[-1] * (~R_in[-2]) * R_in[-1]); 67, comment:#                  - np.log((~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1])) * -0.25); 68, comment:#         * np.exp((np.log((~R_in[-2]) * R_in[-1]); 69, comment:#                  - np.log((~R_in[-2]) * R_in[-1])) * -0.25); 70, expression_statement; 71, comment:# Use the coefficients at the corresponding t_out indices to; 72, comment:# compute the squad interpolant; 73, comment:# R_ip1 = np.array(np.roll(R_in, -1)[i_in_for_out]); 74, comment:# R_ip1[-1] = R_in[-1]*(~R_in[-2])*R_in[-1]; 75, expression_statement; 76, expression_statement; 77, expression_statement; 78, expression_statement; 79, expression_statement; 80, expression_statement; 81, comment:# tau = (t_out - t_in[i_in_for_out]) / ((np.roll(t_in, -1) - t_in)[i_in_for_out]); 82, expression_statement; 83, return_statement; 84, comment:"""Spherical "quadrangular" interpolation of rotors with a cubic spline
  
      This is the best way to interpolate rotations.  It uses the analog
      of a cubic spline, except that the interpolant is confined to the
      rotor manifold in a natural way.  Alternative methods involving
      interpolation of other coordinates on the rotation group or
      normalization of interpolated values give bad results.  The
      results from this method are as natural as any, and are continuous
      in first and second derivatives.
  
      The input `R_in` rotors are assumed to be reasonably continuous
      (no sign flips), and the input `t` arrays are assumed to be
      sorted.  No checking is done for either case, and you may get
      silently bad results if these conditions are violated.
  
      This function simplifies the calling, compared to `squad_evaluate`
      (which takes a set of four quaternions forming the edges of the
      "quadrangle", and the normalized time `tau`) and `squad_vectorized`
      (which takes the same arguments, but in array form, and efficiently
      loops over them).
  
      Parameters
      ----------
      R_in: array of quaternions
          A time-series of rotors (unit quaternions) to be interpolated
      t_in: array of float
          The times corresponding to R_in
      t_out: array of float
          The times to which R_in should be interpolated
  
      """; 85, boolean_operator; 86, block; 87, assignment; 88, assignment; 89, assignment; 90, assignment; 91, assignment; 92, assignment; 93, assignment; 94, assignment; 95, assignment; 96, assignment; 97, assignment; 98, assignment; 99, assignment; 100, assignment; 101, identifier:R_out; 102, comparison_operator:R_in.size == 0; 103, comparison_operator:t_out.size == 0; 104, return_statement; 105, identifier:i_in_for_out; 106, binary_operator:t_in.searchsorted(t_out, side='right')-1; 107, identifier:A; 108, binary_operator:R_in * np.exp((- np.log((~R_in) * np.roll(R_in, -1))
                         + np.log((~np.roll(R_in, 1)) * R_in) * ((np.roll(t_in, -1) - t_in) / (t_in - np.roll(t_in, 1)))
                         ) * 0.25); 109, identifier:B; 110, binary_operator:np.roll(R_in, -1) * np.exp((np.log((~np.roll(R_in, -1)) * np.roll(R_in, -2))
                                      * ((np.roll(t_in, -1) - t_in) / (np.roll(t_in, -2) - np.roll(t_in, -1)))
                                      - np.log((~R_in) * np.roll(R_in, -1))) * -0.25); 111, subscript; 112, subscript; 113, subscript; 114, subscript; 115, subscript; 116, subscript; 117, subscript; 118, binary_operator:R_in[-1] * (~R_in[-2]) * R_in[-1]; 119, identifier:R_ip1; 120, call; 121, subscript; 122, binary_operator:R_in[-1]*(~R_in[-2])*R_in[-1]; 123, identifier:R_ip1; 124, call; 125, identifier:t_inp1; 126, call; 127, subscript; 128, binary_operator:t_in[-1] + (t_in[-1] - t_in[-2]); 129, identifier:tau; 130, binary_operator:(t_out - t_in[i_in_for_out]) / ((t_inp1 - t_in)[i_in_for_out]); 131, identifier:R_out; 132, call; 133, attribute; 134, integer:0; 135, attribute; 136, integer:0; 137, call; 138, call; 139, integer:1; 140, identifier:R_in; 141, call; 142, call; 143, call; 144, identifier:A; 145, integer:0; 146, identifier:R_in; 147, integer:0; 148, identifier:A; 149, unary_operator; 150, identifier:R_in; 151, unary_operator; 152, identifier:B; 153, unary_operator; 154, identifier:R_in; 155, unary_operator; 156, identifier:B; 157, unary_operator; 158, binary_operator:R_in[-1] * (~R_in[-2]); 159, subscript; 160, attribute; 161, argument_list; 162, identifier:R_ip1; 163, unary_operator; 164, binary_operator:R_in[-1]*(~R_in[-2]); 165, subscript; 166, attribute; 167, argument_list; 168, attribute; 169, argument_list; 170, identifier:t_inp1; 171, unary_operator; 172, subscript; 173, parenthesized_expression; 174, parenthesized_expression; 175, parenthesized_expression; 176, attribute; 177, argument_list; 178, identifier:R_in; 179, identifier:size; 180, identifier:t_out; 181, identifier:size; 182, attribute; 183, argument_list; 184, attribute; 185, argument_list; 186, attribute; 187, argument_list; 188, attribute; 189, argument_list; 190, attribute; 191, argument_list; 192, integer:1; 193, integer:1; 194, integer:2; 195, integer:1; 196, integer:1; 197, subscript; 198, parenthesized_expression; 199, identifier:R_in; 200, unary_operator; 201, identifier:np; 202, identifier:roll; 203, identifier:R_in; 204, unary_operator; 205, integer:1; 206, subscript; 207, parenthesized_expression; 208, identifier:R_in; 209, unary_operator; 210, identifier:np; 211, identifier:array; 212, subscript; 213, identifier:np; 214, identifier:roll; 215, identifier:t_in; 216, unary_operator; 217, integer:1; 218, identifier:t_in; 219, unary_operator; 220, binary_operator:t_in[-1] - t_in[-2]; 221, binary_operator:t_out - t_in[i_in_for_out]; 222, subscript; 223, identifier:np; 224, identifier:squad_vectorized; 225, identifier:tau; 226, subscript; 227, subscript; 228, subscript; 229, identifier:R_ip1; 230, identifier:np; 231, identifier:array; 232, tuple; 233, keyword_argument; 234, identifier:t_in; 235, identifier:searchsorted; 236, identifier:t_out; 237, keyword_argument; 238, identifier:np; 239, identifier:exp; 240, binary_operator:(- np.log((~R_in) * np.roll(R_in, -1))
                         + np.log((~np.roll(R_in, 1)) * R_in) * ((np.roll(t_in, -1) - t_in) / (t_in - np.roll(t_in, 1)))
                         ) * 0.25; 241, identifier:np; 242, identifier:roll; 243, identifier:R_in; 244, unary_operator; 245, identifier:np; 246, identifier:exp; 247, binary_operator:(np.log((~np.roll(R_in, -1)) * np.roll(R_in, -2))
                                      * ((np.roll(t_in, -1) - t_in) / (np.roll(t_in, -2) - np.roll(t_in, -1)))
                                      - np.log((~R_in) * np.roll(R_in, -1))) * -0.25; 248, identifier:R_in; 249, unary_operator; 250, unary_operator; 251, integer:1; 252, integer:1; 253, identifier:R_in; 254, unary_operator; 255, unary_operator; 256, integer:1; 257, identifier:R_ip1; 258, identifier:i_in_for_out; 259, integer:1; 260, integer:1; 261, subscript; 262, subscript; 263, identifier:t_out; 264, subscript; 265, parenthesized_expression; 266, identifier:i_in_for_out; 267, identifier:R_in; 268, identifier:i_in_for_out; 269, identifier:A; 270, identifier:i_in_for_out; 271, identifier:B; 272, identifier:i_in_for_out; 273, identifier:dtype; 274, attribute; 275, identifier:side; 276, string; 277, parenthesized_expression; 278, float:0.25; 279, integer:1; 280, parenthesized_expression; 281, unary_operator; 282, integer:1; 283, subscript; 284, integer:1; 285, subscript; 286, identifier:t_in; 287, unary_operator; 288, identifier:t_in; 289, unary_operator; 290, identifier:t_in; 291, identifier:i_in_for_out; 292, binary_operator:t_inp1 - t_in; 293, identifier:np; 294, identifier:quaternion; 295, string_content:right; 296, binary_operator:- np.log((~R_in) * np.roll(R_in, -1))
                         + np.log((~np.roll(R_in, 1)) * R_in) * ((np.roll(t_in, -1) - t_in) / (t_in - np.roll(t_in, 1))); 297, binary_operator:np.log((~np.roll(R_in, -1)) * np.roll(R_in, -2))
                                      * ((np.roll(t_in, -1) - t_in) / (np.roll(t_in, -2) - np.roll(t_in, -1)))
                                      - np.log((~R_in) * np.roll(R_in, -1)); 298, float:0.25; 299, identifier:R_in; 300, unary_operator; 301, identifier:R_in; 302, unary_operator; 303, integer:1; 304, integer:2; 305, identifier:t_inp1; 306, identifier:t_in; 307, unary_operator; 308, binary_operator:np.log((~np.roll(R_in, 1)) * R_in) * ((np.roll(t_in, -1) - t_in) / (t_in - np.roll(t_in, 1))); 309, binary_operator:np.log((~np.roll(R_in, -1)) * np.roll(R_in, -2))
                                      * ((np.roll(t_in, -1) - t_in) / (np.roll(t_in, -2) - np.roll(t_in, -1))); 310, call; 311, integer:2; 312, integer:2; 313, call; 314, call; 315, parenthesized_expression; 316, call; 317, parenthesized_expression; 318, attribute; 319, argument_list; 320, attribute; 321, argument_list; 322, attribute; 323, argument_list; 324, binary_operator:(np.roll(t_in, -1) - t_in) / (t_in - np.roll(t_in, 1)); 325, attribute; 326, argument_list; 327, binary_operator:(np.roll(t_in, -1) - t_in) / (np.roll(t_in, -2) - np.roll(t_in, -1)); 328, identifier:np; 329, identifier:log; 330, binary_operator:(~R_in) * np.roll(R_in, -1); 331, identifier:np; 332, identifier:log; 333, binary_operator:(~R_in) * np.roll(R_in, -1); 334, identifier:np; 335, identifier:log; 336, binary_operator:(~np.roll(R_in, 1)) * R_in; 337, parenthesized_expression; 338, parenthesized_expression; 339, identifier:np; 340, identifier:log; 341, binary_operator:(~np.roll(R_in, -1)) * np.roll(R_in, -2); 342, parenthesized_expression; 343, parenthesized_expression; 344, parenthesized_expression; 345, call; 346, parenthesized_expression; 347, call; 348, parenthesized_expression; 349, identifier:R_in; 350, binary_operator:np.roll(t_in, -1) - t_in; 351, binary_operator:t_in - np.roll(t_in, 1); 352, parenthesized_expression; 353, call; 354, binary_operator:np.roll(t_in, -1) - t_in; 355, binary_operator:np.roll(t_in, -2) - np.roll(t_in, -1); 356, unary_operator; 357, attribute; 358, argument_list; 359, unary_operator; 360, attribute; 361, argument_list; 362, unary_operator; 363, call; 364, identifier:t_in; 365, identifier:t_in; 366, call; 367, unary_operator; 368, attribute; 369, argument_list; 370, call; 371, identifier:t_in; 372, call; 373, call; 374, identifier:R_in; 375, identifier:np; 376, identifier:roll; 377, identifier:R_in; 378, unary_operator; 379, identifier:R_in; 380, identifier:np; 381, identifier:roll; 382, identifier:R_in; 383, unary_operator; 384, call; 385, attribute; 386, argument_list; 387, attribute; 388, argument_list; 389, call; 390, identifier:np; 391, identifier:roll; 392, identifier:R_in; 393, unary_operator; 394, attribute; 395, argument_list; 396, attribute; 397, argument_list; 398, attribute; 399, argument_list; 400, integer:1; 401, integer:1; 402, attribute; 403, argument_list; 404, identifier:np; 405, identifier:roll; 406, identifier:t_in; 407, unary_operator; 408, identifier:np; 409, identifier:roll; 410, identifier:t_in; 411, integer:1; 412, attribute; 413, argument_list; 414, integer:2; 415, identifier:np; 416, identifier:roll; 417, identifier:t_in; 418, unary_operator; 419, identifier:np; 420, identifier:roll; 421, identifier:t_in; 422, unary_operator; 423, identifier:np; 424, identifier:roll; 425, identifier:t_in; 426, unary_operator; 427, identifier:np; 428, identifier:roll; 429, identifier:R_in; 430, integer:1; 431, integer:1; 432, identifier:np; 433, identifier:roll; 434, identifier:R_in; 435, unary_operator; 436, integer:1; 437, integer:2; 438, integer:1; 439, integer:1
code: def squad(R_in, t_in, t_out):
      """Spherical "quadrangular" interpolation of rotors with a cubic spline
  
      This is the best way to interpolate rotations.  It uses the analog
      of a cubic spline, except that the interpolant is confined to the
      rotor manifold in a natural way.  Alternative methods involving
      interpolation of other coordinates on the rotation group or
      normalization of interpolated values give bad results.  The
      results from this method are as natural as any, and are continuous
      in first and second derivatives.
  
      The input `R_in` rotors are assumed to be reasonably continuous
      (no sign flips), and the input `t` arrays are assumed to be
      sorted.  No checking is done for either case, and you may get
      silently bad results if these conditions are violated.
  
      This function simplifies the calling, compared to `squad_evaluate`
      (which takes a set of four quaternions forming the edges of the
      "quadrangle", and the normalized time `tau`) and `squad_vectorized`
      (which takes the same arguments, but in array form, and efficiently
      loops over them).
  
      Parameters
      ----------
      R_in: array of quaternions
          A time-series of rotors (unit quaternions) to be interpolated
      t_in: array of float
          The times corresponding to R_in
      t_out: array of float
          The times to which R_in should be interpolated
  
      """
      if R_in.size == 0 or t_out.size == 0:
          return np.array((), dtype=np.quaternion)
  
      # This list contains an index for each `t_out` such that
      # t_in[i-1] <= t_out < t_in[i]
      # Note that `side='right'` is much faster in my tests
      # i_in_for_out = t_in.searchsorted(t_out, side='left')
      # np.clip(i_in_for_out, 0, len(t_in) - 1, out=i_in_for_out)
      i_in_for_out = t_in.searchsorted(t_out, side='right')-1
  
      # Now, for each index `i` in `i_in`, we need to compute the
      # interpolation "coefficients" (`A_i`, `B_ip1`).
      #
      # I previously tested an explicit version of the loops below,
      # comparing `stride_tricks.as_strided` with explicit
      # implementation via `roll` (as seen here).  I found that the
      # `roll` was significantly more efficient for simple calculations,
      # though the difference is probably totally washed out here.  In
      # any case, it might be useful to test again.
      #
      A = R_in * np.exp((- np.log((~R_in) * np.roll(R_in, -1))
                         + np.log((~np.roll(R_in, 1)) * R_in) * ((np.roll(t_in, -1) - t_in) / (t_in - np.roll(t_in, 1)))
                         ) * 0.25)
      B = np.roll(R_in, -1) * np.exp((np.log((~np.roll(R_in, -1)) * np.roll(R_in, -2))
                                      * ((np.roll(t_in, -1) - t_in) / (np.roll(t_in, -2) - np.roll(t_in, -1)))
                                      - np.log((~R_in) * np.roll(R_in, -1))) * -0.25)
  
      # Correct the first and last A time steps, and last two B time steps.  We extend R_in with the following wrap-around
      # values:
      # R_in[0-1] = R_in[0]*(~R_in[1])*R_in[0]
      # R_in[n+0] = R_in[-1] * (~R_in[-2]) * R_in[-1]
      # R_in[n+1] = R_in[0] * (~R_in[-1]) * R_in[0]
      #           = R_in[-1] * (~R_in[-2]) * R_in[-1] * (~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1]
      #           = R_in[-1] * (~R_in[-2]) * R_in[-1] * (~R_in[-2]) * R_in[-1]
      # A[i] = R_in[i] * np.exp((- np.log((~R_in[i]) * R_in[i+1])
      #                          + np.log((~R_in[i-1]) * R_in[i]) * ((t_in[i+1] - t_in[i]) / (t_in[i] - t_in[i-1]))
      #                          ) * 0.25)
      # A[0] = R_in[0] * np.exp((- np.log((~R_in[0]) * R_in[1]) + np.log((~R_in[0])*R_in[1]*(~R_in[0])) * R_in[0]) * 0.25)
      #      = R_in[0]
      A[0] = R_in[0]
      # A[-1] = R_in[-1] * np.exp((- np.log((~R_in[-1]) * R_in[n+0])
      #                          + np.log((~R_in[-2]) * R_in[-1]) * ((t_in[n+0] - t_in[-1]) / (t_in[-1] - t_in[-2]))
      #                          ) * 0.25)
      #       = R_in[-1] * np.exp((- np.log((~R_in[-1]) * R_in[n+0]) + np.log((~R_in[-2]) * R_in[-1])) * 0.25)
      #       = R_in[-1] * np.exp((- np.log((~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1])
      #                           + np.log((~R_in[-2]) * R_in[-1])) * 0.25)
      #       = R_in[-1] * np.exp((- np.log((~R_in[-2]) * R_in[-1]) + np.log((~R_in[-2]) * R_in[-1])) * 0.25)
      #       = R_in[-1]
      A[-1] = R_in[-1]
      # B[i] = R_in[i+1] * np.exp((np.log((~R_in[i+1]) * R_in[i+2]) * ((t_in[i+1] - t_in[i]) / (t_in[i+2] - t_in[i+1]))
      #                            - np.log((~R_in[i]) * R_in[i+1])) * -0.25)
      # B[-2] = R_in[-1] * np.exp((np.log((~R_in[-1]) * R_in[0]) * ((t_in[-1] - t_in[-2]) / (t_in[0] - t_in[-1]))
      #                            - np.log((~R_in[-2]) * R_in[-1])) * -0.25)
      #       = R_in[-1] * np.exp((np.log((~R_in[-1]) * R_in[0]) - np.log((~R_in[-2]) * R_in[-1])) * -0.25)
      #       = R_in[-1] * np.exp((np.log((~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1])
      #                            - np.log((~R_in[-2]) * R_in[-1])) * -0.25)
      #       = R_in[-1] * np.exp((np.log((~R_in[-2]) * R_in[-1]) - np.log((~R_in[-2]) * R_in[-1])) * -0.25)
      #       = R_in[-1]
      B[-2] = R_in[-1]
      # B[-1] = R_in[0]
      # B[-1] = R_in[0] * np.exp((np.log((~R_in[0]) * R_in[1]) - np.log((~R_in[-1]) * R_in[0])) * -0.25)
      #       = R_in[-1] * (~R_in[-2]) * R_in[-1]
      #         * np.exp((np.log((~(R_in[-1] * (~R_in[-2]) * R_in[-1])) * R_in[-1] * (~R_in[-2]) * R_in[-1] * (~R_in[-2]) * R_in[-1])
      #                  - np.log((~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1])) * -0.25)
      #       = R_in[-1] * (~R_in[-2]) * R_in[-1]
      #         * np.exp((np.log(((~R_in[-1]) * R_in[-2] * (~R_in[-1])) * R_in[-1] * (~R_in[-2]) * R_in[-1] * (~R_in[-2]) * R_in[-1])
      #                  - np.log((~R_in[-1]) * R_in[-1] * (~R_in[-2]) * R_in[-1])) * -0.25)
      #         * np.exp((np.log((~R_in[-2]) * R_in[-1])
      #                  - np.log((~R_in[-2]) * R_in[-1])) * -0.25)
      B[-1] = R_in[-1] * (~R_in[-2]) * R_in[-1]
  
      # Use the coefficients at the corresponding t_out indices to
      # compute the squad interpolant
      # R_ip1 = np.array(np.roll(R_in, -1)[i_in_for_out])
      # R_ip1[-1] = R_in[-1]*(~R_in[-2])*R_in[-1]
      R_ip1 = np.roll(R_in, -1)
      R_ip1[-1] = R_in[-1]*(~R_in[-2])*R_in[-1]
      R_ip1 = np.array(R_ip1[i_in_for_out])
      t_inp1 = np.roll(t_in, -1)
      t_inp1[-1] = t_in[-1] + (t_in[-1] - t_in[-2])
      tau = (t_out - t_in[i_in_for_out]) / ((t_inp1 - t_in)[i_in_for_out])
      # tau = (t_out - t_in[i_in_for_out]) / ((np.roll(t_in, -1) - t_in)[i_in_for_out])
      R_out = np.squad_vectorized(tau, R_in[i_in_for_out], A[i_in_for_out], B[i_in_for_out], R_ip1)
  
      return R_out