Estimating NICAS parameters

This page is a practical guide to running a BUMP parameter estimation job: diagnosing correlation or localization length-scales from an ensemble, building a NICAS operator, and reusing the result. It is model-independent - the examples use a generic geometry, and the places where your model matters are called out explicitly.

For what the underlying drivers do, see General overview. For the SABER blocks that consume the output, see BUMP SABER blocks.

What estimation does

There are two separate pieces of machinery, and it is worth being clear that you can run either one on its own:

• HDIAG samples the ensemble, bins pairs of points by their separation distance, fits a correlation function to the resulting curve, and gives you length-scales. Depending on the drivers you enable, these are correlation length-scales (cor_rh, cor_rv) or localization length-scales (loc_rh, loc_rv).

• NICAS builds the actual convolution operator from a pair of length-scales - whether those came from HDIAG or you specified them yourself.

So there are three useful jobs: diagnostics only, NICAS from length-scales you already know, and the two chained together. All three are shown below.

Note

There is no hdiag SABER block and no separate estimation application. HDIAG is switched on implicitly, by the drivers you enable. Everything runs through the BUMP_NICAS block’s calibration section.

What a model must provide

BUMP takes its vertical coordinate from a geometry field. The field is named by the fields metadata key <variable>.vert_coord if present, and otherwise defaults to a geometry field called vert_coord.

The units of that field are the units of every vertical length-scale you configure or read back. Some model interfaces fill it with a physical depth or height in metres; others fill it with the level index, in which case a vertical length-scale of 2.0 means two model levels and has nothing to do with metres. Check what your model interface puts in vert_coord before interpreting cor_rv or setting a vertical length-scale by hand. Horizontal length-scales are always in metres.

Running the estimation

Parameter estimation is run by the ErrorCovarianceToolbox application. Each model builds its own executable:

mpiexec -n <ntasks> <model>_error_covariance_toolbox.x myconfig.yaml

quench_error_covariance_toolbox.x uses QUENCH simplified model, useful for trying a configuration out before pointing it at real data.

The top-level yaml has three sections:

geometry:
  <...>                            # model-specific geometry

background:
  <...>                            # model-specific state; supplies the valid time
                                   # and, for some models, masks and coordinates

background error:
  covariance model: SABER
  ensemble:                        # NOTE: here, not inside the saber block
    <...>
  saber central block:
    saber block name: BUMP_NICAS
    calibration:
      <...>                        # everything BUMP-specific goes in here

The ensemble lives at the background error level. Four forms are accepted:

Key	Meaning
ensemble	A set of states; perturbations are formed by removing the ensemble mean.
ensemble pert	Increments read directly from disk; used as-is.
ensemble base + ensemble pairs	Perturbations formed as the difference of two sets of states (useful, e.g. for NMC method).
ensemble pert on other geometry + ensemble geometry	Increments read on a different geometry and interpolated.

Add iterative ensemble loading: true at the same indentation level to read members one at a time instead of holding the whole ensemble in memory. On a large grid this is usually necessary.

Note

Within the BUMP-specific section, ensemble sizes: total ensemble size is filled in automatically from the number of members you listed. Set it by hand only when re-reading previously stored moments, where there is no member list to count.

Running HDIAG (diagnostics only)

The simplest useful job, and the one to start with. It answers “what are my ensemble’s length-scales?” without committing to an operator, and it is much cheaper than building NICAS.

saber central block:
  saber block name: BUMP_NICAS
  calibration:
    io:
      data directory: path_to_bump_directory
      files prefix: my_bump_files
    drivers:
      multivariate strategy: univariate  # one diagnostic per variable
      compute moments: true              # required for any HDIAG job
      compute covariance: true
      compute correlation: true
      write diagnostics: true            # -> my_bump_files_diag.nc
      write diagnostics in yaml: true    # -> my_bump_files_diag.yaml (easy to plot)
      compute nicas: false               # diagnostics only
    sampling:
      computation grid size: 1000        # Sc1: where correlations are measured
      diagnostic grid size: 200          # Sc2: where local diagnostics are stored
      distance classes: 15
      distance class width: 200.0e3      # in m; see "Choosing the sampling parameters"
      reduced levels: 11                 # moving window of levels, NOT the level count
      local diagnostic: true             # produce a map rather than one global profile
      averaging length-scale: 1000.0e3   # radius over which local statistics are pooled
      max number of draws: 500000        # the default of 10000 is too small; see pitfalls
    diagnostics:
      target ensemble size: 30           # REQUIRED, and must be > 3
    output model files:
    - parameter: cor_rh
      file:
        <...>                            # model-specific output file configuration
    - parameter: cor_rv
      file:
        <...>
    - parameter: stddev
      file:
        <...>

