Skip to content

Commit

Permalink
Revert test and deprecate analytical jacobian for BezierCurve
Browse files Browse the repository at this point in the history
  • Loading branch information
mikeingold committed Nov 13, 2024
1 parent 0f87431 commit 8c90f3f
Show file tree
Hide file tree
Showing 8 changed files with 18 additions and 65 deletions.
43 changes: 1 addition & 42 deletions src/specializations/BezierCurve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ function integral(
kwargs...
)
paramfunction(t) = _parametric(curve, t; alg = alg)
param_curve = _ParametricGeometry(paramfunction, curve, 1)
param_curve = _ParametricGeometry(paramfunction, 1)
return _integral(f, param_curve, rule; kwargs...)
end

Expand All @@ -50,44 +50,3 @@ end
function _parametric(curve::Meshes.BezierCurve, t; alg::Meshes.BezierEvalMethod)
return curve(t, alg)
end

################################################################################
# jacobian
################################################################################

function jacobian(
curve::_ParametricGeometry{M, C, Meshes.BezierCurve, F, Dim},
ts::Union{AbstractVector{T}, Tuple{T, Vararg{T}}},
::Analytical
) where {M, C, F, Dim, T <: AbstractFloat}
t = only(ts)
# Parameter t restricted to domain [0,1] by definition
if t < 0 || t > 1
throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
end

# Aliases
P = curve.source.controls
N = Meshes.degree(curve.source)

# Ensure that this implementation is tractible: limited by ability to calculate
# binomial(N, N/2) without overflow. It's possible to extend this range by
# converting N to a BigInt, but this results in always returning BigFloat's.
N <= 1028 || error("This algorithm overflows for curves with ⪆1000 control points.")

# Generator for Bernstein polynomial functions
B(i, n) = t -> binomial(Int128(n), i) * t^i * (1 - t)^(n - i)

# Derivative = N Σ_{i=0}^{N-1} sigma(i)
# P indices adjusted for Julia 1-based array indexing
sigma(i) = B(i, N - 1)(t) .* (P[(i + 1) + 1] - P[(i) + 1])
derivative = N .* sum(sigma, 0:(N - 1))

return (derivative,)
end

function _has_analytical(
::Type{_ParametricGeometry{M, C, Meshes.BezierCurve, F, Dim}}
) where {M, C, F, Dim}
return true
end
2 changes: 1 addition & 1 deletion src/specializations/Line.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ function integral(
kwargs...
)
paramfunction(t) = _parametric(line, t)
param_line = _ParametricGeometry(paramfunction, line, 1)
param_line = _ParametricGeometry(paramfunction, 1)
return _integral(f, param_line, rule; kwargs...)
end

Expand Down
2 changes: 1 addition & 1 deletion src/specializations/Plane.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ function integral(
kwargs...
)
paramfunction(t1, t2) = _parametric(plane, t1, t2)
param_plane = _ParametricGeometry(paramfunction, plane, 2)
param_plane = _ParametricGeometry(paramfunction, 2)
return _integral(f, param_plane, rule; kwargs...)
end

Expand Down
2 changes: 1 addition & 1 deletion src/specializations/Ray.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ function integral(
kwargs...
)
paramfunction(t) = _parametric(ray, t)
param_ray = _ParametricGeometry(paramfunction, ray, 1)
param_ray = _ParametricGeometry(paramfunction, 1)
return _integral(f, param_ray, rule; kwargs...)
end

Expand Down
2 changes: 1 addition & 1 deletion src/specializations/Tetrahedron.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ function integral(
kwargs...
) where {F <: Function}
paramfunction(t1, t2, t3) = _parametric(tetrahedron, t1, t2, t3)
tetra = _ParametricGeometry(paramfunction, tetrahedron, 3)
tetra = _ParametricGeometry(paramfunction, 3)
return _integral(f, tetra, rule; kwargs...)
end

