Allow AbstractArray input (#117)

* Allow AbstractArray input

README.md

@@ -70,9 +70,10 @@ Jacobian.
 
 This is the most efficient method, because it minimizes the memory allocations.
 
-In the following, it is assumed that have defined a function `f!(x::Vector,
-fx::Vector)` computing the residual of the system at point `x` and putting it
-into the `fx` argument.
+In the following, it is assumed that you have defined a function
+`f!(x::AbstractVector, fx::AbstractVector)` or, more generally,
+`f!(x::AbstractArray, fx::AbstractArray)` computing the residual of the system
+at point `x` and putting it into the `fx` argument.
 
 In turn, there 3 ways of specifying how the Jacobian should be computed:
 
@@ -93,6 +94,15 @@ df = DifferentiableMultivariateFunction(f!)
 nlsolve(df, initial_x)
 ```
 
+If your function `f!` operates not on `AbstractVector` but on arbitrary
+`AbstractArray` you have to specify `initial_x` in the constructor of the
+`DifferentiableMultivariateFunction`:
+
+```jl
+df = DifferentiableMultivariateFunction(f!, initial_x)
+nlsolve(df, initial_x)
+```
+
 ### Automatic differentiation
 
 Another option if you do not have a function computing the Jacobian is to use
@@ -105,14 +115,18 @@ nlsolve(f!, initial_x, autodiff = true)
 
 ### Jacobian available
 
-If, in addition to `f!`, you have a function `g!(x::Vector, gx::Array)` for
-computing the Jacobian of the system, then the syntax is, as in the example
-above:
+If, in addition to `f!(x::AbstractVector, fx::AbstractVector)`, you have a
+function `g!(x::AbstractVector, gx::AbstractMatrix)` for computing the Jacobian
+of the system, then the syntax is, as in the example above:
 
 ```jl
 nlsolve(f!, g!, initial_x)
 ```
 
+Again it is also possible to specify two functions `f!(x::AbstractArray,
+fx::AbstractArray)` and `g!(x::AbstractArray, gx::AbstractMatrix)` that work on
+arbitrary arrays `x`.
+
 Note that you should not assume that the Jacobian `gx` passed in argument is
 initialized to a zero matrix. You must set all the elements of the matrix in
 the function `g!`.
@@ -125,17 +139,32 @@ df = DifferentiableMultivariateFunction(f!, g!)
 nlsolve(df, initial_x)
 ```
 
+If `f!` and `g!` do not operate on vectors the syntax is:
+
+```jl
+df = DifferentiableMultivariateFunction(f!, g!, initial_x)
+nlsolve(df, initial_x)
+```
+
 ### Optimization of simultaneous residuals and Jacobian
 
-If, in addition to `f!` and `g!`, you have a function `fg!(x::Vector,
-fx::Vector, gx::Array)` that computes both the residual and the Jacobian at
-the same time, you can use the following syntax:
+If, in addition to `f!` and `g!`, you have a function `fg!(x::AbstractVector,
+fx::AbstractVector, gx::AbstractMatrix)` or `fg!(x::AbstractArray,
+fx::AbstractArray, gx::AbstractMatrix)` that computes both the residual and the
+Jacobian at the same time, you can use the following syntax for vectors
 
 ```jl
 df = DifferentiableMultivariateFunction(f!, g!, fg!)
 nlsolve(df, initial_x)
 ```
 
+and similarly for arbitrary arrays
+
+```jl
+df = DifferentiableMultivariateFunction(f!, g!, fg!, initial_x)
+nlsolve(df, initial_x)
+```
+
 If the function `fg!` uses some optimization that make it costless than
 calling `f!` and `g!` successively, then this syntax can possibly improve the
 performance.
@@ -147,45 +176,89 @@ cases are not optimal in terms of memory management.
 
 If only `f!` and `fg!` are available, the helper function `only_f!_and_fg!` can be
 used to construct a `DifferentiableMultivariateFunction` object, that can be
-used as first argument of `nlsolve`. The complete syntax is therefore:
+used as first argument of `nlsolve`. The complete syntax is therefore
 
 ```jl
 nlsolve(only_f!_and_fg!(f!, fg!), initial_x)
 ```
 
+if `f!` and `fg!` operate on vectors, and otherwise
+
+```jl
+nlsolve(only_f!_and_fg!(f!, fg!, initial_x), initial_x)
+```
+
 If only `fg!` is available, the helper function `only_fg!` can be used to
 construct a `DifferentiableMultivariateFunction` object, that can be used as
-first argument of `nlsolve`. The complete syntax is therefore:
+first argument of `nlsolve`. The complete syntax is therefore
 
 ```jl
 nlsolve(only_fg!(fg!), initial_x)
 ```
 
+if `f!` and `fg!` operate on vectors, and otherwise
+
+```jl
+nlsolve(only_fg!(f!, fg!, initial_x), initial_x)
+```
+
 ## With functions returning residuals and Jacobian as output
 
-Here it is assumed that you have a function `f(x::Vector)` that returns a
-newly-allocated vector containing the residuals. The helper function
-`not_in_place` can be used to construct a `DifferentiableMultivariateFunction`
-object, that can be used as first argument of `nlsolve`. The complete syntax is
-therefore:
+Here it is assumed that you have a function `f(x::AbstractVector)` that returns
+a newly-allocated vector containing the residuals. The helper function
+`not_in_place` can be used to construct a function, that can be used as first
+argument of `nlsolve`. The complete syntax is therefore:
 
 ```jl
 nlsolve(not_in_place(f), initial_x)
 ```
 
-Finite-differencing is used to compute the Jacobian in that case.
+If `f` does not operate on vectors but on arbitrary arrays the syntax is
 
-If, in addition, there is a function `g(x::Vector)` returning a newly-allocated
-matrix containing the Jacobian, 
 ```jl
-nlsolve(not_in_place(f), not_in_place(g), initial_x)
+nlsolve(not_in_place(f, initial_x), initial_x)
+```
+
+Via the `autodiff` keyword both finite-differencing and autodifferentiation can
+be used to compute the Jacobian in that case.
+
+If, in addition to `f(x::AbstractVector)`, there is a function
+`g(x::AbstractVector)` returning a newly-allocated matrix containing the
+Jacobian, `not_in_place` can be used to construct an object of type
+`DifferentiableMultivariateFunction` that can be used as first argument of
+`nlsolve`:
+
+```jl
+nlsolve(not_in_place(f, g), initial_x)
+```
+
+For functions `f` and `g` that operate on arbitrary arrays the syntax is:
+
+```jl
+nlsolve(not_in_place(f, g, initial_x), initial_x)
+```
+
+If, in addition to `f` and `g`, there is a function `fg` returning a tuple of a
+newly-allocated vector of residuals and a newly-allocated matrix of the
+Jacobian, `not_in_place` can be used to construct an object of type
+`DifferentiableMultivariateFunction`:
+
+```jl
+nlsolve(not_in_place(f, g, fg), initial_x)
+```
+
+For functions `f`, `g`, and `fg` that operate on arbitrary arrays the syntax is:
+
+```jl
+nlsolve(not_in_place(f, g, fg, initial_x), initial_x)
 ```
 
 ## With functions taking several scalar arguments
 
 If you have a function `f(x::Float64, y::Float64, ...)` that takes the point of
 interest as several scalars and returns a vector or a tuple containing the
-residuals, you can use the helper function `n_ary`. The complete syntax is therefore:
+residuals, you can use the helper function `n_ary`. The complete syntax is
+therefore:
 
 ```jl
 nlsolve(n_ary(f), initial_x)
@@ -198,13 +271,20 @@ Finite-differencing is used to compute the Jacobian.
 If the Jacobian of your function is sparse, it is possible to ask the routines
 to manipulate sparse matrices instead of full ones, in order to increase
 performance on large systems. This can be achieved by constructing an object of
-type `DifferentiableSparseMultivariateFunction`:
+type `DifferentiableSparseMultivariateFunction`. The syntax is therefore
 
 ```jl
 df = DifferentiableSparseMultivariateFunction(f!, g!)
 nlsolve(df, initial_x)
 ```
 
+if `f!` and `g!` operate on vectors, and otherwise
+
+```jl
+df = DifferentiableSparseMultivariateFunction(f!, g!, initial_x)
+nlsolve(df, initial_x)
+```
+
 It is possible to give an optional third function `fg!` to the constructor, as
 for the full Jacobian case.
 
@@ -226,14 +306,24 @@ empty!(gx.rowval)
 empty!(gx.nzval)
 ```
 
-Another solution is to directly pass a Jacobian matrix with a given sparsity. To do so, construct an object of type `DifferentiableGivenSparseMultivariateFunction`
+Another solution is to directly pass a Jacobian matrix with a given sparsity. To
+do so, construct an object of type
+`DifferentiableGivenSparseMultivariateFunction` by
 
 ```jl
 df = DifferentiableGivenSparseMultivariateFunction(f!, g!, J)
 nlsolve(df, initial_x)
 ```
 
-If  `g!` conserves the sparsity structure of `gx`, `gx` will always have the same sparsity as `J`. This sometimes allow to write a faster version of `g!`.
+if `f!` and `g!` operate on vectors, and otherwise
+
+```jl
+df = DifferentiableGivenSparseMultivariateFunction(f!, g!, J, initial_x)
+nlsolve(df, initial_x)
+```
+
+If `g!` conserves the sparsity structure of `gx`, `gx` will always have the same
+sparsity as `J`. This sometimes allow to write a faster version of `g!`.
 
 # Fine tunings
 

REQUIRE

@@ -2,6 +2,6 @@ julia 0.6.0-pre
 Calculus
 Distances
 ForwardDiff 0.5.0
-LineSearches
+LineSearches 2.1.0
 DiffBase
 NLSolversBase 2.1.0

src/anderson.jl

@@ -2,24 +2,24 @@
 # Attempts to accelerate the iteration xn+1 = xn + β f(x)
 
 @views function anderson_{T}(df::AbstractDifferentiableMultivariateFunction,
-                      x0::Vector{T},
-                      xtol::T,
-                      ftol::T,
-                      iterations::Integer,
-                      store_trace::Bool,
-                      show_trace::Bool,
-                      extended_trace::Bool,
-                      m::Integer,
-                      β :: Real)
+                             x0::AbstractArray{T},
+                             xtol::T,
+                             ftol::T,
+                             iterations::Integer,
+                             store_trace::Bool,
+                             show_trace::Bool,
+                             extended_trace::Bool,
+                             m::Integer,
+                             β::Real)
 
     f_calls = 0
     N = length(x0)
     xs = zeros(T,N,m+1) #ring buffer storing the iterates, from newest to oldest
     gs = zeros(T,N,m+1) #ring buffer storing the g of the iterates, from newest to oldest
     residuals = zeros(T, N, m) #matrix of residuals used for the least-squares problem
     alphas = zeros(T, m) #coefficients obtained by least-squares