Add compute localization: true to also diagnose localization length-scales, and request loc_rh / loc_rv in output model files. Localization length-scales are corrected for the sampling noise of a finite ensemble, so they are broader than the raw correlation length-scales; diagnostics: target ensemble size is the ensemble size they are optimal for.

Constructing a NICAS operator from explicit length-scales

This option will build a NICAS operator from a user-provided set of explicit length-scales and does not require an ensemble. See BUMP SABER blocks for the length-scale conventions and the per-level profile: form.

saber central block:
  saber block name: BUMP_NICAS
  calibration:
    io:
      data directory: path_to_bump_directory
      files prefix: my_bump_files
    drivers:
      multivariate strategy: duplicated  # one shared operator; groups are named "common"
      compute nicas: true
      write local nicas: true
    nicas:
      resolution: 8.0
      explicit length-scales: true
      horizontal length-scale:
      - groups: [common]
        value: 1000.0e3                  # metres
      vertical length-scale:
      - groups: [common]
        value: 0.5                       # units of the model's vert_coord field

Running HDIAG and NICAS construction together

The full chain: diagnose the length-scales, then immediately build the operator from them. Take the example in Running HDIAG (diagnostics only) and add the NICAS drivers:

drivers:
  multivariate strategy: univariate
  compute moments: true
  compute covariance: true
  compute correlation: true
  compute localization: true         # for a localization operator
  compute nicas: true
  write local nicas: true
nicas:
  resolution: 8.0                    # no explicit length-scales -- HDIAG supplies them

Running in two steps instead is often better on a large grid: diagnostics are cheap and let you inspect the length-scales before paying for the operator, and the sampling and moments can be reused via overriding sampling file and overriding moments file.

BUMP YAML keys and vocabulary reference

This is a reference for BUMP-specific jargon relevant to estimation. The exhaustive set, with defaults, is in saber/src/saber/bump/BUMPParameters.h.

The sampling subsets

BUMP’s log messages refer to these constantly, so they are worth knowing:

Name	Set by	What it is
Sc0	The model geometry; the geometry section of the YAML	Every model grid point.
Sc1	computation grid size	The points where correlations are actually measured.
Sc2	diagnostic grid size	The points where a local diagnostic is stored - the knots of the length-scale map.
Sc3	distance classes, angular sectors	The partner points paired with each Sc1 point at each separation.
halo A / halo B	The computational grid partitioning	What an MPI task owns / what it owns plus what it borrows from neighbouring tasks.

drivers: keys

Key	Description
multivariate strategy	univariate, duplicated, duplicated and weighted or crossed. Also determines the group names - see the pitfalls below.
compute moments	Accumulate the ensemble moments. Required for any HDIAG job, unless read moments is used instead.
compute covariance, compute correlation	Diagnose covariance / correlation. Either one enables HDIAG.
compute localization	Diagnose localization length-scales. Requires the two above.
compute nicas	Build the NICAS operator.
write local nicas, write global nicas	Write the operator as one file per MPI task, or a single file.
read local nicas, read global nicas	Read a previously written operator instead of building one.
write diagnostics	Write <files prefix>_diag.nc.
write diagnostics in yaml	Also write <files prefix>_diag.yaml - the fitted profiles in plain text.
write diagnostics detail	Per-component curves. Requires fit: number of components greater than 1.
write local sampling, write moments	Store the sampling / moments so a later run can reuse them.

sampling: keys

Key	Description
computation grid size	Number of Sc1 points (\(n_{c1}\)). The main driver of memory and moment-computation cost.
diagnostic grid size	Number of Sc2 points (\(n_{c2}\)), i.e. the resolution of the length-scale map.
distance classes	Number of separation bins.
distance class width	Width of each bin, in metres. Bins run from 0 to (distance classes - 1/2) × distance class width.
angular sectors	1 for isotropic. Greater than 1 (and even) diagnoses anisotropy - see below.
reduced levels	A moving window of levels centred on each level on which vertical diagnostic are performed (not the number of levels, to save CPU and memory).
local diagnostic	true gives a spatially varying map (but increases the cost); false gives one global profile.
averaging length-scale	Radius, in metres, over which Sc1 statistics are pooled into each Sc2 diagnostic. Mutually exclusive with averaging latitude width; exactly one must be set.
max number of draws	Cap on the random draws used to fill the Sc3 pairs. The default of 10000 is far too small on a real grid - see the pitfalls.

diagnostics: and fit: keys

