Computes a perfect foresight simulation using a stacked-time algorithm.

    The initial state is specified either by providing a series of exogenous
    shocks and assuming the model is initially in equilibrium with the first
    value of the shock, or by specifying an initial value for the states.

    Parameters
    ----------
    model : NumericModel
        "dtcscc" model to be solved
    shocks : array-like, dict, or pandas.DataFrame
        A specification of the shocks to the model. Can be any of the
        following (note by "declaration order" below we mean the order
        of `model.symbols["shocks"]`):

        - A 1d numpy array-like specifying a time series for a single
          shock, or all shocks stacked into a single array.
        - A 2d numpy array where each column specifies the time series
          for one of the shocks in declaration order. This must be an
          `N` by number of shocks 2d array.
        - A dict where keys are strings found in
          `model.symbols["shocks"]` and values are a time series of
          values for that shock. For model shocks that do not appear in
          this dict, the shock is set to the calibrated value. Note
          that this interface is the most flexible as it allows the user
          to pass values for only a subset of the model shocks and it
          allows the passed time series to be of different lengths.
        - A DataFrame where columns map shock names into time series.
          The same assumptions and behavior that are used in the dict
          case apply here

        If nothing is given here, `shocks` is set equal to the
        calibrated values found in `model.calibration["shocks"]` for
        all periods.

        If the length of any time-series in shocks is less than `T`
        (see below) it is assumed that that particular shock will
        remain at the final given value for the duration of the
        simulaiton.
    start_states : ndarray or dict
        a vector with the value of initial states, or a calibration
        dictionary with the initial values of states and controls
    T : int
        horizon for the perfect foresight simulation
    maxit : int
        maximum number of iteration for the nonlinear solver
    verbose : boolean
        if True, the solver displays iterations
    tol : float
        stopping criterium for the nonlinear solver
    ignore_constraints : bool
        if True, complementarity constraints are ignored.

    Returns
    -------
    pandas dataframe
        a dataframe with T+1 observations of the model variables along the
        simulation (states, controls, auxiliaries). The first observation is
        the steady-state corresponding to the first value of the shocks. The
        simulation should return to a steady-state corresponding to the last
        value of the exogenous shocks.

