Plotting

Routines for plotting comparisons between model and data.

These can serve as inspiration for custom routines for one's own purposes. Note that all the plotting is done with pylab. To see additional pylab methods: "import pylab; help(pylab)". Pylab's many functions are documented at http://matplotlib.sourceforge.net/contents.html.

plot_1d_comp_Poisson(model, data, fig_num=None, residual='Anscombe', plot_masked=False, show=True)

Poisson comparison between 1D model and data.

Parameters:

Name	Type	Description	Default
model	function	1-dimensional model SFS.	required
data	Spectrum	1-dimensional data SFS.	required
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
residual	str	'Anscombe' for Anscombe residuals, which are more normally distributed for Poisson sampling. 'linear' for the linear residuals, which can be less biased.	'Anscombe'
plot_masked	bool	Additionally plots (in open circles) results for points in the model or data that were masked.	False
show	bool	If True, execute pylab.show command to make sure plot displays.	True

Source code in dadi/Plotting.py

def plot_1d_comp_Poisson(model, data, fig_num=None, residual='Anscombe',
                         plot_masked=False, show=True):
    """
    Poisson comparison between 1D model and data.

    Args:
        model (function): 1-dimensional model SFS.
        data (Spectrum): 1-dimensional data SFS.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        residual (str): 'Anscombe' for Anscombe residuals, which are more normally
            distributed for Poisson sampling. 'linear' for the linear residuals,
            which can be less biased.
        plot_masked (bool): Additionally plots (in open circles) results for points in
            the model or data that were masked.
        show (bool): If True, execute pylab.show command to make sure plot displays.
    """
    if fig_num is None:
        f = pylab.gcf()
    else:
        f = pylab.figure(fig_num, figsize=(7,7))
    pylab.clf()

    if data.folded and not model.folded:
        model = model.fold()

    masked_model, masked_data = Numerics.intersect_masks(model, data)

    ax = pylab.subplot(2,1,1)
    pylab.semilogy(masked_data, '-ob', label='data')
    pylab.semilogy(masked_model, '-or', label='model')

    if plot_masked:
        pylab.semilogy(masked_data.data, '--ob', mfc='w', zorder=-100)
        pylab.semilogy(masked_model.data, '--or', mfc='w', zorder=-100)

    ax.legend(loc='upper right')

    pylab.subplot(2,1,2, sharex = ax)
    if residual == 'Anscombe':
        resid = Inference.Anscombe_Poisson_residual(masked_model, masked_data)
    elif residual == 'linear':
        resid = Inference.linear_Poisson_residual(masked_model, masked_data)
    else:
        raise ValueError("Unknown class of residual '%s'." % residual)
    pylab.plot(resid, '-og')
    if plot_masked:
        pylab.plot(resid.data, '--og', mfc='w', zorder=-100)

    ax.set_xlim(0, data.shape[0]-1)
    if show:
        pylab.show()

plot_1d_comp_multinom(model, data, fig_num=None, residual='Anscombe', plot_masked=False, show=True)

Multinomial comparison between 1D model and data.

Parameters:

Name	Type	Description	Default
model	function	1-dimensional model SFS.	required
data	Spectrum	1-dimensional data SFS.	required
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
residual	str	'Anscombe' for Anscombe residuals, which are more normally distributed for Poisson sampling. 'linear' for the linear residuals, which can be less biased.	'Anscombe'
plot_masked	bool	Additionally plots (in open circles) results for points in the model or data that were masked.	False
show	bool	If True, execute pylab.show command to make sure plot displays.	True

Note

This comparison is multinomial in that it rescales the model to optimally fit the data.

Source code in dadi/Plotting.py

def plot_1d_comp_multinom(model, data, fig_num=None, residual='Anscombe',
                          plot_masked=False, show=True):
    """
    Multinomial comparison between 1D model and data.

    Args:
        model (function): 1-dimensional model SFS.
        data (Spectrum): 1-dimensional data SFS.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        residual (str): 'Anscombe' for Anscombe residuals, which are more normally
            distributed for Poisson sampling. 'linear' for the linear residuals,
            which can be less biased.
        plot_masked (bool): Additionally plots (in open circles) results for points in
            the model or data that were masked.
        show (bool): If True, execute pylab.show command to make sure plot displays.

    Note:
        This comparison is multinomial in that it rescales the model to
        optimally fit the data.
    """
    model = Inference.optimally_scaled_sfs(model, data)

    plot_1d_comp_Poisson(model, data, fig_num, residual,
                         plot_masked, show=show)

plot_1d_fs(fs, fig_num=None, show=True)

Plot a 1-dimensional frequency spectrum.

Parameters:

Name	Type	Description	Default
fs	Spectrum	1-dimensional Spectrum.	required
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
show	bool	If True, execute pylab.show command to make sure plot displays.	True

Note

All the plotting is done with pylab. To see additional pylab methods: "import pylab; help(pylab)". Pylab's many functions are documented at http://matplotlib.sourceforge.net/contents.html.

Source code in dadi/Plotting.py

def plot_1d_fs(fs, fig_num=None, show=True):
    """
    Plot a 1-dimensional frequency spectrum.

    Args:
        fs (Spectrum): 1-dimensional Spectrum.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        show (bool): If True, execute pylab.show command to make sure plot displays.

    Note:
        All the plotting is done with pylab. To see additional pylab methods:
        "import pylab; help(pylab)". Pylab's many functions are documented at
        http://matplotlib.sourceforge.net/contents.html.
    """

    if fig_num is None:
        fig = pylab.gcf()
    else:
        fig = pylab.figure(fig_num, figsize=(7,7))
    fig.clear()

    ax = fig.add_subplot(1,1,1)
    ax.semilogy(fs, '-ob')

    ax.set_xlim(0, fs.sample_sizes[0])
    if show:
        fig.show()

