Skip to content

Instantly share code, notes, and snippets.

@fccoelho
Created February 11, 2010 14:04
Show Gist options
  • Save fccoelho/301530 to your computer and use it in GitHub Desktop.
Save fccoelho/301530 to your computer and use it in GitHub Desktop.
"""
MinMax peak detector
Originally from:
https://fedorahosted.org/pydarn/export/dfbf33a36aab3f576ec585f66b6bc99dfd5750d7/peakdetect/MinMax.py
"""
import sys, copy, math
import numpy
# default parameters
_DEFAULT_ANGLE_PTS = 3
_DEFAULT_ANGLE_THRESHOLD = 1.0
_DEFAULT_WING_WIDTH = 2
_DEFAULT_PEAK_FACTOR = 10
_DEFAULT_AREA_SCALE = 1
_DEFAULT_AREA_THRESHOLD = 1E-5
class peak:
def __init__(self):
self.index=None
self.left=None
self.right=None
self.value=None
self.area=None
self.fwhm=None
def minmax(ia, window,da=_DEFAULT_AREA_SCALE,threshold=_DEFAULT_AREA_THRESHOLD,mode="area"):
"""
@summary: MinMax peak detector
@param ic: An IonChromatogram object
@type ic: pyms.Peak.List.Class.Alignment
@param noise: A number. The noise estimate
@type noise: FloatType
@param window: A string. The width of the window for peak
detection specified as <NUMBER>s or <NUMBER>m, giving a time
in seconds or minutes, respectively
@type window: StringType
@param peak_factor: A number. The value that is multipled to the
noise level to determine the peak threshold for minimum peak
intensities
@type peak_factor: IntType
@return: A list.
@rtype: ListType
---------------------------------------------------------------
@summary: MinMax peak detector
@param ia A numpy object.
@param noise A number. The noise estimate.
@param window The width of the window for peak
detection, giving in points.
@param peak_factor A number. The value that is multipled to the
noise level to determine the peak threshold for minimum peak
intensities.
@return A list.
"""
# if noise < 0 or peak_factor < 0:
# error("negative argument(s) detected.")
# wing_length = window_sele_points(ic, window, half_window=True)
wing_length=window/2
# print " -> Minmax: using window wing of %d point(s)" % (wing_length)
peak_list = _get_peaks(ia, threshold, wing_length,da,mode)
p_max = max_peaks(ia,peak_list,da)
return peak_list,p_max
def _get_peaks(ia, threshold, wing_length,da,mode):
"""
@summary: Peak detection wrapper
@param ia: An numpy object. Ion chromatogram intensity array
@type ia: numpy.ndarray
@param noise: A number. The noise estimate
@type noise: FloatType
@param peak_factor: A number. The value that is multipled to
the noise level to determine the peak threshold for
minimum peak intensities
@type peak_factor: IntType
@param wing_length: An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value
@type wing_length: IntType
@return: A list of maxima points and their corresponding left and right
boundaries. Each maximum is modelled using a list containing threex
elements: the left boundary, maximum point, and right boundary
@rtype: ListType
------------------------------------------------------------------
Peak detection wrapper.
@param ia An numpy object. Ion chromatogram intensity array.
@param noise A number. The noise estimate.
@param peak_factor A number. The value that is multipled to
the noise level to determine the peak threshold for
minimum peak intensities.
@param wing_length An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value.
@return A list of maxima points and their corresponding left and right
boundaries. Each maximum is modelled using a list containing threex
elements: the left boundary, maximum point, and right boundary.
"""
# print " -> Minmax: noise: %lf peak factor: %lf" % (noise,peak_factor)
# peak_threshold = noise * peak_factor
# print " -> Minmax: prob threshold: %lf " % (threshold)
maxima_indices = _identify_maxima(ia, wing_length, threshold,da,mode)
# print " -> Minmax: found local maxima: %d point(s)" % (len(maxima_indices))
peak_list = _peak_boundaries(ia, maxima_indices, wing_length)
# print " -> Minmax: found peak boundaries: %d point(s)" % (len(peak_list))
_resolve_overlaps(ia, peak_list)
#_prune_peaks(peak_list, wing_length)
_adjust_boundaries(ia, peak_list, wing_length, 1.5 )
peak_area(ia, peak_list)
peak_fwhm(ia, peak_list)
return peak_list
def _prune_peaks(peak_list, wing_length):
"""
@summary: Rejects peaks for which the difference between the apex and
left/right boundary is less than wing_length
@param peak_list: A list of pre-peaks
@type peak_list: ListType
@param wing_length: An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value
@type wing_length: IntType
@return: No return value
@rtype: NoneType
-----------------------------------------------------------------
_prune_peaks(peak_list, wing_length)
Rejects peaks for which the difference between the apex and
left/right boundary is less than wing_length.
@param peak_list A list of pre-peaks.
@param wing_length An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value.
@return No return value.
"""
prune_peaks = []
for peak in peak_list:
peak_left = peak[LEFT_BOUNDARY]
peak_apex = peak[PEAK_APEX]
peak_right = peak[RIGHT_BOUNDARY]
if (peak_apex - peak_left) < wing_length or \
(peak_right - peak_apex) < wing_length:
prune_peaks.append(peak)
for peak in prune_peaks:
peak_list.remove(peak)
def _resolve_overlaps(ia, peak_list):
"""
@summary: Resolves overlapped peaks
@param ia: An numpy object. Ion chromatogram intensity array
@type ia: numpy.ndarray
@param peak_list: A list of pre-peaks
@type peak_list: ListType
@return: No return value
@rtype: NoneType
------------------------------------------------------------
_resolve_overlaps(ia, peak_list)
Resolves overlapped peaks.
@param ia An numpy object. Ion chromatogram intensity array.
@param peak_list A list of pre-peaks.
@return No return value.
"""
for ii in range(len(peak_list)-1):
peak_crnt = peak_list[ii]
peak_next = peak_list[ii+1]
if peak_crnt.right > peak_next.left:
_adjust_boundary(ia, peak_crnt, peak_next)
def _adjust_boundary(ia, peak_crnt, peak_next):
"""
@summary: Adjusts the boundary between two overlapping peaks. The new
boundary is set to the minimum between the two peak apexes
@param ia: An numpy object. Ion chromatogram intensity array
@type ia: numpy.ndarray
@param peak_crnt: A list [left,max,right] for peak 1
@type peak_crnt: ListType
@param peak_next: A list [left,max,right] for peak 2
@type peak_next: ListType
@return: No return value
@rtype: NoneType
-----------------------------------------------------------
_adjust_boundary(ia, peak_crnt, peak_next)
Adjusts the boundary between two overlapping peaks. The new
boundary is set to the minimum between the two peak apexes.
@param ia An numpy object. Ion chromatogram intensity array.
@param peak_crnt A list [left,max,right] for peak 1
@param peak_next A list [left,max,right] for peak 2
@return No return value.
"""
crnt_apex = peak_crnt.index
next_apex = peak_next.index
minx = numpy.argmin(ia[crnt_apex+1:next_apex])
minv = ia[crnt_apex+1:next_apex][minx]
split_point = crnt_apex + minx + 1
peak_crnt.right = split_point
peak_next.left = split_point
def _identify_maxima(ia, wing_length, peak_threshold,da,mode):
"""
@summary: Returns an array of indices where each index corresponds to
a local maximum. A point is a local maximum if the following
are satisfied:
1. It is equal or greater than any of 'wing_length' points to
its left and to its right
2. It is greater than at least one of the points to its left,
and one of the points to its right (out of total 'wing_length'
points)
3. A point closer to the edge of 'ia' than 'wing_length' cannot
be a local maximum
@param ia: An numpy object. Ion chromatogram intensity array
@type ia: numpy.ndarray
@param wing_length: An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value
@type wing_length: IntType
@param peak_threshold: A number. The minimum intensity for a peak,
maxima below this value are rejected
@type peak_threshold: FloatType
@return: An array of indices where each index corresponds to a local maximum
@rtype: ListType
--------------------------------------------------------------
_identify_maxima(ia, wing_length, peak_threshold)
Returns an array of indices where each index corresponds to
a local maximum. A point is a local maximum if the following
are satisfied:
1. It is equal or greater than any of 'wing_length' points to
its left and to its right.
2. It is greater than at least one of the points to its left,
and one of the points to its right (out of total 'wing_length'
points).
3. A point closer to the edge of 'ia' than 'wing_length' cannot
be a local maximum.
@param ia An numpy object. Ion chromatogram intensity array.
@param wing_length An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value.
@param peak_threshold A number. The minimum intensity for a peak,
maxima below this value are rejected.
"""
maxima_indices = []
index = 0
for index in xrange(ia.size):
is_maxima = _is_local_maximum(ia, index, wing_length)
if is_maxima:
pt_left = max(0,index-wing_length)
pt_right = min(ia.size-1,index+wing_length)
# print "--> Local maxima: %d %d %d" % (index,pt_left,pt_right)
# print "--> Local maxima: %d %lf %lf %lf" % (index,slice_area,da,peak_threshold)
if mode=="area":
slice_area=ia[pt_left:pt_right+1].sum()
area=slice_area*da
is_large_enough = area >= peak_threshold
else:
is_large_enough = ia[index] >= peak_threshold
if is_large_enough: # local maximum found
maxima_indices.append(index) # append to list
index = index + wing_length + 1 # index skip for wing_length
return maxima_indices
def _peak_boundaries(ia, maxima_indices, wing_length):
"""
@summary: Determines the left and right boundary for each of the maxima.
Returns a list of maxima points and their corresponding left and
right boundaries. Each maxima is modelled using a list of three
elements: the left boundary, maxium point, and right boundary
@param ia: An numpy object. Ion chromatogram intensity array
@type ia: numpy.ndarray
@param maxima_indices: A list contains numpy objects. The list indices of
the maxima points
@type maxima_indices: ListType
@param wing_length: An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value
@type wing_length: IntType
@return: A list of triplets. Giving indices for [left_bound, maximum,
right_bound]
@rtype: ListType
--------------------------------------------------------------------
_peak_boundaries(ia, maxima_indices, wing_length)
Determines the left and right boundary for each of the maxima.
Returns a list of maxima points and their corresponding left and
right boundaries. Each maxima is modelled using a list of three
elements: the left boundary, maxium point, and right boundary.
@param ia An numpy object. Ion chromatogram intensity array.
@param maxima_indices A numpy object. The list indices of
the maxima points.
@param wing_length An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value.
@return A list of triplets. Giving indices for [left_bound, maximum,
right_bound]
"""
peak_list = []
# Catch errors in the list of local maxima.
for crnt_ix in maxima_indices:
if crnt_ix < wing_length or crnt_ix > ia.size-wing_length-1:
error("local maxima too close to the boundary")
for ii in range(len(maxima_indices)-1):
crnt_ix = maxima_indices[ii]
next_ix = maxima_indices[ii+1]
if abs(crnt_ix-next_ix) < wing_length+1:
error("local maxima too close to one another")
for ix in range(len(maxima_indices)):
index = maxima_indices[ix]
# get the previous and next maximum
if ix == 0:
index_prev = 0
else:
index_prev = maxima_indices[ix-1]+1
if ix == len(maxima_indices)-1:
index_next = ia.size
else:
index_next = maxima_indices[ix+1]
left_slice = ia[index_prev:index]
left_slice = left_slice[::-1] # reverse left slice
right_slice = ia[index+1:index_next]
left_min_ix = _get_next_minimum(left_slice, wing_length)
right_min_ix = _get_next_minimum(right_slice, wing_length)
left_boundary = index - left_min_ix - 1
right_boundary = index + right_min_ix + 1
p=peak()
p.left=left_boundary
p.right=right_boundary
p.index=index
p.value=ia[index]
peak_list.append(p)
return peak_list
def _get_next_minimum(ia_slice, wing_length):
"""
@summary: Determines the first next minimum in 'ia_slice' based on 'wing_length',
'ia_slice' is processed as index increases until the first minimum
is found. If no minimum is found, the minimum is set to the first
point after 'wing_length' counting from the end of the slice
@param ia_slice: An numpy object. Ion chromatogram intensity
array slice
@type ia_slice: numpy.ndarray
@param wing_length: An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value
@type wing_length: IntType
@return: An integer. The index of the minimum value
@rtype: IntType
-----------------------------------------------------------------------
_get_next_minimum(ia_slice, wing_length)
Determines the first next minimum in 'ia_slice' based on 'wing_length'.
'ia_slice' is processed as index increases until the first minimum
is found. If no minimum is found, the minimum is set to the first
point after 'wing_length' counting from the end of the slice.
@param ia_slice An numpy object. Ion chromatogram intensity
array slice.
@param wing_length An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value.
@return An integer. The index of the minimum value.
"""
boundary_index = None
for index in range(ia_slice.size):
if _is_local_minimum(ia_slice, index, wing_length):
boundary_index = index
break
# if the minimum was not found, set the boundary to the last
# point in the slice. If the two peaks are exactly 'wing_length'
# points appart (the closest possible), the right boundary of
# first peak will be set to the left adjacent point of the apex
# of the second peak. And the left boundary of the second peak
# will be set to the right adjacent point of the apex of the
# first peak apex. If a peak is at the beginning/end of a signal,
# and no local minima are found when searching towards the boundary
# the boundary will be set to the first/last point of the signal.
if boundary_index == None:
boundary_index = ia_slice.size - 1
return boundary_index
##
# Returns True if the point is a local maximum or False otherwise.
# A point is a local maximum if the following are satisfied:
# 1. It is equal or greater than any of 'wing_length' points to
# its left and to its right.
# 2. It is greater than at least one of the points to its left,
# and one of the points to its right (out of total 'wing_length'
# points).
# 3. A point closer to the edge of 'ia' than 'wing_length' cannot
# be a local maximum.
#
# @param ia An numpy object, ion chromatogram intensity array.
# @param wing_length An integer. The width measured in number of
# points to the left and right of a point to be considered
# for extreme value.
# @return A boolean value.
def _is_local_maximum(ia, index, wing_length):
"""
@summary: Returns True if the point is a local maximum or False otherwise
A point is a local maximum if the following are satisfied:
1. It is equal or greater than any of 'wing_length' points to
its left and to its right
2. It is greater than at least one of the points to its left,
and one of the points to its right (out of total 'wing_length'
points)
3. A point closer to the edge of 'ia' than 'wing_length' cannot
be a local maximum
@param ia: An numpy object. Ion chromatogram intensity array
slice
@type ia: numpy.ndarray
@param index: An integer.Index
@type index: IntType
@param wing_length: An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value
@type wing_length: IntType
@return: If the point is a local maximum
@rtype: BooleanType
--------------------------------------------------------------------
_is_local_maximum(ia, index, wing_length)
Returns True if the point is a local maximum or False otherwise.
A point is a local maximum if the following are satisfied:
1. It is equal or greater than any of 'wing_length' points to
its left and to its right.
2. It is greater than at least one of the points to its left,
and one of the points to its right (out of total 'wing_length'
points).
3. A point closer to the edge of 'ia' than 'wing_length' cannot
be a local maximum.
@param ia An numpy object. Ion chromatogram intensity array
slice.
@param index An integer. Index at which to evaluate if local maximum.
@return A boolean value.
"""
if index > ia.size:
error("The index (%d) is out of bounds (max index: %d)." % \
(index, ia.size - 1))
intensity = ia[index]
valid=False
# Get the index for left/right boundary, fold the edges.
left_index = index - wing_length
if left_index > 0:
left_slice_max = numpy.max(ia[left_index:index])
if left_slice_max < intensity:
right_index = index + wing_length + 1
if right_index < ia.size:
right_slice_max = numpy.max(ia[index+1:right_index])
if right_slice_max < intensity:
valid=True
return valid
# # Test if local maximum. Points closer to the edges than
# # 'wing_length' cannot be a maximum.
# if left_slice.size < wing_length:
# left_wing_valid = False
# else:
# left_wing_valid = _is_greater_than(left_slice, intensity)
#
# if right_slice.size < wing_length:
# right_wing_valid = False
# else:
# right_wing_valid = _is_greater_than(right_slice, intensity)
#
# # Return value: True if both wings passed
# if left_wing_valid and right_wing_valid:
# return True
# else:
# return False
def _is_greater_than(ia, intensity):
"""
@summary: Returns True if 'intensity' is greater than or equal to every
element in 'ia' _and_ is greater than at least one element
@param ia: An numpy object. Ion chromatogram intensity array
slice
@type ia: numpy.ndarray
@param intensity: A number. Intensity to which slice is compared
@type intensity: FloatType
@return: If it is greater
@rtype: BooleanType
--------------------------------------------------------------
_is_greater_than(ia, intensity)
Returns True if 'intensity' is greater than or equal to every
element in 'ia' _and_ is greater than at least one element.
@param ia An numpy object. Ion chromatogram intensity array
slice.
@param intensity A number. Intensity to which slice is compared.
"""
if ia.size < 1:
error("empty array for comparison")
# numpy.where(ia <= intensity)[0] returns an array of indices
# refering to elements of 'ia' that are smaller or equal than
# 'intensity'
lteq = numpy.where(ia <= intensity)[0]
lt = numpy.where(ia < intensity)[0]
if lteq.size == ia.size and lt.size > 0:
return True
else:
return False
def _is_local_minimum(ia, index, wing_length):
"""
@summary: Returns True if the point is a local minimum or False otherwise
A point is a local maximum if the following are satisfied:
1. It is equal or smaller than any of 'wing_length' points to
its left and to its right
2. It is smllaer than at least one of the points to its left,
and one of the points to its right (out of total 'wing_length'
points)
3. A point closer to the edge of 'ia' than 'wing_length' cannot
be a local minimum
@param ia: An numpy object. Ion chromatogram intensity array
slice
@type ia: numpy.ndarray
@param index: An integer. Index at which to evaluate if local maximum
@type index: IntType
@param wing_length: An integer. The width measured in number of
points to the left and right of a point to be considered
for an extreme value
@type wing_length: IntType
@return: A boolean value
@rtype: BooleanType
--------------------------------------------------------------------
_is_local_minimum(ia, index, wing_length)
Returns True if the point is a local minimum or False otherwise.
A point is a local maximum if the following are satisfied:
1. It is equal or smaller than any of 'wing_length' points to
its left and to its right.
2. It is smllaer than at least one of the points to its left,
and one of the points to its right (out of total 'wing_length'
points).
3. A point closer to the edge of 'ia' than 'wing_length' cannot
be a local minimum.
@param ia An numpy object. Ion chromatogram intensity array
slice.
@param index An integer. Index at which to evaluate if local maximum.
@return A boolean value.
"""
if index > ia.size:
error("The index (%d) is out of bounds (max index: %d)." % \
(index, ia.size - 1))
intensity = ia[index]
# Get the index for left/right boundary, fold the edges.
left_index = index - wing_length
if left_index < 0: left_index = 0
left_slice = ia[left_index:index]
right_index = index + wing_length + 1
if right_index > ia.size: right_index = ia.size
right_slice = ia[index+1:right_index]
# Test if local maximum. Points closer to the edges than
# 'wing_length' cannot be a maximum.
if left_slice.size < wing_length:
left_wing_valid = False
else:
left_wing_valid = _is_less_than(left_slice, intensity)
if right_slice.size < wing_length:
right_wing_valid = False
else:
right_wing_valid = _is_less_than(right_slice, intensity)
# Return value: True if both wings passed
if left_wing_valid and right_wing_valid:
return True
else:
return False
def _is_less_than(ia, intensity):
"""
@summary: Returns True if 'intensity' is less than or equal to every
element in 'ia' _and_ is less than at least one element
@param ia: An numpy object. Ion chromatogram intensity array
slice
@type ia: numpy.ndarray
@param intensity: A number. Intensity to which slice is compared
@type intensity: FloatType
@return: a boolean value
@rtype: BooleanType
-----------------------------------------------------------------------
_is_less_than(ia, intensity)
Returns True if 'intensity' is less than or equal to every
element in 'ia' _and_ is less than at least one element.
@param ia An numpy object. Ion chromatogram intensity array
slice.
@param intensity A number. Intensity to which slice is compared.
"""
if ia.size < 1:
error("empty array for comparison")
# numpy.where(ia >= intensity)[0] returns an array of indices
# refering to elements of 'ia' that are greater or equal to
# 'intensity'
lteq = numpy.where(ia >= intensity)[0]
lt = numpy.where(ia > intensity)[0]
if lteq.size == ia.size and lt.size > 0:
return True
else:
return False
def _adjust_boundaries(ia, peak_list, angle_width=_DEFAULT_ANGLE_PTS, \
angle_threshold=_DEFAULT_ANGLE_THRESHOLD):
"""
@summary: Adjusts peak boundaries by fitting the straight line through the
edge points, and trimming point for which the line formes an angle
with the time axis which is less than the threshold
@param ia: An numpy object. Ion chromatogram intensity array
@type ia: numpy.ndarray
@param angle_width: An integer. The number of point for the best fit
@type angle_width: IntType
@param angle_threshold: A number. The angle threshold for peak boundary
points. A line of best fit drawn through an arbitrary number of
points at the peak boundaries must form an angle greater than
the angle threshold. If not, the boundary point will be trimmed
@type angle_threshold: FloatType
@return: No return value
@rtype: NoneType
----------------------------------------------------------------------
_adjust_boundaries(ia, peak_list, angle_width, angle_threshold)
Adjusts peak boundaries by fitting the straight line through the
edge points, and trimming point for which the line formes an angle
with the time axis which is less than the threshold.
@param ia An numpy object. Ion chromatogram intensity array.
@param peak_list A list of pre-peaks.
@param angle_width An integer. The number of point for the best fit.
@param angle_threshold A number. The angle threshold for peak boundary
points. A line of best fit drawn through an arbitrary number of
points at the peak boundaries must form an angle greater than
the angle threshold. If not, the boundary point will be trimmed.
@return No return value.
"""
for peak_item in peak_list:
left_index = peak_item.left
apex_index = peak_item.index
right_index = peak_item.right
apex = ia[apex_index]
left_slice = ia[left_index:apex_index]
right_slice = ia[apex_index+1:right_index+1]
# Reverse the numpy. _trim_boundary() expects the
# numpy from boundary to apex.
right_slice = right_slice[::-1]
left_slice = _trim_boundary(left_slice, apex, angle_width, \
angle_threshold)
right_slice = _trim_boundary(right_slice, apex, angle_width, \
angle_threshold)
peak_item.left = apex_index - left_slice.size
peak_item.right = apex_index + right_slice.size
def _trim_boundary(ia_slice, apex, angle_width, angle_threshold):
"""
@summary: Performs peak boundary trimming
@param ia_slice: An numpy object. Ion chromatogram intensity
array slice representing initial peak boundary
@type ia_slice: numpy.ndarray
@param apex: A number. Intensity at peak apex
@type apex: FloatType
@param angle_width: An integer. The number of point for the best fit
@type angle_width: IntType
@param angle_threshold: A number. The angle threshold for peak boundary
points. A line of best fit drawn through an arbitrary number of
points at the peak boundaries must form an angle greater than
the angle threshold. If not, the boundary point will be trimmed
@type angle_threshold: FloatType
@return: A numpy object. Trimmed peak boundary
@rtype: FloatType
-----------------------------------------------------------------------
_trim_boundary(ia_slice, apex, angle_width, angle_threshold)
Performs peak boundary trimming.
@param ia_slice An numpy object. Ion chromatogram intensity
array slice representing initial peak boundary.
@param apex A number. Intensity at peak apex .
@param angle_width An integer. The number of point for the best fit.
@param angle_threshold A number. The angle threshold for peak boundary
points. A line of best fit drawn through an arbitrary number of
points at the peak boundaries must form an angle greater than
the angle threshold. If not, the boundary point will be trimmed.
@return A numpy object. Trimmed peak boundary.
"""
if ia_slice.size <= angle_width: return ia_slice
index = 0
while index + angle_width < ia_slice.size:
x1 = index
x2 = index + angle_width
x_values = numpy.arange(x1, x2+1)
y_values = ia_slice[x_values]
gradient = _best_fit_line(y_values, x_values)[0]
rise = angle_width * gradient
rise = rise / apex
run = angle_width
angle = math.atan(rise/run)
angle = math.degrees(angle)
index = index + 1
if angle > angle_threshold: break
return ia_slice[index-1:]
def _best_fit_line(y_values, x_values=None):
"""
@summary: Calculates the coefficients to a linear line of best fit
@param y_values: A numpy object representing the y-values
@type y_values: numpy.ndarray
@param x_values: A numpy object representing the x-values
@type x_values: numpy.ndarray
@return: A tuple. Best fit coefficients where the first element is
the gradient, and the second is the y-intercept
@rtype: numpy.float64, numpy.float64
-----------------------------------------------------------------
_best_fit_line(y_values, x_values=None)
Calculates the coefficients to a linear line of best fit. The
@param y_values A numpy object representing the y-values
@param x_values A numpy object representing the x-values
@return A tuple. Best fit coefficients where the first element is
the gradient, and the second is the y-intercept.
"""
if x_values == None:
x_values = numpy.arange(y_values.size)
x_array = numpy.repeat([1], x_values.size)
x_array = numpy.array([x_values, x_array])
x_array = numpy.transpose(x_array)
y_array = copy.deepcopy(y_values)
alpha = numpy.dot(numpy.transpose(x_array), x_array)
beta = numpy.dot(numpy.transpose(x_array), y_array)
alphaM = numpy.asmatrix(alpha)
ialphaA = numpy.asarray(alphaM.I)
coefficients = numpy.dot(ialphaA, beta)
return coefficients[0], coefficients[1]
def peak_area(ia, peaks,da=1):
"""
@summary: Generates a list of peak objects from the list of pre-peaks,
and integrates peak areas
@param ia: A Numpy Array
@type ia: numpy.ndarray
@param peak_list: A list of pre-peaks
@type peak_list: ListType
@param da: area scale
@type peak_list: numpy.float64
@return: A list of peak objects
@rtype: ListType
-----------------------------------------------------------
generate_peaks(ic, peak_list)
Generates a list of peak objects from the list of pre-peaks,
and integrates peak areas.
@param ic An IonChromatogram object.
@param peak_list A list of pre-peaks.
@return A list of peak objects
"""
# ia = ia.get_intensity_array()
for peak in peaks:
pt_left = peak.left
pt_apex = peak.index
pt_right = peak.right
# peak areas is the sum of intensity from peak left
# to its right boundary (boundary points included)
area_slice = ia[pt_left:pt_right+1]
max_slice = numpy.max(ia[pt_left:pt_right+1])
peak.area = area_slice.sum()*da
return True
def peak_fwhm(ia, peaks,da=1):
"""
@summary: Generates a list of peak objects from the list of pre-peaks,
and integrates peak areas
@param ia: A Numpy Array
@type ia: numpy.ndarray
@param peak_list: A list of pre-peaks
@type peak_list: ListType
@param da: area scale
@type peak_list: numpy.float64
@return: A list of peak objects
@rtype: ListType
-----------------------------------------------------------
generate_peaks(ic, peak_list)
Generates a list of peak objects from the list of pre-peaks,
and integrates peak areas.
@param ic An IonChromatogram object.
@param peak_list A list of pre-peaks.
@return A list of peak objects
"""
# ia = ia.get_intensity_array()
widths=[]
for peak in peaks:
pt_left = peak.left
pt_apex = peak.index
pt_right = peak.right
# peak areas is the sum of intensity from peak left
# to its right boundary (boundary points included)
right_slice = ia[pt_apex:pt_right+1]
searcher=ia[pt_apex]
i=0
for i in xrange(len(right_slice)):
if right_slice[i] < float(searcher)/2.:
break
peak.fwhm=i*2*da
return True
def max_peaks(ia, peaks,da=1):
max_peaks=[]
max_value=0
for peak in peaks:
if peak.value > max_value: max_value=peak.value
for peak in peaks:
if peak.value==max_value: max_peaks.append(peak)
return max_peaks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment