ENH: replace evaluate function with methods for call

src/algos/time_iteration.jl

@@ -12,7 +12,7 @@ function residual(model, dprocess, s, x::Array{Array{Float64,2},1}, p, dr)
             for n=1:N
                 S[n,:] = Dolo.transition(model, m, s[n,:], x[i][n,:], M, p)
             end
-            X = evaluate(dr, i, j, S)
+            X = dr(i, j, S)
             for n=1:N
                 res[i][n,:] += w*Dolo.arbitrage(model, m, s[n,:], x[i][n,:], M, S[n,:], X[n,:], p)
             end
@@ -61,7 +61,7 @@ function time_iteration(model, process, init_dr; verbose=true, maxit=100, tol=1e
 
     p = model.calibration[:parameters] :: Vector{Float64}
 
-    x0 = [evaluate(init_dr, i, endo_nodes) for i=1:nsd]
+    x0 = [init_dr(i, endo_nodes) for i=1:nsd]
 
 
     x_lb = Array{Float64,2}[cat(1,[Dolo.controls_lb(model,node(dprocess,i) ,endo_nodes[n,:],p)' for n=1:N]...) for i=1:nsd]

src/algos/time_iteration_direct.jl

@@ -19,7 +19,7 @@ function time_iteration_direct(model, process, init_dr; verbose=true, maxit=100,
     return cat(1,x...)
   end
   # initial guess for controls
-  x0 = [evaluate(init_dr, i, endo_nodes) for i=1:nsd]
+  x0 = [init_dr(i, endo_nodes) for i=1:nsd]
   # set the bound for the controls to check during the iterations not to violate them
   x_lb = Array{Float64,2}[cat(1,[Dolo.controls_lb(model,node(dprocess,i) ,endo_nodes[n,:],p)' for n=1:N]...) for i=1:nsd]
   x_ub = Array{Float64,2}[cat(1,[Dolo.controls_ub(model,node(dprocess,i),endo_nodes[n,:],p)' for n=1:N]...) for i=1:nsd]
@@ -60,7 +60,7 @@ function time_iteration_direct(model, process, init_dr; verbose=true, maxit=100,
                 S[n,:] = Dolo.transition(model, m, s[n,:], x0[i][n,:], M, p)
             end
             # interpolate controles conditional states of tomorrow
-            X = evaluate(dr, i, j, S)
+            X = dr(i, j, S)
             # Compute expectations as a weited average of the exo states w_j
             for n=1:N
                 E_f[i][n,:] += w*Dolo.expectation(model, M, S[n,:], X[n,:], p)

src/algos/value_iteration.jl

@@ -35,7 +35,7 @@ function evaluate_policy(model, dr; verbose=true, maxit=100, )
     S =  copy(endo_nodes)
 
     # Controls at time t
-    x= [evaluate(dr, i, endo_nodes) for i=1:number_of_smooth_drs(dprocess)]
+    x= [dr(i, endo_nodes) for i=1:number_of_smooth_drs(dprocess)]
     s = deepcopy(endo_nodes)
     for i=1:size(res,1)
         m = node(dprocess,i)  ::Vector{Float64}
@@ -70,7 +70,7 @@ function evaluate_policy(model, dr; verbose=true, maxit=100, )
                  end
                  # Update value function
                  for n=1:N
-                     E_V[i][n,1] += w*evaluate(drv,i, j, S[n,:])[1]
+                     E_V[i][n,1] += w*drv(i, j, S[n,:])[1]
                  end
             end
         end
@@ -98,7 +98,7 @@ function update_value(model, β::Float64, dprocess, drv, i, s::Vector{Float64},
         M = inodes(dprocess,i,j) ::Vector{Float64}
         w = iweights(dprocess,i,j) ::Float64
         S = Dolo.transition(model, m, s, x0, M, p)
-        E_V += w*evaluate(drv, i, j, S)[1]
+        E_V += w*drv(i, j, S)[1]
     end
     u = Dolo.felicity(model, m, s, x0, p)[1]
     E_V = u + β.*E_V
@@ -118,7 +118,7 @@ function solve_policy(model, dr; verbose=true, maxit=5000, )
 
     # compute the value function
     absmax(x) = max([maximum(abs(x[i])) for i=1:length(x)]...)
-    p = model.calibration[:parameters] :: Vector{Float64}
+    p = model.calibration[:parameters]
 
     endo_nodes = nodes(grid)
     # Number of endogenous nodes
@@ -143,7 +143,7 @@ function solve_policy(model, dr; verbose=true, maxit=5000, )
     S =  copy(endo_nodes)
 
     # Controls at time t
-    x= [evaluate(dr, i, endo_nodes) for i=1:number_of_smooth_drs(dprocess)]
+    x = [dr(i, endo_nodes) for i=1:number_of_smooth_drs(dprocess)]
     x0 = deepcopy(x)
 
     if verbose
@@ -154,13 +154,13 @@ function solve_policy(model, dr; verbose=true, maxit=5000, )
       println("Evaluating initial policy (done)")
     end
 
-    v0 = [evaluate(drv, i, endo_nodes) for i=1:number_of_smooth_drs(dprocess)]
+    v0 = [drv(i, endo_nodes) for i=1:number_of_smooth_drs(dprocess)]
     v = deepcopy(v0)
     #Preparation for a loop
     tol = 1e-6
-    err=10
+    err = 10
     err_x = 10
-    Err=zeros(maxit)
+    Err = zeros(maxit)
     it = 0
 
     n_eval = 50
@@ -193,7 +193,7 @@ function solve_policy(model, dr; verbose=true, maxit=5000, )
                     lower[lower.<-1000000] = -1000000.0
                     initial_x = x0[i][n,:]
                     # try
-                    results = optimize(DifferentiableFunction(fobj), initial_x, lower, upper, Fminbox(), optimizer = NelderMead)
+                    @time results = optimize(DifferentiableFunction(fobj), initial_x, lower, upper, Fminbox(), optimizer = NelderMead)
                     # results = optimize(DifferentiableFunction(fobj), initial_x, NelderMead())
                     # println(results)
                     xn = Optim.minimizer(results)

src/numeric/decision_rules.jl

@@ -21,10 +21,10 @@ end
 type ConstantDecisionRule <: AbstractDecisionRule
     values::Array{Float64,1}
 end
-evaluate(dr::ConstantDecisionRule, x::Array{Float64,1}) = dr.values
-evaluate(dr::ConstantDecisionRule, x::Array{Float64,2}) = repmat(dr.values',size(x,1),1)
-evaluate(dr::ConstantDecisionRule, i::Int, x::Union{Vector,Matrix}) =  evaluate(dr, x)
-evaluate(dr::ConstantDecisionRule, i::Int, j::Int, x::Union{Vector, Matrix}) = evaluate(dr, x)
+(dr::ConstantDecisionRule)(x::Array{Float64,1}) = dr.values
+(dr::ConstantDecisionRule)(x::Array{Float64,2}) = repmat(dr.values',size(x,1),1)
+(dr::ConstantDecisionRule)(i::Int, x::Union{Vector,Matrix}) =  dr(x)
+(dr::ConstantDecisionRule)(i::Int, j::Int, x::Union{Vector, Matrix}) = dr(x)
 
 
 
@@ -63,7 +63,7 @@ function DecisionRule(mc::DiscreteMarkovProcess, grid::CartesianGrid, values::Ar
     return dr
 end
 
-function evaluate(dr::MCDecisionRule, i::Int, x::Array{Float64,2})
+function (dr::MCDecisionRule)(i::Int, x::Array{Float64,2})
     a = dr.grid.min
     b = dr.grid.max
     n = dr.grid.n
@@ -72,8 +72,8 @@ function evaluate(dr::MCDecisionRule, i::Int, x::Array{Float64,2})
     return res
 end
 
-evaluate(dr::MCDecisionRule, i::Int, x::Array{Float64,1}) =  evaluate(dr, i, x')[:]
-evaluate(dr::MCDecisionRule, i::Int, j::Int, x::Union{Vector, Matrix}) = evaluate(dr, j, x)
+(dr::MCDecisionRule)(i::Int, x::Array{Float64,1}) =  dr(i, x')[:]
+(dr::MCDecisionRule)(i::Int, j::Int, x::Union{Vector, Matrix}) = dr(j, x)
 
 
 # Decision Rule on IID Process
@@ -117,18 +117,18 @@ function DecisionRule(proc::MvNormal, grid::CartesianGrid, values::Array{Float64
     return DecisionRule(proc, grid, [values])
 end
 
-function evaluate(dr::IDecisionRule, x::Array{Float64,2})
+function (dr::IDecisionRule)(x::Array{Float64,2})
     a = dr.grid.min
     b = dr.grid.max
     n = dr.grid.n
     cc = dr.coefficients[1]
     res = splines.eval_UC_multi_spline(a,b,n,cc,x)'
     return res
 end
-evaluate(dr::IDecisionRule, x::Array{Float64,1}) = evaluate(dr, x')[:]
-evaluate(dr::IDecisionRule, i::Int, x::Array{Float64,2}) = evaluate(dr, x)
-evaluate(dr::IDecisionRule, i::Int, x::Array{Float64,1}) =  evaluate(dr, x)
-evaluate(dr::IDecisionRule, i::Int, j::Int, x::Union{Vector, Matrix}) = evaluate(dr, x)
+(dr::IDecisionRule)(x::Array{Float64,1}) = dr(x')[:]
+(dr::IDecisionRule)(i::Int, x::Array{Float64,2}) = dr(x)
+(dr::IDecisionRule)(i::Int, x::Array{Float64,1}) =  dr(x)
+(dr::IDecisionRule)(i::Int, j::Int, x::Union{Vector, Matrix}) = dr(x)
 
 
 
@@ -191,7 +191,7 @@ function DecisionRule(dp::DiscretizedProcess, grid::CartesianGrid, values::Array
     return dr
 end
 
-function evaluate(dr::CPDecisionRule, z::Array{Float64,2})
+function (dr::CPDecisionRule)(z::Array{Float64,2})
     a = dr.grid.min
     b = dr.grid.max
     n = dr.grid.n
@@ -200,10 +200,10 @@ function evaluate(dr::CPDecisionRule, z::Array{Float64,2})
     return res
 end
 
-evaluate(dr::CPDecisionRule, z::Vector) = evaluate(dr,z')[:]
+(dr::CPDecisionRule)(z::Vector) = dr(z')[:]
 
 
-function evaluate(dr::CPDecisionRule, i::Int64, x::Array{Float64,2})
+function (dr::CPDecisionRule)(i::Int64, x::Array{Float64,2})
     a = dr.grid.min
     b = dr.grid.max
     n = dr.grid.n
@@ -217,7 +217,7 @@ end
 
 
 
-function evaluate(dr::CPDecisionRule, i::Int64, j::Int64, x::Array{Float64,2})
+function (dr::CPDecisionRule)(i::Int64, j::Int64, x::Array{Float64,2})
     a = dr.grid.min
     b = dr.grid.max
     n = dr.grid.n
@@ -231,5 +231,5 @@ end
 
 
 
-evaluate(dr::CPDecisionRule, i::Int64, x::Vector) = evaluate(dr,i,x')[:]
-evaluate(dr::CPDecisionRule, i::Int64, j::Int64, x::Vector) = evaluate(dr,i,j,x')[:]
+(dr::CPDecisionRule)(i::Int64, x::Vector) = dr(i,x')[:]
+(dr::CPDecisionRule)(i::Int64, j::Int64, x::Vector) = dr(i,j,x')[:]

src/numeric/simulations.jl

@@ -53,20 +53,20 @@ function QuantEcon.simulate(m::AbstractNumericModel, dr::AbstractDecisionRule,
             epsilons = rand(ϵ_dist, n_exp)'
         end
 
-        s = slice(s_simul, :, :, t)
+        s = view(s_simul, :, :, t)
 
         # extract x and then fill it
-        x = slice(x_simul, :, :, t)
+        x = view(x_simul, :, :, t)
         dr(s, x)
 
         if has_aux
             # extract a and then fill it
-            y = slice(y_simul, :, :, t)
+            y = view(y_simul, :, :, t)
             evaluate!(a, s, x, params, y)
         end
 
         if t < horizon
-            ss = slice(s_simul, :, :, t+1)
+            ss = view(s_simul, :, :, t+1)
             evaluate!(g, s, x, epsilons, params, ss)
         end
     end

test/test_algos.jl

@@ -1,11 +1,9 @@
 path = Pkg.dir("Dolo")
 
-Pkg.build("QuantEcon")
 import Dolo
 
-filename = joinpath(path,"examples","models","rbc_dtcc_mc.yaml")
-# filename = joinpath(path,"examples","models","sudden_stop.yaml")
-model_mc = Dolo.yaml_import(filename)
+fn = joinpath(path,"examples","models","rbc_dtcc_mc.yaml")
+model_mc = Dolo.yaml_import(fn)
 
 drc = Dolo.ConstantDecisionRule(model_mc.calibration[:controls])
 @time dr0, drv0 = Dolo.solve_policy(model_mc, drc, verbose=true, maxit=10000 )
@@ -18,17 +16,17 @@ drc = Dolo.ConstantDecisionRule(model_mc.calibration[:controls])
 
 # compare with prerecorded values
 kvec = linspace(dr.grid.min[1],dr.grid.max[1],10)
-nvec = [Dolo.evaluate(dr,1,[k])[1] for k in kvec]
-ivec = [Dolo.evaluate(dr,1,[k])[2] for k in kvec]
+nvec = [dr(1,[k])[1] for k in kvec]
+ivec = [dr(1,[k])[2] for k in kvec]
 # compare  time_iteration_direct
-nvec_d = [Dolo.evaluate(drd,1,[k])[1] for k in kvec]
-ivec_d = [Dolo.evaluate(drd,1,[k])[2] for k in kvec]
-@assert maximum(abs(nvec_d-nvec))<1e-
+nvec_d = [drd(1,[k])[1] for k in kvec]
+ivec_d = [drd(1,[k])[2] for k in kvec]
+@assert maxabs(nvec_d-nvec)<1e-4
 
 # compare  vfi
-nvec_0 = [Dolo.evaluate(dr0,1,[k])[1] for k in kvec]
-ivec_0 = [Dolo.evaluate(dr0,1,[k])[2] for k in kvec]
-@assert maximum(abs(nvec_0-nvec))<1e-4
+nvec_0 = [dr0(1,[k])[1] for k in kvec]
+ivec_0 = [dr0(1,[k])[2] for k in kvec]
+@assert maxabs(nvec_0-nvec)<1e-4
 
 
 # let's redo when model is stable !
@@ -39,8 +37,8 @@ ivec_0 = [Dolo.evaluate(dr0,1,[k])[2] for k in kvec]
 
 
 # this one needs a lower value of beta or a better initial guess
-filename = joinpath(path,"examples","models","rbc_dtcc_iid.yaml")
-model = Dolo.yaml_import(filename)
+fn = joinpath(path,"examples","models","rbc_dtcc_iid.yaml")
+model = Dolo.yaml_import(fn)
 
 drc = Dolo.ConstantDecisionRule(model.calibration[:controls])
 @time dr0, drv0 = Dolo.solve_policy(model, drc, verbose=true, maxit=1000 )
@@ -52,22 +50,22 @@ drc = Dolo.ConstantDecisionRule(model.calibration[:controls])
 @time drv = Dolo.evaluate_policy(model, dr, verbose=true)
 
 kvec = linspace(dr.grid.min[1],dr.grid.max[1],10)
-nvec = [Dolo.evaluate(dr,1,[k])[1] for k in kvec]
-ivec = [Dolo.evaluate(dr,1,[k])[2] for k in kvec]
-nvec_d = [Dolo.evaluate(drd,1,[k])[1] for k in kvec]
-ivec_d = [Dolo.evaluate(drd,1,[k])[2] for k in kvec]
-nvec_0 = [Dolo.evaluate(dr0,1,[k])[1] for k in kvec]
-ivec_0 = [Dolo.evaluate(dr0,1,[k])[2] for k in kvec]
+nvec = [dr(1,[k])[1] for k in kvec]
+ivec = [dr(1,[k])[2] for k in kvec]
+nvec_d = [drd(1,[k])[1] for k in kvec]
+ivec_d = [drd(1,[k])[2] for k in kvec]
+nvec_0 = [dr0(1,[k])[1] for k in kvec]
+ivec_0 = [dr0(1,[k])[2] for k in kvec]
 
-@assert maximum(abs(nvec_d-nvec))<1e-5
-@assert maximum(abs(nvec_0-nvec))<1e-5 # not satisfied right now (see tol. of optimizer)
+@assert maxabs(nvec_d-nvec)<1e-5
+@assert maxabs(nvec_0-nvec)<1e-5 # not satisfied right now (see tol. of optimizer)
 
 
 
 
 # does not work yet
-filename = joinpath(path,"examples","models","rbc_dtcc_ar1.yaml")
-model = Dolo.yaml_import(filename)
+fn = joinpath(path,"examples","models","rbc_dtcc_ar1.yaml")
+model = Dolo.yaml_import(fn)
 
 Dolo.discretize(model.exogenous)
 

test/test_decision_rules.jl

@@ -14,9 +14,9 @@ vv = [values for i=1:size(mc.values,1)]
 dr = DecisionRule(mc, grid, vv)
 set_values(dr,vv) # filter coefficients
 
-evaluate(dr, 1, [0.1 0.5])
+dr(1, [0.1 0.5])
 s0 = [0.2, 0.5]'
-res =  evaluate(dr, 1, s0)
+res =  dr(1, s0)
 
 
 grid = CartesianGrid([0.0,1.0], [0.1, 0.4], [5,4])
@@ -27,5 +27,5 @@ vv = [values]
 dr = DecisionRule(mvn, grid, vv)
 set_values(dr, values) # filter coefficients
 s0 = [0.2, 0.5]'
-res =  evaluate(dr, 1, s0)
-evaluate(dr, 1, [0.1 0.5])
+res =  dr(1, s0)
+dr(1, [0.1 0.5])

test/test_time_iteration.jl

@@ -28,8 +28,8 @@ steady_state(model, model.calibration)
 
 @time dr = time_iteration(model, process, verbose=true)
 
-ivals = [evaluate(dr, 1, [i])[1] for i=kvec ]
-nvals = [evaluate(dr, 1, [i])[2] for i=kvec ]
+ivals = [dr(1, [i])[1] for i=kvec ]
+nvals = [dr(1, [i])[2] for i=kvec ]
 
 
 using Gadfly
@@ -38,5 +38,5 @@ plot(x=kvec, y=ivals)
 
 plot(x=kvec, y=nvals)
 
-kkvec = [evaluate(dr,1,[k])[1] for k in kvec]
+kkvec = [dr(1,[k])[1] for k in kvec]
 plot(x=kvec, y=yvec, Geom.line)

test/test_value_iteration.jl

@@ -43,5 +43,5 @@ using Gadfly
 kvec = linspace(8,12,100)
 
 drv_answer.grid
-yvec = [evaluate(drv_answer,1,[k])[1] for k in kvec]
+yvec = [drv_answer(1,[k])[1] for k in kvec]
 plot(x=kvec, y=yvec, Geom.line)