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/dev/.

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.

Fields

• hpc_config: Scheduler directives, modules, and environment variables for each ensemble member's job. See SlurmConfig.
• job_records: The jobs submitted with this backend, in submission order.
• failure_rate: The fraction of an iteration's ensemble members that may fail before the calibration is halted [-].
• job_timeout: How long (in seconds) an iteration waits for a running job before giving up [s]. The default is 86400 (24 hours).

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.

Fields

• hpc_config: Scheduler directives, modules, and environment variables for each ensemble member's job. See SlurmConfig.
• job_records: The jobs submitted with this backend, in submission order.
• failure_rate: The fraction of an iteration's ensemble members that may fail before the calibration is halted [-].
• job_timeout: How long (in seconds) an iteration waits for a running job before giving up [s]. The default is 86400 (24 hours).

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.

Fields

• hpc_config: Scheduler directives, modules, and environment variables for each ensemble member's job. See SlurmConfig.
• job_records: The jobs submitted with this backend, in submission order.
• failure_rate: The fraction of an iteration's ensemble members that may fail before the calibration is halted [-].
• job_timeout: How long (in seconds) an iteration waits for a running job before giving up [s]. The default is 86400 (24 hours).

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.

Fields

• hpc_config: Scheduler directives, modules, and environment variables for each ensemble member's job. See PBSConfig.
• job_records: The jobs submitted with this backend, in submission order.
• failure_rate: The fraction of an iteration's ensemble members that may fail before the calibration is halted [-].
• job_timeout: How long (in seconds) an iteration waits for a running job before giving up [s]. The default is 86400 (24 hours).

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.

Fields

• backend: The backend the job was submitted with.
• id: The scheduler's job ID, an Int64 for Slurm and a String for PBS.
• job_script: The script that was submitted. write_job_script writes it to a file, which shows what the scheduler was asked to run.

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)

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 and NaNStatistics 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.

FT is the element type of the samples. Access the samples and metadata with get_samples and get_metadata.

build_samples(vars; FT = Float32, dims = ...)

Build a SampleCollection from ClimaAnalysis.OutputVars.

Accepts a single OutputVar or a Vector of them (one sample made of one or more variables), or a Matrix whose rows are variables and whose columns are samples. The Matrix method also takes ignore_dims, to exclude dimensions from the compatibility checks between samples.

Each variable is flattened in the order given by dims, dropping NaNs, and the same coordinates must be dropped in all samples.

Examples

import ClimaAnalysis
samples = ClimaCalibrate.SampleBuilder.build_samples([ta, hus]; FT = Float64)

See also build_samples_by_times.

build_samples_by_times(vars, time_ranges; FT = Float32, dims = ...)

Build a SampleCollection by windowing vars into one sample per time range.

Each element of time_ranges is a (start, stop) pair of dates or times, and becomes one column of the collection. The time dimension is excluded from the between-sample compatibility checks, since each sample covers a different span.

Windows should not overlap: samples that share time slices are correlated, which biases a covariance estimated from them.

Examples

import ClimaAnalysis, Dates
ranges = [
    (Dates.DateTime(y, 12, 1), Dates.DateTime(y + 1, 9, 1)) for y in 2007:2015
]
samples = ClimaCalibrate.SampleBuilder.build_samples_by_times([ta], ranges)

See also build_samples.

num_samples(sample_collection)

Return the number of samples (columns) in a SampleCollection.

Requires ClimaAnalysis and NaNStatistics to be loaded.

reconstruct_col(sample_collection, i)

Return the ith sample as a vector of ClimaAnalysis.OutputVars.

This undoes the flattening that build_samples applied, so a sample can be inspected or plotted.

Requires ClimaAnalysis and NaNStatistics to be loaded.

get_samples(sample_collection)

Return the matrix of flattened samples, one sample per column.

Requires ClimaAnalysis and NaNStatistics to be loaded.

get_metadata(sample_collection)

Return the matrix of ClimaAnalysis.Var.Metadata, one entry per variable per sample.

Requires ClimaAnalysis and NaNStatistics to be loaded.

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 and NaNStatistics 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 that returns a multiple of the identity.

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

Fields

• scalar: Scalar to multiply the identity matrix by.
• use_latitude_weights: Whether to apply latitude weighting.
• min_cosd_lat: The smallest cosd(lat) used in the latitude weight, which caps the weight at 1 / 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 whose diagonal is the per-season variance across the samples of a SampleCollection.

FT1, FT2, and FT3 are the element types of model_error_scale, regularization, and min_cosd_lat.

Fields

• model_error_scale: A model error scale term added to the diagonal of the covariance matrix.
• regularization: A regularization term added to the diagonal of the covariance matrix.
• use_latitude_weights: Whether to apply latitude weighting.
• min_cosd_lat: The smallest cosd(lat) used in the latitude weight, which caps the weight at 1 / min_cosd_lat [-].

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 dropped when the samples are built, not here: SampleBuilder removes them while flattening and requires the same coordinates to be dropped in all sample, so a NaN whose position varies between samples is an error rather than something silently ignored.

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.

FT1, FT2, and FT3 are the element types of model_error_scale, regularization, and min_cosd_lat; R is the type of rank.

Fields

• model_error_scale: A model error scale term added to the diagonal of the covariance matrix.
• regularization: A regularization term added to the diagonal of the covariance matrix, either a scalar or a QuantileRegularization.
• use_latitude_weights: Whether to apply latitude weighting.
• min_cosd_lat: The smallest cosd(lat) used in the latitude weight, which caps the weight at 1 / min_cosd_lat [-].
• rank: Rank of the singular value decomposition, or nothing to infer it from the data.

