Removing BlockArrays

.travis.yml

@@ -4,8 +4,8 @@ os:
   - linux
   - osx
 julia:
-  - 0.7
   - 1.0
+  - 1.1
   - nightly
 matrix:
   allow_failures:

appveyor.yml

@@ -1,7 +1,7 @@
 environment:
   matrix:
-  - julia_version: 0.7
   - julia_version: 1
+  - julia_version: 1.1
   - julia_version: nightly
 
 platform:
@@ -10,7 +10,6 @@ platform:
 
 matrix:
   allow_failures:
-  - julia_version: 1
   - julia_version: nightly
 
 branches:

src/ProximalAlgorithms.jl

@@ -1,7 +1,6 @@
 module ProximalAlgorithms
 
 using AbstractOperators
-using AbstractOperators.BlockArrays
 using ProximalOperators
 using LinearAlgebra
 using Printf
@@ -11,6 +10,7 @@ import Base: iterate
 include("utilities/identity.jl")
 include("utilities/zero.jl")
 include("utilities/conjugate.jl")
+include("utilities/tomove.jl")
 
 abstract type ProximalAlgorithm{I,T} end
 

src/algorithms/AsymmetricForwardBackwardAdjoint.jl

@@ -5,7 +5,7 @@
 # [2] Condat. "A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms" Journal of Optimization Theory and Applications 158.2 (2013): 460-479.
 # [3] Vũ. "A splitting algorithm for dual monotone inclusions involving cocoercive operators"" Advances in Computational Mathematics, 38(3), pp.667-681.
 
-struct AFBAIterator{I <: Integer, R <: Real, T1 <: BlockArray, T2 <: BlockArray} <: ProximalAlgorithm{I, Tuple{T1, T2}}
+struct AFBAIterator{I <: Integer, R <: Real, T1 <: AbstractArray, T2 <: AbstractArray} <: ProximalAlgorithm{I, Tuple{T1, T2}}
     x::T1
     y::T2
     g
@@ -37,19 +37,19 @@ end
 ################################################################################
 # Constructor(s)
 
-function AFBAIterator(x0::T1, y0::T2; g=IndFree(), h=IndFree(), f=IndFree(), l=IndZero(), L=Identity(blocksize(x0)), theta=1, mu=1, lam=1, betaQ=0, betaR=0, gamma1=-1, gamma2=-1, maxit::I=10000, tol::R=1e-5, verbose=1, verbose_freq = 100) where {I, R, T1, T2}
-    x = blockcopy(x0)
-    xbar = blockcopy(x0)
-    y = blockcopy(y0)
-    ybar = blockcopy(y0)
-    gradf = blockcopy(x0)
-    gradl = blockcopy(y0)
-    FPR_x = blockcopy(x0)
+function AFBAIterator(x0::T1, y0::T2; g=IndFree(), h=IndFree(), f=IndFree(), l=IndZero(), L=Identity(size(x0)), theta=1, mu=1, lam=1, betaQ=0, betaR=0, gamma1=-1, gamma2=-1, maxit::I=10000, tol::R=1e-5, verbose=1, verbose_freq = 100) where {I, R, T1, T2}
+    x = copy(x0)
+    xbar = copy(x0)
+    y = copy(y0)
+    ybar = copy(y0)
+    gradf = copy(x0)
+    gradl = copy(y0)
+    FPR_x = copy(x0)
     FPR_x .= Inf
-    FPR_y = blockcopy(y0)
+    FPR_y = copy(y0)
     FPR_y .= Inf
-    temp_x = blockcopy(x0)
-    temp_y = blockcopy(y0)
+    temp_x = copy(x0)
+    temp_y = copy(y0)
     hconj = Conjugate(h)
     lconj = Conjugate(l)
 

src/algorithms/DouglasRachford.jl

@@ -5,7 +5,7 @@
 # [1] Eckstein, Bertsekas "On the Douglas-Rachford Splitting Method and the Proximal Point Algorithm for Maximal Monotone Operators*", Mathematical Programming, vol. 55, no. 1, pp. 293-318 (1989).
 #
 
-struct DRSIterator{I <: Integer, R <: Real, T <: BlockArray{R}} <: ProximalAlgorithm{I,T}
+struct DRSIterator{I <: Integer, R <: Real, T <: AbstractArray{R}} <: ProximalAlgorithm{I,T}
     x::T
     f
     g
