bug: quantile input equal to 1.0 raises error in transform_to_probit

[mollymwieringa]

Describe the bug

When attempting to regress observation space increments produced by the KQCEF with the KDE_DISTRIBUTION to state variables using probit transforms with the same KDE_DISTRIBUTION, the process crashes. The source of the crash is the inv_cdf function, which reads quantile inputs that are equal to 1.0 as illegal. inv_cdf is called through 2 wrapper functions by probit_or_logit_transform in the to_probit_kde subroutine in transform_to_probit. In the experiments in question, transform_to_probit is called as part of a regression method defined in a user-written program, state_regress_relativefrac.

The expectation for quantiles = 1.0 is that they would be transformed to the upper bound of the system.

Some notes / observations:
The ensemble being transformed at the point of the bug is a posterior ensemble in observation space (values of SIC, updated using the KQCEF). However, incrementing produced increments for every ensemble member equal to 0, so the obs posterior ensemble is identical to the prior, which was transformed immediately before the posterior ensemble using the same function and raised no errors. It is also apparently identical to the obs posterior ensemble produced using the BNRH filter, which was transformed using the BNRH_DISTRIBUTION and raised no errors. As such, I think this is a computational precision error of some kind.

To reproduce the bug: run transform_to_probit on the obs posterior ensemble stored for day_014 of assimilation in test case SC_SIC_aicen, e.g.

transform_to_probit(ens_size, obs_post, dist_for_obs, obs_dist_params, &
 probit_obs_post, .true., obs_bounded_below, obs_bounded_above, &
 obs_lower_bound, obs_upper_bound, ierr)

with
dist_for_obs = KDE_DISTRIBUTION
obs_bounded_below = .true.
obs_bonded_above = .true.
obs_lower_bound = 0.0
obs_upper_bound = 1.0

obs_dist_params needs to be defined as type(distribution_params_type) from the distribution_params_mod

The relevant files can be found in this directory on Derecho:
/glade/u/home/mollyw/work/Projects/fractional-comp/experiment_output/fixed_error_01/SibChuk/SIC_aicen/day_014/

The files themselves are:
(KQCEF obs posterior ensemble): obs_post_ensemble_kde.txt
(obs prior ensemble): obs_prior_ensemble.txt

(BNRHF obs posterior ensemble): obs_post_ensemble_bnrhf.txt
(KQCEF obs increments): obs_increments_kde.txt
(BNRHF obs increments): obs_increments_bnrhf.txt

Error Message

ERROR FROM:
  source : normal_distribution_mod.f90
  routine: inv_cdf
  message:  Illegal Quantile input   1.00000000000000

Which model(s) are you working with?

This is being run in an non-cycling experimental framework with cice-scm (Icepack) output. The model is not advanced or actively run.

Version of DART

Which version of DART are you using?
v11.10.7-339-g5f6be1c1f

Have you modified the DART code?

I have not modified the DART code itself, but I did write a series of programs to explore some new regression approaches for sea ice category assimilation. I also modified a copy of normal_distribution_mod.f90 to use E_MSG instead of E_ERR, which prevented the entire program from exiting, but raised subsequent issues that prevented the program from writing output files (ens_post_*).

The experimental program files can be found in the cice-scm directory of my DART fork (https://github.com/mollymwieringa/DART/tree/fractional_comp/models/cice-scm). They are:

• obs_increments.f90
• state_regression.f90
• fractional_reg_mod.f90

The error arises from the use of the transform_to_probit subroutine in state_regress_probit, called in state_regress_relativefrac (used in state_regression.f90 and defined in fractional_reg_mod.f90).

Build information

I'm running on Derecho with the intel/2024.2.1 compiler

[mollymwieringa]

@hkershaw-brown / @mjs2369 I think Jeff mentioned in standup that this should be assigned to him?

[mjs2369]

@hkershaw-brown / @mjs2369 I think Jeff mentioned in standup that this should be assigned to him?

@mollymwieringa yep, he did say that. I'd go ahead and send him an email pointing him to this issue.

[jlaucar]

The initial problem here does seem to be another in a series of extreme round-off scenarios in which the input value of 1.0 is somehow one bit larger than 1. The current code issues an error for situations where the quantile exceeds 1. Changing the error to a message allows the code to proceed, but presumably it does not pass the next if test. That one says that if the quantile is 1 and the quantity is bounded above, the return should be the upper bound. Instead, it reaches the third if statement which reduces the quantile to the parameter that defines the upper bound value (just less than 1) that will still converge in the inverse cdf.

Is the actual correct thing to do here to treat any value greater than 1 as exactly 1 so that this case would return the upper bound (1 in Molly's case I assume?). Do we need to worry about cases that are much larger than 1. What would be the threshold between roundoff and "much larger".