latools/helpers/plot.py
def calibration_plot(self, analytes=None, datarange=True, loglog=False, ncol=3, srm_group=None, save=True):
    """
    Plot the calibration lines between measured and known SRM values.

    Parameters
    ----------
    analytes : optional, array_like or str
        The analyte(s) to plot. Defaults to all analytes.
    datarange : boolean
        Whether or not to show the distribution of the measured data
        alongside the calibration curve.
    loglog : boolean
        Whether or not to plot the data on a log - log scale. This is
        useful if you have two low standards very close together,
        and want to check whether your data are between them, or
        below them.

    Returns
    -------
    (fig, axes)
    """

    if isinstance(analytes, str):
        analytes = [analytes]

    if analytes is None:
        analytes = [a for a in self.analytes if self.internal_standard not in a]

    if srm_group is not None:
        srm_groups = {int(g): t for g, t in self.stdtab.loc[:, ['group', 'gTime']].values}
        try:
            gTime = srm_groups[srm_group]
        except KeyError:
            text = ('Invalid SRM group selection. Valid options are:\n' +
                    ' Key:  Time Centre\n' + 
                    '\n'.join(['   {:}:  {:.1f}s'.format(k, v) for k, v in srm_groups.items()]))
            print(text)
    else:
        gTime = None

    ncol = int(ncol)
    n = len(analytes)
    nrow = calc_nrow(n + 1, ncol)

    axes = []

    if not datarange:
        fig = plt.figure(figsize=[4.1 * ncol, 3 * nrow])
    else:
        fig = plt.figure(figsize=[4.7 * ncol, 3 * nrow])
        self.get_focus()

    gs = mpl.gridspec.GridSpec(nrows=int(nrow), ncols=int(ncol),
                            hspace=0.35, wspace=0.3)

    mdict = self.srm_mdict

    for g, a in zip(gs, analytes):
        if not datarange:
            ax = fig.add_axes(g.get_position(fig))
            axes.append((ax,))
        else:
            f = 0.8
            p0 = g.get_position(fig)
            p1 = [p0.x0, p0.y0, p0.width * f, p0.height]
            p2 = [p0.x0 + p0.width * f, p0.y0, p0.width * (1 - f), p0.height]
            ax = fig.add_axes(p1)
            axh = fig.add_axes(p2)
            axes.append((ax, axh))
        
        if gTime is None:
            sub = idx[a]
        else:
            sub = idx[a, :, :, gTime]
        x = self.srmtabs.loc[sub, 'meas_mean'].values
        xe = self.srmtabs.loc[sub, 'meas_err'].values
        y = self.srmtabs.loc[sub, 'srm_mean'].values
        ye = self.srmtabs.loc[sub, 'srm_err'].values
        srm = self.srmtabs.loc[sub].index.get_level_values('SRM')
        
        # plot calibration data
        for s, m in mdict.items():
            ind = srm == s
            ax.errorbar(x[ind], y[ind], xerr=xe[ind], yerr=ye[ind],
                        color=self.cmaps[a], alpha=0.6,
                        lw=0, elinewidth=1, marker=m, #'o',
                        capsize=0, markersize=5, label='_')

        # work out axis scaling
        if not loglog:
            xmax = np.nanmax(x + xe)
            ymax = np.nanmax(y + ye)
            if any(x - xe < 0):
                xmin = np.nanmin(x - xe)
                xpad = (xmax - xmin) * 0.05
                xlim = [xmin - xpad, xmax + xpad]
            else:
                xlim = [0, xmax * 1.05]

            if any(y - ye < 0):
                ymin = np.nanmin(y - ye)
                ypad = (ymax - ymin) * 0.05
                ylim = [ymin - ypad, ymax + ypad]
            else:
                ylim = [0, ymax * 1.05]

        else:
            xd = self.srmtabs.loc[a, 'meas_mean'][self.srmtabs.loc[a, 'meas_mean'] > 0].values
            yd = self.srmtabs.loc[a, 'srm_mean'][self.srmtabs.loc[a, 'srm_mean'] > 0].values

            xlim = [10**np.floor(np.log10(np.nanmin(xd))),
                    10**np.ceil(np.log10(np.nanmax(xd)))]
            ylim = [10**np.floor(np.log10(np.nanmin(yd))),
                    10**np.ceil(np.log10(np.nanmax(yd)))]

            # scale sanity checks
            if xlim[0] == xlim[1]:
                xlim[0] = ylim[0]
            if ylim[0] == ylim[1]:
                ylim[0] = xlim[0]

            ax.set_xscale('log')
            ax.set_yscale('log')

        ax.set_xlim(xlim)
        ax.set_ylim(ylim)

        # visual warning if any values < 0
        if xlim[0] < 0:
            ax.axvspan(xlim[0], 0, color=(1,0.8,0.8), zorder=-1)
        if ylim[0] < 0:
            ax.axhspan(ylim[0], 0, color=(1,0.8,0.8), zorder=-1)
        if any(x < 0) or any(y < 0):
            ax.text(.5, .5, 'WARNING: Values below zero.', color='r', weight='bold',
                    ha='center', va='center', rotation=40, transform=ax.transAxes, alpha=0.6)

        # calculate line and R2
        if loglog:
            x = np.logspace(*np.log10(xlim), 100)
        else:
            x = np.array(xlim)
        
        if gTime is None:
            coefs = self.calib_params.loc[:, a]
        else:
            coefs = self.calib_params.loc[gTime, a]
        
        m = np.nanmean(coefs['m'])
        m_nom = nominal_values(m)
        # calculate case-specific paramers
        if 'c' in coefs:
            c = np.nanmean(coefs['c'])
            c_nom = nominal_values(c)
            # calculate R2
            ym = self.srmtabs.loc[a, 'meas_mean'] * m_nom + c_nom
            R2 = R2calc(self.srmtabs.loc[a, 'srm_mean'], ym, force_zero=False)
            # generate line and label
            line = x * m_nom + c_nom
            label = 'y = {:.2e} x'.format(m)
            if c > 0:
                label += '\n+ {:.2e}'.format(c)
            else:
                label += '\n {:.2e}'.format(c)
        else:
            # calculate R2
            ym = self.srmtabs.loc[a, 'meas_mean'] * m_nom
            R2 = R2calc(self.srmtabs.loc[a, 'srm_mean'], ym, force_zero=True)
            # generate line and label
            line = x * m_nom
            label = 'y = {:.2e} x'.format(m)

        # plot line of best fit
        ax.plot(x, line, color=(0, 0, 0, 0.5), ls='dashed')

        # add R2 to label
        if round(R2, 3) == 1:
            label = '$R^2$: >0.999\n' + label
        else:
            label = '$R^2$: {:.3f}\n'.format(R2) + label

        ax.text(.05, .95, pretty_element(a), transform=ax.transAxes,
                weight='bold', va='top', ha='left', size=12)
        ax.set_xlabel('counts/counts ' + self.internal_standard)
        ax.set_ylabel('mol/mol ' + self.internal_standard)
        # write calibration equation on graph happens after data distribution

        # plot data distribution historgram alongside calibration plot
        if datarange:
            # isolate data
            meas = nominal_values(self.focus[a])
            meas = meas[~np.isnan(meas)]

            # check and set y scale
            if np.nanmin(meas) < ylim[0]:
                if loglog:
                    mmeas = meas[meas > 0]
                    ylim[0] = 10**np.floor(np.log10(np.nanmin(mmeas)))
                else:
                    ylim[0] = 0
                ax.set_ylim(ylim)

            m95 = np.percentile(meas[~np.isnan(meas)], 95) * 1.05
            if m95 > ylim[1]:
                if loglog:
                    ylim[1] = 10**np.ceil(np.log10(m95))
                else:
                    ylim[1] = m95

            # hist
            if loglog:
                bins = np.logspace(*np.log10(ylim), 30)
            else:
                bins = np.linspace(*ylim, 30)

            axh.hist(meas, bins=bins, orientation='horizontal',
                        color=self.cmaps[a], lw=0.5, alpha=0.5)

            if loglog:
                axh.set_yscale('log')
            axh.set_ylim(ylim)  # ylim of histogram axis
            ax.set_ylim(ylim)  # ylim of calibration axis
            axh.set_xticks([])
            axh.set_yticklabels([])

        # write calibration equation on graph
        cmax = np.nanmax(y)
        if cmax / ylim[1] > 0.5:
            ax.text(0.98, 0.04, label, transform=ax.transAxes,
                    va='bottom', ha='right')
        else:
            ax.text(0.02, 0.75, label, transform=ax.transAxes,
                    va='top', ha='left')

    if srm_group is None:
        title = 'All SRMs'
    else:
        title = 'SRM Group {:} (centre at {:.1f}s)'.format(srm_group, gTime)
    axes[0][0].set_title(title, loc='left', weight='bold', fontsize=12)
            
    # SRM legend
    ax = fig.add_axes(gs[-1].get_position(fig))
    for lab, m in mdict.items():
        ax.scatter([],[],marker=m, label=lab, color=(0,0,0,0.6))
    ax.legend()
    ax.axis('off')

    if save:
        fig.savefig(self.report_dir + '/calibration.pdf')

    return fig, axes
