GD
Geodynamics

Archives / 2018 / October

Inversion 101 through 201 – Part 2: The inverse problem and deterministic inversion

Inversion 101 through 201 – Part 2: The inverse problem and deterministic inversion

The Geodynamics 101 series serves to showcase the diversity of research topics and methods in the geodynamics community in an understandable manner. We welcome all researchers – PhD students to professors – to introduce their area of expertise in a lighthearted, entertaining manner and touch upon some of the outstanding questions and problems related to their fields. This time, Lars Gebraad, PhD student in seismology at ETH Zürich, Switzerland, talks about inversion and how it can be used in geodynamics! Due to the ambitious scope of this topic, we have a 3-part series on this! You can find the first post here. This week, we have our second post in the series, where Lars discusses the inverse problem and deterministic inversion.

The idea of inversion is to literally invert the forward problem. For the aforementioned problems our knowns become unknowns, and vice versa. We wish to infer physical system parameters from measurements. Note that in our forward problem there are multiple classes of knowns; we have the forcing parameters, the material properties and the boundary conditions. All of these parameters could be inverted for in an inversion, but it is typically one only inverts for one class. Most of the time, medium parameters are the target of inversions. As we will go through the examples, I will gradually introduce some inversion methods. Every new method is marked in blue.

Estimating ingredients from a bread

Now, let’s consider our first example: the recipe for a bread. Let’s say we have a 0.5 kilogram bread. For our first case (case 1a), we will assume that the amount of water is ideal. In this case, we have one free variable to estimate (the used amount of flour), from one observable (the resulting amount of bread). We have an analytical relationship that is both invertible and one-to-one. Therefore, we can use the direct inversion of the relationship to compute the amount of flour. The process would go like this:

Direct inversion

1. Obtain the forward expression G(m) = d;
2. Solve this formula for G^-1 (d) = m;
3. Find m by plugging d into the inverse function.

Applying this direct inversion shows us that 312.5 grams of flour must have been used.

The two properties of the analytical relationship (invertible and one-to one) are very important for our choice of this inversion method. If our relationship was sufficiently complex, we couldn’t solve the formula analytically (though the limitation might then lie by the implementer 😉 ). If a function is not one-to-one, then two different inputs could have the same output, so we cannot analytically invert for a unique solution.

In large linear forward models (as often obtained in finite difference and finite element implementations), we have a matrix formula similar to

A m = d.

In these systems, m is a vector of model parameters, and d a vector of observables. If the matrix is invertible, we can directly obtain the model parameters from our observables. The invertibility of the matrix is a very important concept, which I will come back to in a bit.

Problems in direct inversion, misfits and gradient descents

Let’s see what happens when the first condition is not satisfied. I tried creating a non-linear function that at least I could not invert. Someone might be able to do it, but that is not the aim of this exercise. This new relationship is given here:

Complex flour to bread relationship, another 1D forward model.

The actual forward model is included in the script at the bottom of the post. When we cannot analytically solve the forward model, one approach would be to iteratively optimise the input parameters such that the end result becomes as close as possible to the observed real world. But, before we talk about how to propose inputs, we need a way to compare the outputs of our model to the observables!
When we talk about our bread, we see one end-product. It is very easy to compare the result of our simulation to the real world. But what happens when we make a few breads, or, similarly, when we have a lot of temperature observations in our volcanic model? One (common) way to combine observables and predictions into a single measure of discrepancy is by defining a misfit formula which takes into account all observations and predictions, and returns a single – you guessed it – misfit.
The choice of misfit directly influences which data is deemed more important in an inversion scheme, and many choices exist. One of the most intuitive is the L2-norm (or L2 misfit), which calculates the vector distance between the predicted data (from our forward model) with the observed data as if in Euclidean space. For our bread it would simply be

X = | predicted bread – observed bread |

Let’s try to create a misfit function for our 1D bread model. I use the relationship in the included scripts. Again, we have 500 grams of bread. By calculating the quantity | G(m) – 500 |, we can now make a graph of how the misfit varies as we change the amount of flour. I created a figure which shows the observed and predicted value:

Complex flour to bread relationship with observation superimposed

and a figure which shows the actual misfit at each value of m

Misfit between model and observation

