Source code for csep.core.catalog_evaluations

# Third-Party Imports
import time

import numpy
import scipy.stats

# PyCSEP imports
from csep.core.exceptions import CSEPEvaluationException
from csep.models import (
    CatalogNumberTestResult,
    CatalogSpatialTestResult,
    CatalogMagnitudeTestResult,
    CatalogPseudolikelihoodTestResult,
    CalibrationTestResult
)
from csep.utils.calc import _compute_likelihood
from csep.utils.stats import get_quantiles, cumulative_square_diff


[docs]def number_test(forecast, observed_catalog, verbose=True): """ Performs the number test on a catalog-based forecast. The number test builds an empirical distribution of the event counts for each data. By default, this function does not perform any filtering on the catalogs in the forecast or observation. These should be handled outside of the function. Args: forecast (:class:`csep.core.forecasts.CatalogForecast`): forecast to evaluate observed_catalog (:class:`csep.core.catalogs.AbstractBaseCatalog`): evaluation data Returns: evaluation result (:class:`csep.models.EvaluationResult`): evaluation result """ event_counts = [] t0 = time.time() for i, catalog in enumerate(forecast): # output status if verbose: tens_exp = numpy.floor(numpy.log10(i + 1)) if (i + 1) % 10 ** tens_exp == 0: t1 = time.time() print(f'Processed {i + 1} catalogs in {t1 - t0} seconds', flush=True) event_counts.append(catalog.event_count) obs_count = observed_catalog.event_count delta_1, delta_2 = get_quantiles(event_counts, obs_count) # prepare result result = CatalogNumberTestResult(test_distribution=event_counts, name='Catalog N-Test', observed_statistic=obs_count, quantile=(delta_1, delta_2), status='normal', obs_catalog_repr=str(observed_catalog), sim_name=forecast.name, min_mw=forecast.min_magnitude, obs_name=observed_catalog.name) return result
[docs]def spatial_test(forecast, observed_catalog, verbose=True): """ Performs spatial test for catalog-based forecasts. Args: forecast: CatalogForecast observed_catalog: CSEPCatalog filtered to be consistent with the forecast Returns: CatalogSpatialTestResult """ if forecast.region is None: raise CSEPEvaluationException("Forecast must have region member to perform spatial test.") # get observed likelihood if observed_catalog.event_count == 0: print(f'Spatial test not-invalid because no events in observed catalog.') test_distribution = [] # compute expected rates for forecast if needed if forecast.expected_rates is None: forecast.get_expected_rates(verbose=verbose) expected_cond_count = forecast.expected_rates.sum() forecast_mean_spatial_rates = forecast.expected_rates.spatial_counts() # summing over spatial counts ensures that the correct number of events are used; even through the catalogs should # be filtered before calling this function gridded_obs = observed_catalog.spatial_counts() n_obs = numpy.sum(gridded_obs) # iterate through catalogs in forecast and compute likelihood t0 = time.time() for i, catalog in enumerate(forecast): gridded_cat = catalog.spatial_counts() _, lh_norm = _compute_likelihood(gridded_cat, forecast_mean_spatial_rates, expected_cond_count, n_obs) test_distribution.append(lh_norm) # output status if verbose: tens_exp = numpy.floor(numpy.log10(i + 1)) if (i + 1) % 10 ** tens_exp == 0: t1 = time.time() print(f'Processed {i + 1} catalogs in {t1 - t0} seconds', flush=True) _, obs_lh_norm = _compute_likelihood(gridded_obs, forecast_mean_spatial_rates, expected_cond_count, n_obs) # if obs_lh is -numpy.inf, recompute but only for indexes where obs and simulated are non-zero message = "normal" if obs_lh_norm == -numpy.inf: idx_good_sim = forecast_mean_spatial_rates != 0 new_gridded_obs = gridded_obs[idx_good_sim] new_n_obs = numpy.sum(new_gridded_obs) print(f"Found -inf as the observed likelihood score. " f"Assuming event(s) occurred in undersampled region of forecast.\n" f"Recomputing with {new_n_obs} events after removing {n_obs - new_n_obs} events.") new_ard = forecast_mean_spatial_rates[idx_good_sim] _, obs_lh_norm = _compute_likelihood(new_gridded_obs, new_ard, expected_cond_count, n_obs) message = "undersampled" # check for nans here and remove from spatial distribution test_distribution_spatial_1d = numpy.array(test_distribution) if numpy.isnan(numpy.sum(test_distribution_spatial_1d)): test_distribution_spatial_1d = test_distribution_spatial_1d[~numpy.isnan(test_distribution_spatial_1d)] if n_obs == 0 or numpy.isnan(obs_lh_norm): message = "not-valid" delta_1, delta_2 = -1, -1 else: delta_1, delta_2 = get_quantiles(test_distribution_spatial_1d, obs_lh_norm) result = CatalogSpatialTestResult(test_distribution=test_distribution_spatial_1d, name='S-Test', observed_statistic=obs_lh_norm, quantile=(delta_1, delta_2), status=message, min_mw=forecast.min_magnitude, obs_catalog_repr=str(observed_catalog), sim_name=forecast.name, obs_name=observed_catalog.name) return result
[docs]def magnitude_test(forecast, observed_catalog, verbose=True): """ Performs magnitude test for catalog-based forecasts """ test_distribution = [] if forecast.region.magnitudes is None: raise CSEPEvaluationException("Forecast must have region.magnitudes member to perform magnitude test.") # short-circuit if zero events if observed_catalog.event_count == 0: print("Cannot perform magnitude test when observed event count is zero.") # prepare result result = CatalogMagnitudeTestResult(test_distribution=test_distribution, name='M-Test', observed_statistic=None, quantile=(None, None), status='not-valid', min_mw=forecast.min_magnitude, obs_catalog_repr=str(observed_catalog), obs_name=observed_catalog.name, sim_name=forecast.name) return result # compute expected rates for forecast if needed if forecast.expected_rates is None: forecast.get_expected_rates(verbose=verbose) # returns the average events in the magnitude bins union_histogram = forecast.expected_rates.magnitude_counts() n_union_events = numpy.sum(union_histogram) obs_histogram = observed_catalog.magnitude_counts() n_obs = numpy.sum(obs_histogram) union_scale = n_obs / n_union_events scaled_union_histogram = union_histogram * union_scale # compute the test statistic for each catalog t0 = time.time() for i, catalog in enumerate(forecast): mag_counts = catalog.magnitude_counts() n_events = numpy.sum(mag_counts) if n_events == 0: # print("Skipping to next because catalog contained zero events.") continue scale = n_obs / n_events catalog_histogram = mag_counts * scale # compute magnitude test statistic for the catalog test_distribution.append( cumulative_square_diff(numpy.log10(catalog_histogram + 1), numpy.log10(scaled_union_histogram + 1)) ) # output status if verbose: tens_exp = numpy.floor(numpy.log10(i + 1)) if (i + 1) % 10 ** tens_exp == 0: t1 = time.time() print(f'Processed {i + 1} catalogs in {t1 - t0} seconds', flush=True) # compute observed statistic obs_d_statistic = cumulative_square_diff(numpy.log10(obs_histogram + 1), numpy.log10(scaled_union_histogram + 1)) # score evaluation delta_1, delta_2 = get_quantiles(test_distribution, obs_d_statistic) # prepare result result = CatalogMagnitudeTestResult(test_distribution=test_distribution, name='M-Test', observed_statistic=obs_d_statistic, quantile=(delta_1, delta_2), status='normal', min_mw=forecast.min_magnitude, obs_catalog_repr=str(observed_catalog), obs_name=observed_catalog.name, sim_name=forecast.name) return result
[docs]def pseudolikelihood_test(forecast, observed_catalog, verbose=True): """ Performs the spatial pseudolikelihood test for catalog forecasts. Performs the spatial pseudolikelihood test as described by Savran et al., 2020. The tests uses a pseudolikelihood statistic computed from the expected rates in spatial cells. A pseudolikelihood test based on space-magnitude bins is in a development mode and does not exist currently. Args: forecast: :class:`csep.core.forecasts.CatalogForecast` observed_catalog: :class:`csep.core.catalogs.AbstractBaseCatalog` """ if forecast.region is None: raise CSEPEvaluationException("Forecast must have region member to perform spatial test.") # get observed likelihood if observed_catalog.event_count == 0: print(f'Skipping pseudolikelihood test because no events in observed catalog.') return None test_distribution = [] # compute expected rates for forecast if needed if forecast.expected_rates is None: _ = forecast.get_expected_rates(verbose=verbose) expected_cond_count = forecast.expected_rates.sum() forecast_mean_spatial_rates = forecast.expected_rates.spatial_counts() # summing over spatial counts ensures that the correct number of events are used; even through the catalogs should # be filtered before calling this function gridded_obs = observed_catalog.spatial_counts() n_obs = numpy.sum(gridded_obs) t0 = time.time() for i, catalog in enumerate(forecast): gridded_cat = catalog.spatial_counts() plh, _ = _compute_likelihood(gridded_cat, forecast_mean_spatial_rates, expected_cond_count, n_obs) test_distribution.append(plh) # output status if verbose: tens_exp = numpy.floor(numpy.log10(i + 1)) if (i + 1) % 10 ** tens_exp == 0: t1 = time.time() print(f'Processed {i + 1} catalogs in {t1 - t0} seconds', flush=True) obs_plh, _ = _compute_likelihood(gridded_obs, forecast_mean_spatial_rates, expected_cond_count, n_obs) # if obs_lh is -numpy.inf, recompute but only for indexes where obs and simulated are non-zero message = "normal" if obs_plh == -numpy.inf: idx_good_sim = forecast_mean_spatial_rates != 0 new_gridded_obs = gridded_obs[idx_good_sim] new_n_obs = numpy.sum(new_gridded_obs) print(f"Found -inf as the observed likelihood score. " f"Assuming event(s) occurred in undersampled region of forecast.\n" f"Recomputing with {new_n_obs} events after removing {n_obs - new_n_obs} events.") if new_n_obs == 0: print( f'Skipping pseudo-likelihood based tests for because no events in observed catalog ' f'after correcting for under-sampling in forecast.' ) return None new_ard = forecast_mean_spatial_rates[idx_good_sim] # we need to use the old n_obs here, because if we normalize the ard to a different value the observed # statistic will not be computed correctly. obs_plh, _ = _compute_likelihood(new_gridded_obs, new_ard, expected_cond_count, n_obs) message = "undersampled" # check for nans here test_distribution_1d = numpy.array(test_distribution) if numpy.isnan(numpy.sum(test_distribution_1d)): test_distribution_1d = test_distribution_1d[~numpy.isnan(test_distribution_1d)] if n_obs == 0 or numpy.isnan(obs_plh): message = "not-valid" delta_1, delta_2 = -1, -1 else: delta_1, delta_2 = get_quantiles(test_distribution_1d, obs_plh) # prepare evaluation result result = CatalogPseudolikelihoodTestResult( test_distribution=test_distribution_1d, name='PL-Test', observed_statistic=obs_plh, quantile=(delta_1, delta_2), status=message, min_mw=forecast.min_magnitude, obs_catalog_repr=str(observed_catalog), sim_name=forecast.name, obs_name=observed_catalog.name ) return result
[docs]def calibration_test(evaluation_results, delta_1=False): """ Perform the calibration test by computing a Kilmogorov-Smirnov test of the observed quantiles against a uniform distribution. Args: evaluation_results: iterable of evaluation result objects delta_1 (bool): use delta_1 for quantiles. default false -> use delta_2 quantile score for calibration test """ # this is using "delta_2" which is the cdf value less-equal idx = 0 if delta_1 else 1 quantiles = [] for result in evaluation_results: if result.status == 'not-valid': print(f'evaluation not valid for {result.name}. skipping in calibration test.') else: quantiles.append(result.quantile[idx]) ks, p_value = scipy.stats.kstest(quantiles, 'uniform') result = CalibrationTestResult( test_distribution = quantiles, name=f'{evaluation_results[0].name} Calibration Test', observed_statistic=ks, quantile=p_value, status='normal', min_mw = evaluation_results[0].min_mw, obs_catalog_repr=evaluation_results[0].obs_catalog_repr, sim_name=evaluation_results[0].sim_name, obs_name=evaluation_results[0].obs_name ) return result