Skip to content

Commit

Permalink
Merge branch 'release-1.3.0'
Browse files Browse the repository at this point in the history
  • Loading branch information
@JS3xton
JS3xton committed Jan 25, 2021
2 parents 5e187b0 + 2e01ed0 commit 5a6c1e3
Show file tree
Hide file tree
Showing 98 changed files with 4,255 additions and 832 deletions.
2 changes: 1 addition & 1 deletion FlowCal/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
# Versions should comply with PEP440. For a discussion on single-sourcing
# the version across setup.py and the project code, see
# https://packaging.python.org/en/latest/single_source_version.html
__version__ = '1.2.2'
__version__ = '1.3.0'

from . import io
from . import excel_ui
Expand Down
458 changes: 191 additions & 267 deletions FlowCal/excel_ui.py
View file

Large diffs are not rendered by default.

237 changes: 142 additions & 95 deletions FlowCal/gate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,32 @@
import skimage.measure
import collections

###
# Gate Classes
###

# Output namedtuples returned by gate functions
StartEndGateOutput = collections.namedtuple(
typename='StartEndGateOutput',
field_names=('gated_data',
'mask'))
HighLowGateOutput = collections.namedtuple(
typename='HighLowGateOutput',
field_names=('gated_data',
'mask'))
EllipseGateOutput = collections.namedtuple(
typename='EllipseGateOutput',
field_names=('gated_data',
'mask',
'contour'))
Density2dGateOutput = collections.namedtuple(
typename='Density2dGateOutput',
field_names=('gated_data',
'mask',
'contour',
'bin_edges',
'bin_mask'))

###
# Gate Functions
###
Expand Down Expand Up @@ -78,9 +104,6 @@ def start_end(data, num_start=250, num_end=100, full_output=False):
gated_data = data[mask]

if full_output:
StartEndGateOutput = collections.namedtuple(
'StartEndGateOutput',
['gated_data', 'mask'])
return StartEndGateOutput(gated_data=gated_data, mask=mask)
else:
return gated_data
Expand Down Expand Up @@ -143,9 +166,6 @@ def high_low(data, channels=None, high=None, low=None, full_output=False):
gated_data = data[mask]

