Documentation for derivatives

[Mv77]

Hi,

I've recently been working on wrappers for this repos' interpolators in HARK.

I now want to add the capability of computing derivatives. I have been trying to figure out how to use the derivatives in eval_splines from your test files. I've now run into a question I can't resolve and I was hoping you could provide some guidance.

Consider the following example, which calculates the 0-1-2 order derivatives of f(x) = 2 + 1*x at some points (1.5, 2.5, 3.5, 4.5):

x = np.linspace(0,10,6)
y = 2 + 1.0 * x

grad = eval_spline(
    UCGrid(x),
     y,
     np.array([1.5,2.5,3.5,4.5])[...,None],
     out=None,
     order=1,
     diff=str(((0,), (1,), (2,))),
     extrap_mode="linear"
)

The output is

array([[3.5, 4. , 0. ],
       [4.5, 4. , 0. ],
       [5.5, 4. , 0. ],
       [6.5, 4. , 0. ]])

Which has some properties of the right solution: the first and third column are right, and the second column rightly is a constant. But, in the second column, I'd expect the first order derivative which should be 1.0.

Could you help me figure out if there is something I am misinterpreting?

[aniecon]

Hi Mateo, I tried this as well. And it gave me the correct first derivative only when the grid was set to x = np.linspace(1.5,4.5,4).

[llorracc]

Pablo,

Some context for this: We've embarked on a goal of wiring your interpolation.py into the core of HARK, but can't do that until we figure out this issue. It sure looks like a bug to me; or, if not a bug, then probably there is some place where the documentation should be improved. Is there someone on your team who can look at this, maybe with collaboration/help from @Mv77 and/or @alanlujan91

[albop]

Just saw the thread. I'll have a look at it somewhere next week. It sure looks like a bug, but I need to get back into the logic of it...

…

On Fri, Jul 22, 2022 at 11:05 PM Christopher Llorracc Carroll < ***@***.***> wrote: Pablo, Some context for this: We've embarked on a goal of wiring your interpolation.py into the core of HARK, but can't do that until we figure out this issue. It sure looks like a bug to me; or, if not a bug, then probably there is some place where the documentation should be improved. Is there someone on your team who can look at this, maybe with collaboration/help from @Mv77 <https://github.com/Mv77> and/or @alanlujan91 <https://github.com/alanlujan91> — Reply to this email directly, view it on GitHub <#94 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AACDSKMYZPM4VOWEWSRRUNLVVMEH7ANCNFSM53D7U4ZQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[albop]

Indeed, there was a bug with the evaluation of derivatives for linear interpolators.
I have pushed a fix in master. It should in principle work for uniform and nonuniform grids.

@Mv77 or @aniecon : can you check it works as you expect ? Would you submit a test file ? In principle, it should work for regular and irregular grids, linear and cubic splines. When we are sure it works we can add it to the doc.

NB: a grid is represented as a tuple of 1d grid. A 1d grid can itself be regular, it which case it is given by (min, max, n) or irregular in which case it must be a numpy vector of increasing coordintates. In your example, there is a hidden bug: UCGrid should not accept x as argument. Instead, you should do UCGrid( (0.,10.,11) ) which leads to faster interpolation. And to construct nonuniform grid CGrid( x ) (then x can be either a vector or a min,max,n, tuple). I have pushed another little change so that CGrid, and UCGrid actually do some type check (in the long run my intention was to use python type system for that but so far it will do).

NB2: to debug I used the call code = (get_code_spline(1,1,False,True,allocate=True,extrap_mode="linear", orders=((0,), (1,), (2,)))) to see what code had been generated and noticed 1.0*δ_0 instead of 1.0/δ_0 in the derivative of the basis functions.

[Mv77]

@albop thank you for the fix and for the tips.

I'd be happy to create a test file. I'll submit a PR in the next few days and confirm that it works.

[Mv77]

@albop I'm already working on tests and should have a PR soon. It's looking good. Thank you!

Just a question: I am using order=1 interpolators. First order derivatives look right. Second order derivatives are always coming out as 0. Just wanted to check that this is expected behavior (second derivative of a linear function), and that the interpolators are not expected to get a 2nd order derivative from some other numerical approach.

[llorracc]

So, i guess the tricky thing is that for piecewise linear functions the second derivative is zero everywhere except the kink points, where it is +- infty. ?

On Fri, Jul 29, 2022 at 8:28 PM Mateo Velásquez-Giraldo < ***@***.***> wrote: @albop <https://github.com/albop> I'm already working on tests and should have a PR soon. It's looking good. Thank you! Just a question: I am using order=1 interpolators. First order derivatives look right. Second order derivatives are always coming out as 0. Just wanted to check that this is expected behavior (second derivative of a linear function), and that the interpolators are not expected to get a 2nd order derivative from some other numerical approach. — Reply to this email directly, view it on GitHub <#94 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAKCK75QDPXYVU233BI3BQTVWRZJZANCNFSM53D7U4ZQ> . You are receiving this because you commented.Message ID: ***@***.***>

-- Sent from Gmail Mobile

[albop]

I know that one is closed, but just to comment after death: the constructed splines are differentiable from the right up to order k and to order k-1 on the left. So linear splines are differentiable everywhere from the right, but not differentiable from the left at the knots points. This is actually the expected mathematical result, even though an undesirable consequence is the lack of symmetricity.

…

On Sat, Aug 6, 2022 at 8:21 PM Mateo Velásquez-Giraldo < ***@***.***> wrote: Closed #94 <#94> as completed. — Reply to this email directly, view it on GitHub <#94 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AACDSKN6LYCVK7UNG3INIPLVX2UKZANCNFSM53D7U4ZQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[llorracc]

So linear splines are differentiable everywhere from the right,
but not differentiable from the left at the knots points. This is actually
the expected mathematical result

That's surprising to me. It seems as though it should either be defined both from right and from left with slopes of the corresponding segments, or be undefined from either direction because it is undefined exactly at the knot.

If you do not specify the direction you are coming from, I guess it returns the derivative from the right? To prevent code from breaking when people use it without thinking about the knots point?

In practice, because actual kink points can be so ornery, it has occurred to me to handle this by adding two points to the grid symmetrically slightly to the left and slightly to the right, and connecting them with the unique quadratic that matches the (well-defined) slopes at the left and right. This would yield a function that was differentiable everywhere, and that would preserve the concavity (or convexity) of the function; and with the disadvantage that the quadratic would not yield exactly the right answer exactly at the kink point.