Key	Description
diagnostics: target ensemble size	The ensemble size the localization diagnostics are computed for. By default, it is equal to the input ensemble size.
diagnostics: gaussian approximation	Use the Gaussian approximation for the asymptotic quantities. Worth trying if the localization diagnostics look unstable.
diagnostics: localization option	Select the localization calculation method among several options: default: optimal localization estimation described in equation (165) of this note, might yield unstable results depending on the ensemble distribution. from_squared_correlation: squared correlation-based estimation described in equation (172) of this note, significantly more stable. nice_with_table and nice_without_table based on the NICE method of Vishny et al. (2024), with or without pre-computed lookup tables.
fit: number of components	Number of Gaspari-Cohn functions summed to form the fitted function. Default 1.
fit: horizontal filtering length-scale	Smooths the fitted length-scale field spatially, in metres. Applied after the fit, so it compounds with averaging length-scale, which pools before it.

nicas: keys

Key	Description
resolution	Unitless. The number of points used to discretize the Gaspari-Cohn function from its origin to the limit of its support. Use at least 3; 4-8 if affordable.
max horizontal grid size	Upper bound on the NICAS subgrid size. Cost grows quadratically with resolution and inversely with the length-scale, so if the theoretical grid size exceeds this bound a lower effective resolution is computed and printed in the log. If the effective resolution falls below 3 the run aborts, and you must lower resolution or raise this bound.
explicit length-scales	true to supply length-scales yourself rather than taking them from HDIAG.
horizontal length-scale, vertical length-scale	Per-group value: or per-level profile:.
filter mode, filter resolution	Build the operator as a low-pass filter rather than a covariance square root. In filter mode resolution is ignored in favour of filter resolution. Filter mode can be used e.g. for ensemble scale separation.
same horizontal convolution	Apply one horizontal operator level by level, with no vertical convolution. Much cheaper on a many-level grid, and the right choice for a purely horizontal filter.

Choosing the sampling parameters

computation grid size, diagnostic grid size and averaging length-scale are coupled, and guessing at them independently is the most common way to get a poor length-scale map. The number of Sc1 points pooled into each Sc2 diagnostic is approximately

\[n_{\mathrm{pooled}} \approx n_{c1} \, \frac{\pi r^2}{A}\]

where \(r\) is averaging length-scale and \(A\) is the area covered by the valid points of your domain.

• Aim for \(n_{\mathrm{pooled}}\) between 50 and 100. Below roughly 30 the per-Sc2 fits are noisy and the map looks speckled; above roughly 200 the extra samples no longer improve the fit. Choose computation grid size to land in that range.

• Then choose diagnostic grid size from the map resolution you want. The length-scale field has knots only at Sc2 points and is interpolated between them, and the pooling has already smoothed at scale \(r\), so Sc2 spacing much finer than \(r/2\) buys nothing: \(n_{c2} \approx A / (r/2)^2\).

• distance class width should be roughly twice the mean grid spacing. Bins narrower than the mesh can only draw partners from high latitudes, where converging meridians make the zonal spacing smaller - the result is a latitude-biased sample rather than an empty one, so nothing warns you.

Cost scales as \(n_{c1}\) for memory, \(n_{c1} n_{c2}\) for the averaging, and \(n_{c2}\) for the fitting, which performs a nonlinear minimisation per Sc2 point, per level, per component.

Neither computation grid size nor diagnostic grid size can exceed the number of candidate points, which BUMP reports early in the log:

Decimate full grid, at least    10000 points required,   91600 valid points found

Reusing the output

Switch calibration: to read: and replace the compute drivers with read ones:

saber central block:
  saber block name: BUMP_NICAS
  read:
    io:
      data directory: path_to_bump_directory
      files prefix: my_bump_files
    drivers:
      multivariate strategy: duplicated  # must match the run that wrote the files
      read local nicas: true

Warning

Anything that shaped the stored operator must be restated identically when reading it - multivariate strategy, filter mode and same horizontal convolution in particular. Only the construction-side keys drop away.

Common pitfalls

max number of draws is too small by default, and fails silently

BUMP fills the Sc3 pair slots by drawing uniformly random points over the domain and binning them by distance. The probability that a random point lands in the ring of radius \(r\) and width \(dc\) is \(r \, dc / 2R^2\) - of order \(10^{-5}\) for a first class of a few tens of kilometres. With the default cap of 10000 draws the loop exits with most of the short-distance slots unfilled.

Unfilled slots are quietly masked out. There is no warning, and the progress bar reaches 100% regardless of whether the slots were filled, so the run looks healthy. Because the starved classes are exactly the ones that constrain the fit, the length-scale comes out fitted on the tail of the curve: biased long and noisy.

The best way to confirm the diagnosis of this issue is to check the ..._diag.nc file. If the raw correlation/localization curves are missing for many bins at the beginning of the curve, then max number of draws is indeed too small.