if full_output:
HighLowGateOutput = collections.namedtuple(
'HighLowGateOutput',
['gated_data', 'mask'])
return HighLowGateOutput(gated_data=gated_data, mask=mask)
else:
return gated_data
Expand Down Expand Up @@ -232,9 +252,6 @@ def ellipse(data, channels,
cntr = [ci]

# Build output namedtuple
EllipseGateOutput = collections.namedtuple(
'EllipseGateOutput',
['gated_data', 'mask', 'contour'])
return EllipseGateOutput(
gated_data=data_gated, mask=mask, contour=cntr)
else:
Expand All @@ -247,6 +264,7 @@ def density2d(data,
xscale='logicle',
yscale='logicle',
sigma=10.0,
bin_mask=None,
full_output=False):
"""
Gate that preserves events in the region with highest density.
Expand Down Expand Up @@ -298,6 +316,10 @@ def density2d(data,
Standard deviation for Gaussian kernel used by
`scipy.ndimage.filters.gaussian_filter` to smooth 2D histogram
into a density.
bin_mask : 2D numpy array of bool, optional
A 2D mask array that selects the 2D histogram bins permitted by the
gate. Corresponding bin edges should be specified via `bins`. If
`bin_mask` is specified, `gate_fraction` and `sigma` are ignored.
full_output : bool, optional
Flag specifying to return additional outputs. If true, the outputs
are given as a namedtuple.
Expand All @@ -310,8 +332,14 @@ def density2d(data,
Boolean gate mask used to gate data such that ``gated_data =
data[mask]``.
contour : list of 2D numpy arrays, only if ``full_output==True``
List of 2D numpy array(s) of x-y coordinates tracing out
the edge of the gated region.
List of 2D numpy array(s) of x-y coordinates tracing out the edge of
the gated region. If `bin_mask` is specified, `contour` is None.
bin_edges : 2-tuple of numpy arrays, only if ``full_output==True``
X-axis and y-axis bin edges used by the np.histogram2d() command that
bins events (bin_edges=(x_edges,y_edges)).
bin_mask : 2D numpy array of bool, only if ``full_output==True``
A 2D mask array that selects the 2D histogram bins permitted by the
gate.
Raises
------
Expand Down Expand Up @@ -355,22 +383,12 @@ def density2d(data,
if data_ch.ndim == 1:
data_ch = data_ch.reshape((-1,1))

# Check gating fraction
if gate_fraction < 0 or gate_fraction > 1:
raise ValueError('gate fraction should be between 0 and 1, inclusive')

# Check dimensions
if data_ch.ndim < 2:
raise ValueError('data should have at least 2 dimensions')
if data_ch.shape[0] <= 1:
raise ValueError('data should have more than one event')

# Build output namedtuple if necessary
if full_output:
Density2dGateOutput = collections.namedtuple(
'Density2dGateOutput',
['gated_data', 'mask', 'contour'])

# If ``data_ch.hist_bins()`` exists, obtain bin edges from it if
# necessary.
if hasattr(data_ch, 'hist_bins') and \
Expand Down Expand Up @@ -405,6 +423,8 @@ def density2d(data,

# Make 2D histogram
H,xe,ye = np.histogram2d(data_ch[:,0], data_ch[:,1], bins=bins)
xe = np.array(xe, dtype=float)
ye = np.array(ye, dtype=float)

# Map each event to its histogram bin by sorting events into a 2D array of
# lists which mimics the histogram.
Expand Down Expand Up @@ -447,91 +467,118 @@ def density2d(data,
# Create a 2D array of lists mimicking the histogram to accumulate events
# associated with each bin.
filler = np.frompyfunc(lambda x: list(), 1, 1)
H_events = np.empty_like(H, dtype=np.object)
H_events = np.empty_like(H, dtype=object)
filler(H_events, H_events)

for event_idx, x_bin_idx, y_bin_idx in \
zip(event_indices, x_bin_indices, y_bin_indices):
H_events[x_bin_idx, y_bin_idx].append(event_idx)

# Determine number of events to keep. Only consider events which have not
# been thrown out as outliers.
n = int(np.ceil(gate_fraction*float(len(event_indices))))

# n = 0 edge case (e.g. if gate_fraction = 0.0); incorrectly handled below
if n == 0:
mask = np.zeros(shape=data_ch.shape[0], dtype=bool)
gated_data = data[mask]
# Create bin mask if necessary
contours = None
if bin_mask is None:
# Check gating fraction
if gate_fraction < 0 or gate_fraction > 1:
msg = "gate fraction should be between 0 and 1, inclusive"
raise ValueError(msg)

# Determine number of events to keep. Only consider events which have
# not been thrown out as outliers.
n = int(np.ceil(gate_fraction*float(len(event_indices))))

# n = 0 edge case (e.g. if gate_fraction = 0.0); incorrectly handled
# below
if n == 0:
mask = np.zeros(shape=data_ch.shape[0], dtype=bool)
gated_data = data[mask]
if full_output:
return Density2dGateOutput(
gated_data=gated_data,
mask=mask,
contour=[],
bin_edges=(xe,ye),
bin_mask=np.zeros_like(H, dtype=bool))
else:
return gated_data

# Smooth 2D histogram
sH = scipy.ndimage.filters.gaussian_filter(
H,
sigma=sigma,
order=0,
mode='constant',
cval=0.0,
truncate=6.0)

# Normalize smoothed histogram to make it a valid probability mass
# function
D = sH / np.sum(sH)

# Sort bins by density
vD = D.ravel(order='C')
vH = H.ravel(order='C')
sidx = np.argsort(vD)[::-1]
svH = vH[sidx] # linearized counts array sorted by density

# Find minimum number of accepted bins needed to reach specified
# number of events
csvH = np.cumsum(svH)
Nidx = np.nonzero(csvH >= n)[0][0] # we want to include this index

# Get indices of accepted histogram bins
accepted_bin_indices = sidx[:(Nidx+1)]

# Convert accepted bin indices to bin mask array
bin_mask = np.zeros_like(H, dtype=bool)
v_bin_mask = bin_mask.ravel(order='C')
v_bin_mask[accepted_bin_indices] = True
bin_mask = v_bin_mask.reshape(H.shape, order='C')

# Determine contours if necessary
if full_output:
return Density2dGateOutput(
gated_data=gated_data, mask=mask, contour=[])
else:
return gated_data

# Smooth 2D histogram
sH = scipy.ndimage.filters.gaussian_filter(
H,
sigma=sigma,
order=0,
mode='constant',
cval=0.0,
truncate=6.0)

# Normalize smoothed histogram to make it a valid probability mass function
D = sH / np.sum(sH)

# Sort bins by density
vD = D.ravel()
vH = H.ravel()
sidx = np.argsort(vD)[::-1]
svH = vH[sidx] # linearized counts array sorted by density

# Find minimum number of accepted bins needed to reach specified number
# of events
csvH = np.cumsum(svH)
Nidx = np.nonzero(csvH >= n)[0][0] # we want to include this index

# Get indices of events to keep
vH_events = H_events.ravel()
accepted_indices = vH_events[sidx[:(Nidx+1)]]
accepted_indices = np.array([item # flatten list of lists
for sublist in accepted_indices
for item in sublist])
accepted_indices = np.sort(accepted_indices)

# Convert list of accepted indices to boolean mask array
# Use scikit-image to find the contour of the gated region
#
# To find the contour of the gated region, values in the 2D
# probability mass function ``D`` are used to trace contours at
# the level of the probability associated with the last accepted
# bin, ``vD[sidx[Nidx]]``.

# find_contours() specifies contours as collections of row and
# column indices into the density matrix. The row or column index
# may be interpolated (i.e. non-integer) for greater precision.
contours_ij = skimage.measure.find_contours(D, vD[sidx[Nidx]])

# Map contours from indices into density matrix to histogram x and
# y coordinate spaces (assume values in the density matrix are
# associated with histogram bin centers).
xc = (xe[:-1] + xe[1:]) / 2.0 # x-axis bin centers
yc = (ye[:-1] + ye[1:]) / 2.0 # y-axis bin centers

contours = [np.array([np.interp(contour_ij[:,0],
np.arange(len(xc)),
xc),
np.interp(contour_ij[:,1],
np.arange(len(yc)),
yc)]).T
for contour_ij in contours_ij]

accepted_data_indices = H_events[bin_mask]
accepted_data_indices = np.array([item # flatten list of lists
for sublist in accepted_data_indices
for item in sublist],
dtype=int)

# Convert list of accepted data indices to boolean mask array
mask = np.zeros(shape=data.shape[0], dtype=bool)
mask[accepted_indices] = True
mask[accepted_data_indices] = True

gated_data = data[mask]

if full_output:
# Use scikit-image to find the contour of the gated region
#
# To find the contour of the gated region, values in the 2D probability
# mass function ``D`` are used to trace contours at the level of the
# probability associated with the last accepted bin, ``vD[sidx[Nidx]]``.

# find_contours() specifies contours as collections of row and column
# indices into the density matrix. The row or column index may be
# interpolated (i.e. non-integer) for greater precision.
contours_ij = skimage.measure.find_contours(D, vD[sidx[Nidx]])

# Map contours from indices into density matrix to histogram x and y
# coordinate spaces (assume values in the density matrix are associated
# with histogram bin centers).
xc = (xe[:-1] + xe[1:]) / 2.0 # x-axis bin centers
yc = (ye[:-1] + ye[1:]) / 2.0 # y-axis bin centers

contours = [np.array([np.interp(contour_ij[:,0],
np.arange(len(xc)),
xc),
np.interp(contour_ij[:,1],
np.arange(len(yc)),
yc)]).T
for contour_ij in contours_ij]

return Density2dGateOutput(
gated_data=gated_data, mask=mask, contour=contours)
return Density2dGateOutput(gated_data=gated_data,
mask=mask,
contour=contours,
bin_edges=(xe,ye),
bin_mask=bin_mask)
else:
return gated_data
Loading

0 comments on commit 5a6c1e3

Please sign in to comment.