How to get the gradient of the specific points？

[WillingWeiling]

Dear DanielVandH
I apologize for reaching out to you again. I am writing to seek your guidance on a specific issue regarding obtaining the gradient of a specific point after interpolating a model. Our project has made use of the interpolation method you provided, and we have successfully obtained the model. However, we are currently facing a challenge in retrieving the gradient information of a specific point. It is crucial for us to access the gradient value of that point in order to further refine our algorithms and applications.
The code we have been using is as follows:
itpF(0.0115, 0.1, method=Farin(1), project=false)
this code only returns the Z value of the point,We are wondering if there is any way to retrieve the gradient information with this point.

Thank you for your time and support. We look forward to your response!

Best regards，
Weiling

[DanielVandH]

Hi @WillingWeiling:

I apologize for reaching out to you again.

No problem, happy to help.

We are wondering if there is any way to retrieve the gradient information with this point.

There is no way to obtain this information jointly, but you can obtain a differentiator by differentiating the interpolant. Here is an example, reusing the example from #23:

using NaturalNeighbours 
f = (x, y) -> x^2 + y^2 
df = (x, y) -> (2x, 2y)
x = rand(100)
y = rand(100)
z = f.(x, y)
∇z = df.(x, y) 
itp = interpolate(x, y, z; gradient=∇z)
diffr = differentiate(itp, 1) # or do 2 if you need Hessians also
Natural Neighbour Differentiator
    Order: 1
    z: [1.339810600563677, 0.210456011469002, 0.19118917959543236, 0.5083789692753129, 0.26924336271818383, 0.1514543170757532, 1.2166013926394235, 0.9467794665258464, 1.0268012587999824, 0.35642945216289545  …  0.7322415196082168, 0.22467864978768531, 0.9848781877861323, 0.3013936032790566, 0.07388950833082, 0.07243391092378995, 0.6831778709261185, 0.942251976654646, 1.0163586239719593, 0.4238810850380472]
    ∇: [(1.817110640858952, 1.4343470016463509), (0.6849231218672851, 0.6104951785293478), (0.814659090009161, 0.31794226747504806), (0.8732963723445666, 1.127328400755996), (0.7678882165001901, 0.6980839045794514), (0.17126252767514294, 0.7592670247793813), (1.9300977049817554, 1.0682361254806239), (1.0602368101234436, 1.6318749255266627), (1.9299247656242116, 0.6185429930330322), (0.8942733270628636, 0.7911972100276254)  …  (1.6491675747393493, 0.45739740801757955), (0.8878733483215948, 0.3322582075599356), (1.981582480143811, 0.1133297204251027), (0.8010839801079548, 0.750892049452266), (0.0790826594764864, 0.5378698414053404), (0.3955457844924315, 0.36507420624501363), (0.8939055455681149, 1.3905554139648821), (1.610200373625206, 1.0845564362428683), (0.41906308466449227, 1.972262818936492), (0.9663879838183271, 0.8727076285238589)]
    H: nothing

This diffr is what you need. The relevant information can be found using ?differentiate, which lists four methods:

