'Solve the optimization problem by modified policy iteration. See
the `solve` method.'

def modified_policy_iteration(self, v_init=None, epsilon=None, max_iter=None, k=20):

    if (self.beta == 1):
        raise NotImplementedError(self._error_msg_no_discounting)
    if (max_iter is None):
        max_iter = self.max_iter
    if (epsilon is None):
        epsilon = self.epsilon
    def span(z):
        return (z.max() - z.min())
    def midrange(z):
        return ((z.min() + z.max()) / 2)
    v = np.empty(self.num_states)
    if (v_init is None):
        v[:] = (self.R[(self.R > (- np.inf))].min() / (1 - self.beta))
    else:
        v[:] = v_init
    u = np.empty(self.num_states)
    sigma = np.empty(self.num_states, dtype=int)
    try:
        tol = ((epsilon * (1 - self.beta)) / self.beta)
    except ZeroDivisionError:
        tol = np.inf
    for i in range(max_iter):
        self.bellman_operator(v, Tv=u, sigma=sigma)
        diff = (u - v)
        if (span(diff) < tol):
            v[:] = (u + ((midrange(diff) * self.beta) / (1 - self.beta)))
            break
        self.operator_iteration(T=self.T_sigma(sigma), v=u, max_iter=k)
        v[:] = u
    num_iter = (i + 1)
    res = DPSolveResult(v=v, sigma=sigma, num_iter=num_iter, mc=self.controlled_mc(sigma), method='modified    policy    iteration', epsilon=epsilon, max_iter=max_iter, k=k)
    return res
