def _set_age_values(self, f, include_decay_error=False):
        arc = self.arar_constants
        j = copy(self.j)
        if j is None:
            j = ufloat(1e-4, 1e-7)
        j.tag = 'Position'
        j.std_dev = self.position_jerr or 0
        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage_w_position_err = age

        j = self.j
        if j is None:
            j = ufloat(1e-4, 1e-7, tag='J')

        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage_w_j_err = age

        j = copy(self.j)
        if j is None:
            j = ufloat(1e-4, 1e-7, tag='J')

        j.std_dev = 0
        age = age_equation(j, f, include_decay_error=include_decay_error, arar_constants=arc)
        self.uage = age

        self.age = nominal_value(age)
        self.age_err = std_dev(age)
        self.age_err_wo_j = std_dev(age)

        for iso in self.itervalues():
            iso.age_error_component = self.get_error_component(iso.name)

def test_counter():
    if (sys.version_info < (2, 6, 0)):
        from nose.plugins.skip import SkipTest
        raise SkipTest    
    mcu = Arduino()
    reg = mcu.registers.proxy
    mcu.pins.reset()
    p = mcu.pin(5)
    p.write_mode(OUTPUT)
    p.pwm.write_value(128)
    print 'frequencies_available:', p.pwm.frequencies_available
    for fset in p.pwm.frequencies_available:
        p.pwm.frequency = fset
        assert abs(p.pwm.frequency - fset) <= 1
        print '---------------------------'
        print 'fset=', fset
        print '---------------------------'
        for ms in [10, 20, 50, 100, 200, 500, 1000]:
            for _ in range(1):
                t = ms / 1000.0
                with mcu.counter:
                    mcu.counter.run(t)
                    f = mcu.counter.frequency
                    t = mcu.counter.gate_time
                    err = f - fset

                    print 't=%s  f=%s ' % (t, f)
                    ok_(abs(nominal_value(err)) <= 0.1+std_dev(err), 
                        (abs(nominal_value(err)),std_dev(err)))

    def _load_unknown_computed(self, an, new_list):
        attrs = (('Age', 'uage_w_j_err'),
                 # ('Age', 'age', None, None, 'age_err'),
                 ('w/o J', 'wo_j', '', 'uage', 'age_err_wo_j'),
                 ('K/Ca', 'kca'),
                 ('K/Cl', 'kcl'),
                 ('40Ar*', 'rad40_percent'),
                 ('F', 'uF'),
                 ('w/o Irrad', 'wo_irrad', '', 'uF', 'F_err_wo_irrad'))

        if new_list:
            def comp_factory(n, a, value=None, value_tag=None, error_tag=None):
                if value is None:
                    value = getattr(an, a)

                display_value = True
                if value_tag:
                    value = getattr(an, value_tag)
                    display_value = False

                if error_tag:
                    e = getattr(an, error_tag)
                else:
                    e = std_dev(value)

                return ComputedValue(name=n,
                                     tag=a,
                                     value=nominal_value(value) or 0,
                                     value_tag=value_tag or '',
                                     display_value=display_value,
                                     error=e or 0)

            cv = [comp_factory(*args)
                  for args in attrs]

            self.computed_values = cv
        else:
            age = an.uage
            nage, sage = nominal_value(age), std_dev(age)
            try:
                self.summary_str = u'Age={} {}{}({}%)'.format(floatfmt(nage), PLUSMINUS,
                                                              floatfmt(sage), format_percent_error(nage, sage))
            except:
                pass

            for ci in self.computed_values:
                attr = ci.tag
                if attr == 'wo_j':
                    ci.error = an.age_err_wo_j or 0
                    ci.value = nominal_value(getattr(an, ci.value_tag))
                elif attr == 'wo_irrad':
                    ci.error = an.F_err_wo_irrad or 0
                    ci.value = nominal_value(getattr(an, ci.value_tag))
                else:
                    v = getattr(an, attr)
                    if v is not None:
                        ci.value = nominal_value(v)
                        ci.error = std_dev(v)

    def _value_string(self, t):
        if t == 'uF':
            a, e = self.F, self.F_err
        elif t == 'uage':
            a, e = nominal_value(self.uage), std_dev(self.uage)
        else:
            v = self.get_value(t)
            if isinstance(v, Isotope):
                v = v.get_intensity()
            a, e = nominal_value(v), std_dev(v)

        return a, e

    def _make_intermediate_summary(self, sh, ag, cols, label):
        row = self._current_row

        age_idx = next((i for i, c in enumerate(cols) if c.label == 'Age'), 0)
        cum_idx = next((i for i, c in enumerate(cols) if c.attr == 'cumulative_ar39'), 0)

        fmt = self._get_number_format('summary_age')
        kcafmt = self._get_number_format('summary_kca')

        fmt.set_bottom(1)
        kcafmt.set_bottom(1)

        fmt2 = self._workbook.add_format({'bottom': 1, 'bold': True})
        border = self._workbook.add_format({'bottom': 1})

        for i in range(age_idx + 1):
            sh.write_blank(row, i, '', fmt)

        startcol = 1
        sh.write(row, startcol, '{:02n}'.format(ag.aliquot), fmt2)
        sh.write_rich_string(row, startcol + 1, label, fmt2)
        cols[startcol + 1].calculate_width(label)

        age = ag.uage
        tn = ag.total_n
        if label == 'plateau':
            if not ag.plateau_steps:
                age = None
            else:
                txt = 'n={}/{} steps={}'.format(ag.nsteps, tn, ag.plateau_steps_str)
                sh.write(row, startcol + 2, txt, border)
                sh.write(row, cum_idx + 1, format_mswd(ag.get_plateau_mswd_tuple()), border)

        else:
            txt = 'n={}/{}'.format(ag.nanalyses, tn)
            sh.write(row, startcol + 2, txt, border)
            sh.write(row, cum_idx + 1, format_mswd(ag.get_mswd_tuple()), border)

        if age is not None:
            sh.write_number(row, age_idx, nominal_value(age), fmt)
            sh.write_number(row, age_idx + 1, std_dev(age), fmt)
        else:
            sh.write(row, age_idx, 'No plateau', border)

        sh.write_number(row, age_idx + 2, nominal_value(ag.kca), kcafmt)
        sh.write_number(row, age_idx + 3, std_dev(ag.kca), kcafmt)

        if label == 'plateau':
            sh.write_number(row, cum_idx, ag.plateau_total_ar39(), fmt)
        else:
            sh.write_number(row, cum_idx, ag.valid_total_ar39(), fmt)
        self._current_row += 1

    def _get_error(self, attr, n=3, **kw):
        v = ""
        item = self.item
        if hasattr(item, attr):
            v = getattr(self.item, attr)
            if v:
                v = floatfmt(std_dev(v), n=n, **kw)
        elif hasattr(item, "isotopes"):
            if attr in item.isotopes:
                v = item.isotopes[attr].get_intensity()
                v = floatfmt(std_dev(v), n=n, **kw)

        return v

    def __init__(self, ratios, name, *args, **kw):
        super(Result, self).__init__(*args, **kw)
        vs = array([nominal_value(ri) for ri in ratios])
        es = array([std_dev(ri) for ri in ratios])

        self.name = name
        m = ratios.mean()
        self.value = nominal_value(m)
        self.error = std_dev(m)
        wm, we = calculate_weighted_mean(vs, es)
        self.wm_value = wm
        self.wm_error = we
        self.mswd = calculate_mswd(vs, es, wm=wm)

    def _air_ratio(self):
        a4038 = self.isotope_group.get_ratio('Ar40/Ar38', non_ic_corr=True)
        a4036 = self.isotope_group.get_ratio('Ar40/Ar36', non_ic_corr=True)
        # e4038 = uformat_percent_error(a4038, include_percent_sign=True)
        # e4036 = uformat_percent_error(a4036, include_percent_sign=True)

        lines = [self._make_header('Ratios'),
                 'Ar40/Ar36= {} {}'.format(floatfmt(nominal_value(a4036)), errorfmt(nominal_value(a4036),
                                                                                    std_dev(a4036))),
                 'Ar40/Ar38= {} {}'.format(floatfmt(nominal_value(a4038)), errorfmt(nominal_value(a4038),
                                                                                    std_dev(a4038)))]

        return self._make_lines(lines)

