#cython: profile=False
#cython: boundscheck=False
#cython: wraparound=False
#cython: cdivision=False
"""
State Space Models

Author: Chad Fulton  
License: Simplified-BSD
"""

{{py:

TYPES = {
    "s": ("np.float32_t", "np.float32", "np.NPY_FLOAT32"),
    "d": ("np.float64_t", "float", "np.NPY_FLOAT64"),
    "c": ("np.complex64_t", "np.complex64", "np.NPY_COMPLEX64"),
    "z": ("np.complex128_t", "complex", "np.NPY_COMPLEX128"),
}

}}

# Typical imports
cimport numpy as np
import numpy as np
from statsmodels.src.math cimport *
cimport scipy.linalg.cython_blas as blas
cimport scipy.linalg.cython_lapack as lapack

from statsmodels.tsa.statespace._kalman_filter cimport (
    MEMORY_NO_LIKELIHOOD, MEMORY_NO_FORECAST_COV, MEMORY_NO_STD_FORECAST,
    MEMORY_NO_SMOOTHING, FILTER_CONCENTRATED, FILTER_CHANDRASEKHAR)

{{for prefix, types in TYPES.items()}}
{{py:cython_type, dtype, typenum = types}}
{{py:
combined_prefix = prefix
combined_cython_type = cython_type
if prefix == 'c':
    combined_prefix = 'z'
    combined_cython_type = 'np.complex128_t'
if prefix == 's':
    combined_prefix = 'd'
    combined_cython_type = 'np.float64_t'
combined_suffix = ''
if combined_prefix == 'z':
    combined_suffix = 'u'
}}

# ### Univariate Kalman filter
#
# The following are the routines as defined in the univariate Kalman filter.
#
# See Durbin and Koopman (2012) Chapter 6.4