Looking at the two previous figures may result in some questions and insights to the reader. First and foremost, it should be obvious that this problem does not have a unique solution. Different amounts of flour give exactly the same end result! It is thus impossible to say with certainty how much flour was used.
We have also recast our inversion problem as an optimisation function. Instead of thinking about fitting observations we can now think of our problem as a minimization of some function X (the misfit). This is very important in inversion.
Iterative optimizations schemes such as gradient descent methods (and the plethora of methods derived from it) work in the following fashion:

Gradient descent inversion/optimization

1. Pick a starting point based on the best prior knowledge (e.g., a 500 gram bread could have been logically baked using 500 grams of flour);
2. Calculate the misfit (X) at our initial guess;
3. If X is not low enough (compared to some arbitrary criterion):
• Compute the gradient of X;
• Do a step in the direction of the steepest gradient;
• Recompute X;
• Repeat from 2.
4. If X is low enough:
• Finalise optimisation, with the last point as solution.
These methods are heavily influenced by the starting point of the inversion scheme. If I would have started on the left side of the recipe-domain (e.g. 100 grams of flour), I might well have ended up in a different solution. Gradient descents often get ‘stuck’ in local solutions, which might not even be the optimal one! We will revisit this non-uniqueness in the 2D problem, and give some strategies to mitigate creating more than one solution. Extra material can be found here.

Simplistic gradient descent

Obvious solutions, grid searches, and higher dimensions

One thing that often bugged me about the aforementioned gradient descent methods is the seemingly complicated approach for such simple problems. Anyone looking at the figures could have said

Well, duh, there’s 3 solutions, here, here and here!

Why care about such an complicated way to only get one of them?
The important realisation to make here is that I have precomputed all possible solution for this forward model in the 0 – 700 grams range. This precomputation on a 1D domain was very simple; at a regular interval, compute the predicted value of baked bread. Following this, I could have also programmed my Python routine to extract all the values with a sufficiently low misfit as solutions. This is the basis of a grid search.
Let’s perform a grid search on our second model (1b). Let’s find all predictions with 500 grams of bread as the end result, plus-minus 50 grams. This is the result:

Grid search for the 2D model

The original 312.5 grams of flour as input is part of the solution set. However, the model actually has infinitely many solutions (extending beyond the range of the current grid search)! The reason that a grid search might not be effective is the inherent computational burden. When the forward model is sufficiently expensive in numerical simulation, exploring a model space completely with adequate resolution might take very long. This burden increases with model dimension; if more model parameters are present, the relevant model space to irrelevant model space becomes very small. This is known as the curse of dimensionality (very well explained in Tarantola’s textbook).
Another reason one might want to avoid grid searches is our inability to appropriately process the results. Performing a 5 dimensional or higher grid searches is sometimes possible on computational clusters, but visualizing and interpreting the resulting data is very hard for humans. This is partly why many supercomputing centers have in-house departments for data visualisation, as it is a very involved task to visualise complex data well.

Now: towards solving our physical inversions!

Non-uniqueness, regularization and linear algebra (bleh!)

One big problem in inversion is non-uniqueness: the same result can be obtained from different inputs. The go-to way to combat this is to add extra information of any form to the forward problem. In our bread recipe we could think of adding extra observables to our problem, such as the consistency of the bread, its taste, color, etc. Another option could be to add constraints on the parameters, such as using the minimum amount of ingredients. This is akin to asking the question: given this amount of bread, how much flour and water was minimally used to make it?
Diffusion type problems are notorious for their non-uniqueness. Many different subsurface heat conduction distributions might result in the observations (imagine differently inclined volcanic conduits). An often used method of regularisation (not limited to diffusion type studies!) is spatial smoothing. This method requires that among equally likely solutions, the smoothest solutions are favoured, for it is more ‘natural’ to have smoothly varying parameters. Of course, in many geoscientific settings one would definitely expect sharp contrasts. However, in ‘underdetermined’ problems (i.e., you do not have enough observations to constrain a unique solution), we favour Occam’s Razor and say

The simplest solution must be assumed

When dealing with more parameters than observables (non-uniqueness) in linear models it is interesting to regard the forward problem again. If one would parameterize our volcanic model using 9 parameters for the subsurface and combine that with the 3 measurements at the surface, the result would be an underdetermined inverse problem.

Rough parametrization for the heat conduction model

