[SciPy-dev] value+error type (long)
Gary Ruben
gazzar at email.com
Tue Mar 16 08:19:38 EST 2004
I recently became fed up with performing error calculations on experimental data and decided to make a Python data type which I could use with Numeric. Now I can enter data with error along with the errors as in the following example and all the error calculations are automatically carried through for me. I thought there might be some interest amongst others for inclusion of this in the scipy distribution, although you might prefer it to be in C. I provide it as food for thought, but let me know if there is interest in including it as is and I can look at making it robust by adding testcases etc. and some nicer helper functions for constructors with relative errors etc. Note, the helper functions for extracting the limits would be better implemented as properties of an array type subclassed from the Numeric array type. I gather Numarray would allow this, but I couldn't think of a way to do it in Numeric and I wanted to maintain compatibility. Similarly, I couldn't think of a
nicer way to apply Ufuncs across the prime value and error extrema. It would be nice to not have to call ApplyUfuncToErr() to do this.
Here's a sample of usage of the class, followed by the class itself. Note I've used Gnuplot.py in this but you'll get the idea by looking at the code:
--
R_AB.py:
from Numeric import *
import Gnuplot, Gnuplot.funcutils
from ErrorVal import *
from scipy.interpolate import splrep,splev
# manometer readings [mmHg]
upperPressure = ArrayOfErr([909., 802., 677., 585., 560., 548.], 1.0)
lowerPressure = ArrayOfErr([144., 246., 378., 469., 493., 505.], 1.0)
pressureDiff = upperPressure - lowerPressure
# Temperatures & Pressures from lookup table
lowerTablePressures = array([760. , 540. , 290. , 110. , 60. , 40. ]) # [mmHg]
lowerTableTemps = array([4.210, 3.867, 3.332, 2.685, 2.373, 2.194]) # [K]
upperTablePressures = array([780. , 560. , 300. , 120. , 70. , 50. ]) # [mmHg]
upperTableTemps = array([4.238, 3.902, 3.358, 2.735, 2.447, 2.290]) # [K]
# Linearly interpolate table values to produce temperatures
temps = ( lowerTableTemps +
(pressureDiff - lowerTablePressures) *
(upperTableTemps - lowerTableTemps) /
(upperTablePressures - lowerTablePressures) )
# Voltages across R_AB [V]
V_RStandard = ArrayOfErr([2.016, 2.016, 2.020, 2.017, 2.021, 2.019], 0.001)
R = 100.2 # standard resistor value [Ohm]
I = V_RStandard / R # current [A]
# Voltages across Germanium resistor [V]
V_RAB = array([Err(6.167, 0.002),
Err(6.168, 0.002),
Err(6.170, 0.002),
Err(6.160, 0.002),
Err(6.153, 0.02),
Err(5.894, 0.01)], PyObject)
R_AB = V_RAB / I
logRAB = ApplyUfuncToErr(R_AB, log10) # This means log10(R_AB)
sqrtLogROnT = (logRAB/temps)**0.5
g = Gnuplot.Gnuplot(debug=1)
g.title('Allen-Bradley model log(R)=A+B[log(R)/T]^1/2')
g.ylabel('log(R)')
g.xlabel('[log(R)/T]^1/2')
g('set bar 1') # set bar end length
g('set grid')
# spline fit
x1 = PrimeVals(sqrtLogROnT)[0]
x2 = PrimeVals(sqrtLogROnT)[5]
splineX = arange(x1,x2+.005,.01)
exactSpline = splrep(PrimeVals(sqrtLogROnT), PrimeVals(logRAB), s=0)
exactSplineY = splev(splineX, exactSpline)
g.plot(\
Gnuplot.Data(PrimeVals(sqrtLogROnT), PrimeVals(logRAB), inline=0, with='linesp lt 1 pt 0'),
Gnuplot.Data(splineX, exactSplineY, inline=0, with='linesp lt 2 pt 0'), # exact spline
Gnuplot.Data(PrimeVals(sqrtLogROnT), # x-points
PrimeVals(logRAB), # y-points
MinVals(sqrtLogROnT), MaxVals(sqrtLogROnT), # x-axis error bars
MinVals(logRAB), MaxVals(logRAB), # y-axis error bars
inline=0, with='xyerrorbars -1 0'),
)
raw_input('Gary says copy to clipboard now\n')
--
ErrorVal.py:
from Numeric import *
class Err:
def __init__(self, value, posErr=0., negErr=None):
self.value = float(value)
self.posErr = float(posErr)
if negErr:
self.negErr = float(negErr)
else:
self.negErr = float(posErr)
def __repr__(self):
return "Err(%s,%s,%s)" % (self.value, self.posErr, self.negErr)
def __str__(self):
return "%s +%s/-%s" % (self.value, self.posErr, self.negErr)
def __abs__(self):
value = abs(self.value)
if value == self.value:
posErr = self.posErr
negErr = self.negErr
else:
negErr = self.posErr
posErr = self.negErr
return Err(value, posErr, negErr)
def __neg__(self):
return Err(-self.value, self.negErr, self.posErr)
def __pos__(self):
return Err(self.value, self.posErr, self.negErr)
def __add__(self, other):
value = self.value + other.value
posErr = self.posErr + other.posErr
negErr = self.negErr + other.negErr
return Err(value, posErr, negErr)
def __radd__(self, other):
value = self.value + other.value
posErr = self.posErr + other.posErr
negErr = self.negErr + other.negErr
return Err(value, posErr, negErr)
def __sub__(self, other):
value = self.value - other.value
posErr = self.posErr + other.negErr
negErr = self.negErr + other.posErr
return Err(value, posErr, negErr)
def __rsub__(self, other):
value = other.value - self.value
posErr = other.posErr + self.negErr
negErr = other.negErr + self.posErr
return Err(value, posErr, negErr)
def __mul__(self, other):
value = self.value * other.value
posErr = value * ((self.posErr/self.value) + (other.posErr/other.value))
negErr = value * ((self.negErr/self.value) + (other.negErr/other.value))
return Err(value, posErr, negErr)
def __rmul__(self, other):
value = self.value * other.value
posErr = value * ((self.posErr/self.value) + (other.posErr/other.value))
negErr = value * ((self.negErr/self.value) + (other.negErr/other.value))
return Err(value, posErr, negErr)
def __div__(self, other):
value = self.value / other.value
posErr = value * ((self.posErr/self.value) + (other.posErr/other.value))
negErr = value * ((self.negErr/self.value) + (other.negErr/other.value))
return Err(value, posErr, negErr)
def __rdiv__(self, other):
value = other.value / self.value
posErr = value * ((self.posErr/self.value) + (other.posErr/other.value))
negErr = value * ((self.negErr/self.value) + (other.negErr/other.value))
return Err(value, posErr, negErr)
def __pow__(self, y):
value = self.value**y.value
posErr = (self.value + self.posErr)**(y.value + y.posErr) - value
negErr = value - (self.value - self.negErr)**(y.value - y.negErr)
return Err(value, posErr, negErr)
def __rpow__(self, y):
value = y.value**self.value
posErr = (y.value + y.posErr)**(self.value + self.posErr) - value
negErr = value - (y.value - y.negErr)**(self.value - self.negErr)
return Err(value, posErr, negErr)
def __coerce__(self, other):
if isinstance(other, Err):
return self, other
try:
return self, Err(float(other))
except ValueError:
pass
def ApplyFunc(self, func):
value = func(self.value)
posErr = func(self.value + self.posErr) - value
negErr = value - func(self.value - self.negErr)
return Err(value, posErr, negErr)
def ArrayOfErr(errArray, posErr=0., negErr=None):
'''
Builds a Numeric array of Err() objects given an array and
a +ve or both a +ve and -ve errorval
Example:
egErr = ArrayOfErr([909., 802., 677., 585., 560., 548.], 1.0)
eg2Err = ArrayOfErr([909., 802., 677., 585., 560., 548.], 1.0, 1.5)
'''
valArray = asarray(errArray)
arrayLength = valArray.shape[0]
errs = repeat(array([Err(0.0, posErr, negErr)], PyObject), arrayLength)
return valArray + errs
def PrimeVals(errArray):
'''
Extract a Numeric array with just the prime value of the Err() objects
'''
li = []
for i in errArray:
li.append(i.value)
return asarray(li)
def PosErrs(errArray):
'''
Return a Numeric array with the posErr values of the Err() objects
'''
li = []
for i in errArray:
li.append(i.posErr)
return asarray(li)
def MaxVals(errArray):
'''
Return a Numeric array with the prime+posErr values of the Err() objects
'''
li = []
for i in errArray:
li.append(i.value + i.posErr)
return asarray(li)
def NegErrs(errArray):
'''
Return a Numeric array with the negErr values of the Err() objects
'''
li = []
for i in errArray:
li.append(i.negErr)
return asarray(li)
def MinVals(errArray):
'''
Return a Numeric array with the prime-negArr values of the Err() objects
'''
li = []
for i in errArray:
li.append(i.value - i.negErr)
return asarray(li)
def ErrToFloatArrays(errArray):
'''
Return a Numeric array with 3 rows:
row 1: the prime value of the Err() objects
row 2: the prime+posErr values of the Err() objects
row 3: the prime-negErr values of the Err() objects
'''
r1 = PrimeVals(errArray)
r2 = MaxVals(errArray)
r3 = MinVals(errArray)
return array([r1, r2, r3])
def FloatToErrArrays(floatArray):
'''
Convert a Numeric array of the form produced by ErrToFloatArrays, which contains 3 rows,
to an errArray
'''
r1 = floatArray[0]
r2 = floatArray[1]-r1
r3 = r1-floatArray[2]
return asarray(map(Err,r1,r2,r3), PyObject)
def ApplyUfuncToErr(errArray, func):
'''
Example:
a = ArrayOfErr([99., 82., 67., 58., 50., 54.], 1.0)
b = ApplyUfuncToErr(a, log) # This means log(a)
'''
a = ErrToFloatArrays(errArray)
b = func(a)
c = FloatToErrArrays(b)
return c
--
--
___________________________________________________________
Sign-up for Ads Free at Mail.com
http://promo.mail.com/adsfreejump.htm
More information about the SciPy-Dev
mailing list