Set the value of max number of draws to:

\[f (1/p) \ln(n_{c1} n_{c3} n_{c4})\]

with \(p\) the probability above, and \(f\) is a fudge factor of 2-5. Do not simply set it to a very large value: the loop performs one MPI allgather per draw, so an over-large value becomes the run time.

Note

With angular sectors greater than 1 the loop can never reach its completion condition and always consumes the full max number of draws. Choose the value deliberately in that case.

reduced levels is a window, not a count

It is a moving window of levels centred on each level. Setting it equal to the number of model levels asks BUMP to correlate the top of the model with the bottom, which is rarely meaningful and is expensive: memory scales as reduced levels × number of levels, and with iterative ensemble loading: true four large moment arrays are allocated rather than two.

On a 75-level grid, setting reduced levels to 75 rather than to 11 can cost several GB per MPI task. The symptom is a job that appears to hang - it is thrashing while it allocates and fills the moment arrays. Set it to cover the vertical structure you expect, plus a margin.

Group names might depend on multivariate strategy

Group names can be specified in the model section, for instance like:

saber central block:
  saber block name: BUMP_NICAS
  read:
    model:
      groups:
      - group name: my_first_group
        variables: [var1, var2]
      - group name: my_second_group
        variables: [var3]

However by default, the group names depend on the multivariate strategy:

• For univariate and crossed: one group per variable, the group name being the variable name.

• For duplicated and duplicated and weighted: a single group for all variables, named “common”.

This is important to specify group-specific parameters like explicit length-scales. A mismatch aborts with:

rh is missing for <variable name>

The same naming is used for the netCDF groups inside <files prefix>_nicas, so building with one strategy and reading with another might give NetCDF: No group found. A simple solution is to use the alias key in the io section. For instance if the NICAS file was created for the variable air_horizontal_streamfunction with the univariate strategy, and is read with the duplicated strategy:

saber central block:
  saber block name: BUMP_NICAS
  read:
    io:
      alias:
      - in code: common
        in file: air_horizontal_streamfunction

Vertical length-scales may be auto-filled

With explicit length-scales: true, a group that has a horizontal length-scale but no vertical one gets a vertical length-scale of zero, with a warning in the log. That is usually what you want for a purely horizontal operator. Note the auto-fill only applies once the horizontal length-scale has matched a real group, so a group-name error surfaces as a complaint about rh rather than rv.

Every output model files entry needs its own file: block

BUMP passes the file: sub-configuration straight to the model’s writer. A missing file: key yields an empty configuration rather than an error, so the failure surfaces later and further away, as the model writer complaining about whichever key it happens to read first.

In particular, do not use a YAML merge key to share a common block - <<: *anchor splices the keys in beside parameter:, leaving file: absent. Write each entry out in full, and give each parameter a distinct output file, or they will overwrite each other.

Translating abort messages

BUMP’s namelist checks report the internal Fortran variable rather than the yaml key:

Abort message	What to set
new_mom or load_mom required for new_hdiag	drivers: compute moments (or read moments)
ne should be larger than 3	diagnostics: target ensemble size
fit_ncmp should be larger than 1 for write_hdiag_detail	fit: number of components, or drop drivers: write diagnostics detail
local_rad or local_dlat should be positive	sampling: averaging length-scale (or averaging latitude width)
local_rad and local_dlat cannot be positive at the same time	set exactly one of the two
rh is missing for <name>	the groups: names under nicas: horizontal length-scale do not match the strategy
nc1 should be larger than 2	sampling: computation grid size
nc3 should be positive	sampling: distance classes
dc should be positive	sampling: distance class width
nl0r should be positive	sampling: reduced levels
ens_nsub should be a divider of ens_ne	ensemble sizes: sub-ensembles must divide the ensemble size

Going further

Anisotropy. Set sampling: angular sectors to an even number greater than 1. Pairs are then binned by bearing as well as distance, and the fit produces a horizontal tensor instead of a single radius, available as the output parameters cor_rh1, cor_rh2 and cor_rhc. An equivalent isotropic cor_rh is still written, so results stay comparable with an isotropic run. Note that this multiplies the number of pair slots to fill.

Multi-component fits. Set fit: number of components to 2 or 3 to fit a sum of Gaspari-Cohn functions rather than a single one, which captures a sharp near-origin peak that one function cannot. The components are fitted greedily: the first is fitted to the raw curve, its contribution is subtracted, and the next is fitted to the residual. Outputs become per-component, so each output model files entry needs a component: index, and cor_a gives the amplitudes. NICAS itself consumes only one component, selected by nicas: overriding component in file, so a multi-component fit is primarily a diagnostic tool.