Millisecond Forums

+x

 

Andrew Papale - 11/6/2023

+x

 

Dave - 11/6/2023

+x

 

Andrew Papale - 11/6/2023

+x

 

Andrew Papale - 11/6/2023

+x

 

Andrew Papale - 11/1/2023

+x

 

Dave - 10/30/2023

+x

 

Andrew Papale - 10/28/2023

+x

 

Dave - 10/27/2023

+x

 

Andrew Papale - 10/26/2023

Dear Inquisit Community,

Has anyone tried to implement a cubic spline algorithm through ~30 points in Inquisit?

The problem: I have a list of values (e.g. 1-100) for each trial in an Inquisit list that roughly forms a continuous value landscape with peaks and valleys across 360 possible choices.  At some trial, I insert a new value (e.g. 25) into one of the 360 possible choices (say it was 100).  Therefore, the value at the point that used to be 100 is now 25.  Then, I want to fill in (as smoothly as possible) the values between the inserted value of 25 and the old value landscape that includes points adjacent to the old value (e.g. 90, 80, 70 on the left; and 110, 120, 130 on the right).

Is anything like this even possible? I currently have a parabolic approximation through an endpoint and a vertex working fairly well, but thought I'd ask if there is a better solution.

Thanks,
Andrew

// input: parameter space to construct the splines over
<list x>
/ items = (-1.5, -0.2, 1.0, 5.0, 10.0, 15.0, 20.0)
</list>
<list y>
/ items = (-1.2, 0.0, 0.5, 1.0, 1.2, 2.0, 1.0)
</list>

// temp storage
<list u>
</list>

// output coefficients
<list y2>
</list>

<expressions>
// constructs the cubic splines over the parameter space
/ createSpline = {
    list.y2.reset();
    list.u.reset();
    var n = list.x.itemcount;
    // setting up the vectors with n items
    var i = 1;
    while (i <= n) {
        list.y2.appenditem(null);
        list.u.appenditem(null);
        i += 1;
    };
    // set natural boundary condition
    list.y2.setitem(0.0, 1);
    list.u.setitem(0.0, 1);
    // decomposition
    i = 2;
    while (i <= n-1) {
        //sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
        var sig = (list.x.item(i) - list.x.item(i-1)) / (list.x.item(i+1) - list.x.item(i-1));
        //p=sig*y2[i-1]+2.0
        var p = sig * list.y2.item(i-1) + 2.0;
        //y2[i]=(sig-1.0)/p;
        list.y2.setitem((sig - 1.0) / p, i);
        //u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
        var ui = (list.y.item(i+1) - list.y.item(i)) / (list.x.item(i+1) - list.x.item(i)) - (list.y.item(i) - list.y.item(i-1)) / (list.x.item(i) - list.x.item(i-1));
        list.u.setitem(ui, i);
        //u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
        ui = (6.0 * list.u.item(i) / (list.x.item(i+1) - list.x.item(i-1)) - sig * list.u.item(i-1)) / p;
        list.u.setitem(ui, i);
        i += 1;
    };
    // natural boundary condition
    var qn = 0.0;
    var un = 0.0;
    //y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
    list.y2.setitem((un - qn * list.u.item(n-1)) / (qn * list.y2.item(n-1) + 1.0), n);
    // backsubstitution
    var k = n-1;
    while (k >= 1) {
        //y2[k]=y2[k]*y2[k+1]+u[k];
        var y2k = list.y2.item(k) * list.y2.item(k+1) + list.u.item(k);
        list.y2.setitem(y2k, k);
        k -= 1;
    };
};

// sets up a test vector of x values in the range x1 to xn
/ createTestVector = {
    var start_x = list.x.item(1);
    var end_x = list.x.item(list.x.itemcount);
    var x = start_x;
    // set up a test vector of x values in range x1 to xn
    while (x <= end_x) {
        list.x_test.appenditem(x);
        list.y_test.appenditem(null);
        x += 0.25;
    };
};

/ interpolateSpline = {
    // bisect
    var n = list.x.itemcount;
    var klo = 1;
    var khi = n;
    while (khi - klo > 1) {
        var k = floor((khi + klo) / 2); // same as right bit shift
        if (list.x.item(k) > values.x) {
            khi = k;
        } else {
            klo = k;
        };
    };
    var h = list.x.item(khi) - list.x.item(klo);
    //a=(xa[khi]-x)/h;
    var a = (list.x.item(khi) - values.x) / h;
    //b=(x-xa[klo])/h;
    var b = (values.x - list.x.item(klo)) / h;
    // calculate interpolated y
    // y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0
    values.y = a * list.y.item(klo) + b * list.y.item(khi) + ((a*a*a - a) * list.y2.item(klo) + (b*b*b -b) * list.y2.item(khi)) * (h*h) / 6.0;
    return values.y;
};
</expressions>