This forward model (the Laplace equation) can be discretised by using, for example, finite differences. The resulting matrix equation would be Am = d, with A a 3 x 9 matrix, m a 9 dimensional vector and d a 3 dimensional vector. As one might recall from linear algebra classes, for a matrix to have an inverse, it has to be square. This matrix system is not square, and therefore not invertible!

Aaaaahhh! But don’t panic: there is a solution

By adding either prior information on the parameters, smoothing, or extra datapoints (e.g., taking extra measurements in wells) we can make the 3 x 9 system a perfect 9 x 9 system. By doing this, we condition our system such that it is invertible. However, many times we end up overdetermining our system which could result in a 20 x 9 system, for example. Note that although neither the underdetermined nor the overdetermined systems have an exact matrix inverse, both do have pseudo-inverses. For underdetermined systems, I have not found these to be particularly helpful (but some geophysicists do consider them). Overdetermined matrix systems on the other hand have a very interesting pseudo-inverse: the least squares solution. Finding the least squares solution in linear problems is the same as minimising the L2 norm! Here, two views on inversion come together: solving a specific matrix equation is the same as minimising some objective functional (at least in the linear case). Other concepts from linear algebra play important roles in linear and weakly non-linear inversions. For example, matrix decompositions offer information on how a system is probed with available data, and may provide insights on experimental geophysical survey design to optimise coverage (see “Theory of Model-Based Geophysical Survey and Experimental Design Part A – Linear Problems” by Andrew Curtis).
I would say it is common practice for many geophysicists to pose an inverse problem that is typically underdetermined, and keep adding regularization until the problem is solvable in terms of matrix inversions. I do not necessarily advocate such an approach, but it has its advantages towards more agnostic approaches, as we will see in the post on probabilistic inversion next week!


Summary of deterministic inversions

We’ve seen how the forward model determines our inversion problem, and how many measurements can be combined into a single measure of fit (the misfit). Up to now, three inversion strategies have been introduced:

Direct inversion: analytically find a solution to the forward problem. This method is limited to very specific simple cases, but of course yields near perfect results.
Gradient descent methods: a very widely used class of algorithms that iteratively update solutions based on derivatives of the misfit function. Their drawbacks are mostly getting stuck in local minima, and medium computational cost.
Grid searches: a method that searches the entire parameter space systematically. Although they can map all the features of the inverse problem (by design), they are often much too computationally expensive.

What might be even more important, is that we have seen how to reduce the amount of possible solutions from infinitely many to at least a tractable amount using regularisation. There is only one fundamental piece still missing… Stay tuned for the last blog post in this series for the reveal of this mysterious missing ingredient!

Interested in playing around with inversion yourself? You can find a toy code about baking bread here.


Inversion 101 to 201 – Part 1: The forward problem

Inversion 101 to 201 – Part 1: The forward problem

The Geodynamics 101 series serves to showcase the diversity of research topics and methods in the geodynamics community in an understandable manner. We welcome all researchers – PhD students to professors – to introduce their area of expertise in a lighthearted, entertaining manner and touch upon some of the outstanding questions and problems related to their fields. This time, Lars Gebraad, PhD student in seismology at ETH Zürich, Switzerland, talks about inversion and how it can be used in geodynamics! Due to the ambitious scope of this topic, we have a 3-part series on this with the next part scheduled for publication next week! This week, however, we start with the basics that many geodynamicists will be familiar with: the forward problem.

Inversion methods are at the core of many physical sciences, especially geosciences. The term itself is used mostly by geophysicists, but the techniques employed can be found throughout investigative sciences. Inversion allows us to infer our surroundings, the physical world, from observations. This is especially helpful in geosciences, where one often relies on observations of physical phenomena to infer properties of the (deep) otherwise inaccessible subsurface. Luckily, we have many, many algorithms at our disposal! (…cue computer scientists rejoicing in the distance.)

Inverse theory in a nutshell
Credit: xkcd.

The process of moving between data and physical parameters goes by different names in different fields. In geophysics we mostly use inversion or inference. Inversion is intricately related to optimization in control theory and some branches of machine learning. In this small series of blog posts, I will try to shed a light on how these methods work, and how they can be used in geophysics and in particular geodynamics! In the end, I will focus a little on Bayesian inference, as it is an increasingly popular method in geophysics with growing potential. Be warned though:

Math ahead!

(although as little as possible)

