API

ClimaCalibrate

Calibrate a forward model against observations using EnsembleKalmanProcesses.jl.

Implement forward_model and observation_map for a subtype of AbstractModelInterface, then pass a backend, an EnsembleKalmanProcess, and that interface to calibrate. The backend decides where the ensemble runs: JuliaBackend in the current process, WorkerBackend across Distributed.jl workers, or an HPCBackend as one scheduler job per ensemble member.

See the documentation at https://CliMA.github.io/ClimaCalibrate.jl/stable/.

project_dir()

Return the directory of the currently active Julia project.

This is the default experiment_dir, i.e. what an HPCBackend job script is given as --project unless the model interface overrides it.

AbstractModelInterface

Abstract supertype for user-defined calibration experiments.

Users subtype this to define their experiment-specific configuration and dispatch the calibration interface functions.

Required interface

Subtypes must implement:

• forward_model(interface, iteration, member) which runs the forward model for a single ensemble member.
• observation_map(interface, iteration) which processes model output and returns a G_ensemble matrix.

To use the HPCBackend, the subtypes must also implement:

• model_interface_filepath(interface) which returns the path to the file that defines the model interface. The HPCBackend job script includes this file so that all interface functions defined on the subtype are available in the worker process.

Optional interface

• analyze_iteration(interface, ekp, g_ensemble, prior, output_dir, iteration) which inspects results after each ensemble update. The default implementation logs the mean constrained parameters and covariance-weighted error.
• postprocess_g_ensemble(interface, ekp, g_ensemble, prior, output_dir, iteration) which transforms g_ensemble before the ensemble update. The default implementation returns g_ensemble.

For HPCBackend, the subtypes can also implement:

• experiment_dir(interface) which returns the Julia project directory passed as --project to the job script. The default implementation is to return project_dir().
• exeflags(interface) which returns additional flags (e.g. --threads 4) passed to the Julia executable in the job script. The default implementation is to return the empty string.

Examples

struct MyModelInterface <: ClimaCalibrate.AbstractModelInterface
    config::String
end

function ClimaCalibrate.forward_model(
    interface::MyModelInterface,
    iteration,
    member,
)
    # Run the model using interface.config
end

function ClimaCalibrate.observation_map(interface::MyModelInterface, iteration)
    # Read model outputs and return G_ensemble matrix
end

forward_model(interface::AbstractModelInterface, iteration, member)

Execute the forward model simulation with the given configuration.

This function must be overridden by the user's model interface, dispatching on their subtype of AbstractModelInterface.

The parameters EKP drew for this member are on disk, at parameter_path; write the model's output somewhere observation_map can find it, conventionally under path_to_ensemble_member.

Examples

function ClimaCalibrate.forward_model(
    interface::MyModelInterface,
    iteration,
    member,
)
    member_dir = ClimaCalibrate.path_to_ensemble_member(
        interface.output_dir,
        iteration,
        member,
    )
    parameters = ClimaCalibrate.parameter_path(
        interface.output_dir,
        iteration,
        member,
    )
    run_my_model(; parameters, output_dir = member_dir)
end

See also observation_map.

observation_map(interface::AbstractModelInterface, iteration)

Run the observation map for the specified iteration.

This function must be implemented for each calibration experiment, dispatching on the user's subtype of AbstractModelInterface.

Returns

The G ensemble matrix: column m holds ensemble member m's output, ordered to match the observation entry for entry. A member whose forward model failed should leave a column of NaNs.

Examples

function ClimaCalibrate.observation_map(interface::MyModelInterface, iteration)
    (; output_dir, ensemble_size) = interface
    G_ensemble = fill(NaN, n_observation_entries, ensemble_size)
    for member in 1:ensemble_size
        member_dir = ClimaCalibrate.path_to_ensemble_member(
            output_dir,
            iteration,
            member,
        )
        G_ensemble[:, member] .= process_member_output(member_dir)
    end
    return G_ensemble
end

See also forward_model, postprocess_g_ensemble.

analyze_iteration(interface::AbstractModelInterface, ekp, g_ensemble, prior, output_dir, iteration)

Analyze results after updating the ensemble and before starting the next iteration.

This function is optional to implement.

For example, one may want to print information from the ekp object or plot g_ensemble.

postprocess_g_ensemble(
    interface::AbstractModelInterface,
    ekp,
    g_ensemble,
    prior,
    output_dir,
    iteration
)

Postprocess g_ensemble after evaluating the observation map and before updating the ensemble.

model_interface_filepath(interface::AbstractModelInterface)

Return the path to the file that defines the model interface.

The HPCBackend job script includes this file so that all interface functions defined on the AbstractModelInterface subtype (e.g. forward_model, observation_map, and any optional overrides) are available in the worker process, along with their required packages.

experiment_dir(interface::AbstractModelInterface)

Return the path to the experiment's Julia project directory.

The HPCBackend uses this to construct the job script command:

julia --project=$experiment_dir $exeflags -e '...'

so that each ensemble member's forward model job runs with the correct project environment. By default, returns project_dir() (the currently active project).

You should override this in your AbstractModelInterface subtype if your experiment lives in a separate project directory.

exeflags(::AbstractModelInterface)

Return additional flags passed to the Julia executable in the HPCBackend job script.

The HPCBackend constructs each ensemble member's job command as:

julia --project=$experiment_dir $exeflags -e '...'

Override this in your AbstractModelInterface subtype to pass extra flags such as --threads or -O0. By default, returns "" (no extra flags).

ClimaCalibrate.Calibration

The calibration loop itself.

calibrate writes the ensemble members' parameters, runs the forward model for each member on the chosen backend, evaluates the observation map, and updates the ensemble, repeating until it runs out of iterations or EKP terminates. What it writes goes under a single output directory, which is also what it reads to resume an interrupted run.

calibrate(
    backend,
    ekp::EKP.EnsembleKalmanProcess,
    interface::AbstractModelInterface,
    n_iterations,
    prior,
    output_dir,
)

Run a full calibration with ekp and prior for n_iterations on the given backend, storing the results of the calibration in output_dir.

The loop is the same for all backends: write the ensemble members' parameters, run the forward model for each member, evaluate the observation map, and update the ensemble. Only the middle step differs, and that is what the backend selects. See JuliaBackend, WorkerBackend, and HPCBackend.

If output_dir already contains a calibration, it is resumed: completed iterations are skipped, as are ensemble members that recorded a completed forward model. Pass a fresh output_dir to start over.

Returns

The EnsembleKalmanProcess at the end of the run.

Examples

ekp = ClimaCalibrate.calibrate(
    ClimaCalibrate.JuliaBackend(),
    ekp,
    MyModelInterface(output_dir, ensemble_size),
    10,
    prior,
    output_dir,
)

See also initialize, last_completed_iteration.

AbstractHPCConfig

Scheduler directives, modules, and environment variables for the jobs an HPCBackend submits.

Subtypes:

• SlurmConfig: configuration for the Slurm clusters.
• PBSConfig: configuration for the PBS clusters.

Interface

All subtypes of AbstractHPCConfig must have the following fields:

• directives::OrderedDict{Symbol, Any}: Scheduler directives (e.g., resource requests, time limits, etc.).
• modules::Vector{String}: List of modules to load in the job environment.
• env_vars::OrderedDict{String, Any}: Environment variables to set for the job environment.

Subtypes must also provide the methods:

• generate_directives(config): Return a string of scheduler directives for the job script.
• generate_modules(config): Return a string of module load commands for the job script.
• generate_env_vars(config): Return a string of environment variable export commands for the job script.

SlurmConfig <: AbstractHPCConfig

A configuration holding Slurm directives, modules, and environment variables that will be used when creating job scripts by the SlurmBackends.

SlurmConfig(;
    directives = Pair{Symbol, Any}[],
    modules = String[],
    env_vars = Pair{String, Any}[],
)

Create a SlurmConfig specifying the directives, modules, and env_vars for SlurmBackends.

Defaults

The default directive is

