Solve for the value function and associated Markov decision rule by iterating over
    the value function.

    Parameters:
    -----------
    model :
        "dtmscc" model. Must contain a 'felicity' function.
    grid :
        grid options
    dr :
        decision rule to evaluate

    Returns:
    --------
    mdr : Markov decision rule
        The solved decision rule/policy function
    mdrv: decision rule
        The solved value function

def value_iteration(model,
                    grid={},
                    tol=1e-6,
                    maxit=500,
                    maxit_howard=20,
                    verbose=False,
                    details=True):
    """
    Solve for the value function and associated Markov decision rule by iterating over
    the value function.

    Parameters:
    -----------
    model :
        "dtmscc" model. Must contain a 'felicity' function.
    grid :
        grid options
    dr :
        decision rule to evaluate

    Returns:
    --------
    mdr : Markov decision rule
        The solved decision rule/policy function
    mdrv: decision rule
        The solved value function
    """

    transition = model.functions['transition']
    felicity = model.functions['felicity']
    controls_lb = model.functions['controls_lb']
    controls_ub = model.functions['controls_ub']

    parms = model.calibration['parameters']
    discount = model.calibration['beta']

    x0 = model.calibration['controls']
    m0 = model.calibration['exogenous']
    s0 = model.calibration['states']
    r0 = felicity(m0, s0, x0, parms)

    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

    endo_grid = model.get_grid(**grid)

    exo_grid = dprocess.grid

    mdrv = DecisionRule(exo_grid, endo_grid)

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

    mdr = constant_policy(model)
    
    controls_0 = np.zeros((n_ms, N, n_x))
    for i_ms in range(n_ms):
        controls_0[i_ms, :, :] = mdr.eval_is(i_ms, grid)

    values_0 = np.zeros((n_ms, N, 1))
    # for i_ms in range(n_ms):
    #     values_0[i_ms, :, :] = mdrv(i_ms, grid)

    mdr = DecisionRule(exo_grid, endo_grid)
    # mdr.set_values(controls_0)

    # THIRD: value function iterations until convergence
    it = 0
    err_v = 100
    err_v_0 = 0
    gain_v = 0.0
    err_x = 100
    err_x_0 = 0
    tol_x = 1e-5
    tol_v = 1e-7

    itprint = IterationsPrinter(
        ('N', int), ('Error_V', float), ('Gain_V', float), ('Error_x', float),
        ('Gain_x', float), ('Eval_n', int), ('Time', float),
        verbose=verbose)
    itprint.print_header('Start value function iterations.')

    while (it < maxit) and (err_v > tol or err_x > tol_x):

        t_start = time.time()
        it += 1

        mdr.set_values(controls_0)
        if it > 2:
            ev = evaluate_policy(
                model, mdr, initial_guess=mdrv, verbose=False, details=True)
        else:
            ev = evaluate_policy(model, mdr, verbose=False, details=True)

        mdrv = ev.solution
        for i_ms in range(n_ms):
            values_0[i_ms, :, :] = mdrv.eval_is(i_ms, grid)

        values = values_0.copy()
        controls = controls_0.copy()

        for i_m in range(n_ms):
            m = dprocess.node(i_m)
            for n in range(N):
                s = grid[n, :]
                x = controls[i_m, n, :]
                lb = controls_lb(m, s, parms)
                ub = controls_ub(m, s, parms)
                bnds = [e for e in zip(lb, ub)]

                def valfun(xx):
                    return -choice_value(transition, felicity, i_m, s, xx,
                                         mdrv, dprocess, parms, discount)[0]

                res = scipy.optimize.minimize(valfun, x, bounds=bnds)
                controls[i_m, n, :] = res.x
                values[i_m, n, 0] = -valfun(x)

        # compute error, update value and dr
        err_x = abs(controls - controls_0).max()
        err_v = abs(values - values_0).max()
        t_end = time.time()
        elapsed = t_end - t_start

        values_0 = values
        controls_0 = controls

        gain_x = err_x / err_x_0
        gain_v = err_v / err_v_0

        err_x_0 = err_x
        err_v_0 = err_v

        itprint.print_iteration(
            N=it,
            Error_V=err_v,
            Gain_V=gain_v,
            Error_x=err_x,
            Gain_x=gain_x,
            Eval_n=ev.iterations,
            Time=elapsed)

    itprint.print_finished()

    mdr = DecisionRule(exo_grid, endo_grid)
    mdr.set_values(controls)
    mdrv.set_values(values_0)

    if not details:
        return mdr, mdrv
    else:
        return ValueIterationResult(
            mdr,             #:AbstractDecisionRule
            mdrv,            #:AbstractDecisionRule
            it,              #:Int
            dprocess,        #:AbstractDiscretizedProcess
            err_x<tol_x,     #:Bool
            tol_x,           #:Float64
            err_x,           #:Float64
            err_v<tol_v,     #:Bool
            tol_v,           #:Float64
            err_v,           #:Float64
            None,            #log:     #:ValueIterationLog
            None             #trace:   #:Union{Nothing,IterationTrace
        )