def test_access_to_std_dev():
    "Uniform access to the standard deviation"

    x = ufloat((1, 0.1))
    y = 2*x

    # std_dev for Variable and AffineScalarFunc objects:
    assert uncertainties.std_dev(x) == x.std_dev()
    assert uncertainties.std_dev(y) == y.std_dev()

    # std_dev for other objects:
    assert uncertainties.std_dev([]) == 0
    assert uncertainties.std_dev(None) == 0

	def __residenceTime(self):
		"""
			Calculate the residence time of a single step event. 
		"""
		# if previous errors were detected, the 
		# event is already rejected, don't process it 
		# any further
		if self.mdProcessingStatus != 'normal':
			return

		# set numpy warning handling. 
		# raise divide by zero errors so we 
		# can catch them
		np.seterr(divide='raise')

		ocmu=np.abs(uncertainties.nominal_value(self.mdOpenChCurrent))
		ocsd=np.abs(uncertainties.std_dev(self.mdOpenChCurrent))
		
		bcmu=np.abs(uncertainties.nominal_value(self.mdBlockedCurrent))
		bcsd=np.abs(uncertainties.std_dev(self.mdBlockedCurrent))

		# refine the start estimate
		idx=self.eStartEstimate

		try:
			while np.abs((np.abs(self.eventData[idx])-ocmu)/ocsd) > 5.0:
				idx-=1
			
			# Set the start point
			self.mdEventStart=idx+1

			# Next move the count forward so we are in the blocked channel region of the pulse
			while np.abs((np.abs(self.eventData[idx])-bcmu)/bcsd) > 0.5:
				idx+=1

			# Search for the event end. 7*sigma allows us to prevent false 
			# positives
			while np.abs((np.abs(self.eventData[idx])-bcmu)/bcsd) < 7.0:
				idx+=1

			# Finally backtrack to find the true event end
			while np.abs((np.abs(self.eventData[idx])-bcmu)/bcsd) > 5.0:
				idx-=1
		except ( IndexError, FloatingPointError ):
			self.rejectEvent('eResTime')
			return

		self.mdEventEnd=idx-1

		# residence time in ms
		self.mdResTime=1000.*((self.mdEventEnd-self.mdEventStart)/float(self.Fs))

    def load_measurement(self, an, ar):

        # j = self._get_j(an)
        j = ar.j
        jf = 'NaN'
        if j is not None:
            jj = floatfmt(nominal_value(j), n=7, s=5)
            pe = format_percent_error(nominal_value(j), std_dev(j), include_percent_sign=True)
            jf = u'{} \u00b1{:0.2e}({})'.format(jj, std_dev(j), pe)

        a39 = ar.ar39decayfactor
        a37 = ar.ar37decayfactor
        ms = [MeasurementValue(name='Branch',
                               value=an.branch),
              MeasurementValue(name='DAQ Version',
                               value=an.collection_version),
              MeasurementValue(name='UUID',
                               value=an.uuid),
              MeasurementValue(name='RepositoryID',
                               value=an.repository_identifier),
              MeasurementValue(name='Spectrometer',
                               value=an.mass_spectrometer),
              MeasurementValue(name='Run Date',
                               value=an.rundate.strftime('%Y-%m-%d %H:%M:%S')),
              MeasurementValue(name='Irradiation',
                               value=self._get_irradiation(an)),
              MeasurementValue(name='J',
                               value=jf),
              MeasurementValue(name='Position Error',
                               value=floatfmt(an.position_jerr, use_scientific=True)),
              MeasurementValue(name='Lambda K',
                               value=nominal_value(ar.arar_constants.lambda_k),
                               units='1/a'),
              MeasurementValue(name='Project',
                               value=an.project),
              MeasurementValue(name='Sample',
                               value=an.sample),
              MeasurementValue(name='Material',
                               value=an.material),
              MeasurementValue(name='Comment',
                               value=an.comment),
              MeasurementValue(name='Ar39Decay',
                               value=floatfmt(a39)),
              MeasurementValue(name='Ar37Decay',
                               value=floatfmt(a37)),
              MeasurementValue(name='Sens.',
                               value=floatfmt(an.sensitivity, use_scientific=True),
                               units=an.sensitivity_units)]

        self.measurement_values = ms

    def _update_ratios(self, an):
        for ci in self.computed_values:
            nd = ci.detectors
            n, d = nd.split('/')
            niso, diso = self._get_isotope(n), self._get_isotope(d)

            noncorrected = self._get_non_corrected_ratio(niso, diso)
            corrected, ic = self._get_corrected_ratio(niso, diso)

            ci.trait_set(value=floatfmt(nominal_value(corrected)),
                         error=floatfmt(std_dev(corrected)),
                         noncorrected_value=nominal_value(noncorrected),
                         noncorrected_error=std_dev(noncorrected),
                         ic_factor=nominal_value(ic))

    def _monte_carlo_error_propagation(self, vr, m):
        lambda_total = self._lambda_t
        el = self._lambda_ec
        f = self._f

        vel = nominal_value(el) + std_dev(el) * randn(self._n)
        vt = nominal_value(lambda_total) + std_dev(lambda_total) * randn(self._n)
        vf = nominal_value(f) + std_dev(f) * randn(self._n)

        vt_mc = ones(1, m) * vt
        vf_mc = ones(1, m) * vf
        vel_mc = ones(1, m) * vel

        t_mc = log(vt_mc / vel_mc * vf_mc * vr + 1) / vt_mc
        return mean(t_mc), std(t_mc)

    def _linear_error_propagation(self, age, r, sr):
        """
        age in years
        :param age:
        :param r:
        :return: linear error propagation age error in years
        """

        lambda_total = self._lambda_t
        b = self._lambda_b
        el = self._lambda_ec
        f = self._f

        # partial derivatives
        pd_el = -(1. / lambda_total) * (age + (b * f * r / ((el ** 2) * umath.exp(lambda_total * age))))
        pd_b = (1 / lambda_total) * ((f * r / (el * umath.exp(lambda_total * age))) - age)
        pd_f = r / (el * umath.exp(lambda_total * age))
        pd_r = f / (el * umath.exp(lambda_total * age))

        sel = std_dev(el)
        sb = std_dev(b)
        sf = std_dev(self._f)
        # sr = std_dev(r)

        # (partial derivatives x sigma) ** 2
        pd_el2 = (pd_el * sel) ** 2
        pd_b2 = (pd_b * sb) ** 2
        pd_f2 = (pd_f * sf) ** 2
        pd_r2 = (pd_r * sr) ** 2

        sum_pd = pd_el2 + pd_b2 + pd_f2 + pd_r2

        # covariances
        cov_f_el = 7.1903e-19
        cov_f_b = -6.5839e-19
        cov_el_b = -3.4711e-26

        cov_f_el2 = 2. * cov_f_el * pd_f * pd_el
        cov_f_b2 = 2. * cov_f_b * pd_f * pd_b
        cov_el_b = 2. * cov_el_b * pd_el * pd_b

        sum_cov = cov_f_el2 + cov_f_b2 + cov_el_b

        ss = sum_pd + sum_cov

        # uncertainty in age
        st = ss ** 0.5
        return nominal_value(st)

    def load_isotopes(self):
        isos = self.editor.model.isotopes
        ns = []
        bks = []
        bs = []
        ics = []
        dets = []
        for k in self.editor.model.isotope_keys:
            iso = isos[k]
            iso.use_static = True
            ns.append(iso)
            bks.append(iso.blank)
            bs.append(iso.baseline)

            det = iso.detector
            if not det in dets:
                v, e = nominal_value(iso.ic_factor), std_dev(iso.ic_factor)
                ics.append(ICFactor(value=v, error=e,
                                    ovalue=v, oerror=e,
                                    name=det))
            dets.append(det)

        self.isotopes = ns
        self.blanks = bks
        self.baselines = bs
        self.ic_factors = ics

 def _get_j_err(self):
     j = self.j
     try:
         e = (std_dev(j) / nominal_value(j)) if j is not None else 0
     except ZeroDivisionError:
         e = nan
     return e

    def add_irradiation_production(self, name, pr, ifc):
        kw = {}
        for k, v in ifc.iteritems():
            if k == 'cl3638':
                k = 'P36Cl38Cl'
            else:
                k = k.capitalize()

            kw[k] = float(nominal_value(v))
            kw['{}Er'.format(k)] = float(std_dev(v))

        kw['ClOverKMultiplier'] = pr['Cl_K']
        kw['ClOverKMultiplierEr'] = 0
        kw['CaOverKMultiplier'] = pr['Ca_K']
        kw['CaOverKMultiplierEr'] = 0
        v = binascii.crc32(''.join([str(v) for v in kw.itervalues()]))
        with self.session_ctx() as sess:
            q = sess.query(IrradiationProductionTable)
            q = q.filter(IrradiationProductionTable.ProductionRatiosID == v)
            if not self._query_one(q):
                i = IrradiationProductionTable(Label=name,
                                               ProductionRatiosID=v,
                                               **kw)
                self._add_item(i)
        return v

    def get_mean_raw(self, tau=None):
        vs = []
        corrfunc = self._deadtime_correct
        for r in six.itervalues(self._cp):
            n = int(r['NShots'])
            nv = ufloat(float(r['Ar40']), float(r['Ar40err'])) * 6240
            dv = ufloat(float(r['Ar36']), float(r['Ar36err'])) * 6240
            if tau:
                dv = corrfunc(dv, tau * 1e-9)

            vs.append((n, nv / dv))

        key = lambda x: x[0]
        vs = sorted(vs, key=key)

        mxs = []
        mys = []
        mes = []
        for n, gi in groupby(vs, key=key):
            mxs.append(n)
            ys, es = list(zip(*[(nominal_value(xi[1]), std_dev(xi[1])) for xi in gi]))

            wm, werr = calculate_weighted_mean(ys, es)
            mys.append(wm)
            mes.append(werr)

        return mxs, mys, mes

