Desctrete prediction function

        Input:
            k:int
                Iteration No. Starts at 0. Total number of iterations equal to the
                number of measurements.

            p_m: matrix of size (state_dim, time_series_no)
                Mean value from the previous step. For "multiple time series mode"
                it is matrix, second dimension of which correspond to different
                time series.

            p_P: tuple (Prev_cov, S, V)
                Covariance matrix from the previous step and its SVD decomposition.
                Prev_cov = V * S * V.T The tuple is (Prev_cov, S, V)

            p_dyn_model_callable: object

            calc_grad_log_likelihood: boolean
                Whether to calculate gradient of the marginal likelihood
                of the state-space model. If true then the next parameter must
                provide the extra parameters for gradient calculation.

            p_dm: 3D array (state_dim, time_series_no, parameters_no)
                Mean derivatives from the previous step. For "multiple time series mode"
                it is 3D array, second dimension of which correspond to different
                time series.

            p_dP: 3D array (state_dim, state_dim, parameters_no)
                Mean derivatives from the previous step

        Output:
        ----------------------------
        m_pred, P_pred, dm_pred, dP_pred: metrices, 3D objects
            Results of the prediction steps.

def _kalman_prediction_step_SVD(k, p_m , p_P, p_dyn_model_callable, calc_grad_log_likelihood=False,
                                p_dm = None, p_dP = None):
        """
        Desctrete prediction function

        Input:
            k:int
                Iteration No. Starts at 0. Total number of iterations equal to the
                number of measurements.

            p_m: matrix of size (state_dim, time_series_no)
                Mean value from the previous step. For "multiple time series mode"
                it is matrix, second dimension of which correspond to different
                time series.

            p_P: tuple (Prev_cov, S, V)
                Covariance matrix from the previous step and its SVD decomposition.
                Prev_cov = V * S * V.T The tuple is (Prev_cov, S, V)

            p_dyn_model_callable: object

            calc_grad_log_likelihood: boolean
                Whether to calculate gradient of the marginal likelihood
                of the state-space model. If true then the next parameter must
                provide the extra parameters for gradient calculation.

            p_dm: 3D array (state_dim, time_series_no, parameters_no)
                Mean derivatives from the previous step. For "multiple time series mode"
                it is 3D array, second dimension of which correspond to different
                time series.

            p_dP: 3D array (state_dim, state_dim, parameters_no)
                Mean derivatives from the previous step

        Output:
        ----------------------------
        m_pred, P_pred, dm_pred, dP_pred: metrices, 3D objects
            Results of the prediction steps.

        """

        # covariance from the previous step and its SVD decomposition
        # p_prev_cov = v * S * V.T
        Prev_cov, S_old, V_old = p_P
        #p_prev_cov_tst = np.dot(p_V, (p_S * p_V).T) # reconstructed covariance from the previous step

        # index correspond to values from previous iteration.
        A = p_dyn_model_callable.Ak(k,p_m,Prev_cov) # state transition matrix (or Jacobian)
        Q = p_dyn_model_callable.Qk(k) # state noise matrx. This is necessary for the square root calculation (next step)
        Q_sr = p_dyn_model_callable.Q_srk(k)
        # Prediction step ->
        m_pred = p_dyn_model_callable.f_a(k, p_m, A) # predicted mean

        # coavariance prediction have changed:
        svd_1_matr = np.vstack( ( (np.sqrt(S_old)* np.dot(A,V_old)).T , Q_sr.T) )
        (U,S,Vh) = sp.linalg.svd( svd_1_matr,full_matrices=False, compute_uv=True,
                      overwrite_a=False,check_finite=True)

        # predicted variance computed by the regular method. For testing
        #P_pred_tst = A.dot(Prev_cov).dot(A.T) + Q
        V_new = Vh.T
        S_new = S**2

        P_pred = np.dot(V_new * S_new, V_new.T) # prediction covariance
        P_pred = (P_pred, S_new, Vh.T)
        # Prediction step <-

        # derivatives
        if calc_grad_log_likelihood:
            dA_all_params = p_dyn_model_callable.dAk(k) # derivatives of A wrt parameters
            dQ_all_params = p_dyn_model_callable.dQk(k) # derivatives of Q wrt parameters

            param_number = p_dP.shape[2]

            # p_dm, p_dP - derivatives form the previoius step
            dm_pred = np.empty(p_dm.shape)
            dP_pred = np.empty(p_dP.shape)

            for j in range(param_number):
                dA = dA_all_params[:,:,j]
                dQ = dQ_all_params[:,:,j]

                #dP = p_dP[:,:,j]
                #dm = p_dm[:,:,j]
                dm_pred[:,:,j] = np.dot(dA, p_m) + np.dot(A, p_dm[:,:,j])
                # prediction step derivatives for current parameter:


                dP_pred[:,:,j] = np.dot( dA ,np.dot(Prev_cov, A.T))
                dP_pred[:,:,j] += dP_pred[:,:,j].T
                dP_pred[:,:,j] += np.dot( A ,np.dot(p_dP[:,:,j], A.T)) + dQ

                dP_pred[:,:,j] = 0.5*(dP_pred[:,:,j] + dP_pred[:,:,j].T) #symmetrize
        else:
            dm_pred = None
            dP_pred = None

        return m_pred, P_pred, dm_pred, dP_pred