• :gpus_per_task: 0.

The default environment variables are

• CLIMACOMMS_DEVICE: "CPU" or "GPU" depending on the job directives,
• CLIMACOMMS_CONTEXT: "MPI".

Examples

This example creates a Slurm configuration for a job with a single task, using 12 CPUs and 1 GPU, and a runtime of 720 minutes. It loads the latest version of climacommon and explicitly sets environment variables for ClimaComms.

ClimaCalibrate.SlurmConfig(;
    directives = [
        :ntasks => 1,
        :gpus_per_task => 1,
        :cpus_per_task => 12,
        :time => 720,
    ],
    modules = ["climacommon"],
    env_vars = [
        "CLIMACOMMS_CONTEXT" => "SINGLETON",
        "CLIMACOMMS_DEVICE" => "CUDA",
    ],
)

PBSConfig <: AbstractHPCConfig

A configuration holding PBS directives, modules, and environment variables that will be used when creating job scripts by the DerechoBackend.

PBSConfig(;
    directives = Pair{Symbol, Any}[],
    modules = String[],
    env_vars = Pair{String, Any}[],
)

Create a PBSConfig specifying the directives, modules, and env_vars for the DerechoBackend.

The supported directives are: time, queue, ntasks, cpus_per_task, gpus_per_task, and job_priority. These directive names follow the Slurm naming convention (e.g., time instead of walltime). Any other directives provided will be ignored.

Defaults

The default directives are

• queue: "main@desched1",
• account: "UCIT0011",
• ntasks: 1,
• cpus_per_task: 1,
• gpus_per_task: 0,
• job_priority: "regular".

account is the project the job is charged to, the -A directive. The default is the CliMA allocation on Derecho, so set it to your own.

The default environment variables are

• CLIMACOMMS_DEVICE: "CPU" or "GPU" depending on the job directives,
• CLIMACOMMS_CONTEXT: "MPI".

Examples

This example creates a PBS configuration for a job with a single task, using 12 CPUs and 1 GPU, and a runtime of 720 minutes. It loads the latest version of climacommon and explicitly sets environment variables for ClimaComms.

ClimaCalibrate.PBSConfig(;
    directives = [
        :ntasks => 1,
        :gpus_per_task => 1,
        :cpus_per_task => 12,
        :time => 720,
    ],
    modules = ["climacommon"],
    env_vars = [
        "CLIMACOMMS_CONTEXT" => "SINGLETON",
        "CLIMACOMMS_DEVICE" => "CUDA",
    ],
)

ClimaCalibrate.Backend

Where an ensemble member's forward model runs, and how it is submitted.

A backend is the only thing that changes when a calibration moves from a laptop to a cluster: JuliaBackend runs members one at a time in the current process, WorkerBackend distributes them over Distributed.jl workers, and each HPCBackend submits one scheduler job per member.

This module also holds the scheduler plumbing the HPC backends need: the job configs (SlurmConfig, PBSConfig), job submission and cancellation (submit_job, cancel_job), and job status (JobStatus, job_status).

JuliaBackend(; failure_rate = 0.5)

Run the ensemble members one at a time in the current process.

Keyword Arguments

• failure_rate::Float64: The fraction of an iteration's ensemble members that may fail before the calibration is halted [-]. The default is 0.5.

Examples

backend = ClimaCalibrate.JuliaBackend()
tolerant = ClimaCalibrate.JuliaBackend(; failure_rate = 0.9)

WorkerBackend(; failure_rate, worker_pool, empty_pool_timeout)

Run each ensemble member's forward model on a Distributed.jl worker.

Members are handed to workers as they become free, so a calibration can start before all workers have connected. Add workers with add_workers on a cluster, or with Distributed.addprocs locally; see SlurmManager and PBSManager for the Slurm and PBS cluster managers.

Keyword Arguments

• failure_rate::Float64: The fraction of an iteration's ensemble members that may fail before the calibration is halted [-]. The default is 0.5.
• worker_pool: A worker pool created from the workers available.
• empty_pool_timeout::Int: How long (in seconds) an iteration will wait on an empty worker pool before erroring, so an asynchronous calibration cannot hang forever when no workers ever start. Defaults to 21600.

Examples

wait(ClimaCalibrate.add_workers(4; cluster = :local))
ClimaCalibrate.@worker_setup include("my_model.jl")
backend = ClimaCalibrate.WorkerBackend()

HPCBackend <: AbstractBackend

Backend that submits one scheduler job per ensemble member.

Each job starts a fresh Julia process, includes the file returned by ClimaCalibrate.model_interface_filepath, and runs one member's forward model. Prefer this over WorkerBackend when forward models are long-running, need internal parallelism, or will not all fit in one allocation.

Subtypes:

• SlurmBackend: the Slurm clusters.
• DerechoBackend: NSF NCAR Derecho, which uses PBS.

SlurmBackend <: HPCBackend

Abstract supertype for the clusters that use the Slurm scheduler: CaltechHPCBackend, ClimaGPUBackend, and GCPBackend.

Job submission, status queries, and cancellation are implemented once for this type; the concrete backends differ only in which modules they load and how they launch MPI.

CaltechHPCBackend(config::SlurmConfig)
CaltechHPCBackend(; directives, modules, env_vars, failure_rate, job_timeout)

Submit one scheduler job per ensemble member to Caltech's high-performance computing cluster.

The second form builds the SlurmConfig from its keyword arguments.

Examples

backend = ClimaCalibrate.CaltechHPCBackend(;
    directives = [:time => 60, :ntasks => 1, :cpus_per_task => 8],
    modules = ["climacommon"],
)

See also failure_rate.

ClimaGPUBackend(config::SlurmConfig)
ClimaGPUBackend(; directives, modules, env_vars, failure_rate, job_timeout)

Submit one scheduler job per ensemble member to CliMA's private GPU server.

The second form builds the SlurmConfig from its keyword arguments.

Examples

backend = ClimaCalibrate.ClimaGPUBackend(;
    directives = [:time => 60, :ntasks => 1, :cpus_per_task => 8],
    modules = ["climacommon"],
)

See also failure_rate.

GCPBackend(config::SlurmConfig)
GCPBackend(; directives, modules, env_vars, failure_rate, job_timeout)

Submit one scheduler job per ensemble member to CliMA's private GCP server.

The second form builds the SlurmConfig from its keyword arguments.

Examples

backend = ClimaCalibrate.GCPBackend(;
    directives = [:time => 60, :ntasks => 1, :cpus_per_task => 8],
    modules = ["climacommon"],
)

See also failure_rate.

DerechoBackend(config::PBSConfig)
DerechoBackend(; directives, modules, env_vars, failure_rate, job_timeout)

Submit one scheduler job per ensemble member to NSF NCAR's Derecho supercomputing system.

The second form builds the PBSConfig from its keyword arguments.

Examples

backend = ClimaCalibrate.DerechoBackend(;
    directives = [:time => 60, :ntasks => 1, :cpus_per_task => 8],
    modules = ["climacommon"],
)

See also failure_rate.

failure_rate(backend)

Return the fraction of an iteration's ensemble members that may fail before backend halts the calibration.

job_timeout(backend::HPCBackend)

Return the number of seconds backend waits for a running job before giving up.

The clock starts when a job leaves the queue, so time spent waiting for an allocation does not count against it.

backend_type()

Return the AbstractBackend type that suits the current machine, identified by gethostname(). Defaults to JuliaBackend when the host matches no known cluster.

This returns a type, not a backend that calibrate accepts. Construct a backend by calling the type. When it is an HPCBackend, pass a config with backend_type()(; directives, modules, env_vars); a JuliaBackend takes none of those and is built with JuliaBackend().

get_backend()

Deprecated alias for backend_type.

SlurmManager(ntasks=get(ENV, "SLURM_NTASKS", 1))

The ClusterManager for Slurm clusters, taking in the number of tasks to request with srun.

To execute the srun command, run addprocs(SlurmManager(ntasks)).