plot_2d_comp_Poisson(model, data, vmin=None, vmax=None, resid_range=None, fig_num=None, pop_ids=None, residual='Anscombe', adjust=True, show=True)

Poisson comparison between 2D model and data.

Parameters:

Name	Type	Description	Default
model	Spectrum	2-dimensional model SFS.	required
data	Spectrum	2-dimensional data SFS.	required
vmin	float	Minimum values plotted for sfs.	None
vmax	float	Maximum values plotted for sfs.	None
resid_range	float	Residual plot saturates at +- resid_range.	None
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
residual	str	'Anscombe' for Anscombe residuals, which are more normally distributed for Poisson sampling. 'linear' for the linear residuals, which can be less biased.	'Anscombe'
adjust	bool	Should method use automatic 'subplots_adjust'? For advanced manipulation of plots, it may be useful to make this False.	True

Source code in dadi/Plotting.py

def plot_2d_comp_Poisson(model, data, vmin=None, vmax=None,
                         resid_range=None, fig_num=None,
                         pop_ids=None, residual='Anscombe',
                         adjust=True, show=True):
    """
    Poisson comparison between 2D model and data.

    Args:
        model (Spectrum): 2-dimensional model SFS.
        data (Spectrum): 2-dimensional data SFS.
        vmin (float): Minimum values plotted for sfs.
        vmax (float): Maximum values plotted for sfs.
        resid_range (float): Residual plot saturates at +- resid_range.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        residual (str): 'Anscombe' for Anscombe residuals, which are more normally
            distributed for Poisson sampling. 'linear' for the linear residuals,
            which can be less biased.
        adjust (bool): Should method use automatic 'subplots_adjust'? For advanced
            manipulation of plots, it may be useful to make this False.
    """
    if data.folded and not model.folded:
        model = model.fold()

    masked_model, masked_data = Numerics.intersect_masks(model, data)

    if fig_num is None:
        f = pylab.gcf()
    else:
        f = pylab.figure(fig_num, figsize=(7,7))

    pylab.clf()
    if adjust:
        pylab.subplots_adjust(bottom=0.07, left=0.07, top=0.94, right=0.95, 
                              hspace=0.26, wspace=0.26)

    max_toplot = max(masked_model.max(), masked_data.max())
    min_toplot = min(masked_model.min(), masked_data.min())
    if vmax is None:
        vmax = max_toplot
    if vmin is None:
        vmin = min_toplot
    extend = _extend_mapping[vmin <= min_toplot, vmax >= max_toplot]

    if pop_ids is not None:
        data_pop_ids = model_pop_ids = resid_pop_ids = pop_ids
        if len(pop_ids) != 2:
            raise ValueError('pop_ids must be of length 2.')
    else:
        data_pop_ids = masked_data.pop_ids
        model_pop_ids = masked_model.pop_ids
        if masked_model.pop_ids is None:
            model_pop_ids = data_pop_ids

        if model_pop_ids == data_pop_ids (list[str]):
           resid_pop_ids = model_pop_ids
        else:
            resid_pop_ids = None

    ax = pylab.subplot(2,2,1)
    plot_single_2d_sfs(masked_data, vmin=vmin, vmax=vmax,
                       pop_ids=data_pop_ids, colorbar=False, show=False)
    ax.set_title('data')

    ax2 = pylab.subplot(2,2,2, sharex=ax, sharey=ax)
    plot_single_2d_sfs(masked_model, vmin=vmin, vmax=vmax,
                       pop_ids=model_pop_ids, extend=extend, show=False)
    ax2.set_title('model')

    if residual == 'Anscombe':
        resid = Inference.Anscombe_Poisson_residual(masked_model, masked_data,
                                              mask=vmin)
    elif residual == 'linear':
        resid = Inference.linear_Poisson_residual(masked_model, masked_data,
                                            mask=vmin)
    else:
        raise ValueError("Unknown class of residual '%s'." % residual)

    if resid_range is None:
        resid_range = max((abs(resid.max()), abs(resid.min())))
    resid_extend = _extend_mapping[-resid_range <= resid.min(), 
                                   resid_range >= resid.max()]

    ax3 = pylab.subplot(2,2,3, sharex=ax, sharey=ax)
    plot_2d_resid(resid, resid_range, pop_ids=resid_pop_ids,
                  extend=resid_extend)
    ax3.set_title('residuals')

    ax = pylab.subplot(2,2,4)
    flatresid = numpy.compress(numpy.logical_not(resid.mask.ravel()), 
                               resid.ravel())
    ax.hist(flatresid, bins=20, density=True)
    ax.set_title('residuals')
    ax.set_yticks([])
    if show:
        pylab.show()

plot_2d_comp_multinom(model, data, vmin=None, vmax=None, resid_range=None, fig_num=None, pop_ids=None, residual='Anscombe', adjust=True, show=True)

Multinomial comparison between 2D model and data.

Parameters:

Name	Type	Description	Default
model	Spectrum	2-dimensional model SFS.	required
data	Spectrum	2-dimensional data SFS.	required
vmin	float	Minimum values plotted for sfs.	None
vmax	float	Maximum values plotted for sfs.	None
resid_range	float	Residual plot saturates at +- resid_range.	None
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
residual	str	'Anscombe' for Anscombe residuals, which are more normally distributed for Poisson sampling. 'linear' for the linear residuals, which can be less biased.	'Anscombe'
adjust	bool	Should method use automatic 'subplots_adjust'? For advanced manipulation of plots, it may be useful to make this False.	True
show	bool	Display the figure? False is useful for saving many comparisons in a loop.	True

