The CLI for busy Data Scientists and Engineers

I've been asked to give a talk on the command line interface for a mixed audience of Data Scientists and Engineers. Since the world is becoming ever more container based, I'm focussing on Linux.

Containers

Sometimes you need to diagnose things within the container. Jump into the container with:

docker exec -it CONTAINER_ID /bin/bash

To get the basic diagnostic tools mentioned below, you'll generally need to execute:

apt-get update               # you need to run this first
apt-get install net-tools    # gives you netstat
apt-get install iputils-ping # gives you ping
apt-get install procps       # gives you ps
apt-get install lsof        

Note these installations will all be gone next time you fire up the image as the underlying image does not change.

You can find out which module your favourite command belongs to by running something like this:

$ dpkg -S /bin/netstat
net-tools: /bin/netstat

Formatting

You can use regex expressions in grep with the -P switch. For example, let's search for lines that are strictly composed of two 5-letter words seperated by a space:

$ echo hello world  | grep -P "\w{5}\s\w{5}$"
hello world
$ echo hello wordle | grep -P "\w{5}\s\w{5}$"
$

You can extract elements from a string with an arbitrary delimiter with awk. For example, this takes the first and sixth elements from a line of CSV:

$ echo this,is,a,line,of,csv | awk -F',' '{print $1 " " $6}'
this csv
$

To prettify output, use column like this:

$ echo "hello, world. I love you!
goodbye cruelest world! So sad" | column -t
hello,   world.    I       love  you!
goodbye  cruelest  world!  So    sad
$

To print to standard out as well as to a file, use tee. For example:

$ echo And now for something completely different | tee /tmp/monty_python.txt

And now for something completely different

$ cat /tmp/monty_python.txt 

And now for something completely different

To capture everything typed and output to your terminal (very useful in training), use script:

$ script -f /tmp/my_keystrokes.log

Script started, file is /tmp/my_keystrokes.log

$ echo hello world

hello world

$ cat /tmp/my_keystrokes.log 

Script started on 2022-05-13 16:10:37+0100

$ echo hello world

hello world

$ cat /tmp/my_keystrokes.log 

$

Beware its recursive nature! Anyway, you stop it with an exit.

You can poll an output with watch. For example, this will keep an eye on the Netty threads in a Java application (IntelliJ as it happens):

$ watch "jstack `jps | grep Main | awk '{print \$1}'` | grep -A10 ^\\\"Netty\ Builtin\ Server"

Note that the $1 has been escaped and the quote mark within the quote has been triple escapted. The switch -A10 is just to show the 10 lines After what we pattern matched. Backticks execute a command within a command. Of course, we can avoid this escaping with:

$ watch "jstack $(jps | grep Main | awk '{print $1}') | grep -A10 ^\\\"Netty\ Builtin\ Server"

Note that $(...).

Resources