Keyword arguments can be passed to srun: addprocs(SlurmManager(ntasks), gpus_per_task=1).

By default the workers will inherit the running Julia environment.

To run a calibration, call calibrate(WorkerBackend(), ...).

To run functions on a worker, call remotecall(func, worker_id, args...).

PBSManager(ntasks)

The ClusterManager for PBS/Torque clusters, taking in the number of tasks to request with qsub.

To execute the qsub command, run addprocs(PBSManager(ntasks)). Unlike the SlurmManager, this will not nest scheduled jobs, but will acquire new resources.

Keyword arguments can be passed to qsub: addprocs(PBSManager(ntasks), nodes=2)

By default, the workers will inherit the running Julia environment.

To run a calibration, call calibrate(WorkerBackend(), ...)

To run functions on a worker, call remotecall(func, worker_id, args...)

get_manager(cluster = :auto, nworkers = 1)

Return the ClusterManager for cluster, which is one of :slurm, :pbs, or :auto to pick whichever scheduler's commands are on PATH.

:local workers do not need a manager, so add_workers handles that case before calling this.

add_workers(
    nworkers;
    device = :gpu,
    cluster = :auto,
    time = DEFAULT_WALLTIME,
    kwargs...
)

Add nworkers worker processes to the current Julia session, automatically detecting and configuring for the available computing environment.

This does not wait for the workers to connect. Each worker is submitted as an individual allocation and adds itself to GLOBAL_WORKER_POOL once it has started and loaded its code, so a calibration can begin with an empty pool and pick up workers as they join.

The returned Task runs the (blocking) submission; wait on it to block until all submissions have been processed. Submitted jobs are cancelled automatically when the process exits (via an atexit hook); call cancel_worker_jobs to tear them down earlier.

Use @worker_setup (instead of @everywhere) to load model code so that workers joining later get the same setup.

Arguments

• nworkers::Int: The number of worker processes to add.
• device::Symbol = :gpu: The target compute device type, either :gpu (1 GPU, 4 CPU cores) or :cpu (1 CPU core).
• cluster::Symbol = :auto: The cluster management system to use. Options:
  • :auto: Auto-detect available cluster environment (SLURM, PBS, or local)
  • :slurm: Force use of SLURM scheduler
  • :pbs: Force use of PBS scheduler
  • :local: Force use of local processing (standard addprocs)
• time::Int = DEFAULT_WALLTIME: Walltime in minutes, will be formatted appropriately for the cluster system
• workers_per_node::Int = 1: Number of workers to run per node.
• kwargs: Other kwargs can be passed directly through to addprocs.

Returns

A Task running the submission. wait on it to block until all workers have been submitted; the workers themselves join the pool as they connect.

Examples

# On a cluster: four GPU workers, each its own allocation
wait(ClimaCalibrate.add_workers(4; time = 120))

# Locally, for debugging
wait(ClimaCalibrate.add_workers(2; cluster = :local))

# On a cluster that charges for whole nodes, four workers per allocation
wait(ClimaCalibrate.add_workers(8; workers_per_node = 4))

See also @worker_setup, cancel_worker_jobs, calibration_worker_pool.

@worker_setup expr

Like Distributed.@everywhere, but the expression is also recorded and replayed on any worker that joins later.

using/import statements run on the main process first (to precompile once), and the current source path is propagated so relative include works on workers. As with @everywhere, local variables must be interpolated with $.

calibration_worker_pool()

Return the process-wide GLOBAL_WORKER_POOL, which is what a WorkerBackend draws ensemble members from.

Cluster workers add themselves via the :register hook. Workers added by other means (e.g. plain addprocs/LocalManager or pre-existing workers) are picked up here: each is claimed with _claim_worker and initialized in the background, so it joins the pool once it has the code to run a forward model.

A worker enters workers() when Distributed registers it, which is before the :register hook runs, so pooling ids straight from workers() can hand a member to a worker that has yet to load ClimaCalibrate. Claiming is what keeps the two paths from either racing or initializing the same worker twice.

The name is package-specific because Distributed exports a default_worker_pool of its own, which makes the unqualified name ambiguous under using Distributed, ClimaCalibrate.

cancel_worker_jobs(jobname = worker_jobname())

Cancel all scheduler jobs submitted for workers in this session with scancel (Slurm) or qdel (PBS). This tears down both connected workers (by cancelling their allocation) and any still-pending jobs.

Jobs submitted by add_workers share the job name worker_jobname, so they are cancelled together. Safe to call when no matching jobs exist.

Registered as an atexit hook whenever workers are launched onto a scheduler, so that jobs are not orphaned when the main process exits. It may also be called directly to tear down workers early.

set_worker_logger()

Set the worker's global logger to write to worker_$worker_id.log in its working directory.

Call this from the worker process. add_workers does so for each worker it starts.

Returns

The SimpleLogger that was installed.

set_worker_loggers(workers = workers())

Set the global logger to a simple file logger for the given workers.

map_remotecall_fetch(f::Function, args...; workers = workers())

Call function f from each worker and wait for the results to return.

foreach_remotecall_wait(f::Function, args...; workers = workers())

Call function f from each worker.

JobInfo

A submitted scheduler job: which backend submitted it, its scheduler ID, and the script it runs.

Returned by submit_job, and the argument to job_status, cancel_job, and requeue_job.

JobStatus

An enum representing the current status of a job.

Values

• PENDING: The job is queued and waiting to be scheduled.
• RUNNING: The job is currently executing.
• COMPLETED: The job finished running.
• FAILED: The job terminated with an error as reported by the scheduler.

Use ispending, isrunning, issuccess, isfailed, and iscompleted to query the status of a JobInfo. Each of those queries the scheduler, so ask once and test the result rather than calling several of them on the same job.

Examples

status = ClimaCalibrate.job_status(job)
ClimaCalibrate.iscompleted(status) && ClimaCalibrate.isfailed(status)

See also job_status.

job_status(job::JobInfo)

Return the current job status.

See JobStatus.

job_status(::SlurmBackend, job::JobInfo)

Return the status of job.

squeue only lists jobs that are still queued or running, so a job that has left the queue is looked up with sacct, which is the only way to distinguish a job that succeeded from one that failed, timed out, or was cancelled.

See JobStatus.

job_status(::DerechoBackend, job::JobInfo)

Return the status of job.

See JobStatus.

ispending(job::JobInfo)
ispending(status::JobStatus)

Return true if job is pending (i.e. waiting to be scheduled).

isrunning(job::JobInfo)
isrunning(status::JobStatus)

Return true if job is currently running.

issuccess(job::JobInfo)
issuccess(status::JobStatus)

Return true if job completed successfully.

isfailed(job::JobInfo)
isfailed(status::JobStatus)

Return true if job failed.

iscompleted(job::JobInfo)
iscompleted(status::JobStatus)

Return true if job has finished, either successfully or with a failure.

submit_job(backend::SlurmBackend, job_script::String)

Submit a job that runs job_script with backend.

The job_script should be generated with make_job_script.

submit_job(backend::DerechoBackend, job_script::String)

Submit a job that runs job_script with backend.

The job_script should be generated with make_job_script.

requeue_job(job::JobInfo)

Requeue job by cancelling the job and resubmitting it again.

This function will requeue the job even if the job is completed.

cancel_job(job::JobInfo)

Cancel the job.

cancel_job(::SlurmBackend, job::JobInfo)

Cancel job by running the command scancel.

cancel_job(::DerechoBackend, job::JobInfo)

Cancel job by running the command qdel.

cancel_jobs_at_exit(backend::HPCBackend)

Register an exit hook to cancel all jobs submitted by backend when the Julia process exits.

job_records(backend::HPCBackend)

Return a vector of JobInfos that were requested with backend.

write_job_script(filepath, job::JobInfo)

Write the job scheduler script for job to filepath.

This is useful for debugging the script that was submitted to the backend.

make_job_script(
    backend::SlurmBackend,
    job_body;
    job_name = "slurm_job",
    output = "output.txt",
)

Make a job script with job_body for the backend.

