Python3: How do I make this equation readable?

Question

So for this internship I am doing, I need to use the integral of (x^(9/2))/((1-x)^2) as part of an equation I am graphing. However, the variable that I am graphing along the x axis appears in both limits of integration. Since I am a complete and total novice at python, my code is atrocious, but I ended up copy+pasting the indefinite integral twice, plugging in the limits of integration, and subtracting. How can I make the code better?

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

x = np.arange(0,2.5,0.00001)

zs = 8
zr = 6
rbub = 5
sig = 2.5
XHI = 0.5
sigma = 10**-sig
L = 10**-3
zb = zs - (((1 + zs)/8)**(3/2)) * (0.0275) * (rbub/10)
a = (1+zr)/((1+zs)*(x+1))
b = (1+zb)/((1+zs)*(x+1))

def f(x):
    ans = 0.000140092
    ans = ans * ((1+zs)**(3/2))
    ans = ans * ((x+1)**(3/2))
    ans = ans * XHI
    return ans * ((9/2)*(np.log(1-np.sqrt(b)) - np.log(np.sqrt(b)+1)) + (1/35 * (1/(b-1)) * (10*(b**(9/2)) + 18*(b**(7/2)) + 42*(b**(5/2)) + 210*(b**(3/2)) - 315*(b**(1/2))) - ((9/2)*(np.log(1-np.sqrt(a)) - np.log(np.sqrt(a)+1)) + (1/35 * (1/(a-1)) * (10*(a**(9/2)) + 18*(a**(7/2)) + 42*(a**(5/2)) + 210*(a**(3/2)) - 315*(a**(1/2)))))))

Question on how to improve already working code should be directed to https://codereview.stackexchange.com/ — bphi, May 29 '18 at 17:03
Write a separate assignment statement that evaluates the definite integral and stores it under a name, like `defintval = `. You probably also want to break it down further into smaller steps, so you can use a bit more whitespace so your expression doesn't scroll across multiple screens. — abarnert, May 29 '18 at 17:05
You could try to pit your variables in two seprate lists and iterate on them using a ZIP(). A second thought is to put them in a one list of tuples that have to elements. — Dawid Dave Kosiński, May 29 '18 at 17:06
@CtrlS That directly violates the advice in PEP 8. First, mixing tabs and spaces—even in syntactically insignificant whitespace—is a recipe for disaster in Python; all kinds of tools will detect your tabs and get confused. Second, if you such a huge chain of assignments that you feel the need to line them that way, you probably need to refactor things. Of course PEP 8 is just guidelines, that are sometimes worth violating, but starting off with the idea of violating them without knowing what they're there for doesn't usually lead to good design. — abarnert, May 29 '18 at 17:17
Welcome to Stack Overflow! If you believe that the code works correctly, consider presenting your work (with its unit tests) in a more-complete fashion over at [codereview.se]. You'll likely get some suggestions on making it more efficient, easier to read, and better tested. Before you do that, make sure to read [A guide to Code Review for Stack Overflow users](//codereview.meta.stackexchange.com/a/5778) first, as some things are done differently over there - e.g. question titles should simply say what the code *does*, as the question is always, "How can I improve this?". — Toby Speight, May 29 '18 at 17:32
@CtrlS I've never read Code Complete is, but if you want to argue from authority, it's hard to beat PEP8 as the authority. Besides being an official quasi-standard, it's also the way the standard library is built, and it underlies the rules that every tool from style checkers to formatters to IDEs all follow (although they all extend on it in some ways). — abarnert, May 29 '18 at 17:52
@CtrlS Meanwhile, if you want to ask about why the Python community prefers tabs over spaces, and get references to what people have written, etc., the comments in an unrelated SO question are not the place to do that. But the reality is, you should be using tools that make the question irrelevant. If you can't configure your editor/IDE to reindent code the way you personally prefer on load and reindent back to 4 spaces on save, get a better editor. — abarnert, May 29 '18 at 17:57

score 2 · Answer 1 · answered May 29 '18 at 17:48

Here is one take. Of course, I have no idea what this is doing. You should be in a much better position to add comments / sensible variable names, as others have pointed out. Still, here is what you could do.

First, run a code formatter to make the code more human-friendly.

def f(x):

    ans = 0.000140092
    ans = ans * ((1 + zs) ** (3 / 2))
    ans = ans * ((x + 1) ** (3 / 2))
    ans = ans * XHI

    return ans * (
        (9 / 2) * (np.log(1 - np.sqrt(b)) - np.log(np.sqrt(b) + 1))
        + (
            1
            / 35
            * (1 / (b - 1))
            * (
                10 * (b ** (9 / 2))
                + 18 * (b ** (7 / 2))
                + 42 * (b ** (5 / 2))
                + 210 * (b ** (3 / 2))
                - 315 * (b ** (1 / 2))
            )
            - (
                (9 / 2) * (np.log(1 - np.sqrt(a)) - np.log(np.sqrt(a) + 1))
                + (
                    1
                    / 35
                    * (1 / (a - 1))
                    * (
                        10 * (a ** (9 / 2))
                        + 18 * (a ** (7 / 2))
                        + 42 * (a ** (5 / 2))
                        + 210 * (a ** (3 / 2))
                        - 315 * (a ** (1 / 2))
                    )
                )
            )
        )
    )

Right away you see some symemtry. This chunk

10 * (b ** (9 / 2))
+ 18 * (b ** (7 / 2))
+ 42 * (b ** (5 / 2))
+ 210 * (b ** (3 / 2))
- 315 * (b ** (1 / 2))

is a dot product of some weights and a b raised to a vector of powers. If b were scalar, we could write it as np.dot(weights, np.sqrt(b) ** powers). Maybe we would even score some optimization points from using integral powers.

Putting thigs together, we can get something like this:

weights = np.array([10, 18, 42, 210, -315])
powers = np.array([9, 7, 5, 3, 1])


def log_term(x):
    return (9 / 2) * (np.log(1 - np.sqrt(x)) - np.log(np.sqrt(x) + 1))


def dot_term(x):
    return (1 / 35) * 1 / (x - 1) * np.dot(np.sqrt(x)[..., None] ** powers, weights)


def integrate(x):
    return log_term(x) + dot_term(x)


factor1 = integrate(b) - integrate(a)
factor2 = 0.000140092 * ((1 + zs) ** (3 / 2)) * XHI
factor = factor1 * factor2


def f(x):
    return factor * ((x + 1) ** (3 / 2))

With better variable names and comments this could almost be readable.

Side comment. Both in your original code and in this version, you define x in the body of your script. You also define several variables as a function of x, such as a and b.

Python scoping rules mean that these variables will not change if you pass a different x to f. If you want all of your variables to change with x, you should move the definitions inside the function.

score -2 · Answer 2 · answered May 29 '18 at 17:09

-2

Using good variable names will help a lot for who gonna read the code or even for you, and comments will help too. For the rest, it is a equation, there is no good way of putting it.

answered May 29 '18 at 17:09

Joao Pedro Lopes Mendes

66
10

Downvoted because I disagree with "no good way of putting it". – hilberts_drinking_problem May 29 '18 at 17:15
A equation will be always a equation, probably could be more pythonic but well, i probably did put the words way to hard saying "there is no good way of putting it". :)) – Joao Pedro Lopes Mendes May 29 '18 at 17:16
Good variable names etc. are indeed good practice, but this question is about how to make the equation more readable and that is definitely possible by a variety of means. Otherwise, I would have given a similar answer. – Ctrl S May 29 '18 at 17:27

Python3: How do I make this equation readable?

2 Answers2