Merge branch 'master' of https://github.com/EconForge/Dolo.jl

check_1d.jl

@@ -0,0 +1,40 @@
+using Dolo
+using Test
+
+d = 2
+
+NN = 20
+
+vars =  tuple( (Symbol("d$i") for i in 1:d)... )
+sizes = tuple( (NN+i+1 for i in 1:d)...)
+args = tuple( ( (0.0, 1.0*i, sizes[i]) for i in 1:d)... )
+cg = Dolo.CGrid( args )
+
+@test isbits(cg)
+
+@test Dolo.size(cg) == sizes
+
+
+# check that iterator for is non-allocating
+f(cg,v) = begin e = sum( sum( e for e in cg ) ); v[1] = e ; nothing end
+res = zeros(1)
+f(cg, res)
+@test 0==(@allocated f(cg, res))
+
+# Dolo.enum
+# returns QP object
+# ...
+
+# check that enum for is non-allocating
+function fun(cg::Dolo.CGrid{d}) where d
+    if sum( sum( (sum(loc)*val) for (;loc, val) in Dolo.enum(cg) ) ) < -Inf
+        print("")
+    end
+end
+
+fun(cg)
+@test 0==(@allocated fun(cg))
+@time fun(cg)
+
+
+println((@allocated fun(cg)))

misc.jl

@@ -0,0 +1,45 @@
+using Dolo
+
+root_dir = pkgdir(Dolo)
+model = include("$(root_dir)/examples/ymodels/rbc_iid.jl")
+
+isbits(model)
+
+
+Dolo.time_iteration(model; verbose=false);
+Dolo.time_iteration(model; verbose=false, improve=true, improve_K=1000);
+@time Dolo.time_iteration(model; verbose=false, improve=true, improve_K=1000, interpolation=:cubic);
+
+
+
+dmodel = Dolo.discretize(model)
+
+wksp = Dolo.time_iteration_workspace(dmodel; interp_mode=:linear);
+res = Dolo.time_iteration(dmodel, wksp; verbose=false, improve=false, improve_K=1000, engine=:default);
+
+
+wksp = Dolo.time_iteration_workspace(dmodel; interp_mode=:linear);
+res = Dolo.time_iteration(dmodel, wksp; verbose=false, improve=false, improve_K=1000, engine=:default, trace=true);
+
+
+
+@time Dolo.time_iteration(dmodel, wksp; verbose=false, improve=false, improve_K=1000, engine=:default);
+
+@time Dolo.time_iteration(dmodel, wksp; verbose=false, improve=false, improve_K=1000);
+
+
+@time Dolo.time_iteration(dmodel; verbose=true, improve=true, improve_wait=0, improve_K=500);
+
+@time Dolo.time_iteration(dmodel; verbose=true, improve=true, improve_wait=0, improve_K=500, interpolation=:cubic);
+
+
+
+
+wksp = Dolo.time_iteration_workspace(dmodel; interp_mode=:cubic);
+Dolo.time_iteration(dmodel, wksp; verbose=false, improve=false, improve_K=1000);
+@time Dolo.time_iteration(dmodel, wksp; verbose=false, improve=false, improve_K=1000);
+
+
+@time Dolo.time_iteration(dmodel; verbose=true, improve=true, improve_wait=0, improve_K=500);
+@time Dolo.time_iteration(dmodel; verbose=true, improve=true, improve_wait=0, improve_K=500);
+

self_contained.jl

@@ -0,0 +1,110 @@
+using Test
+using StaticArrays
+import Base: getindex, length
+
+struct CGrid{d}
+    ranges::NTuple{d, Tuple{Float64, Float64, Int64}}
+end
+
+
+getindex(g::CGrid{d}, inds::Vararg{Int64,d}) where d = SVector{d}(
+    ( 
+        ( g.ranges[i][1] + (g.ranges[i][2]-g.ranges[i][1])*( (inds[i]-1)/(g.ranges[i][3]-1)) )
+        for i=1:d
+    )
+)
+
+
+length(g::CGrid{1})= g.ranges[1][3]
+enum(g::CGrid{1})= ( (;loc=(i,), val=g[i]) for i=1:length(g))
+enum(g::CGrid{d}) where d = ((;loc=c, val=g[c...]) for c in Iterators.product( tuple( ((1:r[3]) for r in g.ranges)... )... ) )
+
+
+d = 3
+
+NN = 20
+
+vars =  tuple( (Symbol("d$i") for i in 1:d)... )
+sizes = tuple( (NN+i+1 for i in 1:d)...)
+args = tuple( ( (0.0, 1.0*i, sizes[i]) for i in 1:d)... )
+cg = CGrid( args )
+
+cg.ranges
+
+eeenum(cg::CGrid{d}) where d =  (cg[c...] for c in Iterators.product( tuple( ((1:r[3]) for r in cg.ranges)... )... ) )
+
+
+pit = Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}}((1:22, 1:23, 1:24))
+
+v1, s2 = iterate(pit)
+v2, s3 = iterate(pit, s2)
+
+import Base: iterate
+function iterate(cg::CGrid{d}) where d
+
+    tt = tuple( (1:(e[3]) for e in cg.ranges)...)
+    # pit = Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}}(tt)
+    pit = Iterators.ProductIterator(tt)
+    c1, state = iterate(pit)
+
+    return cg[c1...], (pit,state)
+
+end
+
+function iterate(cg::CGrid{d}, sstate) where d
+
+    pit,state = sstate
+
+    res = iterate(pit, state)
+    if res isa Nothing
+        return nothing
+    else
+        c, newstate = res
+        return cg[c...], (pit,newstate)
+    end
+
+end
+
+eltype(cg::CGrid{d}) where d = SVector{d, Float64}
+
+import Base: length, eltype, isdone
+length(cg::CGrid{d}) where d = prod( (e[3] for e in cg.ranges))
+
+function isdone(cg::CGrid{d}, sstate) where d
+    pit, state= sstate
+    return isdone(pit, state)
+end
+
+function check(cg)
+    sum(sum( cg ))
+end
+
+@code_warntype check(cg)
+check(cg)
+@allocations check(cg) == 1
+
+
+
+
+# function fun2(cg::CGrid{d}) where d
+#     if sum( sum( (sum(loc...)*val) for (;loc, val) in enum(cg) ) ) < -Inf
+#         print("")
+#     end
+# end
+
+function fun3(cg::CGrid{d}) where d
+    if sum( sum( (cg[loc...].*val) for (;loc, val) in enum(cg) ) ) < -Inf
+        print("")
+    end
+end
+
+
+# @code_warntype fun(cg)
+eltype(enum(cg))
+
+# fun2(cg)
+# println( @allocated fun2(cg) )
+
+
+fun3(cg)
+println( @allocated fun3(cg) )

