From 925c5be1e7635a88297234d6c0f8b20e6553977a Mon Sep 17 00:00:00 2001 From: Kanav Gupta Date: Fri, 10 Jul 2020 20:46:51 +0530 Subject: [PATCH 1/3] (fix) hcat non-sparse matrices --- src/jacobians.jl | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/jacobians.jl b/src/jacobians.jl index f85de1e..d4cdc88 100644 --- a/src/jacobians.jl +++ b/src/jacobians.jl @@ -171,7 +171,8 @@ function finite_difference_jacobian( end vecx = _vec(x) vecx1 = _vec(x1) - J = jac_prototype isa Nothing ? (sparsity isa Nothing ? false.*vecfx.*vecx' : zeros(eltype(x),size(sparsity))) : zero(jac_prototype) + J = jac_prototype isa Nothing ? (sparsity isa Nothing ? Array{eltype(x),2}(undef, length(vecfx), 0) : zeros(eltype(x),size(sparsity))) : zero(jac_prototype) + @show J nrows, ncols = size(J) if !(sparsity isa Nothing) @@ -189,7 +190,11 @@ function finite_difference_jacobian( _x1 = reshape(_vecx1,size(x)) vecfx1 = _vec(f(_x1)) dx = (vecfx1-vecfx)/epsilon - J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + if jac_prototype isa Nothing + J = hcat(J, dx) + else + J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + end else tmp = norm(vecx .* (colorvec .== color_i)) epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) @@ -214,7 +219,11 @@ function finite_difference_jacobian( vecfx1 = _vec(f(_x1)) vecfx = _vec(f(_x)) dx = (vecfx1-vecfx)/(2epsilon) - J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + if jac_prototype isa Nothing + J = hcat(J, dx) + else + J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + end else tmp = norm(vecx1 .* (colorvec .== color_i)) epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) @@ -238,7 +247,11 @@ function finite_difference_jacobian( _x = reshape(_vecx,size(x)) vecfx = _vec(f(_x)) dx = imag(vecfx)/epsilon - J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + if jac_prototype isa Nothing + J = hcat(J, dx) + else + J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + end else _vecx = @. vecx + im * epsilon * (colorvec == color_i) _x = reshape(_vecx,size(x)) From e9ad65cab3c209ee5c8a0f5e55711bb535a416e1 Mon Sep 17 00:00:00 2001 From: Kanav Gupta Date: Sat, 11 Jul 2020 16:03:11 +0530 Subject: [PATCH 2/3] use mapreduce --- src/jacobians.jl | 147 ++++++++++++++++++++++++++--------------------- 1 file changed, 82 insertions(+), 65 deletions(-) diff --git a/src/jacobians.jl b/src/jacobians.jl index d4cdc88..4a24760 100644 --- a/src/jacobians.jl +++ b/src/jacobians.jl @@ -172,7 +172,6 @@ function finite_difference_jacobian( vecx = _vec(x) vecx1 = _vec(x1) J = jac_prototype isa Nothing ? (sparsity isa Nothing ? Array{eltype(x),2}(undef, length(vecfx), 0) : zeros(eltype(x),size(sparsity))) : zero(jac_prototype) - @show J nrows, ncols = size(J) if !(sparsity isa Nothing) @@ -182,83 +181,101 @@ function finite_difference_jacobian( end if fdtype == Val{:forward} - @inbounds for color_i ∈ 1:maximum(colorvec) - if sparsity isa Nothing - x_save = ArrayInterface.allowed_getindex(vecx,color_i) - epsilon = compute_epsilon(Val{:forward}, x_save, relstep, absstep, dir) - _vecx1 = Base.setindex(vecx,x_save+epsilon,color_i) - _x1 = reshape(_vecx1,size(x)) - vecfx1 = _vec(f(_x1)) - dx = (vecfx1-vecfx)/epsilon - if jac_prototype isa Nothing - J = hcat(J, dx) - else + + function calculate_Ji_forward(i) + x_save = ArrayInterface.allowed_getindex(vecx, i) + epsilon = compute_epsilon(Val{:forward}, x_save, relstep, absstep, dir) + _vecx1 = Base.setindex(vecx, x_save+epsilon, i) + _x1 = reshape(_vecx1,size(x)) + vecfx1 = _vec(f(_x1)) + dx = (vecfx1-vecfx) / epsilon + return dx + end + + if jac_prototype isa Nothing + J = mapreduce(calculate_Ji_forward, hcat, 1:maximum(colorvec)) + else + @inbounds for color_i ∈ 1:maximum(colorvec) + if sparsity isa Nothing + dx = calculate_Ji(color_i) J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + else + tmp = norm(vecx .* (colorvec .== color_i)) + epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) + _vecx = @. vecx + epsilon * (colorvec == color_i) + _x = reshape(_vecx,size(x)) + vecfx1 = _vec(f(_x)) + dx = (vecfx1-vecfx)/epsilon + Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols) + J = J + Ji end - else - tmp = norm(vecx .* (colorvec .== color_i)) - epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) - _vecx = @. vecx + epsilon * (colorvec == color_i) - _x = reshape(_vecx,size(x)) - vecfx1 = _vec(f(_x)) - dx = (vecfx1-vecfx)/epsilon - Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols) - J = J + Ji end end elseif fdtype == Val{:central} - @inbounds for color_i ∈ 1:maximum(colorvec) - if sparsity isa Nothing - x1_save = ArrayInterface.allowed_getindex(vecx1,color_i) - x_save = ArrayInterface.allowed_getindex(vecx,color_i) - epsilon = compute_epsilon(Val{:forward}, x1_save, relstep, absstep, dir) - _vecx1 = Base.setindex(vecx1,x1_save+epsilon,color_i) - _vecx = Base.setindex(vecx,x_save-epsilon,color_i) - _x1 = reshape(_vecx1,size(x)) - _x = reshape(_vecx,size(x)) - vecfx1 = _vec(f(_x1)) - vecfx = _vec(f(_x)) - dx = (vecfx1-vecfx)/(2epsilon) - if jac_prototype isa Nothing - J = hcat(J, dx) - else + + function calculate_Ji_central(i) + x1_save = ArrayInterface.allowed_getindex(vecx1,i) + x_save = ArrayInterface.allowed_getindex(vecx,i) + epsilon = compute_epsilon(Val{:forward}, x1_save, relstep, absstep, dir) + _vecx1 = Base.setindex(vecx1,x1_save+epsilon,i) + _vecx = Base.setindex(vecx,x_save-epsilon,i) + _x1 = reshape(_vecx1,size(x)) + _x = reshape(_vecx,size(x)) + vecfx1 = _vec(f(_x1)) + vecfx = _vec(f(_x)) + dx = (vecfx1-vecfx)/(2epsilon) + return dx + end + + if jac_prototype isa Nothing + J = mapreduce(calculate_Ji_central, hcat, 1:maximum(colorvec)) + else + @inbounds for color_i ∈ 1:maximum(colorvec) + if sparsity isa Nothing + dx = calculate_Ji(color_i) J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + else + tmp = norm(vecx1 .* (colorvec .== color_i)) + epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) + _vecx1 = @. vecx1 + epsilon * (colorvec == color_i) + _vecx = @. vecx - epsilon * (colorvec == color_i) + _x1 = reshape(_vecx1,size(x)) + _x = reshape(_vecx, size(x)) + vecfx1 = _vec(f(_x1)) + vecfx = _vec(f(_x)) + dx = (vecfx1-vecfx)/(2epsilon) + Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols) + J = J + Ji end - else - tmp = norm(vecx1 .* (colorvec .== color_i)) - epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir) - _vecx1 = @. vecx1 + epsilon * (colorvec == color_i) - _vecx = @. vecx - epsilon * (colorvec == color_i) - _x1 = reshape(_vecx1,size(x)) - _x = reshape(_vecx, size(x)) - vecfx1 = _vec(f(_x1)) - vecfx = _vec(f(_x)) - dx = (vecfx1-vecfx)/(2epsilon) - Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols) - J = J + Ji end end elseif fdtype == Val{:complex} && returntype <: Real epsilon = eps(eltype(x)) - @inbounds for color_i ∈ 1:maximum(colorvec) - if sparsity isa Nothing - x_save = ArrayInterface.allowed_getindex(vecx,color_i) - _vecx = Base.setindex(complex.(vecx),x_save+im*epsilon,color_i) - _x = reshape(_vecx,size(x)) - vecfx = _vec(f(_x)) - dx = imag(vecfx)/epsilon - if jac_prototype isa Nothing - J = hcat(J, dx) - else + + function calculate_Ji_complex(i) + x_save = ArrayInterface.allowed_getindex(vecx,i) + _vecx = Base.setindex(complex.(vecx),x_save+im*epsilon,i) + _x = reshape(_vecx,size(x)) + vecfx = _vec(f(_x)) + dx = imag(vecfx)/epsilon + return dx + end + + if jac_prototype isa Nothing + J = mapreduce(calculate_Ji_complex, hcat, 1:maximum(colorvec)) + else + @inbounds for color_i ∈ 1:maximum(colorvec) + if sparsity isa Nothing + dx = calculate_Ji_complex(color_i) J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) + else + _vecx = @. vecx + im * epsilon * (colorvec == color_i) + _x = reshape(_vecx,size(x)) + vecfx = _vec(f(_x)) + dx = imag(vecfx)/epsilon + Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols) + J = J + Ji end - else - _vecx = @. vecx + im * epsilon * (colorvec == color_i) - _x = reshape(_vecx,size(x)) - vecfx = _vec(f(_x)) - dx = imag(vecfx)/epsilon - Ji = _make_Ji(J,rows_index,cols_index,dx,colorvec,color_i,nrows,ncols) - J = J + Ji end end else From e890b3d37f16ff256af2e5755016292a19962325 Mon Sep 17 00:00:00 2001 From: Kanav Gupta Date: Sat, 11 Jul 2020 21:50:38 +0530 Subject: [PATCH 3/3] fix sparsity --- src/jacobians.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/jacobians.jl b/src/jacobians.jl index 4a24760..e4a60e9 100644 --- a/src/jacobians.jl +++ b/src/jacobians.jl @@ -192,7 +192,7 @@ function finite_difference_jacobian( return dx end - if jac_prototype isa Nothing + if jac_prototype isa Nothing && sparsity isa Nothing J = mapreduce(calculate_Ji_forward, hcat, 1:maximum(colorvec)) else @inbounds for color_i ∈ 1:maximum(colorvec) @@ -227,7 +227,7 @@ function finite_difference_jacobian( return dx end - if jac_prototype isa Nothing + if jac_prototype isa Nothing && sparsity isa Nothing J = mapreduce(calculate_Ji_central, hcat, 1:maximum(colorvec)) else @inbounds for color_i ∈ 1:maximum(colorvec) @@ -261,7 +261,7 @@ function finite_difference_jacobian( return dx end - if jac_prototype isa Nothing + if jac_prototype isa Nothing && sparsity isa Nothing J = mapreduce(calculate_Ji_complex, hcat, 1:maximum(colorvec)) else @inbounds for color_i ∈ 1:maximum(colorvec)