1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
## Copyright 2021-2023 Ian Jauslin
##
## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License.
## You may obtain a copy of the License at
##
## http://www.apache.org/licenses/LICENSE-2.0
##
## Unless required by applicable law or agreed to in writing, software
## distributed under the License is distributed on an "AS IS" BASIS,
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
## See the License for the specific language governing permissions and
## limitations under the License.
# linear interpolation: given vectors x,y, compute a linear interpolation for y(x0)
# assume x is ordered
@everywhere function linear_interpolation(
x0::Float64,
x::Array{Float64,1},
y::Array{Float64,1}
)
# if x0 is beyond all x's, then return the corresponding boundary value.
if x0>x[length(x)]
return y[length(y)]
elseif x0<x[1]
return y[1]
end
# find bracketing interval
i=bracket(x0,x,1,length(x))
# interpolate
return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i])
end
@everywhere function bracket(
x0::Float64,
x::Array{Float64,1},
a::Int64,
b::Int64
)
i=floor(Int64,(a+b)/2)
if x0<x[i]
return bracket(x0,x,a,i)
elseif x0>x[i+1]
return bracket(x0,x,i+1,b)
else
return i
end
end
# polynomial interpolation of a family of points
@everywhere function poly_interpolation(
x::Array{Float64,1},
y::Array{Float64,1}
)
# init for recursion
rec=Array{Polynomial{Float64},1}(undef,length(x))
for i in 1:length(x)
rec[i]=Polynomial([1.])
end
# compute \prod (x-x_i)
poly_interpolation_rec(rec,x,1.,length(x))
# sum terms together
out=0.
for i in 1:length(y)
out+=rec[i]/rec[i](x[i])*y[i]
end
return out
end
# recursive helper function
@everywhere function poly_interpolation_rec(
out::Array{Float64,1},
x::Array{Float64,1},
a::Float64,
b::Float64
)
if a==b
return
end
# mid point
mid=floor(Int64,(a+b)/2)
# multiply left and right of mid
right=Polynomial([1.])
for i in mid+1:b
right*=Polynomial([-x[i],1.])
end
left=Polynomial([1.])
for i in a:mid
left*=Polynomial([-x[i],1.])
end
# multiply into left and right
for i in a:mid
out[i]*=right
end
for i in mid+1:b
out[i]*=left
end
# recurse
poly_interpolation_rec(out,x,a,mid)
poly_interpolation_rec(out,x,mid+1,b)
return
end
## the following does the same, but has complexity N^2, the function above has N*log(N)
#@everywhere function poly_interpolation(
# x::Array{Float64,1},
# y::Array{Float64,1}
#)
# out=Polynomial([0.])
# for i in 1:length(x)
# prod=Polynomial([1.])
# for j in 1:length(x)
# if j!=i
# prod*=Polynomial([-x[j]/(x[i]-x[j]),1/(x[i]-x[j])])
# end
# end
# out+=prod*y[i]
# end
#
# return out
#end
|
