WIP

Project.toml

@@ -12,6 +12,7 @@ Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Dolang = "e5c7262c-e9d2-5620-ad8e-1af14eb8a8e3"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0"

src/algos/simul.jl

@@ -32,7 +32,6 @@ function τ(dmodel::Dolo.DYModel{M}, ss::T, a::SVector) where M<:Union{Dolo.YMod
         )
         for j in 1:size(P, 2)
     )
-
     it
 
 end

src/dolo_model.jl

@@ -6,7 +6,10 @@ struct YModel{C,A,B,D,N,S} <: AModel
     source::S
 end 
 
-YModel(N,A,B,C,D) = YModel{typeof(C),typeof(A),typeof(B),typeof(D),N,Nothing}(A,B,C,D,nothing)
+YModel(N,A,B,C,D) = let
+    println("Who is calling?")
+    YModel{typeof(C),typeof(A),typeof(B),typeof(D),N,Nothing}(A,B,C,D,nothing)
+end
 YModel(N,A,B,C,D,S) = YModel{typeof(C),typeof(A),typeof(B),typeof(D),N,typeof(S)}(A,B,C,D,S)
 
 name(::YModel{C,A,B,D,N}) where C where A where B where D where N = N

src/model.jl

@@ -5,6 +5,13 @@ abstract type AbstractDModel end
 const AModel = AbstractModel
 const ADModel = AbstractDModel
 
+macro akin(model)
+    typ = typeof(eval(model))
+    ntyp = Meta.parse(repr(typ))
+    ntyp.args = ntyp.args[1:end-1]
+    return ntyp
+end
+
 # abstract type AModel{A,B,C,N} end
 # abstract type ADModel{A,B,C,D,N} end
 

src/space.jl

@@ -173,7 +173,7 @@ end
 
 
 
-const DEFAULT_GRID_NPOINTS = 15
+const DEFAULT_GRID_NPOINTS = 50
 
 # discretize(space::CSpace; kwargs...)  = discretize(space, DEFAULT_GRID_NPOINTS)
 # discretize(space::CSpace{d}, n::Int) where d = discretize(space, tuple( (n for i=1:d)...))

src/splines/interp.jl

@@ -14,6 +14,18 @@ getindex(A::Vector{Float64}, i::VectorizationBase.Vec{4,Int64}) = VectorizationB
     return sum( xx .* v[:])
 end
 
+matextract(v::AbstractArray{T,3}, i,j,k) where T = SArray{Tuple{2, 2, 2}, T, 3, 8}(
+    v[i,  j  ,k],
+    v[i+1,j  ,k],
+    v[i  ,j+1,k],
+    v[i+1,j+1,k],
+    v[i  ,j  ,k+1],
+    v[i+1,j  ,k+1],
+    v[i  ,j+1,k+1],
+    v[i+1,j+1,k+1]
+
+)
+
 function interp(ranges::NTuple{d, Tuple{Float64, Float64, Int64}}, values::AbstractArray{T,d}, x::SVector{d, U}) where d where T where U
 
     a = SVector( (e[1] for e in ranges)... )
@@ -30,9 +42,10 @@ function interp(ranges::NTuple{d, Tuple{Float64, Float64, Int64}}, values::Abstr
 
     i_ = floor.(Int, i) .+ 1
 
-    inds = (SVector(e,e+1) for e in i_)
+    # inds = tuple( (SVector(e,e+1) for e in i_)...)
+    # return reduce_tensors( λ, values[inds...] )
 
-    return reduce_tensors( λ, values[inds...] )
+    return reduce_tensors( λ, matextract(values, i_...) )
 
 end
 

src/splines/splines.jl

@@ -104,6 +104,7 @@ end
 
     function SplineInterpolator(ranges::Tuple; values=nothing, k=3)
 
+
         n = [e[3] for e in ranges]
         dims = tuple((i+k-1 for i in n)...)
         θ_ = zeros(eltype(values), dims...)
@@ -121,7 +122,14 @@ end
 
     end
 
-    function (spl::SplineInterpolator{G,C,3})(x) where G where C
+    function (spl::CubicInterpolator{G,C})(x::SVector) where G where C
+        a = tuple( (e[1] for e in spl.grid)...)
+        b = tuple( (e[2] for e in spl.grid)...)
+        n = tuple( (e[3] for e in spl.grid)...)
+        splines.eval_UC_spline(a,b,n, spl.θ, x)
+    end
+
+    function (spl::SplineInterpolator{G,C,3})(x::SVector) where G where C
         a = tuple( (e[1] for e in spl.grid)...)
         b = tuple( (e[2] for e in spl.grid)...)
         n = tuple( (e[3] for e in spl.grid)...)