-    fx = similar(x0) # temp variable to store f!
-    xs[:,1] = x0
+    fx = similar(x0, N) # temp variable to store f!
+    xs[:,1] = vec(x0)
     errs = zeros(iterations)
     n = 1
     tr = SolverTrace()
@@ -36,7 +36,7 @@
 
         # FIXME: How should this flag be set?
         mayterminate = false
-        
+
         if tracing
             dt = Dict()
             if extended_trace
@@ -49,9 +49,9 @@
                     n > 1 ? sqeuclidean(xs[:,1],old_x) : convert(T,NaN),
                     dt,
                     store_trace,
-                    show_trace)            
+                    show_trace)
         end
-        
+
         if converged
             break
         end
@@ -76,14 +76,16 @@
     # returning gs[:,1] rather than xs[:,1] would be better here if
     # xn+1 = xn+beta*f(xn) is convergent, but the convergence
     # criterion is not guaranteed
+    x = similar(x0)
+    copy!(x, xs[:,1])
     return SolverResults("Anderson m=$m β=$β",
-                         x0, copy(xs[:,1]), maximum(abs,fx),
+                         x0, x, maximum(abs,fx),
                          n, x_converged, xtol, f_converged, ftol, tr,
                          f_calls, 0)
 end
 
 function anderson{T}(df::AbstractDifferentiableMultivariateFunction,
-                     initial_x::Vector{T},
+                     initial_x::AbstractArray{T},
                      xtol::Real,
                      ftol::Real,
                      iterations::Integer,

src/autodiff.jl

@@ -1,16 +1,19 @@
-function autodiff(f!, initial_x::Vector)
+function autodiff(f!, initial_x::AbstractArray)
 
-    permf! = (fx, x) -> f!(x, fx)
+    fvec! = reshape_f(f!, initial_x)
+    permf! = (fx::AbstractVector, x::AbstractVector) -> fvec!(x, fx)
 
-    fx2 = copy(initial_x)
-    jac_cfg = ForwardDiff.JacobianConfig(nothing, initial_x, initial_x)
-    g! = (x, gx) -> ForwardDiff.jacobian!(gx, permf!, fx2, x, jac_cfg)
+    fx2 = vec(copy(initial_x))
+    jac_cfg = ForwardDiff.JacobianConfig(nothing, vec(initial_x), vec(initial_x))
+    function g!(x::AbstractVector, gx::AbstractMatrix)
+        ForwardDiff.jacobian!(gx, permf!, fx2, x, jac_cfg)
+    end
 
-    fg! = (x, fx, gx) -> begin
+    fg! = (x::AbstractVector, fx::AbstractVector, gx::AbstractMatrix) -> begin
         jac_res = DiffBase.DiffResult(fx, gx)
         ForwardDiff.jacobian!(jac_res, permf!, fx2, x, jac_cfg)
         DiffBase.value(jac_res)
     end
 
-    return DifferentiableMultivariateFunction(f!, g!, fg!)
+    return DifferentiableMultivariateFunction(fvec!, g!, fg!)
 end