1############################################################################ 2# # 3# File : scidavisUtil.py # 4# Project : SciDAVis # 5# Description : Some Python utilities useful when doing data # 6# analysis with SciDAVis # 7# Copyright : (C) 2006 Knut Franke (knut.franke*gmx.de) # 8# (replace * with @ in the email address) # 9# # 10############################################################################ 11# # 12# This program is free software; you can redistribute it and/or modify # 13# it under the terms of the GNU General Public License as published by # 14# the Free Software Foundation; either version 2 of the License, or # 15# (at your option) any later version. # 16# # 17# This program is distributed in the hope that it will be useful, # 18# but WITHOUT ANY WARRANTY; without even the implied warranty of # 19# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 20# GNU General Public License for more details. # 21# # 22# You should have received a copy of the GNU General Public License # 23# along with this program; if not, write to the Free Software # 24# Foundation, Inc., 51 Franklin Street, Fifth Floor, # 25# Boston, MA 02110-1301 USA # 26# # 27############################################################################ 28 29""" 30 Some classes and functions that are particularly useful when doing calculations in SciDAVis; 31 maybe also outside of it. 32""" 33 34import math, scidavis 35 36class value: 37 """ 38 Scalar variable with associated error estimate. 39 40 This class provides an easy way to do calculations with values subject to (experimental) errors. 41 It is particularly useful for getting a quick idea of how where the main sources of uncertainty 42 are in your calculations and how good the result is. 43 44 However, if you need to do accurate hypthesis tests or have a precise range of possible values 45 at a given clearance level, e.g. for a publication, you will probably want to do an analytical 46 analysis of your final calculation and its error propagation. At the very least, you should be 47 aware of the following shortcomings of this class: 48 49 - it does not handle covariances 50 - it does not handle asymmetric error ranges 51 - it depends exclusively on linearized operations; 52 large errors will be propagated incorrectly in divisions, powers and roots 53 54 The standard constructor is 55 value(val, sigma) 56 where val is the actual value and sigma is the error estimate (standard deviation). You can also 57 construct a value from a SciDAVis table, like this: 58 value(table, valcol, errcol, row) 59 where table is a scidavis.Table instance, valcol the column containing the actual value, errcol the column 60 containing the error estimate and row, if you haven't guessed this, the row from which to read the values. 61 """ 62 def __init__(self,val,var=None,sigma=True,row=None): 63 if val.__class__==scidavis.Table: 64 self.val,self.var = val.cell(var,row),val.cell(sigma,row)**2 65 elif hasattr(val, "inherits") and val.inherits("Fit"): 66 self.val = val.results()[var-1] 67 covy = val.covarianceMatrix("tmp9999") 68 self.var = covy.cell(var, var) 69 covy.close() 70 elif var == None: 71 self.val,self.var = val,val 72 elif sigma: 73 self.val,self.var = val,var**2 74 else: 75 self.val,self.var = val,abs(var) 76 77 # numeric protocol 78 def __abs__(self): return value(abs(self.val), self.var, False) 79 def __neg__(self): return value(-self.val, self.var, False) 80 def __add__(self,o): 81 if o.__class__==value: 82 return value(self.val+o.val, self.var+o.var, False) 83 else: 84 return value(self.val+float(o), self.var, False) 85 def __radd__(self,o): return self.__add__(o) 86 def __sub__(self,o): 87 if o.__class__==value: 88 return value(self.val-o.val,self.var+o.var, False) 89 else: 90 return value(self.val-float(o),self.var, False) 91 def __rsub__(self,o): return -self.__sub__(o) 92 def __mul__(self,o): 93 if o.__class__==value: 94 return value(self.val*o.val,self.var*o.val**2+o.var*self.val**2, False) 95 else: 96 return value(self.val*float(o),self.var*float(o)**2, False) 97 def __rmul__(self,o): return self.__mul__(o) 98 def __div__(self,o): 99 if o.__class__==value: 100 return value(self.val/o.val, self.var/o.val**2+o.var*self.val**2/o.val**4, False) 101 else: 102 return value(self.val/float(o), self.var/float(o)**2, False) 103 def __rdiv__(self,o): 104 return value(float(o)/self.val, self.var*float(o)**2/self.val**4, False) 105 def __pow__(self,o): 106 if o.__class__==value: 107 return value(self.val**o.val, self.var*o.val**2*self.val**(2*(o.val-1))+o.var*math.ln(self.val)*self.val**o.val, False) 108 else: 109 return value(self.val**float(o), self.var*float(o)**2*self.val**(2*(float(o)-1)), False) 110 def __rpow__(self,o): 111 return value(float(o)**self.val, self.var*math.ln(o)*float(o)**self.val, False) 112 113 # comparison protocol 114 def __lt__(self,o): 115 if o.__class__==value: return self.val < o.val 116 else: return self.val < float(o) 117 def __le__(self,o): 118 if o.__class__==value: return self.val <= o.val 119 else: return self.val <= float(o) 120 def __eq__(self,o): 121 if o.__class__==value: return self.val == o.val 122 else: return self.val == float(o) 123 def __ge__(self,o): 124 if o.__class__==value: return self.val >= o.val 125 else: return self.val >= float(o) 126 def __gt__(self,o): 127 if o.__class__==value: return self.val > o.val 128 else: return self.val > float(o) 129 130 # convenience functions 131 def sigma(self): 132 "compute the error estimate (standard deviation)" 133 return math.sqrt(self.var) 134 def __str__(self): 135 return str(self.val) + " +/- " + str(self.sigma()) 136 def sqrt(self): 137 "error estimate aware square root function" 138 if self.val==0: return value(0, self.var, False) 139 return value(math.sqrt(self.val), self.var/(4*self.val), False) 140 141def cartesic(r, phi): 142 "Convert polar coordinates to cartesic ones. Arguments have to instances of class value." 143 x = r*value(math.cos(phi.val), (math.sin(phi.val))**2*phi.var, False) 144 y = r*value(math.sin(phi.val), (math.cos(phi.val))**2*phi.var, False) 145 return (x, y) 146 147def polar(x, y): 148 "Convert (2D) cartesic coordinates to polar ones. Arguments have to instances of class value." 149 r = (x**2+y**2).sqrt() 150 yx = y/x 151 phi = value(math.atan2(y.val, x.val), yx.var/(1+yx.val**2)**2, False) 152 return (r, phi) 153 154def lookup(tab, col, txt): 155 """ 156 lookup(tab, col, txt): 157 search column col of SciDAVis table tab for a cell with text content txt, 158 returning the row of the first match; None on failure. 159 """ 160 for i in range(1,tab.numRows()+1): 161 if tab.text(col, i) == txt: return i 162 163class vec4: 164 """ 165 A relativistic 4-vector. 166 Computations assume a diag(1, -1, -1, -1) metric. 167 Components may be any numbers that can either be cast to float or provide a sqrt() method. 168 """ 169 def __init__(self, one, x=None, y=0, z=0): 170 self.t, self.x, self.y, self.z = 0,0,0,0 171 try: self.t, self.x, self.y, self.z = one[0], one[1], one[2], one[3] 172 except IndexError: pass 173 except (AttributeError, TypeError): 174 try: self.t, self.x, self.y, self.z = one,x[0],x[1],x[2] 175 except IndexError: self.t, self.x, self.y, self.z = one,x[0],x[1],0 176 except (AttributeError,TypeError): self.t, self.x, self.y, self.z = one,x,y,z 177 def __add__(self, o): 178 return vec4(self.t+o.t, self.x+o.x, self.y+o.y, self.z+o.z) 179 def __sub__(self, o): 180 return vec4(self.t-o.t, self.x-o.x, self.y-o.y, self.z-o.z) 181 def __mul__(self, o): 182 return self.t*o.t-self.x*o.x-self.y*o.y-self.z*o.z 183 def __str__(self): 184 return "(%s; %s, %s, %s)" % (self.t, self.x, self.y, self.z) 185 def __abs__(self): 186 s = self*self 187 if s < 0: 188 if (hasattr(s, "sqrt")): return -(-s).sqrt() 189 else: return -math.sqrt(-s) 190 else: 191 if hasattr(s, "sqrt"): return s.sqrt() 192 else: return math.sqrt(s) 193 def abs3(self): 194 "compute the absolute value of the spatial part of this 4-vector" 195 s = self.x**2+self.y**2+self.z**2 196 if hasattr(s, "sqrt"): return s.sqrt() 197 else: return math.sqrt(s) 198 def gamma(self): 199 "assuming this is an energy-momentum four-vector, return the corresponding gamma factor" 200 return self.t/abs(self) 201 def beta(self): 202 "assuming this is an energy-momentum four-vector, return the corresponding beta factor" 203 s = 1-(abs(self)/self.t)**2 204 if hasattr(s, "sqrt"): return s.sqrt() 205 else: return math.sqrt(s) 206 def betagamma(self): 207 "compute self.beta()*self.gamma() more efficiently" 208 s = (self.t/abs(self))**2-1 209 if hasattr(s, "sqrt"): return s.sqrt() 210 else: return math.sqrt(s) 211 212def exportTableToTeX(t, filename=None): 213 """ 214 exportToTeX(table, filename=None): 215 Export table as TeX-tabular to filename. If filename==None, popup a file selection dialog. 216 """ 217 if (scidavis.app.qtVersion() >= 0x050000): 218 from PyQt5.QtWidgets import QFileDialog 219 else: 220 from PyQt4.QtGui import QFileDialog 221 if not filename: 222 filename=QFileDialog.getSaveFileName(scidavis.app,"SciDAVis - Export TeX table","","All files *;;TeX documents (*.tex *.TEX);;"); 223 f=open(filename,'w') 224 f.write('\\begin{tabular}{|' + 'c|'*t.numCols() + '}\\hline\n') 225 for col in range(1,t.numCols()): 226 f.write(t.colName(col) + ' & ') 227 f.write(t.colName(t.numCols()) + ' \\\\\\hline\n') 228 for row in range(1,t.numRows()+1): 229 val = False 230 for col in range(1,t.numCols()+1): 231 if t.text(col,row) != "": val = True 232 if val: 233 for col in range(1,t.numCols()): 234 f.write(t.text(col,row) + ' & ') 235 f.write(t.text(t.numCols(),row) + ' \\\\\n') 236 f.write('\\hline\n\\end{tabular}\n') 237 f.close() 238 239def exportAllTablesToTeX(folder, dir, recurs=True): 240 for w in folder.windows(): 241 if w.isA("Table"): exportTableToTeX(w, dir+w.name()+".tex") 242 if recurs: 243 for f in folder.folders(): 244 exportAllTablesToTeX(f, dir, True) 245 246