First, for people who are serious about learning more about inverse theory, I strongly suggest Albert Tarantola’s “Inverse Problem Theory and Model Parameter Estimation”. You can buy it or download it (for free!) somewhere around here: http://www.ipgp.fr/~tarantola/.

I will discuss the following topics in this blog series (roughly in order of appearance):

1. Direct inversion
2. Gradient descent methods
3. Grid searches
4. Regularization
5. Bayesian inference
6. Markov chain Monte Carlo methods

The forward ‘problem’

Rather backwards, I will start the first blog post with a component of inversion which is needed about halfway through an inversion algorithm. However, this piece determines the physics we are trying to investigate. Choosing the forward ‘problem’ or model is posing the question of how a physical system will evolve, typically in a deterministic setting. This can be any relationship. The forward problem is roughly the same as asking:

I start with these materials and those initial conditions, what will happen in this system?

In this post I will consider two examples of forward problems (one of which is chosen because I have applied it everyday recently. Hint: it is not the second example…)

1. The baking of a bread with known quantities of water and flour, but the mass of the end product is unknown. A recipe is available which predicts how much bread will be produced (let’s say under all ratios of flour to water).

2. Heat conduction through a volcanic zone where geological and thermal properties and base heat influx are known, but surface temperatures are unknown. A partial differential equation is available which predicts the surface temperature (heat diffusion, Laplace’s equation).

Both these systems are forward problems, as we have all the necessary ingredients to simulate the real world. Making a prediction from these ‘input’ variables is a forward model run, a forward simulation or simply a forward solution. These simulations are almost exclusively deterministic in nature, meaning that no randomness is present. Simulations that start exactly the same, will always yield the same, predictable result.

Modelling our bread

For the first system, the recipe might indicate that 500 grams of flour produce 800 grams of bread (leaving aside the water for a moment, making case 1a), and that the relationship is linear in flour. The forward model becomes:

flour * 800 / 500 = bread

This is a plot of that relationship:

Simple flour to bread relationship: a 1D forward model.
I think I am getting hungry now. Does anyone else have the same problem?

Of course, if we add way too much water to our nice bread recipe, we mess it up! So, the correct amount of water is also important for making our bread. Hence, our bread is influenced by both quantities, and the forward model is a function of 2 variables:

bread = g(flour, water)

Let’s call this case 1b. This relationship is depicted here:

Simple flour and water to bread relationship: a 2D forward model.
Yes, I am definitely getting hungry.

It is important to note that the distinction between these two cases creates two different physical systems, in terms of invertibility.

Modelling our volcanic zone (and other PDE’s)

Modelling a specific physical system is a more complicated task. Typically, physical phenomena can be described by set of (partial) differential equations. Other components that are required to predict the end result of a physical deterministic system are

1. The domain extent, suitable discretisation and parameter distribution in the area of investigation, and the relevant material properties (sometimes parameters) at every location in the medium
2. The boundary conditions, and when the system is time varying (transient), additional initial conditions
3. A simulation method

As an example, the relevant property that determines heat conduction is the heat conductivity at every location in the medium. The parameters control much of how the physical system behaves. The boundary and initial conditions, together with the simulation methods are very important parts of our forward model, but are not the focus of this blog post. Typically, simulation methods are numerical methods (e.g., finite differences, finite elements, etc.). The performance and accuracy of the numerical approximations are very important to inversion, as we’ll see later.
The result of the simulation of the volcanic heat diffusion is the temperature field throughout the medium. We might however be only interested in the temperature at real-world measurement points. We consider the following conceptual model:

Conceptual model of heat diffusion in a heterogeneous thermal conductivity model. Domain of interest is the black box. Under homogeneous heat influx, the deterministic simulation will result in a decreasing temperature in the three measurement points from left to right.

The resulting temperatures at the measurement locations will be different for each point, because they are influenced by the heterogeneous subsurface. Real-world observations would for example look something like this

‘Observations’ (made by me of course) from the volcanic model.

Now that you have a good grasp of the essence of the forward problem, I will discuss the inverse problem in the next blog post. Of course, the idea of inversion is to literally invert the forward problem, but if it were as easy as that, I wouldn’t be spending 3 blog posts on it, now would I?

Oceans on Mars: the geodynamic record

Oceans on Mars: the geodynamic record