@@ -23,11 +23,11 @@ end
 ################################################################################
 # Constructor(s)
 
-function DRSIterator(x0::BlockArray{R}; f=Zero(), g=Zero(), gamma::R=1.0, maxit::I=10000, tol::R=1e-4, verbose=1, verbose_freq = 100) where {I, R}
-    y = blockcopy(x0)
-    r = blockcopy(x0)
-    z = blockcopy(x0)
-    FPR_x = blockcopy(x0)
+function DRSIterator(x0::AbstractArray{R}; f=Zero(), g=Zero(), gamma::R=1.0, maxit::I=10000, tol::R=1e-4, verbose=1, verbose_freq = 100) where {I, R}
+    y = copy(x0)
+    r = copy(x0)
+    z = copy(x0)
+    FPR_x = copy(x0)
     FPR_x .= Inf
     return DRSIterator{I, R, typeof(x0)}(x0, f, g, gamma, maxit, tol, verbose, verbose_freq, y, r, z, FPR_x)
 end
@@ -37,7 +37,7 @@ end
 
 maxit(sol::DRSIterator) = sol.maxit
 
-converged(sol::DRSIterator, it) = it > 0 && blockmaxabs(sol.FPR_x)/sol.gamma <= sol.tol
+converged(sol::DRSIterator, it) = it > 0 && maximum(abs,sol.FPR_x)/sol.gamma <= sol.tol
 
 verbose(sol::DRSIterator) = sol.verbose > 0
 verbose(sol::DRSIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0))
@@ -48,12 +48,12 @@ function display(sol::DRSIterator)
 end
 
 function display(sol::DRSIterator, it)
-    @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, blockmaxabs(sol.FPR_x)/sol.gamma)
+    @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, maximum(abs,sol.FPR_x)/sol.gamma)
 end
 
 function Base.show(io::IO, sol::DRSIterator)
     println(io, "Douglas-Rachford Splitting" )
-    println(io, "fpr        : $(blockmaxabs(sol.FPR_x))")
+    println(io, "fpr        : $(maximum(abs,sol.FPR_x))")
     print(  io, "gamma      : $(sol.gamma)")
 end
 

src/algorithms/ForwardBackward.jl

@@ -50,31 +50,31 @@ end
 # Constructor
 
 function FBSIterator(x0::D; 
-                     fs::FS=Zero(), As::AS=Identity(blocksize(x0)), 
-                     fq::FQ=Zero(), Aq::AQ=Identity(blocksize(x0)), 
+                     fs::FS=Zero(), As::AS=Identity(size(x0)), 
+                     fq::FQ=Zero(), Aq::AQ=Identity(size(x0)), 
                      g::G=Zero(), 
                      gamma::R=-1.0, maxit::I=10000, tol::R=1e-4, 
                      adaptive::Bool=false, fast::Bool=false, 
                      verbose::I=1, verbose_freq::I = 100) where {I, R, D, FS, AS, FQ, AQ, G}
-    x = blockcopy(x0)
-    y = blockzeros(x0)
-    z = blockzeros(x0)
-    z_prev = blockzeros(x0)
-    FPR_x = blockzeros(x0)
+    x = copy(x0)
+    y = zero(x0)
+    z = zero(x0)
+    z_prev = zero(x0)
+    FPR_x = zero(x0)
     Aqx = Aq*x
     Asx = As*x
-    Aqz = blockzeros(Aqx)
-    Asz = blockzeros(Asx)
-    gradfq_Aqx = blockzeros(Aqx)
-    gradfs_Asx = blockzeros(Asx)
-    Aqt_gradfq_Aqx = blockzeros(x0)
-    Ast_gradfs_Asx = blockzeros(x0)
-    At_gradf_Ax = blockzeros(x0)
-    Aqz_prev = blockzeros(Aqx)
-    Asz_prev = blockzeros(Asx)
-    gradfq_Aqz = blockzeros(Aqx)
-    gradfs_Asz = blockzeros(Asx)
-    gradfq_Aqz_prev = blockzeros(Aqx)
+    Aqz = zero(Aqx)
+    Asz = zero(Asx)
+    gradfq_Aqx = zero(Aqx)
+    gradfs_Asx = zero(Asx)
+    Aqt_gradfq_Aqx = zero(x0)
+    Ast_gradfs_Asx = zero(x0)
+    At_gradf_Ax = zero(x0)
+    Aqz_prev = zero(Aqx)
+    Asz_prev = zero(Asx)
+    gradfq_Aqz = zero(Aqx)
+    gradfs_Asz = zero(Asx)
+    gradfq_Aqz_prev = zero(Aqx)
     CQ = typeof(Aqx)
     CS = typeof(Asx)
     FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}(x,
@@ -96,7 +96,7 @@ end
 
 maxit(sol::FBSIterator{I}) where {I} = sol.maxit
 
-converged(sol::FBSIterator{I, R, D}, it::I) where {I, R, D}  = it > 0 && blockmaxabs(sol.FPR_x)/sol.gamma <= sol.tol
+converged(sol::FBSIterator{I, R, D}, it::I) where {I, R, D}  = it > 0 && maximum(abs,sol.FPR_x)/sol.gamma <= sol.tol
 
 verbose(sol::FBSIterator) = sol.verbose > 0
 verbose(sol::FBSIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0))
