Doesnt work for inverted arrays

[sharp-trickster]

I've tried this function in C++ and tested 2 arrays a, and b, being:

double a[10] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
double b[10] = { 10, 9, 8, 7, 6, 5, 4, 3, 2, 1 };

Say, for N from -2 to 11

When I run interp(N, a, b) it works:

i=-2 t=10 i=-1 t=10 i=0 t=10 i=1 t=10 i=2 t=9 i=3 t=8 i=4 t=7 i=5 t=6 i=6 t=5 i=7 t=4 i=8 t=3 i=9 t=2 i=10 t=1 i=11 t=1

When I run interp(N, b, a) it doesnt work:

i=-2 t=1 i=-1 t=1 i=0 t=1 i=1 t=1 i=2 t=1 i=3 t=1 i=4 t=1 i=5 t=1 i=6 t=1 i=7 t=1 i=8 t=1 i=9 t=10 i=10 t=10 i=11 t=10

[sharp-trickster]

I made it work changing your function slightly for my inverted cases:

double* invInterp(int N_x, int N_xp, const double* x, const double* xp, const double* yp, double left, double right) {

	double* y = new double[N_x];
	for (int i = 0; i < N_x; i++) y[i] = left;

	int ip = 0;
	int ip_next = 1;
	int i = 0;

	while (i < N_x) {

		double m = ((yp[ip]) - yp[ip_next]) / (xp[ip_next] - xp[ip]);
		
		double q = yp[ip] + m * xp[ip];
		while (x[i] > xp[ip_next]) {
			if (x[i] <= xp[ip])
				y[i] = q - m * x[i];
			i += 1;
			if (i >= N_x) { break; }
			//std::cout << "i:" << i << " ip:" <<ip << "\n";
		}
		ip += 1;
		ip_next += 1;
		if (ip_next == N_xp) {
			while (i < N_x) {
				y[i] = right;
				i++;
			}
			break;
		}
	}

	return y;
}

[stefanoschmidt1995]

Hi! Thanks for spotting this and reaching out.
Indeed, this is a known limitation and I'm glad you've worked to overcome this.

I guess this can be fixed in two ways:

• on the C layer, as you did
• on the python layer, by simply inverting the order of the vector before calling the C code.
  Not sure which one is faster, though.

If you're confident with your solution, feel free to submit a pull request and we can add your function to the public interface.

Cheers!