# Complex scientific formula for nuclear magnetic resonance

I've been reading "clean code" tutorials, and I definitely see the value of using good names that "document themselves" to make intent clear. However, what can be done about complex formulas where describing the meaning of a symbol is difficult?

My solution is to include a reference to the math source, and then use names that relate as closely as possible to the published formula. The code would be unintelligible to another programmer, however, if they didn't have access to the source.

For example, I have a Python program that includes this function factory. The returned function is used elsewhere to plot a graph. The code itself looks like a nightmare:

def dnmr_twosinglet(va, vb, ka, wa, wb, pa):
"""
Accept parameters describing a system of two exchanging, uncoupled nuclei,
and return a function that requires only frequency as an argurment,
that will calculate the intensity of the DNMR lineshape at that frequency.

:param va: The frequency of nucleus 'a' at the slow exchange limit. va > vb
:param vb: The frequency of nucleus 'b' at the slow exchange limit. vb < va
:param ka: The rate constant for state a--> state b
:param wa: The width at half heigh of the signal for nucleus a (at the slow
exchange limit).
:param wb: The width at half heigh of the signal for nucleus b (at the slow
exchange limit).
:param pa: The fraction of the population in state a.
:param pa: fraction of population in state a
wa, wb: peak widths at half height (slow exchange), used to calculate T2s

returns: a function that takes v (x coord or numpy linspace) as an argument
and returns intensity (y).
"""
"""
Formulas for simulating two uncoupled spin-1/2 nuclei are derived from:
Sandstrom, J. "Dynamic NMR Spectroscopy". Academic Press, 1982, p. 15.
"""

pi = np.pi
pi_squared = pi ** 2
T2a = 1 / (pi * wa)
T2b = 1 / (pi * wb)
pb = 1 - pa
tau = pb / ka
dv = va - vb
Dv = (va + vb) / 2
P = tau * (1 / (T2a * T2b) + pi_squared * (dv ** 2)) + (pa / T2a + pb / T2b)
p = 1 + tau * ((pb / T2a) + (pa / T2b))
Q = tau * (- pi * dv * (pa - pb))
R = pi * dv * tau * ((1 / T2b) - (1 / T2a)) + pi * dv * (pa - pb)
r = 2 * pi * (1 + tau * ((1 / T2a) + (1 / T2b)))

def calculate_intensity(v):
"""
:param v: frequency
:return: function that calculates the intensity at v
"""
_Dv = Dv - v
_P = P - tau * 4 * pi_squared * (_Dv ** 2)
_Q = Q + tau * 2 * pi * _Dv
_R = R + _Dv * r
return(_P * p + _Q * _R) / (_P ** 2 + _R ** 2)
return calculate_intensity


However, with the referenced formulae (a subset of which are reproduced below), the meaning of the code should be clear:

$$\textrm{v} = -C_0\frac{\bigg\{P\bigg[1+\tau\Big(\dfrac{p_B}{T_{2A}}+\dfrac{p_B}{T_{2B}}\Big)\bigg]+Q R\bigg\}}{P^2+R^2}$$ $$P=\tau\bigg[\frac{1}{T_{2A} \cdot T_{2B}}-4\pi^2\Delta\nu^2+\pi^2(\delta\nu)^2\bigg]+\frac{p_A}{T_{2A}}+\frac{p_B}{T_{2B}}\\ Q=\tau[2\pi\Delta\nu-\pi\delta\nu(p_A-p_B)]\\ R=2\pi\Delta\nu\bigg[1+\tau\Big(\frac{1}{T_{2A}}+\frac{1}{T_{2B}}\Big)\bigg]+\pi\delta\nu\tau\Big(\frac{1}{T_{2B}}-\frac{1}{T_{2A}}\Big) +\pi\delta\nu(p_A-p_B)$$

etc.

Is this practice considered OK, or could this be improved upon?

• To my eye, aligning more the equals signs (except pi_squared, which is a long name) reduces visual clutter by making it clear where the separation of variable versus assigned value lies. – Richard Apr 11 '17 at 5:35
• Richard: I used to do that, but it violates PEP8. – Geoffrey Sametz Apr 11 '17 at 12:25
• Hmmm. Well it's okay for standards to be wrong, I guess ;-) – Richard Apr 11 '17 at 16:31

I would consider the practice ok in this case.

1. Your variable try appear to be roughly the "ASCII" equivalent of the ones used in the equation.
2. You cite the paper where you got the equations from. (This is probably the most important thing to do) This allows the programmer to understand the context.

Remember readability is subjective, look at something like APL for example and you will see that APL strives for code most closely aligned to math. The creator of this language purposely designed the language to be effective at representing math operations (in particular, matrix manipulation) and it takes a vastly different approach from Python.

One interesting thought experimenet: If we grew up and learned that the word "print" was simply "p", "print" would likely be considered not readable and we would replace each case with "p". However, in our world, we consider the opposite to be true in most cases.

For me it is much easier to read dx instead of something like change_in_x because I have learned to appreciate the former.

The only thing I would consider changing is getting rid of pi_squared and just writting it out explicitly (unless there is some worthwhile performance cost in doing so).

• So I can only upvote this once, and probably cannot improve on it, so I will leave a comment emphasizing that is absolutely correct. Variable names should be clear and descriptive. When it comes to equations, things like dx are the lingua franca of the problem domain, and thus are clear and descriptive. Leaving the equation in a form that is readable (latex?) or a pointer to the paper in question, is the best (and the least) you can do to communicate the problem domain to the next maintainer. – Stephen Rauch Apr 11 '17 at 3:09
• Dair: I put pi_squared in because I wanted to pre-compute everything possible before creating the inner function calculate_intensity with those results. I took this approach over a simple function because then calculate_intensity is applied along an array of hundreds of x coordinates, and I expected that this would speed up the calculations. Speed tests indicated it helped, although it's not the performance bottleneck. – Geoffrey Sametz Apr 11 '17 at 12:33

There's no problem with using short names, as long as they can be clearly understood in context. However, I think that some of the names in the code under review are not clearly understood in context. Comparing the code and the equation for Q or R I see missing terms. I think that the equation for Q actually corresponds to _Q from the inner function, and similarly that R actually corresponds to _R. In my opinion it would improve clarity to give the "outer" variables names like Q_0 and R_0, and probably also add a comment along the lines of

# The constant terms of R and Q


To improve the readability of formulae it is worth considering which common subexpressions to factor out. Every instance of dv which I see is multiplied by pi, and I don't see a numerical-analytical reason not to factor out a pi_dv variable. Similarly for pi_Dv. That would then automatically address the point raised in Dair's review about pi_squared, because instead of pi_squared * (dv ** 2) you would have pi_dv ** 2, and similarly for _Dv.

• There was some chopping done, because the v-dependent terms were separated from the v-independent terms (as you saw) so that the v-independent terms were pre-calculated outside the inner function. Comments and maybe a docstring for the inner function might help. Good point about factoring pi out. – Geoffrey Sametz Apr 11 '17 at 12:45