test_test.py
1 import sys 2 import pytest 3 from sympy import symbols, sin, cos, Rational, expand, collect, simplify, Symbol, Add, S, eye 4 from galgebra.printer import Format, Eprint, latex, GaPrinter 5 from galgebra.ga import Ga 6 from galgebra.mv import Mv, Nga 7 # for backward compatibility 8 from galgebra import ga, metric 9 10 one = S.One 11 12 13 def F(x): 14 global n, nbar 15 Fx = S.Half * ((x * x) * n + 2 * x - nbar) 16 return Fx 17 18 def make_vector(a, n=3, ga=None): 19 if isinstance(a, str): 20 v = S.Zero 21 for i in range(n): 22 a_i = Symbol(a+str(i+1)) 23 v += a_i*ga.basis[i] 24 v = ga.mv(v) 25 return F(v) 26 else: 27 return F(a) 28 29 class TestTest: 30 31 def test_basic_multivector_operations(self): 32 33 g3d = Ga('e*x|y|z') 34 ex, ey, ez = g3d.mv() 35 36 A = g3d.mv('A', 'mv') 37 38 assert str(A) == 'A + A__x*e_x + A__y*e_y + A__z*e_z + A__xy*e_x^e_y + A__xz*e_x^e_z + A__yz*e_y^e_z + A__xyz*e_x^e_y^e_z' 39 assert str(A) == 'A + A__x*e_x + A__y*e_y + A__z*e_z + A__xy*e_x^e_y + A__xz*e_x^e_z + A__yz*e_y^e_z + A__xyz*e_x^e_y^e_z' 40 assert str(A) == 'A + A__x*e_x + A__y*e_y + A__z*e_z + A__xy*e_x^e_y + A__xz*e_x^e_z + A__yz*e_y^e_z + A__xyz*e_x^e_y^e_z' 41 42 X = g3d.mv('X', 'vector') 43 Y = g3d.mv('Y', 'vector') 44 45 assert str(X) == 'X__x*e_x + X__y*e_y + X__z*e_z' 46 assert str(Y) == 'Y__x*e_x + Y__y*e_y + Y__z*e_z' 47 48 assert str((X*Y)) == '(e_x.e_x)*X__x*Y__x + (e_x.e_y)*X__x*Y__y + (e_x.e_y)*X__y*Y__x + (e_x.e_z)*X__x*Y__z + (e_x.e_z)*X__z*Y__x + (e_y.e_y)*X__y*Y__y + (e_y.e_z)*X__y*Y__z + (e_y.e_z)*X__z*Y__y + (e_z.e_z)*X__z*Y__z + (X__x*Y__y - X__y*Y__x)*e_x^e_y + (X__x*Y__z - X__z*Y__x)*e_x^e_z + (X__y*Y__z - X__z*Y__y)*e_y^e_z' 49 assert str((X^Y)) == '(X__x*Y__y - X__y*Y__x)*e_x^e_y + (X__x*Y__z - X__z*Y__x)*e_x^e_z + (X__y*Y__z - X__z*Y__y)*e_y^e_z' 50 assert str((X|Y)) == '(e_x.e_x)*X__x*Y__x + (e_x.e_y)*X__x*Y__y + (e_x.e_y)*X__y*Y__x + (e_x.e_z)*X__x*Y__z + (e_x.e_z)*X__z*Y__x + (e_y.e_y)*X__y*Y__y + (e_y.e_z)*X__y*Y__z + (e_y.e_z)*X__z*Y__y + (e_z.e_z)*X__z*Y__z' 51 52 53 g2d = Ga('e*x|y') 54 ex, ey = g2d.mv() 55 56 X = g2d.mv('X', 'vector') 57 A = g2d.mv('A', 'spinor') 58 59 assert str(X) == 'X__x*e_x + X__y*e_y' 60 assert str(A) == 'A + A__xy*e_x^e_y' 61 62 assert str((X|A)) == 'A__xy*(-(e_x.e_y)*X__x - (e_y.e_y)*X__y)*e_x + A__xy*((e_x.e_x)*X__x + (e_x.e_y)*X__y)*e_y' 63 assert str((X<A)) == 'A__xy*(-(e_x.e_y)*X__x - (e_y.e_y)*X__y)*e_x + A__xy*((e_x.e_x)*X__x + (e_x.e_y)*X__y)*e_y' 64 assert str((A>X)) == 'A__xy*((e_x.e_y)*X__x + (e_y.e_y)*X__y)*e_x + A__xy*(-(e_x.e_x)*X__x - (e_x.e_y)*X__y)*e_y' 65 66 67 o2d = Ga('e*x|y', g=[1, 1]) 68 ex, ey = o2d.mv() 69 70 X = o2d.mv('X', 'vector') 71 A = o2d.mv('A', 'spinor') 72 73 assert str(X) == 'X__x*e_x + X__y*e_y' 74 assert str(A) == 'A + A__xy*e_x^e_y' 75 76 assert str((X*A)) == '(A*X__x - A__xy*X__y)*e_x + (A*X__y + A__xy*X__x)*e_y' 77 assert str((X|A)) == '-A__xy*X__y*e_x + A__xy*X__x*e_y' 78 assert str((X<A)) == '-A__xy*X__y*e_x + A__xy*X__x*e_y' 79 assert str((X>A)) == 'A*X__x*e_x + A*X__y*e_y' 80 81 assert str((A*X)) == '(A*X__x + A__xy*X__y)*e_x + (A*X__y - A__xy*X__x)*e_y' 82 assert str((A|X)) == 'A__xy*X__y*e_x - A__xy*X__x*e_y' 83 assert str((A<X)) == 'A*X__x*e_x + A*X__y*e_y' 84 assert str((A>X)) == 'A__xy*X__y*e_x - A__xy*X__x*e_y' 85 86 def test_check_generalized_BAC_CAB_formulas(self): 87 88 a, b, c, d, e = Ga('a b c d e').mv() 89 90 assert str(a|(b*c)) == '-(a.c)*b + (a.b)*c' 91 assert str(a|(b^c)) == '-(a.c)*b + (a.b)*c' 92 assert str(a|(b^c^d)) == '(a.d)*b^c - (a.c)*b^d + (a.b)*c^d' 93 94 expr = (a|(b^c))+(c|(a^b))+(b|(c^a)) # = (a.b)*c - (b.c)*a - ((a.b)*c - (b.c)*a) 95 assert str(expr.simplify()) == '0' 96 97 assert str(a*(b^c)-b*(a^c)+c*(a^b)) == '3*a^b^c' 98 assert str(a*(b^c^d)-b*(a^c^d)+c*(a^b^d)-d*(a^b^c)) == '4*a^b^c^d' 99 assert str((a^b)|(c^d)) == '-(a.c)*(b.d) + (a.d)*(b.c)' 100 assert str(((a^b)|c)|d) == '-(a.c)*(b.d) + (a.d)*(b.c)' 101 assert str(Ga.com(a^b, c^d)) == '-(b.d)*a^c + (b.c)*a^d + (a.d)*b^c - (a.c)*b^d' 102 assert str((a|(b^c))|(d^e)) == '(-(a.b)*(c.e) + (a.c)*(b.e))*d + ((a.b)*(c.d) - (a.c)*(b.d))*e' 103 104 def test_derivatives_in_rectangular_coordinates(self): 105 106 X = x, y, z = symbols('x y z') 107 o3d = Ga('e_x e_y e_z', g=[1, 1, 1], coords=X) 108 ex, ey, ez = o3d.mv() 109 grad = o3d.grad 110 111 # like sympy.sstr but using our printer 112 sstr = lambda x: GaPrinter().doprint(x) 113 114 f = o3d.mv('f', 'scalar', f=True) 115 A = o3d.mv('A', 'vector', f=True) 116 B = o3d.mv('B', 'bivector', f=True) 117 C = o3d.mv('C', 'mv', f=True) 118 119 assert sstr(f) == 'f' 120 assert sstr(A) == 'A__x*e_x + A__y*e_y + A__z*e_z' 121 assert sstr(B) == 'B__xy*e_x^e_y + B__xz*e_x^e_z + B__yz*e_y^e_z' 122 assert sstr(C) == 'C + C__x*e_x + C__y*e_y + C__z*e_z + C__xy*e_x^e_y + C__xz*e_x^e_z + C__yz*e_y^e_z + C__xyz*e_x^e_y^e_z' 123 124 assert sstr(grad*f) == 'D{x}f*e_x + D{y}f*e_y + D{z}f*e_z' 125 assert sstr(grad|A) == 'D{x}A__x + D{y}A__y + D{z}A__z' 126 assert sstr(grad*A) == 'D{x}A__x + D{y}A__y + D{z}A__z + (-D{y}A__x + D{x}A__y)*e_x^e_y + (-D{z}A__x + D{x}A__z)*e_x^e_z + (-D{z}A__y + D{y}A__z)*e_y^e_z' 127 128 assert sstr(-o3d.I()*(grad^A)) == '(-D{z}A__y + D{y}A__z)*e_x + (D{z}A__x - D{x}A__z)*e_y + (-D{y}A__x + D{x}A__y)*e_z' 129 assert sstr(grad*B) == '(-D{y}B__xy - D{z}B__xz)*e_x + (D{x}B__xy - D{z}B__yz)*e_y + (D{x}B__xz + D{y}B__yz)*e_z + (D{z}B__xy - D{y}B__xz + D{x}B__yz)*e_x^e_y^e_z' 130 assert sstr(grad^B) == '(D{z}B__xy - D{y}B__xz + D{x}B__yz)*e_x^e_y^e_z' 131 assert sstr(grad|B) == '(-D{y}B__xy - D{z}B__xz)*e_x + (D{x}B__xy - D{z}B__yz)*e_y + (D{x}B__xz + D{y}B__yz)*e_z' 132 133 assert sstr(grad<A) == 'D{x}A__x + D{y}A__y + D{z}A__z' 134 assert sstr(grad>A) == 'D{x}A__x + D{y}A__y + D{z}A__z' 135 assert sstr(grad<B) == '(-D{y}B__xy - D{z}B__xz)*e_x + (D{x}B__xy - D{z}B__yz)*e_y + (D{x}B__xz + D{y}B__yz)*e_z' 136 assert sstr(grad>B) == '0' 137 assert sstr(grad<C) == 'D{x}C__x + D{y}C__y + D{z}C__z + (-D{y}C__xy - D{z}C__xz)*e_x + (D{x}C__xy - D{z}C__yz)*e_y + (D{x}C__xz + D{y}C__yz)*e_z + D{z}C__xyz*e_x^e_y - D{y}C__xyz*e_x^e_z + D{x}C__xyz*e_y^e_z' 138 assert sstr(grad>C) == 'D{x}C__x + D{y}C__y + D{z}C__z + D{x}C*e_x + D{y}C*e_y + D{z}C*e_z' 139 140 def test_derivatives_in_spherical_coordinates(self): 141 142 X = r, th, phi = symbols('r theta phi') 143 s3d = Ga('e_r e_theta e_phi', g=[1, r ** 2, r ** 2 * sin(th) ** 2], coords=X, norm=True) 144 er, eth, ephi = s3d.mv() 145 grad = s3d.grad 146 147 # like sympy.sstr but using our printer 148 sstr = lambda x: GaPrinter().doprint(x) 149 150 f = s3d.mv('f', 'scalar', f=True) 151 A = s3d.mv('A', 'vector', f=True) 152 B = s3d.mv('B', 'bivector', f=True) 153 154 assert sstr(f) == 'f' 155 assert sstr(A) == 'A__r*e_r + A__theta*e_theta + A__phi*e_phi' 156 assert sstr(B) == 'B__rtheta*e_r^e_theta + B__rphi*e_r^e_phi + B__thetaphi*e_theta^e_phi' 157 158 assert sstr(grad*f) == 'D{r}f*e_r + D{theta}f*e_theta/r + D{phi}f*e_phi/(r*sin(theta))' 159 assert sstr((grad|A).simplify()) == '(r*D{r}A__r + 2*A__r + A__theta/tan(theta) + D{theta}A__theta + D{phi}A__phi/sin(theta))/r' 160 assert sstr(-s3d.I()*(grad^A)) == '(A__phi/tan(theta) + D{theta}A__phi - D{phi}A__theta/sin(theta))*e_r/r + (-r*D{r}A__phi - A__phi + D{phi}A__r/sin(theta))*e_theta/r + (r*D{r}A__theta + A__theta - D{theta}A__r)*e_phi/r' 161 162 assert latex(grad) == r'\boldsymbol{e}_{r} \frac{\partial}{\partial r} + \boldsymbol{e}_{\theta } \frac{1}{r} \frac{\partial}{\partial \theta } + \boldsymbol{e}_{\phi } \frac{1}{r \sin{\left (\theta \right )}} \frac{\partial}{\partial \phi }' 163 assert latex(B|(eth^ephi)) == r'- B^{\theta \phi } {\left (r,\theta ,\phi \right )}' 164 165 assert sstr(grad^B) == '(r*D{r}B__thetaphi - B__rphi/tan(theta) + 2*B__thetaphi - D{theta}B__rphi + D{phi}B__rtheta/sin(theta))*e_r^e_theta^e_phi/r' 166 167 def test_norm_flag(self): 168 # gh-466 169 rho = symbols('rho', positive=True) 170 theta, phi = symbols('theta phi', real=True) 171 sph_coords = (rho, theta, phi) 172 173 p = (rho*sin(theta)*cos(phi), rho*sin(theta)*sin(phi), rho*cos(theta)) 174 g_sph = Ga('e', coords=sph_coords, X=p, norm=True) 175 assert g_sph.g == eye(3) 176 177 def test_norm_flag_subspace(self): 178 # gh-466 179 # the coordinates here describe a submanifold 180 R = symbols('R', positive=True) 181 theta, z = symbols('theta z', real=True) 182 cyl2_coords = (theta, z) 183 184 p = (R*cos(theta), R*sin(theta), z) 185 g_cyl2 = Ga('e', coords=cyl2_coords, X=p, norm=True) 186 assert g_cyl2.g == eye(2) 187 188 def test_rounding_numerical_components(self): 189 190 o3d = Ga('e_x e_y e_z', g=[1, 1, 1]) 191 ex, ey, ez = o3d.mv() 192 193 X = 1.2*ex+2.34*ey+0.555*ez 194 Y = 0.333*ex+4*ey+5.3*ez 195 196 assert str(X) == '1.2*e_x + 2.34*e_y + 0.555*e_z' 197 assert str(Nga(X, 2)) == '1.2*e_x + 2.3*e_y + 0.55*e_z' 198 assert str(X*Y) == '12.7011 + 4.02078*e_x^e_y + 6.175185*e_x^e_z + 10.182*e_y^e_z' 199 assert str(Nga(X*Y, 2)) == '13.0 + 4.0*e_x^e_y + 6.2*e_x^e_z + 10.0*e_y^e_z' 200 201 def test_noneuclidian_distance_calculation(self): 202 from sympy import solve, sqrt 203 204 g = '0 # #,# 0 #,# # 1' 205 necl = Ga('X Y e',g=g) 206 X, Y, e = necl.mv() 207 208 assert str((X^Y)*(X^Y)) == '(X.Y)**2' 209 210 L = X^Y^e 211 B = (L*e).expand().blade_rep() # D&L 10.152 212 assert str(B) == 'X^Y - (Y.e)*X^e + (X.e)*Y^e' 213 Bsq = B*B 214 assert str(Bsq) == '(X.Y)*((X.Y) - 2*(X.e)*(Y.e))' 215 Bsq = Bsq.scalar() 216 assert str(B) == 'X^Y - (Y.e)*X^e + (X.e)*Y^e' 217 218 BeBr = B*e*B.rev() 219 assert str(BeBr) == '(X.Y)*(-(X.Y) + 2*(X.e)*(Y.e))*e' 220 assert str(B*B) == '(X.Y)*((X.Y) - 2*(X.e)*(Y.e))' 221 assert str(L*L) == '(X.Y)*((X.Y) - 2*(X.e)*(Y.e))' # D&L 10.153 222 223 s, c, Binv, M, S, C, alpha = symbols('s c (1/B) M S C alpha') 224 225 XdotY = necl.g[0, 1] 226 Xdote = necl.g[0, 2] 227 Ydote = necl.g[1, 2] 228 229 Bhat = Binv*B # D&L 10.154 230 R = c+s*Bhat # Rotor R = exp(alpha*Bhat/2) 231 assert str(R) == 'c + (1/B)*s*X^Y - (Y.e)*(1/B)*s*X^e + (X.e)*(1/B)*s*Y^e' 232 233 Z = R*X*R.rev() # D&L 10.155 234 Z.obj = expand(Z.obj) 235 Z.obj = Z.obj.collect([Binv, s, c, XdotY]) 236 assert str(Z) == '((X.Y)**2*(1/B)**2*s**2 - 2*(X.Y)*(X.e)*(Y.e)*(1/B)**2*s**2 + 2*(X.Y)*(1/B)*c*s - 2*(X.e)*(Y.e)*(1/B)*c*s + c**2)*X + 2*(X.e)**2*(1/B)*c*s*Y + 2*(X.Y)*(X.e)*(1/B)*s*(-(X.Y)*(1/B)*s + 2*(X.e)*(Y.e)*(1/B)*s - c)*e' 237 W = Z|Y 238 # From this point forward all calculations are with sympy scalars 239 W = W.scalar() 240 assert str(W) == '(X.Y)**3*(1/B)**2*s**2 - 4*(X.Y)**2*(X.e)*(Y.e)*(1/B)**2*s**2 + 2*(X.Y)**2*(1/B)*c*s + 4*(X.Y)*(X.e)**2*(Y.e)**2*(1/B)**2*s**2 - 4*(X.Y)*(X.e)*(Y.e)*(1/B)*c*s + (X.Y)*c**2' 241 W = expand(W) 242 W = simplify(W) 243 W = W.collect([s*Binv]) 244 245 M = 1/Bsq 246 W = W.subs(Binv**2, M) 247 W = simplify(W) 248 Bmag = sqrt(XdotY**2-2*XdotY*Xdote*Ydote) 249 W = W.collect([Binv*c*s, XdotY]) 250 251 #Double angle substitutions 252 253 W = W.subs(2*XdotY**2-4*XdotY*Xdote*Ydote, 2/(Binv**2)) 254 W = W.subs(2*c*s, S) 255 W = W.subs(c**2, (C+1)/2) 256 W = W.subs(s**2, (C-1)/2) 257 W = simplify(W) 258 W = W.subs(Binv, 1/Bmag) 259 W = expand(W) 260 261 assert str(W.simplify()) == '(X.Y)*C - (X.e)*(Y.e)*C + (X.e)*(Y.e) + S*sqrt((X.Y)*((X.Y) - 2*(X.e)*(Y.e)))' 262 263 Wd = collect(W, [C, S], exact=True, evaluate=False) 264 265 Wd_1 = Wd[one] 266 Wd_C = Wd[C] 267 Wd_S = Wd[S] 268 269 assert str(Wd_1) == '(X.e)*(Y.e)' 270 assert str(Wd_C) == '(X.Y) - (X.e)*(Y.e)' 271 272 assert str(Wd_S) == 'sqrt((X.Y)**2 - 2*(X.Y)*(X.e)*(Y.e))' 273 assert str(Bmag) == 'sqrt((X.Y)**2 - 2*(X.Y)*(X.e)*(Y.e))' 274 275 Wd_1 = Wd_1.subs(Binv, 1/Bmag) 276 Wd_C = Wd_C.subs(Binv, 1/Bmag) 277 Wd_S = Wd_S.subs(Binv, 1/Bmag) 278 279 lhs = Wd_1+Wd_C*C 280 rhs = -Wd_S*S 281 lhs = lhs**2 282 rhs = rhs**2 283 W = expand(lhs-rhs) 284 W = expand(W.subs(1/Binv**2, Bmag**2)) 285 W = expand(W.subs(S**2, C**2-1)) 286 W = W.collect([C, C**2], evaluate=False) 287 288 a = simplify(W[C**2]) 289 b = simplify(W[C]) 290 c = simplify(W[one]) 291 292 assert str(a) == '(X.e)**2*(Y.e)**2' 293 assert str(b) == '2*(X.e)*(Y.e)*((X.Y) - (X.e)*(Y.e))' 294 assert str(c) == '(X.Y)**2 - 2*(X.Y)*(X.e)*(Y.e) + (X.e)**2*(Y.e)**2' 295 296 x = Symbol('x') 297 C = solve(a*x**2+b*x+c, x)[0] 298 assert str(expand(simplify(expand(C)))) == '-(X.Y)/((X.e)*(Y.e)) + 1' 299 300 def test_conformal_representations_of_circles_lines_spheres_and_planes(self): 301 global n, nbar 302 303 g = '1 0 0 0 0,0 1 0 0 0,0 0 1 0 0,0 0 0 0 2,0 0 0 2 0' 304 305 cnfml3d = Ga('e_1 e_2 e_3 n nbar', g=g) 306 307 e1, e2, e3, n, nbar = cnfml3d.mv() 308 309 e = n+nbar 310 311 #conformal representation of points 312 313 A = make_vector(e1, ga=cnfml3d) # point a = (1, 0, 0) A = F(a) 314 B = make_vector(e2, ga=cnfml3d) # point b = (0, 1, 0) B = F(b) 315 C = make_vector(-e1, ga=cnfml3d) # point c = (-1, 0, 0) C = F(c) 316 D = make_vector(e3, ga=cnfml3d) # point d = (0, 0, 1) D = F(d) 317 X = make_vector('x', 3, ga=cnfml3d) 318 319 assert str(A) == 'e_1 + n/2 - nbar/2' 320 assert str(B) == 'e_2 + n/2 - nbar/2' 321 assert str(C) == '-e_1 + n/2 - nbar/2' 322 assert str(D) == 'e_3 + n/2 - nbar/2' 323 assert str(X) == 'x1*e_1 + x2*e_2 + x3*e_3 + (x1**2/2 + x2**2/2 + x3**2/2)*n - nbar/2' 324 325 assert str((A^B^C^X)) == '-x3*e_1^e_2^e_3^n + x3*e_1^e_2^e_3^nbar + (x1**2/2 + x2**2/2 + x3**2/2 - 1/2)*e_1^e_2^n^nbar' 326 assert str((A^B^n^X)) == '-x3*e_1^e_2^e_3^n + (x1/2 + x2/2 - 1/2)*e_1^e_2^n^nbar + x3*e_1^e_3^n^nbar/2 - x3*e_2^e_3^n^nbar/2' 327 assert str((((A^B)^C)^D)^X) == '(-x1**2/2 - x2**2/2 - x3**2/2 + 1/2)*e_1^e_2^e_3^n^nbar' 328 assert str((A^B^n^D^X)) == '(-x1/2 - x2/2 - x3/2 + 1/2)*e_1^e_2^e_3^n^nbar' 329 330 L = (A^B^e)^X 331 332 assert str(L) == '-x3*e_1^e_2^e_3^n - x3*e_1^e_2^e_3^nbar + (-x1**2/2 + x1 - x2**2/2 + x2 - x3**2/2 - 1/2)*e_1^e_2^n^nbar + x3*e_1^e_3^n^nbar - x3*e_2^e_3^n^nbar' 333 334 def test_properties_of_geometric_objects(self): 335 336 global n, nbar 337 338 g = '# # # 0 0,'+ \ 339 '# # # 0 0,'+ \ 340 '# # # 0 0,'+ \ 341 '0 0 0 0 2,'+ \ 342 '0 0 0 2 0' 343 344 c3d = Ga('p1 p2 p3 n nbar', g=g) 345 346 p1, p2, p3, n, nbar = c3d.mv() 347 348 P1 = F(p1) 349 P2 = F(p2) 350 P3 = F(p3) 351 352 L = P1^P2^n 353 delta = (L|n)|nbar 354 assert str(delta) == '2*p1 - 2*p2' 355 356 C = P1^P2^P3 357 delta = ((C^n)|n)|nbar 358 assert str(delta) == '2*p1^p2 - 2*p1^p3 + 2*p2^p3' 359 assert str((p2-p1)^(p3-p1)) == 'p1^p2 - p1^p3 + p2^p3' 360 361 def test_extracting_vectors_from_conformal_2_blade(self): 362 363 g = '0 -1 #,'+ \ 364 '-1 0 #,'+ \ 365 '# # #' 366 367 e2b = Ga('P1 P2 a', g=g) 368 369 P1, P2, a = e2b.mv() 370 371 B = P1^P2 372 Bsq = B*B 373 assert str(Bsq) == '1' 374 ap = a-(a^B)*B 375 assert str(ap) == '-(P2.a)*P1 - (P1.a)*P2' 376 377 Ap = ap+ap*B 378 Am = ap-ap*B 379 380 assert str(Ap) == '-2*(P2.a)*P1' 381 assert str(Am) == '-2*(P1.a)*P2' 382 383 assert str(Ap*Ap) == '0' 384 assert str(Am*Am) == '0' 385 386 aB = a|B 387 assert str(aB) == '-(P2.a)*P1 + (P1.a)*P2' 388 389 def test_ReciprocalFrame(self): 390 ga, *basis = Ga.build('e*u|v|w') 391 392 r_basis = ga.ReciprocalFrame(basis) 393 394 for i, base in enumerate(basis): 395 for r_i, r_base in enumerate(r_basis): 396 if i == r_i: 397 assert (base | r_base).simplify() == 1 398 else: 399 assert (base | r_base).simplify() == 0 400 401 def test_ReciprocalFrame_2d(self): 402 """Test reciprocal frame in 2D to catch even-dimension sign bugs.""" 403 ga2, e1, e2 = Ga.build('e*1|2', g=[1, 1]) 404 basis2 = [e1 + e2, e1 - e2] 405 r_basis2 = ga2.ReciprocalFrame(basis2) 406 for i, base in enumerate(basis2): 407 for j, rbase in enumerate(r_basis2): 408 assert (base | rbase).simplify() == (1 if i == j else 0) 409 410 def test_ReciprocalFrame_append(self): 411 ga, *basis = Ga.build('e*u|v|w') 412 *r_basis, E_sq = ga.ReciprocalFrame(basis, mode='append') 413 414 for i, base in enumerate(basis): 415 for r_i, r_base in enumerate(r_basis): 416 if i == r_i: 417 assert (base | r_base).simplify() == E_sq 418 else: 419 assert (base | r_base).simplify() == 0 420 421 # anything that isn't 'norm' means 'append', but this is deprecated 422 with pytest.warns(DeprecationWarning): 423 assert ga.ReciprocalFrame(basis, mode='nonsense') == (*r_basis, E_sq) 424 425 def test_reciprocal_frame_test(self): 426 427 g = '1 # #,'+ \ 428 '# 1 #,'+ \ 429 '# # 1' 430 431 g3dn = Ga('e1 e2 e3', g=g) 432 433 e1, e2, e3 = g3dn.mv() 434 435 E = e1^e2^e3 436 Esq = (E*E).scalar() 437 assert str(E) == 'e1^e2^e3' 438 assert str(Esq) == '(e1.e2)**2 - 2*(e1.e2)*(e1.e3)*(e2.e3) + (e1.e3)**2 + (e2.e3)**2 - 1' 439 Esq_inv = 1/Esq 440 441 E1 = (e2^e3)*E 442 E2 = (-1)*(e1^e3)*E 443 E3 = (e1^e2)*E 444 445 assert str(E1) == '((e2.e3)**2 - 1)*e1 + ((e1.e2) - (e1.e3)*(e2.e3))*e2 + (-(e1.e2)*(e2.e3) + (e1.e3))*e3' 446 assert str(E2) == '((e1.e2) - (e1.e3)*(e2.e3))*e1 + ((e1.e3)**2 - 1)*e2 + (-(e1.e2)*(e1.e3) + (e2.e3))*e3' 447 assert str(E3) == '(-(e1.e2)*(e2.e3) + (e1.e3))*e1 + (-(e1.e2)*(e1.e3) + (e2.e3))*e2 + ((e1.e2)**2 - 1)*e3' 448 449 w = (E1|e2) 450 w = w.expand() 451 assert str(w) == '0' 452 453 w = (E1|e3) 454 w = w.expand() 455 assert str(w) == '0' 456 457 w = (E2|e1) 458 w = w.expand() 459 assert str(w) == '0' 460 461 w = (E2|e3) 462 w = w.expand() 463 assert str(w) == '0' 464 465 w = (E3|e1) 466 w = w.expand() 467 assert str(w) == '0' 468 469 w = (E3|e2) 470 w = w.expand() 471 assert str(w) == '0' 472 473 w = (E1|e1) 474 w = (w.expand()).scalar() 475 Esq = expand(Esq) 476 assert str(simplify(w/Esq)) == '1' 477 478 w = (E2|e2) 479 w = (w.expand()).scalar() 480 assert str(simplify(w/Esq)) == '1' 481 482 w = (E3|e3) 483 w = (w.expand()).scalar() 484 assert str(simplify(w/Esq)) == '1' 485 486 def test_make_grad(self): 487 ga, e_1, e_2, e_3 = Ga.build('e*1|2|3', g=[1, 1, 1], coords=symbols('x y z')) 488 r = ga.mv(ga.coord_vec) 489 assert ga.make_grad(r) == ga.grad 490 assert ga.make_grad(r, cmpflg=True) == ga.rgrad 491 492 x = ga.mv('x', 'vector') 493 B = ga.mv('B', 'bivector') 494 dx = ga.make_grad(x) 495 dB = ga.make_grad(B) 496 497 # GA4P, eq. (6.29) 498 for a in [ga.mv(1), e_1, e_1^e_2]: 499 r = a.i_grade 500 assert dx * (x ^ a) == (ga.n - r) * a 501 assert dx * (x * a) == ga.n * a 502 503 # derivable via the product rule 504 assert dx * (x*x) == 2*x 505 assert dx * (x*x*x) == (2*x)*x + (x*x)*ga.n 506 507 assert dB * (B*B) == 2*B 508 assert dB * (B*B*B) == (2*B)*B + (B*B)*ga.n 509 510 # an arbitrary chained expression to check we do not crash 511 assert dB * dx * (B * x) == -3 512 assert dx * dB * (x * B) == -3 513 assert dx * dB * (B * x) == 9 514 assert dB * dx * (x * B) == 9 515 516 @pytest.mark.parametrize('g', [ 517 pytest.param(None, id='generic'), 518 pytest.param([1, 1, 1], id='ortho') 519 ]) 520 def test_reciprocal_blades(self, g): 521 ga = Ga('e*1|2|3', g=g) 522 523 for b1 in ga.blades.flat: 524 for b2 in ga.blades.flat: 525 rb2 = ga._reciprocal_blade_dict[b2] 526 527 if b1 == b2: 528 assert ga.scalar_product(b1, rb2).simplify() == S.One 529 else: 530 assert ga.scalar_product(b1, rb2).simplify() == S.Zero 531 532 def test_er_blade(self): 533 """Test er_blade product modes in orthonormal 3D (issue 140).""" 534 ga, e1, e2, e3 = Ga.build('e*1|2|3', g=[1, 1, 1]) 535 er1, er2, er3 = ga.r_basis 536 blade12, blade13, blade23 = ga.blades[2] 537 538 # geometric product: e1*(e1^e2) = e2, (e1^e2)*e1 = -e2 539 assert ga.er_blade(er1, blade12, mode='*', left=True) == ga.blades[1][1] 540 assert ga.er_blade(er1, blade12, mode='*', left=False) == -ga.blades[1][1] 541 # e1*(e2^e3) = e1^e2^e3 542 assert ga.er_blade(er1, blade23, mode='*', left=True) == ga.blades[3][0] 543 544 # wedge product: e1^(e1^e2) = 0, e3^(e1^e2) = e1^e2^e3 545 assert ga.er_blade(er1, blade12, mode='^', left=True) == 0 546 assert ga.er_blade(er3, blade12, mode='^', left=True) == ga.blades[3][0] 547 548 # Hestenes dot: e1|(e1^e2) = e2, e2|(e1^e2) = -e1 549 assert ga.er_blade(er1, blade12, mode='|', left=True) == ga.blades[1][1] 550 assert ga.er_blade(er2, blade12, mode='|', left=True) == -ga.blades[1][0] 551 552 def test_metric_collect(self): 553 ga = Ga('e*1|2', g=[1, 1]) 554 e1, e2 = ga.basis 555 556 assert metric.collect(2*e1 + e2, [e1]) == 2*e1 + e2 557 558 def test_dual_mode(self): 559 ga, e1, e2 = Ga.build('e*1|2', g=[1, 1]) 560 561 default = Ga.dual_mode_value 562 assert default == 'I+' 563 with pytest.raises(ValueError): 564 Ga.dual_mode('illegal') 565 566 d_default = e1.dual() 567 568 # note: this is a global setting, so we have to make sure we put it back 569 try: 570 Ga.dual_mode('I-') 571 d_negated = e1.dual() 572 finally: 573 Ga.dual_mode(default) 574 575 assert d_negated == -d_default 576 577 def test_basis_dict(self): 578 ga = Ga('e*1|2', g=[1, 1]) 579 b = ga.bases_dict() 580 assert b == { 581 'e1': ga.blades[1][0], 582 'e2': ga.blades[1][1], 583 'e12': ga.blades[2][0], 584 } 585 586 def test_single_basis(self): 587 # dual numbers 588 g, delta = Ga.build('delta,', g=[0]) 589 assert delta*delta == 0 590 # which work for automatic differentiation 591 x = Symbol('x') 592 xd = x + delta 593 f = lambda x: x**3 + 2*x*2 + 1 594 assert f(xd) == f(x) + f(x).diff(x) * delta 595 596 def test_no_basis(self): 597 # real numbers! 598 g = Ga('', g=[]) 599 one = g.mv(1) 600 assert (3*one) ^ (2*one) == (6*one) 601 602 def test_deprecations(self): 603 coords = symbols('x y z') 604 ga, e_1, e_2, e_3 = Ga.build('e*1|2|3', coords=coords) 605 606 # none of these have the scalar as their first element, which is why 607 # they're deprecated. 608 with pytest.warns(DeprecationWarning): 609 assert ga.blades_lst[0] == e_1.obj 610 with pytest.warns(DeprecationWarning): 611 assert ga.bases_lst[0] == e_1.obj 612 with pytest.warns(DeprecationWarning): 613 assert ga.indexes_lst[1] == (1,) 614 615 # deprecated to reduce the number of similar members 616 with pytest.warns(DeprecationWarning): 617 ga.blades_to_indexes 618 with pytest.warns(DeprecationWarning): 619 ga.bases_to_indexes 620 with pytest.warns(DeprecationWarning): 621 ga.blades_to_indexes_dict 622 with pytest.warns(DeprecationWarning): 623 ga.bases_to_indexes_dict 624 with pytest.warns(DeprecationWarning): 625 ga.indexes_to_blades 626 with pytest.warns(DeprecationWarning): 627 ga.indexes_to_bases 628 629 # all the above are deprecated in favor of these two, which are _not_ 630 # deprecated 631 ga.indexes_to_blades_dict 632 ga.indexes_to_bases_dict 633 634 # deprecated to reduce the number of similar members 635 with pytest.warns(DeprecationWarning): 636 ga.basic_mul_table 637 with pytest.warns(DeprecationWarning): 638 ga.basic_mul_keys 639 with pytest.warns(DeprecationWarning): 640 ga.basic_mul_values 641 642 # all derived from 643 ga.basic_mul_table_dict 644 645 # deprecated to reduce the number of similar members 646 with pytest.warns(DeprecationWarning): 647 ga.blade_expansion 648 with pytest.warns(DeprecationWarning): 649 ga.base_expansion 650 651 # all derived from 652 ga.blade_expansion_dict 653 ga.base_expansion_dict 654 655 with pytest.warns(DeprecationWarning): 656 import galgebra.utils 657 658 # aliases 659 with pytest.warns(DeprecationWarning): 660 assert ga.X() 661 662 # derived from 663 ga.coord_vec 664 665 # aliases 666 with pytest.warns(DeprecationWarning): 667 ga.lt_x 668 with pytest.warns(DeprecationWarning): 669 ga.lt_coords 670 671 # useless method relating to those aliases 672 from galgebra.lt import Lt 673 with pytest.warns(DeprecationWarning): 674 Lt.setup(ga) 675 676 # more aliases 677 with pytest.warns(DeprecationWarning): 678 ga.mul_table_dict 679 with pytest.warns(DeprecationWarning): 680 ga.wedge_table_dict 681 with pytest.warns(DeprecationWarning): 682 ga.dot_table_dict 683 with pytest.warns(DeprecationWarning): 684 ga.left_contract_table_dict 685 with pytest.warns(DeprecationWarning): 686 ga.right_contract_table_dict 687 with pytest.warns(DeprecationWarning): 688 ga.basic_mul_table_dict 689 690 with pytest.warns(DeprecationWarning): 691 assert ga.geometric_product_basis_blades((e_1.obj, e_2.obj)) == (e_1 * e_2).obj 692 with pytest.warns(DeprecationWarning): 693 assert ga.non_orthogonal_bases_products((e_1.obj, e_2.obj)) == (e_1 * e_2).base_rep().obj 694 with pytest.warns(DeprecationWarning): 695 assert ga.wedge_product_basis_blades((e_1.obj, e_2.obj)) == (e_1 ^ e_2).obj 696 e_12 = e_1 ^ e_2 697 with pytest.warns(DeprecationWarning): 698 assert ga.non_orthogonal_dot_product_basis_blades((e_1.obj, e_12.obj), mode='|') == (e_1 | e_12).obj 699 with pytest.warns(DeprecationWarning): 700 assert ga.non_orthogonal_dot_product_basis_blades((e_1.obj, e_12.obj), mode='<') == (e_1 < e_12).obj 701 with pytest.warns(DeprecationWarning): 702 assert ga.non_orthogonal_dot_product_basis_blades((e_1.obj, e_12.obj), mode='>') == (e_1 > e_12).obj 703 704 # these methods don't do anything 705 with pytest.warns(DeprecationWarning): 706 ga.inverse_metric() 707 with pytest.warns(DeprecationWarning): 708 ga.derivatives_of_g() 709 710 # test the member that is nonsense unless in an orthonormal algebra 711 ga_ortho, e_1, e_2, e_3 = Ga.build('e*1|2|3', g=[1, 1, 1]) 712 e_12 = e_1 ^ e_2 713 with pytest.warns(DeprecationWarning): 714 assert ga_ortho.dot_product_basis_blades((e_1.obj, e_12.obj), mode='|') == (e_1 | e_12).obj 715 with pytest.warns(DeprecationWarning): 716 assert ga_ortho.dot_product_basis_blades((e_1.obj, e_12.obj), mode='<') == (e_1 < e_12).obj 717 with pytest.warns(DeprecationWarning): 718 assert ga_ortho.dot_product_basis_blades((e_1.obj, e_12.obj), mode='>') == (e_1 > e_12).obj 719 720 def test_Cl(self): 721 """Test the Cl(p, q, r) interface (issue 524).""" 722 from galgebra.interop import Cl 723 724 # 3D Euclidean: Cl(3) 725 ga3, e1, e2, e3 = Cl(3) 726 assert e1 * e1 == ga3.mv(1) 727 assert e2 * e2 == ga3.mv(1) 728 assert e3 * e3 == ga3.mv(1) 729 730 # 2D Minkowski: Cl(1, 1) 731 ga11, e1, e2 = Cl(1, 1) 732 assert e1 * e1 == ga11.mv(1) 733 assert e2 * e2 == ga11.mv(-1) 734 735 # Degenerate: Cl(2, 0, 1) 736 ga201, e1, e2, e3 = Cl(2, 0, 1) 737 assert e1 * e1 == ga201.mv(1) 738 assert e2 * e2 == ga201.mv(1) 739 assert e3 * e3 == ga201.mv(0) 740 741 # Import from top-level package 742 from galgebra import Cl as Cl2 743 ga_top, e1_top, e2_top = Cl2(2) 744 assert e1_top * e1_top == ga_top.mv(1) 745 746 # Error on zero dimension 747 with pytest.raises(ValueError): 748 Cl(0) 749 750 def test_Cl_kingdon(self): 751 """Test the kingdon-convention Cl (0-indexed, Iinv+ dual).""" 752 from galgebra.interop.kingdon import Cl as KCl 753 from galgebra.ga import Ga 754 755 # Save and restore global dual mode 756 saved = Ga.dual_mode_value 757 try: 758 # 3D PGA: Cl(3, 0, 1) with 0-indexed basis 759 ga, e0, e1, e2, e3 = KCl(3, 0, 1) 760 assert e0 * e0 == ga.mv(1) 761 assert e3 * e3 == ga.mv(0) 762 # Dual mode should be Iinv+ 763 assert Ga.dual_mode_value == 'Iinv+' 764 finally: 765 Ga.dual_mode(saved) 766 767 def test_Cl_dual_mode_reset(self): 768 """galgebra.interop.Cl must reset dual mode to I+ after kingdon.Cl set Iinv+.""" 769 from galgebra.interop.kingdon import Cl as KCl 770 from galgebra.interop import Cl 771 from galgebra.ga import Ga 772 773 saved = Ga.dual_mode_value 774 try: 775 KCl(3) 776 assert Ga.dual_mode_value == 'Iinv+' 777 Cl(3) 778 assert Ga.dual_mode_value == 'I+' 779 finally: 780 Ga.dual_mode(saved)