# Path: python/pipeline/utils/caiman_stats.py
# def df_percentile(inputData, axis=None):
#     """
#     Extracting the percentile of the data where the mode occurs and its value.
#     Used to determine the filtering level for DF/F extraction. Note that
#     computation can be innacurate for short traces.
#     """
#     if axis is not None:
# 
#         def fnc(x):
#             return df_percentile(x)
# 
#         result = np.apply_along_axis(fnc, axis, inputData)
#         data_prct = result[:, 0]
#         val = result[:, 1]
#     else:
#         # Create the function that we can use for the half-sample mode
#         err = True
#         while err:
#             try:
#                 bandwidth, mesh, density, cdf = kde(inputData)
#                 err = False
#             except:
#                 logging.warning('Percentile computation failed. Duplicating ' + 'and trying again.')
#                 if type(inputData) is not list:
#                     inputData = inputData.tolist()
#                 inputData += inputData
# 
#         data_prct = cdf[np.argmax(density)] * 100
#         val = mesh[np.argmax(density)]
#         if data_prct >= 100 or data_prct < 0:
#             logging.warning('Invalid percentile computed possibly due ' + 'short trace. Duplicating and recomuputing.')
#             if type(inputData) is not list:
#                 inputData = inputData.tolist()
#             inputData *= 2
#             err = True
#         if np.isnan(data_prct):
#             logging.warning('NaN percentile computed. Reverting to median.')
#             data_prct = 50
#             val = np.median(np.array(inputData))
# 
#     return data_prct, val


import numpy as np
import multiprocessing as mp
import glob, os, sys, time
    import matplotlib.pyplot as plt; plt.close()
    import caiman
from caiman import components_evaluation
from caiman.utils import visualization
from caiman.source_extraction.cnmf import map_reduce, initialization, pre_processing, \
                                          merging, spatial, temporal, deconvolution
from .caiman_stats import df_percentile
from scipy.ndimage import percentile_filter
        from scipy.sparse import csr_matrix
    from scipy import ndimage
        from sklearn.decomposition import NMF
    from matplotlib import mlab
    from scipy.sparse import coo_matrix
    from caiman.source_extraction.cnmf import cnmf
    :returns: Deconvolved spike trace.
    :returns: AR coefficients (AR_order) that model the calcium response:
            c(t) = c(t-1) * AR_coeffs[0] + c(t-2) * AR_coeffs[1] + ...
    """
    _, _, _, AR_coeffs, _, spike_trace, _ = deconvolution.constrained_foopsi(trace,
        p=AR_order, method='cvxpy', bas_nonneg=False, fudge_factor=0.96)
        # fudge_factor is a regularization term

    return spike_trace, AR_coeffs


def deconvolve_detrended(trace, scan_fps, detrend_period=600, AR_order=2):
    """Same as the the `deconvolve` method, except that the fluorescence trace is detrended 
    before autoregressive modeling

    :param np.array trace: 1-d array (num_frames) with the fluorescence trace.
    :param float scan_fps: fps of the scan
    :param float detrend_period: number of seconds over which percentiles are computed
    :param int AR_order: Order of the autoregressive process used to model the impulse
        response function, e.g., 0 = no modelling; 2 = model rise plus exponential decay.

    :returns: Deconvolved spike trace.
    :returns: AR coefficients (AR_order) that model the calcium response:
            c(t) = c(t-1) * AR_coeffs[0] + c(t-2) * AR_coeffs[1] + ...
    """
    detrend_window = int(round(detrend_period * scan_fps))
    n_chunks = len(trace) // detrend_window
    if detrend_window > 0 and n_chunks > 0:
        chunks_len = n_chunks * detrend_window
        trace_chunks = trace[:chunks_len].reshape(-1, detrend_window)

# Path: python/pipeline/utils/caiman_stats.py
# def df_percentile(inputData, axis=None):
#     """
#     Extracting the percentile of the data where the mode occurs and its value.
#     Used to determine the filtering level for DF/F extraction. Note that
#     computation can be innacurate for short traces.
#     """
#     if axis is not None:
# 
#         def fnc(x):
#             return df_percentile(x)
# 
#         result = np.apply_along_axis(fnc, axis, inputData)
#         data_prct = result[:, 0]
#         val = result[:, 1]
#     else:
#         # Create the function that we can use for the half-sample mode
#         err = True
#         while err:
#             try:
#                 bandwidth, mesh, density, cdf = kde(inputData)
#                 err = False
#             except:
#                 logging.warning('Percentile computation failed. Duplicating ' + 'and trying again.')
#                 if type(inputData) is not list:
#                     inputData = inputData.tolist()
#                 inputData += inputData
# 
#         data_prct = cdf[np.argmax(density)] * 100
#         val = mesh[np.argmax(density)]
#         if data_prct >= 100 or data_prct < 0:
#             logging.warning('Invalid percentile computed possibly due ' + 'short trace. Duplicating and recomuputing.')
#             if type(inputData) is not list:
#                 inputData = inputData.tolist()
#             inputData *= 2
#             err = True
#         if np.isnan(data_prct):
#             logging.warning('NaN percentile computed. Reverting to median.')
#             data_prct = 50
#             val = np.median(np.array(inputData))
# 
#     return data_prct, val


# Path: python/pipeline/utils/caiman_interface.py
import numpy as np
import multiprocessing as mp
import glob, os, sys, time
    import matplotlib.pyplot as plt; plt.close()
    import caiman
from caiman import components_evaluation
from caiman.utils import visualization
from caiman.source_extraction.cnmf import map_reduce, initialization, pre_processing, \
                                          merging, spatial, temporal, deconvolution
from .caiman_stats import df_percentile
from scipy.ndimage import percentile_filter
        from scipy.sparse import csr_matrix
    from scipy import ndimage
        from sklearn.decomposition import NMF
    from matplotlib import mlab
    from scipy.sparse import coo_matrix
    from caiman.source_extraction.cnmf import cnmf

    :returns: Deconvolved spike trace.
    :returns: AR coefficients (AR_order) that model the calcium response:
            c(t) = c(t-1) * AR_coeffs[0] + c(t-2) * AR_coeffs[1] + ...
    """
    _, _, _, AR_coeffs, _, spike_trace, _ = deconvolution.constrained_foopsi(trace,
        p=AR_order, method='cvxpy', bas_nonneg=False, fudge_factor=0.96)
        # fudge_factor is a regularization term

    return spike_trace, AR_coeffs


def deconvolve_detrended(trace, scan_fps, detrend_period=600, AR_order=2):
    """Same as the the `deconvolve` method, except that the fluorescence trace is detrended 
    before autoregressive modeling

    :param np.array trace: 1-d array (num_frames) with the fluorescence trace.
    :param float scan_fps: fps of the scan
    :param float detrend_period: number of seconds over which percentiles are computed
    :param int AR_order: Order of the autoregressive process used to model the impulse
        response function, e.g., 0 = no modelling; 2 = model rise plus exponential decay.

    :returns: Deconvolved spike trace.
    :returns: AR coefficients (AR_order) that model the calcium response:
            c(t) = c(t-1) * AR_coeffs[0] + c(t-2) * AR_coeffs[1] + ...
    """
    detrend_window = int(round(detrend_period * scan_fps))
    n_chunks = len(trace) // detrend_window
    if detrend_window > 0 and n_chunks > 0:
        chunks_len = n_chunks * detrend_window
        trace_chunks = trace[:chunks_len].reshape(-1, detrend_window)

        data_prct = df_percentile(trace_chunks, axis=1)[0].mean()