@@ -107,12 +107,12 @@ function display(sol::FBSIterator)
 end
 
 function display(sol::FBSIterator, it)
-    @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, blockmaxabs(sol.FPR_x)/sol.gamma)
+    @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, maximum(abs,sol.FPR_x)/sol.gamma)
 end
 
 function Base.show(io::IO, sol::FBSIterator)
     println(io, (sol.fast ? "Fast " : "")*"Forward-Backward Splitting" )
-    println(io, "fpr        : $(blockmaxabs(sol.FPR_x))")
+    println(io, "fpr        : $(maximum(abs,sol.FPR_x))")
     print(  io, "gamma      : $(sol.gamma)")
 end
 
@@ -131,30 +131,29 @@ function initialize!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}) where
     sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asx)
     mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx)
     mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx)
-    blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx)
+    sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx
     sol.f_Ax = sol.fs_Asx + sol.fq_Aqx
 
     if sol.gamma <= 0.0 # estimate L in this case, and set gamma = 1/L
         # this part should be as follows:
         # 1) if adaptive = false and only fq is present then L is "accurate"
         # 2) otherwise L is "inaccurate" and set adaptive = true
         # TODO: implement case 1), now 2) is always performed
-        # xeps = sol.x .+ sqrt(eps())
-        xeps = (x -> x .+ sqrt(eps())).(sol.x)
+        xeps = sol.x .+ sqrt(eps())
         Aqxeps = sol.Aq*xeps
         gradfq_Aqxeps, = gradient(sol.fq, Aqxeps)
         Asxeps = sol.As*xeps
         gradfs_Asxeps, = gradient(sol.fs, Asxeps)
         At_gradf_Axeps = sol.As'*gradfs_Asxeps .+ sol.Aq'*gradfq_Aqxeps
-        L = blockvecnorm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*blocklength(xeps)))
+        L = norm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*length(xeps)))
         sol.adaptive = true
         # in both cases set gamma = 1/L
         sol.gamma = 1.0/L
     end
 
-    blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax)
+    sol.y .= sol.x .- sol.gamma .* sol.At_gradf_Ax
     prox!(sol.z, sol.g, sol.y, sol.gamma)
-    blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.z)
+    sol.FPR_x .= sol.x .- sol.z
 
     return sol.z
 end
@@ -169,18 +168,18 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w
 
     if sol.adaptive
         for it_gam = 1:100 # TODO: replace/complement with lower bound on gamma
-            normFPR_x = blockvecnorm(sol.FPR_x)
-            uppbnd = sol.f_Ax - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_x^2
+            normFPR_x = norm(sol.FPR_x)
+            uppbnd = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_x^2
             mul!(sol.Aqz, sol.Aq, sol.z)
             fq_Aqz = gradient!(sol.gradfq_Aqz, sol.fq, sol.Aqz)
             mul!(sol.Asz, sol.As, sol.z)
             fs_Asz = gradient!(sol.gradfs_Asz, sol.fs, sol.Asz)
             f_Az = fs_Asz + fq_Aqz
             if f_Az > uppbnd + 1e-6*abs(sol.f_Ax)
                 sol.gamma = 0.5*sol.gamma
-                blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax)
+                sol.y .= sol.x - sol.gamma .* sol.At_gradf_Ax
                 prox!(sol.z, sol.g, sol.y, sol.gamma)
