WIXLIB

Appendix A: Object-Oriented Programming in a Nutshell231In the for loop, we make use of the class method two body step() to updatepositions and velocities using the symplectic Euler scheme. In contrast to an instancemethod, it does not act on a particular object. In the above listing of the class, you cansee that the definition of two body step() is marked by the @classmethoddecorator. The argument named cls is analogous to self, except that it refers tothe class name instead of the current object. The two objects for which we want tocompute the Euler step are passed just like arguments of an ordinary function. But,hold on, did we not state that functions only receive input through their arguments,without changing them? Well, two body step() can change the states of thetwo bodies because Python uses a mechanism that is known as call by object reference. You can easily convince yourself that theses changes persist after callingtwo body step().Like functions, class methods can change public as well as protected attributesof their arguments, provided that they are mutable objects (call by objectreference).The gravitational force between the two bodies is computed with the instance methodgravity(). Inside two body step(), it is called in the same fashion asdistance() in the example above (see line 34), with cls pointing to the class.As a result, the states of body1 and body2 are iterated in the loop through all timesteps and successive positions are recorded, like ephemerides, in the arrays orbit1and orbit2. It is left as an exercise to plot the orbits and to compare them withFig. 4.8.Why not split two body step() into separate updates for the two bodies,which could be implemented as instance methods? In that case, the gravitationalforce between the bodies would be evaluated twice rather than applying Newton’sthird law to compute the their accelerations at once. Since force computation involvesthe most expensive arithmetic operations, this should be avoided. You can try to writea version for three-body interactions and apply it, for instance, to the Algol system(see Sect. 4.3). In systems consisting of many gravitating bodies, such as the stellarcluster discussed in Sect. 4.5, performance becomes the main issue. Direct summation of gravitational forces over all pairs of bodies becomes quickly intractable andapproximate algorithms have to be applied. If you want to test it, see Exercise 4.18.Although defining objects for pieces of data helps to write well structured code,the overhead can become problematic in numerical applications. As a consequence, itis preferable to collect data in objects rather than representing each item individuallyby an object. You can do an OOP project of your own in Exercise 4.15. The task is towrite a class for a large number of test particles orbiting a gravitating central mass.

Appendix BMaking Python FasterB.1Using ArraysThe native data structure for ordered collections of data items in Python is a list.At first glance, lists are very similar to NumPy arrays. For example, the days ofequinoxes and solstices in 2020 could be defined as a list:1N [79, 171, 265, 355]Compare this to the definition of the array N in Sect. 2.1. The list on the right-handside is just what is passed as argument to np.array(). What does this functiondo? Well, it converts a given list into an array. If you enter2N[1], type(N[1]), type(N)you will see the output(171, int, list)So N[1] is an integer, which is the second element of the list N (indexing works inthe same way as with arrays). Now let us convert this list into an array and assignedit to the same name as before. The original list will be deleted in this case:345678import numpy as np# convert to arrayN np.array(N)print(N[1], N[1].dtype, N.dtype)Now the data type of the element N[1] is given by the NumPy attribute dtype:171 int64 int64 Springer Nature Switzerland AG 2021W. Schmidt and M. Völschow, Numerical Python in Astronomy and Astrophysics,Undergraduate Lecture Notes in 33