Note

This comparison is multinomial in that it rescales the model to optimally fit the data.

Source code in dadi/Plotting.py

def plot_2d_comp_multinom(model, data, vmin=None, vmax=None,
                          resid_range=None, fig_num=None,
                          pop_ids=None, residual='Anscombe',
                          adjust=True,show=True):
    """
    Multinomial comparison between 2D model and data.

    Args:
        model (Spectrum): 2-dimensional model SFS.
        data (Spectrum): 2-dimensional data SFS.
        vmin (float): Minimum values plotted for sfs.
        vmax (float): Maximum values plotted for sfs.
        resid_range (float): Residual plot saturates at +- resid_range.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        residual (str): 'Anscombe' for Anscombe residuals, which are more normally
            distributed for Poisson sampling. 'linear' for the linear residuals,
            which can be less biased.
        adjust (bool): Should method use automatic 'subplots_adjust'? For advanced
            manipulation of plots, it may be useful to make this False.
        show (bool): Display the figure? False is useful for saving many comparisons
            in a loop.

    Note:
        This comparison is multinomial in that it rescales the model to
        optimally fit the data.
    """
    model = Inference.optimally_scaled_sfs(model, data)

    plot_2d_comp_Poisson(model, data, vmin=vmin, vmax=vmax,
                         resid_range=resid_range, fig_num=fig_num,
                         pop_ids=pop_ids, residual=residual,
                         adjust=adjust,show=show)

plot_2d_meta_resid(s_resid, ns_resid, resid_range=None, fig_num=None, pop_ids=None, adjust=True, show=True)

Comparison between 2D nonsynonymous residual and 2D synonymous residual.

Parameters:

Name	Type	Description	Default
s_resid	array - like	residual SFS from synonymous data.	required
ns_resid	array - like	residual SFS from nonsynonymous data.	required
resid_range	float	Residual plot saturates at +- resid_range. This range applies to both the residual SFS's supplied as well as the meta-residual plot.	None
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
adjust	bool	Should method use automatic 'subplots_adjust'? For advanced manipulation of plots, it may be useful to make this False.	True
show	bool	Display the plot? False can be useful when plotting many in a loop.	True

Source code in dadi/Plotting.py

def plot_2d_meta_resid(s_resid,ns_resid,resid_range=None,
                       fig_num=None, pop_ids=None, 
                       adjust=True, show=True):

    """
    Comparison between 2D nonsynonymous residual and 2D synonymous residual.

    Args:
        s_resid (array-like): residual SFS from synonymous data.
        ns_resid (array-like): residual SFS from nonsynonymous data.
        resid_range (float): Residual plot saturates at +- resid_range. This range applies to both
            the residual SFS's supplied as well as the meta-residual plot.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        adjust (bool): Should method use automatic 'subplots_adjust'? For advanced
            manipulation of plots, it may be useful to make this False.
        show (bool): Display the plot? False can be useful when plotting many in a loop.
    """

    if ns_resid.folded and not s_resid.folded:
        s_resid = s_resid.fold()

    masked_s, masked_ns = Numerics.intersect_masks(s_resid,ns_resid)

    if fig_num is None:
        f = pylab.gcf()
    else:
        f = pylab.figure(fig_num, figsize=(7,7))

    pylab.clf()
    if adjust:
        pylab.subplots_adjust(bottom=0.07, left=0.07, top=0.94, right=0.95, 
                              hspace=0.26, wspace=0.26) 

    max_toplot = max(masked_s.max(), masked_ns.max())
    min_toplot = min(masked_s.min(), masked_ns.min())

    if pop_ids is not None:
        ns_pop_ids = s_pop_ids = resid_pop_ids = pop_ids
        if len(pop_ids) != 2:
            raise ValueError('pop_ids must be of length 2.')
    else:
        ns_pop_ids = masked_ns.pop_ids
        s_pop_ids = masked_s.pop_ids
        if masked_s.pop_ids is None:
            s_pop_ids = ns_pop_ids

        if s_pop_ids == ns_pop_ids (list[str]):
           resid_pop_ids = s_pop_ids
        else:
            resid_pop_ids = None

    if resid_range is None:
        resid_range = max((abs(masked_s.max()), abs(masked_s.min())))
    resid_extend = _extend_mapping[-resid_range <= masked_s.min(), 
                                   resid_range >= masked_s.max()]

    ax = pylab.subplot(2,2,1)
    plot_2d_resid(masked_s, resid_range=resid_range, pop_ids=resid_pop_ids,
                  extend=resid_extend)
    ax.set_title('Synonymous Residuals')

    if resid_range is None:
        resid_range = max((abs(masked_ns.max()), abs(masked_ns.min())))
    resid_extend = _extend_mapping[-resid_range <= masked_ns.min(), 
                                   resid_range >= masked_ns.max()]

    ax2 = pylab.subplot(2,2,2, sharex=ax, sharey=ax)
    plot_2d_resid(masked_ns, resid_range=resid_range, pop_ids=resid_pop_ids,
                  extend=resid_extend)
    ax2.set_title('Nonsynonymous Residuals')

    resid = masked_s-masked_ns

    if resid_range is None:
        resid_range = max((abs(resid.max()), abs(resid.min())))
    resid_extend = _extend_mapping[-resid_range <= resid.min(), 
                                   resid_range >= resid.max()]

    ax3 = pylab.subplot(2,2,3, sharex=ax, sharey=ax)
    plot_2d_resid(resid, resid_range, pop_ids=resid_pop_ids,
                  extend=resid_extend,cmap=pylab.cm.PuOr_r)
    ax3.set_title('Meta-residuals')

    ax = pylab.subplot(2,2,4)
    flatresid = numpy.compress(numpy.logical_not(resid.mask.ravel()), 
                               resid.ravel())
    ax.hist(flatresid, bins=20, density=True,color='purple')

    resid.data[resid.mask==True]=0
    sum_squares=numpy.sum(resid.data**2)
    ax.set_title(r'$res^2$ = '+'{0:.3f}'.format(sum_squares))
    ax.set_yticks([])

    if show:
        pylab.show()

