5
\$\begingroup\$

Please be gentle with me-- I just began learning to code a few weeks ago as a hobby in order to support my other hobby of learning general relativity, so this is the first bit of code I've ever written. I would love to keep getting better at coding, though, so any feedback on ways to improve my code would be much appreciated.

What this code does is take as input the number of dimensions of a manifold, the coordinate labels being used, and the components of a metric, and outputs the non-zero components of the metric (exactly what was input, it just looks prettier), and also of the inverse metric, derivatives of the metric, the Christoffel symbols, the derivatives of the Christoffel symbols, the Riemann curvature tensor, the Ricci curvature tensor, the Ricci scalar, and the Einstein tensor (with 2 covariant indices, but also with 1 contravariant and 1 covariant).

For those of you who run the code, here are some useful tips on the user inputs (I will also include an example at the bottom): When inputting metric components, you can use '^' instead of '**' for exponents, and when multiplying a number by a symbol or something in parentheses, you don't need to include a '*'. Feel free to include undefined functions- just make sure you include its arguments if you want it to be differentiated correctly (e.g. if you want a function that will have a non-zero derivative of x, then type 'f(x)' in your expression instead of just 'f'). Also feel free to use greek letters (spelled out in English) when inputting coordinate labels and/or functions in your metric components.

Here is the code:

from sympy import *
from dataclasses import dataclass
from IPython.display import display as Idisplay
from IPython.display import Math

greek = ['alpha', 'beta', 'gamma', 'Gamma', 'delta', 'Delta', 'epsilon',
         'varepsilon', 'zeta', 'eta', 'theta', 'vartheta', 'Theta', 'iota',
         'kappa', 'lambda', 'Lambda', 'mu', 'nu', 'xi', 'Xi', 'pi', 'Pi',
         'rho', 'varrho', 'sigma', 'Sigma', 'tau', 'upsilon', 'Upsilon',
         'phi', 'varphi', 'Phi', 'chi', 'psi', 'Psi', 'omega', 'Omega']

n = int(input('Enter the number of dimensions:\n'))
coords = []
for i in range(n):
    coords.append(Symbol(str(input('Enter coordinate label %d:\n' % i))))


@dataclass(frozen=False, order=True)
class Tensor:
    name: str
    symbol: str
    key: str
    components: list
    
    def rank(self):
        return self.key.count('*')
    
    def tensor_zeros(self, t=0):
        for i in range(self.rank()):
            t = [t,] * n
        return MutableDenseNDimArray(t)
    
    def coord_id(self, o):
        a = []
        for i in range(self.rank()):
            c = int(o/(n**(self.rank() - i - 1)))
            a.append(str(coords[c]))
            if any(letter in a[i] for letter in greek) is True:
                a[i] = '\\' + a[i] + ' '
            o -= c * (n**(self.rank() - i - 1))  
        x = self.key
        w = 0
        for i in x:
            if i == '*':
                x = x.replace('*', a[w], 1)
                w += 1
        return self.symbol + x
    
    def print_tensor(self):
        for o in range(len(self.components)):
            if self.components[o] != 0:
                Idisplay(Math(latex(Eq(Symbol(self.coord_id(o)),
                                       self.components[o]))))
        print('\n\n')


def assign(instance, thing):
    instance.components = thing.reshape(len(thing)).tolist()


def fix_input(expr):
    expr = expr.replace('^', '**')
    for i in range(len(expr)-1):
        if expr[i].isnumeric() and (expr[i+1].isalpha() or
                                    expr[i+1] == '('):
            expr = expr[:i+1] + '*' + expr[i+1:]
    return expr


metric = Tensor('metric tensor', 'g', '_**', [])
metric_inv = Tensor('inverse of metric tensor', 'g', '__**', [])
metric_d = Tensor('partial derivative of metric tensor', 'g', '_**,*', [])
Christoffel = Tensor('Christoffel symbol - 2nd kind', 'Gamma', '__*_**', [])
Christoffel_d = Tensor('partial derivative of Christoffel symbol',
                       'Gamma', '__*_**,*', [])