The command htop gives lots of information on a running system. Pressing P or M orders the resources by processor or memory usage respectively. VIRT and RES are your virtual memory (how much your application has asked for) and resident memory (how much it's actually using) the latter is normally the most important. The load average tells you how much work is backing up. Anything over the number of processors you have is suboptimal. How many processors do you have?

$ grep -c ^processor /proc/cpuinfo 

16

$

The top command also lists zombie tasks. I'm told that these are threads that are irretrievable stuck, probably due to some hardware driver issue.

File handles can be seen using lsof. This can be useful to see, for example, where something is logging. For instance, guessing that IntelliJ logs to a file that has log in its name, we can run:

$ lsof -p 12610 2>/dev/null | grep log$

java    12610 henryp   12w      REG              259,3    7393039 41035613 /home/henryp/.cache/JetBrains/IntelliJIdea2021.2/log/idea.log

$

The 2>/dev/null pipes errors (the 2) to a dark pit that is ignore.

To see what your filewall is dropping (useful when you've misconfigured a node), run:

$ sudo iptables -L -n -v -x

To see current network connections, run:

$ netstat -nap 

You might want to pipe that to grep LISTEN to see what processes are listening and on which port. Very useful if something already has control of port 8080.

For threads, you can see what's going on by accessing the /proc directory. While threads are easy to see in Java (jstack), Python is a little more opaque, not least because the Global Interpretter Lock (GIL) only really allows one physical thread of execution (even if Python can allow logical threads). To utilise more processors, you must start a heavyweight thread (see "multiprocessing" here). Anyway, find the process ID you're interested in and run something like:

$ sudo cat /proc/YOUR_PROCESS_ID/stack
[<0>] do_wait+0x1cb/0x230
[<0>] kernel_wait4+0x89/0x130
[<0>] __do_sys_wait4+0x95/0xa0
[<0>] __x64_sys_wait4+0x1e/0x20
[<0>] do_syscall_64+0x57/0x190
[<0>] entry_SYSCALL_64_after_hwframe+0x44/0xa9

as everything is a file in Unix, right? This happens to be the stack in Python code that is time.sleeping. You'll see a similar looking stack for a Java thread that happens to be waiting.

If you want to pin work to certain cores, use something like taskset. For example, if I wanted to run COMMAND on all but one of my 16 cores, I run:

taskset 0xFFFE COMMAND

This is very useful if some data munging is so intense it's bringing my system down. Using this, at least one thread is left for the OS.

Finally, vmstat gives you lots of information about the health of the box such as blocks being read/written from/to disk (bo/bi), the number of processes runnable (not necessarily running) and the number blcoked (r/b) and the number of context switches per second (cs)

Tips for effective MLOps

Some miscellaneous tips I've discovered after over a year of being hands-on with a clinical ETL pipeline.

Technology

Set up a local dev environment (Git, Pip, IDE of choice, Python environment etc). Being able to test locally cannot be more important. For instance, I've been using Palantir's Foundry but since it's just a wrapper around Spark, you can have a full, locally run, test suite using PySpark.

If you can't get the data on your laptop (GDPR etc), use simulated data. If your model cannot find a signal in the noise you control, it likely won't when it faces real data. Make these tests automated so it can be part of your CI/CD pipeline.

Keep everything as deterministic as you can. Using commands like SQL's first means you have no guarantees what will come out of the sausage machine. How can your QAs test that?

Prefer to write intermediary data sets to persistent storage and load them again in the next stage over processing everything in memory with a single output. It might make the pipeline take longer to run but it gives you something like an "audit log" of the steps in your model.

Have diagnositic transforms in your pipeline. Typically, they make assertions to ensure self-consistency (eg, all patients' hospital length-of-stay must be non-negative; the aggregated total of overnight patients cannot exceed the total number of beds). These transforms output any violations into data sets that will never be seen by the end user, just the MLOps team. Diagnostic transforms have the advantage over the local test in that they can use real data as they are embedded in the production pipeline.

Process

Testing data is harder than testing code so spare a thought for your QAs. Talk to the them long before you cut a line of code and establish whether they can test what you're building. If not, find out what you need to give them.

Talk to data owners, tell them what you're trying to do and ask them how they'd do it. They often tell you things you'd never know a priori - like which data is, ahem, "less reliable". 

Spend some time (a week?) just playing with the data (EDA) so things become more apparent. It's important to know, for instance, if a column called "region" in A&E data set is the region of the patient or the care provider. Since a patient in A&E may be treated the other side of the country (perhaps they were taken ill while on holiday), the two regions might have nothing to do with one another.

Make sure everybody is singing from the same hymn sheet with a ubiquitous language. Calculating the number of weeks a patient waits couldn't be easier, right? Take the number of days between the start and end date and divide by 7, right? But do we round up or down? And MS SQLServer has its own weird way of doing this calculation.

Be iterative. Don't clean data, build features and then try models in a waterfall manner. 

Data

Don't forget the effect of time on your data (data drift). For instance, care provider networks often get restructured so does hospital trust X in 2020 go by a different name today? Patients move around the country so the same person may be associated with more than one locale. Did more data sources come online over time?

How often is the data updated and how stale is it? These are not the same thing. For instance, the data may be updated daily but its content lags behind by several months.

Were dodgy assumptions made on upstream data sets that you need to be aware of? For instance, removing dead people from the the patients data set might have made sense upstream when dealing with vaccination data but not for historical accident and emergency data.

When unifying data sets, watch out for semantically equivalent values that may be lexically different (eg, if one data set uses an ethnicity of "British white" and another that says "White British", then a union of the two will look like two different ethnicities). This problem becomes harder when one data set uses, say, different age bands. How would you reconcile them?

More MLOps: bug hunting

We have several ML models for which we calculate the odds (the coefficients in a linear regression model) and the risks (the difference in probabilities between the factual and the counterfactual). 

Previously, we saw that for one of our six models, the sign of the odds and risks differed when they should always be the same.

Fortunately, we caught this by putting a transform into our pipeline that takes the output of both the linear regression model and factual/counterfactual probabilities and simply checks them.

Part of our ML pipeline. Some transforms create production data, others diagnostic.

In the graph above, anything with _regression_tests appended is a diagnostic transform while everthing else is production data that will ultimately be seen by the customer.

Recently, a diagnostic test was reporting:

difference in risk (0.14285714285714285) does not agree with the sign of the coefficient (-4182.031458784066) for STP QR1 when p-value = 0.988606800326618

Transform containing pipeline diagnostic data.

Aside: there are other violations in this file but they're not unexpected. For instance, there are 42 health regions in England but the data has some random rubbish from the upstream data capture systems. The last row is due to us expecting at least 10 statistically signifcant outputs from our model when we only found 9. That's fine, we'll turn up the squelch on that one. A red flag would be if the number of insights dropped to zero. This happened when we needed the outcome to be binary but accidentally applied a function of type ℕ => {0,1} twice. Oops.

In a previous post, we found the problem was the "risks" and the "odds" transforms, although using the same input file, were consuming that file at significantly different times. The problem manifested itself when most regions for that one model had different signs. But this single diagnostic we see in the violations above is different. Only one region (QR1) differs. All the others in the data set were fine. 

Looking at the offending region, the first things to note was that we had very little data for it. In fact, of the 42 health regions, it had the least by a large margin.

Number of rows of data per health region for a particular model

So, the next step was to build a Python test harness and using these 127 rows and run the data against production code locally. To run this on my laptop is so much easier as I can, say, create a loop and run the model 20 times to see if there was some non-determinsm happening. Apart from tiny changes in numbers due to numerical instability, I could not recreate the bug.

So, like in the previous post, could my assumption that the data is the same for both risk and odds calculations be flawed? The pipeline was saying that the two descendant data sets were built seconds after the parent.

Both risk and odds data descend from the same parent

Could there be some bizarre AWS eventual consistency issue? This seems not the case as S3 read/write consistency now has strong consistency (see this AWS post from December 2020).

Conclusion

Why just one region in one model has a different sign is still a mystery. Multiple attempts to see the same bug have failed. However, since there is so little data for this region, the p-values are very high (c. 0.99) we disregard any insights anyway. We'll just have to keep an eye out for it happening again.

The Scala 3 Type System - part 2

This is the successor to my first post about Scala 3's great new type system.

Peano Axioms (again)

For fun, I'm modelling the natural numbers thus:

   sealed trait Nat
   case object Zero extends Nat
   case class Succ[N <: Nat](n: N) extends Nat 

   transparent inline def toNat(inline x: Int): Nat =
     inline x match
       case 0 => Zero
       case 1 => Succ(Zero)
       case _ => Succ(toNat(x - 1))

But despite me using transparent, the function toNat just returns the super-type toNat not Succ[Succ[Succ[....

scala> toNat(3)
val res0: Succ[Nat] = Succ(Succ(Succ(Zero)))

So, with a bit of help from the guys on the Scala metaprogramming Discord channel, I came up with this:

  transparent inline def toNat[I <: Int]: ToNat[I] =
    inline erasedValue[I] match
      case z: 0    => Zero
      case s: S[t] => Succ(toNat[t])

which worked:

scala> toNat[5]

val res0: Succ[Succ[Succ[Succ[Succ[Zero.type]]]]] = Succ(Succ(Succ(Succ(Succ(Zero)))))

You'll notice two major differences. First, toNat no longer takes any arguments. We use the parameterized type in lieu of a parameter

The second is this bizarre S[_] (which is imported from scala.compiletime.ops.int). It gives us the notion of a successor. Let's look more closely at it.

S Express(ion)

From the Scala source code:

  /** Successor of a natural number where zero is the type 0 and successors are reduced as if the definition was:

   *
   *  ```scala
   *  type S[N <: Int] <: Int = N match {
   *    case 0 => 1
   *    case 1 => 2
   *    case 2 => 3
   *    // ...
   *    case 2147483646 => 2147483647
   *  }
   *  ```
   *  @syntax markdown
   */
  type S[N <: Int] <: Int

What sorcery is this? I see no implementation of that huge match statement!

PhillHenry
The bit I didn't understand was S[t] => ... I assume the compiler enforces that t must be preceeding natural number of S. Is this compiler magic? Or can anybody enforce that constraint in plain Scala?

ragnar
Its compiler magic in the same way that the literal 3 is compiler magic. Basically val x: S[3] = 4 and val y: S[x.type] = 5 work as you would expect, the same holds true in reverse for pattern matching.

(the Scala/metaprogramming Discord server, 3 April 2022).

Compiler Options

By default, the compiler will limit the depth of recursion to 32 levels (after all, you need to tell the compiler when things have become pathological). This is far too small for us and 1024 sounds a bit more sensible. However, I use SBT and wasted an hour telling SBT exactly how to use this value. You need to add:

     scalacOptions ++= Seq("-Xmax-inlines",  "1024")

Note -X is broken into 2 elements of the Seq. If you foolishly try to add it as a single element as it would appear in the CLI, the compiler/SBT will ignore it.

The takeaway point

A good approach to metaprogramming is to imagine that the functions are being called but the arguments correspond to the parameterized types, not the function's parameters.  

Python/Pandas Tips

A few miscellaneous tips that I keep referring to so I thought I'd make a post.

Flatmap in Python

Given a list of lists t,

flat_list = [item for sublist in t for item in sublist]

This is a more Pythonic solution [StackOverflow] than 

from itertools import chain

where you might do something like this:

spark_map = F.create_map([F.lit(x) for x in chain(*python_map.items())])

Here, the chain is basically flattening a list (tuples) in a list (.items).

Pandas

Even if your main code does not use Pandas, it can be very useful for writing assertions in your tests.

Say, you have data in the form of lists of lists in variable rows. You can organise this into a Pandas data frame with something like [SO]:

df = pd.DataFrame(rows, range(len(rows)), ["COLUMN_1", "COLUMN_2"])

Now, if I want to get the COLUMN_1 value for the highest COLUMN_2, I do:

df.iloc[df['PREDICTION_COLUMN'].idxmax()]["COHORT"]

One can sort Pandas data frames with something like this [StackOverflow]:

df = df.reindex(df.coefficients.abs().sort_values().index)

where the absolute value of the column coefficients is what I want to sort on.

One can filter with something like:

cleaned = df[(df["p_values"] < 0.05) & ((df["coefficients"] > 0.1) | (df["coefficients"] < -0.1))]

Or, if you want to filter something that does not contain a string, say:

fun[~fun["feature"].str.contains("intercept")].head(50)

Pandas doesn't try to guess the types in the file so if you have a date for instance, you need to do something like:

df['time'] = pd.to_datetime(df['time'])

to convert a string to a timestamp type.

To create a new column from old columns, the syntax is simply:

df["z_score"] = df["coefficients"] / df["standard_error"]

You can concat dataframes with something like 

pd.concat([df, df2], axis=1)

where the axis is that of the concatenation - 0 for rows (vertically), 1 for columns (horizontally) etc.

But to join, you need something like:

pd.merge(df1, df2, on=cols, how="inner")

where the parameters are obvious.

Tests

Useful tip when running tests: if you want only a subset to run, you can use a glob ignore. For example:

python -m pytest integration_tests/ --ignore-glob='integration_tests/test_actionable*.py'

will ignore tests that start with test_actionable.

Riddles of Sphinx

Sphinx is a tool to autogenerate documentation. It's generic but is often used in the Python ecosystem. 

Once configured, you run it via an old school Makefile. It will then do to your Python code was Javadoc does to Java. It doesn't appear to have the ability to mix code and documentation like ScalaTest.

How to use it

I found this StackOverflow answer a lot clearer than the official documentation. The upshot of it is:

1. Create a docs folder. In this, you can run sphinx-quickstart if you like (separating source from build - see SOURCEDIR and BUILDDIR in the Makefile) but one way or another, your conf.py file will need to know where to find the code. If you've followed the instruction in the SO answer, it's two levels up so conf.py needs to have:
   
   import os
   import sys
   sys.path.insert(0, os.path.abspath(os.path.join('..', '..', 'src')))
   
   (note those two '..') in it. Also, note that this is the top-level directory of your code. If there are relative imports in your code, these need to be relative to this directory.
2. You'll also need the extensions to parse Python.
3. The index.rst file needs to be told what to pull in. You do this by adding a reference under toctree (let's arbitrarily call it module) that corresponds to another file (module.rst) in the same directory. 
4. In turn, our arbitrarily names module.rst needs in its toctree a reference to our final sub-module. We'll call it my_docs and it will necessitate a file called my_docs.rst in the same directory.
5. Here, in the automodule block, we'll tell Sphinx which file want documented. Note that you don't need the .py extension.
6. Note that the .rst files do not need to conform to your directory structure.

PDFs

Sphinx natively creates Latex files but to create PDF, I needed to run:

sudo apt-get install latexmk

which allowed me to run:

make latexpdf

Dependencies

I had a problem with third-party libraries on which my code depended. Having found my code, Sphinx was telling me no module named pyspark when it hit the import statement in my code. This SO answer helped out but it left me with a nasty hard-coding in my conf.py. This may be more down to the way Python manages its dependencies (which is horrible compared to Java).

Writing effective Sphinx docstring

You can reference constants in your docstring with something like :py:const:`FQN.CONSTANT_NAME` (where FQN is the fully qualified path to your Python file). If you also put something like this in your .rst file:

.. autoattribute:: FQN.CONSTANT_NAME

then Sphinx will not only provide a hyperlink to the constant but also print out the value of that constant.

Executing the code in comments

For this, try doctest. It not only executes code in the comments but checks the results.

MLOps debugging: an example

In our ML pipeline, we use generalized linear models to calculate the odds of certain clinical outcomes. We showed this to the client but odds were hard for them to understand. "Can we have probabilities instead?" they asked. 

So, having trained the GLMs, we fit the same data and calculate the average probabilities for each cohort. We then bastardise the data to create the counterfactuals For example, socioeconomic status is of interest so let's make everybody the same level (counterfactual) then make our predictions again and compare with the factual results.

Now, one would imagine that if the corresponding GLM coefficient were positive, then the average probability from the counterfactuals would be less than that of the factual data. This is simple maths, right? You dot product your feature vector with the coefficients and plug it into the sigmoid function. A positive coefficient would lead to a bigger product and so a smaller denominator.

We were mostly seeing this but about 10% of our predicitions violated this rule.

Investigation

A transform compares the two approaches (risks and odds) and with a bit of Pandas:

> import pandas as pd

> df = pd.read_csv("/home/henryp/Documents/CandF/fitted_vs_coeffs.txt", delimiter='\t')

> sub = df[["coefficients", "cohort_risk", "counterfactual_cohort_risk", "p_values"]]

> sub["diff"] = sub["cohort_risk"] - sub["counterfactual_cohort_risk"]

> sub = sub.reindex(sub["diff"].abs().sort_values().index)

> sub[-5:]

    coefficients  cohort_risk  counterfactual_cohort_risk  p_values  feature                       model_id               diff

25     -0.520122     0.239911                    0.184961  0.121023  ethnic_category_one_hot_Other a_and_e_cancer         0.054950

23     -0.561247     0.277648                    0.219477  0.006748  ethnic_category_one_hot_Asian a_and_e_cancer         0.058171

44     -0.292537     0.545455                    0.480899  0.821720  ethnic_category_one_hot_Asian wait_18_ophthalmology  0.064556

50     -0.224494     0.307358                    0.229469  0.480723  ethnic_category_one_hot_Black a_and_e_cancer         0.077889

8      -5.340222     0.636364                    0.551042  0.957624  ethnic_category_one_hot_Asian wait_18_ophthalmology  0.085322

we note that the p-values indicate a lack of statistical significance in the coefficient even if the difference between factual and counterfactual probabilities can be pretty large.

But given the same data should produce the same p-values and coefficients each time (modulo floating point issues), how could this happen?

In Generalized Logistic Regression models, the coefficients don't change much on each run, even when shuffling the data. This is true irrespective of p-values. Let's demonstrate.

Shuffling

The Scala code taken from the official Spark docs for a GLR (here) gives this:

glr.fit(dataset.orderBy(rand())).coefficients.toArray.zip(glr.fit(dataset.orderBy(rand())).coefficients.toArray).map{ case (x, y) => x - y }

res24: Array[Double] = Array(1.231653667943533E-16, -2.220446049250313E-16, 2.220446049250313E-16, 0.0, 1.1102230246251565E-16, -2.220446049250313E-16, -1.6653345369377348E-16, 2.220446049250313E-16, 2.220446049250313E-16, -2.220446049250313E-16)

Running this a few times shows the differences are miniscule even though the p-values are pretty large:

glr.fit(dataset.orderBy(rand())).summary.pValues

res25: Array[Double] = Array(0.989426199114056, 0.32060241580811044, 0.3257943227369877, 0.004575078538306521, 0.5286281628105467, 0.16752945248679119, 0.7118614002322872, 0.5322327097421431, 0.467486325282384, 0.3872259825794293, 0.753249430501097)

Any differences in the coefficients between runs are almost certainly due to floating point arithmetic.

Modification

However, things are very different when we modify the data; even a small perturbation can radically change coefficients with high p-values.

Note that there are only 500 rows in this test data set. But let's fit two models where on average we drop one row at random. Then we see things like this:

def compare2(): Unit = {

  val tol         = 0.002 

  val fitted      = glr.fit(dataset.where(rand() > tol))

  val pValues     = fitted.summary.pValues 

  val stdErrors   = fitted.summary.coefficientStandardErrors

  val coeffs_rand = fitted.coefficients.toArray  

  val other       = glr.fit(dataset.where(rand() > tol))

  val diffs       = other.coefficients.toArray.zip(coeffs_rand).map{ case (x, y) => x - y }

  val diffs_pc    = diffs.zip(coeffs_rand).map{ case (x, y) => (x / y) * 100 } ;  

  val meta        = pValues.zip(diffs).zip(diffs_pc).zip(stdErrors).map{ case(((p, d), pc), s) => (p, d, pc, s) }

  println("%-15s%-15s%-15s%-15s".format("p-value", "difference", "% difference", "std error"))

  meta.sortBy(x => -math.abs(x._1)).foreach { case (p, d, pc, s) => println(s"%-15.3f%-15.3f%-15.3f%-15.3f".format(p, d, pc, s)) }

}

scala> compare2()

p-value        difference     % difference   std error      

0.977          0.051          -224.562       0.793          

0.602          0.101          -24.475        0.792          

0.574          0.043          9.546          0.793          

0.551          -0.035         7.432          0.795          

0.459          0.044          -7.145         0.828          

0.456          0.085          14.590         0.776          

0.283          -0.068         -7.897         0.803          

0.275          0.133          -15.230        0.796          

0.134          -0.067         -5.484         0.811          

0.006          0.119          5.188          0.830        

We need run it only a dozen or so times to see that the coefficients with the largest p-value tend to have the largest percentage discrepancies between the two fits and that the sign can change.

The Culprit

These CLI antics lead me to question the "given the same data" assumption. A closer look at the timestamps of the input data sets revealed they were built a week apart (Palantir's Foundry will by default silently fall back to data sets on a different branch if they have not been built on the current branch. This data may be stale). Well, our data slowly grows over time. So, could these small differences be responsible for big swings in the GLM coefficients? I created in the pipeline a single snapshot of the data from which both the risk and odds values derived and there were zero discrepancies. Sweet.