234Appendix B: Making Python FasterIn contrast to a list, this data type is also an attribute of the whole array because thedata type must be uniform, i.e. all elements must have an identical data type.This does not apply to lists:9101112# redefine listN [79, "summer solstice", 265, "winter solstice"]N[1], type(N[1]), type(N)Now some of the elements are strings in place of integers:(’summer solstice’, str, list)You can easily convince yourself that the first and third elements are still of typeint. Although the meaning of the expressions in line 10 might be the same for ahuman, the meaning that pertains to these data types is fundamentally different inPython. What happens if N is again converted into an array? Try for yourself. Youare probably up for a surprise.Uniformity of the data type is important for the memory layout of NumPy arrays.The speed of modern CPUs is fast compared to the time required to load data fromthe computer’s main memory (RAM). If data cannot be loaded more or less in bigchunks, the CPU will inevitably be idle for a while, before it receives the next piece ofdata. This situation is typically encountered when working with Python lists becausedifferent elements of the list can be stored at random positions in memory. NumPyarrays, on the other hand, alleviate the problem by arranging elements consecutivelyin memory.2 The total memory required to store an array is given by the numberof elements times the size of each element in bytes and there is a simple mappingbetween array elements and their positions in memory.3It turns out that this has a substantial impact on the efficiency of numerical computations using NumPy functions or array operations. As an example, let us computethe Planck spectrum for a given temperature (see Sect. 3.1.2). First, we do this byusing lists. For convenience, we use physical constants from scipy.constants.12345678import mathfrom scipy.constants import h,c,k,sigma# list of wavenumbersn 1000lambda max 2e-6lambda step lambda max/nwavelength [i*lambda step for i in range(1,n 1)]2 This is possiblethrough a hierarchy of caches which provide access to blocks of memory at muchhigher speed than the main memory.3 This is also true for multi-dimensional arrays, although the mapping is slightly more complicated.

Appendix B: Making Python Faster235The last line shows how to iteratively create a Python list of uniformly spacedwavenumbers by means of an inline for loop. The following function calculates the intensity for a list of wavelengths. The effective temperature of the Sun isset as default.910111213141516171819def planck spectrum(wavelength, T 5778):# create empty listspectrum []# loop through wavelengths and append flux valuesfor val in wavelength:spectrum.append(2*h*c**2 /(val**5 * (math.exp(min(700, h*c/(val*k*T))) - 1)))return spectrumHere, we have an explict loop that runs through all wavelengths and computes thecorresponding intensity using the exponential function from the math module. Theresult is appended to the list spectrum, which is initialized as an empty list withoutany elements in line 12. This is an important difference between lists and arrays.The latter always have to be initialized with non-zero length. Moreover, the methodappend() in the example above modifies an existing object. As a result, appendedelements can be located anywhere in memory. In contrast, NumPy’s append()function returns a newly allocated array that is a copy of the original array plus oneor more new elements.A tool for performance measurement that is particularly easy to use is the magiccommand timeit in interactive Python (alternatively use -m timeit as command line option if Python is executed on the shell):20%timeit planck spectrum(wavelength)which outputs (numbers will differ depending on your computer architecture):834 ms 3.44 ms per loop (mean std. dev. of 7 runs, 1000 loops each)This tells us that execution of planck spectrum(wavelength) took 834 ms(a little less than one second) averaged over 7 runs. Since timeit is tailored tosmall code snippets that are executed in a very short time, the measurement periodis artificially increased by a large number of loops for higher precision (the numberof loops is automatically adjusted).It is left as an exercise for you to time the implementation ofplanck spectrum() based on NumPy from Sect. 3.1.2 (do not forget to convertwavelength to an array). We get:50.3 ms 1.28 ms per loop (mean std. dev. of 7 runs, 10000 loops each)With NumPy, the computation is more than ten times faster! We can assign the arrayreturned by planck spectrum() to explore its properties:

2362122Appendix B: Making Python Fastersolar planck spectrum(wavelength)solar.flagsThe flags attribute gives us the following information (abridged):C CONTIGUOUS : TrueF CONTIGUOUS : TrueIndeed, the resulting array is contiguous in memory.4 You may also check the arraywavelength.However, this is not the only reason for higher efficiency. Array operations alsoavoid function calls inside loops, which are particularly costly. Especially beginnersmight prefer an explicit loop as it is easy to read. For example, we can apply thetrapezoidal rule to integrate the Planck spectrum (see Sect. 3.2.2) using a for loop:123456789101112131415def integr trapez(f, a, b, n):# integration steph (b - a)/n# initialisationtmp 0.5*f(a)# loop through subintervals between a h and b-hfor i in range(1,n):tmp f(a i*h)tmp 0.5*f(b)return h*tmpTiming the integration of the Planck spectrum,16%timeit integr trapez(planck spectrum, 1e-9, 364.7e-9, 100)yields a mean execution time of474 µs 13.4 µs per loop (mean std. dev. of 7 runs, 1000 loops each)For comparability, we used the Numpy version of planck spectrum(). Nowlet us see what we get if the integrator from Sect. 3.2.2 is used:76.7 µs 893 ns per loop (mean std. dev. of 7 runs, 10000 loops each)Once more, the speed up is palpable. Passing an array as function argument andapplying np.sum() turns out to be much faster than calling the function for singlevalues inside the body of a for loop, as in the code example above. So the lesson4 The two types are relevantfor two-dimensional arrays, where the layout can be row-major as in Cor column-major as in Fortran.

