latools/latools.py
def calibrate(self, analytes=None, drift_correct=True,
                  srms_used=['NIST610', 'NIST612', 'NIST614'],
                  zero_intercept=True, n_min=10, reload_srm_database=False):
        """
        Calibrates the data to measured SRM values.

        Assumes that y intercept is zero.

        Parameters
        ----------  
        analytes : str or iterable
            Which analytes you'd like to calibrate. Defaults to all.
        drift_correct : bool
            Whether to pool all SRM measurements into a single calibration,
            or vary the calibration through the run, interpolating
            coefficients between measured SRMs.
        srms_used : str or iterable
            Which SRMs have been measured. Must match names given in
            SRM data file *exactly*.
        n_min : int
            The minimum number of data points an SRM measurement
            must have to be included.

        Returns
        -------
        None
        """
        if analytes is None:
            analytes = self.analytes[self.analytes != self.internal_standard]
        elif isinstance(analytes, str):
            analytes = [analytes]

        if not hasattr(self, 'srmtabs'):
            self.srm_id_auto(srms_used=srms_used, n_min=n_min, reload_srm_database=reload_srm_database)

        # make container for calibration params
        if not hasattr(self, 'calib_params'):
            gTime = self.stdtab.gTime.unique()
            self.calib_params = pd.DataFrame(columns=pd.MultiIndex.from_product([analytes, ['m']]),
                                            index=gTime)

        calib_analytes = self.srmtabs.index.get_level_values(0).unique()

        if zero_intercept:
            fn  = lambda x, m: x * m
        else:
            fn = lambda x, m, c: x * m + c

        for a in calib_analytes:
            if zero_intercept:
                if (a, 'c') in self.calib_params:
                    self.calib_params.drop((a, 'c'), 1, inplace=True)
            if drift_correct:
                for g in self.stdtab.gTime.unique():
                    ind = idx[a, :, :, g]
                    if self.srmtabs.loc[ind].size == 0:
                        continue
                    # try:
                    meas = self.srmtabs.loc[ind, 'meas_mean']
                    srm = self.srmtabs.loc[ind, 'srm_mean']
                    # TODO: replace curve_fit with Sambridge's 2D likelihood function for better uncertainty incorporation.
                    merr = self.srmtabs.loc[ind, 'meas_err']
                    serr = self.srmtabs.loc[ind, 'srm_err']
                    sigma = np.sqrt(merr**2 + serr**2)

                    if len(meas) > 1:
                        # multiple SRMs - do a regression
                        p, cov = curve_fit(fn, meas, srm, sigma=sigma)
                        pe = unc.correlated_values(p, cov)                
                        self.calib_params.loc[g, (a, 'm')] = pe[0]
                        if not zero_intercept:
                            self.calib_params.loc[g, (a, 'c')] = pe[1]
                    else:
                        # deal with case where there's only one datum
                        self.calib_params.loc[g, (a, 'm')] = (un.uarray(srm, serr) / 
                                                              un.uarray(meas, merr))[0]
                        if not zero_intercept:
                            self.calib_params.loc[g, (a, 'c')] = 0

                    # This should be obsolete, because no-longer sourcing locator from calib_params index.
                    # except KeyError:
                    #     # If the calibration is being recalculated, calib_params
                    #     # will have t=0 and t=max(uTime) values that are outside
                    #     # the srmtabs index.
                    #     # If this happens, drop them, and re-fill them at the end.
                    #     self.calib_params.drop(g, inplace=True)
            else:
                ind = idx[a, :, :, :]
                meas = self.srmtabs.loc[ind, 'meas_mean']
                srm = self.srmtabs.loc[ind, 'srm_mean']
                merr = self.srmtabs.loc[ind, 'meas_err']
                serr = self.srmtabs.loc[ind, 'srm_err']
                sigma = np.sqrt(merr**2 + serr**2)
                
                if len(meas) > 1:
                    p, cov = curve_fit(fn, meas, srm, sigma=sigma)
                    pe = unc.correlated_values(p, cov)                
                    self.calib_params.loc[:, (a, 'm')] = pe[0]
                    if not zero_intercept:
                        self.calib_params.loc[:, (a, 'c')] = pe[1]
                else:
                    self.calib_params.loc[:, (a, 'm')] = (un.uarray(srm, serr) / 
                                                          un.uarray(meas, merr))[0]
                    if not zero_intercept:
                        self.calib_params.loc[:, (a, 'c')] = 0

        # if fill:
        # fill in uTime=0 and uTime = max cases for interpolation
        if self.calib_params.index.min() == 0:
            self.calib_params.drop(0, inplace=True)
            self.calib_params.drop(self.calib_params.index.max(), inplace=True)
        self.calib_params.loc[0, :] = self.calib_params.loc[self.calib_params.index.min(), :]
        maxuT = np.max([d.uTime.max() for d in self.data.values()])  # calculate max uTime
        self.calib_params.loc[maxuT, :] = self.calib_params.loc[self.calib_params.index.max(), :]
        # sort indices for slice access
        self.calib_params.sort_index(1, inplace=True)
        self.calib_params.sort_index(0, inplace=True)

        # calculcate interpolators for applying calibrations
        self.calib_ps = Bunch()
        for a in analytes:
            # TODO: revisit un_interp1d to see whether it plays well with correlated values. 
            # Possible re-write to deal with covariance matrices?
            self.calib_ps[a] = {'m': un_interp1d(self.calib_params.index.values,
                                                self.calib_params.loc[:, (a, 'm')].values)}
            if not zero_intercept:
                self.calib_ps[a]['c'] = un_interp1d(self.calib_params.index.values,
                                                    self.calib_params.loc[:, (a, 'c')].values)

        with self.pbar.set(total=len(self.data), desc='Applying Calibrations') as prog:
            for d in self.data.values():
                d.calibrate(self.calib_ps, analytes)
                prog.update()

        # record SRMs used for plotting
        markers = 'osDsv<>PX'  # for future implementation of SRM-specific markers.
        if not hasattr(self, 'srms_used'):
            self.srms_used = set(srms_used)
        else:
            self.srms_used.update(srms_used)
        self.srm_mdict = {k: markers[i] for i, k in enumerate(self.srms_used)}

        self.stages_complete.update(['calibrated'])
        self.focus_stage = 'calibrated'

        return
