Skip to content

Commit

Permalink
update blockbandedmatrices
Browse files Browse the repository at this point in the history
  • Loading branch information
@matthieugomez
matthieugomez committed Jun 17, 2024
1 parent 98e0018 commit efea6f4
Show file tree
Hide file tree
Showing 3 changed files with 3 additions and 74 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "EconPDEs"
uuid = "a3315474-fad9-5060-8696-cee5f38a87b7"
version = "1.1.0"
version = "1.1.1"

[deps]
BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
Expand All @@ -13,7 +13,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"

[compat]
BlockBandedMatrices = "0.7, 0.8, 0.9, 0.10, 0.11, 0.12"
BlockBandedMatrices = "0.7, 0.8, 0.9, 0.10, 0.11, 0.12, 0.13"
FiniteDiff = "2"
LinearAlgebra = "1"
NLsolve = "1, 2, 3, 4"
Expand Down
1 change: 1 addition & 0 deletions examples/ConsumptionProblem/WangWangYang.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ function solve!(pts, m, ws, ps)
end
return pts
end

m = WangWangYangModel()
ws = range(m.wmin^(1/2), m.wmax^(1/2), length = 100).^2
ps = 1 .+ stategrid[:w]
Expand Down
72 changes: 0 additions & 72 deletions src/pdesolve.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,75 +174,3 @@ function _setindex!(@nospecialize(a), apm, stategrid::StateGrid, Tsolution, y_M:
end





function pdesolve2(apm, @nospecialize(grid), @nospecialize(yend); is_algebraic = OrderedDict(k => false for k in keys(yend)), bc = nothing, verbose = true, kwargs...)
stategrid = StateGrid(NamedTuple(grid))
S = size(stategrid)
all(size(v) == S for v in values(yend)) || throw(ArgumentError("The length of initial guess (e.g. terminal value) does not equal the length of the state space"))
all(keys(is_algebraic) .== keys(yend)) || throw(ArgumentError("the terminal guess yend and the is_algebric keyword argument must have the same names"))
is_algebraic = OrderedDict(first(p) => fill(last(p), S) for p in pairs(is_algebraic))
y = OrderedDict(first(p) => collect(last(p)) for p in pairs(yend))
# convert to Matrix
yend_M = catlast(values(yend))
is_algebraic_M = catlast(values(is_algebraic))
bc_M = _Array_bc(bc, yend, grid)
Tsolution = Type{tuple(keys(yend)...)}

a = get_a2(apm, stategrid, Tsolution, yend_M, bc_M)
y_M, residual_norm = finiteschemesolve((ydot, y) -> hjb2!(apm, stategrid, Tsolution, ydot, y, bc_M, size(yend_M)), vec(yend_M); is_algebraic = vec(is_algebraic_M), J0c = J0c, verbose = verbose, kwargs... )
y_M = reshape(y_M, size(yend_M)...)
_setindex!(y, y_M)
if a !== nothing
_setindex2!(a, apm, stategrid, Tsolution, y_M, bc_M)
a = merge(y, a)
end
return EconPDEResult(y, residual_norm, a)
end

function get_a2(apm, stategrid::StateGrid, Tsolution, y_M::AbstractArray, bc_M::AbstractArray)
derivatives = differentiate2(Tsolution, stategrid, y_M, bc_M)
result = apm(stategrid, derivatives)
if length(result) == 1
return nothing
else
return OrderedDict(a_key => Array{Float64}(undef, size(stategrid)) for a_key in keys(result[2]))
end
end

# create hjb! that accepts and returns AbstractVector rather than AbstractArrays
function hjb2!(apm, stategrid::StateGrid, Tsolution, ydot::AbstractVector, y::AbstractVector, bc_M::AbstractArray, ysize::NTuple)
y_M = reshape(y, ysize...)
ydot_M = reshape(ydot, ysize...)
vec(hjb!(apm, stategrid, Tsolution, ydot_M, y_M, bc_M))
end

function hjb2!(apm, stategrid::StateGrid, Tsolution, ydot_M::AbstractArray, y_M::AbstractArray, bc_M::AbstractArray)
solution = differentiate2(Tsolution, stategrid, y_M, bc_M)
out = apm(stategrid, solution)
if isa(outi[1], Vector)
_setindex2!(ydot_M, Tsolution, outi, i)
else
_setindex2!(ydot_M, Tsolution, outi[1], i)
end
return ydot_M
end


@generated function _setindex2!(ydot_M::AbstractArray, ::Type{Tsolution}, outi::NamedTuple) where {Tsolution}
N = length(Tsolution.parameters[1])
quote
$(Expr(:meta, :inline))
$(Expr(:block, [Expr(:call, :setindex!, :ydot_M, Expr(:call, :getproperty, :outi, Meta.quot(Symbol(Tsolution.parameters[1][k], :t))), :Colon(), k) for k in 1:N]...))
end
end


function _setindex2!(@nospecialize(a), apm, stategrid::StateGrid, Tsolution, y_M::AbstractArray, bc_M::AbstractArray)
solution = differentiate2(Tsolution, stategrid, y_M, bc_M)
out = apm(stategrid, solution)
for (k, v) in zip(values(a), values(outi))
copyto!(k, v)
end
end

2 comments on commit efea6f4

@matthieugomez
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/109198

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v1.1.1 -m "<description of version>" efea6f4d5fd7e3c7e45c1e4436e8a6ddeb2fa53c
git push origin v1.1.1

Please sign in to comment.