-
Notifications
You must be signed in to change notification settings - Fork 105
/
constant.jl
102 lines (88 loc) · 3.26 KB
/
constant.jl
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
export Nearest, Previous, Next
abstract type ConstantInterpType end
"""
Default parameter for `Constant` that performs *nearest-neighbor* interpolation.
Can optionally be specified as
```
Constant() === Constant{Nearest}()
```
"""
struct Nearest <: ConstantInterpType end
"""
Parameter for `Constant` that performs *previous-neighbor* interpolations.
Applied through
```
Constant{Previous}()
```
"""
struct Previous <: ConstantInterpType end
"""
Parameter for `Constant` that performs *next-neighbor* interpolations.
Applied through
```
Constant{Next}()
```
"""
struct Next <: ConstantInterpType end
struct Constant{T<:ConstantInterpType,BC<:Union{Throw{OnGrid},Periodic{OnCell}}} <: DegreeBC{0}
bc::BC
function Constant{T, BC}(bc::BC=BC()) where {T<:ConstantInterpType, BC<:Union{Throw{OnGrid},Periodic{OnCell}}}
new{T, BC}(bc)
end
end
# Default to Nearest and Throw{OnGrid}
Constant() = Constant{Nearest}()
Constant(bc::BC) where {BC <: BoundaryCondition} = Constant{Nearest}(bc)
Constant(::Type{T}) where T <: ConstantInterpType = Constant{T}()
Constant{T}() where {T<:ConstantInterpType} = Constant{T,Throw{OnGrid}}(Throw(OnGrid()))
Constant{T}(bc::BC) where {T<:ConstantInterpType,BC<:BoundaryCondition} = Constant{T,BC}(bc)
Constant{T}(::Periodic{Nothing}) where {T<:ConstantInterpType} = Constant{T,Periodic{OnCell}}(Periodic(OnCell()))
function Base.show(io::IO, deg::Constant)
print(io, nameof(typeof(deg)), '{', typeof(deg).parameters[1], '}', '(')
show(io, deg.bc)
print(io, ')')
end
function Base.show(io::IO, deg::Constant{T,Throw{OnGrid}}) where {T <: ConstantInterpType}
print(io, nameof(typeof(deg)), '{', typeof(deg).parameters[1], '}', '(', ')')
end
"""
Constant b-splines are *nearest-neighbor* interpolations, and effectively
return `A[round(Int,x)]` when interpolating without scaling.
Scaling can lead to inaccurate position of the *neighbors* due to limited numerical precision.
`Constant{Previous}` interpolates to the previous value and is thus equivalent
to `A[floor(Int,x)]` without scaling.
`Constant{Next}` interpolates to the next value and is thus equivalent
to `A[ceil(Int,x)]` without scaling.
"""
Constant
function positions(c::Constant{Previous}, ax, x) # discontinuity occurs at integer locations
xm = floorbounds(x, ax)
δx = x - xm
fast_trunc(Int, xm), δx
end
function positions(c::Constant{Next}, ax, x) # discontinuity occurs at integer locations
xm = ceilbounds(x, ax)
δx = x - xm
fast_trunc(Int, xm), δx
end
function positions(c::Constant{Nearest}, ax, x) # discontinuity occurs at half-integer locations
xm = roundbounds(x, ax)
δx = x - xm
fast_trunc(Int, xm), δx
end
function positions(c::Constant{Previous,Periodic{OnCell}}, ax, x)
# We do not use floorbounds because we do not want to add a half at
# the lowerbound to round up.
xm = floor(x)
δx = x - xm
modrange(fast_trunc(Int, xm), ax), δx
end
function positions(c::Constant{Next,Periodic{OnCell}}, ax, x) # discontinuity occurs at integer locations
xm = ceilbounds(x, ax)
δx = x - xm
modrange(fast_trunc(Int, xm), ax), δx
end
value_weights(::Constant, δx) = (1,)
gradient_weights(::Constant, δx) = (0,)
hessian_weights(::Constant, δx) = (0,)
padded_axis(ax::AbstractUnitRange, ::BSpline{<:Constant}) = ax