The job body must be a single Julia command.

make_job_script(
    backend::DerechoBackend,
    job_body;
    job_name = "pbs_job.txt",
    output = "output.txt",
)

Make a job script with job_body for the backend.

The job body must be a single Julia command.

initialize(ekp::EKP.EnsembleKalmanProcess, prior, output_dir)

Initialize a calibration, saving the initial parameter ensemble to a folder within output_dir.

If output_dir already holds a calibration, the stored first-iteration EnsembleKalmanProcess is returned instead and nothing is written. Overwriting it would replace the parameters that the completed forward models were run with: ekp is typically rebuilt on restart, and unless the caller seeded the RNG, its initial ensemble is a fresh random draw. The ensemble update would then pair G(u_old) from the checkpointed members with u_new.

last_completed_iteration(output_dir)

Return the last completed iteration of the calibration in output_dir, or 0 if none has completed.

An iteration counts as complete once its G_ensemble.jld2 exists and the next iteration's eki_file.jld2 has been written, i.e. once the ensemble update that consumed it finished. calibrate resumes from the iteration after this one.

Examples

ClimaCalibrate.last_completed_iteration(output_dir)

See also load_latest_ekp, model_completed.

terminated_iteration(output_dir)

Return the iteration at which the EnsembleKalmanProcess scheduler terminated the calibration in output_dir, or nothing if it ran to the end.

update_ensemble! saves the next iteration's eki_file.jld2 whether or not the scheduler terminated, and the ensemble it saves is the one that was already there. last_completed_iteration alone would send a restart off to run those same parameters again.

Examples

julia> ClimaCalibrate.terminated_iteration(output_dir)
4

save_G_ensemble(output_dir::AbstractString, iteration, G_ensemble)

Save the ensemble's observation map output to the correct directory. Takes an output directory, iteration number, and the ensemble output to save.

update_ensemble(output_dir::AbstractString, iteration, prior)

Update the EnsembleKalmanProcess object and save the parameters for the next iteration.

update_ensemble!(ekp, G_ens, output_dir, iteration, prior)

Update an EKP object with data G_ens, saving the object and parameters for the next iteration to disk.

observation_map_and_update!(
    ekp,
    output_dir,
    iteration,
    prior,
    interface,
)

Compute the observation map and update the given EKP object.

get_prior(param_dict::AbstractDict; names = nothing)
get_prior(prior_path::AbstractString; names = nothing)

Construct the combined prior distribution from a param_dict or a TOML configuration file specified by prior_path.

If names is provided, only those parameters are used.

Examples

prior = ClimaCalibrate.get_prior("prior.toml")
subset = ClimaCalibrate.get_prior("prior.toml"; names = ["coefficient_a"])

get_param_dict(distribution; names)

Generate a dictionary for parameters based on the specified distribution, assumed to be of floating-point type. If names is not provided, the distribution's names will be used.

path_to_iteration(output_dir, iteration)

Return the path to the directory for a given iteration within the specified output directory.

path_to_ensemble_member(output_dir, iteration, member)

Return the path to an ensemble member's directory for a given iteration and member number.

This is where a forward model should write its output.

Examples

ClimaCalibrate.path_to_ensemble_member("output", 3, 7)
# "output/iteration_003/member_007"

See also parameter_path, checkpoint_path, path_to_model_log, path_to_iteration.

path_to_model_log(output_dir, iteration, member)

Return the path to an ensemble member's forward model log for a given iteration and member number.

parameter_path(output_dir, iteration, member)

Return the path to an ensemble member's parameter file.

ClimaCalibrate writes this file before the forward model runs. It is TOML in the format ClimaParams.jl reads, so it can be parsed with TOML.parsefile or passed straight to ClimaParams.create_toml_dict.

Examples

ClimaCalibrate.parameter_path("output", 3, 7)
# "output/iteration_003/member_007/parameters.toml"

checkpoint_path(output_dir, iteration, member)

Return the path to an ensemble member's checkpoint file.

load_latest_ekp(output_dir)

Return the most recent EnsembleKalmanProcess struct from the given output directory.

Returns nothing if no EKP structs are found.

load_ekp_struct(output_dir, iteration)

Return the EnsembleKalmanProcess struct for a completed iteration.

A run that was killed while writing this file leaves it truncated, which JLD2 reports as a missing superblock without naming the file, so the error says which one it is and what to do about it.

ekp_path(output_dir, iteration)

Return the path to the serialized EnsembleKalmanProcess struct file for a given iteration.

save_eki_and_parameters(ekp, output_dir, iteration, prior)

Save the EnsembleKalmanProcess state and each ensemble member's parameters for iteration.

Helper for initialize and update_ensemble.

model_started(output_dir, iteration, member)

Return true if the ensemble member's forward model started but did not finish.

This is how an interrupted run is distinguished from one that never began.

model_completed(output_dir, iteration, member)

Return true if the ensemble member's forward model finished successfully.

Returns false when no checkpoint exists, which is also the case for a member that has not been run yet.

write_model_started(output_dir, iteration, member)

Record that an ensemble member's forward model is about to run.

The checkpoint is overwritten by write_model_completed once the model finishes, so a member left in the "started" state is one that was interrupted.

write_model_completed(output_dir, iteration, member)

Record that an ensemble member's forward model finished successfully, so that a restart skips it.

The forward model itself calls this, which is why it is exported: an HPCBackend job script runs it as its last statement.

ClimaCalibrate.EKPUtils

Helpers for working with EnsembleKalmanProcesses.jl objects that do not depend on the rest of ClimaCalibrate.

Covers building minibatchers and ObservationSeries from a vector of samples (minibatcher_over_samples, observation_series_from_samples), allocating a G ensemble matrix of the right shape (g_ens_matrix), and looking up which observations an iteration is being scored against (get_observations_for_nth_iteration).

minibatcher_over_samples(n_samples, batch_size)

Create a FixedMinibatcher that divides n_samples into batches of size batch_size.

If n_samples is not divisible by batch_size, the remaining samples are dropped and a warning is emitted. batch_size larger than n_samples leaves no minibatch at all, which is an error: EKP divides by the number of minibatches in an epoch.

Examples

minibatcher = ClimaCalibrate.minibatcher_over_samples(10, 5)

See also observation_series_from_samples.

minibatcher_over_samples(samples, batch_size)

Create a FixedMinibatcher that divides a vector of samples into batches of size batch_size.

If the number of samples is not divisible by batch_size, the remaining samples will be dropped.

observation_series_from_samples(samples, batch_size, names = nothing)

Create an EKP.ObservationSeries from a vector of EKP.Observation samples.

If the number of samples is not divisible by batch_size, the remaining samples are dropped.

Examples

obs_series = ClimaCalibrate.observation_series_from_samples(observations, 5)

See also minibatcher_over_samples.

get_observations_for_nth_iteration(obs_series::EKP.ObservationSeries, N)

For the Nth iteration, return a vector of the observation(s) being processed.

get_metadata_for_nth_iteration(obs_series::EKP.ObservationSeries, N)

For the Nth iteration, return a vector of the metadata of the observation(s) being processed.

g_ens_matrix(eki::EKP.EnsembleKalmanProcess{FT}) where {FT <: AbstractFloat}

Construct a G ensemble matrix of type FT, filled with NaN, sized for the current iteration's observation and ensemble.

Starting from NaN means a member whose forward model failed is left as NaN, which is how EKP is told to ignore it.

Examples

G_ensemble = ClimaCalibrate.g_ens_matrix(ekp)

residual(
    ekp::EKP.EnsembleKalmanProcess;
    N = EKP.get_N_iterations(ekp),
    ignore_nan = true,
)

Return the normalized residual (mean(G) - obs) / σ for the Nth iteration, where σ is the square root of the diagonal of the observation noise covariance.

If ignore_nan = true, then the mean of the G ensemble at each index is computed over the ensemble members that are not NaN.

See the Visualization documentation for how to interpret the residual.

ClimaCalibrate.SampleBuilder

