Compute value function corresponding to policy ``dr``

    Parameters:
    -----------

    model:
        "dtcscc" model. Must contain a 'value' function.

    mdr:
        decision rule to evaluate

    Returns:
    --------

    decision rule:
        value function (a function of the space similar to a decision rule
        object)

def evaluate_policy(model,
                    mdr,
                    tol=1e-8,
                    maxit=2000,
                    grid={},
                    verbose=True,
                    initial_guess=None,
                    hook=None,
                    integration_orders=None,
                    details=False,
                    interp_type='cubic'):
    """Compute value function corresponding to policy ``dr``

    Parameters:
    -----------

    model:
        "dtcscc" model. Must contain a 'value' function.

    mdr:
        decision rule to evaluate

    Returns:
    --------

    decision rule:
        value function (a function of the space similar to a decision rule
        object)

    """

    process = model.exogenous
    dprocess = process.discretize()

    n_ms = dprocess.n_nodes()  # number of exogenous states
    n_mv = dprocess.n_inodes(
        0)  # this assume number of integration nodes is constant

    x0 = model.calibration['controls']
    v0 = model.calibration['values']
    parms = model.calibration['parameters']
    n_x = len(x0)
    n_v = len(v0)
    n_s = len(model.symbols['states'])

    endo_grid = model.get_grid(**grid)
    exo_grid = dprocess.grid

    if initial_guess is not None:
        mdrv = initial_guess
    else:
        mdrv = DecisionRule(exo_grid, endo_grid, interp_type=interp_type)

    grid = mdrv.endo_grid.nodes()
    N = grid.shape[0]

    if isinstance(mdr, np.ndarray):
        controls = mdr
    else:
        controls = np.zeros((n_ms, N, n_x))
        for i_m in range(n_ms):
            controls[i_m, :, :] = mdr.eval_is(i_m, grid)

    values_0 = np.zeros((n_ms, N, n_v))
    if initial_guess is None:
        for i_m in range(n_ms):
            values_0[i_m, :, :] = v0[None, :]
    else:
        for i_m in range(n_ms):
            values_0[i_m, :, :] = initial_guess.eval_is(i_m, grid)

    val = model.functions['value']
    g = model.functions['transition']

    sh_v = values_0.shape

    err = 10
    inner_maxit = 50
    it = 0

    if verbose:
        headline = '|{0:^4} | {1:10} | {2:8} | {3:8} |'.format(
            'N', ' Error', 'Gain', 'Time')
        stars = '-' * len(headline)
        print(stars)
        print(headline)
        print(stars)

    t1 = time.time()

    err_0 = np.nan

    verbit = (verbose == 'full')

    while err > tol and it < maxit:

        it += 1

        t_start = time.time()

        mdrv.set_values(values_0.reshape(sh_v))
        values = update_value(val, g, grid, controls, values_0, mdr, mdrv,
                              dprocess, parms).reshape((-1, n_v))
        err = abs(values.reshape(sh_v) - values_0).max()

        err_SA = err / err_0
        err_0 = err

        values_0 = values.reshape(sh_v)

        t_finish = time.time()
        elapsed = t_finish - t_start

        if verbose:
            print('|{0:4} | {1:10.3e} | {2:8.3f} | {3:8.3f} |'.format(
                it, err, err_SA, elapsed))

    # values_0 = values.reshape(sh_v)

    t2 = time.time()

    if verbose:
        print(stars)
        print("Elapsed: {} seconds.".format(t2 - t1))
        print(stars)

    if not details:
        return mdrv
    else:
        return EvaluationResult(mdrv, it, tol, err)