plot_2d_resid(resid, resid_range=None, ax=None, pop_ids=None, extend='neither', colorbar=True, cmap=pylab.cm.RdBu_r)

Linear heatmap of 2D residual array.

Parameters:

Name	Type	Description	Default
resid	array - like	Residual array to plot.	required
resid_range	float	Values > resid_range or < -resid_range saturate the color spectrum.	None
ax	int	Axes object to plot into. If None, the result of pylab.gca() is used.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
extend	str	Whether the colorbar should have 'extension' arrows. See help(pylab.colorbar) for more details.	'neither'
colorbar	bool	Should we plot a colorbar?	True
cmap	cm	Pylab colormap to use for plotting.	RdBu_r

Returns:

Name	Type	Description
cb	figure element	The created colorbar.

Source code in dadi/Plotting.py

def plot_2d_resid(resid, resid_range=None, ax=None, pop_ids=None,
                  extend='neither', colorbar=True, cmap=pylab.cm.RdBu_r):
    """
    Linear heatmap of 2D residual array.

    Args:
        resid (array-like): Residual array to plot.
        resid_range (float): Values > resid_range or < -resid_range saturate the color
            spectrum.
        ax (int): Axes object to plot into. If None, the result of pylab.gca() is used.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        extend (str): Whether the colorbar should have 'extension' arrows. See
            help(pylab.colorbar) for more details.
        colorbar (bool): Should we plot a colorbar?
        cmap (pylab.cm): Pylab colormap to use for plotting.

    Returns:
        cb (figure element): The created colorbar.
    """
    if ax is None:
        ax = pylab.gca()

    if resid_range is None:
        resid_range = abs(resid).max()

    mappable=ax.pcolor(resid, cmap=cmap, vmin=-resid_range, 
                       vmax=resid_range, edgecolors='none')

    cbticks = [-resid_range, 0, resid_range]
    format = matplotlib.ticker.FormatStrFormatter('%.2g')
    cb = ax.figure.colorbar(mappable, ticks=cbticks, format=format,
                            extend=extend)
    if not colorbar:
        ax.figure.delaxes(ax.figure.axes[-1])
    else:
        try:
            ax.figure.dadi_colorbars.append(cb)
        except AttributeError:
            ax.figure.dadi_colorbars = [cb]

    ax.plot([0,resid.shape[1]],[0, resid.shape[0]], '-k', lw=0.2)

    if pop_ids is None:
        if resid.pop_ids is not None:
            pop_ids = resid.pop_ids
        else:
            pop_ids = ['pop0','pop1']
    ax.set_ylabel(pop_ids[0], verticalalignment='top')
    ax.set_xlabel(pop_ids[1], verticalalignment='bottom')

    ax.xaxis.set_major_formatter(_ctf)
    ax.xaxis.set_major_locator(_sfsTickLocator())
    ax.yaxis.set_major_formatter(_ctf)
    ax.yaxis.set_major_locator(_sfsTickLocator())
    for tick in ax.xaxis.get_ticklines() + ax.yaxis.get_ticklines():
        tick.set_visible(False)

    ax.set_xlim(0, resid.shape[1])
    ax.set_ylim(0, resid.shape[0])

    return cb

plot_3d_comp_Poisson(model, data, vmin=None, vmax=None, resid_range=None, fig_num=None, pop_ids=None, residual='Anscombe', adjust=True, show=True)

Poisson comparison between 3D model and data.

Parameters:

Name	Type	Description	Default
model	Spectrum	3-dimensional model SFS.	required
data	Spectrum	3-dimensional data SFS.	required
vmin	float	Minimum values plotted for sfs.	None
vmax	float	Maximum values plotted for sfs.	None
resid_range	float	Residual plot saturates at +- resid_range.	None
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
residual	str	'Anscombe' for Anscombe residuals, which are more normally distributed for Poisson sampling. 'linear' for the linear residuals, which can be less biased.	'Anscombe'
adjust	bool	Should method use automatic 'subplots_adjust'? For advanced manipulation of plots, it may be useful to make this False.	True
show	bool	If True, execute pylab.show command to make sure plot displays.	True

Source code in dadi/Plotting.py