Turn ClimaAnalysis.OutputVars into a matrix of flattened samples.

This is the first of the two steps in building an observation: SampleBuilder produces a SampleCollection, and ClimaCalibrate.ObservationRecipe then estimates a noise covariance from it and assembles the EKP.Observation.

Each column of the collection is one sample, and each carries the metadata needed to reconstruct the OutputVars later.

Requires ClimaAnalysis to be loaded.

SampleCollection

An object for storing a collection of samples and their associated metadata as matrices.

The collection of samples is represented as a matrix. Each column of the matrix of samples represents one sample which is a vertical concatenation of one or more flattened ClimaAnalysis.OutputVars. Each column of metadata holds the corresponding ClimaAnalysis.Var.Metadata, one for each ClimaAnalysis.OutputVar.

For each row of the matrix of the metadata, for dimensions that are not ignored, it is guaranteed that

1. the short names are the same,
2. the flattened vector size are the same,
3. the units are the same,
4. the dimensions are the same,
5. the number of dimensions are the same,
6. the dimension units are the same,
7. the dimension values are the same,
8. the coordinates where the NaNs are dropped are the same.

build_samples(var::OutputVar; FT = Float32, dims = ("longitude", "latitude", "pressure_level", "x", "y", "z", "time"))

Return a SampleCollection with a single sample consisting of var.

The matrix of samples has element type FT (defaults to Float32). OutputVars are flattened in the order of dims.

build_samples(
    var_sample::Vector;
    FT = Float32,
    dims = ("longitude", "latitude", "pressure_level", "x", "y", "z", "time"),
)

Return a SampleCollection with a single sample consisting of the OutputVars in var_sample.

The matrix of samples has element type FT (defaults to Float32). OutputVars are flattened in the order of dims.

build_samples(
    var_samples::Matrix;
    FT = Float32,
    dims = ("longitude", "latitude", "pressure_level", "x", "y", "z", "time"),
    ignore_dims = (),
)

Return a SampleCollection from var_samples.

The matrix of samples has element type FT (defaults to Float32).

The OutputVars are flattened in the order of dims. When validating between OutputVars, checks for dimensions in ignore_dims are skipped.

It is the user's responsibility to ensure that ignoring those dimensions is appropriate. For example, build_samples_by_times ignores the time dimension because each sample is windowed to a different time range.

build_samples_by_times(
    var::OutputVar,
    time_ranges;
    FT = Float32,
    dims = ("longitude", "latitude", "pressure_level", "x", "y", "z", "time"),
)

Generate samples from a single OutputVar by windowing its time dimension with time_ranges (one sample per time range).

The matrix of samples has element type FT (defaults to Float32).

build_samples_by_times(
    vars::Vector,
    time_ranges;
    FT = Float32,
    dims = ("longitude", "latitude", "pressure_level", "x", "y", "z", "time"),
)

Generate samples from a vector of OutputVars by windowing the times of the OutputVars in vars using time_ranges.

Each sample has all OutputVars in vars, but the times are windowed according to time_ranges. The matrix of samples has element type FT (defaults to Float32).

num_samples(sample_collection::SampleCollection)

Return the number of samples in sample_collection.

reconstruct_col(sample_collection::SampleCollection, i::Integer)

Reconstruct the ith column of sample_collection as a vector of OutputVars.

get_samples(sample_collection::SampleCollection)

Return the matrix of samples stored in sample_collection.

Mutating this matrix also mutates the matrix in sample_collection.

get_metadata(sample_collection::SampleCollection)

Return the matrix of ClimaAnalysis.Var.Metadata in sample_collection.

Mutating this matrix also mutates the matrix in sample_collection.

var_indices(sample_collection::SampleCollection)

Return a vector of ranges, where each range contains the indices of the rows of the matrix of samples that belong to each OutputVar.

The ranges are computed from the metadata of the first sample, which is valid for every sample, since all samples share the same variables.

ClimaCalibrate.ObservationRecipe

Estimate a noise covariance from a SampleCollection and build an EKP.Observation from it.

Three estimators are available: ScalarCovariance for a multiple of the identity, SeasonalDiagonalCovariance for the per-season variance across years, and SVDplusDCovariance for a low-rank sample covariance plus a diagonal term. All of them take samples built by ClimaCalibrate.SampleBuilder.

Also reconstructs the flattened vectors back into OutputVars (reconstruct_vars, reconstruct_g), so a calibration's observations and forward map output can be inspected.

Requires ClimaAnalysis to be loaded.

AbstractCovarianceEstimator

An object that estimates the noise covariance matrix from the samples in a SampleCollection.

AbstractCovarianceEstimator have to provide one function, ObservationRecipe.covariance.

The function has to have the signature

ObservationRecipe.covariance(
    covar_estimator::AbstractCovarianceEstimator,
    sample_collection,
)

and return a noise covariance matrix. The SampleCollection carries the matrix of flattened samples and their metadata. The covariance matrix does not depend on which sample is chosen as the observation.

Subtypes:

• ScalarCovariance: a multiple of the identity.
• SeasonalDiagonalCovariance: the per-season variance across samples.
• SVDplusDCovariance: a low-rank sample covariance plus a diagonal term.

ScalarCovariance <: AbstractCovarianceEstimator

Covariance estimator contain the necessary information to construct the scalar covariance matrix.

FT1 and FT2 are the element types of scalar and min_cosd_lat.

ScalarCovariance(;
    scalar = 1.0,
    use_latitude_weights = false,
    min_cosd_lat = 0.1,
)

Create a ScalarCovariance which specifies how the covariance matrix should be formed. When used with ObservationRecipe.observation or ObservationRecipe.covariance, return a Diagonal matrix.

Keyword Arguments

• scalar: Scalar value to multiply the identity matrix by.

• use_latitude_weights: If true, then latitude weighting is applied to the covariance matrix. Latitude weighting is multiplying the values along the diagonal of the covariance matrix by (1 / max(cosd(lat), min_cosd_lat)). See the keyword argument min_cosd_lat for more information.

• min_cosd_lat: Control the minimum latitude weight when use_latitude_weights is true. The weight is 1 / max(cosd(lat), min_cosd_lat), so this is the largest weight any point can be given, 1 / min_cosd_lat. Without it the weight grows without bound toward the poles, where cosd(lat) reaches zero, and the diagonal entries span so many orders of magnitude that the covariance is badly conditioned.

SeasonalDiagonalCovariance <: AbstractCovarianceEstimator

Covariance estimator that contain the necessary information to construct a diagonal matrix whose diagonal is the per-season variance across the samples of a SampleCollection.

SeasonalDiagonalCovariance(;
    model_error_scale = 0.0,
    regularization = 0.0,
    use_latitude_weights = false,
    min_cosd_lat = 0.1,
)

Create a SeasonalDiagonalCovariance which specifies how the covariance matrix should be formed. When used with ObservationRecipe.observation or ObservationRecipe.covariance, return a Diagonal matrix.

The samples used to compute the covariance matrix come from the SampleCollection, where each sample is one year of seasonal statistics. NaNs are ignored when computing the seasonal variance.

Keyword Arguments

• model_error_scale: Noise from the model error added to the covariance matrix. This is (model_error_scale * seasonal_mean).^2, where seasonal_mean is the seasonal mean for each of the quantity for each of the season (DJF, MAM, JJA, SON).

• regularization: A diagonal matrix of the form regularization * I is added to the covariance matrix. It is added before latitude weighting, so with use_latitude_weights = true the effective regularization varies with latitude, unlike SVDplusDCovariance, which adds it afterwards.

• use_latitude_weights: If true, then latitude weighting is applied to the covariance matrix. Latitude weighting is multiplying the values along the diagonal of the covariance matrix by (1 / max(cosd(lat), min_cosd_lat)). See the keyword argument min_cosd_lat for more information.

