#!/usr/bin/env python3.4
# coding: latin-1
# (c) Massachusetts Institute of Technology 2015-2018
# (c) Brian Teague 2018-2021
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
cytoflow.views.kde_2d
---------------------
"""
from traits.api import provides, Constant
import matplotlib.pyplot as plt
import numpy as np
from sklearn.neighbors import KernelDensity
from statsmodels.nonparametric.bandwidths import bw_scott, bw_silverman
import cytoflow.utility as util
from .i_view import IView
from .base_views import Base2DView
[docs]@provides(IView)
class Kde2DView(Base2DView):
"""
Plots a 2-d kernel-density estimate. Sort of like a smoothed histogram.
The density is visualized with a set of isolines.
Attributes
----------
Examples
--------
Make a little data set.
.. plot::
:context: close-figs
>>> import cytoflow as flow
>>> import_op = flow.ImportOp()
>>> import_op.tubes = [flow.Tube(file = "Plate01/RFP_Well_A3.fcs",
... conditions = {'Dox' : 10.0}),
... flow.Tube(file = "Plate01/CFP_Well_A4.fcs",
... conditions = {'Dox' : 1.0})]
>>> import_op.conditions = {'Dox' : 'float'}
>>> ex = import_op.apply()
Plot a density plot
.. plot::
:context: close-figs
>>> flow.Kde2DView(xchannel = 'V2-A',
... xscale = 'log',
... ychannel = 'Y2-A',
... yscale = 'log',
... huefacet = 'Dox').plot(ex)
"""
id = Constant('edu.mit.synbio.cytoflow.view.kde2d')
friend_id = Constant("2D Kernel Density Estimate")
[docs] def plot(self, experiment, **kwargs):
"""
Plot a faceted 2d kernel density estimate
Parameters
----------
shade : bool
Shade the interior of the isoplot? (default = `False`)
min_alpha, max_alpha : float
The minimum and maximum alpha blending values of the isolines,
between 0 (transparent) and 1 (opaque).
n_levels : int
How many isolines to draw? (default = 10)
bw : str or float
The bandwidth for the gaussian kernel, controls how lumpy or smooth the
kernel estimate is. Choices are:
- ``scott`` (the default) - ``1.059 * A * nobs ** (-1/5.)``, where ``A`` is ``min(std(X),IQR/1.34)``
- ``silverman`` - ``.9 * A * nobs ** (-1/5.)``, where ``A`` is ``min(std(X),IQR/1.34)``
If a float is given, it is the bandwidth. Note, this is in
scaled units, not data units.
gridsize : int
How many times to compute the kernel on each axis? (default: 100)
Notes
-----
Other ``kwargs`` are passed to `matplotlib.axes.Axes.contour <https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.contour.html>`_
"""
super().plot(experiment, **kwargs)
def _grid_plot(self, experiment, grid, **kwargs):
kwargs.setdefault('shade', False)
kwargs.setdefault('min_alpha', 0.2)
kwargs.setdefault('max_alpha', 0.9)
kwargs.setdefault('n_levels', 10)
lim = kwargs.pop('lim')
xlim = lim[self.xchannel]
ylim = lim[self.ychannel]
scale = kwargs.pop('scale')
xscale = scale[self.xchannel]
yscale = scale[self.ychannel]
legend_data = {}
grid.map(_bivariate_kdeplot,
self.xchannel,
self.ychannel,
xscale = xscale,
yscale = yscale,
legend_data = legend_data,
**kwargs)
return dict(xlim = xlim,
xscale = xscale,
ylim = ylim,
yscale = yscale,
legend_data = legend_data)
# yoinked from seaborn/distributions.py, with modifications for scaling.
def _bivariate_kdeplot(x, y, xscale=None, yscale=None, shade=False,
bw="scott", gridsize=50, cut=3, clip=None, legend=True,
legend_data = None, **kwargs):
ax = plt.gca()
label = kwargs.pop('label', None)
# Determine the clipping
clip = [(-np.inf, np.inf), (-np.inf, np.inf)]
x = xscale(x)
y = yscale(y)
x_nan = np.isnan(x)
y_nan = np.isnan(y)
x = x[~(x_nan | y_nan)]
y = y[~(x_nan | y_nan)]
if bw == 'scott':
bw_x = bw_scott(x)
bw_y = bw_scott(y)
bw = (bw_x + bw_y) / 2
elif bw == 'silverman':
bw_x = bw_silverman(x)
bw_y = bw_silverman(y)
bw = (bw_x + bw_y) / 2
elif isinstance(bw, float):
bw_x = bw_y = bw
else:
raise util.CytoflowViewError(None,
"Bandwith must be 'scott', 'silverman' or a float")
kde = KernelDensity(bandwidth = bw, kernel = 'gaussian').fit(np.column_stack((x, y)))
x_support = _kde_support(x, bw_x, gridsize, cut, clip[0])
y_support = _kde_support(y, bw_y, gridsize, cut, clip[1])
xx, yy = np.meshgrid(x_support, y_support)
z = kde.score_samples(np.column_stack((xx.ravel(), yy.ravel())))
z = z.reshape(xx.shape)
z = np.exp(z)
n_levels = kwargs.pop("n_levels", 10)
color = kwargs.pop("color")
kwargs['colors'] = (color, )
min_alpha = kwargs.pop("min_alpha", 0.2)
if shade:
min_alpha = 0
max_alpha = kwargs.pop("max_alpha", 0.9)
x_support = xscale.inverse(x_support)
y_support = yscale.inverse(y_support)
xx, yy = np.meshgrid(x_support, y_support)
contour_func = ax.contourf if shade else ax.contour
try:
cset = contour_func(xx, yy, z, n_levels, **kwargs)
except ValueError as e:
raise util.CytoflowViewError(None,
"Something went wrong in {}, bandwidth = {}. "
.format(contour_func.__name__, bw)) from e
num_collections = len(cset.collections)
alpha = np.linspace(min_alpha, max_alpha, num = num_collections)
for el in range(num_collections):
cset.collections[el].set_alpha(alpha[el])
# Label the axes
if hasattr(x, "name") and legend:
ax.set_xlabel(x.name)
if hasattr(y, "name") and legend:
ax.set_ylabel(y.name)
if label is not None:
ax.set_title(label)
# Add legend data
if 'label' in kwargs:
legend_data[kwargs['label']] = plt.Rectangle((0, 0), 1, 1, fc = color)
return ax
def _kde_support(data, bw, gridsize, cut, clip):
"""Establish support for a kernel density estimate."""
support_min = max(data.min() - bw * cut, clip[0])
support_max = min(data.max() + bw * cut, clip[1])
return np.linspace(support_min, support_max, gridsize)
util.expand_class_attributes(Kde2DView)
util.expand_method_parameters(Kde2DView, Kde2DView.plot)
if __name__ == '__main__':
import cytoflow as flow
tube1 = flow.Tube(file = '../../cytoflow/tests/data/Plate01/RFP_Well_A3.fcs',
conditions = {"Dox" : 10.0})
tube2 = flow.Tube(file = '../../cytoflow/tests/data/Plate01/CFP_Well_A4.fcs',
conditions = {"Dox" : 1.0})
ex = flow.ImportOp(conditions = {"Dox" : "float"}, tubes = [tube1, tube2])
thresh = flow.ThresholdOp()
thresh.name = "Y2-A+"
thresh.channel = 'Y2-A'
thresh.threshold = 200.0
ex2 = thresh.apply(ex)
scatter = flow.ScatterplotView()
scatter.name = "Scatter"
scatter.xchannel = "FSC-A"
scatter.ychannel = "SSC-A"
scatter.xscale = "logicle"
scatter.yscale = "logicle"
scatter.huefacet = 'Dox'
plt.ioff()
scatter.plot(ex2)
plt.show()