/
utils.jl
156 lines (134 loc) · 6.78 KB
/
utils.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
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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
"""
adapted_grid(f, anchors;
max_recursions::Real = 7, max_change = 5°, rand_factor::Real = 0.05,
grid_step = 1.0μm)
Computes a resampled `grid` given anchor points so that `f.(grid)` is sufficiently
smooth. The method used is to create an initial grid around the anchor points and refine
intervals. When an interval becomes "straight enough" it is no longer divided.
Adapted from a contribution to PlotUtils.jl from Kristoffer Carlsson.
- `max_recursions`: how many times each interval is allowed to be refined.
- `max_change`: specifies acceptable change between evaluations of `f` on subsequent grid
points, as estimated by the derivative times the distance between grid points.
Typically, `f` is the angle of a path in the plane, so this is often an angle threshold.
This condition is approximately valid in the end result, but may be weakly violated.
This condition may be grossly violated if `max_recursions` is too low.
- `rand_factor`: between anchor points, `adapted_grid` will wiggle initial grid points
a bit to prevent aliasing. The wiggling is sampled uniformly from the interval:
`[-rand_factor, rand_factor]`, times the distance between three grid points
(e.g. `i+1` and `i-1`). A random number generator is given a fixed seed every time
`adapted_grid` is called, so the rendered results are deterministic.
- `grid_step`: Step size for initial grid points. If you set this to be larger than
the maximum anchor point, then the lowest resolution consistent with `max_change` is used
(unless `f` has some fast variations that the algorithm might miss).
"""
function adapted_grid end
adapted_grid(f, anchors::Tuple{S,T}; kwargs...) where {S <: Coordinate, T <: Coordinate} =
adapted_grid(f, StaticArrays.SVector(anchors); kwargs...)
adapted_grid(f, anchors::AbstractVector{T};
max_recursions::Real = 7,
max_change = 5°,
rand_factor::Real = 0.05,
grid_step::Coordinate = 1.0 * Unitful.ContextUnits(
ifelse(T<:Length, Unitful.μm, Unitful.NoUnits),
Unitful.unit(Unitful.upreferred(zero(T))))) where {T <: Coordinate} =
assemble_grids(f, anchors, max_recursions, max_change, rand_factor, grid_step)
function assemble_grids(f, anchors::AbstractVector,
max_recursions, max_change, rand_factor, grid_step)
dimension(f(anchors[1])) != dimension(max_change) &&
throw(ArgumentError("max_change must have dimensions of f($(anchors[1]))"))
g = i->begin
# Want a range that specially excludes the last point; use linspace with npts.
@inbounds npts = Int(ceil(NoUnits((anchors[i+1] - anchors[i]) / grid_step))) - 1
npts = ifelse(npts < 5, 5, ifelse(iseven(npts), npts+1, npts))
@inbounds grid = make_grid(f, range(anchors[i], stop = anchors[i+1], length = npts),
max_recursions, max_change, rand_factor)
grid[1:(end-1)]
end
grid = vcat(g.(1:(length(anchors)-1))...)
push!(grid, anchors[end])
end
function make_grid(f, initial_grid, max_recursions, max_change, rand_factor)
# Initial assertions and checks
@assert max_recursions >= 0
@assert rand_factor >= 0
n_points = length(initial_grid)
@assert isodd(n_points) && n_points >= 3
!issorted(initial_grid) && throw(ArgumentError("initial grid must be sorted."))
# Initial number of points
n_intervals = div(n_points, 2)
xs = collect(initial_grid)
if rand_factor != 0
rng = MersenneTwister(1337)
for i in 2:length(xs)-1
xs[i] += rand_factor * 2 * (rand(rng) - 0.5) * (xs[i+1] - xs[i-1])
end
end
n_tot_refinements = zeros(Int, n_intervals)
while true
results = zeros(typeof(f(initial_grid[1])), n_intervals)
active = Vector{Bool}(undef, n_intervals)
# derivs = [ForwardDiff.derivative(f, xs[i]) for i in 1:(2*n_intervals + 1)]
for interval in 1:n_intervals
p = 2 * interval
tot_w = 0.0
# Do a small convolution
for (q,w) in ((-1, 0.25), (0, 0.5), (1, 0.25))
interval == 1 && q == -1 && continue
interval == n_intervals && q == 1 && continue
tot_w += w
i = p + q
results[interval] += abs(ForwardDiff.derivative(f, xs[i]) * (xs[i+1]-xs[i-1])/2) * w
end
results[interval] /= tot_w
# Only consider intervals that have not been refined too much and have a high enough curvature
active[interval] = n_tot_refinements[interval] < max_recursions && results[interval] > max_change
end
if all(x -> x >= max_recursions, n_tot_refinements[active])
break
end
n_target_refinements = div(n_intervals, 2)
interval_candidates = collect(1:n_intervals)[active]
n_refinements = min(n_target_refinements, length(interval_candidates))
perm = sortperm(results[active])
intervals_to_refine = sort(interval_candidates[perm[length(perm) - n_refinements + 1:end]])
n_intervals_to_refine = length(intervals_to_refine)
n_new_points = 2*length(intervals_to_refine)
# Do division of the intervals
new_xs = similar(xs, n_points + n_new_points)
new_tot_refinements = similar(n_tot_refinements, n_intervals + n_intervals_to_refine)
k = 0
kk = 0
for i in 1:n_points
if iseven(i) # This is a point in an interval
interval = div(i, 2)
if interval in intervals_to_refine
# The refined interval will become two intervals in the next round.
# Both new intervals take on the old intervals refinement number + 1.
kk += 1
new_tot_refinements[interval - 1 + kk] = n_tot_refinements[interval] + 1
new_tot_refinements[interval + kk] = n_tot_refinements[interval] + 1
# Insert new x taken between i and i-1 of interval.
k += 1
new_xs[i-1 + k] = (xs[i] + xs[i-1]) / 2
# Insert old x
new_xs[i + k] = xs[i]
# Insert new x taken between i and i+1 of interval.
new_xs[i+1 + k] = (xs[i+1] + xs[i]) / 2
# Increment k again since we added two new points
k += 1
else
# Just keep the old x and refinement count to maintain the interval
new_tot_refinements[interval + kk] = n_tot_refinements[interval]
new_xs[i + k] = xs[i]
end
else
new_xs[i + k] = xs[i]
end
end
xs = new_xs
n_tot_refinements = new_tot_refinements
n_points = n_points + n_new_points
n_intervals = div(n_points, 2)
end
return xs
end