• min_cosd_lat: Control the minimum latitude weight when use_latitude_weights is true. The weight is 1 / max(cosd(lat), min_cosd_lat), so this is the largest weight any point can be given, 1 / min_cosd_lat. Without it the weight grows without bound toward the poles, where cosd(lat) reaches zero, and the diagonal entries span so many orders of magnitude that the covariance is badly conditioned.

SVDplusDCovariance <: AbstractCovarianceEstimator

Covariance estimator that returns an EKP.SVDplusD: a low-rank sample covariance plus a diagonal term.

SVDplusDCovariance(;
    model_error_scale = 0.0,
    regularization = 0.0,
    use_latitude_weights = false,
    use_weighted_samples_for_diagonal = true,
    min_cosd_lat = 0.1,
    rank = nothing
)

Create a SVDplusDCovariance which specifies how the covariance matrix should be formed. When used with ObservationRecipe.observation or ObservationRecipe.covariance, return a EKP.SVDplusD covariance matrix.

The samples used to compute the covariance matrix come from the SampleCollection, where each sample is one column.

Keyword Arguments

• model_error_scale: Noise from the model error added to the covariance matrix. This is (model_error_scale * mean(samples, dims = 2)).^2, where mean(samples, dims = 2) is the mean of the samples.

• regularization: If a scalar is used, a diagonal matrix of the form regularization * I is added to the covariance matrix. See QuantileRegularization for another option for regularization.

• use_latitude_weights: If true, then latitude weighting is applied to the covariance matrix. Latitude weighting is multiplying the columns of the matrix of samples by 1 / sqrt(max(cosd(lat), 0.1)). See the keyword argument min_cosd_lat for more information.

• use_weighted_samples_for_diagonal: If true and use_latitude_weights is true, then the diagonal term is computed from the latitude weighted samples. Otherwise, the diagonal term is computed from the samples without latitude weighting. This has no effect when use_latitude_weights is false.

• min_cosd_lat: Control the minimum latitude weight when use_latitude_weights is true. The weight is 1 / max(cosd(lat), min_cosd_lat), so this is the largest weight any point can be given, 1 / min_cosd_lat. Without it the weight grows without bound toward the poles, where cosd(lat) reaches zero, and the diagonal entries span so many orders of magnitude that the covariance is badly conditioned.

• rank: Rank of the singular value decomposition (SVD). If nothing is passed in, then the rank is automatically inferred from the data.

SVDplusDCovariance(
    diagonal::AbstractDiagonalTerm;
    use_latitude_weights = false,
    use_weighted_samples_for_diagonal = true,
    min_cosd_lat = 0.1,
    rank = nothing,
)

Create a SVDplusDCovariance whose diagonal matrix is described by the diagonal term diagonal. When used with ObservationRecipe.observation or ObservationRecipe.covariance, return a EKP.SVDplusD covariance matrix.

Passing model_error_scale = x and regularization = y to the keyword constructor is the same as passing diagonal = ModelErrorScaleDiagonal(x) .+ ScalarDiagonal(y). Passing regularization = QuantileRegularization(q) instead is the same as passing diagonal = ModelErrorScaleDiagonal(x) .+ QuantileDiagonal(q, ModelErrorScaleDiagonal(x)).

Keyword Arguments

• use_latitude_weights: If true, then latitude weighting is applied to the covariance matrix. Latitude weighting is multiplying the columns of the matrix of samples by 1 / sqrt(max(cosd(lat), 0.1)). See the keyword argument min_cosd_lat for more information.

• use_weighted_samples_for_diagonal: If true and use_latitude_weights is true, then the diagonal term is computed from the latitude weighted samples. Otherwise, the diagonal term is computed from the samples without latitude weighting. This has no effect when use_latitude_weights is false.

• min_cosd_lat: Control the minimum latitude weight when use_latitude_weights is true. The weight is 1 / max(cosd(lat), min_cosd_lat), so this is the largest weight any point can be given, 1 / min_cosd_lat. Without it the weight grows without bound toward the poles, where cosd(lat) reaches zero, and the diagonal entries span so many orders of magnitude that the covariance is badly conditioned.

• rank: Rank of the singular value decomposition (SVD). If nothing is passed in, then the rank is automatically inferred from the data.

QuantileRegularization

Regularization using the quantile of the model error scale for each OutputVar.

The same quantile is used for each OutputVar when making the observation.

This is used for the SVDplusDCovariance matrix.

Examples

In the example below, a regularization using the 0.05 quantile of the model error scale for each variable is initialized.

qtl_regularization = QuantileRegularization(0.05)

covariance(
    covar_estimator::ScalarCovariance,
    sample_collection::SampleCollection,
)

Compute the scalar covariance matrix.

The data in the matrix of samples in sample_collection is ignored.

covariance(
    covar_estimator::SeasonalDiagonalCovariance,
    sample_collection::SampleCollection,
)

Compute the diagonal covariance matrix of seasonal quantities from the samples in sample_collection.

The diagonal entries are the per-entry variance across the samples of sample_collection. This is the variance of each season across years with NaNs ignored. Each sample must represent the same sequence of seasons with one time slice per season. At least two samples (years) are required to estimate a variance.

covariance(
    covar_estimator::SVDplusDCovariance,
    sample_collection::SampleCollection,
)

Compute the EKP.SVDplusD covariance matrix from the samples in sample_collection.

observation(
    covar_estimator::AbstractCovarianceEstimator,
    sample_collection::SampleCollection,
    i::Integer;
    name = nothing,
)

Return an EKP.Observation with the ith sample of sample_collection as the observation, a covariance matrix defined by covar_estimator, name determined from the short names of the observation, and metadata.

short_names(obs::EKP.Observation)

Get the short names of the variables from the metadata in the EKP.Observation.

If the short name is not available, then nothing is returned instead.

reconstruct_g(ekp::EKP.EnsembleKalmanProcess, it::Integer)

Reconstruct the G ensemble matrix of the itth iteration as a matrix of OutputVars.

reconstruct_g_mean(ekp::EKP.EnsembleKalmanProcess, it::Integer)

Reconstruct the mean forward model evaluation at the itth iteration as a vector of OutputVars.

reconstruct_g_mean_final(ekp::EKP.EnsembleKalmanProcess)

Reconstruct the mean forward model evaluation at the last iteration as a vector of OutputVars.

reconstruct_diag_cov(obs::EKP.Observation)

Reconstruct the diagonal of the covariance matrix in obs as a vector of OutputVars.

This function only supports observations that contain diagonal covariance matrices.

The units of the reconstructed OutputVars are squared, since the diagonal of the covariance matrix contains variances.

reconstruct_vars(obs::EKP.Observation)

Reconstruct the OutputVars from the samples in obs.

reconstruct_residual(
    ekp::EKP.EnsembleKalmanProcess,
    it::Integer;
    ignore_nan = true,
)

Reconstruct the normalized residual (mean(G) - obs) / σ of the itth iteration as a vector of OutputVars, where σ is the square root of the diagonal of the observation noise covariance.

If ignore_nan = true, then the mean of the G ensemble at each index is computed over the ensemble members that are not NaN.

The units of the reconstructed OutputVars are empty, since the residual is normalized by σ.

seasonally_aligned_yearly_sample_date_ranges(var::OutputVar)

Generate sample dates that conform to a seasonally aligned year from dates(var).

A seasonally aligned year is defined to be from December to November of the following year.

This function is useful for finding the sample dates of samples consisting of all four seasons in a single year. For example, one can pass these date ranges to SampleBuilder.build_samples_by_times to build the samples for a SVDplusDCovariance or SeasonalDiagonalCovariance.

AbstractDiagonalTerm

A description of how to build the diagonal matrix from a SampleCollection. It is not a matrix itself.

A diagonal term is lazy because the covariance recipe may transform the samples before the diagonal is computed. The term only records what to compute, and compute_diagonal builds the matrix once the samples are available.

To define a custom diagonal term, subtype AbstractDiagonalTerm and implement a method of compute_diagonal for it.

compute_diagonal(diagonal_term::AbstractDiagonalTerm, sample_collection)

Compute the diagonal matrix described by diagonal_term from the samples in sample_collection.

