Skip to content

Commit a13cbe4

Browse files
MWillettMMWillettM
authored
Add files via upload
1 parent d6c0145 commit a13cbe4

File tree

3 files changed

+319
-0
lines changed

3 files changed

+319
-0
lines changed

‎DEMO.py‎

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
import numpy as np
2+
import random
3+
import radial_functions as rf
4+
import time
5+
6+
#An implementation of the Faul - Goodsen - Powell algorithm. Given data to interpolate,
7+
#the algorithm returns the coefficents (the lambda_i) and the constant term (alpha)
8+
#for a multiquadric interpolant s(x) of the form
9+
10+
#s(x) = sum_i^n lambda_i phi(x-x_i) + alpha
11+
12+
#This interpolant is accurate to within your prescribed error
13+
#The x_i are the data centers, phi(x) is the multiquadric radial basis fuction (x^2+c^2)^0.5 and
14+
#alpha is a constant term.
15+
16+
#INPUT VARIABLES:
17+
#n is number of data points
18+
#c is the multiquadric parameter >0. Using a smaller value is generally quicker
19+
#q is a variable controlling the size of sets for the cardinal functions generated as part of the algorithm. Usually we use 5<q<50 - q=30 is a safe value to use. The smaller the value of q, the quicker and dirtier the algorithm will be.
20+
#d is the dimension of the data to interpolate.
21+
#error is the error you wish to be within.
22+
#We generate our own random data for this implementation.
23+
24+
# x_i are the centers
25+
#f_i are the function values,
26+
#c is the multiquadric shape parameter
27+
#In this case we use randomly generated points in the unit ball - d controls the dimension of the data.
28+
#q is a parameter for the algorithm
29+
30+
31+
def FGP_DEMO(n,c,q,d,error):
32+
start_time = time.time()
33+
#set up intial interpolation data - here we use randomly generated points in the unit d-ball.
34+
x_i = rf.points_in_unit_ball(n,d)
35+
f_i = [random.random() for _ in range(n)]
36+
#set up intital values
37+
lambdas = np.zeros(n)
38+
alpha = 0.5*(np.max(f_i)+np.min(f_i))
39+
r = f_i - np.ones(n)*alpha
40+
41+
#randomly shuffling the order of the datapoints to avoid combinations that result in slow convergence.
42+
omeg = list(range(1,n+1))
43+
random.shuffle(omeg)
44+
data = []
45+
46+
#setting up the interpolation matrix
47+
phi = np.ones((n,n))
48+
for i, point in enumerate(x_i):
49+
for j, center in enumerate(x_i):
50+
phi[i,j] = rf.MQ(point-center,c)
51+
52+
#step3 returns the 'lsets' - approximations for the q nearest neighbours for each point (apart from 1)
53+
for m in range(1, n):
54+
newomeg, lset, lvalue = rf.step3(omeg,phi,q,m, n)
55+
omeg = newomeg
56+
data.append([lvalue,lset,rf.step4(lset,x_i,c)])
57+
data = sorted(data, key = lambda x:x[0])
58+
59+
#setup complete, we now look to find the coeffificents for the interpolant to within the prescribed error.
60+
k = 0
61+
err = np.max(np.abs(r))
62+
while err > error:
63+
k+=1
64+
tau = np.zeros(n)
65+
for m in range(len(data)):
66+
taudummy = rf.step5(data[m][1],data[m][2],r)
67+
tau += taudummy
68+
69+
if k == 1:
70+
delta = tau.copy()
71+
72+
else:
73+
beta = np.dot(tau, d)/ np.dot(delta,d)
74+
delta = tau - beta * delta
75+
76+
d = phi @ delta
77+
78+
lambdas, alpha, r, err = rf.step8(lambdas, alpha, r, d, delta)
79+
if k > 500:
80+
break
81+
end_time = time.time()
82+
print(f"Time taken: {end_time - start_time} seconds")
83+
84+
#k is the iteration count
85+
#lambdas are the coefficients of the interpolant we generate
86+
#alpha is the constant in the interpolant we generate
87+
#err is the error in our interpolant.
88+
return k, lambdas, alpha, err