SVDplusDCovariance(;
    model_error_scale = 0.0,
    regularization = 0.0,
    use_latitude_weights = false,
    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.

• 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, sample_collection)

Estimate the observational noise covariance from sample_collection.

The result does not depend on which sample is used as the observation. See ScalarCovariance, SeasonalDiagonalCovariance, and SVDplusDCovariance.

Examples

import ClimaAnalysis
estimator = ClimaCalibrate.ObservationRecipe.SVDplusDCovariance(;
    regularization = 1e-3,
)
covar = ClimaCalibrate.ObservationRecipe.covariance(estimator, samples)

See also observation.

observation(covar_estimator, sample_collection, i; name, covariance)

Build an EKP.Observation from the ith sample of sample_collection, with a noise covariance estimated by covar_estimator.

The observation carries the metadata of its samples, which is what ClimaCalibrate.EnsembleBuilder uses to line model output up with it, and what the reconstruct_* functions use to turn the flattened vectors back into OutputVars.

Examples

import ClimaAnalysis
obs = ClimaCalibrate.ObservationRecipe.observation(estimator, samples, 1)

See also covariance, reconstruct_vars.

short_names(obs)

Return the short names of the variables in an EKP.Observation, in the order they were stacked.

Requires ClimaAnalysis and NaNStatistics to be loaded.

reconstruct_g(ekp, iter)

Return the G ensemble matrix of iteration iter as a matrix of ClimaAnalysis.OutputVars, one row per variable and one column per ensemble member.

Requires ClimaAnalysis and NaNStatistics to be loaded, and observations built by this module.

reconstruct_g_mean(ekp, iter)

Return the mean forward map evaluation of iteration iter as a vector of ClimaAnalysis.OutputVars.

Requires ClimaAnalysis and NaNStatistics to be loaded, and observations built by this module.

reconstruct_g_mean_final(ekp)

Return the mean forward map evaluation of the last completed iteration as a vector of ClimaAnalysis.OutputVars.

Requires ClimaAnalysis and NaNStatistics to be loaded, and observations built by this module.

reconstruct_diag_cov(obs)

Return the diagonal of an observation's noise covariance as a vector of ClimaAnalysis.OutputVars, so the noise can be plotted alongside the data.

Only meaningful for a diagonal covariance. Requires ClimaAnalysis and NaNStatistics to be loaded.

reconstruct_vars(obs)

Return the observation itself as a vector of ClimaAnalysis.OutputVars.

This undoes the flattening that ClimaCalibrate.SampleBuilder applied, so an observation can be plotted or compared against model output.

Requires ClimaAnalysis and NaNStatistics to be loaded.

seasonally_aligned_yearly_sample_date_ranges(var)

Return the (start, stop) date ranges that split var into one sample per seasonal year, starting at December.

Pass the result to SampleBuilder.build_samples_by_times to build the samples that SeasonalDiagonalCovariance expects.

Requires ClimaAnalysis and NaNStatistics to be loaded.

analyze_residual(ekp, iter; n_eigenvectors = 3)

Analyze the model-data residual y - G(u) 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 and NaNStatistics 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)

Create a builder for the G ensemble matrix of ekp.

The builder reads the metadata off the observations that ekp is being scored against, so it knows which rows of the matrix each variable occupies. Fill it column by column with fill_g_ens_col!, then hand the matrix to EKP with get_g_ensemble.

Examples

import ClimaAnalysis
builder = ClimaCalibrate.EnsembleBuilder.GEnsembleBuilder(ekp)
for member in 1:EKP.get_N_ens(ekp)
    for var in preprocess_member_output(member)
        ClimaCalibrate.EnsembleBuilder.fill_g_ens_col!(builder, member, var)
    end
end
G_ensemble = ClimaCalibrate.EnsembleBuilder.get_g_ensemble(builder)

See also is_complete, missing_short_names.

fill_g_ens_col!(builder, col_idx, var; checkers = ..., verbose = false)
fill_g_ens_col!(builder, col_idx, value::AbstractFloat)

Fill the part of ensemble member col_idx's column that var corresponds to, and return whether var was used.

var is matched against the observation metadata by short name, and validated against it by the checkers before anything is written, so model output that does not line up with the observation is rejected instead of being calibrated against silently. Pass verbose = true to have each failed check say why.

The second form fills the member's whole column with a single value, which is how a failed forward model is marked (with NaN).

Requires ClimaAnalysis and NaNStatistics to be loaded.

is_complete(builder)

Return true once all entries of the G ensemble matrix have been filled.

Requires ClimaAnalysis and NaNStatistics to be loaded.

get_g_ensemble(builder)

Return the G ensemble matrix, with NaN wherever nothing was filled in.

Requires ClimaAnalysis and NaNStatistics to be loaded.

ranges_by_short_name(builder, short_name)

Return the row ranges of the G ensemble matrix that short_name occupies, as a vector with one range per observation carrying that short name.

Requires ClimaAnalysis and NaNStatistics to be loaded.

metadata_by_short_name(builder, short_name)

Return the observation metadata the builder matches short_name against, as a vector with one entry per observation carrying that short name.

Requires ClimaAnalysis and NaNStatistics to be loaded.

missing_short_names(builder, col_idx)

Return the short names that ensemble member col_idx is still missing.

This is the first thing to check when is_complete returns false.

Requires ClimaAnalysis and NaNStatistics to be loaded.

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 and NaNStatistics 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 metadata the way the observation was, is within the threshold defined in SignChecker of the proportion of positive values in data, false otherwise.

This check assumes var can be flattened with metadata, which the default checkers establish before it runs.

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.