Fix array input.

src/differentiable_vectors/autodiff.jl

@@ -1,35 +1,30 @@
-# if only f! and vector shaped x is given, construct j! or fj! according to autodiff
 function OnceDifferentiable(f!, F::AbstractArray, x::AbstractArray, autodiff = :central)
-    fv!(fxv::AbstractVector, xv::AbstractVector)  = f!(reshape(fxv, size(F)...), reshape(xv, size(x)...))
-    OnceDifferentiable(fv!, vec(F), vec(x), autodiff)
-end
-function OnceDifferentiable(f!, F::AbstractVector, x::AbstractVector, autodiff = :central)
     if autodiff == :central
-        function fj!(fx, jx, x)
-            f!(fx, x)
+        function fj!(F, J, x)
+            f!(F, x)
             function f(x)
-                fx = similar(x)
-                f!(fx, x)
-                return fx
+                F = similar(x)
+                f!(F, x)
+                return F
             end
-            finite_difference_jacobian!(f, x, fx, jx, autodiff)
+            finite_difference_jacobian!(f, x, F, J, autodiff)
         end
-        function j!(jx, x)
-            fx = similar(x)
-            fj!(fx, jx, x)
+        function j!(J, x)
+            F = similar(x)
+            fj!(F, J, x)
         end
         return OnceDifferentiable(f!, j!, fj!, similar(x), similar(x))
     elseif autodiff == :forward || autodiff == true
         jac_cfg = ForwardDiff.JacobianConfig(f!, x, x)
         ForwardDiff.checktag(jac_cfg, f!, x)
 
-        fx2 = copy(x)   
-        function g!(gx, x)
-            ForwardDiff.jacobian!(gx, f!, fx2, x, jac_cfg, Val{false}())
+        F2 = copy(x)   
+        function g!(J, x)
+            ForwardDiff.jacobian!(J, f!, F2, x, jac_cfg, Val{false}())
         end
-        function fg!(fx, gx, x)
-            jac_res = DiffBase.DiffResult(fx, gx)
-            ForwardDiff.jacobian!(jac_res, f!, fx2, x, jac_cfg, Val{false}())
+        function fg!(F, J, x)
+            jac_res = DiffBase.DiffResult(F, J)
+            ForwardDiff.jacobian!(jac_res, f!, F2, x, jac_cfg, Val{false}())
             DiffBase.value(jac_res)
         end
 

src/nlsolve/utils.jl

@@ -1,13 +1,13 @@
-function wdot{T}(wx::AbstractVector{T}, x::AbstractVector{T},
-                 wy::AbstractVector{T}, y::AbstractVector{T})
+function wdot{T}(wx::AbstractArray{T}, x::AbstractArray{T},
+                 wy::AbstractArray{T}, y::AbstractArray{T})
     out = zero(T)
     @inbounds @simd for i in 1:length(x)
         out += wx[i]*x[i] * wy[i]*y[i]
     end
     return out
 end
 