def small_sep13 ( mx , mxErr, lag=0):
	""" 
	Given a frequency matrix and errors calculates the (scaled) small 
	separation d13 and propagates errors.
	
	Notes
	-----
	The parameter lag is the difference between the radial orders of
	the first modes of l=1 and l=3.  
	"""

	(Nn, Nl) = np.shape(mx)
	d13    = np.zeros((1,Nn))
	d13Err = np.zeros((1,Nn))
	d13.fill(None)      # values that can't be calculated remain NaN

	for n in range(1-lag,Nn):
		if (mx[n,1] != 0. and mx[n-1+lag,3] != 0.):     
        		a = un.ufloat( (mx[n,1], mxErr[n,1]) )
        		b = un.ufloat( (mx[n-1+lag,3], mxErr[n-1+lag,3]) )
        		result = (a-b) / 5.
        		d13[0,n-1+lag]    = un.nominal_value(result)
        		d13Err[0,n-1+lag] = un.std_dev(result)

	return d13, d13Err

    def _calculate_integrated_age(self, ans, weighting):
        ret = ufloat(0, 0)
        if ans and all((not isinstance(a, InterpretedAgeGroup) for a in ans)):

            rs = array([a.get_computed_value('rad40') for a in ans])
            ks = array([a.get_computed_value('k39') for a in ans])
            sks = ks.sum()

            weights = None
            if weighting == 'Volume':
                weights = ks / sks
            elif weighting == 'Variance':
                weights = [1 / std_dev(k) ** 2 for k in rs/ks]

            if weights is not None:
                wmean, sum_weights = average([nominal_value(fi) for fi in rs / ks], weights=weights, returned=True)
                werr = sum_weights ** -0.5
                f = ufloat(wmean, werr)
            else:
                f = rs.sum() / sks

            a = ans[0]
            j = a.j
            try:
                ret = age_equation(f, j, a.arar_constants)  # / self.age_scalar
            except ZeroDivisionError:
                pass

        return ret