Riemann = Tensor('Riemann curvature tensor', 'R', '__*_***', [])
Ricci = Tensor('Ricci curvature tensor', 'R', '_**', [])
Einstein = Tensor('Einstein tensor', 'G', '_**', [])
Einstein_alt = Tensor('Einstein tensor', 'G', '__*_*', [])

# user inputs metric:
g = eye(n)
while True:
    diag = str(input('Is metric diagonal? y for yes, n for no\n')).lower()
    if diag == 'y':
        for i in range(n):
            g[i, i] = sympify(fix_input(str(input(
                'What is g_[%s%s]?\n' % (str(coords[i]), str(coords[i])
                )))))
    else:
        for i in range(n):
            for j in range(i, n):
                g[i, j] = sympify(fix_input(str(input(
                    'What is g_[%s%s]?\n' % (str(coords[i]), str(coords[j])
                    )))))
                g[j, i] = g[i, j]
    if g.det() == 0:
        print('\nMetric is singular, try again!\n')
        continue
    else:
        break

# calculate everything:
# inverse metric:
g_inv = MutableDenseNDimArray(g.inv())
assign(metric_inv, g_inv)
g = MutableDenseNDimArray(g)
assign(metric, g)
# first derivatives of metric components:
g_d = metric_d.tensor_zeros()
for i in range(n):
    for j in range(i):
        for d in range(n):
            g_d[i, j, d] = g_d[j, i, d]
    for j in range(i, n):
        for d in range(n):
            g_d[i, j, d] = diff(g[i, j], coords[d])
assign(metric_d, g_d)
# Christoffel symbols for Levi-Civita connection (Gam^i_jk):
Gamma = Christoffel.tensor_zeros()
for i in range(n):
    for j in range(n):
        for k in range(j):
            Gamma[i, j, k] = Gamma[i, k, j]
        for k in range(j, n):
            for l in range(n):
                Gamma[i, j, k] += S(1)/2 * g_inv[i, l] * (
                    -g_d[j, k, l] + g_d[k, l, j] + g_d[l, j, k]
                    )
assign(Christoffel, Gamma)
# first derivatives of Christoffel symbols (Gam^i_jk,d):
Gamma_d = Christoffel_d.tensor_zeros()
for i in range(n):
    for j in range(n):
        for k in range(j):
            for d in range(n):
                Gamma_d[i, j, k, d] = Gamma_d[i, k, j, d]
        for k in range(j, n):
            for d in range(n):
                Gamma_d[i, j, k, d] = simplify(diff(Gamma[i, j, k],
                                                    coords[d]))
assign(Christoffel_d, Gamma_d)
# Riemann curvature tensor (R^i_jkl):
Rie = Riemann.tensor_zeros()
for i in range(n):
    for j in range(n):
        for k in range(n):
            for l in range(k):
                Rie[i, j, k, l] = -Rie[i, j, l, k]
            for l in range(k, n):
                Rie[i, j, k, l] = Gamma_d[i, j, l, k] - Gamma_d[i, j, k, l]
                for h in range(n):
                    Rie[i, j, k, l] += (Gamma[h, j, l] * Gamma[i, h, k]
                                    - Gamma[h, j, k] * Gamma[i, h, l])
                    Rie[i, j, k, l] = simplify(Rie[i, j, k, l])
assign(Riemann, Rie)
# Ricci curvature tensor (R_jl):
Ric = simplify(tensorcontraction(Rie, (0, 2)))
assign(Ricci, Ric)
# Ricci curvature scalar:
R = 0
for i in range(n):
    for j in range(n):
        R += g_inv[i, j] * Ric[i, j]