help?> differentiate
search: differentiate

  differentiate(itp::NaturalNeighboursInterpolant, order)

  Differentiate a given interpolant itp up to degree order (1 or 2). The returned object is a
  NaturalNeighboursDifferentiator struct, which is callable.

  For calling the resulting struct, we define the following methods:

  (∂::NaturalNeighboursDifferentiator)(x, y, zᵢ, nc, id::Integer=1; parallel=false, method=default_diff_method(∂), kwargs...)
  (∂::NaturalNeighboursDifferentiator)(x, y, id::Integer=1; parallel=false, method=default_diff_method(∂), interpolant_method=Sibson(), rng=Random.default_rng(), project = true, kwargs...)
  (∂::NaturalNeighboursDifferentiator)(vals::AbstractVector, x::AbstractVector, y::AbstractVector; parallel=true, method=default_diff_method(∂), interpolant_method=Sibson(), kwargs...)
  (∂::NaturalNeighboursDifferentiator{I, O})(x::AbstractVector, y::AbstractVector; parallel=true, method=default_diff_method(∂), interpolant_method=Sibson(), kwargs...) where {I, O}

    1. This method is useful if you already have an estimate for the function value, zᵢ, at the data site, (x,
       y), provided you also provide the NaturalCoordinates nc. id is the thread id.

    2. This method is for scalars, with id referring to a thread id.

    3. This method is an in-place method for vectors, storing ∂(x[i], y[i], 1) into vals[i].

    4. This method is similar to (3), but vals is constructed and returned.

  The available keyword arguments are:

    •  parallel=true: Whether to use multithreading. Ignored for the first two methods.

    •  method=default_diff_method(∂): Default method for evaluating the interpolant. default_diff_method(∂)
       returns Direct(). The method must be a AbstractDifferentiator.

    •  interpolant_method=Sibson(): The method used for evaluating the interpolant to estimate zᵢ for the latter
       three methods. See AbstractInterpolator for the available methods.

    •  rng=Random.default_rng(): The random number generator used for estimating zᵢ for the latter three methods,
       or for constructing the natural coordinates.

    •  project=false: Whether to project any extrapolated points onto the boundary of the convex hull of the data
       sites and perform two-point interpolation, or to simply replace any extrapolated values with Inf, when
       evaluating the interpolant in the latter three methods.

    •  use_cubic_terms=true: If estimating second order derivatives, whether to use cubic terms. Only relevant
       for method == Direct().

    •  alpha=0.1: If estimating second order derivatives, the weighting parameter used for estimating the second
       order derivatives. Only relevant for method == Iterative().

    •  use_sibson_weight=true: Whether to weight the residuals in the associated least squares problems by the
       associated Sibson coordinates. Only relevant for method == Iterative() if order == 2.

  The outputs are:

    •  order == 1: The scalar methods return a Tuple of the form (∂x, ∂y), while the vector methods return a
       vector of Tuples of the form (∂x, ∂y).

    •  order == 2: The scalar methods return a (∇, ℋ), where ∇ is a Tuple of the form (∂x, ∂y) and ℋ is a Tuple
       of the form (∂xx, ∂yy, ∂xy). The vector methods return a vector of (∇, ℋ)s.

  │ Warning
  │
  │  Until we implement ghost point extrapolation, behaviour near the convex hull of your data sites may in
  │  some cases be undesirable, despite the extrapolation method we describe above, even for points that are
  │  inside the convex hull. If you want to control this behaviour so that you discard any points that are very
  │  close to the convex hull, see identify_exterior_points and the tol keyword argument.

Here are some examples:

julia> df(0.521, 0.2001) # exact
(1.042, 0.4002)

julia> diffr(0.521, 0.2001)
(1.111846841430759, 0.49196370869830647)

julia> diffr(0.521, 0.2001, interpolant_method = Farin(1)) # this one is much better than above!
(1.042, 0.40020000000000017)

julia> diffr(0.521, 0.2001, interpolant_method = Farin(1), project = false)
(1.042, 0.40020000000000017)

julia> diffr(0.521, 0.2001, interpolant_method = Farin(1), project = false, method = Direct())
(1.042, 0.40020000000000017)

julia> diffr(0.521, 0.2001, interpolant_method = Farin(1), project = false, method = Iterative())
(1.042, 0.40020000000000017)

I hope that helps. Feel free to let me know if you've any more questions.

[DanielVandH]

(Note that, just like for interpolation, it is much faster to evaluate the differentiator at vector inputs rather than scalar inputs, so that multithreading can be used to evaluate in batches.)

[WillingWeiling]

@DanielVandH Hi,
Thank you so much for your prompt reply and helpful assistance!! You've really cleared up my confusion! And I'm grateful for your kind advice on the computational efficiency of numerical differentiation through multithreading. I appreciate the reminder.
Wishing you all the best！：）

Best regards，
Weiling

[WillingWeiling]

Following your instructions, I utilized the following code snippet:
Gsx = Gs[:,1]
Gsy = Gs[:,2]
Gsz = Gs[:,3]
dz = [(row[4], row[5]) for row in eachrow(Gs)] #Gs is a matrix containing data information.
itpF = NaturalNeighbours.interpolate(model.Gsx, model.Gsy, model.Gsz; gradient=model.dz)
diffr = differentiate(itpF, 1)
F = itpF(ϕA, ϕB, method=Nearest(), project=false)
mu = diffr(ϕA, ϕB, interpolant_method=Nearest(), project=false, method=Direct()) /mu = diffr(ϕA, ϕB, interpolant_method=Nearest(), project=false, method=Iterative()) The results obtained by selecting method = Iterative() or Direct() are the same.
My objective is to obtain an interpolated convex surface and analyze a cross-section by selecting ten points. However, upon plotting the results, I observed a certain shift in the interpolated gradient positions compared to the true data points.

The horizontal axis in the graph represents x-values, while the vertical axis corresponds to the interpolated gradient data obtained from diffr.

I am wondering if this discrepancy may be attributed to the fact that the provided gradient=model.dz was not effectively utilized during the gradient calculation. I would greatly appreciate your insights into this matter and any suggestions you may have to address this issue.