From 850a2fdad81b1c24a0ab57766f5048e363108541 Mon Sep 17 00:00:00 2001 From: JakHolzer <53743814+JakHolzer@users.noreply.github.com> Date: Thu, 10 Sep 2020 15:14:55 +0200 Subject: [PATCH 1/4] Add files via upload function to find peaks (more detailed description in email) --- pyzebra/ccl_findpeaks.py | 81 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) create mode 100644 pyzebra/ccl_findpeaks.py diff --git a/pyzebra/ccl_findpeaks.py b/pyzebra/ccl_findpeaks.py new file mode 100644 index 0000000..a0d2fd6 --- /dev/null +++ b/pyzebra/ccl_findpeaks.py @@ -0,0 +1,81 @@ +import numpy as np +import scipy as sc +from scipy.interpolate import interp1d +from scipy.signal import savgol_filter + + +def ccl_findpeaks(data, keys, int_threshold=None, prominence=None, smooth=True, window_size=None, poly_order=None): + + """ function iterates through the dictionary created by load_cclv2 and locates peaks for each measurement + args: data (dictionary from load_cclv2), + + int_threshold - fraction of threshold_intensity/max_intensity, must be positive num between 0 and 1 + i.e. will only detect peaks above 75% of max intensity + + prominence - defines a drop of values that must be between two peaks, must be positive number + i.e. if promimence is 20, it will detect two neigbouring peaks of 300 and 310 intesities, if none of the itermediate values are lower that 290 + + smooth - if true, smooths data by savitzky golay filter, if false - no smoothing + + window_size - window size for savgol filter, must be odd positive integer + + poly_order = order of the polynomial used in savgol filter, must be positive integer smaller than window_size + returns: dictionary with following structure: + D{M34{ 'num_of_peaks': 1, #num of peaks + 'peak_indexes': [20], # index of peaks in omega array + 'peak_heights': [90.], # height of the peaks (if data vere smoothed its the heigh of the peaks in smoothed data) + + + """ + + if type(data) is dict and data["file_type"] == 'ccl': + int_threshold = 0.75 if int_threshold is None else int_threshold + prominence = 50 if prominence is None else prominence + smooth = False if smooth is None else smooth + window_size = 7 if window_size is None else window_size + poly_order = 3 if poly_order is None else poly_order + + if 0 <= int_threshold <= 1: + pass + else: + int_threshold = 0.75 + print('Invalid value for int_threshold, select value between 0 and 1, new value set to:', int_threshold) + if isinstance(window_size, int) is True and (window_size % 2) != 0 and window_size >= 1: + pass + else: + window_size = 7 + print('Invalid value for window_size, select positive odd integer, new value set to:', window_size) + if isinstance(poly_order, int) is True and window_size > poly_order >= 1: + pass + else: + poly_order = 3 + print('Invalid value for poly_order, select positive integer smaller than window_size, new value set to:', poly_order) + if isinstance(prominence, (int, float)) is True and prominence > 0: + pass + else: + prominence = 50 + print('Invalid value for prominence, select positive number, new value set to:', + prominence) + + omega = data["Measurements"][str(keys)]["omega"] + counts = np.array(data["Measurements"][str(keys)]["counts"]) + if smooth is True: + itp = interp1d(omega, counts, kind='linear') + absintensity = [abs(number) for number in counts] + lowest_intensity = min(absintensity) + counts[counts < 0] = lowest_intensity + smooth_peaks = savgol_filter(itp(omega), window_size, poly_order) + + else: + smooth_peaks = counts + + indexes = sc.signal.find_peaks(smooth_peaks, height=int_threshold*max(smooth_peaks), prominence=prominence) + data["Measurements"][str(keys)]["num_of_peaks"] = len(indexes[0]) + data["Measurements"][str(keys)]["peak_indexes"] = indexes[0] + data["Measurements"][str(keys)]["peak_heights"] = indexes[1]["peak_heights"] + data["Measurements"][str(keys)]["smooth_peaks"] = smooth_peaks # smoothed curve + + return data + + else: + return print('Data is not a dictionary or was not made from ccl file') \ No newline at end of file -- 2.47.2 From f9ced1cebc6e1ca7f1ec71aff281179c6fde337c Mon Sep 17 00:00:00 2001 From: JakHolzer <53743814+JakHolzer@users.noreply.github.com> Date: Fri, 11 Sep 2020 15:30:30 +0200 Subject: [PATCH 2/4] guard close changed --- pyzebra/ccl_findpeaks.py | 121 +++++++++++++++++++++------------------ 1 file changed, 66 insertions(+), 55 deletions(-) diff --git a/pyzebra/ccl_findpeaks.py b/pyzebra/ccl_findpeaks.py index a0d2fd6..a8df1d4 100644 --- a/pyzebra/ccl_findpeaks.py +++ b/pyzebra/ccl_findpeaks.py @@ -4,78 +4,89 @@ from scipy.interpolate import interp1d from scipy.signal import savgol_filter -def ccl_findpeaks(data, keys, int_threshold=None, prominence=None, smooth=True, window_size=None, poly_order=None): +def ccl_findpeaks( + data, keys, int_threshold=None, prominence=None, smooth=True, window_size=None, poly_order=None +): - """ function iterates through the dictionary created by load_cclv2 and locates peaks for each measurement + """function iterates through the dictionary created by load_cclv2 and locates peaks for each measurement args: data (dictionary from load_cclv2), int_threshold - fraction of threshold_intensity/max_intensity, must be positive num between 0 and 1 i.e. will only detect peaks above 75% of max intensity prominence - defines a drop of values that must be between two peaks, must be positive number - i.e. if promimence is 20, it will detect two neigbouring peaks of 300 and 310 intesities, if none of the itermediate values are lower that 290 + i.e. if promimence is 20, it will detect two neigbouring peaks of 300 and 310 intesities, + if none of the itermediate values are lower that 290 smooth - if true, smooths data by savitzky golay filter, if false - no smoothing window_size - window size for savgol filter, must be odd positive integer - poly_order = order of the polynomial used in savgol filter, must be positive integer smaller than window_size - returns: dictionary with following structure: + poly_order = order of the polynomial used in savgol filter, must be positive integer smaller than + window_size returns: dictionary with following structure: D{M34{ 'num_of_peaks': 1, #num of peaks 'peak_indexes': [20], # index of peaks in omega array - 'peak_heights': [90.], # height of the peaks (if data vere smoothed its the heigh of the peaks in smoothed data) + 'peak_heights': [90.], # height of the peaks (if data vere smoothed + its the heigh of the peaks in smoothed data) + """ + if type(data) is not dict and data["file_type"] != "ccl": + print("Data is not a dictionary or was not made from ccl file") - """ + int_threshold = 0.75 if int_threshold is None else int_threshold + prominence = 50 if prominence is None else prominence + smooth = False if smooth is None else smooth + window_size = 7 if window_size is None else window_size + poly_order = 3 if poly_order is None else poly_order - if type(data) is dict and data["file_type"] == 'ccl': - int_threshold = 0.75 if int_threshold is None else int_threshold - prominence = 50 if prominence is None else prominence - smooth = False if smooth is None else smooth - window_size = 7 if window_size is None else window_size - poly_order = 3 if poly_order is None else poly_order + if 0 <= int_threshold <= 1: + pass + else: + int_threshold = 0.75 + print( + "Invalid value for int_threshold, select value between 0 and 1, new value set to:", + int_threshold, + ) + if isinstance(window_size, int) is True and (window_size % 2) != 0 and window_size >= 1: + pass + else: + window_size = 7 + print( + "Invalid value for window_size, select positive odd integer, new value set to:", + window_size, + ) + if isinstance(poly_order, int) is True and window_size > poly_order >= 1: + pass + else: + poly_order = 3 + print( + "Invalid value for poly_order, select positive integer smaller than window_size, new value set to:", + poly_order, + ) + if isinstance(prominence, (int, float)) is True and prominence > 0: + pass + else: + prominence = 50 + print("Invalid value for prominence, select positive number, new value set to:", prominence) - if 0 <= int_threshold <= 1: - pass - else: - int_threshold = 0.75 - print('Invalid value for int_threshold, select value between 0 and 1, new value set to:', int_threshold) - if isinstance(window_size, int) is True and (window_size % 2) != 0 and window_size >= 1: - pass - else: - window_size = 7 - print('Invalid value for window_size, select positive odd integer, new value set to:', window_size) - if isinstance(poly_order, int) is True and window_size > poly_order >= 1: - pass - else: - poly_order = 3 - print('Invalid value for poly_order, select positive integer smaller than window_size, new value set to:', poly_order) - if isinstance(prominence, (int, float)) is True and prominence > 0: - pass - else: - prominence = 50 - print('Invalid value for prominence, select positive number, new value set to:', - prominence) - - omega = data["Measurements"][str(keys)]["omega"] - counts = np.array(data["Measurements"][str(keys)]["counts"]) - if smooth is True: - itp = interp1d(omega, counts, kind='linear') - absintensity = [abs(number) for number in counts] - lowest_intensity = min(absintensity) - counts[counts < 0] = lowest_intensity - smooth_peaks = savgol_filter(itp(omega), window_size, poly_order) - - else: - smooth_peaks = counts - - indexes = sc.signal.find_peaks(smooth_peaks, height=int_threshold*max(smooth_peaks), prominence=prominence) - data["Measurements"][str(keys)]["num_of_peaks"] = len(indexes[0]) - data["Measurements"][str(keys)]["peak_indexes"] = indexes[0] - data["Measurements"][str(keys)]["peak_heights"] = indexes[1]["peak_heights"] - data["Measurements"][str(keys)]["smooth_peaks"] = smooth_peaks # smoothed curve - - return data + omega = data["Measurements"][str(keys)]["omega"] + counts = np.array(data["Measurements"][str(keys)]["counts"]) + if smooth is True: + itp = interp1d(omega, counts, kind="linear") + absintensity = [abs(number) for number in counts] + lowest_intensity = min(absintensity) + counts[counts < 0] = lowest_intensity + smooth_peaks = savgol_filter(itp(omega), window_size, poly_order) else: - return print('Data is not a dictionary or was not made from ccl file') \ No newline at end of file + smooth_peaks = counts + + indexes = sc.signal.find_peaks( + smooth_peaks, height=int_threshold * max(smooth_peaks), prominence=prominence + ) + data["Measurements"][str(keys)]["num_of_peaks"] = len(indexes[0]) + data["Measurements"][str(keys)]["peak_indexes"] = indexes[0] + data["Measurements"][str(keys)]["peak_heights"] = indexes[1]["peak_heights"] + data["Measurements"][str(keys)]["smooth_peaks"] = smooth_peaks # smoothed curve + + return data -- 2.47.2 From 90c02565957a23a58135c00fe6a96a71a43b37ad Mon Sep 17 00:00:00 2001 From: JakHolzer <53743814+JakHolzer@users.noreply.github.com> Date: Fri, 11 Sep 2020 16:45:31 +0200 Subject: [PATCH 3/4] Update ccl_findpeaks.py --- pyzebra/ccl_findpeaks.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/pyzebra/ccl_findpeaks.py b/pyzebra/ccl_findpeaks.py index a8df1d4..6e5ba56 100644 --- a/pyzebra/ccl_findpeaks.py +++ b/pyzebra/ccl_findpeaks.py @@ -5,7 +5,7 @@ from scipy.signal import savgol_filter def ccl_findpeaks( - data, keys, int_threshold=None, prominence=None, smooth=True, window_size=None, poly_order=None + data, keys, int_threshold=0.8, prominence=50, smooth=False, window_size=7, poly_order=3 ): """function iterates through the dictionary created by load_cclv2 and locates peaks for each measurement @@ -33,12 +33,6 @@ def ccl_findpeaks( if type(data) is not dict and data["file_type"] != "ccl": print("Data is not a dictionary or was not made from ccl file") - int_threshold = 0.75 if int_threshold is None else int_threshold - prominence = 50 if prominence is None else prominence - smooth = False if smooth is None else smooth - window_size = 7 if window_size is None else window_size - poly_order = 3 if poly_order is None else poly_order - if 0 <= int_threshold <= 1: pass else: -- 2.47.2 From d30b348425e259c8d814fa1e784da6575062a25c Mon Sep 17 00:00:00 2001 From: JakHolzer <53743814+JakHolzer@users.noreply.github.com> Date: Mon, 14 Sep 2020 09:42:14 +0200 Subject: [PATCH 4/4] Update ccl_findpeaks.py --- pyzebra/ccl_findpeaks.py | 33 ++++++++++++++------------------- 1 file changed, 14 insertions(+), 19 deletions(-) diff --git a/pyzebra/ccl_findpeaks.py b/pyzebra/ccl_findpeaks.py index 6e5ba56..fc52425 100644 --- a/pyzebra/ccl_findpeaks.py +++ b/pyzebra/ccl_findpeaks.py @@ -15,54 +15,49 @@ def ccl_findpeaks( i.e. will only detect peaks above 75% of max intensity prominence - defines a drop of values that must be between two peaks, must be positive number - i.e. if promimence is 20, it will detect two neigbouring peaks of 300 and 310 intesities, + i.e. if promimence is 20, it will detect two neigbouring peaks of 300 and 310 intesities, if none of the itermediate values are lower that 290 smooth - if true, smooths data by savitzky golay filter, if false - no smoothing window_size - window size for savgol filter, must be odd positive integer - poly_order = order of the polynomial used in savgol filter, must be positive integer smaller than + poly_order = order of the polynomial used in savgol filter, must be positive integer smaller than window_size returns: dictionary with following structure: D{M34{ 'num_of_peaks': 1, #num of peaks 'peak_indexes': [20], # index of peaks in omega array - 'peak_heights': [90.], # height of the peaks (if data vere smoothed + 'peak_heights': [90.], # height of the peaks (if data vere smoothed its the heigh of the peaks in smoothed data) """ if type(data) is not dict and data["file_type"] != "ccl": print("Data is not a dictionary or was not made from ccl file") - if 0 <= int_threshold <= 1: - pass - else: - int_threshold = 0.75 + if not 0 <= int_threshold <= 1: + int_threshold = 0.8 print( "Invalid value for int_threshold, select value between 0 and 1, new value set to:", int_threshold, ) - if isinstance(window_size, int) is True and (window_size % 2) != 0 and window_size >= 1: - pass - else: + + if isinstance(window_size, int) is False or (window_size % 2) == 0 or window_size <= 1: window_size = 7 print( - "Invalid value for window_size, select positive odd integer, new value set to:", - window_size, - ) - if isinstance(poly_order, int) is True and window_size > poly_order >= 1: - pass - else: + "Invalid value for window_size, select positive odd integer, new value set to!:", + window_size) + + if isinstance(poly_order, int) is False or window_size < poly_order: poly_order = 3 print( "Invalid value for poly_order, select positive integer smaller than window_size, new value set to:", poly_order, ) - if isinstance(prominence, (int, float)) is True and prominence > 0: - pass - else: + + if isinstance(prominence, (int, float)) is False and prominence < 0: prominence = 50 print("Invalid value for prominence, select positive number, new value set to:", prominence) + omega = data["Measurements"][str(keys)]["omega"] counts = np.array(data["Measurements"][str(keys)]["counts"]) if smooth is True: -- 2.47.2