TreeMesh over a 2D function?

[psteinb]

Hi, thanks for offering discretize to the open-source community. I am currently looking into optimally sampling a 2D function. I had hopes that discretize might help here. I already sat down and went through the tutorial on TreeMesh. As I am new to the field, I apologize for my potentially naive question. Let's assume I have an arbitrary 2D function on a cartesian grid like so:

import numpy as np

x = np.linspace(-4,4,100)
y = np.linspace(-4,4,100)

xm, ym = np.meshgrid(x,y)

z = xm*np.sin(xm**2 + ym**2) / (xm**2 + ym**2)

How would I use TreeMesh to find an optimal discretization of z with less sample points than the default cartesian sampling used here with np.meshgrid? My goal is to minimize the number of samples over the plane [-4,4]x[-4,4]. The examples speak of 2D, but effectively they deal with a function y = f(x).

[domfournier]

Hi @psteinb,

You could refine based on the z value, mesh.refine(f(x, y))

I would suggest that you probably want to refine on the gradient of the function instead. Here is a notebook with a slightly simpler function (just cause I didn't feel like taking long derivatives by hand).

Refine2DTreeMesh.zip

You get the idea.

[psteinb]

Thanks for providing the solution so quickly and explicitely. In my use case, I don't have an analytical description of the function in question. ;) So I'd have to approximate the gradient along z within the cell, i.e. your zfunGrad function.

Btw, I was confused by this part of the documentation a bit. What is x0 (bottom west corner as I found in the examples)? Why is the width denoted as h as in height (is that from the usual gradient formula?)?

Related to my question, could I approximate the gradient then by the difference between cell.center and cell.x0 are? E.g. like so:

def refineZgradient2(cell):
    delta = cell.center - cell.x0
    dz = (zfun(*(cell.x0)) - zfun(*(cell.center + delta)))/ 2*norm(cell.h)
    return int(mesh.max_level * np.abs(dz))

This gives decent results, but they are far from perfect.

[domfournier]

Guess you could yes. I would probably be faster (and safer) to compute your gradients over a fixed length distance. You might also have to normalize your gradients on the overall max though so that the return index always falls between 0 and mesh.max_level. It worked in the example above because the trig functions are between [0, 1]. Yours to experiment with.

[psteinb]

thanks for your excellent help. I'll have to understand better how the algorithm works. I'll start digging the source code. If there is any other source of information, please let me know.