<list x_test>
/ selectionmode = sequence
</list>
<list y_test>
</list>

<values>
/ x = null
/ y = null
</values>

<block exampleblock>
/ onblockbegin = [
    expressions.createSpline;
    expressions.createTestVector;
]
/ trials = [1=interpolateSpline]
</block>

<trial interpolateSpline>
/ ontrialbegin = [
    values.x = list.x_test.nextvalue; // get our x value
    values.y = null;
    expressions.interpolateSpline; // interpolate y
    list.y_test.setitem(values.y, list.x_test.currentindex);
]
/ trialduration = 0
/ branch = [
    if (list.x_test.selectedcount < list.x_test.itemcount) {
        return trial.interpolateSpline;
    };
]
</trial>

<data>
/ columns = (trialnum values.x values.y)
/ separatefiles = true
</data>

The above is essentially a straight Inquisit implementation of the cubic spline interpolation algorithm from Numerical Recipes in C.

Plotting the output, i.e. interpolated x / y values, gives



with





being the original points the splines were constructed over (shown with & without linear interpolation).

Or here's the interpolation result over a very different set of inputs

<list x>
/ items = (-10, -6.5, -3.5, 0, 3.5, 6.5, 10)
</list>
<list y>
/ items = (0, 0.1, -0.1, 1, -0.1, 0.1, 0)
</list>

Dave,

This is absolutely fantastic.  Thank you so much on behalf of the entire lab.

Andre