ScalarDiagonal{FT <: AbstractFloat} <: AbstractDiagonalTerm

A diagonal term whose entries are constant for each corresponding variable.

Example

scalars = [1e-6, 1e-4]
ScalarDiagonal(scalars)

In the example, scalars[i] fills every diagonal entry belonging to the ith variable. If scalars is length 1, then the same constant is used for all variables.

ScalarDiagonal(scalar::AbstractFloat)

A diagonal term where every entry is scalar.

ModelErrorScaleDiagonal{FT <: AbstractFloat} <: AbstractDiagonalTerm

A diagonal term whose entries are (scale * mean)^2, where mean is the mean of each entry across the samples, ignoring NaNs, and scale is the model error scale of the variable that the entry belongs to.

Example

model_error_scales = [0.05, 0.1]
ModelErrorScaleDiagonal(model_error_scales)

In the example, model_error_scales[i] scales the mean of every entry belonging to the ith variable. Scales must not be negative. If model_error_scales is length 1, then the same model error scale is used for all variables.

ModelErrorScaleDiagonal(model_error_scale::AbstractFloat)

A diagonal term whose entries are (model_error_scale * mean)^2, using the same scale for every variable.

VarianceDiagonal <: AbstractDiagonalTerm

A diagonal term whose entries are the variance of each row of the sample matrix, taken across the samples and ignoring NaNs.

Example

VarianceDiagonal()

QuantileDiagonal <: AbstractDiagonalTerm

A diagonal term whose entries are constant for each variable, where the constant is a quantile of the entries that another diagonal term produces for that variable.

Example

quantiles = [0.5, 0.05]
QuantileDiagonal(quantiles, VarianceDiagonal())

In the example, the variances are computed first. Then, for the ith variable, the quantiles[i] quantile of that variable's variances becomes the value of every diagonal entry belonging to that variable. Quantiles must be in (0, 1]. This is useful for smoothing out or flooring a diagonal term computed from the samples. If quantiles is length 1, then the same quantile is used for all variables.

QuantileDiagonal(
    quantile::AbstractFloat,
    diag_term::AbstractDiagonalTerm,
)

A diagonal term whose entries are constant for each variable, where the constant is the quantile of the entries that diag_term produces for that variable. The same quantile is used for every variable.

BroadcastedDiagonal <: AbstractDiagonalTerm

A diagonal term whose entries are a broadcast expression over other diagonal terms, Diagonal matrices, vectors, and scalars.

Example

2.0 .* ModelErrorScaleDiagonal(0.05) .+ ScalarDiagonal(1e-6)

Once every term is computed, the expression is broadcast over the resulting Diagonal matrices and follows the same broadcasting rules. Operations that keep a Diagonal diagonal, such as 2.0 .* term, term .^ 2, and term1 .+ term2, are supported. Operations that would not, such as term .+ 1e-6, error when the diagonal is computed.

You should not directly construct a BroadcastedDiagonal since it is constructed by any dotted expression containing a diagonal term.

analyze_residual(ekp, iter; n_eigenvectors = 3)

Analyze the model-data residual G(u) - y at iteration iter using the leading eigenvectors of the observational noise covariance.

Projecting the residual onto those eigenvectors and normalizing by the corresponding eigenvalues turns it into z-scores, so a value much larger than one means the residual has structure the noise model does not account for. A high structured energy in one variable points at that variable's observation or its part of the observation map.

The noise covariance comes from EKP.get_obs_noise_cov with build = false, so this works with any StructuredMatrix EKP supports (SVD, Diagonal, SVDplusD), not only SVDplusD.

Arguments

• ekp: The EKP.EnsembleKalmanProcess to analyze.
• iter: Which iteration's residual to analyze.

Keyword Arguments

• n_eigenvectors = 3: How many leading eigenvectors to project onto.

Returns

A NamedTuple with:

• normalized_projections: (n_eigenvectors × n_variables) matrix of z-scores per variable [-]. Values much greater than one indicate mismatch beyond the noise model.
• structured_energy: Normalized whitened energy across all variables [-], approximately one under the noise model.
• structured_energy_by_variable: How that energy is split between variables [-]. Under the noise model the entries sum to structured_energy, so compare them against each other rather than against one.
• residual_norm_by_variable: norm(diff[rᵥ]) for each variable.
• metadata: A ClimaAnalysis.Var.Metadata per variable, in the same order as the columns of normalized_projections and the elements of structured_energy_by_variable and residual_norm_by_variable.

Examples

import ClimaAnalysis   # required
result = ClimaCalibrate.analyze_residual(ekp, 3; n_eigenvectors = 3)
result.structured_energy

compute_structured_energy(projections)

Given the matrix of normalized projections from compute_normalized_projections, compute the total structured energy in the whitened space:

energy = (1/n_eig) * ∑ᵢ zᵢ²,   where zᵢ = ∑ᵥ projections[i, v] = aᵢ / √λᵢ

zᵢ is the global whitened projection onto eigenvector i. Under the noise model, zᵢ ~ N(0, 1), so energy ≈ 1 is consistent with noise. Values >> 1 indicate mismatch beyond what the structured noise explains. Values << 1 suggest overfitting to noise or an overestimated noise covariance.

compute_structured_energy_by_variable(projections)

Given the matrix of normalized projections from compute_normalized_projections, compute the per-variable structured energy in the whitened space:

energy_v = (1/n_eig) * ∑ᵢ projections[i, v]²

Returns a vector of length n_variables. These are shares of the total: under the noise model they sum to the same value that compute_structured_energy reports, which is approximately one, so each entry is well below one on its own. A variable holding most of the total is the one carrying the structure, and a variable near zero is one the leading eigenvectors barely touch. Compare variables against each other rather than against one.

See also analyze_residual.

compute_normalized_projections(diff, eigvectors, eigvalues, ranges)

For each eigenvector i and variable v, compute the variable-specific contribution to the normalized projection aᵢᵛ / √λᵢ, where aᵢᵛ = vᵢ[rᵥ]ᵀ diff[rᵥ].

This is a z-score: values >> 1 indicate model-data mismatch beyond noise; values ≤ 1 are consistent with structured noise.

Returns a Matrix of shape (n_eigenvectors, n_variables). ranges is a vector of index ranges, one per variable, giving the observation indices for that variable. Pass the result to compute_structured_energy or compute_structured_energy_by_variable for scalar summaries.

See also analyze_residual.

ClimaCalibrate.EnsembleBuilder

Assemble the G ensemble matrix from ClimaAnalysis.OutputVars.

GEnsembleBuilder reads the metadata off the observations in an EnsembleKalmanProcess and works out where each variable belongs in the matrix, so index ranges do not have to be tracked by hand. It validates each OutputVar against the observation it is filling in, checking short name, units, dimension names, dimension units, and dimension values. A mismatch between model output and observations raises an error instead of being calibrated against silently.

Requires ClimaAnalysis to be loaded.

GEnsembleBuilder{FT <: AbstractFloat}

An object to help build G ensemble matrix by using the metadata stored in the EKP.EnsembleKalmanProcess object. Metadata must come from ClimaAnalysis.

GEnsembleBuilder takes in preprocessed OutputVars and automatically constructs the corresponding G ensemble matrix for the current iteration of the calibration.

FT is the element type of the G ensemble matrix.

GEnsembleBuilder(ekp::EKP.EnsembleKalmanProcess{FT})
    where {FT <: AbstractFloat}

Construct a GEnsembleBuilder where the element type of the G ensemble matrix is FT.

EnsembleBuilder.fill_g_ens_col!(
    g_ens_builder::GEnsembleBuilder,
    col_idx,
    var::OutputVar;
    checkers = (),
    verbose = false
)

Fill the col_idxth column of the G ensemble matrix from var using g_ens_builder. If it was successful, return true, otherwise, return false.

It is assumed that the times or dates of a single OutputVar is a superset of the times or dates of one or more metadata in the minibatch.

