diff --git a/cfr/reconjob.py b/cfr/reconjob.py index 0ae3e829..6bbeeeda 100644 --- a/cfr/reconjob.py +++ b/cfr/reconjob.py @@ -443,7 +443,7 @@ def crop_clim(self, tag, lat_min=None, lat_max=None, lon_min=None, lon_max=None, def calib_psms(self, ptype_psm_dict=None, ptype_season_dict=None, ptype_clim_dict=None, calib_period=None, - use_predefined_R=False, verbose=False, **kwargs): + use_predefined_R=False, R_floor=None, R_floor_frac=None, verbose=False, **kwargs): ''' Calibrate the PSMs. Args: @@ -452,6 +452,10 @@ def calib_psms(self, ptype_psm_dict=None, ptype_season_dict=None, ptype_clim_dic ptype_clim_dict (dict): the dictionary to denote the required climate variables for each proxy type; ['tas'] for all by default. calib_period (tuple or list): the time period for calibration. use_predefined_R (bool): use the predefined observation error covariance instead of by calibration. + R_floor (float): the minimum observation error variance assigned from calibration. Defaults to 0 (no floor). + R_floor_frac (float): the minimum observation error variance assigned from calibration, as a fraction of + each record's own variance (e.g., 0.01 caps the signal-to-noise ratio at 10 whatever the record's units). + Where both are set, the larger floor applies. Defaults to 0 (no floor). verbose (bool, optional): print verbose information. Defaults to False. ''' ptype_psm_dict_default = {ptype: 'Linear' for ptype in set(self.proxydb.type_list)} @@ -479,6 +483,9 @@ def calib_psms(self, ptype_psm_dict=None, ptype_season_dict=None, ptype_clim_dic default=[1850, 2015], verbose=verbose) + R_floor = self.io_cfg('psm_R_floor', R_floor, default=0, verbose=verbose) + R_floor_frac = self.io_cfg('psm_R_floor_frac', R_floor_frac, default=0, verbose=verbose) + for pid, pobj in tqdm(self.proxydb.records.items(), total=self.proxydb.nrec, desc='Calibrating the PSMs'): psm_name = ptype_psm_dict[pobj.ptype] @@ -506,16 +513,24 @@ def calib_psms(self, ptype_psm_dict=None, ptype_season_dict=None, ptype_clim_dic pobj.psm.calibrate(season_list=ptype_season_dict[pobj.ptype], calib_period=calib_period, **kwargs) # give the calibrated records a tag + n_floored = 0 for pid, pobj in self.proxydb.records.items(): if pobj.psm.calib_details is None: if verbose: p_warning(f'>>> PSM for {pid} failed to be calibrated.') else: self.proxydb.records[pid].tags.add('calibrated') if not use_predefined_R: - self.proxydb.records[pid].R = pobj.psm.calib_details['PSMmse'] # assign obs err variance + R = pobj.psm.calib_details['PSMmse'] # assign obs err variance + # a near-zero R (an almost perfect calibration fit) can make the Kalman gain blow up + floor = max(R_floor, R_floor_frac * np.nanvar(np.asarray(pobj.value, dtype=float))) + if R < floor: + R = floor + n_floored += 1 + self.proxydb.records[pid].R = R if verbose: p_success(f'>>> {self.proxydb.nrec_tags("calibrated")} records tagged "calibrated" with ProxyRecord.psm created') + if n_floored > 0: p_success(f'>>> R raised to the floor for {n_floored} records') def forward_psms(self, verbose=False, ptype_forward_dict=None): ''' Forward the PSMs. diff --git a/tests/test_reconjob.py b/tests/test_reconjob.py index e32862c9..1d350a3f 100644 --- a/tests/test_reconjob.py +++ b/tests/test_reconjob.py @@ -122,6 +122,32 @@ def test_calibrated_records_have_R(self, calibrated_job): assert hasattr(pobj, 'R'), f'{pid} missing R' assert np.isfinite(pobj.R), f'{pid}.R is not finite' + def _calib(self, **kwargs): + job = ReconJob() + job.prior = {'tas': make_prior()} + job.proxydb = make_proxy_db() + job.calib_psms( + ptype_psm_dict={'tree.TRW': 'TempPlusNoise'}, + ptype_clim_dict={'tree.TRW': ['tas']}, + **kwargs, + ) + return job + + def test_R_floor(self): + job = self._calib(R_floor=5.0) + for pid, pobj in job.proxydb.records.items(): + assert pobj.R == pytest.approx(5.0) + + def test_R_floor_frac(self): + job = self._calib(R_floor_frac=2.0) + for pid, pobj in job.proxydb.records.items(): + assert pobj.R == pytest.approx(2.0 * np.nanvar(pobj.value)) + + def test_no_R_floor_by_default(self): + job = self._calib() + for pid, pobj in job.proxydb.records.items(): + assert pobj.R == pytest.approx(pobj.psm.calib_details['PSMmse']) + def test_calibrated_records_have_psm(self, calibrated_job): for pid, pobj in calibrated_job.proxydb.records.items(): if 'calibrated' in pobj.tags: