common.py
from numba import njit
from numpy import empty
from math import sqrt
from random import uniform
@njit
def vec2(x,y):
'''Convert two numbers to a 2D numpy array '''
v = empty((2,),dtype=float)
v[0] = x
v[1] = y
return v
@njit
def vec3(x,y,z):
'''Convert three numbers to a 2D numpy array '''
v = empty((3,),dtype=float)
v[0] = x
v[1] = y
v[2] = z
return v
@njit
def vec4(x,y,z,w):
'''Convert four numbers to a 2D numpy array '''
v = empty((4,),dtype=float)
v[0] = x
v[1] = y
v[2] = z
v[3] = w
return v
# @njit
# def V(x, inv_noise_scale, offset):
# ''' V returns a curl noise vector field'''
# p = empty(3)
# p[0:2] = x*inv_noise_scale+offset
# return scnoise_curl(p)[0:2]
@njit
def V(x, inv_noise_scale, offset, noise_grad):
''' V returns a curl noise vector field'''
p = empty(3)
p[0:2] = x*inv_noise_scale+offset
p[2] = 98.45
g = noise_grad(p)
gh = vec2(-g[1], g[0])
return gh
@njit
def Euler(x, dt, inv_noise_scale, offset, noise_grad):
'''Perform one Euler step of length dt from x along the curl noise vector field'''
return x + dt * V(x, inv_noise_scale, offset, noise_grad)
@njit
def RK4(x, dt, inv_noise_scale, offset, noise_grad):
'''Perform one Runge-Kutta 4th order step of length dt from x along the curl noise vector field'''
a = dt * V(x, inv_noise_scale, offset, noise_grad)
b = dt * V(x+a/2, inv_noise_scale, offset, noise_grad)
c = dt * V(x+b/2, inv_noise_scale, offset, noise_grad)
d = dt * V(x+c, inv_noise_scale, offset, noise_grad)
return x + (a+2*b+2*c+d)/6
@njit
def rand_vec2():
'''Convert two numbers to a 2D numpy array '''
v = empty((2,),dtype=float)
v[0] = uniform(0,1)
v[1] = uniform(0,1)
return v
## UNUSED
@njit
def smoothstep(a,b,x):
t = x-a/b-a
if t < 0: return 0
if 1 < t: return 1
return 3*t**2-2*t**3
@njit
def Rand(x, dt, inv_noise_scale, offset):
'''Perform a step of length dt from x in a random direction.'''
return x + dt * rand_vec2()
@njit
def NoiseV(x, dt, inv_noise_scale, offset, noise_grad):
'''Perform a step of length dt from x in a noise vector direction.'''
p = vec3(0,0,offset)
p[0:2] = x*inv_noise_scale
g = noise_grad(p)
return x + dt * g[0:2]
def L2_star_disc(P):
n = P.shape[0]
l2sd = (1/3)**2
l2sd -= (2/n) * sum([ 0.25 * (1 - x[0]**2) * (1 - x[1]**2) for x in P])
s = 0
for xi in P:
for xk in P:
s += (1-max(xi[0], xk[0]))*(1-max(xi[1], xk[1]))
l2sd += s/(n**2)
return sqrt(l2sd)