cdef int {{prefix}}forecast_univariate({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):

    # Constants
    cdef:
        int i, j, k
        int inc = 1
        int forecast_cov_t = kfilter.t
        int check
        {{cython_type}} forecast_error_cov
        {{cython_type}} forecast_error_cov_inv
        {{cython_type}} forecast_error_cov_inv_prev = 0.0

    if kfilter.conserve_memory & MEMORY_NO_FORECAST_COV > 0:
        forecast_cov_t = 1

    # Initialize the filtered states
    {{prefix}}initialize_filtered(kfilter, model)

    # Make sure the loglikelihood is set to zero if necessary

    # Iterate over the observations at time t
    for i in range(model._k_endog):

        # #### Forecast for time t
        # `forecast` $= Z_{t,i} a_{t,i} + d_{t,i}$
        # Note: $Z_{t,i}$ is a row vector starting at [i,0,t] and ending at
        # [i,k_states,t]
        # Note: zdot and cdot are broken, so have to use gemv for those

        # #### Forecast error for time t
        # `forecast_error` $\equiv v_t = y_t -$ `forecast`
        {{prefix}}forecast_error(kfilter, model, i)

        # #### Forecast error covariance matrix for time t
        # $F_{t,i} \equiv Z_{t,i} P_{t,i} Z_{t,i}' + H_{t,i}$
        # TODO what about Kalman convergence?
        # Note: zdot and cdot are broken, so have to use gemv for those
        if not kfilter.converged:
            forecast_error_cov = {{prefix}}forecast_error_cov(kfilter, model, i)
        else:
            forecast_error_cov = kfilter._forecast_error_cov[i + i*kfilter.k_endog]

        # Handle numerical issues that can cause a very small negative
        # forecast_error_cov
        check = {{prefix}}check1(kfilter, forecast_error_cov)
        if check:
            kfilter._forecast_error_cov[i + i*kfilter.k_endog] = 0
            forecast_error_cov = 0

        # If we have a non-zero variance
        # (if we have a zero-variance then we are done with this iteration)
        check = {{prefix}}check2(kfilter, i, forecast_error_cov)
        if check:
            forecast_error_cov_inv = {{prefix}}forecast_error_cov_inv(kfilter, model, i, forecast_error_cov)

            {{prefix}}standardized_forecast_error(kfilter, model, i, forecast_error_cov_inv)

            # Save temporary array data
            {{prefix}}temp_arrays(kfilter, model, i, forecast_error_cov_inv)

            # #### Filtered state for time t
            # $a_{t,i+1} = a_{t,i} + P_{t,i} Z_{t,i}' F_{t,i}^{-1} v_{t,i}$  
            # Make a new temporary array  
            # K_{t,i} = P_{t,i} Z_{t,i}' F_{t,i}^{-1}
            {{prefix}}filtered_state(kfilter, model, i, forecast_error_cov_inv)

            # Chandrasekhar arrays
            if not kfilter.converged and (kfilter.filter_method & FILTER_CHANDRASEKHAR > 0):
                if kfilter.t > 0:
                    forecast_error_cov_inv_prev = 1.0 / kfilter.forecast_error_cov[i, i, forecast_cov_t - 1]
                {{prefix}}chandrasekhar_recursion(kfilter, model, i, forecast_error_cov, forecast_error_cov_inv, forecast_error_cov_inv_prev)

            # #### Filtered state covariance for time t
            # $P_{t,i+1} = P_{t,i} - P_{t,i} Z_{t,i}' F_{t,i}^{-1} Z_{t,i} P_{t,i}'$
            if not kfilter.converged:
                {{prefix}}filtered_state_cov(kfilter, model, i, forecast_error_cov_inv)

            # #### Loglikelihood
            {{prefix}}loglikelihood(kfilter, model, i, forecast_error_cov, forecast_error_cov_inv)
        else:
            # Otherwise, we need to record that this observation is not associated
            # with a loglikelihood step (so that it can be excluded in the denominator
            # when computing the scale)
            kfilter.nobs_kendog_univariate_singular = kfilter.nobs_kendog_univariate_singular + 1

    {{prefix}}symmetry(kfilter, model)

    return 0

cdef void {{prefix}}initialize_filtered({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    cdef int inc = 1
    blas.{{prefix}}copy(&kfilter.k_states, kfilter._input_state, &inc,
                                           kfilter._filtered_state, &inc)
    if not kfilter.converged:
        blas.{{prefix}}copy(&kfilter.k_states2, kfilter._input_state_cov, &inc,
                                                kfilter._filtered_state_cov, &inc)

cdef int {{prefix}}check1({{prefix}}KalmanFilter kfilter, {{cython_type}} forecast_error_cov):
    if not kfilter.converged:
        return forecast_error_cov{{if combined_prefix == 'z'}}.real{{endif}} < 0
    else:
        return False

cdef int {{prefix}}check2({{prefix}}KalmanFilter kfilter, int i, {{cython_type}} forecast_error_cov):
    if not kfilter.converged:
        kfilter.forecast_error_ipiv[i] = forecast_error_cov{{if combined_prefix == 'z'}}.real{{endif}} > kfilter.tolerance_diffuse
    return kfilter.forecast_error_ipiv[i]

cdef {{cython_type}} {{prefix}}forecast_error_cov_inv({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i, {{cython_type}} forecast_error_cov):
    if not kfilter.converged:
        kfilter.forecast_error_fac[i, i] = 1.0 / forecast_error_cov

    return kfilter.forecast_error_fac[i, i]

cdef void {{prefix}}standardized_forecast_error({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i, {{cython_type}} forecast_error_cov_inv):
    if not (kfilter.conserve_memory & MEMORY_NO_STD_FORECAST > 0):
            kfilter._standardized_forecast_error[i] = (
                kfilter._forecast_error[i] * forecast_error_cov_inv**0.5)

cdef void {{prefix}}symmetry({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    cdef int j, k
    # Make final filtered_state_cov symmetric (is not currently symmetric
    # due to use of ?syr or ?her)
    if not kfilter.converged:
        for j in range(model._k_states):      # columns
            for k in range(model._k_states):  # rows
                if k > j: # row > column => in lower triangle
                    kfilter._filtered_state_cov[j + k*kfilter.k_states] = kfilter._filtered_state_cov[k + j*kfilter.k_states]

cdef void {{prefix}}forecast_error({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i):
    cdef:
        int inc = 1
        {{cython_type}} alpha = 1
        {{cython_type}} beta = 0
        int k_states = model._k_states
    # Adjust for a VAR transition (i.e. design = [#, 0], where the zeros
    # correspond to all states except the first k_posdef states)
    if model.subset_design:
        k_states = model._k_posdef

    # `forecast` $= Z_{t,i} a_{t,i} + d_{t,i}$
    {{if combined_prefix == 'd'}}
    kfilter._forecast[i] = (
        model._obs_intercept[i] +
        blas.{{prefix}}dot(&k_states, &model._design[i], &model._k_endog,
                                      kfilter._filtered_state, &inc)
    )
    {{else}}
    blas.{{prefix}}gemv("N", &inc, &k_states,
                   &alpha, kfilter._filtered_state, &inc,
                           &model._design[i], &model._k_endog,
                   &beta, kfilter._tmp0, &inc)
    kfilter._forecast[i] = model._obs_intercept[i] + kfilter._tmp0[0]
    {{endif}}

    # `forecast_error` $\equiv v_t = y_t -$ `forecast`
    kfilter._forecast_error[i] = model._obs[i] - kfilter._forecast[i]

cdef {{cython_type}} {{prefix}}forecast_error_cov({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i):
    cdef:
        int inc = 1
        {{cython_type}} alpha = 1
        {{cython_type}} beta = 0
        {{cython_type}} forecast_error_cov
        int k_states = model._k_states

    # Adjust for a VAR transition (i.e. design = [#, 0], where the zeros
    # correspond to all states except the first k_posdef states)
    if model.subset_design:
        k_states = model._k_posdef

    # *Intermediate calculation* (used just below and then once more)  
    # $M_{t,i} = P_{t,i} Z_{t,i}'$  
    # $(m \times 1) = (m \times m) (1 \times m)'$
    # blas.{{prefix}}gemv("N", &model._k_states, &k_states,
    #       &alpha, kfilter._filtered_state_cov, &kfilter.k_states,
    #               &model._design[i], &model._k_endog,
    #       &beta, &kfilter._M[i*kfilter.k_states], &inc)

    # $F_{t,i} \equiv Z_{t,i} P_{t,i} Z_{t,i}' + H_{t,i}$
    {{if combined_prefix == 'd'}}
    blas.{{prefix}}symv("L", &model._k_states,
          &alpha, kfilter._filtered_state_cov, &kfilter.k_states,
                  &model._design[i], &model._k_endog,
          &beta, &kfilter._M[i*kfilter.k_states], &inc)

    forecast_error_cov = (
        model._obs_cov[i + i*model._k_endog] +
        blas.{{prefix}}dot(&k_states, &model._design[i], &model._k_endog,
                                      &kfilter._M[i*kfilter.k_states], &inc)
    )
    {{else}}
    blas.{{prefix}}symm("R", "L", &inc, &model._k_states,
          &alpha, kfilter._filtered_state_cov, &kfilter.k_states,
                  &model._design[i], &model._k_endog,
          &beta, &kfilter._M[i*kfilter.k_states], &inc)

    blas.{{prefix}}gemv("N", &inc, &k_states,
                   &alpha, &kfilter._M[i*kfilter.k_states], &inc,
                           &model._design[i], &model._k_endog,
                   &beta, kfilter._tmp0, &inc)
    forecast_error_cov = model._obs_cov[i + i*model._k_endog] + kfilter._tmp0[0]
    {{endif}}
    kfilter._forecast_error_cov[i + i*kfilter.k_endog] = forecast_error_cov
    return forecast_error_cov

cdef void {{prefix}}temp_arrays({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i, {{cython_type}} forecast_error_cov_inv):
    cdef:
        int k_states = model._k_states
        tmp_1 = 0

    # Adjust for a VAR transition (i.e. design = [#, 0], where the zeros
    # correspond to all states except the first k_posdef states)
    if model.subset_design:
        k_states = model._k_posdef

    # $\\#_1 = P_{t,i} Z_{t,i}'$ - set above
    # $\\#_2 = v_{t,i} / F_{t,i}$
    kfilter._tmp2[i] = kfilter._forecast_error[i] * forecast_error_cov_inv
    # $\\#_3 = Z_{t,i} / F_{t,i}$
    # $\\#_4 = H_{t,i} / F_{t,i}$
    if not kfilter.converged:
        blas.{{prefix}}copy(&k_states, &model._design[i], &model._k_endog,
                                       &kfilter._tmp3[i], &kfilter.k_endog)
        blas.{{prefix}}scal(&k_states, &forecast_error_cov_inv, &kfilter._tmp3[i], &kfilter.k_endog)

        kfilter._tmp4[i + i*kfilter.k_endog] = model._obs_cov[i + i*model._k_endog] * forecast_error_cov_inv
    elif kfilter.conserve_memory & MEMORY_NO_SMOOTHING > 0:
        # If we're converged and we're not storing these arrays, then we
        # already have the converged values and there's nothing more to do
        pass
    else:
        # If we're converged and we are storing these arrays, then we
        # just need to copy them from the previous iteration
        blas.{{prefix}}copy(&k_states, &kfilter.tmp3[i, 0, kfilter.t - 1], &kfilter.k_endog,
                                       &kfilter._tmp3[i], &kfilter.k_endog)
        kfilter._tmp4[i + i*kfilter.k_endog] = kfilter.tmp4[i, i, kfilter.t - 1]

cdef void {{prefix}}filtered_state({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i, {{cython_type}} forecast_error_cov_inv):
    cdef int j
    # $a_{t,i+1} = a_{t,i} + P_{t,i} Z_{t,i}' F_{t,i}^{-1} v_{t,i}$  
    for j in range(model._k_states):
        if not kfilter.converged:
            kfilter._kalman_gain[j + i*kfilter.k_states] = kfilter._M[j + i*kfilter.k_states] * forecast_error_cov_inv
        kfilter._filtered_state[j] = (
            kfilter._filtered_state[j] +
            kfilter._forecast_error[i] * kfilter._kalman_gain[j + i*kfilter.k_states]
        )

cdef void {{prefix}}filtered_state_cov({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i, {{cython_type}} forecast_error_cov_inv):
    cdef:
        int inc = 1, j, k
        {{cython_type}} scalar = -1.0 * forecast_error_cov_inv
        {{cython_type}} alpha = 1.0
        {{cython_type}} gamma = -1.0
        int k_states = model._k_states
        int k_states1 = model._k_states

    # Adjust for a VAR transition (i.e. design = [#, 0], where the zeros
    # correspond to all states except the first k_posdef states)
    if model.subset_design:
        k_states = model._k_posdef
        if model._k_posdef > model._k_states:
            k_states1 = model._k_posdef + 1

    # $P_{t,i+1} = P_{t,i} - P_{t,i} Z_{t,i}' F_{t,i}^{-1} Z_{t,i} P_{t,i}'$
    # blas.{{prefix}}ger{{combined_suffix}}(&model._k_states, &model._k_states,
    #     &gamma, &kfilter._M[i*kfilter.k_states], &inc,
    #             &kfilter._kalman_gain[i*kfilter.k_states], &inc,
    #     kfilter._filtered_state_cov, &kfilter.k_states
    # )

    {{if combined_prefix == 'd'}}
    blas.{{prefix}}syr("L", &model._k_states,
        &scalar, &kfilter._M[i*kfilter.k_states], &inc,
                 kfilter._filtered_state_cov, &kfilter.k_states
    )
    {{else}}
    blas.{{prefix}}syrk("L", "N", &model._k_states, &inc,
        &scalar, &kfilter._M[i*kfilter.k_states], &kfilter.k_states,
        &alpha, kfilter._filtered_state_cov, &kfilter.k_states)
    {{endif}}

cdef void {{prefix}}chandrasekhar_recursion({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i, {{cython_type}} forecast_error_cov, {{cython_type}} forecast_error_cov_inv, {{cython_type}} forecast_error_cov_inv_prev):
    # Constants
    cdef:
        int inc = 1
        int j, k
        {{cython_type}} alpha = 1.0
        {{cython_type}} beta = 0.0
        {{cython_type}} gamma = -1.0

    # Initialization
    if kfilter.t == 0:
        if i == 0:
            kfilter.CM[:] = 0

        # W[:, i:i+1] = T @ (P @ Z[i].T)
        # W[:, i:i+1] = T @ K @ F[i, i]
        # Note: we scale by forecast error cov here b/c kalman_gain was
        # computed above as K = P @ Z[i].T @ (1 / F[i, i])
        blas.{{prefix}}gemv("N", &model._k_states, &model._k_states,
          &forecast_error_cov, model._transition, &model._k_states,
                  &kfilter._kalman_gain[i * kfilter.k_states], &inc,
          &beta, &kfilter.CW[0, i], &inc)
        # M[i, i] = Finv[i, i]
        kfilter.CM[i, i] = -forecast_error_cov_inv
    # Updating
    else:
        # M.T @ W.T. (p x p) (p x m)
        blas.{{prefix}}gemm("T", "T", &model._k_endog, &model._k_states, &model._k_endog,
          &alpha, &kfilter.CM[0, 0], &kfilter.k_endog,
                  &kfilter.CW[0, 0], &kfilter.k_states,
          &beta, &kfilter.CMW[0, 0], &kfilter.k_endog)
        # MW @ Z[i].T (p x m) (m x 1) -> (p x 1)
        blas.{{prefix}}gemv("N", &model._k_endog, &model._k_states,
          &alpha, &kfilter.CMW[0, 0], &kfilter.k_endog,
                  &model._design[i], &model._k_endog,
          &beta, &kfilter.CMWZ[0, 0], &inc)

        # M = M + MWZ @ MWZ.T / F_prev[i, i]
        # Note: syr / syrk only fills in lower triangle here
        {{if combined_prefix == 'd'}}
        blas.{{prefix}}syr("L", &model._k_endog,
            &forecast_error_cov_inv_prev, &kfilter.CMWZ[0, 0], &inc,
                     &kfilter.CM[0, 0], &kfilter.k_endog)
        {{else}}
        blas.{{prefix}}syrk("L", "N", &model._k_endog, &inc,
            &forecast_error_cov_inv_prev, &kfilter.CMWZ[0, 0], &kfilter.k_endog,
            &alpha, &kfilter.CM[0, 0], &kfilter.k_endog)
        {{endif}}

        # Fill in the upper triangle
        for j in range(model._k_endog):      # columns
            for k in range(model._k_endog):  # rows
                if k > j: # row > column => in lower triangle
                    kfilter.CM[j, k] = kfilter.CM[k, j]

        # Compute W
        # W -> tmpW
        blas.{{prefix}}copy(&model._k_endogstates, &kfilter.CW[0, 0], &inc, &kfilter.CtmpW[0, 0], &inc)

        if i == model.k_endog - 1:
            # W = (T - T @ K @ Z[i]) @ W

            # Compute T @ K: (m x m) (m x 1) -> (m x 1)
            # Note: we previously copied CW -> CtmpW, so overwriting CW is okay
            blas.{{prefix}}gemv("N", &model._k_states, &model._k_states,
              &alpha, model._transition, &model._k_states,
                      &kfilter._kalman_gain[i*kfilter.k_states], &inc,
              &beta, &kfilter.CW[0, 0], &inc)
            # T -> tmp00
            blas.{{prefix}}copy(&model._k_states2, model._transition, &inc, kfilter._tmp00, &inc)
            # T - (T @ K) @ Z[i]: (m x 1) (1 x m) -> (m x m)
            blas.{{prefix}}ger{{combined_suffix}}(&model._k_states, &model._k_states,
                &gamma, &kfilter.CW[0, 0], &inc,
                       &model._design[i], &model._k_endog,
                kfilter._tmp00, &kfilter.k_states)

            # (T - T @ K @ Z[i]) @ tmpW -> W
            blas.{{prefix}}gemm("N", "N", &model._k_states, &model._k_endog, &model._k_states,
                  &alpha, kfilter._tmp00, &kfilter.k_states,
                          &kfilter.CtmpW[0, 0], &kfilter.k_states,
                  &beta, &kfilter.CW[0, 0], &kfilter.k_states)
        else:
            # W = (I - I @ K @ Z[i]) @ W
            # K @ Z[i] (m x 1) (1 x m) -> (m x m)
            kfilter.tmp0[:] = 0
            blas.{{prefix}}ger{{combined_suffix}}(&model._k_states, &model._k_states,
                &alpha, &kfilter._kalman_gain[i*kfilter.k_states], &inc,
                        &model._design[i], &model._k_endog,
                kfilter._tmp0, &kfilter.k_states)

            # W = - K @ Z[i] @ W + W
            blas.{{prefix}}gemm("N", "N", &model._k_states, &model._k_endog, &model._k_states,
                  &gamma, kfilter._tmp0, &kfilter.k_states,
                          &kfilter.CtmpW[0, 0], &kfilter.k_states,
                  &alpha, &kfilter.CW[0, 0], &kfilter.k_states)


cdef void {{prefix}}loglikelihood({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, int i, {{cython_type}} forecast_error_cov, {{cython_type}} forecast_error_cov_inv):
    kfilter._loglikelihood[0] = (
        kfilter._loglikelihood[0] - 0.5*(
            {{combined_prefix}}log(2 * NPY_PI * forecast_error_cov)
        )
    )
    if kfilter.filter_method & FILTER_CONCENTRATED:
        kfilter._scale[0] = kfilter._scale[0] + kfilter._forecast_error[i]**2 * forecast_error_cov_inv
    else:
        kfilter._loglikelihood[0] = kfilter._loglikelihood[0] - 0.5 * (kfilter._forecast_error[i]**2 * forecast_error_cov_inv)

cdef int {{prefix}}updating_univariate({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    # the updating step was performed in the forecast_univariate step
    return 0

cdef int {{prefix}}prediction_univariate({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    # Constants
    cdef:
        int inc = 1
        int i, j
        {{cython_type}} alpha = 1.0
        {{cython_type}} beta = 0.0
        {{cython_type}} gamma = -1.0

    # #### Predicted state for time t+1
    # $a_{t+1} = T_t a_{t,n} + c_t$

    # #### Predicted state covariance matrix for time t+1
    # $P_{t+1} = T_t P_{t,n} T_t' + Q_t^*$
    #
    # TODO check behavior during convergence
    if not model.companion_transition:
        {{prefix}}predicted_state(kfilter, model)
        if not kfilter.converged:
            if kfilter.filter_method & FILTER_CHANDRASEKHAR > 0:
                {{prefix}}predicted_state_cov_chandrasekhar(kfilter, model)
            else:
                {{prefix}}predicted_state_cov(kfilter, model)
    else:
        {{prefix}}companion_predicted_state(kfilter, model)
        if not kfilter.converged:
            {{prefix}}companion_predicted_state_cov(kfilter, model)

    # #### Kalman gain for time t
    # $K_t = T_t P_t Z_t' F_t^{-1}$  
    # Kalman gain calculation done in forecasting step.

    return 0

cdef void {{prefix}}predicted_state({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    cdef:
        int inc = 1
        {{cython_type}} alpha = 1.0

    # $a_{t+1} = T_t a_{t,n} + c_t$
    blas.{{prefix}}copy(&model._k_states, model._state_intercept, &inc, kfilter._predicted_state, &inc)
    blas.{{prefix}}gemv("N", &model._k_states, &model._k_states,
          &alpha, model._transition, &model._k_states,
                  kfilter._filtered_state, &inc,
          &alpha, kfilter._predicted_state, &inc)

cdef void {{prefix}}predicted_state_cov({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    cdef:
        int inc = 1
        {{cython_type}} alpha = 1.0
        {{cython_type}} beta = 0.0

    # $P_{t+1} = T_t P_{t,n} T_t' + Q_t^*$
    blas.{{prefix}}copy(&model._k_states2, model._selected_state_cov, &inc, kfilter._predicted_state_cov, &inc)
    # `tmp0` array used here, dimension $(m \times m)$  

    # $\\#_0 = T_t P_{t|t} $

    # $(m \times m) = (m \times m) (m \times m)$
    # blas.{{prefix}}gemm("N", "N", &model._k_states, &model._k_states, &model._k_states,
    #       &alpha, model._transition, &model._k_states,
    #               kfilter._filtered_state_cov, &kfilter.k_states,
    #       &beta, kfilter._tmp0, &kfilter.k_states)
    blas.{{prefix}}symm("R", "L", &model._k_states, &model._k_states,
          &alpha, kfilter._filtered_state_cov, &kfilter.k_states,
                  model._transition, &model._k_states,
          &beta, kfilter._tmp0, &kfilter.k_states)
    # $P_{t+1} = 1.0 \\#_0 T_t' + 1.0 \\#$  
    # $(m \times m) = (m \times m) (m \times m) + (m \times m)$
    blas.{{prefix}}gemm("N", "T", &model._k_states, &model._k_states, &model._k_states,
          &alpha, kfilter._tmp0, &kfilter.k_states,
                  model._transition, &model._k_states,
          &alpha, kfilter._predicted_state_cov, &kfilter.k_states)

cdef void {{prefix}}predicted_state_cov_chandrasekhar({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    # Constants
    cdef:
        int inc = 1
        {{cython_type}} alpha = 1.0
        {{cython_type}} beta = 0.0
    blas.{{prefix}}copy(&model._k_states2, kfilter._input_state_cov, &inc, kfilter._predicted_state_cov, &inc)
    # M @ W.T. (p x p) (p x m)
    blas.{{prefix}}gemm("N", "T", &model._k_endog, &model._k_states, &model._k_endog,
      &alpha, &kfilter.CM[0, 0], &kfilter.k_endog,
              &kfilter.CW[0, 0], &kfilter.k_states,
      &beta, &kfilter.CMW[0, 0], &kfilter.k_endog)
    # P = P + W M W.T
    blas.{{prefix}}gemm("N", "N", &model._k_states, &model._k_states, &model._k_endog,
      &alpha, &kfilter.CW[0, 0], &kfilter.k_states,
              &kfilter.CMW[0, 0], &kfilter.k_endog,
      &alpha, kfilter._predicted_state_cov, &kfilter.k_states)

cdef void {{prefix}}companion_predicted_state({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    cdef:
        int i
        int inc = 1
        {{cython_type}} alpha = 1.0

    # $a_{t+1} = T_t a_{t,n} + c_t$
    blas.{{prefix}}copy(&model._k_states, model._state_intercept, &inc, kfilter._predicted_state, &inc)
    blas.{{prefix}}gemv("N", &model._k_posdef, &model._k_states,
          &alpha, model._transition, &model._k_states,
                  kfilter._filtered_state, &inc,
          &alpha, kfilter._predicted_state, &inc)

    for i in range(model._k_posdef, model._k_states):
        kfilter._predicted_state[i] = kfilter._predicted_state[i] + kfilter._filtered_state[i - model._k_posdef]

cdef void {{prefix}}companion_predicted_state_cov({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    cdef:
        int i, j, idx
        int inc = 1
        {{cython_type}} alpha = 1.0
        {{cython_type}} beta = 0.0
        {{cython_type}} tmp

    # $P_{t+1} = T_t P_{t,n} T_t' + Q_t^*$
    
    # `tmp0` array used here, dimension $(p \times m)$  
    # $\\#_0 = \phi_t P_{t|t} $

    # $(p \times m) = (p \times m) (m \times m)$
    # TODO: symm?
    blas.{{prefix}}gemm("N", "N", &model._k_posdef, &model._k_states, &model._k_states,
          &alpha, model._transition, &model._k_states,
                  kfilter._filtered_state_cov, &kfilter.k_states,
          &beta, kfilter._tmp0, &kfilter.k_states)
                
    # $P_{t+1} = 1.0 \\#_0 \phi_t' + 1.0 \\#$  
    # $(m \times m) = (p \times m) (m \times p) + (m \times m)$
    blas.{{prefix}}gemm("N", "T", &model._k_posdef, &model._k_posdef, &model._k_states,
          &alpha, kfilter._tmp0, &kfilter.k_states,
                  model._transition, &model._k_states,
          &beta, kfilter._predicted_state_cov, &kfilter.k_states)

    # Fill in the basic matrix blocks
    for i in range(kfilter.k_states):      # columns
        for j in range(kfilter.k_states):  # rows
            idx = j + i*kfilter.k_states

            # Add the Q matrix to the upper-left block
            if i < model._k_posdef and j < model._k_posdef:
                kfilter._predicted_state_cov[idx] = (
                    kfilter._predicted_state_cov[idx] + 
                    model._state_cov[j + i*model._k_posdef]
                )

            # Set the upper-right block to be the first m-p columns of
            # \phi _t P_{t|t}, and the lower-left block to the its transpose
            elif i >= model._k_posdef and j < model._k_posdef:
                tmp = kfilter._tmp0[j + (i-model._k_posdef)*kfilter.k_states]
                kfilter._predicted_state_cov[idx] = tmp
                kfilter._predicted_state_cov[i + j*model._k_states] = tmp

            # Set the lower-right block 
            elif i >= model._k_posdef and j >= model._k_posdef:
                kfilter._predicted_state_cov[idx] = (
                    kfilter._filtered_state_cov[(j - model._k_posdef) + (i - model._k_posdef)*kfilter.k_states]
                )

cdef {{cython_type}} {{prefix}}inverse_noop_univariate({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, {{cython_type}} determinant) except *:
    return -np.inf

cdef {{cython_type}} {{prefix}}loglikelihood_univariate({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model, {{cython_type}} determinant):
    return 0

cdef {{cython_type}} {{prefix}}scale_univariate({{prefix}}KalmanFilter kfilter, {{prefix}}Statespace model):
    return 0

{{endfor}}