1
\$\begingroup\$

I previously posted some code I've been working on and had a fantastic answer, but as I didn't post the full code I was then unable to bring it all together using the much faster numpy implementation, so a review of the full code is much appreciated. I now realise that using vectorised operations would be much, much faster, but have little experience.

import numpy as np
from numpy import linalg as LA

def g(x, y):
    gx = -x * (x ** 2 + y ** 2)
    gy = -y * (x ** 2 + y ** 2)
    return np.array((gx, gy))

def wendland(r):
    if r < 1:
        return ((1-r)**4)*(4*r+1)
    else:
        return 0

def make_collocation_points(h, x1, x2, y1, y2):
    i = 0
    x = []
    y = []

    for j in range(x1,x2+1):
        for k in np.arange(y1, y2, h):
            if LA.norm(g(j*h, k*h) - (j*h, k*h)) > 0.00001: # Don't accept values on the chain recurrent set!
                x.append(j*h)
                y.append(k*h)

                i += 1

    return x, y, i

def get_aij(x, y, n):
    a = np.zeros((n,n))
    for j in range(0, n):
        for k in range(0, n):
            a[j, k] = (
                wendland(LA.norm(g(x[j], y[j]) - g(x[k], y[k])))
                - wendland(LA.norm(g(x[j], y[j]) - np.array([x[k], y[k]])))
                - wendland(LA.norm([x[j], y[j]] - g(x[k], y[k])))
                + wendland(LA.norm(np.array([x[j], y[j]]) - np.array([x[k], y[k]])))
            )
    return a

def solve_for_B(a, n):
    r = -1 * np.ones(n)
    return LA.solve(a, r)


def find_v_and_OD(x, y, n):
    x_range = np.arange(x1, x2 + h, h)
    y_range = np.arange(y1, y2 + h, h)
    v_val = np.zeros((len(x_range), len(y_range)))
    delV = np.zeros((len(x_range), len(y_range)))

    for i, r in enumerate(x_range):
        for j, q in enumerate(y_range):
            v_val[i, j] = v(r * h, q * h, x, y)
            delV[i, j] = orbital_derivative(r * h, q * h, x, y)
    return v_val, delV


def v(s, t, x, y):
    out = np.multiply(
        b[0],
        (
            wendland(LA.norm(np.array([s, t]) - g(x[0], y[0])))
            - wendland(LA.norm(np.array([s, t]) - np.array(x[0], y[0])))
        ),
    )
    for k in range(1, len(b)):
        out += np.multiply(
            b[k],
            (
                wendland(LA.norm(np.array([s, t]) - g(x[k], y[k])))
                - wendland(LA.norm(np.array([s, t]) - np.array((x[k], y[k]))))
            ),
        )
    return out


def orbital_derivative(s, t, x, y):
    out = np.multiply(
        b[0],
        (
            wendland(LA.norm(g(s, t) - g(x[0], y[0])))
            - wendland(LA.norm(g(s, t) - np.array((x[0], y[0]))))
            - wendland(LA.norm(np.array([s, t]) - g(x[0], y[0])))
            + wendland(LA.norm(np.array([s, t]) - np.array([x[0], y[0]])))
        ),
    )
    for k in range(1, len(b)):
        out += np.multiply(
            b[k],
            (
                wendland(LA.norm(g(s, t) - g(x[k], y[k])))
                - wendland(LA.norm(g(s, t) - np.array((x[k], y[k]))))
                - wendland(LA.norm(np.array([s, t]) - g(x[k], y[k])))
                + wendland(LA.norm(np.array([s, t]) - np.array([x[k], y[k]])))
            ),
        )
    return out

if __name__ == "__main__":
    h = 0.11
    x1, x2 = -15, 15
    y1, y2 = -15, 15

    x, y, n = make_collocation_points(h, x1, x2, y1, y2)
    a = get_aij(x, y, n)
    b = solve_for_B(a, n)
    v_vals, delV = find_v_and_OD(x, y, n)
\$\endgroup\$
3
  • 1
    \$\begingroup\$ Is C on the table? If not, why not? \$\endgroup\$
    – Reinderien
    Sep 16, 2021 at 17:28
  • \$\begingroup\$ I'm using the Streamlit Python package to create an interactive version of this code, hence the switch from MATLAB to Python (Plus the old MATLAB code was also garbage). I've also never used C before, so thats another reason I guess! \$\endgroup\$
    – Zac L
    Sep 16, 2021 at 18:05
  • 1
    \$\begingroup\$ I'm struggling with this question. You already know what's wrong: the code is barely vectorized at all. You know what to do: vectorize the code. You even have sample code from the previous question showing the requisite techniques for half of the functions. I think it's time for you to invest some effort in learning how to use Numpy correctly, and if you encounter specific issues, those are good questions for Stack Overflow. Said another way: you need to meet us in the middle and show some effort. \$\endgroup\$
    – Reinderien
    Sep 16, 2021 at 23:59

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.