0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022"""bring mathematical functions from the Numeric and cvisual modles into the PyGeo namespace,
0023and create convenience functions and classes for analytical operations"""
0024
0025import Numeric as N
0026import cmath
0027
0028import LinearAlgebra
0029from pygeo.base.support.pygeoopts import EPS
0030
0031
0032try:
0033 import cvisual
0034except ImportError:
0035 import visual.cvisual as cvisual
0036
0037PyGeoMath=['PI','math_E','absolute','add','subtract','multiply','divide', 'sqrt',
0038'arctan', 'arctan2', 'cos','sin','arccos','arcsin','tan','identity','hypot',
0039'matrixmultiply','transpose', 'arrayrange','array','dot', 'reshape', 'ones','conjugate',
0040'trace', 'asarray','shape', 'determinant','solve','SVD','inverse', 'eigenvectors',
0041'AlgebraError','cross','vector','toXY','cross3','quadratic','homogenous','cross_ratio',
0042 'mod','mod2','Hermitian','Mobius']
0043
0044__all__=PyGeoMath
0045
0046ArrayType=type(N.array([]))
0047
0048
0049PI=N.pi
0050math_E=N.e
0051absolute=N.absolute
0052add=N.add
0053subtract=N.subtract
0054multiply=N.multiply
0055divide=N.divide
0056sqrt=N.sqrt
0057arctan=N.arctan
0058arctan2=N.arctan2
0059cos=N.cos
0060sin=N.sin
0061arccos=N.arccos
0062arcsin=N.arcsin
0063tan=N.tan
0064identity=N.identity
0065hypot=N.hypot
0066matrixmultiply=N.matrixmultiply
0067transpose=N.transpose
0068arrayrange=N.arrayrange
0069array=N.array
0070dot=N.dot
0071reshape=N.reshape
0072ones=N.ones
0073conjugate=N.conjugate
0074trace=N.trace
0075asarray=N.asarray
0076shape=N.shape
0077determinant=LinearAlgebra.determinant
0078solve=LinearAlgebra.solve_linear_equations
0079SVD= LinearAlgebra.singular_value_decomposition
0080inverse= LinearAlgebra.inverse
0081eigenvectors=LinearAlgebra.eigenvectors
0082AlgebraError = LinearAlgebra.LinAlgError
0083vector=cvisual.vector
0084cross=cvisual.cross
0085
0086def toXY(point):
0087 """ return a 3 element array representing the origin centered projection
0088 of the given point to the plane parallel with the XY plane at unit distance"""
0089
0090 __mat2=array([[-1.,0.,0.,0.],[0.,-1.,0.,0.],[0.,0.,0.,-1.],[0.,0.,0.,-1]])
0091 tf= matrixmultiply(homogenous(point),__mat2)
0092 return array([tf[0]/tf[3],tf[1]/tf[3],1])
0093
0094def quadratic(a,b,c):
0095 """the return the roots of the quadratic equation with
0096 the three given values as coefficients - i.e. a*x**2 + b&x +c """
0097
0098 try:
0099 x1=sqrt(b**2-4*a*c)
0100 try:
0101 h1=(-b-x1)/(2*a)
0102 h2=(-b+x1)/(2*a)
0103 except ZeroDivisionError:
0104 h1=None
0105 h2=None
0106 except ValueError:
0107 h1=None
0108 h2=None
0109 return h1,h2
0110
0111def cross3(v1,v2,v3):
0112 """return the crossproduct of the vectors connecting the 1st given
0113 position with that of the 2nd and 3rd given positions"""
0114
0115 return cross(v2-v1,v3-v1)
0116
0117def homogenous(narray):
0118 """return a 4 element array as a homogenous coordinate of the given
0119 3 element vector"""
0120
0121 return array((narray[0],narray[1],narray[2],1.),'d')
0122
0123def cross_ratio(a,b,c,d):
0124 """ return the crossratio of 4 ordered points"""
0125
0126 ab=a.distance(b)
0127 bc=b.distance(c)
0128 ad=a.distance(d)
0129 dc=-c.distance(d)
0130 return (ab/bc)/(ad/dc)
0131
0132def mod(c):
0133 """return the modulus (absolute value) of the given complex number"""
0134
0135 return absolute(c)
0136
0137def mod2(c):
0138 """return the square of the modulus of the given complex number"""
0139
0140 return absolute(c)**2
0141
0142
0143
0144
0145class base2by2(object):
0146 """the Numeric Usearray coce, but derived from object, so that inherited classes
0147 can use properties"""
0148
0149 def __init__(self):
0150 self.array=N.zeros((2,2),'D')
0151 self.shape = (2,2)
0152 self._typecode = self.array.typecode()
0153 self.name = ''.split(str(self.__class__))[0]
0154
0155 def __repr__(self):
0156 if len(self.array.shape) > 0:
0157 return self.__class__.__name__+repr(self.array)[len("array"):]
0158 else:
0159 return self.__class__.__name__+"("+repr(self.array)+")"
0160
0161 def __array__(self,t=None):
0162 if t: return asarray(self.array,t)
0163 return asarray(self.array)
0164
0165 def __float__(self):
0166 return float(asarray(self.array))
0167
0168
0169 def __len__(self): return len(self.array)
0170
0171 def __getitem__(self, index):
0172 return self.array[index]
0173
0174 def __getslice__(self, i, j):
0175 return self.array[i:j]
0176
0177
0178 def __setitem__(self, index, value): self.array[index] = asarray(value,self._typecode)
0179 def __setslice__(self, i, j, value): self.array[i:j] = asarray(value,self._typecode)
0180
0181
0182
0183
0184
0185
0186 def __abs__(self): return self._rc(absolute(self.array))
0187 def __neg__(self): return self._rc(-self.array)
0188
0189 def __add__(self, other):
0190 return self._rc(self.array+asarray(other))
0191 __radd__ = __add__
0192
0193 def __iadd__(self, other):
0194 add(self.array, other, self.array)
0195 return self
0196
0197 def __sub__(self, other):
0198 return self._rc(self.array-asarray(other))
0199 def __rsub__(self, other):
0200 return self._rc(asarray(other)-self.array)
0201 def __isub__(self, other):
0202 subtract(self.array, other, self.array)
0203 return self
0204
0205 def __mul__(self, other):
0206 return self._rc(multiply(self.array,asarray(other)))
0207 __rmul__ = __mul__
0208 def __imul__(self, other):
0209 multiply(self.array, other, self.array)
0210 return self
0211
0212 def __div__(self, other):
0213 return self._rc(divide(self.array,asarray(other)))
0214 def __rdiv__(self, other):
0215 return self._rc(divide(asarray(other),self.array))
0216 def __idiv__(self, other):
0217 divide(self.array, other, self.array)
0218 return self
0219
0220 def __mod__(self, other):
0221 return self._rc(remainder(self.array, other))
0222 def __rmod__(self, other):
0223 return self._rc(remainder(other, self.array))
0224 def __imod__(self, other):
0225 remainder(self.array, other, self.array)
0226 return self
0227
0228 def __divmod__(self, other):
0229 return (self._rc(divide(self.array,other)),
0230 self._rc(remainder(self.array, other)))
0231 def __rdivmod__(self, other):
0232 return (self._rc(divide(other, self.array)),
0233 self._rc(remainder(other, self.array)))
0234
0235 def __pow__(self,other):
0236 return self._rc(power(self.array,asarray(other)))
0237 def __rpow__(self,other):
0238 return self._rc(power(asarray(other),self.array))
0239 def __ipow__(self,other):
0240 power(self.array, other, self.array)
0241 return self
0242
0243 def __lshift__(self,other):
0244 return self._rc(left_shift(self.array, other))
0245 def __rshift__(self,other):
0246 return self._rc(right_shift(self.array, other))
0247 def __rlshift__(self,other):
0248 return self._rc(left_shift(other, self.array))
0249 def __rrshift__(self,other):
0250 return self._rc(right_shift(other, self.array))
0251 def __ilshift__(self,other):
0252 left_shift(self.array, other, self.array)
0253 return self
0254 def __irshift__(self,other):
0255 right_shift(self.array, other, self.array)
0256 return self
0257
0258 def __and__(self, other):
0259 return self._rc(bitwise_and(self.array, other))
0260 def __rand__(self, other):
0261 return self._rc(bitwise_and(other, self.array))
0262 def __iand__(self, other):
0263 bitwise_and(self.array, other, self.array)
0264 return self
0265
0266 def __xor__(self, other):
0267 return self._rc(bitwise_xor(self.array, other))
0268 def __rxor__(self, other):
0269 return self._rc(bitwise_xor(other, self.array))
0270 def __ixor__(self, other):
0271 bitwise_xor(self.array, other, self.array)
0272 return self
0273
0274 def __or__(self, other):
0275 return self._rc(bitwise_or(self.array, other))
0276 def __ror__(self, other):
0277 return self._rc(bitwise_or(other, self.array))
0278 def __ior__(self, other):
0279 bitwise_or(self.array, other, self.array)
0280 return self
0281
0282 def __neg__(self):
0283 return self._rc(-self.array)
0284 def __pos__(self):
0285 return self._rc(self.array)
0286 def __abs__(self):
0287 return self._rc(abs(self.array))
0288 def __invert__(self):
0289 return self._rc(invert(self.array))
0290
0291 def _scalarfunc(a, func):
0292 if len(a.shape) == 0:
0293 return func(a[0])
0294 else:
0295 raise TypeError, "only rank-0 arrays can be converted to Python scalars."
0296
0297 def __complex__(self): return self._scalarfunc(complex)
0298 def __float__(self): return self._scalarfunc(float)
0299 def __int__(self): return self._scalarfunc(int)
0300 def __long__(self): return self._scalarfunc(long)
0301 def __hex__(self): return self._scalarfunc(hex)
0302 def __oct__(self): return self._scalarfunc(oct)
0303
0304 def __lt__(self,other): return self._rc(less(self.array,other))
0305 def __le__(self,other): return self._rc(less_equal(self.array,other))
0306 def __eq__(self,other): return self._rc(equal(self.array,other))
0307 def __ne__(self,other): return self._rc(not_equal(self.array,other))
0308 def __gt__(self,other): return self._rc(greater(self.array,other))
0309 def __ge__(self,other): return self._rc(greater_equal(self.array,other))
0310
0311 def copy(self): return self._rc(self.array.copy())
0312
0313 def tostring(self): return self.array.tostring()
0314
0315 def byteswapped(self): return self._rc(self.array.byteswapped())
0316
0317 def astype(self, typecode): return self._rc(self.array.astype(typecode))
0318
0319 def typecode(self): return self._typecode
0320
0321 def itemsize(self): return self.array.itemsize()
0322
0323 def iscontiguous(self): return self.array.iscontiguous()
0324
0325 def _rc(self, a):
0326
0327 return self.__class__(a)
0328
0329 def determinant(self):
0330 det = self.array[0][0]*self.array[1][1] - self.array[0][1]*self.array[1][0]
0331 return det
0332
0333 def normalize(self):
0334 """normalize to matrix of unit determinant"""
0335 sdet=sqrt(self.determinant())
0336 self.array/=sdet
0337 return self
0338
0339 def inverse(self):
0340 """ invert the matrix in place"""
0341
0342 return self._rc([[self.array[1][1],-self.array[0][1]],[-self.array[1][0],self.array[0][0]]])/self.determinant()
0343
0344 def transpose(self):
0345
0346 return self._rc([[self.array[0][0],self.array[1][0]],[self.array[0][1],self.array[1][1]]])
0347
0348 def conjugate(self):
0349 return self._rc(conjugate(self.array))
0350
0351class Hermitian(base2by2):
0352 """ the Hermitian matrix [[A,B],[C,D]] """
0353 def __init__(self,mat):
0354 base2by2.__init__(self)
0355 self.A=mat[0][0]
0356 self.B=mat[0][1]
0357 self.C=mat[1][0]
0358 self.D=mat[1][1]
0359
0360 def setA(self,value):
0361 try:
0362 if N.absolute(value.imag) > EPS:
0363 raise TypeError, "malformed Hermitian matrix"
0364 else:
0365 self.array[0][0] = complex(value.real,value.imag)
0366 except AttributeError:
0367 self.array[0][0] = complex(value)
0368
0369 def getA(self):
0370 return self.array[0][0].real
0371
0372 def setB(self,value):
0373 try:
0374 self.array[0][1] = complex(value.real,value.imag)
0375 except AttributeError:
0376 self.array[0][1] = complex(value)
0377
0378 def getB(self):
0379 return self.array[0][1]
0380
0381 def setC(self,value):
0382 try:
0383 self.array[1][0] = complex(value.real,value.imag)
0384 except AttributeError:
0385 self.array[1][0] = complex(value)
0386
0387 def getC(self):
0388 return self.array[1][0]
0389
0390 def setD(self,value):
0391 try:
0392 if N.absolute(value.imag) > EPS:
0393 print self, value.imag
0394 raise TypeError, "malformed Hermitian matrix"
0395 else:
0396 self.array[1][1] = complex(value.real,value.imag)
0397 except AttributeError:
0398 self.array[1][1] = complex(value)
0399
0400 def getD(self):
0401 return self.array[1][1].real
0402
0403 A = property(getA,setA)
0404 B = property(getB,setB)
0405 C = property(getC,setC)
0406 D = property(getD,setD)
0407
0408 def __repr__(self):
0409 return "Hermitian"+repr(self.array)[5:]
0410
0411 def transform(self,mat):
0412 t1=mat.inverse()
0413 return Hermitian(matrixmultiply(t1.transpose(),matrixmultiply(self.array,t1.conjugate())))
0414
0415 conjugate = transpose
0416
0417class Mobius(base2by2):
0418 """ the complex valued Mobius transformation matrix [[a,b],[c,d]]"""
0419
0420 def __init__(self,mat):
0421 base2by2.__init__(self)
0422 self.a=mat[0][0]
0423 self.b=mat[0][1]
0424 self.c=mat[1][0]
0425 self.d=mat[1][1]
0426
0427 def setA(self,value):
0428 try:
0429 self.array[0][0] = complex(value.real,value.imag)
0430 except AttributeError:
0431 self.array[0][0] = complex(value)
0432
0433 def getA(self):
0434 return self.array[0][0]
0435
0436 def setB(self,value):
0437 try:
0438 self.array[0][1] = complex(value.real,value.imag)
0439 except AttributeError:
0440 self.array[0][1] = complex(value)
0441
0442 def getB(self):
0443 return self.array[0][1]
0444
0445 def setC(self,value):
0446 try:
0447 self.array[1][0] = complex(value.real,value.imag)
0448 except AttributeError:
0449 self.array[1][0] = complex(value)
0450
0451 def getC(self):
0452 return self.array[1][0]
0453
0454 def setD(self,value):
0455 try:
0456 self.array[1][1] = complex(value.real,value.imag)
0457 except AttributeError:
0458 self.array[1][1] = complex(value)
0459
0460 def getD(self):
0461 return self.array[1][1]
0462
0463 a = property(getA,setA)
0464 b = property(getB,setB)
0465 c = property(getC,setC)
0466 d = property(getD,setD)
0467
0468 def normal_form(self):
0469 """convert in place to the normal form (fixed points at origin and infinity) of the
0470 transformation represented by the Mobius matrix"""
0471
0472 f1,f2=self.getFixed()
0473 try:
0474 mat=array([[1.,-f2],[1.,-f1]])
0475 self.conjugate(mat)
0476 except TypeError:
0477 print "fixed point at infinity"
0478
0479 def getFixed(self):
0480 """return the fixed points of the transformation represented by the
0481 Mobius matrix"""
0482
0483 if absolute(self.determinant() -1) > EPS:
0484 self.normalize()
0485 evalues, evectors = eigenvectors(self.array)
0486 return evectors
0487
0488 def getMultiplier(self):
0489 """return the multipliers of the transformation represented by the
0490 Mobius matrix and associated with the fixed points"""
0491
0492 if absolute(self.determinant()-1)> EPS:
0493 self.normalize()
0494 evalues, evectors = eigenvectors(self.array)
0495 return 1./(evalues[0]**2),1./(evalues[1]**2)
0496
0497 def getPole(self):
0498 """return the pole - i.e. the point transformed to the point at infinity - of the
0499 transformation represented by the Mobius matrix"""
0500
0501 return -self.d/self.c
0502
0503 def getInversePole(self):
0504 """return the inverse pole - i.e. the point transformed to the point at infinity - by
0505 the transformation by the invserse of Mobius matrix"""
0506
0507 return self.a/self.c
0508
0509 def getTrace(self):
0510 """return trace - the sum of the diagonal elements - of the Mboius matrix"""
0511 return self.a + self.d
0512
0513
0514
0515
0516
0517 def getConstant(self):
0518 """ return the characteristic constant associated with the Mobius matrix"""
0519 determinant=self.getDeterminant()
0520 if absolute(self.determinant()-1)> EPS:
0521 self.normalize()
0522 trace=self.getTrace()
0523 invariant=(self.a-self.d)**2 +4 * self.b*self.c
0524 return (trace+sqrt(invariant*determinant))/(trace-sqrt(invariant*determinant))
0525
0526 def getInvariant(self):
0527
0528 if absolute(self.determinant()-1)> EPS:
0529 self.normalize()
0530 return (self.a-self.d)**2 +4 * self.b*self.c
0531
0532
0533__author__ = "Arthur Siegel <ajsiegel@optonline.com>"
0534__date__ = "Date: 2006-02-02"
0535__revision__ = "Revision: a1"
0536__url__ = "http://source.net/projects/pygeo"
0537__license__ ="GPL <http://www.opensource.org/licenses/gpl-license.php>"