From 96ae0e281f4913005cdc559a62e99f69eab45ee8 Mon Sep 17 00:00:00 2001 From: Winant Pablo Date: Fri, 6 Mar 2020 23:30:25 +0100 Subject: [PATCH] Work on d.r. --- dolo/algos/improved_time_iteration.py | 14 ++-- dolo/algos/time_iteration.py | 8 +-- dolo/algos/value_iteration.py | 4 +- dolo/compiler/model.py | 17 +++-- dolo/compiler/model_numeric.py | 2 +- dolo/compiler/objects.py | 21 +----- dolo/numeric/decision_rule.py | 68 +++++++++---------- dolo/numeric/grids.py | 44 +++++++++++- examples/models/rbc_iid.yaml | 2 +- trash/dolo/algos/dtcscc/accuracy.py | 2 +- trash/dolo/algos/dtcscc/nonlinearsystem.py | 2 +- .../dtcscc/parameterized_expectations.py | 4 +- trash/dolo/algos/dtcscc/time_iteration_2.py | 4 +- 13 files changed, 108 insertions(+), 84 deletions(-) diff --git a/dolo/algos/improved_time_iteration.py b/dolo/algos/improved_time_iteration.py index ec5eeaf2..d89643b1 100644 --- a/dolo/algos/improved_time_iteration.py +++ b/dolo/algos/improved_time_iteration.py @@ -230,7 +230,7 @@ def radius_jac(res,dres,jres,fut_S,dumdr,tol=1e-10,maxit=1000,verbose=False): from .results import AlgoResult, ImprovedTimeIterationResult def improved_time_iteration(model, method='jac', dr0=None, dprocess=None, - interp_type='spline', mu=2, maxbsteps=10, verbose=False, + interp_method='spline', mu=2, maxbsteps=10, verbose=False, tol=1e-8, smaxit=500, maxit=1000, complementarities=True, compute_radius=False, invmethod='iti', details=True): @@ -260,7 +260,7 @@ def vprint(*args, **kwargs): parms = model.calibration['parameters'] dp = dprocess - + if dp is None: dp = model.exogenous.discretize() @@ -268,12 +268,12 @@ def vprint(*args, **kwargs): n_s = len(model.symbols['states']) - grid = model.get_grid() + grid = model.get_endo_grid() - if interp_type == 'spline': - ddr = DecisionRule(dp.grid, grid, dprocess=dp) - ddr_filt = DecisionRule(dp.grid, grid, dprocess=dp) - elif interp_type == 'smolyak': + if interp_method in ('spline', 'linear'): + ddr = DecisionRule(dp.grid, grid, dprocess=dp, interp_method=interp_method) + ddr_filt = DecisionRule(dp.grid, grid, dprocess=dp, interp_method=interp_method) + elif interp_method == 'smolyak': ddr = SmolyakDecisionRule(n_m, grid.min, grid.max, mu) ddr_filt = SmolyakDecisionRule(n_m, grid.min, grid.max, mu) derivative_type = 'numerical' diff --git a/dolo/algos/time_iteration.py b/dolo/algos/time_iteration.py index 8e811df6..0f2d9a3e 100644 --- a/dolo/algos/time_iteration.py +++ b/dolo/algos/time_iteration.py @@ -33,7 +33,7 @@ def residuals_simple(f, g, s, x, dr, dprocess, parms): def time_iteration(model, dr0=None, dprocess=None, with_complementarities=True, verbose=True, grid={}, - maxit=1000, inner_maxit=10, tol=1e-6, hook=None, details=False): + maxit=1000, inner_maxit=10, tol=1e-6, hook=None, details=False, interp_method='cubic'): ''' Finds a global solution for ``model`` using backward time-iteration. @@ -83,11 +83,11 @@ def vprint(t): n_x = len(x0) n_s = len(model.symbols['states']) - endo_grid = model.get_grid(**grid) + endo_grid = model.get_endo_grid(**grid) exo_grid = dprocess.grid - mdr = DecisionRule(exo_grid, endo_grid, dprocess=dprocess) + mdr = DecisionRule(exo_grid, endo_grid, dprocess=dprocess, interp_method=interp_method) grid = mdr.endo_grid.nodes N = grid.shape[0] @@ -105,7 +105,7 @@ def vprint(t): for i_m in range(n_ms): m = dprocess.node(i_m) controls_0[i_m, :, :] = dr0(m, grid) - + f = model.functions['arbitrage'] g = model.functions['transition'] diff --git a/dolo/algos/value_iteration.py b/dolo/algos/value_iteration.py index de8737f7..6816a6b3 100644 --- a/dolo/algos/value_iteration.py +++ b/dolo/algos/value_iteration.py @@ -62,7 +62,7 @@ def value_iteration(model, n_mv = dprocess.n_inodes( 0) # this assume number of integration nodes is constant - endo_grid = model.get_grid(**grid) + endo_grid = model.get_endo_grid(**grid) exo_grid = dprocess.grid @@ -252,7 +252,7 @@ def evaluate_policy(model, n_v = len(v0) n_s = len(model.symbols['states']) - endo_grid = model.get_grid(**grid) + endo_grid = model.get_endo_grid(**grid) exo_grid = dprocess.grid if dr0 is not None: diff --git a/dolo/compiler/model.py b/dolo/compiler/model.py index df834524..5bc52429 100644 --- a/dolo/compiler/model.py +++ b/dolo/compiler/model.py @@ -183,9 +183,7 @@ def get_exogenous(self): return ProductProcess(*exogenous.values()) - def get_grid(self): - - states = self.symbols['states'] + def get_endo_grid(self): # determine bounds: domain = self.get_domain() @@ -204,13 +202,20 @@ def get_grid(self): args = {'min': min, 'max': max} if grid_type.lower() in ('cartesian', 'cartesiangrid'): - # from dolo.numerid.grids import CartesianGrid - from dolo.numeric.grids import CartesianGrid + from dolo.numeric.grids import UniformCartesianGrid orders = get_address(self.data, ["options:grid:n", "options:grid:orders"]) if orders is None: orders = [20] * len(min) - grid = CartesianGrid(min=min, max=max, n=orders) + grid = UniformCartesianGrid(min=min, max=max, n=orders) + elif grid_type.lower() in ('nonuniformcartesian', 'nonuniformcartesiangrid'): + from dolo.compiler.language import eval_data + from dolo.numeric.grids import NonUniformCartesianGrid + calibration = self.get_calibration() + nodes = [eval_data(e, calibration) for e in self.data['options']['grid']] + print(nodes) + # each element of nodes should be a vector + return NonUniformCartesianGrid(nodes) elif grid_type.lower() in ('smolyak', 'smolyakgrid'): from dolo.numeric.grids import SmolyakGrid mu = get_address(self.data, ["options:grid:mu"]) diff --git a/dolo/compiler/model_numeric.py b/dolo/compiler/model_numeric.py index 13be0083..f3350903 100644 --- a/dolo/compiler/model_numeric.py +++ b/dolo/compiler/model_numeric.py @@ -265,7 +265,7 @@ def get_domain(model): domain.pop(k) return domain - def get_grid(model, **dis_opts): + def get_endo_grid(model, **dis_opts): import copy domain = model.get_domain() a = np.array([e[0] for e in domain.values()]) diff --git a/dolo/compiler/objects.py b/dolo/compiler/objects.py index 8e0261f5..fef247bf 100644 --- a/dolo/compiler/objects.py +++ b/dolo/compiler/objects.py @@ -1,17 +1,10 @@ -# import numpy as np -# -# from dolo.numeric.processes import MvNormal, DiscreteMarkovProcess, VAR1, MarkovProduct from dolo.numeric.processes_iid import * -# from dolo.numeric.grids import CartesianGrid, SmolyakGrid -# Normal = MvNormal -# MarkovChain = DiscreteMarkovProcess -# AR1 = VAR1 + # from dataclasses import dataclass from dolo.compiler.language import language_element # not sure we'll keep that import numpy as np -import typing from typing import List, Union Scalar = Union[int, float] @@ -52,18 +45,6 @@ class MvNormal: signature = {'Mu': 'list(float)', 'Sigma': 'Matrix'} -@language_element -def Matrix(*lines): - vec = np.array(lines) - assert(vec.ndim==2) - return vec - - -@language_element -def Vector(*els): - vec = np.array(els) - assert(vec.ndim==1) - return vec #%% diff --git a/dolo/numeric/decision_rule.py b/dolo/numeric/decision_rule.py index f5c93754..4ef91779 100644 --- a/dolo/numeric/decision_rule.py +++ b/dolo/numeric/decision_rule.py @@ -155,7 +155,7 @@ def eval_ijs(self, i, j, s): @multimethod def get_coefficients(itp, exo_grid: CartesianGrid, endo_grid: CartesianGrid, interp_type: Linear, x): - grid = cat_grids(exo_grid, endo_grid) # one single CartesianGrid + grid = exo_grid + endo_grid xx = x.reshape(tuple(grid.n)+(-1,)) return xx @@ -165,10 +165,11 @@ def eval_ms(itp, exo_grid: CartesianGrid, endo_grid: CartesianGrid, interp_type: assert(m.ndim==s.ndim==2) - grid = cat_grids(exo_grid, endo_grid) # one single CartesianGrid + grid = exo_grid + endo_grid # one single CartesianGrid + coeffs = itp.coefficients - d = len(grid.n) - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) + + gg = grid.__numba_repr__() from interpolation.splines import eval_linear x = np.concatenate([m, s], axis=1) @@ -188,11 +189,9 @@ def eval_is(itp, exo_grid: CartesianGrid, endo_grid: CartesianGrid, interp_type: def get_coefficients(itp, exo_grid: CartesianGrid, endo_grid: CartesianGrid, interp_type: Cubic, x): from interpolation.splines.prefilter_cubic import prefilter_cubic - grid = cat_grids(exo_grid, endo_grid) # one single CartesianGrid + grid = exo_grid + endo_grid # one single CartesianGrid x = x.reshape(tuple(grid.n)+(-1,)) - d = len(grid.n) - # this gg could be stored as a member of itp - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) + gg = grid.__numba_repr__() return prefilter_cubic(gg, x) @multimethod @@ -202,10 +201,10 @@ def eval_ms(itp, exo_grid: CartesianGrid, endo_grid: CartesianGrid, interp_type: assert(m.ndim==s.ndim==2) - grid = cat_grids(exo_grid, endo_grid) # one single CartesianGrid + grid = exo_grid + endo_grid # one single CartesianGrid coeffs = itp.coefficients - d = len(grid.n) - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) + + gg = grid.__numba_repr__() x = np.concatenate([m, s], axis=1) @@ -230,11 +229,8 @@ def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_ty from interpolation.splines import eval_linear assert(s.ndim==2) - - grid = endo_grid # one single CartesianGrid coeffs = itp.coefficients[i] - d = len(grid.n) - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) + gg = endo_grid.__numba_repr__() return eval_linear(gg, coeffs, s) @@ -244,10 +240,7 @@ def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_ty @multimethod def get_coefficients(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_type: Cubic, x): from interpolation.splines.prefilter_cubic import prefilter_cubic - grid = endo_grid # one single CartesianGrid - d = len(grid.n) - # this gg could be stored as a member of itp - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) + gg = endo_grid.__numba_repr__() return [prefilter_cubic(gg, x[i].reshape( tuple(grid.n) + (-1,))) for i in range(x.shape[0])] @@ -256,12 +249,8 @@ def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_ty from interpolation.splines import eval_cubic assert(s.ndim==2) - - grid = endo_grid # one single CartesianGrid coeffs = itp.coefficients[i] - d = len(grid.n) - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) - + gg = endo_grid.__numba_repr__() return eval_cubic(gg, coeffs, s) @@ -277,13 +266,26 @@ def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_ty from interpolation.splines import eval_linear assert(s.ndim==2) - grid = endo_grid # one single CartesianGrid coeffs = itp.coefficients[i] - d = len(grid.n) - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) + gg = endo_grid.__numba_repr__() return eval_linear(gg, coeffs, s) +# Empty x Cartesian x Linear + +@multimethod +def get_coefficients(itp, exo_grid: EmptyGrid, endo_grid: CartesianGrid, interp_type: Linear, x): + grid = exo_grid + endo_grid + xx = x.reshape(tuple(grid.n)+(-1,)) + return xx + +@multimethod +def eval_s(itp, exo_grid: EmptyGrid, endo_grid: CartesianGrid, interp_type: Linear, s): + from interpolation.splines import eval_linear + assert(s.ndim==2) + coeffs = itp.coefficients + gg = endo_grid.__numba_repr__() + return eval_linear(gg, coeffs, s) # Empty x Cartesian x Cubic @@ -291,21 +293,17 @@ def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_ty def get_coefficients(itp, exo_grid: EmptyGrid, endo_grid: CartesianGrid, interp_type: Cubic, x): from interpolation.splines.prefilter_cubic import prefilter_cubic grid = endo_grid # one single CartesianGrid - d = len(grid.n) - # this gg could be stored as a member of itp - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) + gg = endo_grid.__numba_repr__() return prefilter_cubic(gg, x[0].reshape( tuple(grid.n) + (-1,))) + @multimethod def eval_s(itp, exo_grid: EmptyGrid, endo_grid: CartesianGrid, interp_type: Cubic, s): from interpolation.splines import eval_cubic assert(s.ndim==2) - grid = endo_grid # one single CartesianGrid coeffs = itp.coefficients - d = len(grid.n) - gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) - + gg = endo_grid.__numba_repr__() return eval_cubic(gg, coeffs, s) @multimethod @@ -401,7 +399,7 @@ def __init__(self, values, model=None): self.p = model.calibration['parameters'] self.exo_grid = model.exogenous.discretize() # this is never used - self.endo_grid = model.get_grid() + self.endo_grid = model.get_endo_grid() self.gufun = gufun def eval_ms(self, m, s): diff --git a/dolo/numeric/grids.py b/dolo/numeric/grids.py index ca201366..7944cb99 100644 --- a/dolo/numeric/grids.py +++ b/dolo/numeric/grids.py @@ -41,6 +41,9 @@ def n_nodes(self): def node(self, i): return None + def __add__(self, g): + return g + class PointGrid(Grid): type = 'point' @@ -72,22 +75,47 @@ def __init__(self, nodes): class CartesianGrid(Grid): - type = 'cartesian' + pass + +class UniformCartesianGrid(CartesianGrid): + + type = 'UniformCartesian' def __init__(self, min, max, n=[]): self.d = len(min) + + # this should be a tuple self.min = np.array(min, dtype=float) self.max = np.array(max, dtype=float) if len(n) == 0: self.n = np.zeros(n, dtype=int) + 20 else: self.n = np.array(n, dtype=int) + + # this should be done only on request. self.__nodes__ = mlinspace(self.min, self.max, self.n) + # def node(i:) + # pass + + def __add__(self, g): + + if not isinstance(g, UniformCartesianGrid): + raise Exception("Not implemented.") + n = np.array( tuple(self.n) + tuple(g.n)) + min = np.array( tuple(self.min) + tuple(self.min) ) + max = np.array( tuple(self.max) + tuple(self.max) ) -class NonUniformCartesianGrid(Grid): + return UniformCartesianGrid(min, max, n) + + def __numba_repr__(self): + return tuple([(self.min[i], self.max[i], self.n[i]) for i in range(self.d)]) + + + +class NonUniformCartesianGrid(CartesianGrid): type = "NonUniformCartesian" @@ -95,7 +123,19 @@ def __init__(self, list_of_nodes): list_of_nodes = [np.array(l) for l in list_of_nodes] self.min = [min(l) for l in list_of_nodes] self.max = [max(l) for l in list_of_nodes] + self.n = np.array([(len(e)) for e in list_of_nodes]) + # this should be done only on request. self.__nodes__ = cartesian(list_of_nodes) + self.list_of_nodes = list_of_nodes # think of a better name + + def __add__(self, g): + + if not isinstance(g, NonUniformCartesianGrid): + raise Exception("Not implemented.") + return NonUniformCartesianGrid( self.list_of_nodes + g.list_of_nodes ) + + def __numba_repr__(self): + return tuple([np.array(e) for e in self.list_of_nodes]) diff --git a/examples/models/rbc_iid.yaml b/examples/models/rbc_iid.yaml index 680f7270..eb071534 100644 --- a/examples/models/rbc_iid.yaml +++ b/examples/models/rbc_iid.yaml @@ -76,7 +76,7 @@ exogenous: !Normal domain: z: [-2*sig_z/(1-rho^2)^0.5, 2*sig_z/(1-rho^2)^0.5] k: [k*0.5, k*1.5] - + options: grid: !Cartesian orders: [5, 50] diff --git a/trash/dolo/algos/dtcscc/accuracy.py b/trash/dolo/algos/dtcscc/accuracy.py index 78db39a1..0426fdd7 100644 --- a/trash/dolo/algos/dtcscc/accuracy.py +++ b/trash/dolo/algos/dtcscc/accuracy.py @@ -69,7 +69,7 @@ def omega(model, dr, n_exp=10000, grid={}, bounds=None, epsilons = np.random.multivariate_normal(mean, sigma, n_draws) weights = np.ones(epsilons.shape[0])/n_draws - approx = model.get_grid(**grid) + approx = model.get_endo_grid(**grid) a = approx.a b = approx.b orders = approx.orders diff --git a/trash/dolo/algos/dtcscc/nonlinearsystem.py b/trash/dolo/algos/dtcscc/nonlinearsystem.py index e307f18d..6bfe48ef 100644 --- a/trash/dolo/algos/dtcscc/nonlinearsystem.py +++ b/trash/dolo/algos/dtcscc/nonlinearsystem.py @@ -51,7 +51,7 @@ def nonlinear_system(model, dr0=None, maxit=10, tol=1e-8, grid={}, distribution distrib = model.get_distribution(**distribution) nodes, weights = distrib.discretize() - approx = model.get_grid(**grid) + approx = model.get_endo_grid(**grid) ms = create_interpolator(approx, approx.interpolation) grid = ms.grid diff --git a/trash/dolo/algos/dtcscc/parameterized_expectations.py b/trash/dolo/algos/dtcscc/parameterized_expectations.py index 245a4686..2b9ad35e 100644 --- a/trash/dolo/algos/dtcscc/parameterized_expectations.py +++ b/trash/dolo/algos/dtcscc/parameterized_expectations.py @@ -65,7 +65,7 @@ def vprint(t): if direct is True: d = model.functions['direct_response'] - approx = model.get_grid(**grid) + approx = model.get_endo_grid(**grid) grid = approx.grid interp_type = approx.interpolation dr = create_interpolator(approx, interp_type) # Interp for control @@ -263,7 +263,7 @@ def parameterized_expectations_direct(model, verbose=False, dr0=None, if pert_order > 1: raise Exception("Perturbation order > 1 not supported (yet).") - approx = model.get_grid(**grid) + approx = model.get_endo_grid(**grid) grid = approx.grid interp_type = approx.interpolation dr = create_interpolator(approx, interp_type) diff --git a/trash/dolo/algos/dtcscc/time_iteration_2.py b/trash/dolo/algos/dtcscc/time_iteration_2.py index a2f0c380..6c8604dc 100644 --- a/trash/dolo/algos/dtcscc/time_iteration_2.py +++ b/trash/dolo/algos/dtcscc/time_iteration_2.py @@ -69,7 +69,7 @@ def parameterized_expectations(model, verbose=False, dr0=None, if pert_order > 1: raise Exception("Perturbation order > 1 not supported (yet).") - approx = model.get_grid(**grid) + approx = model.get_endo_grid(**grid) grid = approx.grid interp_type = approx.interpolation dr = create_interpolator(approx, interp_type) @@ -230,7 +230,7 @@ def parameterized_expectations_direct(model, verbose=False, dr0=None, if pert_order > 1: raise Exception("Perturbation order > 1 not supported (yet).") - approx = model.get_grid(**grid) + approx = model.get_endo_grid(**grid) grid = approx.grid interp_type = approx.interpolation dr = create_interpolator(approx, interp_type)