Picture by me of the blood moon!

This quarter I’m taking a physics/programming class where we learn numerical methods in the context of physics(high energy more or less) with python. We’ve covered some interesting tools and numerical approaches such as Markov Chain Monte Carlo, Bayesian inference etc.

Naturally, we would use numpy in python extensively since the vectorized nature makes it very easy to manipulate arrays. To people who’re less familiar, this is what I mean:

import numpy as np
xs = np.linspace(1,10,10) #[1,2,3,...,10]
ys = np.power(x,2) #[1,4,9....,100]

The problem occurs when a numpy function can handle both number and array (also overide*). Continuing the example above, now image I have a function:

def integrate(u):
    ary = np.exp(-ys) + u # this is array because ys is array
    return np.sum(ary) # this will sum the above array
    
s = [integrate(x) for x in xs]
s = integrate(xs) # these two will yield different result

This is because when passing xs to integrate, python has no way of knowing if I wans it to treat it as int or array. Of course, this is not too nasty because expand list comprehension is easy. But it took me too long to realize this is the bug.

In contrast, in julia I could have made clear in the definiton of the integrate(u) for example:

function integrate(u::Number)
    ...
end

Initially I tried to google how to force typing in python and then realize that’s the opposite of python philosophy. This is also reflected in the ugly __mul__(self, other) magic function in making a class in python. In that case, you would have to write multiple if isinstance(other, <something>): clauses to make things workd, where in julia I could have the blessed super/sub type to figure out what method to use when executing the function.