def plot_3d_comp_Poisson(model, data, vmin=None, vmax=None,
                         resid_range=None, fig_num=None, pop_ids=None, 
                         residual='Anscombe', adjust=True, show=True):
    """
    Poisson comparison between 3D model and data.

    Args:
        model (Spectrum): 3-dimensional model SFS.
        data (Spectrum): 3-dimensional data SFS.
        vmin (float): Minimum values plotted for sfs.
        vmax (float): Maximum values plotted for sfs.
        resid_range (float): Residual plot saturates at +- resid_range.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        residual (str): 'Anscombe' for Anscombe residuals, which are more normally
            distributed for Poisson sampling. 'linear' for the linear residuals,
            which can be less biased.
        adjust (bool): Should method use automatic 'subplots_adjust'? For advanced
            manipulation of plots, it may be useful to make this False.
        show (bool): If True, execute pylab.show command to make sure plot displays.
    """
    if data.folded and not model.folded:
        model = model.fold()

    masked_model, masked_data = Numerics.intersect_masks(model, data)

    if fig_num is None:
        f = pylab.gcf()
    else:
        f = pylab.figure(fig_num, figsize=(8,10))

    pylab.clf()
    if adjust:
        pylab.subplots_adjust(bottom=0.07, left=0.07, top=0.95, right=0.95)

    modelmax = max(masked_model.sum(axis=sax).max() for sax in range(3))
    datamax = max(masked_data.sum(axis=sax).max() for sax in range(3))
    modelmin = min(masked_model.sum(axis=sax).min() for sax in range(3))
    datamin = min(masked_data.sum(axis=sax).min() for sax in range(3))
    max_toplot = max(modelmax, datamax)
    min_toplot = min(modelmin, datamin)

    if vmax is None:
        vmax = max_toplot
    if vmin is None:
        vmin = min_toplot
    extend = _extend_mapping[vmin <= min_toplot, vmax >= max_toplot]

    # Calculate the residuals
    if residual == 'Anscombe':
        resids = [Inference.\
                  Anscombe_Poisson_residual(masked_model.sum(axis=2-sax), 
                                            masked_data.sum(axis=2-sax), 
                                            mask=vmin) for sax in range(3)]
    elif residual == 'linear':
        resids =[Inference.\
                 linear_Poisson_residual(masked_model.sum(axis=2-sax), 
                                         masked_data.sum(axis=2-sax), 
                                         mask=vmin) for sax in range(3)]
    else:
        raise ValueError("Unknown class of residual '%s'." % residual)


    min_resid = min([r.min() for r in resids])
    max_resid = max([r.max() for r in resids])
    if resid_range is None:
        resid_range = max((abs(max_resid), abs(min_resid)))
    resid_extend = _extend_mapping[-resid_range <= min_resid, 
                                   resid_range >= max_resid]

    if pop_ids is not None:
        if len(pop_ids) != 3:
            raise ValueError('pop_ids must be of length 3.')
        data_ids = model_ids = resid_ids = pop_ids
    else:
        data_ids = masked_data.pop_ids
        model_ids = masked_model.pop_ids

        if model_ids is None:
            model_ids = data_ids

        if model_ids == data_ids:
           resid_ids = model_ids
        else:
            resid_ids = None

    for sax in range(3):
        marg_data = masked_data.sum(axis=2-sax)
        marg_model = masked_model.sum(axis=2-sax)

        curr_ids = []
        for ids in [data_ids, model_ids, resid_ids]:
            if ids is None:
                ids = ['pop0', 'pop1', 'pop2']

            if ids is not None:
                ids = list(ids)
                del ids[2-sax]

            curr_ids.append(ids)

        ax = pylab.subplot(4,3,sax+1)
        plot_colorbar = (sax == 2)
        plot_single_2d_sfs(marg_data, vmin=vmin, vmax=vmax, pop_ids=curr_ids[0],
                           extend=extend, colorbar=plot_colorbar, show=False)

        pylab.subplot(4,3,sax+4, sharex=ax, sharey=ax)
        plot_single_2d_sfs(marg_model, vmin=vmin, vmax=vmax, 
                           pop_ids=curr_ids[1], extend=extend, colorbar=False, show=False)

        resid = resids[sax]
        pylab.subplot(4,3,sax+7, sharex=ax, sharey=ax)
        plot_2d_resid(resid, resid_range, pop_ids=curr_ids[2],
                      extend=resid_extend, colorbar=plot_colorbar)

        ax = pylab.subplot(4,3,sax+10)
        flatresid = numpy.compress(numpy.logical_not(resid.mask.ravel()), 
                                   resid.ravel())
        ax.hist(flatresid, bins=20, density=True)
        ax.set_yticks([])
    if show:
        pylab.show()

plot_3d_comp_multinom(model, data, vmin=None, vmax=None, resid_range=None, fig_num=None, pop_ids=None, residual='Anscombe', adjust=True, show=True)

Multinomial comparison between 3D model and data.

Parameters:

Name	Type	Description	Default
model	Spectrum	3-dimensional model SFS.	required
data	Spectrum	3-dimensional data SFS.	required
vmin	float	Minimum values plotted for sfs.	None
vmax	float	Maximum values plotted for sfs.	None
resid_range	float	Residual plot saturates at +- resid_range.	None
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
residual	str	'Anscombe' for Anscombe residuals, which are more normally distributed for Poisson sampling. 'linear' for the linear residuals, which can be less biased.	'Anscombe'
adjust	bool	Should method use automatic 'subplots_adjust'? For advanced manipulation of plots, it may be useful to make this False.	True
show	bool	If True, execute pylab.show command to make sure plot displays.	True

Note

This comparison is multinomial in that it rescales the model to optimally fit the data.

Source code in dadi/Plotting.py

def plot_3d_comp_multinom(model, data, vmin=None, vmax=None,
                          resid_range=None, fig_num=None,
                          pop_ids=None, residual='Anscombe', adjust=True, show=True):
    """
    Multinomial comparison between 3D model and data.

    Args:
        model (Spectrum): 3-dimensional model SFS.
        data (Spectrum): 3-dimensional data SFS.
        vmin (float): Minimum values plotted for sfs.
        vmax (float): Maximum values plotted for sfs.
        resid_range (float): Residual plot saturates at +- resid_range.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        residual (str): 'Anscombe' for Anscombe residuals, which are more normally
            distributed for Poisson sampling. 'linear' for the linear residuals,
            which can be less biased.
        adjust (bool): Should method use automatic 'subplots_adjust'? For advanced
            manipulation of plots, it may be useful to make this False.
        show (bool): If True, execute pylab.show command to make sure plot displays.

    Note:
        This comparison is multinomial in that it rescales the model to
        optimally fit the data.
    """
    model = Inference.optimally_scaled_sfs(model, data)

    plot_3d_comp_Poisson(model, data, vmin=vmin, vmax=vmax,
                         resid_range=resid_range, fig_num=fig_num,
                         pop_ids=pop_ids, residual=residual,
                         adjust=adjust, show=show)

