Page MenuHomec4science
Files F62921719

measure.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 16, 13:16

measure.py
View Options

import math
import numpy
import statistics
import scipy
import scipy.optimize
def error_variance(*errs):
if len(errs) == 1:
return math.fabs(*errs)
out = 0
for err in errs:
out += err**2
return math.sqrt(out)
def error_abs(*errs):
out = 0
for err in errs:
out += math.fabs(err)
return out
class Measure:
def __init__(self, v, e = 0.0, sig_digit = 3, **args):
self.v = float(v)
self.e = math.fabs(float(e))
self.sig_digit = sig_digit
self.scient = False
# format
if "format" in args:
self.format = args["format"] # !!! no control
elif "scient" in args and \
isinstance(args["scient"], bool) and \
args["scient"]:
self.scient = True
self.format = "%%.%d%s" % (sig_digit, 'e')
else:
self.format = "%%.%d%s" % (sig_digit, 'g')
# std deviation, !! no control
self.errop = error_variance
if "stddev" in args:
self.errop = args["stddev"]
def __str__(self):
return ("%s +- %s" % (self.format, self.format)) % (self.v, self.e)
def __tuple__(self):
return (self.v, self.e)
def __iadd__(self, m):
if isinstance(m, Measure):
self.v += m.v
self.e = self.errop(self.e, m.e)
else:
self.v += m
return self
def __isub__(self, m):
if isinstance(m, Measure):
self.v -= m.v
self.e = self.errop(self.e, m.e)
else:
self.v -= m
return self
def __imul__(self, m):
if isinstance(m, Measure):
self.e = self.errop(self.v * m.e, self.e * m.v)
self.v *= m.v
else:
self.v *= m
self.e = self.errop(self.e * m)
return self
def __itruediv__(self, m):
if isinstance(m, Measure):
self.e = self.errop(self.e / m.v, self.v * m.e / m.v**2)
self.v /= m.v
else:
self.v /= m
self.e = self.errop(self.e / m)
return self
def __radd__(self, k):
return Measure(self.v + k, self.e, self.sig_digit, format = self.format)
def __rsub__(self, k):
return Measure(k - self.v, self.e, self.sig_digit, format = self.format)
def __rmul__(self, k):
return Measure(self.v * k, self.errop(self.e * k), self.sig_digit, format = self.format)
def __rtruediv__(self, k):
return Measure(k / self.v, self.errop(k * self.e / self.v**2), self.sig_digit, format = self.format)
def __add__(self, m):
if not isinstance(m, Measure):
return self.__radd__(m)
return Measure(self.v + m.v, self.errop(self.e, m.e), min(self.sig_digit, m.sig_digit), scient = self.scient)
def __sub__(self, m):
if not isinstance(m, Measure):
return Measure(self.v - m, self.e, self.sig_digit, scient = self.scient)
return Measure(self.v - m.v, self.errop(self.e, m.e), min(self.sig_digit, m.sig_digit), scient = self.scient)
def __mul__(self, m):
if not isinstance(m, Measure):
return self.__rmul__(m)
return Measure(self.v * m.v, self.errop(self.v * m.e, self.e * m.v), min(self.sig_digit, m.sig_digit), scient = self.scient)
def __truediv__(self, m):
if not isinstance(m, Measure):
return Measure(self.v / m, self.errop(self.e / m), self.sig_digit, scient = self.scient)
return Measure(self.v / m.v, self.errop(self.e / m.v, self.v * m.e / m.v**2), min(self.sig_digit, m.sig_digit), scient = self.scient)
def __pow__(self, r):
if not isinstance(r, int) and not isinstance(r, float):
raise TypeError("Measure: pow is only for int or float args")
err = r * self.v**(r-1) * self.e
return Measure(self.v**r, self.errop(err), self.sig_digit, scient = self.scient)
def __neg__(self):
return Measure(-self.v, self.e, self.sig_digit, scient = self.scient)
def __eq__(self, m):
if isinstance(m, Measure):
return m.v == self.v and m.e == self.e
return self.v == m
def __ne__(self, m):
if isinstance(m, Measure):
return m.v != self.v or m.e != self.e
return self.v != m
def __lt__(self, m):
if isinstance(m, Measure):
return self.v < m.v
return self.v < m
def __le__(self, m):
if isinstance(m, Measure):
return self.v <= m.v
return self.v <= m
def __gt__(self, m):
if isinstance(m, Measure):
return self.v > m.v
return self.v > m
def __ge__(self, m):
if isinstance(m, Measure):
return self.v >= m.v
return self.v >= m
def square(self):
return Measure(self.v**2, self.errop(2 * self.v * self.e), self.sig_digit, scient = self.scient)
def sqrt(self):
return Measure(math.sqrt(self.v), self.e / (2 * math.sqrt(self.v)), self.sig_digit, scient = self.scient)
def abs(self):
return Measure(math.fabs(self.v), self.e, self.sig_digit, scient = self.scient)
def exp(self):
ex = math.exp(self.v)
return Measure(ex, self.errop(ex * self.e), self.sig_digit, scient = self.scient)
def log(self):
return Measure(math.log(self.v), self.errop(self.e / self.v), self.sig_digit, scient = self.scient)
def sin(self):
return Measure(math.sin(self.v), self.errop(math.cos(self.v) * self.e), self.sig_digit, scient = self.scient)
def cos(self):
return Measure(math.cos(self.v), self.errop(math.sin(self.v) * self.e), self.sig_digit, scient = self.scient)
def __float__(self):
return self.v
def extr(self, coeff = 1.0):
return (self.v - self.e * coeff, self.v + self.e * coeff)
#
# Utils
#
# print a line into stream
def print_line(stream, *measures, **kargs):
M = []
for m in measures:
M.extend((m.format % m.v, m.format % m.e))
print(*M, file=stream, sep=" ", end="\n")
# give list of proper max and min, see Measure.extr
def extremums(measures, coeff = 1.0):
down = []
up = []
for m in measures:
e = m.extr(coeff)
down.append(e[0])
up.append(e[1])
return (tuple(up), tuple(down))
def values(measures):
v = []
for m in measures:
v.append(m.v)
return tuple(v)
# generate a list of measures
def genMeasures(values, error, sig_digit = 3):
out = []
for value in values:
out.append(Measure(value, error, sig_digit))
return out
# generate a measure basing on min and max
def minMaxMeasure(m, M):
return Measure((float(M) + float(m))/2, (float(M) - float(m))/2)
# return the maximal error
def maxError(measures):
M = 0
for m in measures:
if (m.e > M):
M = m.e
return M
# TODO, don't remember
def computeRange(y, Y, x, X):
dy = Y - y
dx = (X - x).v
rm = (dy.v-dy.e)/dx
rp = (dy.v+dy.e)/dx
return ( (rm, Y.v - Y.e - rm * x.v), (rp, Y.v + Y.e - rp * x.v) )
#
# analysis instruments
#
def polyfit(X, Y, deg = 1, **args):
coeff = 1
if "coeff" in args:
coeff = args["coeff"]
del args["coeff"]
sig_digit = 8
if "sig_digit" in args and isinstance(args["sig_digit"], int):
sig_digit = args["sig_digit"]
del args["sig_digit"]
extX = extremums(X, coeff)
extY = extremums(Y, coeff)
res = [None] * 4
# x down, y down
res[0] = numpy.polyfit(extX[0], extY[0], deg)
# x down, y up
res[1] = numpy.polyfit(extX[0], extY[1], deg)
# x up, y down
res[2] = numpy.polyfit(extX[1], extY[0], deg)
# x up, y up
res[3] = numpy.polyfit(extX[1], extY[1], deg)
(centered, cov) = (None, None)
try:
(centered, cov) = numpy.polyfit(values(X), values(Y), deg, None, False, None, True)
except:
centered = numpy.polyfit(values(X), values(Y), deg)
# compute the nearest coefficients
coeff = []
for i in range(deg + 1):
sample = (res[0][i], res[1][i], res[2][i], res[3][i])
stddev = statistics.variance(sample)
if cov is not None:
stddev += cov[i, i]
coeff.append(Measure( \
centered[i], \
math.sqrt(stddev), \
sig_digit \
))
return coeff
# fit a set of data
# basing on a general expression
# expr = expr(x, ...)
# where ... stands for the set of coefficients to deduce
def fit(X, Y, expr, p0 = None, **args):
coeff = 1
if "coeff" in args:
coeff = args["coeff"]
del args["coeff"]
sig_digit = 8
if "sig_digit" in args and isinstance(args["sig_digit"], int):
sig_digit = args["sig_digit"]
del args["sig_digit"]
extX = extremums(X, coeff)
extY = extremums(Y, coeff)
res = [None] * 4
# x down, y down
res[0] = scipy.optimize.curve_fit(expr, extX[0], extY[0], p0, **args)[0]
# x down, y up
res[1] = scipy.optimize.curve_fit(expr, extX[0], extY[1], p0, **args)[0]
# x up, y down
res[2] = scipy.optimize.curve_fit(expr, extX[1], extY[0], p0, **args)[0]
# x up, y up
res[3] = scipy.optimize.curve_fit(expr, extX[1], extY[1], p0, **args)[0]
# centered solution
centered, cov = scipy.optimize.curve_fit(expr, values(X), values(Y), p0, **args)
# compute the nearest coefficients
coeff = []
for i in range(len(res[0])):
sample = (res[0][i], res[1][i], res[2][i], res[3][i], centered[i])
coeff.append(Measure( \
centered[i], \
math.sqrt(statistics.variance(sample) + cov[i,i]), \
sig_digit \
))
return coeff
"""
def fit(X, Y, expr, p0 = None):
_X = values(X)
_Y = values(Y)
(res, pcov) = scipy.optimize.curve_fit(expr, _X, _Y, p0)
var = maxError(X)**2 + + maxError(Y)**2
out = []
for i in range(len(res)):
out.append(Measure(res[i], \
#math.sqrt(pcov[i][i] + var), \
math.sqrt(var), \
8))
return out
"""
def rawfit(X, Y, expr, p0 = None):
res = scipy.optimize.curve_fit(expr, X, Y, p0)[0]
return genMeasures(res, 0, 8)
#
# file loading instruments
#
def splitline(line):
out = (line[:-1]).replace('\t', ' ').split(" ")
flag = False
while(not flag):
i = 0
flag = False
while(i < len(out)):
if (out[i] == ''):
del out[i]
flag = True
else:
i += 1
return out
def loadMeasures(filename, values, errors = 0.0, gain = 1.0, offset = 0.0, verbose = False, **args):
f = open(filename, 'r')
sig_digit = 3
if "sig_digit" in args and isinstance(args["sig_digit"], int):
sig_digit = int(args["sig_digit"])
del args["sig_digit"]
firstline = None
if isinstance(values, str):
firstline = splitline(f.readline())
values = firstline.index(values)
errStr = False
if isinstance(errors, str):
if firstline is None:
firstline = splitline(f.readline())
errors = firstline.index(errors)
errStr = True
gainStr = False
if isinstance(gain, str):
if firstline is None:
firstline = splitline(f.readline())
gain = firstline.index(gain)
gainStr = True
offStr = False
if isinstance(offset, str):
if firstline is None:
firstline = splitline(f.readline())
offset = firstline.index(offset)
offStr = True
out = []
if verbose:
print("File %s: loading columns (%d, %d)" % (filename, values, errors))
for line in f:
l = splitline(line)
if line.startswith('#'):
continue # ignore
try:
# evaluate error
err = 0
if (errStr):
err = float(l[errors])
else:
err = errors
# evaluate gain scaling
g = 1.0
if (gainStr):
g = float(l[gain])
else:
g = gain
# evaluate offset
off = 1.0
if (offStr):
off = float(l[offset])
else:
off = offset
if values < len(l):
out.append(Measure(l[values], err, sig_digit, **args) * g + off)
except ValueError:
print("Could not parse line", line, "Ignoring it")
f.close()
return out

Event Timeline

c4science · Help