<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def eda_scr(signal, sampling_rate=1000, treshold=0.1, method="fast"):
 """ Skin-Conductance Responses extraction algorithm. Parameters signal : list or array EDA signal array. sampling_rate : int Sampling rate (samples/second). treshold : float SCR minimum treshold (in terms of signal standart deviation). method : str "fast" or "slow". Either use a gradient-based approach or a local extrema one. Returns onsets, peaks, amplitudes, recoveries : lists SCRs features. Example Notes *Authors* - `Dominique Makowski <https://dominiquemakowski.github.io/>`_ *Dependencies* - biosppy - numpy - pandas *See Also* - BioSPPy: https://github.com/PIA-Group/BioSPPy References - Kim, K. H., Bang, S. W., & Kim, S. R. (2004). Emotion recognition system using short-term monitoring of physiological signals. Medical and biological engineering and computing, 42(3), 419-427. - Gamboa, H. (2008). Multi-Modal Behavioral Biometrics Based on HCI and Electrophysiology (Doctoral dissertation, PhD thesis, Universidade Técnica de Lisboa, Instituto Superior Técnico). """

    # Processing
    # ===========
    if method == "slow":
        # Compute gradient (sort of derivative)
        gradient = np.gradient(signal)
        # Smoothing
        size = int(0.1 * sampling_rate)
        smooth, _ = biosppy.tools.smoother(signal=gradient, kernel='bartlett', size=size, mirror=True)
        # Find zero-crossings
        zeros, = biosppy.tools.zero_cross(signal=smooth, detrend=True)

        # Separate onsets and peaks
        onsets = []
        peaks = []
        for i in zeros:
            if smooth[i+1] > smooth[i-1]:
                onsets.append(i)
            else:
                peaks.append(i)
        peaks = np.array(peaks)
        onsets = np.array(onsets)
    else:
        # find extrema
        peaks, _ = biosppy.tools.find_extrema(signal=signal, mode='max')
        onsets, _ = biosppy.tools.find_extrema(signal=signal, mode='min')

    # Keep only pairs
    peaks = peaks[peaks > onsets[0]]
    onsets = onsets[onsets < peaks[-1]]

    # Artifact Treatment
    # ====================
    # Compute rising times
    risingtimes = peaks-onsets
    risingtimes = risingtimes/sampling_rate*1000

    peaks = peaks[risingtimes > 100]
    onsets = onsets[risingtimes > 100]

    # Compute amplitudes
    amplitudes = signal[peaks]-signal[onsets]

    # Remove low amplitude variations
    mask = amplitudes > np.std(signal)*treshold
    peaks = peaks[mask]
    onsets = onsets[mask]
    amplitudes = amplitudes[mask]

    # Recovery moments
    recoveries = []
    for x, peak in enumerate(peaks):
        try:
            window = signal[peak:onsets[x+1]]
        except IndexError:
            window = signal[peak:]
        recovery_amp = signal[peak]-amplitudes[x]/2
        try:
            smaller = find_closest_in_list(recovery_amp, window, "smaller")
            recovery_pos = peak + list(window).index(smaller)
            recoveries.append(recovery_pos)
        except ValueError:
            recoveries.append(np.nan)
    recoveries = np.array(recoveries)

    return(onsets, peaks, amplitudes, recoveries)