plot_3d_pairwise(data, vmin=None, vmax=None, fig_num=None, pop_ids=None, adjust=True, show=True)

Poisson comparison between 3D model and data.

Parameters:

Name	Type	Description	Default
data	Spectrum	3-dimensional data SFS.	required
vmin	float	Minimum values plotted for sfs.	None
vmax	float	Maximum values plotted for sfs.	None
fig_num	int	Clear and use figure fig_num for display. If None, a new figure window is created.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
adjust	bool	Should method use automatic 'subplots_adjust'? For advanced manipulation of plots, it may be useful to make this False.	True
show	bool	If True, execute pylab.show command to make sure plot displays.	True

Source code in dadi/Plotting.py

def plot_3d_pairwise(data, vmin=None, vmax=None,
                     fig_num=None, pop_ids=None, 
                     adjust=True, show=True):
    """
    Poisson comparison between 3D model and data.

    Args:
        data (Spectrum): 3-dimensional data SFS.
        vmin (float): Minimum values plotted for sfs.
        vmax (float): Maximum values plotted for sfs.
        fig_num (int): Clear and use figure fig_num for display. If None, a new figure
            window is created.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        adjust (bool): Should method use automatic 'subplots_adjust'? For advanced
            manipulation of plots, it may be useful to make this False.
        show (bool): If True, execute pylab.show command to make sure plot displays.
    """

    if fig_num is None:
        f = pylab.gcf()
    else:
        f = pylab.figure(fig_num, figsize=(10,4))

    pylab.clf()
    # pylab.tight_layout()
    # if adjust:
    #     pylab.subplots_adjust(bottom=0.07, left=0.07, top=0.95, right=0.95)
    max_toplot = max(data.sum(axis=sax).max() for sax in range(3))
    min_toplot = min(data.sum(axis=sax).min() for sax in range(3))

    if vmax is None:
        vmax = max_toplot
    if vmin is None:
        vmin = min_toplot
    extend = _extend_mapping[vmin <= min_toplot, vmax >= max_toplot]

    if pop_ids is not None:
        if len(pop_ids) != 3:
            raise ValueError('pop_ids must be of length 3.')
        data_ids = model_ids = resid_ids = pop_ids
    else:
        data_ids = data.pop_ids

    for sax in range(3):
        marg_data = data.sum(axis=2-sax)
        marg_model = data.sum(axis=2-sax)

        if data_ids is None:
            ids = ['pop0', 'pop1', 'pop2']
        else:
            ids = list(data_ids)
        del ids[2-sax]

        ax = pylab.subplot(1,3,sax+1)
        plot_colorbar = (sax == 2)
        plot_single_2d_sfs(marg_data, vmin=vmin, vmax=vmax, pop_ids=ids,
                           extend=extend, colorbar=plot_colorbar, show=False)

        # ax.set_yticks([])
    pylab.tight_layout()
    if show:
        pylab.show()

plot_3d_spectrum(fs, fignum=None, vmin=None, vmax=None, pop_ids=None, show=True)

Logarithmic heatmap of single 3D FS.

Note that this method is slow, because it relies on matplotlib's software rendering. For faster and better looking plots, use plot_3d_spectrum_mayavi.

Parameters:

Name	Type	Description	Default
fs	Spectrum	FS to plot.	required
vmin	float	Values in fs below vmin are masked in plot.	None
vmax	float	Values in fs above vmax saturate the color spectrum.	None
fignum	int	Figure number to plot into. If None, a new figure will be created.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
show	bool	If True, execute pylab.show command to make sure plot displays.	True

Source code in dadi/Plotting.py