-                blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.z)
+                sol.FPR_x .= sol.x .- sol.z
             else
                 break
             end
@@ -198,7 +197,7 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w
             sol.gradfs_Asx, sol.gradfs_Asz = sol.gradfs_Asz, sol.gradfs_Asx
             mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx)
             mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx)
-            blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx)
+            sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx
             sol.f_Ax = sol.fs_Asx + sol.fq_Aqx
         end
     else
@@ -207,27 +206,19 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w
         extr = (sol.theta - 1.0)/theta1
         sol.theta = theta1
         # perform extrapolation
-        #sol.x .= sol.z .+ extr.*(sol.z .- sol.z_prev)
-        blockaxpy!(sol.x, sol.z, -1, sol.z_prev)
-        blockaxpy!(sol.x, sol.z, extr, sol.x)
+        sol.x .= sol.z .+ extr.*(sol.z .- sol.z_prev)
 
         sol.z, sol.z_prev = sol.z_prev, sol.z
         if sol.adaptive == true
             # extrapolate other extrapolable quantities
-            # diff_Aqx = (sol.Aqz .- sol.Aqz_prev)
-            # diff_gradfq_Aqx = (sol.gradfq_Aqz .- sol.gradfq_Aqz_prev)
-            blockaxpy!(sol.Aqx, sol.Aqz, -1.0, sol.Aqz_prev)
-            blockaxpy!(sol.gradfq_Aqx, sol.gradfq_Aqz, -1.0, sol.gradfq_Aqz_prev)
-            coeff_linear = real(blockvecdot(sol.gradfq_Aqz, sol.Aqx))
-            coeff_quadr = real(blockvecdot(sol.Aqx, sol.gradfq_Aqx))
+            sol.Aqx .= sol.Aqz .- sol.Aqz_prev
+            sol.gradfq_Aqx .= sol.gradfq_Aqz .- sol.gradfq_Aqz_prev
+            coeff_linear = real(dot(sol.gradfq_Aqz, sol.Aqx))
+            coeff_quadr = real(dot(sol.Aqx, sol.gradfq_Aqx))
             sol.fq_Aqx = fq_Aqz + extr*coeff_linear + 0.5*coeff_quadr*extr^2
-            # sol.Aqx .= sol.Aqz .+ extr.*diff_Aqx
-            blockaxpy!(sol.Aqx, sol.Aqz, extr, sol.Aqx)
-            # sol.gradfq_Aqx .= sol.gradfq_Aqz .+ extr.*diff_gradfq_Aqx
-            blockaxpy!(sol.gradfq_Aqx, sol.gradfq_Aqz, extr, sol.gradfq_Aqx)
-            # sol.Asx .= sol.Asz .+ extr.*(sol.Asz .- sol.Asz_prev)
-            blockaxpy!(sol.Asx, sol.Asz, -1.0, sol.Asz_prev)
-            blockaxpy!(sol.Asx, sol.Asz, extr, sol.Asx     )
+            sol.Aqx .= sol.Aqz .+ extr.*(sol.Aqz .- sol.Aqz_prev)
+            sol.gradfq_Aqx .= sol.gradfq_Aqz .+ extr.*(sol.gradfq_Aqz .- sol.gradfq_Aqz_prev)
+            sol.Asx .= sol.Asz .+ extr .*(sol.Asz .- sol.Asz_prev)
             # store the z-quantities for future extrapolation
             sol.Aqz_prev, sol.Aqz = sol.Aqz, sol.Aqz_prev
             sol.gradfq_Aqz_prev, sol.gradfq_Aqz = sol.gradfq_Aqz, sol.gradfq_Aqz_prev
@@ -237,7 +228,7 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w
             mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx)
             # TODO: we can probably save the MATVEC by Aq' in the next line
             mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx)
-            blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx)
+            sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx
             sol.f_Ax = sol.fs_Asx + sol.fq_Aqx
         end
     end
@@ -248,12 +239,12 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w
         sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asx)
         mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx)
         mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx)
-        blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx)
+        sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx
         sol.f_Ax = sol.fs_Asx + sol.fq_Aqx
     end
-    blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax)
+    sol.y .= sol.x .- sol.gamma .* sol.At_gradf_Ax
     prox!(sol.z, sol.g, sol.y, sol.gamma)
-    blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.z)
+    sol.FPR_x .= sol.x .- sol.z
 
     return sol.z