adding IndHyperslab, fixing IndHalfspace and Linear

docs/src/functions.md

@@ -28,6 +28,7 @@ IndBinary
 IndBox  
 IndGraph     
 IndHalfspace  
+IndHyperslab
 IndPoint
 IndPolyhedral              
 IndSimplex    

src/ProximalOperators.jl

@@ -46,6 +46,7 @@ include("functions/indBox.jl")
 include("functions/indFree.jl")
 include("functions/indGraph.jl")
 include("functions/indHalfspace.jl")
+include("functions/indHyperslab.jl")
 include("functions/indNonnegative.jl")
 include("functions/indNonpositive.jl")
 include("functions/indPoint.jl")

src/functions/indHalfspace.jl

@@ -12,39 +12,40 @@ For an array `a` and a scalar `b`, returns the indicator of set
 S = \\{x : \\langle a,x \\rangle \\leq b \\}.
 ```
 """
-struct IndHalfspace{R <: Real, T <: AbstractVector{R}} <: ProximableFunction
+struct IndHalfspace{R <: Real, T <: AbstractArray{R}} <: ProximableFunction
   a::T
   b::R
-  function IndHalfspace{R,T}(a::T, b::R) where {R <: Real, T <: AbstractVector{R}}
-    norma = norm(a)
-    new(a/norma, b/norma)
+  norm_a::R
+  function IndHalfspace{R, T}(a::T, b::R) where {R <: Real, T <: AbstractArray{R}}
+    norm_a = norm(a)
+    if norm_a == 0 && b < 0
+        error("function is improper")
+    end
+    new(a, b, norm_a)
   end
 end
 
+IndHalfspace(a::T, b::R) where {R <: Real, T <: AbstractArray{R}} = IndHalfspace{R, T}(a, b)
+
 is_convex(f::IndHalfspace) = true
 is_set(f::IndHalfspace) = true
-is_cone(f::IndHalfspace) = (f.b == 0)
-
-IndHalfspace(a::T, b::R) where {R <: Real, T <: AbstractVector{R}} = IndHalfspace{R, T}(a, b)
+is_cone(f::IndHalfspace) = f.b == 0 || f.b == Inf
 
-function (f::IndHalfspace)(x::AbstractArray{T}) where T <: Real
-  s = dot(f.a, x) - f.b
-  if s <= 1e-14
-    return zero(T)
+function (f::IndHalfspace{R})(x::AbstractArray{R}) where R
+  if dot(f.a, x) - f.b <= eps(R)*f.norm_a*(1 + abs(f.b))
+    return zero(R)
   end
-  return +Inf
+  return R(Inf)
 end
 
-function prox!(y::AbstractArray{T}, f::IndHalfspace, x::AbstractArray{T}, gamma::Real=1.0) where T <: Real
-  s = dot(f.a, x) - f.b
-  if s <= 0
-    y .= x
-    return zero(T)
-  end
-  for k in eachindex(x)
-    y[k] = x[k] - s*f.a[k]
+function prox!(y::AbstractArray{R}, f::IndHalfspace{R}, x::AbstractArray{R}, gamma::R=one(R)) where R
+  s = dot(f.a, x)
+  if s > f.b
+    y .= x .- ((s - f.b)/f.norm_a^2) .* f.a
+  else
+    copyto!(y, x)
   end
-  return zero(T)
+  return zero(R)
 end
 
 fun_name(f::IndHalfspace) = "indicator of a halfspace"
@@ -54,10 +55,10 @@ fun_params(f::IndHalfspace) =
   string( "a = ", typeof(f.a), " of size ", size(f.a), ", ",
           "b = $(f.b)")
 
-function prox_naive(f::IndHalfspace, x::AbstractArray{T}, gamma::Real=1.0) where T <: Real
+function prox_naive(f::IndHalfspace{R}, x::AbstractArray{R}, gamma::R=one(R)) where R
   s = dot(f.a, x) - f.b
   if s <= 0
     return x, 0.0
   end
-  return x - s*f.a, 0.0
+  return x - (s/norm(f.a)^2)*f.a, 0.0
 end

src/functions/indHyperslab.jl

@@ -0,0 +1,70 @@
+# indicator of a hyperslab
+
+export IndHyperslab
+
+"""
+**Indicator of a hyperslab**
+
+    IndHalfspace(low, a, upp)
+
+For an array `a` and scalars `low` and `upp`, returns the indicator of set
+```math
+S = \\{x : low \\leq \\langle a,x \\rangle \\leq upp \\}.
+```
+"""
+struct IndHyperslab{R <: Real, T <: AbstractArray{R}} <: ProximableFunction
+  low::R
+  a::T
+  upp::R
+  norm_a::R
+  function IndHyperslab{R, T}(low::R, a::T, upp::R) where {R <: Real, T <: AbstractArray{R}}
+    norm_a = norm(a)
+    if (norm_a == 0 && (upp < 0 || low > 0)) || upp < low
+        error("function is improper")
+    end
+    new(low, a, upp, norm_a)
+  end
+end
+
+IndHyperslab(low::R, a::T, upp::R) where {R <: Real, T <: AbstractArray{R}} = IndHyperslab{R, T}(low, a, upp)
+
+is_convex(f::IndHyperslab) = true
+is_set(f::IndHyperslab) = true
+is_cone(f::IndHyperslab{R}) where R =
+    (f.low == f.upp == 0) ||
+    (f.low == 0 && f.upp == Inf) ||
+    (f.low == -Inf && f.upp == 0) ||
+    (f.low == -Inf && f.upp == Inf)
+
+function (f::IndHyperslab{R})(x::AbstractArray{R}) where R
+  s = dot(f.a, x)
+  # if f.low <= s <= f.upp
+  # tol = (R(1) + abs(s))*eps(R)
+  if f.low - s <= eps(R)*f.norm_a*(1 + abs(f.low)) && s - f.upp <= eps(R)*f.norm_a*(1 + abs(f.upp))
+    return zero(R)
+  end
+  return R(Inf)
+end
+
+function prox!(y::AbstractArray{R}, f::IndHyperslab{R}, x::AbstractArray{R}, gamma::R=one(R)) where R
+  s = dot(f.a, x)
+  if s < f.low && f.norm_a > 0
+    y .= x .- ((s - f.low)/f.norm_a^2) .* f.a
+  elseif s > f.upp && f.norm_a > 0
+    y .= x .- ((s - f.upp)/f.norm_a^2) .* f.a
+  else
+    copyto!(y, x)
+  end
+  return zero(R)
+end
+
+function prox_naive(f::IndHyperslab{R}, x::AbstractArray{R}, gamma::R=one(R)) where R
+  s = dot(f.a, x)
+  if s < f.low && f.norm_a > 0
+    return x - ((s - f.low)/norm(f.a)^2) * f.a, 0
+  elseif s > f.upp && f.norm_a > 0
+    return x - ((s - f.upp)/norm(f.a)^2) * f.a, 0
+  else
+    return x, 0
+  end
+end

src/functions/linear.jl

@@ -10,25 +10,25 @@ Returns the function
 f(x) = \\langle c, x \\rangle.
 ```
 """
-struct Linear{R <: RealOrComplex, A <: AbstractArray{R}} <: ProximableFunction
+struct Linear{R <: Real, A <: AbstractArray{R}} <: ProximableFunction
   c::A
 end
 
 is_separable(f::Linear) = true
 is_convex(f::Linear) = true
 is_smooth(f::Linear) = true
 
-function (f::Linear{RC, A})(x::AbstractArray{RC}) where {R <: Real, RC <: Union{R, Complex{R}}, A <: AbstractArray{RC}}
+function (f::Linear{R})(x::AbstractArray{R}) where R
   return dot(f.c, x)
 end
 
-function prox!(y::AbstractArray{RC}, f::Linear{RC, A}, x::AbstractArray{RC}, gamma::Union{R, AbstractArray{R}}=1.0) where {R <: Real, RC <: Union{R, Complex{R}}, A <: AbstractArray{RC}}
+function prox!(y::AbstractArray{R}, f::Linear{R}, x::AbstractArray{R}, gamma::Union{R, AbstractArray{R}}=1.0) where R
   y .= x .- gamma.*(f.c)
   fy = dot(f.c, y)
   return fy
 end
 
-function gradient!(y::AbstractArray{RC}, f::Linear{RC, A}, x::AbstractArray{RC}) where {R <: Real, RC <: Union{R, Complex{R}}, A <: AbstractArray{RC}}
+function gradient!(y::AbstractArray{R}, f::Linear{R}, x::AbstractArray{R}) where R
   y .= f.c
   return dot(f.c, x)
 end

test/runtests.jl

@@ -78,6 +78,7 @@ end
   include("test_logisticLoss.jl")
   include("test_quadratic.jl")
   include("test_linear.jl")
+  include("test_indhyperslab.jl")
   include("test_calls.jl")
   include("test_graph.jl")
 end

test/test_indhyperslab.jl

@@ -0,0 +1,25 @@
+using Test
+using Random
+using LinearAlgebra
+using ProximalOperators
+
+Random.seed!(0)
+
+@testset "IndHyperslab" begin
+
+for R in [Float16, Float32, Float64]
+    for shape in [(5,), (3, 5), (3, 4, 5)]
+        c = randn(R, shape)
+        x = randn(R, shape)
+        cx = dot(c, x)
+
+        for (low, upp) in [(cx-R(1), cx+R(1)), (cx-R(2), cx-R(1)), (cx+R(1), cx+R(2))]
+            f = IndHyperslab(low, c, upp)
+            predicates_test(f)
+            call_test(f, x)
+            prox_test(f, x, R(0.5)+rand(R))
+        end
+    end
+end
+
+end

test/test_linear.jl

@@ -3,17 +3,18 @@ using Random
 
 Random.seed!(0)
 
+@testset "Linear" begin
+
 for R in [Float16, Float32, Float64]
-    for T in [R, Complex{R}]
-        for shape in [(5,), (3, 5), (3, 4, 5)]
-            # Real
-            c = randn(T, shape)
-            f = Linear(c)
-            predicates_test(f)
-            x = randn(T, shape)
-            @test gradient(f, x) == (c, f(x))
-            call_test(f, x)
-            prox_test(f, x, R(0.5)+rand(R))
-        end
+    for shape in [(5,), (3, 5), (3, 4, 5)]
+        c = randn(R, shape)
+        f = Linear(c)
+        predicates_test(f)
+        x = randn(R, shape)
+        @test gradient(f, x) == (c, f(x))
+        call_test(f, x)
+        prox_test(f, x, R(0.5)+rand(R))
     end
 end
+
+end