def plot_3d_spectrum(fs, fignum=None, vmin=None, vmax=None, pop_ids=None,
                     show=True):
    """
    Logarithmic heatmap of single 3D FS.

    Note that this method is slow, because it relies on matplotlib's software
    rendering. For faster and better looking plots, use plot_3d_spectrum_mayavi.

    Args:
        fs (Spectrum): FS to plot.
        vmin (float): Values in fs below vmin are masked in plot.
        vmax (float): Values in fs above vmax saturate the color spectrum.
        fignum (int): Figure number to plot into. If None, a new figure will be created.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        show (bool): If True, execute pylab.show command to make sure plot displays.
    """
    import mpl_toolkits.mplot3d as mplot3d

    fig = pylab.figure(fignum)
    ax = mplot3d.Axes3D(fig)

    if vmin is None:
        vmin = fs.min()
    if vmax is None:
        vmax = fs.max()

    # Which entries should I plot?
    toplot = numpy.logical_not(fs.mask)
    toplot = numpy.logical_and(toplot, fs.data >= vmin)

    # Figure out the color mapping.
    normalized = (numpy.log(fs)-numpy.log(vmin))\
            /(numpy.log(vmax)-numpy.log(vmin))
    normalized = numpy.minimum(normalized, 1)
    colors = pylab.cm.hsv(normalized)

    # We draw by calculating which faces are visible and including each as a
    # polygon.
    polys, polycolors = [],[]
    for ii in range(fs.shape[0]):
        for jj in range(fs.shape[1]):
            for kk in range(fs.shape[2]):
                if not toplot[ii,jj,kk]:
                    continue
                if kk < fs.shape[2]-1 and toplot[ii,jj,kk+1]:
                    pass
                else:
                    polys.append([[ii-0.5,jj+0.5,kk+0.5],[ii+0.5,jj+0.5,kk+0.5],
                                  [ii+0.5,jj-0.5,kk+0.5],[ii-0.5,jj-0.5,kk+0.5]]
                                 )
                    polycolors.append(colors[ii,jj,kk])
                if kk > 0 and toplot[ii,jj,kk-1]:
                    pass
                else:
                    polys.append([[ii-0.5,jj+0.5,kk-0.5],[ii+0.5,jj+0.5,kk-0.5],
                                  [ii+0.5,jj-0.5,kk-0.5],[ii-0.5,jj-0.5,kk-0.5]]
                                 )
                    polycolors.append(colors[ii,jj,kk])
                if jj < fs.shape[1]-1 and toplot[ii,jj+1,kk]:
                    pass
                else:
                    polys.append([[ii-0.5,jj+0.5,kk+0.5],[ii+0.5,jj+0.5,kk+0.5],
                                  [ii+0.5,jj+0.5,kk-0.5],[ii-0.5,jj+0.5,kk-0.5]]
                                 )
                    polycolors.append(colors[ii,jj,kk])
                if jj > 0 and toplot[ii,jj-1,kk]:
                    pass
                else:
                    polys.append([[ii-0.5,jj-0.5,kk+0.5],[ii+0.5,jj-0.5,kk+0.5],
                                  [ii+0.5,jj-0.5,kk-0.5],[ii-0.5,jj-0.5,kk-0.5]]
                                 )
                    polycolors.append(colors[ii,jj,kk])
                if ii < fs.shape[0]-1 and toplot[ii+1,jj,kk]:
                    pass
                else:
                    polys.append([[ii+0.5,jj-0.5,kk+0.5],[ii+0.5,jj+0.5,kk+0.5],
                                  [ii+0.5,jj+0.5,kk-0.5],[ii+0.5,jj-0.5,kk-0.5]]
                                 )
                    polycolors.append(colors[ii,jj,kk])
                if ii > 0 and toplot[ii-1,jj,kk]:
                    pass
                else:
                    polys.append([[ii-0.5,jj-0.5,kk+0.5],[ii-0.5,jj+0.5,kk+0.5],
                                  [ii-0.5,jj+0.5,kk-0.5],[ii-0.5,jj-0.5,kk-0.5]]
                                 )
                    polycolors.append(colors[ii,jj,kk])


    polycoll = mplot3d.art3d.Poly3DCollection(polys, facecolor=polycolors, 
                                              edgecolor='k', linewidths=0.5)
    ax.add_collection(polycoll)

    # Set the limits
    ax.set_xlim3d(-0.5,fs.shape[0]-0.5)
    ax.set_ylim3d(-0.5,fs.shape[1]-0.5)
    ax.set_zlim3d(-0.5,fs.shape[2]-0.5)

    if pop_ids is None:
        if fs.pop_ids is not None:
            pop_ids = fs.pop_ids
        else:
            pop_ids = ['pop0','pop1','pop2']
    ax.set_xlabel(pop_ids[0], horizontalalignment='left')
    ax.set_ylabel(pop_ids[1], verticalalignment='bottom')
    ax.set_zlabel(pop_ids[2], verticalalignment='bottom')

    # XXX: I can't set the axis ticks to be just the endpoints.

    if show:
        pylab.show()

plot_3d_spectrum_mayavi(fs, fignum=None, vmin=None, vmax=None, pop_ids=None, show=True)

Logarithmic heatmap of single 3D FS.

This method relies on MayaVi2's mlab interface. See http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/mlab.html . To edit plot properties, click leftmost icon in the toolbar.

If you get an ImportError upon calling this function, it is likely that you don't have mayavi installed.

Parameters:

Name	Type	Description	Default
fs	Spectrum	FS to plot.	required
vmin	float	Values in fs below vmin are masked in plot.	None
vmax	float	Values in fs above vmax saturate the color spectrum.	None
fignum	int	Figure number to plot into. If None, a new figure will be created. Note that these are MayaVi figures, which are separate from matplotlib figures.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
show	bool	If True, execute mlab.show command to make sure plot displays.	True

Source code in dadi/Plotting.py