def calibrate(self, analytes=None, drift_correct=True,
                  srms_used=['NIST610', 'NIST612', 'NIST614'],
                  zero_intercept=True, n_min=10, reload_srm_database=False):
        """
        Calibrates the data to measured SRM values.

        Assumes that y intercept is zero.

        Parameters
        ----------  
        analytes : str or iterable
            Which analytes you'd like to calibrate. Defaults to all.
        drift_correct : bool
            Whether to pool all SRM measurements into a single calibration,
            or vary the calibration through the run, interpolating
            coefficients between measured SRMs.
        srms_used : str or iterable
            Which SRMs have been measured. Must match names given in
            SRM data file *exactly*.
        n_min : int
            The minimum number of data points an SRM measurement
            must have to be included.

        Returns
        -------
        None
        """
        if analytes is None:
            analytes = self.analytes[self.analytes != self.internal_standard]
        elif isinstance(analytes, str):
            analytes = [analytes]

        if not hasattr(self, 'srmtabs'):
            self.srm_id_auto(srms_used=srms_used, n_min=n_min, reload_srm_database=reload_srm_database)

        # make container for calibration params
        if not hasattr(self, 'calib_params'):
            gTime = self.stdtab.gTime.unique()
            self.calib_params = pd.DataFrame(columns=pd.MultiIndex.from_product([analytes, ['m']]),
                                            index=gTime)

        calib_analytes = self.srmtabs.index.get_level_values(0).unique()

        if zero_intercept:
            fn  = lambda x, m: x * m
        else:
            fn = lambda x, m, c: x * m + c

        for a in calib_analytes:
            if zero_intercept:
                if (a, 'c') in self.calib_params:
                    self.calib_params.drop((a, 'c'), 1, inplace=True)
            if drift_correct:
                for g in self.stdtab.gTime.unique():
                    ind = idx[a, :, :, g]
                    if self.srmtabs.loc[ind].size == 0:
                        continue
                    # try:
                    meas = self.srmtabs.loc[ind, 'meas_mean']
                    srm = self.srmtabs.loc[ind, 'srm_mean']
                    # TODO: replace curve_fit with Sambridge's 2D likelihood function for better uncertainty incorporation.
                    merr = self.srmtabs.loc[ind, 'meas_err']
                    serr = self.srmtabs.loc[ind, 'srm_err']
                    sigma = np.sqrt(merr**2 + serr**2)

                    if len(meas) > 1:
                        # multiple SRMs - do a regression
                        p, cov = curve_fit(fn, meas, srm, sigma=sigma)
                        pe = unc.correlated_values(p, cov)                
                        self.calib_params.loc[g, (a, 'm')] = pe[0]
                        if not zero_intercept:
                            self.calib_params.loc[g, (a, 'c')] = pe[1]
                    else:
                        # deal with case where there's only one datum
                        self.calib_params.loc[g, (a, 'm')] = (un.uarray(srm, serr) / 
                                                              un.uarray(meas, merr))[0]
                        if not zero_intercept:
                            self.calib_params.loc[g, (a, 'c')] = 0

                    # This should be obsolete, because no-longer sourcing locator from calib_params index.
                    # except KeyError:
                    #     # If the calibration is being recalculated, calib_params
                    #     # will have t=0 and t=max(uTime) values that are outside
                    #     # the srmtabs index.
                    #     # If this happens, drop them, and re-fill them at the end.
                    #     self.calib_params.drop(g, inplace=True)
            else:
                ind = idx[a, :, :, :]
                meas = self.srmtabs.loc[ind, 'meas_mean']
                srm = self.srmtabs.loc[ind, 'srm_mean']
                merr = self.srmtabs.loc[ind, 'meas_err']
                serr = self.srmtabs.loc[ind, 'srm_err']
                sigma = np.sqrt(merr**2 + serr**2)
                
                if len(meas) > 1:
                    p, cov = curve_fit(fn, meas, srm, sigma=sigma)
                    pe = unc.correlated_values(p, cov)                
                    self.calib_params.loc[:, (a, 'm')] = pe[0]
                    if not zero_intercept:
                        self.calib_params.loc[:, (a, 'c')] = pe[1]
                else:
                    self.calib_params.loc[:, (a, 'm')] = (un.uarray(srm, serr) / 
                                                          un.uarray(meas, merr))[0]
                    if not zero_intercept:
                        self.calib_params.loc[:, (a, 'c')] = 0

        # if fill:
        # fill in uTime=0 and uTime = max cases for interpolation
        if self.calib_params.index.min() == 0:
            self.calib_params.drop(0, inplace=True)
            self.calib_params.drop(self.calib_params.index.max(), inplace=True)
        self.calib_params.loc[0, :] = self.calib_params.loc[self.calib_params.index.min(), :]
        maxuT = np.max([d.uTime.max() for d in self.data.values()])  # calculate max uTime
        self.calib_params.loc[maxuT, :] = self.calib_params.loc[self.calib_params.index.max(), :]
        # sort indices for slice access
        self.calib_params.sort_index(1, inplace=True)
        self.calib_params.sort_index(0, inplace=True)

        # calculcate interpolators for applying calibrations
        self.calib_ps = Bunch()
        for a in analytes:
            # TODO: revisit un_interp1d to see whether it plays well with correlated values. 
            # Possible re-write to deal with covariance matrices?
            self.calib_ps[a] = {'m': un_interp1d(self.calib_params.index.values,
                                                self.calib_params.loc[:, (a, 'm')].values)}
            if not zero_intercept:
                self.calib_ps[a]['c'] = un_interp1d(self.calib_params.index.values,
                                                    self.calib_params.loc[:, (a, 'c')].values)

        with self.pbar.set(total=len(self.data), desc='Applying Calibrations') as prog:
            for d in self.data.values():
                d.calibrate(self.calib_ps, analytes)
                prog.update()

        # record SRMs used for plotting
        markers = 'osDsv<>PX'  # for future implementation of SRM-specific markers.
        if not hasattr(self, 'srms_used'):
            self.srms_used = set(srms_used)
        else:
            self.srms_used.update(srms_used)
        self.srm_mdict = {k: markers[i] for i, k in enumerate(self.srms_used)}

        self.stages_complete.update(['calibrated'])
        self.focus_stage = 'calibrated'

        return
Calibrates the data to measured SRM values.

        Assumes that y intercept is zero.

        Parameters
        ----------  
        analytes : str or iterable
            Which analytes you'd like to calibrate. Defaults to all.
        drift_correct : bool
            Whether to pool all SRM measurements into a single calibration,
            or vary the calibration through the run, interpolating
            coefficients between measured SRMs.
        srms_used : str or iterable
            Which SRMs have been measured. Must match names given in
            SRM data file *exactly*.
        n_min : int
            The minimum number of data points an SRM measurement
            must have to be included.

        Returns
        -------
        None