test/test_interpolation.jl

@@ -0,0 +1,78 @@
+using Test
+using StaticArrays
+
+@testset verbose=true "Cartesian Interpolation" begin
+
+@testset verbose=true "Interpolation" begin
+
+    for d in 1:3 
+
+        @testset verbose=true "Dimension $d" begin
+
+            NN = 20
+
+            vars =  tuple( (Symbol("d$i") for i in 1:d)... )
+            pairs = (vars[i]=>[0.0, 1.0*i] for i in 1:d)
+            dvars = Dict( pairs )
+            cs = Dolo.CSpace(;
+                pairs...
+            )
+
+            sg = Dolo.discretize(cs)
+
+            coeffs = [(i+1)^2 for i=1:d]
+            coeffs_2 = [(i+1)^3 for i=1:d]
+
+            fl(x) = SVector( sum(x .* coeffs), sum(x.* coeffs_2) )
+            fc(x) = SVector( sum(x .* coeffs + x.^2 .* coeffs), sum(x.* coeffs_2  + x.^3 .* coeffs) )
+
+
+            function testfun(dfun, sg)
+                if sum(sum(dfun(e) for e in sg)) <0.0
+                    println("Oups")
+                else
+                    nothing
+                end
+            end
+
+            @testset verbose=true "Linear" begin
+
+                values = [fl(e) for e in sg]
+                gvec = Dolo.GVector(sg, values)
+
+                dfun = Dolo.DFun(cs, gvec)
+
+                _values = [dfun(e) for e in sg]
+
+                d = _values - values
+                @test isapprox(values, _values)
+                testfun(dfun, sg)
+                @test 0== @allocated testfun(dfun, sg)
+
+            end
+
+            @testset verbose=true "Cubic" begin
+
+                values = [fc(e) for e in sg]
+                gvec = Dolo.GVector(sg, values)
+
+                dfun = Dolo.DFun(cs, gvec; interp_mode=:cubic)
+
+                _values = [dfun(e) for e in sg]
+
+                d = _values - values
+                @test isapprox(values, _values)
+
+                testfun(dfun, sg)
+                @test 0== @allocated testfun(dfun, sg)
+
+            end
+
+        end
+
+    end
+
+end
+
+
+end

test/time_iteration.jl

@@ -0,0 +1,33 @@
+root_dir = pkgdir(Dolo)
+
+@testset "Benchmark Models" verbose=true begin
+
+    model_files = ["rbc_iid", "rbc_mc", "rbc_ar1"]
+
+    for fn in model_files
+
+        @testset "$fn" verbose=true begin
+
+            @testset "Julia" verbose=true begin
+
+                model = include("$(root_dir)/examples/ymodels/$fn.jl")
+                @test isbits(model)
+                sol = Dolo.time_iteration(model; verbose=false)
+
+            end
+
+
+            @testset "YAML" verbose=true begin
+
+                model = include("$(root_dir)/examples/ymodels/$fn.jl")
+                @test isbits(model)
+                sol = Dolo.time_iteration(model; verbose=false)
+
+            end
+
+        end
+
+
+    end
+
+end

test_dev.jl

@@ -0,0 +1,39 @@
+d = 1
+
+NN = 20
+
+vars =  tuple( (Symbol("d$i") for i in 1:d)... )
+sizes = tuple( (NN+i+1 for i in 1:d)...)
+args = tuple( ( (0.0, 1.0*i, sizes[i]) for i in 1:d)... )
+cg = Dolo.CGrid( args )
+
+@test isbits(cg)
+
+@test Dolo.size(cg) == sizes
+
+
+# check that iterator for is non-allocating
+
+
+# check that 'iterate' is non-allocating
+function fun1(cg::Dolo.CGrid{d}) where d
+    if sum( sum( e for e in cg) ) < -Inf
+        print("This should never be ran.")
+    end
+end
+fun1(cg)
+@test 0==(@allocated fun1(cg))
+
+# Dolo.enum
+# returns QP object
+# ...
+
+# check that 'enum' for is non-allocating
+function fun(cg::Dolo.CGrid{d}) where d
+    if sum( sum( sum(q.loc...)*(q.val) for q in Dolo.enum(cg) ) ) < 0.0
+        print("")
+    end
+end
+
+fun(cg)
+@test 0==(@allocated fun(cg))