def deterministic_solve(model, shocks=None, start_states=None, T=100,
                        ignore_constraints=False, maxit=100,
                        initial_guess=None, verbose=False, tol=1e-6):
    """
    Computes a perfect foresight simulation using a stacked-time algorithm.

    The initial state is specified either by providing a series of exogenous
    shocks and assuming the model is initially in equilibrium with the first
    value of the shock, or by specifying an initial value for the states.

    Parameters
    ----------
    model : NumericModel
        "dtcscc" model to be solved
    shocks : array-like, dict, or pandas.DataFrame
        A specification of the shocks to the model. Can be any of the
        following (note by "declaration order" below we mean the order
        of `model.symbols["shocks"]`):

        - A 1d numpy array-like specifying a time series for a single
          shock, or all shocks stacked into a single array.
        - A 2d numpy array where each column specifies the time series
          for one of the shocks in declaration order. This must be an
          `N` by number of shocks 2d array.
        - A dict where keys are strings found in
          `model.symbols["shocks"]` and values are a time series of
          values for that shock. For model shocks that do not appear in
          this dict, the shock is set to the calibrated value. Note
          that this interface is the most flexible as it allows the user
          to pass values for only a subset of the model shocks and it
          allows the passed time series to be of different lengths.
        - A DataFrame where columns map shock names into time series.
          The same assumptions and behavior that are used in the dict
          case apply here

        If nothing is given here, `shocks` is set equal to the
        calibrated values found in `model.calibration["shocks"]` for
        all periods.

        If the length of any time-series in shocks is less than `T`
        (see below) it is assumed that that particular shock will
        remain at the final given value for the duration of the
        simulaiton.
    start_states : ndarray or dict
        a vector with the value of initial states, or a calibration
        dictionary with the initial values of states and controls
    T : int
        horizon for the perfect foresight simulation
    maxit : int
        maximum number of iteration for the nonlinear solver
    verbose : boolean
        if True, the solver displays iterations
    tol : float
        stopping criterium for the nonlinear solver
    ignore_constraints : bool
        if True, complementarity constraints are ignored.

    Returns
    -------
    pandas dataframe
        a dataframe with T+1 observations of the model variables along the
        simulation (states, controls, auxiliaries). The first observation is
        the steady-state corresponding to the first value of the shocks. The
        simulation should return to a steady-state corresponding to the last
        value of the exogenous shocks.

    """

    # TODO:

    # if model.model_spec == 'fga':
    #     from dolo.compiler.converter import GModel_fg_from_fga
    #     model = GModel_fg_from_fga(model)

    # definitions
    n_s = len(model.calibration['states'])
    n_x = len(model.calibration['controls'])

    epsilons = _shocks_to_epsilons(model, shocks, T)

    # final initial and final steady-states consistent with exogenous shocks
    if start_states is None:
        start_states = model.calibration

    if isinstance(start_states, dict):
        # at least that part is clear
        start_equilibrium = start_states
        start_s = start_equilibrium['states']
        start_x = start_equilibrium['controls']
        final_s = start_equilibrium['states']
        final_x = start_equilibrium['controls']
    elif isinstance(start_states, np.ndarray):
        start_s = start_states
        start_x = model.calibration['controls']
        final_s = model.calibration['states']
        final_x = model.calibration['controls']

    # if start_constraints:
    #     # we ignore start_constraints
    #     start_dict.update(start_constraints)
    #     final_equilibrium = start_constraints.copy()
    # else:
    #     final_eqm = find_deterministic_equilibrium(model,
    #                                                constraints=final_dict)
    # final_s = final_eqm['states']
    # final_x = final_eqm['controls']
    #
    # start_s = start_states
    # start_x = final_x

    # TODO: for start_x, it should be possible to use first order guess

    final = np.concatenate([final_s, final_x])
    start = np.concatenate([start_s, start_x])

    if verbose is True:
        print("Initial states : {}".format(start_s))
        print("Final controls : {}".format(final_x))

    p = model.calibration['parameters']

    if initial_guess is None:
        initial_guess = np.row_stack([start*(1-l) + final*l
                                         for l in linspace(0.0, 1.0, T+1)])

    else:
        if isinstance(initial_guess, pd.DataFrame):
            initial_guess = np.array(initial_guess).T.copy()
        initial_guess = initial_guess[:, :n_s+n_x]
        initial_guess[0, :n_s] = start_s
        initial_guess[-1, n_s:] = final_x

    sh = initial_guess.shape

    if model.x_bounds and not ignore_constraints:
        initial_states = initial_guess[:, :n_s]
        [lb, ub] = [u(initial_states, p) for u in model.x_bounds]
        lower_bound = initial_guess*0 - np.inf
        lower_bound[:, n_s:] = lb
        upper_bound = initial_guess*0 + np.inf
        upper_bound[:, n_s:] = ub
        test1 = max(lb.max(axis=0) - lb.min(axis=0))
        test2 = max(ub.max(axis=0) - ub.min(axis=0))
        if test1 > 0.00000001 or test2 > 0.00000001:
            msg = "Not implemented: perfect foresight solution requires that "
            msg += "controls have constant bounds."
            raise Exception(msg)
    else:
        ignore_constraints = True
        lower_bound = None
        upper_bound = None

    nn = sh[0]*sh[1]

    def fobj(vec):
        o = det_residual(model, vec.reshape(sh), start_s, final_x, epsilons)[0]
        return o.ravel()

    if not ignore_constraints:
        def ff(vec):
            return det_residual(model, vec.reshape(sh), start_s, final_x,
                                epsilons, jactype='sparse')

        x0 = initial_guess.ravel()
        sol, nit = ncpsolve(ff, lower_bound.ravel(), upper_bound.ravel(),
                            initial_guess.ravel(), verbose=verbose,
                            maxit=maxit, tol=tol, jactype='sparse')

        sol = sol.reshape(sh)

    else:

        def ff(vec):
            return det_residual(model, vec.reshape(sh), start_s, final_x,
                                epsilons, diff=False).ravel()

        x0 = initial_guess.ravel()
        sol = root(ff, x0, jac=False)
        res = ff(sol.x)
        sol = sol.x.reshape(sh)

    if 'auxiliary' in model.functions:
        colnames = (model.symbols['states'] + model.symbols['controls'] +
                    model.symbols['auxiliaries'])
        # compute auxiliaries
        y = model.functions['auxiliary'](sol[:, :n_s], sol[:, n_s:], p)
        sol = np.column_stack([sol, y])
    else:
        colnames = model.symbols['states'] + model.symbols['controls']

    sol = np.column_stack([sol, epsilons])
    colnames = colnames + model.symbols['shocks']

    ts = pd.DataFrame(sol, columns=colnames)
    return ts