From 0bb85a55da0184ed846bb0fa41f9a0d29432b797 Mon Sep 17 00:00:00 2001 From: Patrick Kofod Mogensen Date: Wed, 13 Dec 2017 08:48:59 +0100 Subject: [PATCH] Fix array input. --- src/differentiable_vectors/autodiff.jl | 35 ++++++++--------- src/nlsolve/utils.jl | 8 ++-- src/solvers/mcp.jl | 52 +++++++++++++------------- src/solvers/mcp_func_defs.jl | 6 +-- src/solvers/newton.jl | 16 ++++---- src/solvers/trust_region.jl | 24 ++++++------ test/runtests.jl | 2 + 7 files changed, 69 insertions(+), 74 deletions(-) diff --git a/src/differentiable_vectors/autodiff.jl b/src/differentiable_vectors/autodiff.jl index 7b407bd..a80598e 100644 --- a/src/differentiable_vectors/autodiff.jl +++ b/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 diff --git a/src/nlsolve/utils.jl b/src/nlsolve/utils.jl index b21cf74..cf89f9b 100644 --- a/src/nlsolve/utils.jl +++ b/src/nlsolve/utils.jl @@ -1,5 +1,5 @@ -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] @@ -7,7 +7,7 @@ function wdot{T}(wx::AbstractVector{T}, x::AbstractVector{T}, 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)) diff --git a/src/solvers/mcp.jl b/src/solvers/mcp.jl index c38eb4b..e476eb1 100644 --- a/src/solvers/mcp.jl +++ b/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 diff --git a/src/solvers/mcp_func_defs.jl b/src/solvers/mcp_func_defs.jl index efa9894..3137b77 100644 --- a/src/solvers/mcp_func_defs.jl +++ b/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, diff --git a/src/solvers/newton.jl b/src/solvers/newton.jl index 4579641..764296b 100644 --- a/src/solvers/newton.jl +++ b/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,8 +96,8 @@ 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) @@ -105,7 +105,7 @@ function newton_{T}(df::OnceDifferentiable, # 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 diff --git a/src/solvers/trust_region.jl b/src/solvers/trust_region.jl index 2a9e8a5..a562bb4 100644 --- a/src/solvers/trust_region.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index b48d2ba..185741e 100644 --- a/test/runtests.jl +++ b/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)