This function implements the continuous-discrete Kalman Filter algorithm
        These notations for the State-Space model are assumed:
            d/dt x(t) =  F * x(t) + L * w(t);         w(t) ~ N(0, Qc)
            y_{k} = H_{k} * x_{k} + r_{k};            r_{k-1} ~ N(0, R_{k})

        Returns estimated filter distributions x_{k} ~ N(m_{k}, P(k))

        Current Features:
        ----------------------------------------
        1) The function generaly do not modify the passed parameters. If
        it happens then it is an error. There are several exeprions: scalars
        can be modified into a matrix, in some rare cases shapes of
        the derivatives matrices may be changed, it is ignored for now.

        2) Copies of F,L,Qc are created in memory because they may be used later
        in smoother. References to copies are kept in "AQcomp" object
        return parameter.

        3) Function support "multiple time series mode" which means that exactly
        the same State-Space model is used to filter several sets of measurements.
        In this case third dimension of Y should include these state-space measurements
        Log_likelihood and Grad_log_likelihood have the corresponding dimensions then.

        4) Calculation of Grad_log_likelihood is not supported if matrices
        H, or R changes overf time (with index k). (later may be changed)

        5) Measurement may include missing values. In this case update step is
        not done for this measurement. (later may be changed)

        Input:
        -----------------

        F: (state_dim, state_dim) matrix
            F in the model.

        L: (state_dim, noise_dim) matrix
            L in the model.

        Qc: (noise_dim, noise_dim) matrix
            Q_c in the model.

        p_H: scalar, matrix (measurement_dim, state_dim) , 3D array
            H_{k} in the model. If matrix then H_{k} = H - constant.
            If it is 3D array then H_{k} = p_Q[:,:, index[2,k]]

        p_R: scalar, square symmetric matrix, 3D array
            R_{k} in the model. If matrix then R_{k} = R - constant.
            If it is 3D array then R_{k} = p_R[:,:, index[3,k]]

        P_inf: (state_dim, state_dim) matrix
            State varince matrix on infinity.

        X: 1D array
            Time points of measurements. Needed for converting continuos
            problem to the discrete one.

        Y: matrix or vector or 3D array
            Data. If Y is matrix then samples are along 0-th dimension and
            features along the 1-st. If 3D array then third dimension
            correspond to "multiple time series mode".

        index: vector
            Which indices (on 3-rd dimension) from arrays p_H, p_R to use
            on every time step. If this parameter is None then it is assumed
            that p_H, p_R do not change over time and indices are not needed.
            index[0,:] - correspond to H, index[1,:] - correspond to R
            If index.shape[0] == 1, it is assumed that indides for all matrices
            are the same.

        m_init: vector or matrix
            Initial distribution mean. If None it is assumed to be zero.
            For "multiple time series mode" it is matrix, second dimension of
            which correspond to different time series. In regular case ("one
            time series mode") it is a vector.

        P_init: square symmetric matrix or scalar
            Initial covariance of the states. If the parameter is scalar
            then it is assumed that initial covariance matrix is unit matrix
            multiplied by this scalar. If None the unit matrix is used instead.
            "multiple time series mode" does not affect it, since it does not
            affect anything related to state variaces.

        p_kalman_filter_type: string, one of ('regular', 'svd')
            Which Kalman Filter is used. Regular or SVD. SVD is more numerically
            stable, in particular, Covariace matrices are guarantied to be
            positive semi-definite. However, 'svd' works slower, especially for
            small data due to SVD call overhead.

        calc_log_likelihood: boolean
            Whether to calculate marginal likelihood of the state-space model.

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

        grad_params_no: int
            If previous parameter is true, then this parameters gives the
            total number of parameters in the gradient.

        grad_calc_params: dictionary
            Dictionary with derivatives of model matrices with respect
            to parameters "dF", "dL", "dQc", "dH", "dR", "dm_init", "dP_init".
            They can be None, in this case zero matrices (no dependence on parameters)
            is assumed. If there is only one parameter then third dimension is
            automatically added.

        Output:
        --------------

        M: (no_steps+1,state_dim) matrix or (no_steps+1,state_dim, time_series_no) 3D array
            Filter estimates of the state means. In the extra step the initial
            value is included. In the "multiple time series mode" third dimension
            correspond to different timeseries.

        P: (no_steps+1, state_dim, state_dim) 3D array
            Filter estimates of the state covariances. In the extra step the initial
            value is included.

        log_likelihood: double or (1, time_series_no) 3D array.

            If the parameter calc_log_likelihood was set to true, return
            logarithm of marginal likelihood of the state-space model. If
            the parameter was false, return None. In the "multiple time series mode" it is a vector
            providing log_likelihood for each time series.

        grad_log_likelihood: column vector or (grad_params_no, time_series_no) matrix
            If calc_grad_log_likelihood is true, return gradient of log likelihood
            with respect to parameters. It returns it column wise, so in
            "multiple time series mode" gradients for each time series is in the
            corresponding column.

        AQcomp: object
            Contains some pre-computed values for converting continuos model into
            discrete one. It can be used later in the smoothing pahse.

def cont_discr_kalman_filter(cls, F, L, Qc, p_H, p_R, P_inf, X, Y, index = None,
                                 m_init=None, P_init=None,
                                 p_kalman_filter_type='regular',
                                 calc_log_likelihood=False,
                                 calc_grad_log_likelihood=False,
                                 grad_params_no=0, grad_calc_params=None):
        """
        This function implements the continuous-discrete Kalman Filter algorithm
        These notations for the State-Space model are assumed:
            d/dt x(t) =  F * x(t) + L * w(t);         w(t) ~ N(0, Qc)
            y_{k} = H_{k} * x_{k} + r_{k};            r_{k-1} ~ N(0, R_{k})

        Returns estimated filter distributions x_{k} ~ N(m_{k}, P(k))

        Current Features:
        ----------------------------------------
        1) The function generaly do not modify the passed parameters. If
        it happens then it is an error. There are several exeprions: scalars
        can be modified into a matrix, in some rare cases shapes of
        the derivatives matrices may be changed, it is ignored for now.

        2) Copies of F,L,Qc are created in memory because they may be used later
        in smoother. References to copies are kept in "AQcomp" object
        return parameter.

        3) Function support "multiple time series mode" which means that exactly
        the same State-Space model is used to filter several sets of measurements.
        In this case third dimension of Y should include these state-space measurements
        Log_likelihood and Grad_log_likelihood have the corresponding dimensions then.

        4) Calculation of Grad_log_likelihood is not supported if matrices
        H, or R changes overf time (with index k). (later may be changed)

        5) Measurement may include missing values. In this case update step is
        not done for this measurement. (later may be changed)

        Input:
        -----------------

        F: (state_dim, state_dim) matrix
            F in the model.

        L: (state_dim, noise_dim) matrix
            L in the model.

        Qc: (noise_dim, noise_dim) matrix
            Q_c in the model.

        p_H: scalar, matrix (measurement_dim, state_dim) , 3D array
            H_{k} in the model. If matrix then H_{k} = H - constant.
            If it is 3D array then H_{k} = p_Q[:,:, index[2,k]]

        p_R: scalar, square symmetric matrix, 3D array
            R_{k} in the model. If matrix then R_{k} = R - constant.
            If it is 3D array then R_{k} = p_R[:,:, index[3,k]]

        P_inf: (state_dim, state_dim) matrix
            State varince matrix on infinity.

        X: 1D array
            Time points of measurements. Needed for converting continuos
            problem to the discrete one.

        Y: matrix or vector or 3D array
            Data. If Y is matrix then samples are along 0-th dimension and
            features along the 1-st. If 3D array then third dimension
            correspond to "multiple time series mode".

        index: vector
            Which indices (on 3-rd dimension) from arrays p_H, p_R to use
            on every time step. If this parameter is None then it is assumed
            that p_H, p_R do not change over time and indices are not needed.
            index[0,:] - correspond to H, index[1,:] - correspond to R
            If index.shape[0] == 1, it is assumed that indides for all matrices
            are the same.

        m_init: vector or matrix
            Initial distribution mean. If None it is assumed to be zero.
            For "multiple time series mode" it is matrix, second dimension of
            which correspond to different time series. In regular case ("one
            time series mode") it is a vector.

        P_init: square symmetric matrix or scalar
            Initial covariance of the states. If the parameter is scalar
            then it is assumed that initial covariance matrix is unit matrix
            multiplied by this scalar. If None the unit matrix is used instead.
            "multiple time series mode" does not affect it, since it does not
            affect anything related to state variaces.

        p_kalman_filter_type: string, one of ('regular', 'svd')
            Which Kalman Filter is used. Regular or SVD. SVD is more numerically
            stable, in particular, Covariace matrices are guarantied to be
            positive semi-definite. However, 'svd' works slower, especially for
            small data due to SVD call overhead.

        calc_log_likelihood: boolean
            Whether to calculate marginal likelihood of the state-space model.

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

        grad_params_no: int
            If previous parameter is true, then this parameters gives the
            total number of parameters in the gradient.

        grad_calc_params: dictionary
            Dictionary with derivatives of model matrices with respect
            to parameters "dF", "dL", "dQc", "dH", "dR", "dm_init", "dP_init".
            They can be None, in this case zero matrices (no dependence on parameters)
            is assumed. If there is only one parameter then third dimension is
            automatically added.

        Output:
        --------------

        M: (no_steps+1,state_dim) matrix or (no_steps+1,state_dim, time_series_no) 3D array
            Filter estimates of the state means. In the extra step the initial
            value is included. In the "multiple time series mode" third dimension
            correspond to different timeseries.

        P: (no_steps+1, state_dim, state_dim) 3D array
            Filter estimates of the state covariances. In the extra step the initial
            value is included.

        log_likelihood: double or (1, time_series_no) 3D array.

            If the parameter calc_log_likelihood was set to true, return
            logarithm of marginal likelihood of the state-space model. If
            the parameter was false, return None. In the "multiple time series mode" it is a vector
            providing log_likelihood for each time series.

        grad_log_likelihood: column vector or (grad_params_no, time_series_no) matrix
            If calc_grad_log_likelihood is true, return gradient of log likelihood
            with respect to parameters. It returns it column wise, so in
            "multiple time series mode" gradients for each time series is in the
            corresponding column.

        AQcomp: object
            Contains some pre-computed values for converting continuos model into
            discrete one. It can be used later in the smoothing pahse.
        """

        p_H = np.atleast_1d(p_H)
        p_R = np.atleast_1d(p_R)

        X.shape, old_X_shape  = cls._reshape_input_data(X.shape, 2) # represent as column
        if (X.shape[1] != 1):
            raise ValueError("Only one dimensional X data is supported.")

        Y.shape, old_Y_shape  = cls._reshape_input_data(Y.shape) # represent as column

        state_dim = F.shape[0]
        measurement_dim = Y.shape[1]
        time_series_no = Y.shape[2] # multiple time series mode

        if  ((len(p_H.shape) == 3) and (len(p_H.shape[2]) != 1)) or\
            ((len(p_R.shape) == 3) and (len(p_R.shape[2]) != 1)):
            model_matrices_chage_with_time = True
        else:
            model_matrices_chage_with_time = False

        # Check index
        old_index_shape = None
        if index is None:
            if (len(p_H.shape) == 3) or (len(p_R.shape) == 3):
                raise ValueError("Parameter index can not be None for time varying matrices (third dimension is present)")
            else: # matrices do not change in time, so form dummy zero indices.
                index = np.zeros((1,Y.shape[0]))
        else:
            if len(index.shape) == 1:
                index.shape = (1,index.shape[0])
                old_index_shape = (index.shape[0],)

            if (index.shape[1] != Y.shape[0]):
                raise ValueError("Number of measurements must be equal the number of H_{k}, R_{k}")

        if (index.shape[0] == 1):
            H_time_var_index = 0; R_time_var_index = 0
        elif (index.shape[0] == 4):
            H_time_var_index = 0; R_time_var_index = 1
        else:
            raise ValueError("First Dimension of index must be either 1 or 2.")

        (p_H, old_H_shape) = cls._check_SS_matrix(p_H, state_dim, measurement_dim, which='H')
        (p_R, old_R_shape) = cls._check_SS_matrix(p_R, state_dim, measurement_dim, which='R')

        if m_init is None:
            m_init = np.zeros((state_dim, time_series_no))
        else:
            m_init = np.atleast_2d(m_init).T

        if P_init is None:
            P_init = P_inf.copy()

        if p_kalman_filter_type not in ('regular', 'svd'):
            raise ValueError("Kalman filer type neither 'regular nor 'svd'.")

        # Functions to pass to the kalman_filter algorithm:
        # Parameters:
        # k - number of Kalman filter iteration
        # m - vector for calculating matrices. Required for EKF. Not used here.
        # f_hl = lambda k,m,H: np.dot(H, m)
        # f_H = lambda k,m,P: p_H[:,:, index[H_time_var_index, k]]
        #f_R = lambda k: p_R[:,:, index[R_time_var_index, k]]
        #o_R = R_handling( p_R, index, R_time_var_index, 20)

        if calc_grad_log_likelihood:

            dF = cls._check_grad_state_matrices(grad_calc_params.get('dF'), state_dim, grad_params_no, which = 'dA')
            dQc = cls._check_grad_state_matrices(grad_calc_params.get('dQc'), state_dim, grad_params_no, which = 'dQ')
            dP_inf = cls._check_grad_state_matrices(grad_calc_params.get('dP_inf'), state_dim, grad_params_no, which = 'dA')

            dH = cls._check_grad_measurement_matrices(grad_calc_params.get('dH'), state_dim, grad_params_no, measurement_dim, which = 'dH')
            dR = cls._check_grad_measurement_matrices(grad_calc_params.get('dR'), state_dim, grad_params_no, measurement_dim, which = 'dR')

            dm_init = grad_calc_params.get('dm_init') # Initial values for the Kalman Filter
            if dm_init is None:
                # multiple time series mode. Keep grad_params always as a last dimension
                dm_init = np.zeros( (state_dim, time_series_no, grad_params_no) )

            dP_init = grad_calc_params.get('dP_init') # Initial values for the Kalman Filter
            if dP_init is None:
                dP_init = dP_inf(0).copy() # get the dP_init matrix, because now it is a function

        else:
            dP_inf = None
            dF = None
            dQc = None
            dH = None
            dR = None
            dm_init = None
            dP_init = None

        measurement_callables = Std_Measurement_Callables_Class(p_H, H_time_var_index, p_R, index, R_time_var_index, 20, dH, dR)
        #import pdb; pdb.set_trace()

        dynamic_callables = cls._cont_to_discrete_object(X, F, L, Qc, compute_derivatives=calc_grad_log_likelihood,
                                              grad_params_no=grad_params_no,
                                              P_inf=P_inf, dP_inf=dP_inf, dF = dF, dQc=dQc)

        if print_verbose:
            print("General: run Continuos-Discrete Kalman Filter")
        # Also for dH, dR and probably for all derivatives
        (M, P, log_likelihood, grad_log_likelihood, AQcomp) = cls._cont_discr_kalman_filter_raw(state_dim,
                        dynamic_callables, measurement_callables,
                        X, Y, m_init=m_init, P_init=P_init,
                        p_kalman_filter_type=p_kalman_filter_type,
                        calc_log_likelihood=calc_log_likelihood,
                        calc_grad_log_likelihood=calc_grad_log_likelihood, grad_params_no=grad_params_no,
                        dm_init=dm_init, dP_init=dP_init)

        if old_index_shape is not None:
            index.shape = old_index_shape

        if old_X_shape is not None:
            X.shape = old_X_shape

        if old_Y_shape is not None:
            Y.shape = old_Y_shape

        if old_H_shape is not None:
            p_H.shape = old_H_shape

        if old_R_shape is not None:
            p_R.shape = old_R_shape

        return (M, P, log_likelihood, grad_log_likelihood, AQcomp)