def plot_3d_spectrum_mayavi(fs, fignum=None, vmin=None, vmax=None, 
                            pop_ids=None, show=True):
    """
    Logarithmic heatmap of single 3D FS.

    This method relies on MayaVi2's mlab interface. See http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/mlab.html . To edit plot
    properties, click leftmost icon in the toolbar.

    If you get an ImportError upon calling this function, it is likely that you
    don't have mayavi installed.

    Args:
        fs (Spectrum): FS to plot.
        vmin (float): Values in fs below vmin are masked in plot.
        vmax (float): Values in fs above vmax saturate the color spectrum.
        fignum (int): Figure number to plot into. If None, a new figure will be created.
            Note that these are MayaVi figures, which are separate from
            matplotlib figures.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        show (bool): If True, execute mlab.show command to make sure plot displays.
    """
    from enthought.mayavi import mlab

    fig = mlab.figure(fignum, bgcolor=(1,1,1))
    mlab.clf(fig)

    if vmin is None:
        vmin = fs.min()
    if vmax is None:
        vmax = fs.max()

    # Which entries should I plot?
    toplot = numpy.logical_not(fs.mask)
    toplot = numpy.logical_and(toplot, fs.data >= vmin)

    # For the color mapping
    normalized = (numpy.log(fs)-numpy.log(vmin))\
            /(numpy.log(vmax)-numpy.log(vmin))
    normalized = numpy.minimum(normalized, 1)

    xs,ys,zs = numpy.indices(fs.shape)
    flat_xs = xs.flatten()
    flat_ys = ys.flatten()
    flat_zs = zs.flatten()
    flat_toplot = toplot.flatten()

    mlab.barchart(flat_xs[flat_toplot], flat_ys[flat_toplot], 
                  flat_zs[flat_toplot], normalized.flatten()[flat_toplot], 
                  colormap='hsv', scale_mode='none', lateral_scale=1, 
                  figure=fig)

    if pop_ids is None:
        if fs.pop_ids is not None:
            pop_ids = fs.pop_ids
        else:
            pop_ids = ['pop0','pop1','pop2']

    a = mlab.axes(xlabel=pop_ids[0],ylabel=pop_ids[1],zlabel=pop_ids[2], 
                  figure=fig, color=(0,0,0))
    a.axes.label_format = ""
    a.title_text_property.color = (0,0,0)
    mlab.text3d(fs.sample_sizes[0],fs.sample_sizes[1],fs.sample_sizes[2]+1, 
                '(%i,%i,%i)'%tuple(fs.sample_sizes), scale=0.75, figure=fig,
                color=(0,0,0))
    mlab.view(azimuth=-40, elevation=65, distance='auto', focalpoint='auto')

    if show:
        mlab.show()

plot_single_2d_sfs(sfs, vmin=None, vmax=None, ax=None, pop_ids=None, extend='neither', colorbar=True, cmap=pylab.cm.viridis_r, show=True)

Heatmap of single 2D SFS.

Parameters:

Name	Type	Description	Default
sfs	Spectrum	SFS to plot.	required
vmin	flaot	Values in sfs below vmin are masked in plot.	None
vmax	float	Values in sfs above vmax saturate the color spectrum.	None
ax	int	Axes object to plot into. If None, the result of pylab.gca() is used.	None
pop_ids	list[str]	If not None, override pop_ids stored in Spectrum.	None
extend	str	Whether the colorbar should have 'extension' arrows. See help(pylab.colorbar) for more details.	'neither'
colorbar	bool	Should we plot a colorbar?	True
cmap	pylab.cm function	Pylab colormap to use for plotting.	viridis_r
show	bool	If True, execute pylab.show command to make sure plot displays.	True

Returns:

Name	Type	Description
cb	cm	The created colorbar.

Source code in dadi/Plotting.py

def plot_single_2d_sfs(sfs, vmin=None, vmax=None, ax=None, 
                       pop_ids=None, extend='neither', colorbar=True,
                       cmap=pylab.cm.viridis_r, show=True):
    """
    Heatmap of single 2D SFS.

    Args:
        sfs (Spectrum): SFS to plot.
        vmin (flaot): Values in sfs below vmin are masked in plot.
        vmax (float): Values in sfs above vmax saturate the color spectrum.
        ax (int): Axes object to plot into. If None, the result of pylab.gca() is used.
        pop_ids (list[str]): If not None, override pop_ids stored in Spectrum.
        extend (str): Whether the colorbar should have 'extension' arrows. See
            help(pylab.colorbar) for more details.
        colorbar (bool): Should we plot a colorbar?
        cmap (pylab.cm function): Pylab colormap to use for plotting.
        show (bool): If True, execute pylab.show command to make sure plot displays.

    Returns:
        cb (pylab.cm): The created colorbar.
    """
    if ax is None:
        ax = pylab.gca()

    if vmin is None:
        vmin = sfs.min()
    if vmax is None:
        vmax = sfs.max()

    if vmax / vmin > 10:
        # Under matplotlib 1.0.1, default LogFormatter omits some tick lines.
        # This works more consistently.
        norm = matplotlib.colors.LogNorm(vmin=vmin*(1-1e-3), vmax=vmax*(1+1e-3))
        format = matplotlib.ticker.LogFormatterMathtext()
    else:
        norm = matplotlib.colors.Normalize(vmin=vmin*(1-1e-3), 
                                           vmax=vmax*(1+1e-3))
        format = None
    mappable=ax.pcolor(numpy.ma.masked_where(sfs<vmin, sfs), 
                       cmap=cmap, edgecolors='none',
                       norm=norm)
    cb = ax.figure.colorbar(mappable, extend=extend, format=format)
    if not colorbar:
        ax.figure.delaxes(ax.figure.axes[-1])
    else:
        # A hack so we can manually work around weird ticks in some colorbars
        try:
            ax.figure.dadi_colorbars.append(cb)
        except AttributeError:
            ax.figure.dadi_colorbars = [cb]

    ax.plot([0,sfs.shape[1]],[0, sfs.shape[0]], '-k', lw=0.2)

    if pop_ids is None:
        if sfs.pop_ids is not None:
            pop_ids = sfs.pop_ids
        else:
            pop_ids = ['pop0','pop1']
    ax.set_ylabel(pop_ids[0], verticalalignment='top')
    ax.set_xlabel(pop_ids[1], verticalalignment='bottom')

    ax.xaxis.set_major_formatter(_ctf)
    ax.xaxis.set_major_locator(_sfsTickLocator())
    ax.yaxis.set_major_formatter(_ctf)
    ax.yaxis.set_major_locator(_sfsTickLocator())
    for tick in ax.xaxis.get_ticklines() + ax.yaxis.get_ticklines():
        tick.set_visible(False)

    ax.set_xlim(0, sfs.shape[1])
    ax.set_ylim(0, sfs.shape[0])

    if show:
        pylab.show()

    return cb