Apart from our own planet Earth, there are a lot of Peculiar Planets out there! In this series we take a look at a planetary body or system worthy of our geodynamic attention, and this week we are back to our own solar system, more precisely to our neighbour Mars. In this post, Robert Citron, PhD student at the University of California, Berkeley, writes about the links between oceans, shorelines, and volcanism on Mars. A David Bowie approved blog post!

Robert Citron

Whether Mars was warm enough in its early history to support large oceans remains controversial. Although today Mars is extremely cold and dry, several lines of geological evidence suggest early Mars was periodically warm and wet. Evidence for ancient liquid water includes river channels, deltas and alluvial fans, lakes, and even shorelines of an extensive ocean.

Such features are carved into Mars’ ancient crust, which contains a remarkable geologic record spanning from over 4 billion years ago to present. How and where fluvial erosion takes place is highly dependent on topography. However, Mars is a dynamic planet and the topography observed today does not necessarily represent the planetary surface billions of years ago. Geological markers that seem misaligned today, such as river flow directions and sea levels, may be more consistent with ancient topography. Geodynamic models of how the planet has changed shape over time can therefore be used to test and constrain evidence of water on early Mars.

Shorelines are one example where geodynamic models have helped interpret the geological record. Perhaps the most compelling evidence for ancient Martian oceans are the hypothetical palaeo-shorelines that border Mars’ northern lowland basin. However, the shorelines fail to follow an equipotential surface, or contours of constant elevation, which would be expected if they formed via a standing body of liquid water. One explanation is that the shorelines used to follow an equipotential surface, but subsequent changes to the planet’s shape warped them to their present-day elevations, which vary by up to several kilometers.

Two geodynamic processes that likely changed the global shape of Mars are surface loading and true polar wander. True polar wander occurs when a planet reorients relative to its spin axis, which reshapes the planet because it changes the location of the equatorial bulge produced by the planet’s rotation. Large scale true polar wander on Mars was examined by Perron et al. (2007), which found that it could have warped past shoreline profiles to their present-day topographic profiles.

Another possibility is that Martian shoreline markers were deflected by flexure associated with surface loads. The emplacement or removal of material on a planet’s surface can cause flexure and displacement of nearby crust. This is observed on Earth, where melting of glaciers has unburdened underlying crust, allowing for rebound and displacement of past shoreline markers. Similar processes could be at work on Mars, but on a global scale.

Mars topography: The massive Tharsis volcanic province (red) is situated on the boundary between the southern highlands (orange) and northern lowlands (blue). The lowlands may have been covered by one or more ancient oceans, and are bordered by palaeo-shorelines. Two of the most prominent shorelines are the older Arabia shoreline (dashed line) and younger Deuteronilus shoreline (solid line). Image constructed by R. Citron using MOLA data. Shoreline data from Ivanov et al. (2017) and Perron et al. (2007).

The largest load on Mars is the Tharsis rise, a volcanic province that dominates the topography and gravity of the planet. Tharsis was built by volcanic activity over hundreds of millions to billions of years. Its emplacement changed Mars’ shape on a global scale; in addition to the Tharsis rise, there is a circum-Tharsis depression and an antipodal bulge.

In recent work (Citron et al. 2018), we found that the present-day variations in shoreline elevations can be explained by flexure from Tharsis and its associated loading. Of the two most prominent Mars shorelines, the older (~ 4 billion years old) Arabia shoreline corresponds to pre-Tharsis topography, deformed by almost all of the flexure associated with Tharsis. The younger (~ 3.6 billion years old) Deuteronilus shoreline corresponds to late-Tharsis topography, requiring only ~17% of Tharsis loading to explain its variations in elevation. This suggests that the Arabia shoreline formed before or during the early stages of Tharsis, and the Deuteronilus shoreline formed during the late stages of Tharsis growth. The match between the present-day shoreline markers and ancient equipotential surfaces supports the hypothesis that the markers do indicate shorelines formed by an ancient ocean.

The timing of ancient Martian oceans is consistent with recent work by Bouley et al. (2016), which found that the Mars valley networks (ancient river channels) also better fit Mars’ pre-Tharsis topography. In the topography of Mars prior to Tharsis, the flow direction of the channels are more consistent with the topographic gradient, and the channels occur at latitudes and elevations where climate models predict water ice (resulting in ice melt) to form.