def calibration_plot(self, analytes=None, datarange=True, loglog=False, ncol=3, srm_group=None, save=True):
    """
    Plot the calibration lines between measured and known SRM values.

    Parameters
    ----------
    analytes : optional, array_like or str
        The analyte(s) to plot. Defaults to all analytes.
    datarange : boolean
        Whether or not to show the distribution of the measured data
        alongside the calibration curve.
    loglog : boolean
        Whether or not to plot the data on a log - log scale. This is
        useful if you have two low standards very close together,
        and want to check whether your data are between them, or
        below them.

    Returns
    -------
    (fig, axes)
    """

    if isinstance(analytes, str):
        analytes = [analytes]

    if analytes is None:
        analytes = [a for a in self.analytes if self.internal_standard not in a]

    if srm_group is not None:
        srm_groups = {int(g): t for g, t in self.stdtab.loc[:, ['group', 'gTime']].values}
        try:
            gTime = srm_groups[srm_group]
        except KeyError:
            text = ('Invalid SRM group selection. Valid options are:\n' +
                    ' Key:  Time Centre\n' + 
                    '\n'.join(['   {:}:  {:.1f}s'.format(k, v) for k, v in srm_groups.items()]))
            print(text)
    else:
        gTime = None

    ncol = int(ncol)
    n = len(analytes)
    nrow = calc_nrow(n + 1, ncol)

    axes = []

    if not datarange:
        fig = plt.figure(figsize=[4.1 * ncol, 3 * nrow])
    else:
        fig = plt.figure(figsize=[4.7 * ncol, 3 * nrow])
        self.get_focus()

    gs = mpl.gridspec.GridSpec(nrows=int(nrow), ncols=int(ncol),
                            hspace=0.35, wspace=0.3)

    mdict = self.srm_mdict

    for g, a in zip(gs, analytes):
        if not datarange:
            ax = fig.add_axes(g.get_position(fig))
            axes.append((ax,))
        else:
            f = 0.8
            p0 = g.get_position(fig)
            p1 = [p0.x0, p0.y0, p0.width * f, p0.height]
            p2 = [p0.x0 + p0.width * f, p0.y0, p0.width * (1 - f), p0.height]
            ax = fig.add_axes(p1)
            axh = fig.add_axes(p2)
            axes.append((ax, axh))
        
        if gTime is None:
            sub = idx[a]
        else:
            sub = idx[a, :, :, gTime]
        x = self.srmtabs.loc[sub, 'meas_mean'].values
        xe = self.srmtabs.loc[sub, 'meas_err'].values
        y = self.srmtabs.loc[sub, 'srm_mean'].values
        ye = self.srmtabs.loc[sub, 'srm_err'].values
        srm = self.srmtabs.loc[sub].index.get_level_values('SRM')
        
        # plot calibration data
        for s, m in mdict.items():
            ind = srm == s
            ax.errorbar(x[ind], y[ind], xerr=xe[ind], yerr=ye[ind],
                        color=self.cmaps[a], alpha=0.6,
                        lw=0, elinewidth=1, marker=m, #'o',
                        capsize=0, markersize=5, label='_')

        # work out axis scaling
        if not loglog:
            xmax = np.nanmax(x + xe)
            ymax = np.nanmax(y + ye)
            if any(x - xe < 0):
                xmin = np.nanmin(x - xe)
                xpad = (xmax - xmin) * 0.05
                xlim = [xmin - xpad, xmax + xpad]
            else:
                xlim = [0, xmax * 1.05]

            if any(y - ye < 0):
                ymin = np.nanmin(y - ye)
                ypad = (ymax - ymin) * 0.05
                ylim = [ymin - ypad, ymax + ypad]
            else:
                ylim = [0, ymax * 1.05]

        else:
            xd = self.srmtabs.loc[a, 'meas_mean'][self.srmtabs.loc[a, 'meas_mean'] > 0].values
            yd = self.srmtabs.loc[a, 'srm_mean'][self.srmtabs.loc[a, 'srm_mean'] > 0].values

            xlim = [10**np.floor(np.log10(np.nanmin(xd))),
                    10**np.ceil(np.log10(np.nanmax(xd)))]
            ylim = [10**np.floor(np.log10(np.nanmin(yd))),
                    10**np.ceil(np.log10(np.nanmax(yd)))]

            # scale sanity checks
            if xlim[0] == xlim[1]:
                xlim[0] = ylim[0]
            if ylim[0] == ylim[1]:
                ylim[0] = xlim[0]

            ax.set_xscale('log')
            ax.set_yscale('log')

        ax.set_xlim(xlim)
        ax.set_ylim(ylim)

        # visual warning if any values < 0
        if xlim[0] < 0:
            ax.axvspan(xlim[0], 0, color=(1,0.8,0.8), zorder=-1)
        if ylim[0] < 0:
            ax.axhspan(ylim[0], 0, color=(1,0.8,0.8), zorder=-1)
        if any(x < 0) or any(y < 0):
            ax.text(.5, .5, 'WARNING: Values below zero.', color='r', weight='bold',
                    ha='center', va='center', rotation=40, transform=ax.transAxes, alpha=0.6)

        # calculate line and R2
        if loglog:
            x = np.logspace(*np.log10(xlim), 100)
        else:
            x = np.array(xlim)
        
        if gTime is None:
            coefs = self.calib_params.loc[:, a]
        else:
            coefs = self.calib_params.loc[gTime, a]
        
        m = np.nanmean(coefs['m'])
        m_nom = nominal_values(m)
        # calculate case-specific paramers
        if 'c' in coefs:
            c = np.nanmean(coefs['c'])
            c_nom = nominal_values(c)
            # calculate R2
            ym = self.srmtabs.loc[a, 'meas_mean'] * m_nom + c_nom
            R2 = R2calc(self.srmtabs.loc[a, 'srm_mean'], ym, force_zero=False)
            # generate line and label
            line = x * m_nom + c_nom
            label = 'y = {:.2e} x'.format(m)
            if c > 0:
                label += '\n+ {:.2e}'.format(c)
            else:
                label += '\n {:.2e}'.format(c)
        else:
            # calculate R2
            ym = self.srmtabs.loc[a, 'meas_mean'] * m_nom
            R2 = R2calc(self.srmtabs.loc[a, 'srm_mean'], ym, force_zero=True)
            # generate line and label
            line = x * m_nom
            label = 'y = {:.2e} x'.format(m)

        # plot line of best fit
        ax.plot(x, line, color=(0, 0, 0, 0.5), ls='dashed')

        # add R2 to label
        if round(R2, 3) == 1:
            label = '$R^2$: >0.999\n' + label
        else:
            label = '$R^2$: {:.3f}\n'.format(R2) + label

        ax.text(.05, .95, pretty_element(a), transform=ax.transAxes,
                weight='bold', va='top', ha='left', size=12)
        ax.set_xlabel('counts/counts ' + self.internal_standard)
        ax.set_ylabel('mol/mol ' + self.internal_standard)
        # write calibration equation on graph happens after data distribution

        # plot data distribution historgram alongside calibration plot
        if datarange:
            # isolate data
            meas = nominal_values(self.focus[a])
            meas = meas[~np.isnan(meas)]

            # check and set y scale
            if np.nanmin(meas) < ylim[0]:
                if loglog:
                    mmeas = meas[meas > 0]
                    ylim[0] = 10**np.floor(np.log10(np.nanmin(mmeas)))
                else:
                    ylim[0] = 0
                ax.set_ylim(ylim)

            m95 = np.percentile(meas[~np.isnan(meas)], 95) * 1.05
            if m95 > ylim[1]:
                if loglog:
                    ylim[1] = 10**np.ceil(np.log10(m95))
                else:
                    ylim[1] = m95

            # hist
            if loglog:
                bins = np.logspace(*np.log10(ylim), 30)
            else:
                bins = np.linspace(*ylim, 30)

            axh.hist(meas, bins=bins, orientation='horizontal',
                        color=self.cmaps[a], lw=0.5, alpha=0.5)

            if loglog:
                axh.set_yscale('log')
            axh.set_ylim(ylim)  # ylim of histogram axis
            ax.set_ylim(ylim)  # ylim of calibration axis
            axh.set_xticks([])
            axh.set_yticklabels([])

        # write calibration equation on graph
        cmax = np.nanmax(y)
        if cmax / ylim[1] > 0.5:
            ax.text(0.98, 0.04, label, transform=ax.transAxes,
                    va='bottom', ha='right')
        else:
            ax.text(0.02, 0.75, label, transform=ax.transAxes,
                    va='top', ha='left')

    if srm_group is None:
        title = 'All SRMs'
    else:
        title = 'SRM Group {:} (centre at {:.1f}s)'.format(srm_group, gTime)
    axes[0][0].set_title(title, loc='left', weight='bold', fontsize=12)
            
    # SRM legend
    ax = fig.add_axes(gs[-1].get_position(fig))
    for lab, m in mdict.items():
        ax.scatter([],[],marker=m, label=lab, color=(0,0,0,0.6))
    ax.legend()
    ax.axis('off')

    if save:
        fig.savefig(self.report_dir + '/calibration.pdf')

    return fig, axes
Plot the calibration lines between measured and known SRM values.

    Parameters
    ----------
    analytes : optional, array_like or str
        The analyte(s) to plot. Defaults to all analytes.
    datarange : boolean
        Whether or not to show the distribution of the measured data
        alongside the calibration curve.
    loglog : boolean
        Whether or not to plot the data on a log - log scale. This is
        useful if you have two low standards very close together,
        and want to check whether your data are between them, or
        below them.

    Returns
    -------
    (fig, axes)