In both Matlab and Numpy, your indexed loops are not the right way to express vectorized operations. You should adopt a strong bias against loops and individual element indexing, learn how broadcasting works, and pass around ndarrays instead of lists.
Based on the above, a potential vectorized solution - with a regression test to ensure that the results have not changed - is
from numbers import Real
from typing import Sequence, Union
import numpy as np
from numpy.linalg import norm
from scipy.stats import describe
EPSILON = 1e-12
def g(xy: np.ndarray) -> np.ndarray:
return -xy * np.sum(xy**2, axis=0)
def wendland(r: np.ndarray) -> np.ndarray:
r = norm(r, axis=0)
out = (1 - r)**4 * (1 + 4*r)
out[r >= 1] = 0
return out
def make_collocation_points(h: Real, x1: int, x2: int, y1: int, y2: int) -> np.ndarray:
j = np.arange(x1, x2+1)
k = np.arange(y1, y2, h)
jk = np.stack(np.meshgrid(j, k)) * h
# Don't accept values on the chain recurrent set!
g_norm = norm(g(jk) - jk, axis=0)
jk = jk[:, g_norm > 1e-5]
return jk
def get_aij(xy: np.ndarray) -> np.ndarray:
# broadcast to get the effect of xj, xk, yj, yk
xyj, xyk = np.broadcast_arrays(
xy[:, np.newaxis, :],
xy[:, :, np.newaxis],
)
gj = g(xyj)
gk = g(xyk)
return (
+ wendland(gj - gk) - wendland(xyj - gk)
- wendland(gj - xyk) + wendland(xyj - xyk)
)
def assert_close(
a: Union[Real, np.ndarray],
b: Union[Real, Sequence[Real]],
) -> None:
assert np.all(np.isclose(a, b, rtol=0, atol=EPSILON))
def regression_test() -> None:
h = 0.11
n = 2 # 15 is way too slow
x1, x2 = -n, n
y1, y2 = -n, n
xy = make_collocation_points(h, x1, x2, y1, y2)
assert xy.shape == (2, 185)
stats = describe(xy.T)
assert_close(stats.minmax[0], (-0.22, -0.22))
assert_close(stats.minmax[1], (+0.22, +0.2156))
assert_close(stats.mean, (0, -2.2e-3))
assert_close(stats.variance, (0.02433152173913048, 0.016781450543478273))
a = get_aij(xy)
assert a.shape == (185, 185)
stats = describe(a.flatten())
assert_close(stats.minmax[0], -0.17266964737638246)
assert_close(stats.minmax[1], +1.10925185890418380)
assert_close(stats.mean, 0.11771392866408713)
assert_close(stats.variance, 0.061445578720171375)
def demo() -> None:
h = 0.11
n = 10
x1, x2 = -n, n
y1, y2 = -n, n
xy = make_collocation_points(h, x1, x2, y1, y2)
print(describe(xy.T))
a = get_aij(xy)
print(describe(a.flatten()))
if __name__ == "__main__":
regression_test()
demo()
I don't have a whole lot of patience for long execution runs, but for n=10 the demo finishes in about five seconds even though it produces 14,607,684 elements.
