This manuscript describes an approach to computer model emulation and sensitivity analysis using parallel and distributed computing on high performance computing (HPC). Specifically, they fit scalar Gaussian process (GP) emulators to grid-cell level Earth system model output in an embarassingly-parallel fashion, and use the fitted emulators to conduct main effects-style sensitivity analysis using generalized additive models (GAMs). Each model fit occurs on a core independently of the other fits and speedup is realized by running many such fits simultaneously on many cores, with no between-node communication required.

The resulting parallel speedups seem reasonable and are significant in magnitude.  It might be easier to see the parallel efficiency if plotted differently, e.g. on a log-log plot, or with a theoretical maximum speedup curve (equal to number tasks) superimposed, or speedup as a fraction of theoretical.  I am not sure I would say that optimized walltime is nearly constant up to three orders of magnitude in task count; we'd expect it to be constant when the tasks are less than the concurrency cap of 200 cores, but it should (and does) increase after that.  The speedup graphs are not as informative as they could be, for that reason, because the baseline was never run (and the speedup calculated) for task sizes exceeding the number of cores, where different speedup behavior is expected (of course this is because of the increasing expense of running the baseline) - perhaps at least one long baseline analysis at (say) 400 cores could be run?

Other than that, a few other points to consider:



Could the GAM have been fit directly to the data instead of to GP emulator samples? The GAM essentially is a simple emulator (neglecting parameter interactions).

It may be worth pointing out that the GAM sensitivity analysis is like a (unnormalized) first-order Sobol' sensitivity analysis with the assumption that there are no between-parameter interactions.

Could this software be extended to analyze time series at a grid cell? This would require a different form of emulator and different distribution of data onto tasks.  (Obviously a spatial analysis would be harder since communication across boundaries would become inevitable, although spatially-partitioned Gaussian processes can be applied.)

This paper discusses how the speed of training of Gaussian Process emulators and associated fitting of Generalised Additive Models can be be improved by an optimised workflow, use of in-memory I/O, and (I think) an alternative algorithm (as opposed to an alternative implementation) of the necessary prediction matrix.

The overall outcome is a significant improvement in performance, although it it is not clear to me whether the scale up goes beyond 200 cores. Is Figure 5, showing the results of running 10K tasks through 200 cores? In which case, how does the speedup go beyond what can be achieved in 200 cores alone? There is no point showing more than 200 tasks if the rest is linear?

This little example shows the problem I have with the paper. I think the work done is very interesting, and worthy of publication, but the write-up doesn't quite do it justification. The authors assume a lot of the readers and the material is not always well ordered.

I would like to see this appear as it is a good piece of work, and a nice example of optimisation a workflow, but would recommend a substantial revision of the material in terms of order and clarity (I say substantial in that i think it would be of substantial benefit, but probably wouldn't take that long and I would expect the editor could decide on this without coming back to me). 

In particular, I'd like to see

1. Some clarity in section 3.1 as to whether figure 4 is necessary at all? It appears that the innovation here was to take out a loop over grid points and use slurm job arrays? Is that it? The comment about "spawning competing rscript processes" is presumably the key part of this? But if this is just making better use of the parallelism and isn't changing the algorithm, why do we need figure 4? And if it is changing the algorithm, how is it doing it, there are no details?

2. A better explanation in section 3.2 of of what is going on, that doesn't mix together algorithm and implementation. I found this really hard to understand. Pseudo function calls are mixed with algorithmic notation. (I said "I think" in my summary, because I have vaguely understood what was going on here, but I would like to actually understand it based on the text.)

3. The tables in B1-B4 need some thinking about. I understand that these are exemplars showing the difference between the two implementations at four points, but these sort of things are notorious. Why not show scatter plots (or whatever) of the difference across the domain (of however many points it was), for some of the parameters, not all, so we can have confidence that these are not the "best" results, and are typical?

I should say that Appendix A is excellent and very useful.

Minor comments:

P 2, line 28; the sentence beginning "However" implies that somehow the size of PPEs is linked the complexity of the model. Might be cleaner to say something like: "However, complex models are expensive, and costs are prohibitive beyond o(100) members"

(There is also an assumption in this sentence about how big PPEs need to be, which is surely problem related).

The paragraph then carries on and conflates PPES with the assumption that they are only useful if they have to be expanded in size ... this paragraph probably just needs reordering, since the justification for the expansion in size follows the assumption.

P2 ,line 39/40 bottlenecks are or bottleneck is

P 6 On page 5 we were told that each task had 16 GB of memory, but here we are told that peak memory use exceeded 80 GB. I've gotten confused as to what is going on. Is this an example of the "otherwise noted"?

P 7. It is not necessary to keep citing Lawrence et al after the first time JASMIN was mentioned.