Fix a few easier bugs and typos.

src/plans/augmented_lagrangian_plan.jl

@@ -44,8 +44,8 @@ function set_manopt_parameter!(alc::AugmentedLagrangianCost, ::Val{:λ}, λ)
     return alc
 end
 function (L::AugmentedLagrangianCost)(M::AbstractManifold, p)
-    gp = get_inequality_constraints(M, L.co, p)
-    hp = get_equality_constraints(M, L.co, p)
+    gp = get_inequality_constraint(M, L.co, p, :)
+    hp = get_equality_constraint(M, L.co, p, :)
     m = length(gp)
     n = length(hp)
     c = get_cost(M, L.co, p)
@@ -124,7 +124,9 @@ function (LG::AugmentedLagrangianGrad)(
         end
     end
     if n > 0
-        X .+= sum((hp .* LG.ρ .+ LG.λ) .* get_grad_equality_constraints(M, LG.co, p, range))
+        X .+= sum(
+            (hp .* LG.ρ .+ LG.λ) .* get_grad_equality_constraint(M, LG.co, p, :, range)
+        )
     end
     return X
 end

src/plans/count.jl

@@ -307,7 +307,7 @@ end
 # Constraint
 function get_constraints(M::AbstractManifold, co::ManifoldCountObjective, p)
     _count_if_exists(co, :Constraints)
-    return get_constraint(M, co.objective, p)
+    return get_constraints(M, co.objective, p)
 end
 function get_equality_constraint(
     M::AbstractManifold, co::ManifoldCountObjective, p, c::Colon
@@ -340,7 +340,7 @@ function get_inequality_constraint(
     M::AbstractManifold, co::ManifoldCountObjective, p, i::Colon
 )
     _count_if_exists(co, :InequalityConstraints)
-    return get_inequality_constraints(M, co.objective, p)
+    return get_inequality_constraint(M, co.objective, p, i)
 end
 function get_inequality_constraint(
     M::AbstractManifold, co::ManifoldCountObjective, p, i::Integer
@@ -390,7 +390,7 @@ function get_grad_inequality_constraint(
     M::AbstractManifold, co::ManifoldCountObjective, p, i::Colon
 )
     _count_if_exists(co, :GradInequalityConstraints)
-    return get_grad_inequality_constraints(M, co.objective, p, i)
+    return get_grad_inequality_constraint(M, co.objective, p, i)
 end
 function get_grad_inequality_constraint!(
     M::AbstractManifold, X, co::ManifoldCountObjective, p, i::Colon

src/plans/exact_penalty_method_plan.jl

@@ -71,17 +71,17 @@ function set_manopt_parameter!(epc::ExactPenaltyCost, ::Val{:u}, u)
     return epc
 end
 function (L::ExactPenaltyCost{<:LogarithmicSumOfExponentials})(M::AbstractManifold, p)
-    gp = get_inequality_constraints(M, L.co, p)
-    hp = get_equality_constraints(M, L.co, p)
+    gp = get_inequality_constraint(M, L.co, p, :)
+    hp = get_equality_constraint(M, L.co, p, :)
     m = length(gp)
     n = length(hp)
     cost_ineq = (m > 0) ? sum(L.u .* log.(1 .+ exp.(gp ./ L.u))) : 0.0
     cost_eq = (n > 0) ? sum(L.u .* log.(exp.(hp ./ L.u) .+ exp.(-hp ./ L.u))) : 0.0
     return get_cost(M, L.co, p) + (L.ρ) * (cost_ineq + cost_eq)
 end
 function (L::ExactPenaltyCost{<:LinearQuadraticHuber})(M::AbstractManifold, p)
-    gp = get_inequality_constraints(M, L.co, p)
-    hp = get_equality_constraints(M, L.co, p)
+    gp = get_inequality_constraint(M, L.co, p, :)
+    hp = get_equality_constraint(M, L.co, p, :)
     m = length(gp)
     n = length(hp)
     cost_eq_greater_u = (m > 0) ? sum((gp .- L.u / 2) .* (gp .> L.u)) : 0.0
@@ -135,107 +135,45 @@ function (EG::ExactPenaltyGrad)(M::AbstractManifold, p)
     return EG(M, X, p)
 end
 function (EG::ExactPenaltyGrad{<:LogarithmicSumOfExponentials})(M::AbstractManifold, X, p)
-    gp = get_inequality_constraints(M, EG.co, p)
-    hp = get_equality_constraints(M, EG.co, p)
+    gp = get_inequality_constraint(M, EG.co, p, :)
+    hp = get_equality_constraint(M, EG.co, p, :)
     m = length(gp)
     n = length(hp)
     # start with `gradf`
     get_gradient!(M, X, EG.co, p)
     c = 0
     # add gradient of the components of g
     (m > 0) && (c = EG.ρ .* exp.(gp ./ EG.u) ./ (1 .+ exp.(gp ./ EG.u)))
-    (m > 0) && (X .+= sum(get_grad_inequality_constraints(M, EG.co, p) .* c))
+    (m > 0) && (X .+= sum(get_grad_inequality_constraint(M, EG.co, p, :) .* c))
     # add gradient of the components of h
+    # TODO: This could be reduced and only active gradients be evaluated
     (n > 0) && (
         c =
             EG.ρ .* (exp.(hp ./ EG.u) .- exp.(-hp ./ EG.u)) ./
             (exp.(hp ./ EG.u) .+ exp.(-hp ./ EG.u))
     )
-    (n > 0) && (X .+= sum(get_grad_equality_constraints(M, EG.co, p) .* c))
+    (n > 0) && (X .+= sum(get_grad_equality_constraint(M, EG.co, p, :) .* c))
     return X
 end
 
 # Default (functions constraints): evaluate all gradients
 function (EG::ExactPenaltyGrad{<:LinearQuadraticHuber})(
     M::AbstractManifold, X, p::P
 ) where {P}
-    gp = get_inequality_constraints(M, EG.co, p)
-    hp = get_equality_constraints(M, EG.co, p)
+    gp = get_inequality_constraint(M, EG.co, p, :)
+    hp = get_equality_constraint(M, EG.co, p, :)
     m = length(gp)
     n = length(hp)
     get_gradient!(M, X, EG.co, p)
     if m > 0
-        gradgp = get_grad_inequality_constraints(M, EG.co, p)
+        gradgp = get_grad_inequality_constraint(M, EG.co, p, :)
         X .+= sum(gradgp .* (gp .>= EG.u) .* EG.ρ) # add the ones >= u
         X .+= sum(gradgp .* (gp ./ EG.u .* (0 .<= gp .< EG.u)) .* EG.ρ) # add < u
     end
     if n > 0
+        # TODO: This could be reduced and only active gradients be evaluated
         c = (hp ./ sqrt.(hp .^ 2 .+ EG.u^2)) .* EG.ρ
-        X .+= sum(get_grad_equality_constraints(M, EG.co, p) .* c)
+        X .+= sum(get_grad_equality_constraint(M, EG.co, p, :) .* c)
     end
     return X
 end
-# Variant 2: vectors of allocating gradients
-#= TODO: rework to easier access with substructures
-function (
-    EG::ExactPenaltyGrad{
-        <:LinearQuadraticHuber,
-        <:ConstrainedManifoldObjective{AllocatingEvaluation,<:VectorConstraint},
-    }
-)(
-    M::AbstractManifold, X, p::P
-) where {P}
-    m = length(EG.co.g)
-    n = length(EG.co.h)
-    get_gradient!(M, X, EG.co, p)
-    for i in 1:m
-        gpi = get_inequality_constraint(M, EG.co, p, i)
-        if (gpi >= 0) # these are the only necessary allocations.
-            (gpi .>= EG.u) && (X .+= gpi .* EG.ρ)
-            (0 < gpi < EG.u) && (X .+= gpi .* (gpi / EG.u) * EG.ρ)
-        end
-    end
-    for j in 1:n
-        hpj = get_equality_constraint(M, EG.co, p, j)
-        if hpj > 0
-            c = hpj / sqrt(hpj^2 + EG.u^2)
-            X .+= get_grad_equality_constraint(M, EG.co, p, j) .* (c * EG.ρ)
-        end
-    end
-    return X
-end
-
-# Variant 3: vectors of mutating gradients
-function (
-    EG::ExactPenaltyGrad{
-        <:LinearQuadraticHuber,
-        <:ConstrainedManifoldObjective{InplaceEvaluation,<:VectorConstraint},
-    }
-)(
-    M::AbstractManifold, X, p::P
-) where {P}
-    m = length(EG.co.g)
-    n = length(EG.co.h)
-    get_gradient!(M, X, EG.co, p)
-    Y = zero_vector(M, p)
-    for i in 1:m
-        gpi = get_inequality_constraint(M, EG.co, p, i)
-        if (gpi >= 0) # the cases where to evaluate the gradient
-            # only evaluate the gradient if `gpi > 0`
-            get_grad_inequality_constraint!(M, Y, EG.co, p, i)
-            # just add the gradient scaled by ρ
-            (gpi >= EG.u) && (X .+= EG.ρ .* Y)
-            # use a different factor, but exclude the case `g = 0` as well
-            (0 < gpi < EG.u) && (X .+= ((gpi / EG.u) * EG.ρ) .* Y)
-        end
-    end
-    for j in 1:n
-        hpj = get_equality_constraint(M, EG.co, p, j)
-        if hpj > 0
-            get_grad_equality_constraint!(M, Y, EG.co, p, j)
-            X .+= ((hpj / sqrt(hpj^2 + EG.u^2)) * EG.ρ) .* Y
-        end
-    end
-    return X
-end
-=#

src/plans/vectorial_plan.jl

@@ -153,11 +153,11 @@ Since `i` is assumed to be a linear index, you can provide
 get_cost(M::AbstractManifold, vgf::VectorGradientFunction, p, i, range=nothing)
 function get_cost(
     M::AbstractManifold,
-    vgf::VectorGradientFunction{<:AllocatingEvaluation,<:FunctionVectorialType},
+    vgf::VectorGradientFunction{E,<:FunctionVectorialType},
     p,
     i,
     range=nothing,
-)
+) where {E}
     c = vgf.costs!!(M, p)
     if isa(c, Number)
         return c
@@ -167,20 +167,20 @@ function get_cost(
 end
 function get_cost(
     M::AbstractManifold,
-    vgf::VectorGradientFunction{<:AllocatingEvaluation,<:ComponentVectorialType},
+    vgf::VectorGradientFunction{E,<:ComponentVectorialType},
     p,
     i::Integer,
     range=nothing,
-)
+) where {E}
     return vgf.costs!![i](M, p)
 end
 function get_cost(
     M::AbstractManifold,
-    vgf::VectorGradientFunction{<:AllocatingEvaluation,<:ComponentVectorialType},
+    vgf::VectorGradientFunction{E,<:ComponentVectorialType},
     p,
     i,
     range=nothing,
-)
+) where {E}
     return [f(M, p) for f in vgf.costs!![i]]
 end
 
@@ -351,6 +351,18 @@ function get_gradient!(
     copyto!(mP, X, vgf.jacobian!!(M, p)[mP, i])
     return X
 end
+function get_gradient!(
+    M::AbstractManifold,
+    X,
+    vgf::VectorGradientFunction{<:AllocatingEvaluation,FT,<:FunctionVectorialType},
+    p,
+    i::Integer,
+    range::Union{AbstractPowerRepresentation,Nothing}=NestedPowerRepresentation(),
+) where {FT}
+    mP = PowerManifold(M, range, vgf.range_dimension)
+    copyto!(mP, X, vgf.jacobian!!(M, p)[mP, i])
+    return X
+end
 #
 #
 # Part II: In-place evaluations
@@ -419,7 +431,7 @@ function get_gradient!(
 ) where {FT}
     n = _vgf_index_to_length(i, vgf.range_dimension)
     pM = PowerManifold(M, range, n)
-    for (j, f) in zip(j, vgf.jacobian!![i])
+    for (j, f) in zip(i, vgf.jacobian!![i])
         f(M, X[pM, j], p)
     end
     return X
@@ -437,7 +449,7 @@ function get_gradient!(
     P = fill(p, pM)
     x = zero_vector(pM, P)
     vgf.jacobian!!(M, x, p)
-    copyto!(M, X, p, x[mP, i])
+    copyto!(M, X, p, x[pM, i])
     return X
 end
 function get_gradient!(

src/solvers/augmented_Lagrangian_method.jl

@@ -76,8 +76,8 @@ mutable struct AugmentedLagrangianMethodState{
         λ_max::R=20.0,
         λ_min::R=-λ_max,
         μ_max::R=20.0,
-        μ::V=ones(length(get_inequality_constraints(M, co, p))),
-        λ::V=ones(length(get_equality_constraints(M, co, p))),
+        μ::V=ones(length(get_inequality_constraint(M, co, p, :))),
+        λ::V=ones(length(get_equality_constraint(M, co, p, :))),
         ρ::R=1.0,
         τ::R=0.8,
         θ_ρ::R=0.3,
@@ -391,9 +391,9 @@ function augmented_Lagrangian_method!(
     ϵ_min::Real=1e-6,
     ϵ_exponent::Real=1 / 100,
     θ_ϵ::Real=(ϵ_min / ϵ)^(ϵ_exponent),
-    μ::Vector=ones(length(get_inequality_constraints(M, cmo, p))),
+    μ::Vector=ones(length(get_inequality_constraint(M, cmo, p, :))),
     μ_max::Real=20.0,
-    λ::Vector=ones(length(get_equality_constraints(M, cmo, p))),
+    λ::Vector=ones(length(get_equality_constraint(M, cmo, p, :))),
     λ_max::Real=20.0,
     λ_min::Real=-λ_max,
     τ::Real=0.8,