Skip to content

Commit

Permalink
Test the new variants.
Browse files Browse the repository at this point in the history
  • Loading branch information
@kellertuer
kellertuer committed Dec 30, 2024
1 parent 85c1a64 commit 2002ef5
Show file tree
Hide file tree
Showing 3 changed files with 165 additions and 32 deletions.
35 changes: 18 additions & 17 deletions src/plans/nonlinear_least_squares_plan.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,8 +95,11 @@ function get_cost(
) where {E<:AbstractEvaluationType,H<:AbstractManifoldHessianObjective}
v = 0.0
for i in 1:length(nlso.objective)
v += get_cost(vector_space(1), nlso.smoothing, get_value(nlso.objective, p, i)^2)
v += get_cost(
vector_space(1), nlso.smoothing, abs(get_value(M, nlso.objective, p, i))^2
)
end
v *= 0.5
return v
end
function get_cost(
Expand All @@ -108,8 +111,8 @@ function get_cost(
vector_space=Rn,
value_cache=get_value(M, nlso.objective, p),
) where {E<:AbstractEvaluationType,H<:AbstractManifoldHessianObjective}
return sum(
get_cost(vector_space(1), nlso.smoothing, value_cache[i]) for
return 0.5 * sum(
get_cost(vector_space(1), nlso.smoothing, abs(value_cache[i])^2) for
i in 1:length(value_cache)
)
end
Expand All @@ -125,8 +128,11 @@ function get_cost(
) where {E<:AbstractEvaluationType}
v = 0.0
for i in 1:length(nlso.objective)
v += get_cost(vector_space(1), nlso.smoothing, get_value(nlso.objective, p, i)^2, i)
v += get_value(
vector_space(1), nlso.smoothing, abs(get_value(M, nlso.objective, p, i))^2, i
)
end
v *= 0.5
return v
end
function get_cost(
Expand All @@ -138,8 +144,8 @@ function get_cost(
vector_space=Rn,
value_cache=get_value(M, nlso.objective, p),
) where {E<:AbstractEvaluationType}
return sum(
get_cost(vector_space(1), nlso.smoothing, value_cache[i], i) for
return 0.5 * sum(
get_value(vector_space(1), nlso.smoothing, abs(value_cache[i])^2, i) for
i in 1:length(value_cache)
)
end
Expand Down Expand Up @@ -191,6 +197,7 @@ function get_jacobian!(
J,
nlso::NonlinearLeastSquaresObjective{E,AHVF,<:AbstractVectorGradientFunction},
p;
vector_space=Rn,
basis::AbstractBasis=get_basis(nlso.objective.jacobian_type),
value_cache=get_value(M, nlso.objective, p),
) where {E,AHVF}
Expand Down Expand Up @@ -226,7 +233,7 @@ _doc_get_residuals_nlso = """
get_residuals(M::AbstractManifold, nlso::NonlinearLeastSquaresObjective, p)
get_residuals!(M::AbstractManifold, V, nlso::NonlinearLeastSquaresObjective, p)
Compute the vector of residuals ``s_i(f_i(p))``, ``i=1,…,n`` given the manifold `M`,
Compute the vector of residuals ``f_i(p)``, ``i=1,…,n`` given the manifold `M`,
the [`NonlinearLeastSquaresObjective`](@ref) `nlso` and a current point ``p`` on `M`.
# Keyword arguments
Expand Down Expand Up @@ -258,11 +265,10 @@ function get_residuals!(
E,<:AbstractVectorFunction{E,<:ComponentVectorialType},H
},
p;
vector_space=Rn,
kwargs...,
) where {E<:AbstractEvaluationType,H<:AbstractManifoldHessianObjective}
for i in 1:length(nlso.objective)
V[i] = get_cost(vector_space(1), nlso.smoothing, get_value(nlso.objective, p, i)^2)
V[i] = get_value(M, nlso.objective, p, i)
end
return V
end
Expand All @@ -273,11 +279,10 @@ function get_residuals!(
E,<:AbstractVectorFunction{E,<:FunctionVectorialType},H
},
p;
vector_space=Rn,
value_cache=get_value(M, nlso.objective, p),
) where {E<:AbstractEvaluationType,H<:AbstractManifoldHessianObjective}
for i in 1:length(value_cache)
V[i] = get_cost(vector_space(1), nlso.smoothing, value_cache[i])
V[i] = value_cache[i]
end
return V
end
Expand All @@ -289,13 +294,10 @@ function get_residuals!(
E,<:AbstractVectorFunction{E,<:ComponentVectorialType},<:AbstractVectorFunction
},
p;
vector_space=Rn,
kwargs...,
) where {E<:AbstractEvaluationType}
for i in 1:length(nlso.objective)
V[i] = get_cost(
vector_space(1), nlso.smoothing, get_value(nlso.objective, p, i)^2, i
)
V[i] = get_value(M, nlso.objective, p, i)
end
return V
end
Expand All @@ -306,11 +308,10 @@ function get_residuals!(
E,<:AbstractVectorFunction{E,<:FunctionVectorialType},<:AbstractVectorFunction
},
p;
vector_space=Rn,
value_cache=get_value(M, nlso.objective, p),
) where {E<:AbstractEvaluationType}
for i in 1:length(value_cache)
V[i] = get_cost(vector_space(1), nlso.smoothing, value_cache[i], i)
V[i] = value_cache[i]
end
return V
end
Expand Down
29 changes: 17 additions & 12 deletions src/plans/vectorial_plan.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,11 +21,14 @@ for example ``g_i(p) ∈ ℝ^m`` or ``\operatorname{grad} g_i(p) ∈ T_p\mathcal
# Fields
* `basis` an [`AbstractBasis`](@extref `ManifoldsBase.AbstractBasis`) to indicate the default representation.
# Constructor
CoordinateVectorialType(basis=DefaultOrthonormalBasis())
"""
struct CoordinateVectorialType{B<:AbstractBasis} <: AbstractVectorialType
basis::B
end

CoordinateVectorialType() = CoordinateVectorialType(DefaultOrthonormalBasis())

Check warning on line 31 in src/plans/vectorial_plan.jl

View check run for this annotation

Codecov / codecov/patch

src/plans/vectorial_plan.jl#L31 

Added line #L31 was not covered by tests
"""
get_basis(::AbstractVectorialType)
Expand Down Expand Up @@ -641,7 +644,7 @@ function get_jacobian(
n = vgf.range_dimension
d = manifold_dimension(M, p)
mP = PowerManifold(M, range, vgf.range_dimension)
gradients = zero_vector(mP, fill(p, pM))
gradients = zero_vector(mP, fill(p, mP))
vgf.jacobian!!(M, gradients, p)

Check warning on line 648 in src/plans/vectorial_plan.jl

View check run for this annotation

Codecov / codecov/patch

src/plans/vectorial_plan.jl#L644-L648 

Added lines #L644 - L648 were not covered by tests
# generate the first row to get an eltype
c1 = get_coordinates(M, p, gradients[mP, 1], basis)
Expand Down Expand Up @@ -711,10 +714,8 @@ function get_jacobian!(
n = vgf.range_dimension
gradients = vgf.jacobian!!(M, p)
mP = PowerManifold(M, range, vgf.range_dimension)
println("range ", range, "mP", mP)
for i in 1:n
c = @view J[i, :]
println("G ", gradients)
get_coordinates!(M, c, p, gradients[mP, i], basis)
end
return J
Expand Down Expand Up @@ -748,15 +749,21 @@ end
# Part II: mutating vgf – allocating jacobian
# (a) We have a single gradient function
function get_jacobian!(
M::AbstractManifold, JF, vgf::VGF, p; basis::B=DefaultOrthonormalBasis()
M::AbstractManifold,
JF,
vgf::VGF,
p;
basis::B=DefaultOrthonormalBasis(),
range::AbstractPowerRepresentation=get_range(vgf.jacobian_type),
) where {
FT,
VGF<:AbstractVectorGradientFunction{<:InplaceEvaluation,FT,<:FunctionVectorialType},
B<:AbstractBasis,
}
gradients = zero_vector(mP, fill(p, pM))
mP = PowerManifold(M, range, vgf.range_dimension)
gradients = zero_vector(mP, fill(p, mP))
vgf.jacobian!!(M, gradients, p)
for i in 1:n
for i in 1:(vgf.range_dimension)
JF[i, :] .= get_coordinates(M, p, gradients[mP, i], basis)
end
return JF
Expand All @@ -772,7 +779,7 @@ function get_jacobian!(
) where {
FT,VGF<:AbstractVectorGradientFunction{<:InplaceEvaluation,FT,<:ComponentVectorialType}
}
for i in 1:n
for i in 1:(vgf.range_dimension)
vgf.jacobian!![i](M, X, p)
JF[i, :] .= get_coordinates(M, p, X, basis)
end
Expand Down Expand Up @@ -847,8 +854,7 @@ function get_gradient(
)
n = _vgf_index_to_length(i, vgf.range_dimension)
pM = PowerManifold(M, range, n)
P = fill(p, pM)
X = zero_vector(pM, P)
X = zero_vector(pM, fill(p, pM))
return get_gradient!(M, X, vgf, p, i, range)
end
# (c) Special cases where allocations can be skipped
Expand Down Expand Up @@ -994,8 +1000,7 @@ function get_gradient!(
) where {FT<:AbstractVectorialType}
# a type wise safe way to allocate what usually should yield a n-times-d matrix
pM = PowerManifold(M, range, vgf.range_dimension...)
P = fill(p, pM)
Y = zero_vector(pM, P)
Y = zero_vector(pM, fill(p, pM))
JF = reshape(
get_coordinates(pM, fill(p, pM), Y, vgf.jacobian_type.basis),
power_dimensions(pM)...,
Expand Down
133 changes: 130 additions & 3 deletions test/plans/test_nonlinear_least_squares_plan.jl
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,141 @@
using Manifolds, Manopt, Test

@testset "Nonlinear lest squares plane" begin
@testset "Test cost/residual/jacobian cases with smoothing" begin end
@testset "Nonlinear lest squares plan" begin
@testset "Test cost/residual/jacobian cases with smoothing" begin
# a simple nlso objetive on R2
M = Euclidean(2)
d1 = [1, 0]
d2 = [0, 1]
f1(M, x) = norm(x - d1)
f2(M, x) = norm(x - d2)
f(M, x) = [f1(M, x), f2(M, x)]
# Components
f!(M, V, x) = (V .= [f1(M, x), f2(M, x)])
j1(M, x) = (x - d1) / norm(x - d1)
j1!(M, X, x) = (X .= (x - d1) / norm(x - d1))
j2(M, x) = (x - d2) / norm(x - d2)
j2!(M, X, x) = (X .= (x - d2) / norm(x - d2))
# Function
JF(M, x) = [j1(M, x), j2(M, x)]
JF!(M, JF, x) = (JF .= [j1(M, x), j2(M, x)])
# Jacobi matrix
J(M, x) = cat(j1(M, x), j2(M, x); dims=2)
J!(M, J, x) = (J .= cat(j1(M, x), j2(M, x); dims=2))
# Smoothing types
s1 = Manopt.smoothing_factory(:Identity)
s2 = Manopt.smoothing_factory((:Identity, 2))

# Test all (new) possible combinations of vectorial cost and Jacobian
# (1) [F]unction (Gradient), [C]omponent (Gradients), [J] Coordinate (Jacobian in Basis)
# (2) [a]llocating [i] inplace
# (3) [s] single smoothing [v] vector smoothing
nlsoFas = NonlinearLeastSquaresObjective(
f, JF, 2; jacobian_type=FunctionVectorialType(), smoothing=s1
)
nlsoFav = NonlinearLeastSquaresObjective(
f, JF, 2; jacobian_type=FunctionVectorialType(), smoothing=s2
)
nlsoFis = NonlinearLeastSquaresObjective(
f!,
JF!,
2;
evaluation=InplaceEvaluation(),
jacobian_type=FunctionVectorialType(),
smoothing=s1,
)
nlsoFiv = NonlinearLeastSquaresObjective(
f!,
JF!,
2;
evaluation=InplaceEvaluation(),
jacobian_type=FunctionVectorialType(),
smoothing=s2,
)

nlsoCas = NonlinearLeastSquaresObjective(
[f1, f2],
[j1, j2],
2;
function_type=ComponentVectorialType(),
jacobian_type=ComponentVectorialType(),
smoothing=s1,
)
nlsoCav = NonlinearLeastSquaresObjective(
[f1, f2],
[j1, j2],
2;
function_type=ComponentVectorialType(),
jacobian_type=ComponentVectorialType(),
smoothing=s2,
)
nlsoCis = NonlinearLeastSquaresObjective(
[f1, f2],
[j1!, j2!],
2;
function_type=ComponentVectorialType(),
jacobian_type=ComponentVectorialType(),
evaluation=InplaceEvaluation(),
smoothing=s1,
)
nlsoCiv = NonlinearLeastSquaresObjective(
[f1, f2],
[j1!, j2!],
2;
function_type=ComponentVectorialType(),
jacobian_type=ComponentVectorialType(),
evaluation=InplaceEvaluation(),
smoothing=s2,
)

nlsoJas = NonlinearLeastSquaresObjective(f, J, 2; smoothing=s1)
nlsoJav = NonlinearLeastSquaresObjective(f, J, 2; smoothing=s2)
nlsoJis = NonlinearLeastSquaresObjective(
f!, J!, 2; evaluation=InplaceEvaluation(), smoothing=s1
)
nlsoJiv = NonlinearLeastSquaresObjective(
f!, J!, 2; evaluation=InplaceEvaluation(), smoothing=s2
)

p = [0.5, 0.5]
V = [0.0, 0.0]
Vt = [1 / sqrt(2), 1 / sqrt(2)]
G = zeros(2, 2)
Gt = 1 / sqrt(2) .* [-1.0 1.0; 1.0 -1.0]
for nlso in [
nlsoFas,
nlsoFav,
nlsoFis,
nlsoFiv,
nlsoCas,
nlsoCav,
nlsoCis,
nlsoCiv,
nlsoJas,
nlsoJav,
nlsoJis,
nlsoJiv,
]
c = get_cost(M, nlso, p)
@test c ≈ 0.5
fill!(V, 0.0)
get_residuals!(M, V, nlso, p)
@test V == get_residuals(M, nlso, p)
@test V ≈ Vt
@test 0.5 * sum(abs.(V) .^ 2) ≈ c
fill!(G, 0.0)
get_jacobian!(M, G, nlso, p)
@test G == get_jacobian(M, nlso, p)
@test G == Gt
end
end
@testset "Smootthing factory" begin
s1 = Manopt.smoothing_factory(:Identity)
@test s1 isa ManifoldHessianObjective

@test Manopt.smoothing_factory(s1) === s1 # Passthrough for mhos
s2 = Manopt.smoothing_factory((:Identity, 2))
@test s2 isa VectorHessianFunction
@test length(s2) == 2
@test Manopt.smoothing_factory(s2) === s2 # Passthrough for vhfs

s3 = Manopt.smoothing_factory((:Identity, 3.0))
@test s3 isa ManifoldHessianObjective
Expand Down

0 comments on commit 2002ef5

Please sign in to comment.