Sambridge's Neighborhood Algorithm (Coffee Break Scripting)

neighborhood.py

import matplotlib.pyplot as plt
import numpy as np

def voronoi(nx, ny, dh, xp, yp, val, xrng, yrng):
    """
    Simple Voronoi 2D interpolation.

    :param n1: number of points in the first dimension of the output model
    :param n2: number of points in the second dimension of the output model
    :param dh: spatial sampling
    :param xp: x-coordinates of points to interpolate
    :param yp: y-coordinates of points to interpolate
    :param val: value of the points to interpolate
    :param xrng: tuple (xmin, xmax)
    :param yrng: tuple (ymin, ymax )
    """

    # Get the number of points
    npts = np.size(xp)

    # Declare output array
    model = np.zeros((nx, ny), dtype=np.float32)

    # Loop over dimensions
    for ix in range(0, nx):
        x = xrng[0]+float(ix)*dh
        for iy in range(0, ny):
            y = yrng[0]+float(iy)*dh
            d = np.sqrt((x-xp[:])**2+(y-yp[:])**2)
            imin = np.argmin(d)
            dmin = d[imin]
            model[ix, iy] = val[imin]

    return model

def z_function(x,y):
    return np.sin(x)*(np.exp(1-np.cos(y))**2)+np.cos(y)*(np.exp(1-np.sin(x))**2)+(x-y)**2

def firstTrial(ns, xmin, xmax, ymin, ymax):
    # First trial
    for i in range(ns):
        x = np.random.random_sample()*(xmax-xmin) + xmin
        y = np.random.random_sample()*(xmax-xmin) + xmin
        v = z_function(x, y)
        if i == 0:
            points = np.asarray([[v, x, y]])
        else:
            points = np.append(points, [[v, x, y]], axis=0)
    return points

# NA parameters
ns = 15 # Number of sample per iteration
nv = 5  # Number of selected voronoi 
niter = 9 # number of iterations
xmin = 0.; xmax = 8. # x range
ymin = 0.; ymax = 8. # y range
nx = 401; ny = 401; dh = 0.02 # Grid

# First Trial
points = firstTrial(ns, xmin, xmax, ymin, ymax)

# Sort
points = points[points[:, 0].argsort()]

# Create model
model = voronoi(nx, ny, dh, points[:,1], points[:, 2], points[:,0], (0., 8.), (0., 8.))

# Calculate model value for each grid point
x= np.linspace(xmin, xmax, nx)
y= np.linspace(ymin, ymax, ny)
X,Y= np.meshgrid(x,y)
Z= z_function(X,Y)

# Plot
plt.subplot(3, 5, 1)
plt.xlabel("x")
plt.ylabel("y")
plt.imshow(model.T, aspect="auto", origin="lower", interpolation="gaussian", cmap="gnuplot_r", extent=[xmin, xmax, ymin, ymax], vmin=-100, vmax=60)
contour = plt.contour(Z, 10, colors = 'white', linewidths = 0.5, extent=[xmin, xmax, ymin, ymax])
plt.clabel(contour, inline=1, fontsize=6)
plt.text(0.5, 7.5, "iter:: 1", fontfamily="sans-serif", color="white", fontweight="bold")

# Loop over iterations
for iiter in range(niter):
    print(iiter)
    # number of sample per Voronoi cells
    nsv = int(ns/nv)
    # Loop over voronoi
    for iv in range(nv):
        #Loop over nsv
        count = 0
        while count != nsv:
            # Trial for the nv best results
            x = np.random.random_sample()*(xmax-xmin) + xmin
            y = np.random.random_sample()*(xmax-xmin) + xmin
            # Calc distance
            d = np.sqrt((x-points[:,1])**2+(y-points[:,2])**2)
            if np.argmin(d) == iv:
                count += 1
                v = z_function(x, y)
                if iv == 0:
                    pointstmp = np.asarray([[v, x, y]])
                else:
                    pointstmp = np.append(pointstmp, [[v, x, y]], axis=0)
    
    points = np.append(points, pointstmp, axis=0)
    
    # Sort
    points = points[points[:, 0].argsort()]

    model = voronoi(nx, ny, dh, points[:,1], points[:, 2], points[:,0], (0., 8.), (0., 8.))

    plt.subplot(3, 5, 2+iiter)
    plt.xlabel("x")
    plt.ylabel("y")
    plt.imshow(model.T, aspect="auto", origin="lower", interpolation="gaussian", cmap="gnuplot_r", extent=[xmin, xmax, ymin, ymax], vmin=-100, vmax=60)
    contour = plt.contour(Z, 10, colors = 'white', linewidths = 0.5, extent=[xmin, xmax, ymin, ymax])
    plt.clabel(contour, inline=1, fontsize=6)
    plt.text(0.5, 7.5, "iter:: "+str(2+iiter), fontfamily="sans-serif", color="white", fontweight="bold")
    
plt.show()