‎FULL_IMPLEMENTATION.py‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
import numpy as np
2+
import random
3+
import radial_functions as rf
4+
import time
5+
6+
#An implementation of the Faul - Goodsen - Powell algorithm. Given data to interpolate,
7+
#the algorithm returns the coefficents (the lambda_i) and the constant term (alpha)
8+
#for a multiquadric interpolant s(x) of the form
9+
10+
#s(x) = sum_i^n lambda_i phi(x-x_i) + alpha
11+
12+
#This interpolant is accurate to within your prescribed error
13+
#The x_i are the data centers, phi(x) is the multiquadric radial basis fuction (x^2+c^2)^0.5 and
14+
#alpha is a constant term.
15+
16+
#INPUT VARIABLES:
17+
#data are the data centers
18+
#values are the data values at the centers
19+
#c is the multiquadric parameter >0. Using a smaller value is generally quicker
20+
#q is a variable controlling the size of sets for the cardinal functions generated as part of the algorithm. Usually we use 5<q<50 - q=30 is a safe value to use. The smaller the value of q, the quicker and dirtier the algorithm will be.
21+
#error is the error you wish to be within.
22+
#We generate our own random data for this implementation.
23+
24+
def FGP(data,values,c,q,error):
25+
26+
n = len(data)
27+
x_i = [i for i in data]
28+
f_i = [i for i in values]
29+
start_time = time.time()
30+
31+
#set up intital values
32+
lambdas = np.zeros(n)
33+
alpha = 0.5*(np.max(f_i)+np.min(f_i))
34+
r = f_i - np.ones(n)*alpha
35+
36+
#randomly shuffling the order of the datapoints to avoid combinations that result in slow convergence.
37+
omeg = list(range(1,n+1))
38+
random.shuffle(omeg)
39+
data = []
40+
41+
#setting up the interpolation matrix
42+
phi = np.ones((n,n))
43+
for i, point in enumerate(x_i):
44+
for j, center in enumerate(x_i):
45+
phi[i,j] = rf.MQ(point-center,c)
46+
47+
#step3 returns the 'lsets' - approximations for the q nearest neighbours for each point (apart from 1)
48+
for m in range(1, n):
49+
newomeg, lset, lvalue = rf.step3(omeg,phi,q,m, n)
50+
omeg = newomeg
51+
data.append([lvalue,lset,rf.step4(lset,x_i,c)])
52+
data = sorted(data, key = lambda x:x[0])
53+
54+
#setup complete, we now look to find the coeffificents for the interpolant to within the prescribed error.
55+
k = 0
56+
err = np.max(np.abs(r))
57+
while err > error:
58+
k+=1
59+
tau = np.zeros(n)
60+
for m in range(len(data)):
61+
taudummy = rf.step5(data[m][1],data[m][2],r)
62+
tau += taudummy
63+
64+
if k == 1:
65+
delta = tau.copy()
66+
67+
else:
68+
beta = np.dot(tau, d)/ np.dot(delta,d)
69+
delta = tau - beta * delta
70+
71+
d = phi @ delta
72+
73+
lambdas, alpha, r, err = rf.step8(lambdas, alpha, r, d, delta)
74+
if k > 500:
75+
break
76+
end_time = time.time()
77+
print(f"Time taken: {end_time - start_time} seconds")
78+
79+
#k is the iteration count
80+
#lambdas are the coefficients of the interpolant we generate
81+
#alpha is the constant in the interpolant we generate
82+
#err is the error in our interpolant.
83+
return k, lambdas, alpha, err

