interpolation stuff

[spinkney]

• suggested in https://github.com/stan-dev/design-docs/pull/18#issuecomment-628717246 coded up in https://discourse.mc-stan.org/t/linear-interpolation-and-searchsorted-in-stan/13318/6?u=spinkney

real linear_interpolation_v(real x, vector x_pred, vector y_pred){
    int K = rows(x_pred);
    vector[K] deltas = x - x_pred;
    real ans;
    real t;
    real w;
    int i;
    if(x < x_pred[1] || x > x_pred[K]) reject("x is outside of the x_pred grid!");
    if(rows(y_pred) != K) reject("x_pred and y_pred aren't of the same size");
    //this is which.min()
    i = sort_indices_asc(fabs(deltas))[1];
    if(deltas[i] <= 0) i -= 1;
    ans = y_pred[i];
    real x1 = x_pred[i];
    real x2 = x_pred[i + 1];
    real y1 = y_pred[i];
    real y2 = y_pred[i + 1];
    ans = y1 + (y2 - y1) * (x - x1) / (x2 - x1);
    t = (x - x1) / (x2 - x1);
   // 3 weight functions here to handle
    w = 1 - t;
    //w = 1 / (1 + exp(1 / (1 - t) - 1 / t));
    //w = 1 - 3 * pow(t, 2) + 2 * pow(t, 3);
    ans = w * y1 + (1 - w) * y2;
    return ans;
  }

• Quadratic envelop method in https://github.com/stan-dev/design-docs/pull/18#issuecomment-812690206. This method leaves out the bookend bins. Maybe do linear interpolation there?

vector get_quad_coeffs(vector x3, vector y3){
  matrix[3, 3] xm = rep_matrix(1.0, 3, 3) ;
  for (i in 2:3)
    xm[, i] = xm[, i - 1] .* x3;
  return xm \ y3;
  }

matrix grid_to_quad(vector xs, vector ys){
  int N = rows(ys);
  int M = N - 1;
  matrix[N, 3] grdm = rep_matrix(not_a_number(), N, 3);
  vector[3] as;
  int idx[3];
  if (N<3) reject("data grid should have at least 3 rows. Found N=",N);
  for (i in 2:M){
    idx = {i - 1, i, i + 1};
    as = get_quad_coeffs(xs[idx], ys[idx]);
    grdm[i, ] = to_row_vector(as);
  }
  return grdm;
}

real get_w(real t){
  return 1 / (1 + exp(1 / (1 - t) - 1 / t));
 // return 1 - t;
 // return 1 - 3 * pow(t, 2) + 2 * pow(t, 3);
}


real get_smooth(real x, vector xs, vector ys){
  int N =rows(xs);
  real t;
  real w;
  real res;
  matrix[N, 1] v;
  matrix[N, 3] grdm = grid_to_quad(xs, ys);
  matrix[3, 1] xp = [[1], [x], [square(x)]];
  for (i in 1:N){
    if (x <= xs[i] && x >= xs[i - 1]) {
      t = (x - xs[i - 1]) / (xs[i] - xs[i-1]);
      w = get_w(t);
      v = grdm * xp;
      res = w * v[i - 1, 1] + (1 - w) * v[i, 1];
    }
  }
  return res;
}

[spinkney]

monotone cubic interpolation https://en.wikipedia.org/wiki/Monotone_cubic_interpolation