R = simplify(R)
# Einstein tensor (G_ij):
G = Einstein.tensor_zeros()
for i in range(n):
    for j in range(i):
        G[i, j] = G[j, i]
    for j in range(i, n):
        G[i, j] = simplify(Ric[i, j] - S(1)/2 * R * g[i, j])
assign(Einstein, G)
# G^i_j:
G_alt = Einstein_alt.tensor_zeros()
for i in range(n):
    for j in range(n):
        for k in range(n):
            G_alt[i, j] += g_inv[i, k] * G[k, j]
        G_alt[i, j] = simplify(G_alt[i, j])
assign(Einstein_alt, G_alt)

# print it all
print()
metric.print_tensor()
metric_inv.print_tensor()
metric_d.print_tensor()
Christoffel.print_tensor()
Christoffel_d.print_tensor()
Riemann.print_tensor()
Ricci.print_tensor()
if R != 0:
    Idisplay(Math(latex(Eq(Symbol('R'), R))))
    print('\n\n')
Einstein.print_tensor()
Einstein_alt.print_tensor()

EDIT: this code should be executed in Jupyter

Example input:

  • number of dimensions: 4
  • coordinate 0: t
  • coordinate 1: l
  • coordinate 2: theta
  • coordinate 3: phi
  • metric diagonal?: y
  • g_tt: -1
  • g_ll: 1
  • g_thetatheta: r(l)^2
  • g_phiphi: r(l)^2sin(theta)^2
\$\endgroup\$
2
  • \$\begingroup\$ Please replace your output screenshot with a text block. \$\endgroup\$
    – Reinderien
    Commented Sep 6, 2021 at 15:03
  • \$\begingroup\$ Also please indicate how this needs to be executed (i.e. ipython or Jupyter) \$\endgroup\$
    – Reinderien
    Commented Sep 6, 2021 at 15:18

2 Answers 2

4
\$\begingroup\$
  • In general, have your classes on top of the file, and then all the code. Don't have them inbetween.
  • In Python, if we know a variable won't change, and we just have it as a reference or only read from it, we call it a constant. By convention, we use all uppercase names for those: greek -> GREEK.
  • Names should be representative of what they contain. Does greek contain Greeks?!?? Maybe GREEK_CHARACTERS or GREEK_SYMBOLS or GREEK_LETTERS would represent better what's inside it.
  • What happens if the user enters an invalid coordinate? You should account and check that the user input is a valid one.
  • You should use the if __name__ == '__main__' pattern (see more here).
  • IMPORTANT: Divide your code into functions you can reuse. This will avoid code repetition and make your code more modular and easy to understand.

This could become one function def clean_coordinates(coordinates). And same for the rest of the code.

g = eye(n)
while True:
    diag = str(input('Is metric diagonal? y for yes, n for no\n')).lower()
    if diag == 'y':
        for i in range(n):
            g[i, i] = sympify(fix_input(str(input(
                'What is g_[%s%s]?\n' % (str(coords[i]), str(coords[i])
                )))))
    else:
        for i in range(n):
            for j in range(i, n):
                g[i, j] = sympify(fix_input(str(input(
                    'What is g_[%s%s]?\n' % (str(coords[i]), str(coords[j])
                    )))))
                g[j, i] = g[i, j]
    if g.det() == 0:
        print('\nMetric is singular, try again!\n')
        continue
    else:
        break
  • Try to avoid mixing logic with input. E.g., in the example above, you have logic (fixing the input, etc.), but also user input (asking for the diagonal). Instead, have logic functions which handle the logic and take as a parameter whatever they need (e.g. def clean_coordinates(coordinates, diag):) and call them with the user input. This will make your code more modular, testable, clean, organized, reusable, etc.
  • 200 lines of code is a moderately large file... Maybe you can split your code into multiple files which makes it easier to read/understand?

There surely is a lot more, we are just scratching the surface, but I think this is enough for one CR. If you fix all of this and get back, ping me and we can look into more issues!

