Created
January 19, 2022 12:06
-
-
Save PageotD/c79a1e93da8fa917f807100478a6ebd2 to your computer and use it in GitHub Desktop.
Sambridge's Neighborhood Algorithm (Coffee Break Scripting)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment