Replace spline function

[dricki]

• All source code must be compatible with the GNU Lesser General Public License version 3, as that is the target license for release. In particular, the GPL-spline implementation in current evalresp releases or any other GPL-licensed code cannot be used.

• There should be no dependencies on purchasing COTS software to satisfy solutions.

[he-hesce]

Resolved in pull request #17

[he-hesce]

I would recommend test building the new code on Windows.

[he-hesce]

Future work:

• rename function from evr_spline() evalresp_spline() (c file, header file, resp_fctns.c)
• rename evr_spline.h to evalresp_splineP.h (header file, resp_fctns.c)
• remove k and tension arguments from function declaration (header file), function definition (c file), and function call (resp_fctns.c)

[dricki]

Henrik - we found a case of not so good spline.

I am attaching the AMPITUDE resposnse BEFORE and AFTER spline modification with X=counts and X=Freqs

The results are generated by the 01_args.robot (test is errorneously named B62 interpolation)

Please see if you can explain/fix the bad jerks at first 10 samples there

[dricki]

Here are the screenshots

[he-hesce]

Old routine has some functionality which ignored "out of range" values. COuld be it.

[dricki]

Ok, interesting. Are you going to look further how it worked and if we can implement it? Indeed, I think the problem could be because interpolation was done on the log-spaced dataset

[he-hesce]

However, this is still done in the wrapper interpolate_list_blockette() function.

[he-hesce]

I will look if the old routine did something extra special than the replacement spline interpolation is not doing.

[he-hesce]

I am finishing documentation first.

[andrewcooke-isti]

i've updated master so that test now has a tolerance of 0.1 (10%) rather than the 1e-8 of other tests.

feel free to modify to a tighter value if appropriate (it still fails).

[andrewcooke-isti]

when is the "future work" above going to happen? it's going to clash with refactoring (i am about to look at command line args, including tension, for example). can i just go ahead and do it?

[andrewcooke-isti]

or maybe i am going about this refactoring wrong?

[he-hesce]

Yepp. The future work is for the refactoring or general clean-up of naming which has been requested. You can just go ahead and do it.

[he-hesce]

It was requested that everything be named evalresp and not evr or evresp.

[he-hesce]

Hmm, the first 8 points looks horrible with the new code. So it could be that the algorithm needs to look at more points initially.

[he-hesce]

The issue is the boundary condition flags when calculating the second derivatives.

    Input, int IBCBEG, left boundary condition flag:
    0: the cubic spline should be a quadratic over the first interval;
    1: the first derivative at the left endpoint should be YBCBEG;
    2: the second derivative at the left endpoint should be YBCBEG;
    3: Not-a-knot: the third derivative is continuous at T(2).

    Input, double YBCBEG, the values to be used in the boundary
    conditions if IBCBEG is equal to 1 or 2.

    Input, int IBCEND, right boundary condition flag:
    0: the cubic spline should be a quadratic over the last interval;
    1: the first derivative at the right endpoint should be YBCEND;
    2: the second derivative at the right endpoint should be YBCEND;
    3: Not-a-knot: the third derivative is continuous at T(N-1).

    Input, double YBCEND, the values to be used in the boundary
    conditions if IBCEND is equal to 1 or 2.

When using IBCBEG = 0 YBCBEG = 0.0 [not used] IBCEND = 0 YBCEND = 0.0 [not used]

When using IBCBEG = 1 YBCBEG = 0.0 IBCEND = 1 YBCEND = 0.0

When using IBCBEG = 2 YBCBEG = 0.0 IBCEND = 2 YBCEND = 0.0

When using IBCBEG = 3 YBCBEG = 0.0 [not used] IBCEND = 3 YBCEND = 0.0 [not used]

---

So IBCBEG = 2 is the best result. Maybe one could tweak YBCBEG = 0.0 to improve it further.

[he-hesce]

IBCBEG = 2 but with logarithmic scale on x axis only.

[he-hesce]

Clearly even with IBCBEG = 2, the spline is pretty horrible. It is more of a straight line. Particularly at the first 10 data points.

[he-hesce]

Using spline_cubic_val() instead of spline_cubic_val2() appears to resolve the issue.

[he-hesce]

Pushed use of spline_cubic_val() and IBCBEG = 0, IBCEND = 0 as commit 146f7094c1400ddcce7645fd3fc2e7b3bf0a3c86

[he-hesce]

Using spline_cubic_val() instead of spline_cubic_val2()

When using IBCBEG = 0 IBCEND = 0

When using IBCBEG = 1 YBCBEG = 0.0 IBCEND = 1 YBCEND = 0.0

When using IBCBEG = 2 YBCBEG = 0.0 IBCEND = 2 YBCEND = 0.0

When using IBCBEG = 3 IBCEND = 3

[he-hesce]

You may want to look at if you want to use IBCBEG = 1 or IBCBEG = 2 instead of IBCBEG = 0. IBCBEG = 3 seems too optimistic on that spline curvature.

[he-hesce]

IBCBEG = 2 being closest to target (old cubic spline). For this example.

[dricki]

OK, I see that the solution is IBCEG =2 and spline_cubic_val(). Please also check how the phase looks and if the current solution pass the test, please close the issue. Thanks. Really good investigation

