Work on d.r.

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,20 +260,20 @@ def vprint(*args, **kwargs):
     parms = model.calibration['parameters']
 
     dp = dprocess
-    
+
     if dp is None:
         dp = model.exogenous.discretize()
 
     n_m = max(dp.n_nodes,1)
 
     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'

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']
 

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:

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"])

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()])

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
 
 #%%
 

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,35 +266,44 @@ 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
 
 @multimethod
 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):