Ian Jauslin
summaryrefslogtreecommitdiff
blob: 066bc20acc723a870ea640b730bae7da827f8e74 (plain)
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