This function relies on the short names in the metadata. This function will not behave as intended if the short names are mislabeled or not present.

Furthermore, this function assumes that all observations are generated using ObservationRecipe.observation which guarantees that the metadata exists and the correct placement of metadata.

EnsembleBuilder.fill_g_ens_col!(
    g_ens_builder::GEnsembleBuilder,
    col_idx,
    val::AbstractFloat
)

Fill the col_idxth column of the G ensemble matrix with val.

This returns true.

This is useful if you want to completely fill a column of a G ensemble matrix with NaNs if a simulation crashed.

EnsembleBuilder.is_complete(g_ens_builder::GEnsembleBuilder)

Return true if all the entries of the G ensemble matrix are filled out and false otherwise.

EnsembleBuilder.get_g_ensemble(g_ens_builder::GEnsembleBuilder)

Return the G ensemble matrix from g_ens_builder.

If the G ensemble matrix is not completed, then a warning is thrown. See ClimaCalibrate.EnsembleBuilder.is_complete to check if the G ensemble matrix is completely filled out.

ranges_by_short_name(g_ens_builder::GEnsembleBuilder, short_name)

Return a vector of ranges for the G ensemble matrix that correspond with the short name.

metadata_by_short_name(g_ens_builder::GEnsembleBuilder, short_name)

Return a vector of metadata that correspond with short_name.

missing_short_names(g_ens_builder::GEnsembleBuilder, col_idx)

Return a set of the short names of the metadata that are not filled out for the col_idxth column of g_ens_builder.

ClimaCalibrate.Checker

The validation checks that ClimaCalibrate.EnsembleBuilder runs when matching model output against an observation.

Each checker answers one question about an OutputVar and the Metadata it is being matched to: do the short names agree, the units, the dimension names, their units, their values. GEnsembleBuilder runs a default set of them; pass others with the checkers keyword argument, or define your own by subtyping AbstractChecker and implementing check.

Requires ClimaAnalysis to be loaded.

AbstractChecker

An object that performs validation checks between the simulation data and metadata from observational data. This is used by GEnsembleBuilder to validate OutputVars from simulation data against the Metadata in the observations in the EnsembleKalmanProcess object.

An AbstractChecker must implement the Checker.check function.

The function must have the signature:

import ClimaCalibrate.Checker
Checker.check(::YourChecker,
              var::OutputVar,
              metadata::Metadata;
              data = nothing,
              verbose = false)

and return true or false.

Subtypes:

• ShortNameChecker: the short names agree.
• DimNameChecker: the dimension names agree.
• DimUnitsChecker: the dimension units agree.
• UnitsChecker: the units agree.
• DimValuesChecker: the dimension values agree.
• SequentialIndicesChecker: the matched dates are sequential.
• SignChecker: the proportion of positive values agrees.

struct ShortNameChecker <: AbstractChecker end

A struct that checks the short name between simulation data and metadata.

struct DimNameChecker <: AbstractChecker end

A struct that checks the dimension names between simulation data and metadata.

struct DimUnitsChecker <: AbstractChecker end

A struct that checks the units of the dimensions between simulation data and metadata.

struct UnitsChecker <: AbstractChecker end

A struct that checks the units between the simulation data and metadata.

struct DimValuesChecker <: AbstractChecker end

A struct that checks the values of the dimensions between the simulation data and metadata.

struct SequentialIndicesChecker <: AbstractChecker end

A struct that checks that the indices of the dates of the simulation data corresponding to the dates of the metadata is sequential.

struct SignChecker{FT <: AbstractFloat} <: AbstractChecker

A struct that checks that the proportion of positive values in the simulation data and observational data is roughly the same.

To change the default threshold of 0.05, you can pass a float to SignChecker.

import ClimaCalibrate
sign_checker = ClimaCalibrate.Checker.SignChecker(0.01)

check(checker::AbstractChecker,
      var,
      metadata;
      data = nothing,
      verbose = false)

Return true if the check passes, false otherwise.

If verbose=true, then provides information for why a check did not succeed.

Checker.check(
    ::ShortNameChecker,
    var::OutputVar,
    metadata::Metadata;
    data = nothing,
    verbose = false,
)

Return true if var and metadata have the same short name, false otherwise.

Checker.check(
    ::DimNameChecker,
    var::OutputVar,
    metadata::Metadata;
    data = nothing,
    verbose = false,
)

Return true if var and metadata have the same dimensions, false otherwise.

Checker.check(
    ::DimUnitsChecker,
    var::OutputVar,
    metadata::Metadata;
    data = nothing,
    verbose = false,
)

Return true if the units of the dimensions in var and metadata are the same, false otherwise. This function assumes var and metadata have the same dimensions.

Checker.check(
    ::UnitsChecker,
    var::OutputVar,
    metadata::Metadata;
    data = nothing,
    verbose = false,
)

Return true if var and metadata have the same units, false otherwise.

Checker.check(
    ::DimValuesChecker,
    var::OutputVar,
    metadata::Metadata;
    data = nothing,
    verbose = false,
)

Return true if the values of the dimensions in var and metadata are compatible for the purpose of filling out the G ensemble matrix, false otherwise.

The nontemporal dimensions are compatible if the values are approximately the same. The temporal dimensions are compatible if the temporal dimension of metadata is a subset of the temporal dimension of var.

Checker.check(
    ::SequentialIndicesChecker,
    var::OutputVar,
    metadata::Metadata;
    data = nothing,
    verbose = false,
)

Return true if the dates of metadata map to consecutive indices within the dates of var, false otherwise.

Checker.check(
    ::SignChecker,
    var::OutputVar,
    metadata::Metadata;
    data,
    verbose = false,
)

Return true if the proportion of positive values in var, flattened with observation metadata, is within the threshold defined in SignChecker of the proportion of positive values in data, false otherwise.

ClimaCalibrate.Visualization

Makie plots of a calibration's forward map output against its observations.

plot_g draws all ensemble members, plot_g_mean the ensemble mean, and plot_obs the observations, so the three can be overlaid to see whether the ensemble is bracketing the data.

Requires a Makie backend such as CairoMakie to be loaded.

plot_g

Plot members of the G ensemble matrix as line plots.

If the iter keyword argument is not passed, then this plots the last completed G ensemble matrix. Otherwise, it plots the iterth G ensemble matrix.

All Makie keyword arguments compatible with Makie.lines are also compatible with this function.

Examples

import CairoMakie
fig, ax, _ = ClimaCalibrate.Visualization.plot_g(ekp; color = (:grey, 0.3))
ClimaCalibrate.Visualization.plot_obs!(ax, ekp; color = :black)

See also plot_g_mean, plot_obs.

plot_g!

This is the mutating variant of the plotting function plot_g.

plot_g_mean

Plot mean forward map evaluation as a line plot.

If the iter keyword argument is not passed, then this plots the last mean forward map evaluation. Otherwise, it plots the iterth mean forward map evaluation.

All Makie keyword arguments compatible with Makie.lines are also compatible with this function.

plot_g_mean!

This is the mutating variant of the plotting function plot_g_mean.

plot_obs

Plot the observations as a line plot.

If the iter keyword argument is not passed, then this plots the observations of the last iteration. Otherwise, it plots the observations of the iterth iteration.

All Makie keyword arguments compatible with Makie.lines are also compatible with this function.

plot_obs!

This is the mutating variant of the plotting function plot_obs.

plot_residual

Plot the normalized residual (mean(G) - obs) / σ as a scatter plot, where σ is the square root of the diagonal of the observation noise covariance. Reference lines are drawn at zero and at plus and minus one and two σ; pass reference_lines = false to not draw them.

If the iter keyword argument is not passed, then this plots the residual of the last iteration. Otherwise, it plots the residual of the iterth iteration.

If the keyword argument ignore_nan = true, then the mean of the G ensemble at each index is computed over the ensemble members that are not NaN.

All Makie keyword arguments compatible with Makie.scatter are also compatible with this function.

plot_residual!

This is the mutating variant of the plotting function plot_residual.