The timing of the shorelines and valley networks relative to Tharsis volcanism suggests a close link between the stability of water on Mars and volcanic activity. Atmospheric models predict a cold and icy early Mars, however it is possible that oceans may be more sustainable during periods of heightened volcanism. Tharsis activity has also been associated with outflow channels indicative of catastrophic flooding that may have inundated the northern plains with water. Further examination of the link between Tharsis volcanism and oceans could increase our understanding of early Mars habitability.

Further reading:

Bouley, S., Baratoux, D., Matsuyama, I., Forget, F., Séjourné, A., Turbet, M., & Costard, F. (2016). Late Tharsis formation and implications for early Mars. Nature531(7594), 344.

Citron, R. I., Manga, M., & Hemingway, D. J. (2018). Timing of oceans on Mars from shoreline deformation. Nature555(7698), 643.

Ivanov, M. A., Erkeling, G., Hiesinger, H., Bernhardt, H., & Reiss, D. (2017). Topography of the Deuteronilus contact on Mars: Evidence for an ancient water/mud ocean and long-wavelength topographic readjustments. Planetary and Space Science144, 49-70.

Matsuyama, I., & Manga, M. (2010). Mars without the equilibrium rotational figure, Tharsis, and the remnant rotational figure. Journal of Geophysical Research: Planets115(E12).

Perron, J. T., Mitrovica, J. X., Manga, M., Matsuyama, I., & Richards, M. A. (2007). Evidence for an ancient martian ocean in the topography of deformed shorelines. Nature447(7146), 840.

Ramirez, R. M., & Craddock, R. A. (2018). The geological and climatological case for a warmer and wetter early Mars. Nature Geoscience11(4), 230.

It’s just coding … – Scientific software development in geodynamics

The Spaghetti code challenge. Source: Wikimedia Commons, Plamen petkov 92, CC-BY-SA 4.0

As big software packages become a commonplace in geodynamics, which skills should a geodynamicist aim at having in software development? Which techniques should be considered a minimum standard for our software? This week Rene Gassmöller, project scientist at UC Davis, Computational Infrastructure for Geodynamics, shares his insights on the best practices to make scientific software better, and how we can work to translate these into our field. Enjoy the read!

Rene  Gassmöller

Nowadays we often equate geodynamics with computational geodynamics. While there are still interesting analytical studies to be made, and important data to be gathered, it is increasingly common that PhD students in geodynamics are expected to work exclusively on data interpretation, computational models, and in particular the accompanying development of geodynamic software packages. But as it turns out, letting an unprepared PhD student (or unprepared postdoc or faculty member for that matter) work on a big software package is a near guarantee for the project to develop into a sizeable bowl of spaghetti code (see figure above for a representative illustration).

Note, that I intentionally write about ‘software packages’ instead of ‘code’, as many of these packages — think of Gplates (Müller et al, 2018), ObsPy (Krischer et al, 2015), FeniCS (Alneas et al, 2015) , or the project I am working on, ASPECT (Heister et al, 2017) — have necessarily left the stage of a quickly written ‘code’ for a single purpose, and developed into multi-purpose tools with a complex internal structure. With this growing complexity, the activity of scientific ‘coding’ evolved into ‘developing software’. However, when students enter the field of geophysics, they are rarely prepared for this challenge. Hannay et al. (2009) report that while researchers typically spend 30% or more of their time developing software, 90% of them are primarily self-taught, and only few of them received formal training for writing software, including tests and documentation. Nobody told them: Programming and engineering software are two very different things. Many undergraduate and graduate geoscience curricula today include classes about the basics of programming (e.g. in Python, R, or Matlab), and also discuss numerical and computational methods. While these concepts are crucial for solving scientific problems, they are not sufficient for managing the complexity of growing scientific software. Writing a 50-line script is a very different task from contributing to an inherited and poorly documented PhD project of 1,000 lines, which again is very different from managing a multi-developer project of 100,000 lines of source code. A recurring theme is that these differences are only discovered when damage has already been done. Hannay et al. (2009) note:

Codes often start out small and only grow large with time as the software proves its usefulness in scientific investigations. The demand for proper software engineering is therefore seldom visible until it is “too late”.

But what are these ‘proper software engineering techniques’?

Best practices vs. Best techniques in practice