\$\endgroup\$
8
  • 3
    \$\begingroup\$ "In general, have your classes on top of the file, and then all the code. Don't have them inbetween." If multiple classes are being used, separating all class-definitions into a separate file should be even better. I'm a bit on the fence whether it should already be preferred in this case. Could some of the remaining functions be put into a class for improved stability and maintainability? \$\endgroup\$
    – Mast
    Commented Sep 5, 2021 at 10:08
  • \$\begingroup\$ I agree! I think however first refactoring into functions is key. Then we can start with the file separation if needed, etc.! \$\endgroup\$
    – miquelvir
    Commented Sep 5, 2021 at 10:35
  • 1
    \$\begingroup\$ Overall good feedback. I'm not convinced that 200 lines is long enough to justify file separation for its own sake. Instead, file separation should be done on logical, structural and dependency boundaries, which @Mast hints at. \$\endgroup\$
    – Reinderien
    Commented Sep 6, 2021 at 15:08
  • \$\begingroup\$ True. This is why I said "maybe you can split". Basically now it is a bit hard to figure out, but once divided into functions, we might see some separation become evident! \$\endgroup\$
    – miquelvir
    Commented Sep 6, 2021 at 15:42
  • \$\begingroup\$ @kettle-corn-with-a-q-u please consider marking the question as answered if this is all! if you have any other question let us know \$\endgroup\$
    – miquelvir
    Commented Sep 6, 2021 at 19:52
0
\$\begingroup\$

Interesting project!

Your greek letter escaping and fix_input are unnecessary gymnastics that fight against the norm. You should simply let the user know that they should be using Sympy-compatible format and leave it at that. Replacing ^ with ** is not "fixing" the input; it's presenting an ad-hoc, non-standard, undocumented equation interface to the user that fairly departs from what Python users would expect (^ means xor).

Overall Tensor being a class is a good idea. However, components being a list mutated via a class-external assign method is ungood. Don't accept a list [] at all; declare via components: List[Any] = field(default_factory=list). The Any here is a placeholder because I don't know what actually goes here. But why degrade this to a list - why not keep the NDimArray? Also, assign should be a method on the class.

thing, a, c, x, w, o are pretty foot-directed firearms. They're utterly non-descriptive and will not help us review your code, nor will they likely mean anything to you in 9+ months.

if any(...) is True can just be if any(...). However, this will go away if you abandon the equation conversion effort, which you should.

You need to indicate that this code is to be executed from Jupyter, which you've left as a mystery to your readers. There should be an associated .ipynb that calls into this code. Also consider adding a requirements.txt which is parseable by pip, in which you'd include sympy, notebook and IPython.

\$\endgroup\$
4
  • \$\begingroup\$ Thanks for your time and your input! So I’ve been laboring under a vague idea of eventually making this into a program that is easy to use for everyone regardless of familiarity with Python. Plus I personally was not familiar with Python before a month or two ago (but very familiar with LaTeX), so I still have a strong preference for using carrots over double asterisks. And if I had known that I could use NDimArray as a data type, I absolutely would have done that! So thanks for making my life easier! \$\endgroup\$ Commented Sep 8, 2021 at 4:09
  • \$\begingroup\$ And the nondescript variable names—yeah, I know. I’m a mathematician at heart and new to coding, so it’s taking a lot of effort for me to prioritize descriptiveness and clarity over succinctness. I’m working on it though! \$\endgroup\$ Commented Sep 8, 2021 at 4:10
  • \$\begingroup\$ All very fair. If you want to write a WYSIWYG equation editor frontend consider a MathJax web frontend browser component which can still have Sympy in the backend using something like flask. \$\endgroup\$
    – Reinderien
    Commented Sep 8, 2021 at 11:38
  • \$\begingroup\$ That looks very promising! Thanks for the tip! \$\endgroup\$ Commented Sep 9, 2021 at 2:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.