Expand Down
2 changes: 1 addition & 1 deletion src/specializations/Triangle.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ function integral(
kwargs...
) where {F <: Function}
paramfunction(t1, t2) = _parametric(triangle, t1, t2)
tri = _ParametricGeometry(paramfunction, triangle, 2)
tri = _ParametricGeometry(paramfunction, 2)
return _integral(f, tri, rule; kwargs...)
end

Expand Down
22 changes: 8 additions & 14 deletions src/specializations/_ParametricGeometry.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,53 +12,47 @@ no longer be required.
# Fields
- `fun::Function` - a parametric function: (ts...) -> Meshes.Point
- `source::Meshes.Geometry` - the source geometry being transformed
# Type Structure
- `M <: Meshes.Manifold` - same usage as in `Meshes.Geometry{M, C}`
- `C <: CoordRefSystems.CRS` - same usage as in `Meshes.Geometry{M, C}`
- `G <: Meshes.Geometry` - the `Meshes.Geometry` type being transformed
- `F` - type of the callable integrand function
- `Dim` - number of parametric dimensions
"""
struct _ParametricGeometry{M <: Meshes.Manifold, C <: CRS, G <: Geometry, F, Dim} <:
struct _ParametricGeometry{M <: Meshes.Manifold, C <: CRS, F, Dim} <:
Meshes.Primitive{M, C}
fun::F
source::G

function _ParametricGeometry{M, C}(
fun::F,
source::G,
Dim::Int64
) where {M <: Meshes.Manifold, C <: CRS, G <: Geometry, F}
return new{M, C, G, F, Dim}(fun, source)
) where {M <: Meshes.Manifold, C <: CRS, F}
return new{M, C, F, Dim}(fun)
end
end

"""
_ParametricGeometry(fun, type, dims)
_ParametricGeometry(fun, dims)
Construct a `_ParametricGeometry` using a provided parametric function `fun` for a geometry
with `dims` parametric dimensions.
# Arguments
- `fun::Function` - parametric function mapping `(ts...) -> Meshes.Point`
- `type::Type{G} where {G <: Geometry}` - the `Meshes.Geometry` type being transformed
- `dims::Int64` - number of parametric dimensions, i.e. `length(ts)`
"""
function _ParametricGeometry(
fun::F,
source::G,
dims::Int64
) where {F <: Function, G <: Geometry}
Dim::Int64
) where {F <: Function}
p = fun(_zeros(dims)...)
return _ParametricGeometry{Meshes.manifold(p), Meshes.crs(p)}(fun, source, dims)
return _ParametricGeometry{Meshes.manifold(p), Meshes.crs(p)}(fun, Dim)
end

# Allow a _ParametricGeometry to be called like a Geometry
(g::_ParametricGeometry)(ts...) = g.fun(ts...)

Meshes.paramdim(::_ParametricGeometry{M, C, F, G, Dim}) where {M, C, F, G, Dim} = Dim
Meshes.paramdim(::_ParametricGeometry{M, C, F, Dim}) where {M, C, F, Dim} = Dim

"""
_parametric(geometry::G, ts...) where {G <: Meshes.Geometry}
Expand Down
8 changes: 4 additions & 4 deletions test/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,17 +24,17 @@ end

# _has_analytical of instances
bezier = BezierCurve([Point(t, sin(t), 0.0) for t in range(-π, π, length = 361)])
@test _has_analytical(bezier) == true
@test _has_analytical(bezier) == false
sphere = Sphere(Point(0, 0, 0), 1.0)
@test _has_analytical(sphere) == false

# _has_analytical of types
@test _has_analytical(Meshes.BezierCurve) == true
@test _has_analytical(Meshes.BezierCurve) == false
@test _has_analytical(Meshes.Sphere) == false

# _default_method
@test _default_method(Meshes.BezierCurve) isa Analytical
@test _default_method(bezier) isa Analytical
@test _default_method(Meshes.BezierCurve) isa FiniteDifference
@test _default_method(bezier) isa FiniteDifference
@test _default_method(Meshes.Sphere) isa FiniteDifference
@test _default_method(sphere) isa FiniteDifference

Expand Down

0 comments on commit 8c90f3f

Please sign in to comment.