FYI, I cleaned up the code in the original response a little bit and added a few more test inputs and comments. (No bug fixes or major changes to the code itself, so it's fine to keep using the previous version.)

Thanks Dave, I've got this working basically how we want it to now (see attached).  We have these values mapped onto 360 degrees of a circle and I need to figure out 1) how to handle the cases where these insertions appear on the edges near 0/360 degrees.  Then, I need to 2) iteratively insert all these erasures (which will occur at some frequency set by the experimental design) into all future trials, which occur in multiples of 360 in the values-5888.txt and RTs-5888.txt vectors.  I have a general idea of how to proceed with 2) using while loops and modular division and 1) shifting the RT values, performing the spline interpolation, and then shifting them back into place.

Any advice is appreciated, though I feel like you've done enough for us here already.

Dear Dave,

I would like some algorithm advice.  I have successfully iterated the spline generation over all future trials of my task (yay for small victories).  However, the spline computation takes about 80ms per trial, which is prohibitively long when trying to do this over 100+ trials (~8s).  My thought was to instead compute the spline only once on trial N, save the values, and then paste them into trials N+1, N+2, ... ,N+M up to 300.  This will be slightly different than my original plan but I think it will save computation time.

The copying and pasting approach takes about 2.5s for 100+ trials.  This is still too slow.  The target is about 500ms.  My next thought was to copy and paste only 10 or so trials in the future, like a buffer.  This will be slightly more tricky to code I think.

I really don't know what you refer to as the "copying and pasting" approach, i.e. what exactly you implemented, why and how

I apologize for the confusion.  I have created a spline on trial 1 between points (for example) 135 - 165 (as in the picture above).  There is an underlying value distribution that is changing across trials 1-300 over top of which I need to insert this spline.  Now, I have copied that spline into trials 2 through 300 for the same x values.  I have encoded all values for all trials in one really long list <values-5888.txt> of size (360*300) It is taking too long to compute for a real time task and I was asking for advice on how to do this algorithm faster.  Sorry this is a fairly complicated task.

+x

 

Dave - 11/6/2023

+x

 

Andrew Papale - 11/6/2023

+x

 

Dave - 11/6/2023

+x

 

Andrew Papale - 11/6/2023

+x

 

Andrew Papale - 11/6/2023

+x

 

Andrew Papale - 11/1/2023

+x

 

Dave - 10/30/2023

+x

 

Andrew Papale - 10/28/2023

+x

 

Dave - 10/27/2023

+x

 

Andrew Papale - 10/26/2023

Dear Inquisit Community,

Has anyone tried to implement a cubic spline algorithm through ~30 points in Inquisit?

The problem: I have a list of values (e.g. 1-100) for each trial in an Inquisit list that roughly forms a continuous value landscape with peaks and valleys across 360 possible choices.  At some trial, I insert a new value (e.g. 25) into one of the 360 possible choices (say it was 100).  Therefore, the value at the point that used to be 100 is now 25.  Then, I want to fill in (as smoothly as possible) the values between the inserted value of 25 and the old value landscape that includes points adjacent to the old value (e.g. 90, 80, 70 on the left; and 110, 120, 130 on the right).

Is anything like this even possible? I currently have a parabolic approximation through an endpoint and a vertex working fairly well, but thought I'd ask if there is a better solution.

Thanks,
Andrew

// input: parameter space to construct the splines over
<list x>
/ items = (-1.5, -0.2, 1.0, 5.0, 10.0, 15.0, 20.0)
</list>
<list y>
/ items = (-1.2, 0.0, 0.5, 1.0, 1.2, 2.0, 1.0)
</list>

// temp storage
<list u>
</list>

// output coefficients
<list y2>
</list>

<expressions>
// constructs the cubic splines over the parameter space
/ createSpline = {
    list.y2.reset();
    list.u.reset();
    var n = list.x.itemcount;
    // setting up the vectors with n items
    var i = 1;
    while (i <= n) {
        list.y2.appenditem(null);
        list.u.appenditem(null);
        i += 1;
    };
    // set natural boundary condition
    list.y2.setitem(0.0, 1);
    list.u.setitem(0.0, 1);
    // decomposition
    i = 2;
    while (i <= n-1) {
        //sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
        var sig = (list.x.item(i) - list.x.item(i-1)) / (list.x.item(i+1) - list.x.item(i-1));
        //p=sig*y2[i-1]+2.0
        var p = sig * list.y2.item(i-1) + 2.0;
        //y2[i]=(sig-1.0)/p;
        list.y2.setitem((sig - 1.0) / p, i);
        //u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
        var ui = (list.y.item(i+1) - list.y.item(i)) / (list.x.item(i+1) - list.x.item(i)) - (list.y.item(i) - list.y.item(i-1)) / (list.x.item(i) - list.x.item(i-1));
        list.u.setitem(ui, i);
        //u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
        ui = (6.0 * list.u.item(i) / (list.x.item(i+1) - list.x.item(i-1)) - sig * list.u.item(i-1)) / p;
        list.u.setitem(ui, i);
        i += 1;
    };
    // natural boundary condition
    var qn = 0.0;
    var un = 0.0;
    //y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
    list.y2.setitem((un - qn * list.u.item(n-1)) / (qn * list.y2.item(n-1) + 1.0), n);
    // backsubstitution
    var k = n-1;
    while (k >= 1) {
        //y2[k]=y2[k]*y2[k+1]+u[k];
        var y2k = list.y2.item(k) * list.y2.item(k+1) + list.u.item(k);
        list.y2.setitem(y2k, k);
        k -= 1;
    };
};

// sets up a test vector of x values in the range x1 to xn
/ createTestVector = {
    var start_x = list.x.item(1);
    var end_x = list.x.item(list.x.itemcount);
    var x = start_x;
    // set up a test vector of x values in range x1 to xn
    while (x <= end_x) {
        list.x_test.appenditem(x);
        list.y_test.appenditem(null);
        x += 0.25;
    };
};

/ interpolateSpline = {
    // bisect
    var n = list.x.itemcount;
    var klo = 1;
    var khi = n;
    while (khi - klo > 1) {
        var k = floor((khi + klo) / 2); // same as right bit shift
        if (list.x.item(k) > values.x) {
            khi = k;
        } else {
            klo = k;
        };
    };
    var h = list.x.item(khi) - list.x.item(klo);
    //a=(xa[khi]-x)/h;
    var a = (list.x.item(khi) - values.x) / h;
    //b=(x-xa[klo])/h;
    var b = (values.x - list.x.item(klo)) / h;
    // calculate interpolated y
    // y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0
    values.y = a * list.y.item(klo) + b * list.y.item(khi) + ((a*a*a - a) * list.y2.item(klo) + (b*b*b -b) * list.y2.item(khi)) * (h*h) / 6.0;
    return values.y;
};
</expressions>

<list x_test>
/ selectionmode = sequence
</list>
<list y_test>
</list>

<values>
/ x = null
/ y = null
</values>

<block exampleblock>
/ onblockbegin = [
    expressions.createSpline;
    expressions.createTestVector;
]
/ trials = [1=interpolateSpline]
</block>

<trial interpolateSpline>
/ ontrialbegin = [
    values.x = list.x_test.nextvalue; // get our x value
    values.y = null;
    expressions.interpolateSpline; // interpolate y
    list.y_test.setitem(values.y, list.x_test.currentindex);
]
/ trialduration = 0
/ branch = [
    if (list.x_test.selectedcount < list.x_test.itemcount) {
        return trial.interpolateSpline;
    };
]
</trial>

<data>
/ columns = (trialnum values.x values.y)
/ separatefiles = true
</data>

The above is essentially a straight Inquisit implementation of the cubic spline interpolation algorithm from Numerical Recipes in C.

Plotting the output, i.e. interpolated x / y values, gives



with





being the original points the splines were constructed over (shown with & without linear interpolation).

Or here's the interpolation result over a very different set of inputs

<list x>
/ items = (-10, -6.5, -3.5, 0, 3.5, 6.5, 10)
</list>
<list y>
/ items = (0, 0.1, -0.1, 1, -0.1, 0.1, 0)
</list>

Dave,

This is absolutely fantastic.  Thank you so much on behalf of the entire lab.

Andre

FYI, I cleaned up the code in the original response a little bit and added a few more test inputs and comments. (No bug fixes or major changes to the code itself, so it's fine to keep using the previous version.)

Thanks Dave, I've got this working basically how we want it to now (see attached).  We have these values mapped onto 360 degrees of a circle and I need to figure out 1) how to handle the cases where these insertions appear on the edges near 0/360 degrees.  Then, I need to 2) iteratively insert all these erasures (which will occur at some frequency set by the experimental design) into all future trials, which occur in multiples of 360 in the values-5888.txt and RTs-5888.txt vectors.  I have a general idea of how to proceed with 2) using while loops and modular division and 1) shifting the RT values, performing the spline interpolation, and then shifting them back into place.

Any advice is appreciated, though I feel like you've done enough for us here already.

Dear Dave,

I would like some algorithm advice.  I have successfully iterated the spline generation over all future trials of my task (yay for small victories).  However, the spline computation takes about 80ms per trial, which is prohibitively long when trying to do this over 100+ trials (~8s).  My thought was to instead compute the spline only once on trial N, save the values, and then paste them into trials N+1, N+2, ... ,N+M up to 300.  This will be slightly different than my original plan but I think it will save computation time.

The copying and pasting approach takes about 2.5s for 100+ trials.  This is still too slow.  The target is about 500ms.  My next thought was to copy and paste only 10 or so trials in the future, like a buffer.  This will be slightly more tricky to code I think.

I really don't know what you refer to as the "copying and pasting" approach, i.e. what exactly you implemented, why and how

I apologize for the confusion.  I have created a spline on trial 1 between points (for example) 135 - 165 (as in the picture above).  There is an underlying value distribution that is changing across trials 1-300 over top of which I need to insert this spline.  Now, I have copied that spline into trials 2 through 300 for the same x values.  I have encoded all values for all trials in one really long list <values-5888.txt> of size (360*300) It is taking too long to compute for a real time task and I was asking for advice on how to do this algorithm faster.  Sorry this is a fairly complicated task.

I'm sure you've done some profiling and measurements, so *what* exact part or parts of your code are taking too long to compute? As for the algorithm in the original post, it's not going to get any faster. However, given that -- as I noted -- spline construction over 15 points takes about 40 ms, and obtaining an interpolated y value takes less than 1 ms, I doubt that the spline creation or interpolation routines are the bottleneck here.

+x

 

Dave - 11/6/2023

+x

 

Dave - 11/6/2023

+x

 

Andrew Papale - 11/6/2023

+x

 

Dave - 11/6/2023

+x

 

Andrew Papale - 11/6/2023

+x

 

Andrew Papale - 11/6/2023

+x

 

Andrew Papale - 11/1/2023

+x

 

Dave - 10/30/2023

+x

 

Andrew Papale - 10/28/2023

+x

 

Dave - 10/27/2023

+x

 

Andrew Papale - 10/26/2023

Dear Inquisit Community,

Has anyone tried to implement a cubic spline algorithm through ~30 points in Inquisit?

The problem: I have a list of values (e.g. 1-100) for each trial in an Inquisit list that roughly forms a continuous value landscape with peaks and valleys across 360 possible choices.  At some trial, I insert a new value (e.g. 25) into one of the 360 possible choices (say it was 100).  Therefore, the value at the point that used to be 100 is now 25.  Then, I want to fill in (as smoothly as possible) the values between the inserted value of 25 and the old value landscape that includes points adjacent to the old value (e.g. 90, 80, 70 on the left; and 110, 120, 130 on the right).

Is anything like this even possible? I currently have a parabolic approximation through an endpoint and a vertex working fairly well, but thought I'd ask if there is a better solution.

Thanks,
Andrew

// input: parameter space to construct the splines over
<list x>
/ items = (-1.5, -0.2, 1.0, 5.0, 10.0, 15.0, 20.0)
</list>
<list y>
/ items = (-1.2, 0.0, 0.5, 1.0, 1.2, 2.0, 1.0)
</list>

// temp storage
<list u>
</list>

// output coefficients
<list y2>
</list>

<expressions>
// constructs the cubic splines over the parameter space
/ createSpline = {
    list.y2.reset();
    list.u.reset();
    var n = list.x.itemcount;
    // setting up the vectors with n items
    var i = 1;
    while (i <= n) {
        list.y2.appenditem(null);
        list.u.appenditem(null);
        i += 1;
    };
    // set natural boundary condition
    list.y2.setitem(0.0, 1);
    list.u.setitem(0.0, 1);
    // decomposition
    i = 2;
    while (i <= n-1) {
        //sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
        var sig = (list.x.item(i) - list.x.item(i-1)) / (list.x.item(i+1) - list.x.item(i-1));
        //p=sig*y2[i-1]+2.0
        var p = sig * list.y2.item(i-1) + 2.0;
        //y2[i]=(sig-1.0)/p;
        list.y2.setitem((sig - 1.0) / p, i);
        //u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
        var ui = (list.y.item(i+1) - list.y.item(i)) / (list.x.item(i+1) - list.x.item(i)) - (list.y.item(i) - list.y.item(i-1)) / (list.x.item(i) - list.x.item(i-1));
        list.u.setitem(ui, i);
        //u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
        ui = (6.0 * list.u.item(i) / (list.x.item(i+1) - list.x.item(i-1)) - sig * list.u.item(i-1)) / p;
        list.u.setitem(ui, i);
        i += 1;
    };
    // natural boundary condition
    var qn = 0.0;
    var un = 0.0;
    //y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
    list.y2.setitem((un - qn * list.u.item(n-1)) / (qn * list.y2.item(n-1) + 1.0), n);
    // backsubstitution
    var k = n-1;
    while (k >= 1) {
        //y2[k]=y2[k]*y2[k+1]+u[k];
        var y2k = list.y2.item(k) * list.y2.item(k+1) + list.u.item(k);
        list.y2.setitem(y2k, k);
        k -= 1;
    };
};

// sets up a test vector of x values in the range x1 to xn
/ createTestVector = {
    var start_x = list.x.item(1);
    var end_x = list.x.item(list.x.itemcount);
    var x = start_x;
    // set up a test vector of x values in range x1 to xn
    while (x <= end_x) {
        list.x_test.appenditem(x);
        list.y_test.appenditem(null);
        x += 0.25;
    };
};

/ interpolateSpline = {
    // bisect
    var n = list.x.itemcount;
    var klo = 1;
    var khi = n;
    while (khi - klo > 1) {
        var k = floor((khi + klo) / 2); // same as right bit shift
        if (list.x.item(k) > values.x) {
            khi = k;
        } else {
            klo = k;
        };
    };
    var h = list.x.item(khi) - list.x.item(klo);
    //a=(xa[khi]-x)/h;
    var a = (list.x.item(khi) - values.x) / h;
    //b=(x-xa[klo])/h;
    var b = (values.x - list.x.item(klo)) / h;
    // calculate interpolated y
    // y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0
    values.y = a * list.y.item(klo) + b * list.y.item(khi) + ((a*a*a - a) * list.y2.item(klo) + (b*b*b -b) * list.y2.item(khi)) * (h*h) / 6.0;
    return values.y;
};
</expressions>

<list x_test>
/ selectionmode = sequence
</list>
<list y_test>
</list>

<values>
/ x = null
/ y = null
</values>

<block exampleblock>
/ onblockbegin = [
    expressions.createSpline;
    expressions.createTestVector;
]
/ trials = [1=interpolateSpline]
</block>

<trial interpolateSpline>
/ ontrialbegin = [
    values.x = list.x_test.nextvalue; // get our x value
    values.y = null;
    expressions.interpolateSpline; // interpolate y
    list.y_test.setitem(values.y, list.x_test.currentindex);
]
/ trialduration = 0
/ branch = [
    if (list.x_test.selectedcount < list.x_test.itemcount) {
        return trial.interpolateSpline;
    };
]
</trial>

<data>
/ columns = (trialnum values.x values.y)
/ separatefiles = true
</data>

The above is essentially a straight Inquisit implementation of the cubic spline interpolation algorithm from Numerical Recipes in C.

Plotting the output, i.e. interpolated x / y values, gives



with





being the original points the splines were constructed over (shown with & without linear interpolation).

Or here's the interpolation result over a very different set of inputs

<list x>
/ items = (-10, -6.5, -3.5, 0, 3.5, 6.5, 10)
</list>
<list y>
/ items = (0, 0.1, -0.1, 1, -0.1, 0.1, 0)
</list>

Dave,

This is absolutely fantastic.  Thank you so much on behalf of the entire lab.

Andre

FYI, I cleaned up the code in the original response a little bit and added a few more test inputs and comments. (No bug fixes or major changes to the code itself, so it's fine to keep using the previous version.)

Thanks Dave, I've got this working basically how we want it to now (see attached).  We have these values mapped onto 360 degrees of a circle and I need to figure out 1) how to handle the cases where these insertions appear on the edges near 0/360 degrees.  Then, I need to 2) iteratively insert all these erasures (which will occur at some frequency set by the experimental design) into all future trials, which occur in multiples of 360 in the values-5888.txt and RTs-5888.txt vectors.  I have a general idea of how to proceed with 2) using while loops and modular division and 1) shifting the RT values, performing the spline interpolation, and then shifting them back into place.

Any advice is appreciated, though I feel like you've done enough for us here already.

Dear Dave,

I would like some algorithm advice.  I have successfully iterated the spline generation over all future trials of my task (yay for small victories).  However, the spline computation takes about 80ms per trial, which is prohibitively long when trying to do this over 100+ trials (~8s).  My thought was to instead compute the spline only once on trial N, save the values, and then paste them into trials N+1, N+2, ... ,N+M up to 300.  This will be slightly different than my original plan but I think it will save computation time.

The copying and pasting approach takes about 2.5s for 100+ trials.  This is still too slow.  The target is about 500ms.  My next thought was to copy and paste only 10 or so trials in the future, like a buffer.  This will be slightly more tricky to code I think.

I really don't know what you refer to as the "copying and pasting" approach, i.e. what exactly you implemented, why and how

I apologize for the confusion.  I have created a spline on trial 1 between points (for example) 135 - 165 (as in the picture above).  There is an underlying value distribution that is changing across trials 1-300 over top of which I need to insert this spline.  Now, I have copied that spline into trials 2 through 300 for the same x values.  I have encoded all values for all trials in one really long list <values-5888.txt> of size (360*300) It is taking too long to compute for a real time task and I was asking for advice on how to do this algorithm faster.  Sorry this is a fairly complicated task.

I'm sure you've done some profiling and measurements, so *what* exact part or parts of your code are taking too long to compute? As for the algorithm in the original post, it's not going to get any faster. However, given that -- as I noted -- spline construction over 15 points takes about 40 ms, and obtaining an interpolated y value takes less than 1 ms, I doubt that the spline creation or interpolation routines are the bottleneck here.

Here's the original code with some runtime measurements plugged in.

<usermanual>
    For background, see:
    Press et al., Numerical Recipes in C: The Art of Scientific Computing, 2nd Ed. (1992),
    Chapter 3.3 Cubic Spline Interpolation, pg. 113 - 116.
</usermanual>

// input: parameter space to construct the splines over
// x values must be in ascending order, x values must be unique
// at least 3 points are required, otherwise interpolation
// obviously isn't possible
<list x>
/ items = (
    -10, -6.5, -3.5, -1.5, -0.2, 0, 1, 3.5, 5, 6.5, 10, 15, 20, 23, 30,
    32, 36.5, 39.5, 40.5, 41.2, 42, 43, 43.5, 45, 46.5, 50, 55, 60, 63, 70)
// items = (0,1,2)
// items = (-10, -6.5, -3.5)
// items = (-10, -6.5, -3.5, -1.5, -0.2, 0, 1, 3.5, 5, 6.5, 10, 15, 20, 23, 30) // x values of 1st example input space
// items = (-10, -6.5, -3.5, 0, 3.5, 6.5, 10) // x values of 2nd example input space
// items = (-1.5, -0.2, 1.0, 5.0, 10.0, 15.0, 20.0) // x values of 3rd example input space
/ selectionmode = sequence
</list>
<list y>
/ items = (
    0, 0.1, -0.1, -1.2, 0, -0.4, 0.5, -0.1, 0.1, 0, 2, 1, 5, 3, -1,
    1, 0.3, -0.2, -1.4, -0.9, -0.7, 0.25, -0.41, 0.3, 0, 4, 0, 3, 3, -1)
// items = (1,3,2)
// items = (0, 0.1, -0.1)
// items = (0, 0.1, -0.1, -1.2, 0, -0.4, 0.5, -0.1, 0.1, 0, 2, 1, 5, 3, -1) // y values of 1st example input space
// items = (0, 0.1, -0.1, 1, -0.1, 0.1, 0) // y values of 2nd example input space
// items = (-1.2, 0.0, 0.5, 1.0, 1.2, 2.0, 1.0) // y values of 3rd example input space
/ selectionmode = sequence
</list>

// temp storage
<list u>
</list>
// output y of 2nd derivative
<list y2>
</list>

<expressions>
// constructs the cubic splines over the input parameter space
/ createSpline = {
    values.tStartConstructionRoutine = script.elapsedtime;
    list.y2.reset();
    list.u.reset();
    var n = list.x.itemcount;
    // setting up the vectors with n items
    var i = 1;
    while (i <= n) {
        list.y2.appenditem(null);
        list.u.appenditem(null);
        i += 1;
    };
    i = 1;
    // natural boundary condition on leftmost point
    list.y2.setitem(0.0, i);
    list.u.setitem(0.0, i);
    // decomposition
    i = 2;
    while (i <= n-1) {
        //sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
        var sig = (list.x.item(i) - list.x.item(i-1)) / (list.x.item(i+1) - list.x.item(i-1));
        //p=sig*y2[i-1]+2.0
        var p = sig * list.y2.item(i-1) + 2.0;
        //y2[i]=(sig-1.0)/p;
        var y2i = (sig - 1.0) / p;
        list.y2.setitem(y2i, i);
        //u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
        var ui = (list.y.item(i+1) - list.y.item(i)) / (list.x.item(i+1) - list.x.item(i)) - (list.y.item(i) - list.y.item(i-1)) / (list.x.item(i) - list.x.item(i-1));
        list.u.setitem(ui, i);
        //u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
        ui = (6.0 * list.u.item(i) / (list.x.item(i+1) - list.x.item(i-1)) - sig * list.u.item(i-1)) / p;
        list.u.setitem(ui, i);
        i += 1;
    };
    // natural boundary condition on rightmost point
    var qn = 0.0;
    var un = 0.0;
    //y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
    var y2n = (un - qn * list.u.item(n-1)) / (qn * list.y2.item(n-1) + 1.0);
    list.y2.setitem(y2n, n);
    // backsubstitution
    var k = n-1;
    while (k >= 1) {
        //y2[k]=y2[k]*y2[k+1]+u[k];
        var y2k = list.y2.item(k) * list.y2.item(k+1) + list.u.item(k);
        list.y2.setitem(y2k, k);
        k -= 1;
    };
    values.tEndConstructionRoutine = script.elapsedtime;
};

// sets up a test vector of x values in the range x1 to xn
/ createTestVector = {
    values.tStartTestVectorSetup = script.elapsedtime;
    var x1 = list.x.item(1);
    var xn = list.x.item(list.x.itemcount);
    var x = x1;
    // set up a test vector of x values in range x1 to xn
    while (x <= xn) {
        list.x_test.appenditem(x);
        list.y_test.appenditem(null);
        x += 0.03125; // decrease for denser sampling (more points), increase for sparser sampling (fewer points)
    };
    values.tEndTestVectorSetup = script.elapsedtime;
};

// given an x value, finds the insertion point and returns interpolated y value
/ interpolateSpline = {
    values.tStartInterpolationRoutine = script.elapsedtime;
    // bisect
    var x = values.x;
    var n = list.x.itemcount;
    var klo = 1;
    var khi = n;
    while (khi - klo > 1) {
        var k = floor((khi + klo) / 2); // same as right bit shift
        if (list.x.item(k) > x) {
            khi = k;
        } else {
            klo = k;
        };
    };
    var h = list.x.item(khi) - list.x.item(klo);
    //a=(xa[khi]-x)/h;
    var a = (list.x.item(khi) - x) / h;
    //b=(x-xa[klo])/h;
    var b = (x - list.x.item(klo)) / h;
    // calculate interpolated y
    // y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0
    var y = a * list.y.item(klo) + b * list.y.item(khi) + ((a*a*a - a) * list.y2.item(klo) + (b*b*b - b) * list.y2.item(khi)) * (h*h) / 6.0;
    values.tEndInterpolationRoutine = script.elapsedtime;
    return y;
};
</expressions>

// test vectors holding interpolated x / y values
<list x_test>
/ selectionmode = sequence
/ selectionrate = always
</list>
<list y_test>
</list>

<values>
/ x = null
/ y = null
</values>

<values>
/ tStartTestVectorSetup = null
/ tEndTestVectorSetup = null
/ tStartConstructionRoutine = null
/ tEndConstructionRoutine = null
/ tStartInterpolationRoutine = null
/ tEndInterpolationRoutine = null
/ tStartValueGeneration = null
/ tEndValueGeneration = null
</values>

<expressions>
/ runtimeTestVectorSetup = values.tEndTestVectorSetup - values.tStartTestVectorSetup;
/ runtimeConstructionRoutine = values.tEndConstructionRoutine - values.tStartConstructionRoutine;
/ runtimeInterpolationRoutine = values.tEndInterpolationRoutine - values.tStartInterpolationRoutine;
/ runtimeValueGeneration = values.tEndValueGeneration - values.tStartValueGeneration;
</expressions>

<block exampleBlock>
/ onblockbegin = [
    expressions.createSpline; // call the spline construction routine
    expressions.createTestVector; // call the test vector generation routine
]
/ trials = [1=interpolateSplineLoop;] // to generate the interpolated values in a single loop
// trials = [1=interpolateSpline;] // to write out the interpolated values one per row
</block>

// pass x values to the interpolation routine, get interpolated y values
<trial interpolateSpline>
/ ontrialbegin = [
    values.x = list.x_test.nextvalue; // get next x value from the test vector
    values.y = expressions.interpolateSpline; // call the interpolation routine to get our interpolated y value
    list.y_test.setitem(values.y, list.x_test.currentindex); // add interpolated y to the test vector
]
/ trialduration = 0
/ branch = [
    if (list.x_test.selectedcount < list.x_test.itemcount) {    // run as often as we have x values in our test vector
        return trial.interpolateSpline;                         // this is just for ease of data output and plotting, it's not a problem to
    };                                                            // wrap everything in a single while() loop; -> see trial.interpolateSplineLoop
]
/ recorddata = true
</trial>

<trial interpolateSplineLoop>
/ ontrialbegin = [
    values.tStartValueGeneration = script.elapsedtime;
    while (list.x_test.selectedcount < list.x_test.itemcount) {
        values.x = list.x_test.nextvalue;
        values.y = expressions.interpolateSpline;
        list.y_test.setitem(values.y, list.x_test.currentindex);
    };
    values.tEndValueGeneration = script.elapsedtime;
]
/ validresponse = (0)
/ trialduration = 0
</trial>

<data>
/ columns = (script.elapsedtime list.x.itemcount list.y.itemcount list.x_test.itemcount list.y_test.itemcount
    expressions.runtimeTestVectorSetup, expressions.runtimeConstructionRoutine expressions.runtimeInterpolationRoutine expressions.runtimeValueGeneration
    trialnum values.x values.y) // write out x and y one value pair per row for easy plotting
/ separatefiles = true
</data>

Interpolation over 30 input points (small circles in the plot below), and obtaining a total of 2561 interpolated y values for the x range -10 to 70, in 0.03125 increments (the tiny grey dots in the plot below).


Setting up the test vector (2561 x values for which to obtain interpolated y values) typically takes less than 10ms on my system.
Spline creation over the 30-point input space takes around 40ms.
Obtaining the 2561 interpolated y values when done in a single loop takes around 80ms, that's about 0.03ms per y value.

Here's one typical result:

elapsedtime:158
list.x.itemcount: 30
list.y.itemcoun: 30
list.x_test.itemcount: 2561
list.y_test.itemcount: 2561
runtimeTestVectorSetup: 7
runtimeConstructionRoutine: 43
runtimeValueGeneration: 78

Overall, I'd say that's pretty fast.