Skip to content

Commit 1ce1340

Browse files
utiberiousutiberious
committed
fix: deferred-simp approach for curvi_linear_latex (pygae#576)
Replaces the `Simp.profile([simplify])` default with a two-phase strategy to avoid the SymPy 1.13 TR3 O(N·M) traversal (PR #26390): 1. **Build phase** — `_no_simp_build()` context manager sets an identity simplifier (`Simp.profile([lambda e: e])`) around each `Ga.build()` call. The ~70 intermediate `Simp.apply` calls during metric normalisation now do nothing, cutting the build time from ~25 s/call to < 0.01 s/call. 2. **Output phase** — each result expression (gradient, divergence, curl, Laplacian, grad-of-bivector) is explicitly simplified once with `_ts()` = `Mv.simplify(modes=lambda e: trigsimp(e, method='old'))` before formatting. `trigsimp(method='old')` avoids `fu.py`/TR3 entirely. Alternative considered: replacing the explicit `trigsimp(method='old')` in `_ts()` with plain `simplify()` to recover the pre-SymPy-1.13 canonical output form. This does **not** work — the un-simplified metric components folded into the gradient/divergence/curl expressions are large enough that TR3 is just as slow on the final result as it was on the intermediate metric components, causing the same > 10 min hang. Closes pygae#576 (alongside pygae#590)
1 parent 444fdd5 commit 1ce1340

1 file changed

Lines changed: 86 additions & 41 deletions

File tree

‎examples/LaTeX/curvi_linear_latex.py‎

Lines changed: 86 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,53 @@
11
from __future__ import print_function
22
import sys
3+
from contextlib import contextmanager
34

4-
from sympy import symbols,sin,cos,sinh,cosh
5+
from sympy import symbols,sin,cos,sinh,cosh,trigsimp
56
from galgebra.ga import Ga
7+
from galgebra.metric import Simp
68
from galgebra.printer import Format, xpdf, Print_Function, Eprint
79

10+
11+
@contextmanager
12+
def _no_simp_build():
13+
"""Suppress Simp during Ga.build; restore trigsimp(method='old') afterwards.
14+
15+
``Ga.build(norm=True)`` calls ``Simp.apply`` ~70 times during metric
16+
normalisation. Since SymPy 1.13 (PR #26390) each call triggers a slow
17+
O(N·M) ``.replace()`` traversal inside ``TR3``/``futrig`` that is a no-op
18+
for symbolic trig arguments but costs ~25 s per invocation on these
19+
curvilinear-coordinate expressions.
20+
21+
Strategy: use an identity simplifier during the build phase (so the metric
22+
components are stored in unsimplified form), then simplify each final output
23+
expression once with ``trigsimp(method='old')`` via ``_ts()`` before
24+
formatting. ``trigsimp(method='old')`` uses the ``_trigsimp`` code path
25+
which avoids ``fu.py`` / TR3 entirely and runs in < 0.1 s per expression.
26+
27+
Note: calling plain ``simplify()`` on the final results does **not** help —
28+
the unsimplified metric components make the gradient/divergence/curl
29+
expressions large enough that TR3 is just as slow on them as it was on the
30+
intermediate metric components. ``trigsimp(method='old')`` is required to
31+
avoid the SymPy 1.13 slow path entirely.
32+
"""
33+
orig = Simp.modes[:]
34+
Simp.profile([lambda e: e]) # identity — no simplification during build
35+
try:
36+
yield
37+
finally:
38+
Simp.profile(orig)
39+
40+
41+
def _ts(mv):
42+
"""Apply trigsimp(method='old') to each coefficient of a multivector."""
43+
return mv.simplify(modes=lambda e: trigsimp(e, method='old'))
44+
45+
846
def derivatives_in_spherical_coordinates():
947
#Print_Function()
1048
coords = (r,th,phi) = symbols('r theta phi', real=True)
11-
(sp3d,er,eth,ephi) = Ga.build('e_r e_theta e_phi',g=[1,r**2,r**2*sin(th)**2],coords=coords)
49+
with _no_simp_build():
50+
(sp3d,er,eth,ephi) = Ga.build('e_r e_theta e_phi',g=[1,r**2,r**2*sin(th)**2],coords=coords)
1251
grad = sp3d.grad
1352

1453
f = sp3d.mv('f','scalar',f=True)
@@ -21,11 +60,11 @@ def derivatives_in_spherical_coordinates():
2160
print('A =',A)
2261
print('B =',B)
2362

24-
print('grad*f =',grad*f)
25-
print('grad|A =',grad|A)
26-
print('grad\\times A = -I*(grad^A) =',-sp3d.i*(grad^A))
27-
print('%\\nabla^{2}f =',grad|(grad*f))
28-
print('grad^B =',grad^B)
63+
print('grad*f =',_ts(grad*f))
64+
print('grad|A =',_ts(grad|A))
65+
print('grad\\times A = -I*(grad^A) =',_ts(-sp3d.i*(grad^A)))
66+
print('%\\nabla^{2}f =',_ts(grad|(grad*f)))
67+
print('grad^B =',_ts(grad^B))
2968

3069
"""
3170
print '( \\nabla\\W\\nabla )\\bm{e}_{r} =',((grad^grad)*er).trigsimp()
@@ -39,7 +78,8 @@ def derivatives_in_spherical_coordinates():
3978
def derivatives_in_paraboloidal_coordinates():
4079
#Print_Function()
4180
coords = (u,v,phi) = symbols('u v phi', real=True)
42-
(par3d,er,eth,ephi) = Ga.build('e_u e_v e_phi',X=[u*v*cos(phi),u*v*sin(phi),(u**2-v**2)/2],coords=coords,norm=True)
81+
with _no_simp_build():
82+
(par3d,er,eth,ephi) = Ga.build('e_u e_v e_phi',X=[u*v*cos(phi),u*v*sin(phi),(u**2-v**2)/2],coords=coords,norm=True)
4383
grad = par3d.grad
4484

4585
f = par3d.mv('f','scalar',f=True)
@@ -52,10 +92,10 @@ def derivatives_in_paraboloidal_coordinates():
5292
print('A =',A)
5393
print('B =',B)
5494

55-
print('grad*f =',grad*f)
56-
print('grad|A =',grad|A)
57-
(-par3d.i*(grad^A)).Fmt(3,'grad\\times A = -I*(grad^A)')
58-
print('grad^B =',grad^B)
95+
print('grad*f =',_ts(grad*f))
96+
print('grad|A =',_ts(grad|A))
97+
_ts(-par3d.i*(grad^A)).Fmt(3,'grad\\times A = -I*(grad^A)')
98+
print('grad^B =',_ts(grad^B))
5999

60100
return
61101

@@ -64,7 +104,8 @@ def derivatives_in_elliptic_cylindrical_coordinates():
64104
#Print_Function()
65105
a = symbols('a', real=True)
66106
coords = (u,v,z) = symbols('u v z', real=True)
67-
(elip3d,er,eth,ephi) = Ga.build('e_u e_v e_z',X=[a*cosh(u)*cos(v),a*sinh(u)*sin(v),z],coords=coords,norm=True)
107+
with _no_simp_build():
108+
(elip3d,er,eth,ephi) = Ga.build('e_u e_v e_z',X=[a*cosh(u)*cos(v),a*sinh(u)*sin(v),z],coords=coords,norm=True)
68109
grad = elip3d.grad
69110

70111
f = elip3d.mv('f','scalar',f=True)
@@ -77,19 +118,20 @@ def derivatives_in_elliptic_cylindrical_coordinates():
77118
print('A =',A)
78119
print('B =',B)
79120

80-
print('grad*f =',grad*f)
81-
print('grad|A =',grad|A)
82-
print('-I*(grad^A) =',-elip3d.i*(grad^A))
83-
print('grad^B =',grad^B)
121+
print('grad*f =',_ts(grad*f))
122+
print('grad|A =',_ts(grad|A))
123+
print('-I*(grad^A) =',_ts(-elip3d.i*(grad^A)))
124+
print('grad^B =',_ts(grad^B))
84125
return
85126

86127

87128
def derivatives_in_prolate_spheroidal_coordinates():
88129
#Print_Function()
89130
a = symbols('a', real=True)
90131
coords = (xi,eta,phi) = symbols('xi eta phi', real=True)
91-
(ps3d,er,eth,ephi) = Ga.build('e_xi e_eta e_phi',X=[a*sinh(xi)*sin(eta)*cos(phi),a*sinh(xi)*sin(eta)*sin(phi),
92-
a*cosh(xi)*cos(eta)],coords=coords,norm=True)
132+
with _no_simp_build():
133+
(ps3d,er,eth,ephi) = Ga.build('e_xi e_eta e_phi',X=[a*sinh(xi)*sin(eta)*cos(phi),a*sinh(xi)*sin(eta)*sin(phi),
134+
a*cosh(xi)*cos(eta)],coords=coords,norm=True)
93135
grad = ps3d.grad
94136

95137
f = ps3d.mv('f','scalar',f=True)
@@ -102,19 +144,20 @@ def derivatives_in_prolate_spheroidal_coordinates():
102144
print('A =',A)
103145
print('B =',B)
104146

105-
print('grad*f =',grad*f)
106-
print('grad|A =',grad|A)
107-
(-ps3d.i*(grad^A)).Fmt(3,'-I*(grad^A)')
108-
(grad^B).Fmt(3,'grad^B')
147+
print('grad*f =',_ts(grad*f))
148+
print('grad|A =',_ts(grad|A))
149+
_ts(-ps3d.i*(grad^A)).Fmt(3,'-I*(grad^A)')
150+
_ts(grad^B).Fmt(3,'grad^B')
109151
return
110152

111153

112154
def derivatives_in_oblate_spheroidal_coordinates():
113155
Print_Function()
114156
a = symbols('a', real=True)
115157
coords = (xi,eta,phi) = symbols('xi eta phi', real=True)
116-
(os3d,er,eth,ephi) = Ga.build('e_xi e_eta e_phi',X=[a*cosh(xi)*cos(eta)*cos(phi),a*cosh(xi)*cos(eta)*sin(phi),
117-
a*sinh(xi)*sin(eta)],coords=coords,norm=True)
158+
with _no_simp_build():
159+
(os3d,er,eth,ephi) = Ga.build('e_xi e_eta e_phi',X=[a*cosh(xi)*cos(eta)*cos(phi),a*cosh(xi)*cos(eta)*sin(phi),
160+
a*sinh(xi)*sin(eta)],coords=coords,norm=True)
118161
grad = os3d.grad
119162

120163
f = os3d.mv('f','scalar',f=True)
@@ -125,18 +168,19 @@ def derivatives_in_oblate_spheroidal_coordinates():
125168
print('A =',A)
126169
print('B =',B)
127170

128-
print('grad*f =',grad*f)
129-
print('grad|A =',grad|A)
130-
print('-I*(grad^A) =',-os3d.i*(grad^A))
131-
print('grad^B =',grad^B)
171+
print('grad*f =',_ts(grad*f))
172+
print('grad|A =',_ts(grad|A))
173+
print('-I*(grad^A) =',_ts(-os3d.i*(grad^A)))
174+
print('grad^B =',_ts(grad^B))
132175
return
133176

134177

135178
def derivatives_in_bipolar_coordinates():
136179
Print_Function()
137180
a = symbols('a', real=True)
138181
coords = (u,v,z) = symbols('u v z', real=True)
139-
(bp3d,eu,ev,ez) = Ga.build('e_u e_v e_z',X=[a*sinh(v)/(cosh(v)-cos(u)),a*sin(u)/(cosh(v)-cos(u)),z],coords=coords,norm=True)
182+
with _no_simp_build():
183+
(bp3d,eu,ev,ez) = Ga.build('e_u e_v e_z',X=[a*sinh(v)/(cosh(v)-cos(u)),a*sin(u)/(cosh(v)-cos(u)),z],coords=coords,norm=True)
140184
grad = bp3d.grad
141185

142186
f = bp3d.mv('f','scalar',f=True)
@@ -147,20 +191,21 @@ def derivatives_in_bipolar_coordinates():
147191
print('A =',A)
148192
print('B =',B)
149193

150-
print('grad*f =',grad*f)
151-
print('grad|A =',grad|A)
152-
print('-I*(grad^A) =',-bp3d.i*(grad^A))
153-
print('grad^B =',grad^B)
194+
print('grad*f =',_ts(grad*f))
195+
print('grad|A =',_ts(grad|A))
196+
print('-I*(grad^A) =',_ts(-bp3d.i*(grad^A)))
197+
print('grad^B =',_ts(grad^B))
154198
return
155199

156200

157201
def derivatives_in_toroidal_coordinates():
158202
Print_Function()
159203
a = symbols('a', real=True)
160204
coords = (u,v,phi) = symbols('u v phi', real=True)
161-
(t3d,eu,ev,ephi) = Ga.build('e_u e_v e_phi',X=[a*sinh(v)*cos(phi)/(cosh(v)-cos(u)),
162-
a*sinh(v)*sin(phi)/(cosh(v)-cos(u)),
163-
a*sin(u)/(cosh(v)-cos(u))],coords=coords,norm=True)
205+
with _no_simp_build():
206+
(t3d,eu,ev,ephi) = Ga.build('e_u e_v e_phi',X=[a*sinh(v)*cos(phi)/(cosh(v)-cos(u)),
207+
a*sinh(v)*sin(phi)/(cosh(v)-cos(u)),
208+
a*sin(u)/(cosh(v)-cos(u))],coords=coords,norm=True)
164209
grad = t3d.grad
165210

166211
f = t3d.mv('f','scalar',f=True)
@@ -171,10 +216,10 @@ def derivatives_in_toroidal_coordinates():
171216
print('A =',A)
172217
print('B =',B)
173218

174-
print('grad*f =',grad*f)
175-
print('grad|A =',grad|A)
176-
print('-I*(grad^A) =',-t3d.i*(grad^A))
177-
print('grad^B =',grad^B)
219+
print('grad*f =',_ts(grad*f))
220+
print('grad|A =',_ts(grad|A))
221+
print('-I*(grad^A) =',_ts(-t3d.i*(grad^A)))
222+
print('grad^B =',_ts(grad^B))
178223
return
179224

180225

0 commit comments

Comments
 (0)