Appendix B: Making Python Faster237learned is to avoid explicit loops and especially function calls inside loops wheneverpossible. Even in cases where this is not feasible, there are means of improvingperformance, as shown in the following section.B.2Cythonizing CodeAs a case study, we consider the Strömgren sphere discussed in Sect. 4.1.1. Tosolve the initial value problem with the fourth-order Runge-Kutta method (RK4), weiteratively call the Python function rk4 step(). As demonstrated in the previoussection, function calls inside a loop have a negative impact on performance. We cantime the complete RK4 integration by defining a wrapper function:123456789101112import numpy as npdef solve stroemgren(r0, dt, n steps):t np.linspace(0, n steps*dt, n steps 1)r np.zeros(n steps 1)r[0] r0for n in range(n steps):r[n 1] rk4 step(lambda t, r: (1 - r**3)/(3*r**2),t[n], r[n], dt)return (t,r)This allows us to utilize the %timeit command introduced in Sect. B.1:13%timeit solve stroemgren(0.01, 1e-3, 10000)The measured execution time (subject to the computer in use) is94.3 ms 997 µs per loop (mean std. dev. of 7 runs, 10 loops each)which is about a tenth of a second.If this is so, why do we not implement the RK4 scheme for a specific differentialequation rather than using a function for arbitrary initial value problems? Withouta generic, reusable implementation, it would be necessary to rewrite the code everytime a new problem is to be solved. This is time consuming and prone to errors.A better solution is to translate the Python code for the RK4 scheme into thecompilable language C. A compiler produces machine code that is much faster thancode that is executed by the Python interpreter. This is made possible by Cython andrequires some preparation.5 To begin with, we need to put the functions we want toturn into C code in a module with file extension .pyx. In our example, this meanscopying the definition of rk4 step() into a file named stroemgren.pyx (or5 Formore information, see cython.readthedocs.io/en/latest/index.html.

238Appendix B: Making Python Fasterextract the file from the zip-archive for this section). An important modification totake advantage of Cython is static typing. Remember that variables in Python areversatile objects without a fixed data type. For example, you can initially assign afloat to a variable, then change it to a string and later to something entirely different.As a consquence, there is no way of telling what kind of data is passed as actualargument to a function before it is called. This is nice and extremely flexible, butnot favourable in terms of efficiency. In a low-level language such as C, you need toexactly specify the type of every variable, function argument, and return value. Soif you want to translate Python code into C via Cython, you should make data typesexplicit in your Python source code:123456789# excerpt from stroemgren.pyxcpdef double crk4 step(f, double t, double x, double k4 dtdtdtdt****f(t, x)f(t 0.5*dt, x 0.5*k1)f(t 0.5*dt, x 0.5*k2)f(t dt, x k3)return x (k1 2*(k2 k3) k4)/6.0The C type double corresponds to a floating point number in Python. For localvariables you need to use the Cython keyword cdef followed by the type. Inaddition, the function will return a static type if it is declared with cpdef instead ofdef. Formal arguments of a function are simply prefixed with the type. You mightnotice that first argument has no type because f is the name of a Python function.As preparation for using the module stroemgren.pyx, we create a smallPython script setup.py:from setuptools import setupfrom Cython.Build import cythonizesetup(ext modules cythonize("stroemgren.pyx"))The central directive is cythonize("stroemgren.pyx"). It instructs Pythonto cynthonize the code contained in the module stroemgren.pyx. For the nextstep, you need a C compiler, such as the GNU compiler gcc, on your system.6 Thecommandpython setup.py build ext --inplace6 OnLinux and Mac systems you can try to type gcc -v on the command line. If the compiler isinstalled, you will get some information, otherwise an error message will result. Anaconda comeswith integrated compiler tools. If there is no compiler on your computer, search the web for availablecompilers for your system.

