code: 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:  # Raised if beta = 0
              tol = np.inf
  
          for i in range(max_iter):
              # Policy improvement
              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
              # Partial policy evaluation with k iterations
              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
docs: Solve the optimization problem by modified policy iteration. See
          the `solve` method.