"""some more complex examples of actual use cases than found in unit tests"""
from __future__ import absolute_import, division, print_function, unicode_literals
from builtins import *
import matplotlib.pyplot as plt
from numpy_indexed import *
__author__ = "Eelco Hoogendoorn"
__license__ = "LGPL"
__email__ = "hoogendoorn.eelco@gmail.com"
def test_radial_reduction():
"""radial reduction example"""
x = np.linspace(-2,2, 64)
y = x[:, None]
x = x[None, :]
R = np.sqrt(x**2+y**2)
def airy(r, sigma):
from scipy.special import j1
r = r / sigma * np.sqrt(2)
a = (2*j1(r)/r)**2
a[r==0] = 1
return a
def gauss(r, sigma):
return np.exp(-(r/sigma)**2)
distribution = np.random.choice([gauss, airy])(R, 0.3)
sample = np.random.poisson(distribution*200+10).astype(np.float)
#is this an airy or gaussian function? hard to tell with all this noise!
plt.imshow(sample, interpolation='nearest', cmap='gray')
plt.show()
#radial reduction to the rescue!
#if we are sampling an airy function, you will see a small but significant rise around x=1
g = group_by(np.round(R, 5).flatten())
plt.errorbar(
g.unique,
g.mean(sample.flatten())[1],
g.std (sample.flatten())[1] / np.sqrt(g.count))
plt.xlim(0, 2)
plt.show()
def test_mesh_solver():
"""
meshing example
demonstrates the use of multiplicity, and group.median
"""
#set up some random points, and get their delaunay triangulation
points = np.random.random((20000,2))*2-1
points = points[np.linalg.norm(points,axis=1) < 1]
from scipy.spatial.qhull import Delaunay
d = Delaunay(points)
tris = d.simplices
#the operations provided in this module allow us to express potentially complex
#computational geometry questions elegantly in numpy
#Delaunay.neighbors could be used as well,
#but the point is to express this in pure numpy, without additional library functionaoty
edges = tris[:,[[0,1],[1,2],[2,0]]].reshape(-1,2)
sorted_edges = np.where(edges[:,0:1]<edges[:,1:2], edges, edges[:,::-1])
#we test to see how often each edge occurs, or how many indicent simplices it has
#this is a very general method of finding the boundary of any topology
#and we can do so here with only one simple and readable command, multiplicity == 1
if semantics.backwards_compatible:
boundary_edges = edges[multiplicity(sorted_edges, axis=0)==1]
else:
boundary_edges = edges[multiplicity(sorted_edges)==1]
boundary_points = unique(boundary_edges)
if False:
print(boundary_edges)
print(incidence(boundary_edges))
#create some random values on faces
#we want to smooth them over the mesh to create a nice hilly landscape
face_values = np.random.normal(size=d.nsimplex)
#add some salt and pepper noise, to make our problem more interesting
face_values[np.random.randint(d.nsimplex, size=10)] += 1000
#start with a median step, to remove salt-and-pepper noise
#toggle to mean to see the effect of the median filter
g = group_by(tris.flatten())
prestep = g.median if True else g.mean
vertex_values = prestep(np.repeat(face_values, 3))[1]
vertex_values[boundary_points] = 0
#actually, we can compute the mean without grouping
tris_per_vert = g.count
def scatter(x):
r = np.zeros(d.npoints, x.dtype)
for idx in tris.T: np.add.at(r, idx, x)
return r / tris_per_vert
def gather(x):
return x[tris].mean(axis=1)
#iterate a little
for i in range(100):
face_values = gather(vertex_values)
vertex_values = scatter(face_values)
vertex_values[boundary_points] = 0
#display our nicely rolling hills and their boundary
x, y = points.T
plt.tripcolor(x, y, triangles = tris, facecolors = face_values)
plt.scatter(x[boundary_points], y[boundary_points])
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.axis('equal')
plt.show()