In a previous blog post, Krister Karlsen already discussed the value of version control systems for reproducibility of computational research. It is needless to say that these systems (originally also termed source code control systems, e.g. Rochkind, 1975) are just as valuable for scientific software development as they are for reproducibility of results. However, they are not sufficient for developing reliable scientific software. Wilson et al. (2014) summarize a list of 8 best practices that make scientific software better:

  1. Write programs for people, not computers.
    • A program should not require its readers to hold more than a handful of facts in memory at once.
    • Make names consistent, distinctive, and meaningful.
    • Make code style and formatting consistent.
  2. Let the computer do the work.
    • Make the computer repeat tasks.
    • Save recent commands in a file for re-use.
    • Use a build tool to automate workflows.
  3. Make incremental changes.
    • Work in small steps with frequent feedback and course correction.
    • Use a version control system.
    • Put everything that has been created manually in version control.
  4. Don’t repeat yourself (or others).
    • Every piece of data must have a single authoritative representation in the system.
    • Modularize code rather than copying and pasting.
    • Re-use code instead of rewriting it.
  5. Plan for mistakes.
    • Add assertions to programs to check their operation.
    • Use an off-the-shelf unit testing library.
    • Turn bugs into test cases.
    • Use a symbolic debugger.
  6. Optimize software only after it works correctly.
    • Use a profiler to identify bottlenecks.
    • Write code in the highest-level language possible.
  7. Document design and purpose, not mechanics.
    • Document interfaces and reasons, not implementations.
    • Refactor code in preference to explaining how it works.
    • Embed the documentation for a piece of software in that software.
  8. Collaborate.
    • Use pre-merge code reviews.
    • Use pair programming when bringing someone new up to speed and when tackling particularly tricky problems.
    • Use an issue tracking tool.

There is a lot to be said about each of these techniques, but that would be beyond the scope of this blog post (please see Wilson et al.’s excellent and concise paper if you are interested). What I would like to emphasize here is that these techniques are often requested, but rarely taught. What are peer code reviews? How do I gradually introduce tests and refactor a legacy code? Who knows if it is better to use unit testing, integration testing, regression testing, or benchmarking for a given change of the code? And do I really need to know the difference? After all, a common argument against using software development techniques in applied computational science disciplines boils down to:

  • We can not expect these software development techniques from geodynamicists.
  • We should not employ the same best practices as Google, Amazon, Apple, because they do not apply to us.
  • There is no time to learn/apply these techniques, because we have to conduct our research, write our publications, secure our funding.

While from a philosophical standpoint it is easy to dismiss these statements as not adhering to best practices, and possibly impacting the reliability of the created software, it is harder to tackle them from a practical perspective. Of course it is true that implementing a sophisticated testing infrastructure for a one-line shell command is neither useful nor necessary. Maybe the same is true for a 20 line script that is written to specifically convert one dataset into another, but in this case putting it under version control would already be useful in order to record your process and apply it to other datasets. And from my own experience it is extraordinarily easy to miss the threshold at 40-100 lines at which writing documentation and implementing first testing procedures become crucial to avoid cursing yourself in the future for not explaining what you did and why you did it. So why are there detailed instructions for lab notes and experimental procedures, but not for geodynamic software design and reliability of scientific software? Geoscience, chemistry, and physics have established multi-semester lab and field exercises, to drill students towards careful scientific analysis. Should we develop comparable exercises for scientific software development (beyond numerical methods and basic programming)? How would an equivalent of these classes look like for computational methods? And is there a point where the skills of software development and geodynamics research grow so far apart we have to consider them separately and establish a unique career track, such as the Research Software Engineer that is becoming more popular in the UK?