[he-hesce]

I will change to IBCBEG=2 and IBCEND=2, and use YBCBEG = 0.0 and YBCEND = 0.0.

[he-hesce]

spline_cubic_val()

IBCBEG = 2 YBCBEG = 0.0 IBCEND = 2 YBCEND = 0.0

AMPLITUDE logarithmic on x axis:

PHASE logarithmic on x axis:

[he-hesce]

On the tests, I get the following for the C tests:

check_response.c: In function ‘test_response’:
check_response.c:77: error: too few arguments to function ‘print_resp’

So failure during compilation of tests. Ignoring that.

And for Robot tests I get on the relevant test:

------------------------------------------------------------------------------
B62 interpolation                                                     | FAIL |
Keyword 'Support.Count And Compare Target Files Two Float Cols' expected 0 to 1 arguments, got 2.
------------------------------------------------------------------------------

master branch.

[andrewcooke-isti]

i'll have a look at those in a moment.

[andrewcooke-isti]

ok, first is dylan not changing tests to match his new logging. second is probably because you are re-using an existing env, but the robot test library was changed (sorry - i don't rebuild every time because it means making a new env, downloading stuff etc). please "rm -fr env" in your evalresp directory and run the tests again. i see more failures than before, some of which are likely related to logging changes. (is there a good reason why dylan is committing changes before the work is complete?)

[he-hesce]

Ok, test runs... but probably the new numbers are diffing too much:

B62 interpolation                                                     | FAIL |
166705.000000 and 149728.900000 differ at /home/xxx/git/evalresp/tests/robot/run/base/args/b62/AMP.IM.ATTU..BHE and /home/xxx/git/evalresp/tests/robot/target/base/args/b62/AMP.IM.ATTU..BHE at line 1

[andrewcooke-isti]

close. you can adjust the tolerance in the test. it's currently 10% (0.1). maybe 20% would be ok...

[he-hesce]

The first 8 points (lines) still differ more than the rest. It's the new spline algorithm finding its "spline" footing.

[he-hesce]

Basically, after the first 8 it is much much closer.

[he-hesce]

From 12 points/line, it is the same value. So the issue is how the spline algorithm handles the beginning (and end) of the spline.

[andrewcooke-isti]

do we need an extra param to skip lines?! that's a little clunky. maybe call out to a shell task and run head or tail? also clunky but at least not enshrined in the library.

[andrewcooke-isti]

typically splines have boundary conditions that are effectively gradients. sounds like those are wrong.

[he-hesce]

See discussion above what I have tried concerning the beg/end conditions.

[he-hesce]

    Input, int IBCBEG, left boundary condition flag:
    0: the cubic spline should be a quadratic over the first interval;
    1: the first derivative at the left endpoint should be YBCBEG;
    2: the second derivative at the left endpoint should be YBCBEG;
    3: Not-a-knot: the third derivative is continuous at T(2).

    Input, double YBCBEG, the values to be used in the boundary
    conditions if IBCBEG is equal to 1 or 2.

    Input, int IBCEND, right boundary condition flag:
    0: the cubic spline should be a quadratic over the last interval;
    1: the first derivative at the right endpoint should be YBCEND;
    2: the second derivative at the right endpoint should be YBCEND;
    3: Not-a-knot: the third derivative is continuous at T(N-1).

    Input, double YBCEND, the values to be used in the boundary
    conditions if IBCEND is equal to 1 or 2.

[he-hesce]

We are currently using IBCBEG = 2 YBCBEG = 0.0 IBCEND = 2 YBCEND = 0.0

As that produces the "best" (closest to old spline) results of IBCBEG 1,2,3,4.

But perhaps YBCBEG (and YBCEND) can be changed from 0.0 to something else.

[andrewcooke-isti]

maybe log tolerance rather than linear?

[andrewcooke-isti]

did you try 20%? from those numbers it should be ok.

[he-hesce]

No, didn't try 20% as I think that is a bit much of a margin. Would be better if we could fix the boundary condition. But I am not a spline nor math expert, so I would need help with that.

[andrewcooke-isti]

tbh i think the problem may be with what we had, rather than what you added. it was a slightly non-standard approach that included "tension", which was set at some arbitrary value. i don't think a "normal" spline (which i guess you used) is going to reproduce that exactly.

when i looked at this i thought i found the code that eric used. traced it back to some old fortran paper? i can't remember any details, but i suspect you'd need to duplicate that to make it match.

[andrewcooke-isti]

sorry i should perhaps have explained all that before you started, but i doubt anyone told me you would be doing it.

[he-hesce]

Sure, could be it. So Ilya can sign off on updating the target results if he is happy with the results above. I guess we could plot the input values (x,y) and see how they compare to the new (x2,y2) interpolated values.

[dricki]

I am suprised that 10% is not enough. What I've seen on figures was really good and close to the spline produced by Eric's code within a couple of percents.....

[he-hesce]

@andrewcooke-isti the old code was from GNU plotutils.

[he-hesce]

@andrewcooke-isti a very very old version from the 90s. I had to dig in the GNU archives to find it. Has been replaced in plotutils since 15+ years or so.

[andrewcooke-isti]

before that it was in a paper. i traced it back when quoting for the work.