Appendix B: Making Python Faster239executed on the command line instructs Cython to produce a shared object file (alsoknow as shared library).7 Shared libraries are extensions that can be used by programswhen they are executed.Now we can import and use crk4 step() just like any other Python functiondefined in a module:12345678910111213import numpy as npfrom stroemgren import crk4 stepdef solve stroemgren(r0, dt, n steps):t np.linspace(0, n steps*dt, n steps 1)r np.zeros(n steps 1)r[0] r0for n in range(n steps):r[n 1] crk4 step(lambda t, r: (1 - r**3)/(3*r**2),t[n], r[n], dt)return (t,r)A measurement of the execution time confirms that using crk4 step() insteadof rk4 step() makes a difference:14%timeit solve stroemgren(0.01, 1e-3, 10000)The measured execution time (subject to the computer in use) is21.5 ms 1.14 ms per loop (mean std. dev. of 7 runs, 10 loops each)The gain in speed is about a factor four.It turns out that we can do even better if we refrain from passing a Python functionas argument (something that does not translate well into C). In the example above, aPython lambda is used to specify the derivate dr̃ /dt defined by Eq. (4.11). A moreefficient computation can be accomplished by using a pure C function in the modulestroemgren.pyx:151617# excerpt from stroemgren.pyxcdef double rdot(double t, double r):return (1.0 - r**3)/(3.0*r**2)The Cython keyword cdef indicates that this function will be compiled as a Cfunction and cannot be called directly by a Python program. For this reason, it hasto be referenced explicitly from within the RK4 integrator:7 Ifyou list the files in your work directory, you will see that also a file named stroemgren.c isproduced. This file contains the C code generated by Cython from which the C compiler producesthe shared object file. It should be mentioned that the intermediate C code is not intended to be readby humans.

240181920212223242526Appendix B: Making Python Faster# excerpt from stroemgren.pyxcpdef double stroemgren step(double t, double r, double k4 dtdtdtdt****rdot(t, r)rdot(t 0.5*dt, r 0.5*k1)rdot(t 0.5*dt, r 0.5*k2)rdot(t dt, r k3)return x (k1 2*(k2 k3) k4)/6.0Of course, we give away versatility since the integrator stroemgren step()requires a particular function name. However, we could easily change its definition tosolve other first-order differential equations. Let us see how well the fully cythonizedfunctions perform in the following program:1234567891011from stroemgren import stroemgren stepdef solve stroemgren(r0, dt, n steps):t np.linspace(0, n steps*dt, n steps 1)r np.zeros(n steps 1)r[0] r0for n in range(n steps):r[n 1] stroemgren step(t[n], r[n], dt)return (t,r)Executing12%timeit solve stroemgren(0.01, 1e-3, 10000)we finally get7.97 ms 260 µs per loop (mean std. dev. of 7 runs, 100 loops each)Compared to the plain Python version, we have achieved a speed-up by more thana factor of ten. You might still think that is of no concern whether your Pythonprogram finishes in a somewhat shorter or longer fraction of a second. However, ifyou increase the number of steps to 100.000—the convergence study in Sect. 4.1.1suggests that this is the recommended number of steps—you will begin to noticethe difference without using timeit. If it comes to more demanding problems,for example, the colliding disk simulations in Sect. 4.1.1, performance optimizationbecomes an issue, and even more so in scientific data analysis in research projects.You need to find a balance between the versatility offered by Python and theefficiency of languages such as C and Fortran that serves your purpose. There ismuch more to Cython than what we can cover here. We have given you only a firsttaste and your are invited to explore the capabilities on your own.8 Since Cython8Agood starting point are the online tutorials: cython.readthedocs.io/en/latest/src/tutorial.

Ap

Object-Oriented Programming in a Nutshell Classical programming languages, such as Fortran or C, are procedural languages: . 226 Appendix A: Object-Oriented Programming in a Nutshell An essential part of a Python class is the definition of the _init_method, which is also called constructor. It defines the attributes of an object belonging to