In my personal opinion we have made great progress over the last years in defining best practices for scientific software (see e.g. https://software.ac.uk/resources/online-sustainability-evaluation, or https://geodynamics.org/cig/dev/best-practices/). However, it is still considered a personal task to acquire the necessary skills and to find the correct balance between careful engineering and overdesigning software. Establishing courses and resources that discuss these questions could greatly benefit our community, and allow for a more reliable scientific progress in geodynamics.

Collaborative software development – The overlooked social challenge

The contributor funnel. The atmosphere and usability of a project influence how many users will join a project, how long they stick around, and if they will take responsibility for the project by contributing to it or eventually become maintainers. Credit: https://opensource.guide/

Now that we covered every topic a scientist can learn about scientific software development in a single blog post, what can go wrong when you put several of them together to work on a software package? Needless to say, a lot. No matter if your software project is a closed-source, intra-workgroup project, or an open-source project with users and developers spread over different continents, things are going to get exponentially more complicated the more people work on your software. Not only does discussion and interaction take more time, there will also be conflicting ideas about computational methods, software design, or implementation. Using state-of-the-art tools like collaborative development platforms (Github, Gitlab, Bitbucket, pick your favourite) and modern discussion channels like chats (Slack, Gitter), forums (Discourse), or video conferences (Skype, Hangouts, Zoom) can alleviate a part of the communication barriers. But ultimately, the social challenges remain. How does a project decide between competing goals of flexibility and performance? Who is going to enforce a code of conduct in a project to keep the development environment open and friendly? Does a project create a welcoming atmosphere that invites new contributions, or does it repel newcomers by unrealistic standards and inappropriate behavior? How should maintainers of scientific software deal with unrealistic feature requests by users? How to encourage new users to become contributors and take responsibility for the software they benefit from? How to compromise or combine providing improvements to the upstream project versus publishing them as scientific papers? How to provide credit to contributors?

In my opinion it is unfortunate that these questions about scientific software projects are even less discussed than the (now increasing) awareness of reproducibility. On the bright side, there is already a trove of experiences in the open-source community. The same questions about attribution and credit, collaboration and community-management, and correctness and security have been discussed over the past decades in open-source projects all over the world, and nowadays a good number of resources provide guidance, such as https://opensource.guide/, or the excellent book  ‘How to Run a Successful Free Software Project’ (Fogel, 2017). Not all of it can be transferred to science, but we would waste time and energy to dismiss these experiences and instead repeat their mistakes.

Let us talk about engineering scientific software

I realize that in this blog post I opened more questions than I answered. Maybe that is because I am not aware of the answers that are already out there. But maybe it is also caused by a lack of attention that these questions receive. I feel that there are no established guidelines for which software development skills a geodynamicist should have, and what techniques should be considered a minimum standard for our software. If that is the case, I would invite you to have a discussion about it. Maybe we can agree on a set of guidelines and improve the state of software in geodynamics. But at the very least I hope I inspired some thought about the topic, and provided some resources to learn more about a discussion that will likely grow more important over the coming years.

References:

M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes and G. N. Wells. The FEniCS Project Version 1.5. Archive of Numerical Software, vol. 3, 2015, http://dx.doi.org/10.11588/ans.2015.100.20553.

Fogel, K. (2017). Producing Open Source Software: How to Run a Successful Free Software Project. O'Reilly Media, 2nd edition.

Hannay, J. E., MacLeod, C., Singer, J., Langtangen, H. P., Pfahl, D., & Wilson, G. (2009). How do scientists develop and use scientific software?. In Proceedings of the 2009 ICSE workshop on Software Engineering for Computational Science and Engineering (pp. 1-8). IEEE Computer Society.

Heister, T., Dannberg, J., Gassmöller, R., & Bangerth, W. (2017). High accuracy mantle convection simulation through modern numerical methods–II: realistic models and problems. Geophysical Journal International, 210(2), 833-851.

Krischer, L., Megies, T., Barsch, R., Beyreuther, M., Lecocq, T., Caudron, C., & Wassermann, J. (2015). ObsPy: A bridge for seismology into the scientific Python ecosystem. Computational Science & Discovery, 8(1), 014003.

Müller, R.D., Cannon, J., Qin, X., Watson, R.J., Gurnis, M., Williams, S., Pfaffelmoser, T., Seton, M., Russell, S.H. & Zahirovic, S. (2018). GPlates–Building a Virtual Earth Through Deep Time. Geochemistry, Geophysics, Geosystems.

Open Source Guides. https://opensource.guide/. Oct, 2018.

Rochkind, M. J. (1975). The source code control system. IEEE transactions on Software Engineering, (4), 364-370.

Wilson, G., Aruliah, D.A., Brown, C.T., Hong, N.P.C., Davis, M., Guy, R.T., Haddock, S.H., Huff, K.D., Mitchell, I.M., Plumbley, M.D. and Waugh, B. (2014). Best practices for scientific computing. PLoS biology, 12(1), e1001745.