‎radial_functions.py‎

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
1+
import numpy as np
2+
3+
def MQ(x,c):
4+
k = np.sqrt(c**2+np.linalg.norm(x)**2)
5+
return k
6+
7+
#x is the input, l is the vector of values of lambda, a is the value of alpha, X is the vector of x_i, c is the shape parameter
8+
def interp(x, l, a, X, c):
9+
radial_shifts = np.array([MQ(np.abs(x - i),c) for i in X])
10+
l = np.array(l)
11+
interp = np.dot(l, radial_shifts) + a
12+
return interp
13+
14+
def smallest_distance(x):
15+
min_distance = float('inf')
16+
for i in range(len(x)-1):
17+
dist = x[i+1][1]-x[i][1]
18+
if dist < min_distance:
19+
min_distance = dist
20+
Beta = i
21+
Gamma = i+1
22+
return min_distance, Beta, Gamma
23+
24+
#X is the input vector, x is the value with which to evaluate d(x), d is delta, c is shape parameter - this returns d^k(x)
25+
def search_direction(X,x,d,c):
26+
values = [MQ(x-i,c) for i in X]
27+
k = np.dot(values, d)
28+
return k
29+
30+
31+
def step3(omeg, phi, q, m, n):
32+
omeg1 = omeg.copy()
33+
ell = omeg1[m-1]
34+
loopcount = 0
35+
if n - m + 1 > q:
36+
flag = True
37+
38+
while flag:
39+
loopcount+=1
40+
41+
dist2ell = np.zeros(n-m+1)
42+
for j in range(m, n+1):
43+
jj = omeg1[j-1]
44+
dist2ell[j-m] = phi[ell - 1, jj-1]
45+
46+
Sorted_distances = np.argsort(dist2ell)
47+
Lset = [omeg1[i+m-1] for i in Sorted_distances[:q]]
48+
mindist2 = 0.5 * min(dist2ell[1:])
49+
minbeta = 0
50+
mingamma = 1
51+
mindist = np.inf
52+
53+
for beta in range(q):
54+
jbeta = Lset[beta]
55+
56+
for gamma in range(beta+1, q):
57+
jgamma = Lset[gamma]
58+
dist2betagamma = phi[jbeta-1, jgamma-1]
59+
mindistnew = min(mindist, dist2betagamma)
60+
61+
if mindistnew != mindist:
62+
minbeta = jbeta
63+
mingamma = jgamma
64+
mindist = mindistnew
65+
if mindist < mindist2:
66+
dist2gammaell = phi[mingamma - 1, ell -1 ]
67+
dist2betaell = phi[minbeta -1, ell -1]
68+
69+
if dist2gammaell < dist2betaell:
70+
minbeta, mingamma = mingamma, minbeta
71+
72+
mhat = omeg1.index(minbeta)
73+
ell = minbeta
74+
dummy = omeg1[m-1]
75+
omeg1[m-1] = omeg1[mhat]
76+
omeg1[mhat] = dummy
77+
else:
78+
flag = False
79+
else:
80+
Lset = omeg1[m-1:n]
81+
82+
return omeg1, Lset, ell
83+
84+
85+
def step4(lset, x_values,c):
86+
size = len(lset)
87+
x_ell = [x_values[i-1] for i in lset]
88+
z_matrix = np.ones((size+1, size+1))
89+
z_matrix[size][size] = 0
90+
dirac_vector = np.zeros(size+1)
91+
dirac_vector[0] = 1
92+
93+
for i, point in enumerate(x_ell):
94+
for j, center in enumerate(x_ell):
95+
z_matrix[i,j] = MQ(point-center,c)
96+
97+
zeta = np.linalg.solve(z_matrix,dirac_vector)
98+
zeta = zeta[:-1]
99+
zeta[-1] = -sum(zeta[:-1])
100+
return zeta
101+
102+
def step5(lset, zeta, r):
103+
n = len(lset)
104+
taudummy = np.zeros(len(r))
105+
sum = 0
106+
for i in range(n):
107+
sum+= zeta[i]*r[lset[i]-1]
108+
myuell = sum/zeta[0]
109+
110+
for j in range(n):
111+
taudummy[lset[j]-1] = myuell*zeta[j]
112+
113+
return taudummy
114+
115+
116+
def step8(lambdas, alpha, r, d, delta):
117+
gamma = np.dot(delta, r)/np.dot(delta, d)
118+
r1 = r - gamma * d
119+
err = np.max(np.abs(r1))
120+
c = 0.5*(np.max(r1) + np.min(r1))
121+
alpha1 = alpha + c
122+
r2 = r1 - c*np.ones(len(r))
123+
lambdas1 = lambdas + gamma*delta
124+
125+
return lambdas1, alpha1, r2, err
126+
127+
128+
#Different underlying distributions to draw from in the DEMO version.
129+
130+
def points_in_unit_ball(n,d):
131+
cube = np.random.standard_normal(size=(n, d))
132+
norms = np.linalg.norm(cube,axis=1)
133+
surface_sphere = cube/norms[:,np.newaxis]
134+
scales = np.random.uniform(0,1, size= n)
135+
points = surface_sphere * (scales[:, np.newaxis])**(1/d)
136+
return points
137+
138+
def points_in_cube(n,d):
139+
points = np.random.uniform(0,1,size=(n,d))
140+
return points
141+
142+
def points_normal(n,d):
143+
points = np.random.normal(0,1,size=(n,d))
144+
return points
145+
146+
def points_integer_grid(n,d, s_min, s_max):
147+
points = np.random.randint(s_min,s_max, size = (n,d))
148+
return points

0 commit comments

Comments
 (0)