-wnorm{T}(w::AbstractVector{T}, x::AbstractVector{T}) = sqrt(wdot(w, x, w, x))
+wnorm{T}(w::AbstractArray{T}, x::AbstractArray{T}) = sqrt(wdot(w, x, w, x))
 assess_convergence(f, ftol) = assess_convergence(NaN, NaN, f, NaN, ftol)
 function assess_convergence(x,
                             x_previous,
@@ -29,7 +29,7 @@ function assess_convergence(x,
     return x_converged, f_converged, converged
 end
 
-function check_isfinite(x::AbstractVector)
+function check_isfinite(x::AbstractArray)
     i = find((!).(isfinite.(x)))
     if !isempty(i)
         throw(IsFiniteException(i))

src/solvers/mcp.jl

@@ -15,26 +15,26 @@ end
 function mcp_smooth(df::OnceDifferentiable,
                     lower::Vector, upper::Vector)
 
-    function f!(fx, x)
-        value!(df, fx, x)
+    function f!(F, x)
+        value!(df, F, x)
         for i = 1:length(x)
             if  isfinite.(upper[i])
-                fx[i] += (x[i]-upper[i]) + sqrt(fx[i]^2+(x[i]-upper[i])^2)
+                F[i] += (x[i]-upper[i]) + sqrt(F[i]^2+(x[i]-upper[i])^2)
             end
             if  isfinite.(lower[i])
-                fx[i] += (x[i]-lower[i]) - sqrt(fx[i]^2+(x[i]-lower[i])^2)
+                F[i] += (x[i]-lower[i]) - sqrt(F[i]^2+(x[i]-lower[i])^2)
             end
         end
     end
 
-    function j!(gx, x)
-        fx = similar(x)
-        value_jacobian!(df, fx, gx, x)
+    function j!(J, x)
+        F = similar(x)
+        value_jacobian!(df, F, J, x)
 
         # Derivatives of phiplus
-        sqplus = sqrt.(fx.^2 .+ (x .- upper).^2)
+        sqplus = sqrt.(F.^2 .+ (x .- upper).^2)
 
-        dplus_du = 1 + fx./sqplus
+        dplus_du = 1 .+ F./sqplus
 
         dplus_dv = similar(x)
         for i = 1:length(x)
@@ -46,7 +46,7 @@ function mcp_smooth(df::OnceDifferentiable,
         end
 
         # Derivatives of phiminus
-        phiplus = copy(fx)
+        phiplus = copy(F)
         for i = 1:length(x)
             if isfinite(upper[i])
                 phiplus[i] += (x[i]-upper[i]) + sqplus[i]
@@ -55,7 +55,7 @@ function mcp_smooth(df::OnceDifferentiable,
 
         sqminus = sqrt.(phiplus.^2 .+ (x .- lower).^2)
 
-        dminus_du = 1-phiplus./sqminus
+        dminus_du = 1.-phiplus./sqminus
 
         dminus_dv = similar(x)
         for i = 1:length(x)
@@ -69,12 +69,12 @@ function mcp_smooth(df::OnceDifferentiable,
         # Final computations
         for i = 1:length(x)
             for j = 1:length(x)
-                gx[i,j] *= dminus_du[i]*dplus_du[i]
+                J[i,j] *= dminus_du[i]*dplus_du[i]
             end
-            gx[i,i] += dminus_dv[i] + dminus_du[i]*dplus_dv[i]
+            J[i,i] += dminus_dv[i] + dminus_du[i]*dplus_dv[i]
         end
     end
-    return OnceDifferentiable(f!, j!, similar(lower), jacobian(df), similar(lower))
+    return OnceDifferentiable(f!, j!, similar(df.F), similar(df.x_f))
 end
 
 # Generate a function whose roots are the solutions of the MCP.
@@ -86,28 +86,28 @@ end
 # hence the difference in the function.
 function mcp_minmax(df::OnceDifferentiable,
                     lower::Vector, upper::Vector)
-    function f!(fx, x)
-        value!(df, fx, x)
+    function f!(F, x)
+        value!(df, F, x)
         for i = 1:length(x)
-            if fx[i] < x[i]-upper[i]
-                fx[i] = x[i]-upper[i]
+            if F[i] < x[i]-upper[i]
+                F[i] = x[i]-upper[i]
             end
-            if fx[i] > x[i]-lower[i]
-                fx[i] = x[i]-lower[i]
+            if F[i] > x[i]-lower[i]
+                F[i] = x[i]-lower[i]
             end
         end
     end
 
-    function j!(gx, x)
-        fx = similar(x)
-        value_jacobian!(df, fx, gx, x)
+    function j!(J, x)
+        F = similar(x)
+        value_jacobian!(df, F, J, x)
         for i = 1:length(x)
-            if fx[i] < x[i]-upper[i] || fx[i] > x[i]-lower[i]
+            if F[i] < x[i]-upper[i] || F[i] > x[i]-lower[i]
                 for j = 1:length(x)
-                    gx[i,j] = (i == j ? 1.0 : 0.0)
+                    J[i,j] = (i == j ? 1.0 : 0.0)
                 end
             end
         end
     end
-    return OnceDifferentiable(f!, j!, similar(lower), similar(lower))
+    return OnceDifferentiable(f!, j!, similar(df.F), similar(df.x_f))
 end

src/solvers/mcp_func_defs.jl

@@ -76,11 +76,9 @@ function mcpsolve{T}(f!,
                   factor::Real = one(T),
                   autoscale::Bool = true,
                   autodiff::Bool = false)
-    if !autodiff
-        df = OnceDifferentiable(f!, similar(initial_x), initial_x)
-    else
+#    if !autodiff
         df = OnceDifferentiable(f!, initial_x, initial_x)
-    end
+#    end
     @reformulate df
     nlsolve(rf,
             initial_x, method = method, xtol = xtol, ftol = ftol,

src/solvers/newton.jl

@@ -67,7 +67,7 @@ function newton_{T}(df::OnceDifferentiable,
         if xlin != xold
             value!(df, xlin)
         end
-        dot(value(df), value(df)) / 2
+        vecdot(value(df), value(df)) / 2
     end
 
     # The line search algorithm will want to first compute ∇fo(xₖ).
@@ -79,11 +79,11 @@ function newton_{T}(df::OnceDifferentiable,
         if xlin != xold
             value_jacobian!(df, xlin)
         end
-        At_mul_B!(storage, jacobian(df), value(df))
+        At_mul_B!(storage, jacobian(df), vec(value(df)))
     end
     function fgo!(storage::AbstractVector, xlin::AbstractVector)
         go!(storage, xlin)
-        dot(value(df), value(df)) / 2
+        vecdot(value(df), value(df)) / 2
     end
     dfo = OnceDifferentiable(fo, go!, fgo!, real(zero(T)), x)
 
@@ -96,16 +96,16 @@ function newton_{T}(df::OnceDifferentiable,
         end
 
         try
-            At_mul_B!(g, jacobian(df), value(df))
-            p = jacobian(df)\value(df)
+            At_mul_B!(g, jacobian(df), vec(value(df)))
+            p = jacobian(df)\vec(value(df))
             scale!(p, -1)
         catch e
             if isa(e, Base.LinAlg.LAPACKException) || isa(e, Base.LinAlg.SingularException)
                 # Modify the search direction if the jacobian is singular
                 # FIXME: better selection for lambda, see Nocedal & Wright p. 289
                 fjac2 = jacobian(df)'*jacobian(df)
                 lambda = convert(T,1e6)*sqrt(n*eps())*norm(fjac2, 1)
-                p = -(fjac2 + lambda*eye(n))\g
+                p = -(fjac2 + lambda*eye(n))\vec(g)
             else
                 throw(e)
             end
@@ -114,7 +114,7 @@ function newton_{T}(df::OnceDifferentiable,
         copy!(xold, x)
 
         LineSearches.clear!(lsr)
-        push!(lsr, zero(T), dot(value(df),value(df))/2, dot(g, p))
+        push!(lsr, zero(T), vecdot(value(df),value(df))/2, vecdot(g, p))
 
         alpha = linesearch!(dfo, xold, p, x, lsr, one(T), mayterminate)
 
@@ -126,7 +126,7 @@ function newton_{T}(df::OnceDifferentiable,
     end
 
     return SolverResults("Newton with line-search",
-                         initial_x, reshape(x, size(initial_x)...), norm(value(df), Inf),
+                         initial_x, reshape(x, size(initial_x)...), vecnorm(value(df), Inf),
                          it, x_converged, xtol, f_converged, ftol, tr,
                          first(df.f_calls), first(df.df_calls))
 end

src/solvers/trust_region.jl

@@ -22,18 +22,18 @@ macro trustregiontrace(stepnorm)
     end)
 end
 
-function dogleg!{T}(p::AbstractVector{T}, r::AbstractVector{T}, d::AbstractVector{T},
+function dogleg!{T}(p::AbstractArray{T}, r::AbstractArray{T}, d::AbstractArray{T},
                     J::AbstractMatrix{T}, delta::Real)
     local p_i
     try
-        p_i = J \ r # Gauss-Newton step
+        p_i = J \ vec(r) # Gauss-Newton step
     catch e
         if isa(e, Base.LinAlg.LAPACKException) || isa(e, Base.LinAlg.SingularException)
             # If Jacobian is singular, compute a least-squares solution to J*x+r=0
             U, S, V = svd(full(J)) # Convert to full matrix because sparse SVD not implemented as of Julia 0.3
             k = sum(S .> eps())
             mrinv = V * diagm([1./S[1:k]; zeros(eltype(S), length(S)-k)]) * U' # Moore-Penrose generalized inverse of J
-            p_i = mrinv * r
+            p_i = mrinv * vec(r)
         else
             throw(e)
         end
@@ -51,11 +51,11 @@ function dogleg!{T}(p::AbstractVector{T}, r::AbstractVector{T}, d::AbstractVecto
 
         # compute g = J'r ./ (d .^ 2)
         g = p
-        At_mul_B!(g, J, r)
+        At_mul_B!(vec(g), J, vec(r))
         g .= g ./ d.^2
 
         # compute Cauchy point
-        p_c = - wnorm(d, g)^2 / sum(abs2, J*g) .* g
+        p_c = - wnorm(d, g)^2 / sum(abs2, J*vec(g)) .* vec(g)
 
         if wnorm(d, p_c) >= delta
             # Cauchy point is out of the region, take the largest step along
@@ -68,13 +68,13 @@ function dogleg!{T}(p::AbstractVector{T}, r::AbstractVector{T}, d::AbstractVecto
             # from this point on we will only need p_i in the term p_i-p_c.
             # so we reuse the vector p_i by computing p_i = p_i - p_c and then
             # just so we aren't confused we name that p_diff
-            p_i .-= p_c
+            p_i .-= vec(p_c)
             p_diff = p_i
 
             # Compute the optimal point on dogleg path
             b = 2 * wdot(d, p_c, d, p_diff)
             a = wnorm(d, p_diff)^2
-            tau = (-b+sqrt(b^2 - 4a*(wnorm(d, p_c)^2 - delta^2)))/(2a)
+            tau = (-b + sqrt(b^2 - 4a*(wnorm(d, p_c)^2 - delta^2)))/(2a)
             p_c .+= tau .* p_diff
             copy!(p, p_c)
         end
@@ -91,10 +91,10 @@ function trust_region_{T}(df::OnceDifferentiable,
                           extended_trace::Bool,
                           factor::T,
                           autoscale::Bool)
-    x = vec(copy(initial_x)) # Current point
+    x = copy(initial_x) # Current point
     nn = length(x)
     xold = similar(x) # Old point
-    r = similar(x)          # Current residual
+    r = similar(df.F)       # Current residual
     r_predict = similar(x)  # predicted residual
     p = similar(x)          # Step
     d = similar(x)          # Scaling vector
@@ -142,9 +142,9 @@ function trust_region_{T}(df::OnceDifferentiable,
         value!(df, x)
 
         # Ratio of actual to predicted reduction (equation 11.47 in N&W)
-        A_mul_B!(r_predict, jacobian(df), p)
+        A_mul_B!(vec(r_predict), jacobian(df), vec(p))
         r_predict .+= r
-        rho = (sum(abs2,r) - sum(abs2, value(df))) / (sum(abs2,r) - sum(abs2,r_predict))
+        rho = (sum(abs2, r) - sum(abs2, value(df))) / (sum(abs2, r) - sum(abs2, r_predict))
 
         if rho > eta
             # Successful iteration
@@ -181,7 +181,7 @@ function trust_region_{T}(df::OnceDifferentiable,
         name *= " and autoscaling"
     end
     return SolverResults(name,
-                         initial_x, reshape(x, size(initial_x)...), norm(r, Inf),
+                         initial_x, reshape(x, size(initial_x)...), vecnorm(r, Inf),
                          it, x_converged, xtol, f_converged, ftol, tr,
                          first(df.f_calls), first(df.df_calls))
 end

test/runtests.jl

@@ -1,6 +1,7 @@
 using NLsolve
 using DiffEqBase
 using Base.Test
+import Base.convert
 
 add_jl(x) = endswith(x, ".jl") ? x : x*".jl"
 
@@ -9,6 +10,7 @@ type WrappedArray{T,N} <: DEDataArray{T,N}
 end
 
 Base.reshape(A::WrappedArray, dims::Int...) = WrappedArray(reshape(A.x, dims...))
+Base.convert(A::Type{WrappedArray{T,N}}, B::Array{T,N}) where {T,N} = WrappedArray(copy(B))
 
 if length(ARGS) > 0
     tests = map(add_jl, ARGS)