From e36bc42dec70188e7acf184d844de8d6240bd0ea Mon Sep 17 00:00:00 2001 From: pciturri Date: Tue, 6 Jun 2023 18:32:40 +0200 Subject: [PATCH 1/3] ft: added preliminary implementation of the brier score test, and an unit_test that only evaluates a forecast and plot the result. --- csep/core/brier_evaluations.py | 156 +++++++++++++++++++++++++++++++++ tests/test_brier.py | 16 ++++ 2 files changed, 172 insertions(+) create mode 100644 csep/core/brier_evaluations.py create mode 100644 tests/test_brier.py diff --git a/csep/core/brier_evaluations.py b/csep/core/brier_evaluations.py new file mode 100644 index 00000000..38483276 --- /dev/null +++ b/csep/core/brier_evaluations.py @@ -0,0 +1,156 @@ +import numpy +import scipy.stats +import scipy.spatial + +from csep.models import EvaluationResult +from csep.core.exceptions import CSEPCatalogException + + +def bier_score_ndarray(forecast, observations): + """ Computes the brier score + """ + + observations = (observations >= 1).astype(int) + prob_success = 1 - scipy.stats.poisson.cdf(0, forecast) + brier = [] + + for p, o in zip(prob_success.ravel(), observations.ravel()): + + if o == 0: + brier.append(-2 * p ** 2) + else: + brier.append(-2 * (p - 1) ** 2) + brier = numpy.sum(brier) + + for n_dim in observations.shape: + brier /= n_dim + + return brier + + +def _simulate_catalog(sim_cells, sampling_weights, sim_fore, random_numbers=None): + # Modified this code to generate simulations in a way that every cell gets one earthquake + # Generate uniformly distributed random numbers in [0,1), this + if random_numbers is None: + # Reset simulation array to zero, but don't reallocate + sim_fore.fill(0) + num_active_cells = 0 + while num_active_cells < sim_cells: + random_num = numpy.random.uniform(0,1) + loc = numpy.searchsorted(sampling_weights, random_num, side='right') + if sim_fore[loc] == 0: + sim_fore[loc] = 1 + num_active_cells = num_active_cells + 1 + else: + # Find insertion points using binary search inserting to satisfy a[i-1] <= v < a[i] + pnts = numpy.searchsorted(sampling_weights, random_numbers, side='right') + # Create simulated catalog by adding to the original locations + numpy.add.at(sim_fore, pnts, 1) + + assert sim_fore.sum() == sim_cells, "simulated the wrong number of events!" + return sim_fore + + +def _brier_score_test(forecast_data, observed_data, num_simulations=1000, random_numbers=None, + seed=None, use_observed_counts=True, verbose=True): + """ Computes binary conditional-likelihood test from CSEP using an efficient simulation based approach. + + Args: + + """ + + # Array-masking that avoids log singularities: + forecast_data = numpy.ma.masked_where(forecast_data <= 0.0, forecast_data) + + # set seed for the likelihood test + if seed is not None: + numpy.random.seed(seed) + + # used to determine where simulated earthquake should be placed, by definition of cumsum these are sorted + sampling_weights = numpy.cumsum(forecast_data.ravel()) / numpy.sum(forecast_data) + + # data structures to store results + sim_fore = numpy.zeros(sampling_weights.shape) + simulated_brier = [] + n_active_cells = len(numpy.unique(numpy.nonzero(observed_data.ravel()))) + + # main simulation step in this loop + for idx in range(num_simulations): + if use_observed_counts: + num_cells_to_simulate = int(n_active_cells) + + if random_numbers is None: + sim_fore = _simulate_catalog(num_cells_to_simulate, + sampling_weights, + sim_fore) + else: + sim_fore = _simulate_catalog(num_cells_to_simulate, + sampling_weights, + sim_fore, + random_numbers=random_numbers[idx, :]) + + # compute Brier score + current_brier = bier_score_ndarray(forecast_data.data, sim_fore) + + # append to list of simulated Brier score + simulated_brier.append(current_brier) + + # just be verbose + if verbose: + if (idx + 1) % 100 == 0: + print(f'... {idx + 1} catalogs simulated.') + + # observed Brier score + obs_brier = bier_score_ndarray(forecast_data.data, observed_data) + + # quantile score + qs = numpy.sum(simulated_brier <= obs_brier) / num_simulations + + # float, float, list + return qs, obs_brier, simulated_brier + + +def brier_score_test(gridded_forecast, + observed_catalog, + num_simulations=1000, + seed=None, + random_numbers=None, + verbose=False): + """ Performs the Brier conditional test on Gridded Forecast using an Observed Catalog. + + Normalizes the forecast so the forecasted rate are consistent with the observations. This modification + eliminates the strong impact differences in the number distribution have on the forecasted rates. + + """ + + # grid catalog onto spatial grid + try: + _ = observed_catalog.region.magnitudes + except CSEPCatalogException: + observed_catalog.region = gridded_forecast.region + gridded_catalog_data = observed_catalog.spatial_magnitude_counts() + + # simply call likelihood test on catalog data and forecast + qs, obs_brier, simulated_brier = _brier_score_test( + gridded_forecast.data, + gridded_catalog_data, + num_simulations=num_simulations, + seed=seed, + random_numbers=random_numbers, + use_observed_counts=True, + verbose=verbose + ) + + # populate result data structure + result = EvaluationResult() + result.test_distribution = simulated_brier + result.name = 'Brier score-Test' + result.observed_statistic = obs_brier + result.quantile = qs + result.sim_name = gridded_forecast.name + result.obs_name = observed_catalog.name + result.status = 'normal' + result.min_mw = numpy.min(gridded_forecast.magnitudes) + + return result + diff --git a/tests/test_brier.py b/tests/test_brier.py new file mode 100644 index 00000000..14ed114a --- /dev/null +++ b/tests/test_brier.py @@ -0,0 +1,16 @@ +import csep +import numpy +from csep.core.brier_evaluations import brier_score_test +from csep.utils.plots import plot_consistency_test +from datetime import datetime + +forecast = csep.load_gridded_forecast(csep.datasets.hires_ssm_italy_fname) + +start = datetime(2010, 1, 1) +end = datetime(2015, 1, 1) +cat = csep.query_bsi(start, end, min_magnitude=5.0) +cat.filter_spatial(region=csep.regions.italy_csep_region( + magnitudes=numpy.arange(5.0, 9.0, 0.1)), in_place=True) + +result = brier_score_test(forecast, cat) +plot_consistency_test(result, show=True) \ No newline at end of file From 5fb55cabe541c5f4458bf44a1bcd36939877658e Mon Sep 17 00:00:00 2001 From: pciturri Date: Mon, 15 Jan 2024 23:00:38 +0100 Subject: [PATCH 2/3] ft: significantly increased brier test speed --- csep/core/brier_evaluations.py | 80 +++++++++++++++------------------- 1 file changed, 36 insertions(+), 44 deletions(-) diff --git a/csep/core/brier_evaluations.py b/csep/core/brier_evaluations.py index 38483276..ae86f0bb 100644 --- a/csep/core/brier_evaluations.py +++ b/csep/core/brier_evaluations.py @@ -1,49 +1,42 @@ import numpy -import scipy.stats -import scipy.spatial +import numpy as np +from scipy.stats import poisson from csep.models import EvaluationResult from csep.core.exceptions import CSEPCatalogException -def bier_score_ndarray(forecast, observations): +def _brier_score_ndarray(forecast, observations): """ Computes the brier score """ - - observations = (observations >= 1).astype(int) - prob_success = 1 - scipy.stats.poisson.cdf(0, forecast) - brier = [] - - for p, o in zip(prob_success.ravel(), observations.ravel()): - - if o == 0: - brier.append(-2 * p ** 2) - else: - brier.append(-2 * (p - 1) ** 2) - brier = numpy.sum(brier) + prob_success = 1 - poisson.cdf(0, forecast) + brier_cell = np.square(prob_success.ravel() - (observations.ravel() > 0)) + brier = -2 * brier_cell.sum() for n_dim in observations.shape: brier /= n_dim - return brier -def _simulate_catalog(sim_cells, sampling_weights, sim_fore, random_numbers=None): - # Modified this code to generate simulations in a way that every cell gets one earthquake - # Generate uniformly distributed random numbers in [0,1), this +def _simulate_catalog(sim_cells, sampling_weights, random_numbers=None): + sim_fore = numpy.zeros(sampling_weights.shape) + if random_numbers is None: # Reset simulation array to zero, but don't reallocate sim_fore.fill(0) num_active_cells = 0 while num_active_cells < sim_cells: random_num = numpy.random.uniform(0,1) - loc = numpy.searchsorted(sampling_weights, random_num, side='right') + loc = numpy.searchsorted(sampling_weights, random_num, + side='right') if sim_fore[loc] == 0: sim_fore[loc] = 1 - num_active_cells = num_active_cells + 1 + num_active_cells += 1 else: - # Find insertion points using binary search inserting to satisfy a[i-1] <= v < a[i] - pnts = numpy.searchsorted(sampling_weights, random_numbers, side='right') + # Find insertion points using binary search inserting + # to satisfy a[i-1] <= v < a[i] + pnts = numpy.searchsorted(sampling_weights, random_numbers, + side='right') # Create simulated catalog by adding to the original locations numpy.add.at(sim_fore, pnts, 1) @@ -51,46 +44,45 @@ def _simulate_catalog(sim_cells, sampling_weights, sim_fore, random_numbers=None return sim_fore -def _brier_score_test(forecast_data, observed_data, num_simulations=1000, random_numbers=None, - seed=None, use_observed_counts=True, verbose=True): - """ Computes binary conditional-likelihood test from CSEP using an efficient simulation based approach. +def _brier_score_test(forecast_data, observed_data, num_simulations=1000, + random_numbers=None, seed=None, verbose=True): + """ Computes binary conditional-likelihood test from CSEP using an + efficient simulation based approach. Args: """ - # Array-masking that avoids log singularities: forecast_data = numpy.ma.masked_where(forecast_data <= 0.0, forecast_data) # set seed for the likelihood test if seed is not None: numpy.random.seed(seed) + import time + start = time.process_time() - # used to determine where simulated earthquake should be placed, by definition of cumsum these are sorted - sampling_weights = numpy.cumsum(forecast_data.ravel()) / numpy.sum(forecast_data) + # used to determine where simulated earthquake should + # be placed, by definition of cumsum these are sorted + sampling_weights = (numpy.cumsum(forecast_data.ravel()) / + numpy.sum(forecast_data)) # data structures to store results - sim_fore = numpy.zeros(sampling_weights.shape) simulated_brier = [] n_active_cells = len(numpy.unique(numpy.nonzero(observed_data.ravel()))) + num_cells_to_simulate = int(n_active_cells) # main simulation step in this loop for idx in range(num_simulations): - if use_observed_counts: - num_cells_to_simulate = int(n_active_cells) - if random_numbers is None: sim_fore = _simulate_catalog(num_cells_to_simulate, - sampling_weights, - sim_fore) + sampling_weights) else: sim_fore = _simulate_catalog(num_cells_to_simulate, sampling_weights, - sim_fore, random_numbers=random_numbers[idx, :]) # compute Brier score - current_brier = bier_score_ndarray(forecast_data.data, sim_fore) + current_brier = _brier_score_ndarray(forecast_data.data, sim_fore) # append to list of simulated Brier score simulated_brier.append(current_brier) @@ -100,8 +92,7 @@ def _brier_score_test(forecast_data, observed_data, num_simulations=1000, random if (idx + 1) % 100 == 0: print(f'... {idx + 1} catalogs simulated.') - # observed Brier score - obs_brier = bier_score_ndarray(forecast_data.data, observed_data) + obs_brier = _brier_score_ndarray(forecast_data.data, observed_data) # quantile score qs = numpy.sum(simulated_brier <= obs_brier) / num_simulations @@ -116,10 +107,12 @@ def brier_score_test(gridded_forecast, seed=None, random_numbers=None, verbose=False): - """ Performs the Brier conditional test on Gridded Forecast using an Observed Catalog. - - Normalizes the forecast so the forecasted rate are consistent with the observations. This modification - eliminates the strong impact differences in the number distribution have on the forecasted rates. + """ Performs the Brier conditional test on Gridded Forecast using an + Observed Catalog. + Normalizes the forecast so the forecasted rate are consistent with the + observations. This modification + eliminates the strong impact differences in the number distribution + have on the forecasted rates. """ @@ -137,7 +130,6 @@ def brier_score_test(gridded_forecast, num_simulations=num_simulations, seed=seed, random_numbers=random_numbers, - use_observed_counts=True, verbose=verbose ) From a50d1e3fb2a43778c29c30e41949a171ab2d7dcc Mon Sep 17 00:00:00 2001 From: pciturri Date: Mon, 22 Jan 2024 20:21:18 +0100 Subject: [PATCH 3/3] ft: finalized brier score implementation as well as brier consistency test. It is implemented only using a Binary calculation for spatial-magnitude bins. tests: Implemented brier score tests. --- csep/core/brier_evaluations.py | 56 +++++++++++++++++++++++-------- tests/test_brier.py | 16 --------- tests/test_evaluations.py | 61 +++++++++++++++++++++++++++++++++- 3 files changed, 103 insertions(+), 30 deletions(-) delete mode 100644 tests/test_brier.py diff --git a/csep/core/brier_evaluations.py b/csep/core/brier_evaluations.py index ae86f0bb..cbe30cf1 100644 --- a/csep/core/brier_evaluations.py +++ b/csep/core/brier_evaluations.py @@ -7,8 +7,22 @@ def _brier_score_ndarray(forecast, observations): - """ Computes the brier score + """ Computes the brier (binary) score for spatial-magnitude cells + using the formula: + + Q(Lambda, Sigma) = 1/N sum_{i=1}^N (Lambda_i - Ind(Sigma_i > 0 ))^2 + + where Lambda is the forecast array, Sigma is the observed catalog, N the + number of spatial-magnitude cells and Ind is the indicator function, which + is 1 if Sigma_i > 0 and 0 otherwise. + + Args: + forecast: 2d array of forecasted rates + observations: 2d array of observed counts + Returns + brier: float, brier score """ + prob_success = 1 - poisson.cdf(0, forecast) brier_cell = np.square(prob_success.ravel() - (observations.ravel() > 0)) brier = -2 * brier_cell.sum() @@ -19,6 +33,18 @@ def _brier_score_ndarray(forecast, observations): def _simulate_catalog(sim_cells, sampling_weights, random_numbers=None): + """ + Simulates a catalog by sampling from the sampling_weights array. + Identical to binomial_evaluations._simulate_catalog + + Args: + sim_cells: + sampling_weights: + random_numbers: + + Returns: + + """ sim_fore = numpy.zeros(sampling_weights.shape) if random_numbers is None: @@ -46,10 +72,19 @@ def _simulate_catalog(sim_cells, sampling_weights, random_numbers=None): def _brier_score_test(forecast_data, observed_data, num_simulations=1000, random_numbers=None, seed=None, verbose=True): - """ Computes binary conditional-likelihood test from CSEP using an - efficient simulation based approach. + """ Computes the Brier consistency test conditional on the total observed + number of events Args: + forecast_data: 2d array of forecasted rates for spatial_magnitude cells + observed_data: 2d array of a catalog resampled to spatial_magnitude + cells + num_simulations: number of synthetic catalog simulations + random_numbers: numpy array of random numbers to use for simulation + seed: seed for random number generator + verbose: print status updates + + """ # Array-masking that avoids log singularities: @@ -58,8 +93,6 @@ def _brier_score_test(forecast_data, observed_data, num_simulations=1000, # set seed for the likelihood test if seed is not None: numpy.random.seed(seed) - import time - start = time.process_time() # used to determine where simulated earthquake should # be placed, by definition of cumsum these are sorted @@ -93,7 +126,6 @@ def _brier_score_test(forecast_data, observed_data, num_simulations=1000, print(f'... {idx + 1} catalogs simulated.') obs_brier = _brier_score_ndarray(forecast_data.data, observed_data) - # quantile score qs = numpy.sum(simulated_brier <= obs_brier) / num_simulations @@ -107,13 +139,11 @@ def brier_score_test(gridded_forecast, seed=None, random_numbers=None, verbose=False): - """ Performs the Brier conditional test on Gridded Forecast using an - Observed Catalog. - Normalizes the forecast so the forecasted rate are consistent with the - observations. This modification - eliminates the strong impact differences in the number distribution - have on the forecasted rates. - + """ + Performs the Brier conditional test on a Gridded Forecast using an + Observed Catalog. Normalizes the forecast so the forecasted rate are + consistent with the observations. This modification eliminates the strong + impact differences in the number distribution have on the forecasted rates. """ # grid catalog onto spatial grid diff --git a/tests/test_brier.py b/tests/test_brier.py deleted file mode 100644 index 14ed114a..00000000 --- a/tests/test_brier.py +++ /dev/null @@ -1,16 +0,0 @@ -import csep -import numpy -from csep.core.brier_evaluations import brier_score_test -from csep.utils.plots import plot_consistency_test -from datetime import datetime - -forecast = csep.load_gridded_forecast(csep.datasets.hires_ssm_italy_fname) - -start = datetime(2010, 1, 1) -end = datetime(2015, 1, 1) -cat = csep.query_bsi(start, end, min_magnitude=5.0) -cat.filter_spatial(region=csep.regions.italy_csep_region( - magnitudes=numpy.arange(5.0, 9.0, 0.1)), in_place=True) - -result = brier_score_test(forecast, cat) -plot_consistency_test(result, show=True) \ No newline at end of file diff --git a/tests/test_evaluations.py b/tests/test_evaluations.py index c7468528..6f8791c7 100644 --- a/tests/test_evaluations.py +++ b/tests/test_evaluations.py @@ -1,10 +1,12 @@ import os import numpy +import numpy as np +import scipy import unittest import csep.core.poisson_evaluations as poisson import csep.core.binomial_evaluations as binary - +import csep.core.brier_evaluations as brier def get_datadir(): root_dir = os.path.dirname(os.path.abspath(__file__)) data_dir = os.path.join(root_dir, 'artifacts', 'Comcat') @@ -115,5 +117,62 @@ def test_binomial_likelihood(self): numpy.testing.assert_allclose(simulated_ll[0], -7.921741654647629) +class TestPoissonBrier(unittest.TestCase): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.seed = 1 + self.forecast_data = numpy.array([[0.3, 0.2, 0.1], [0.2, 0.1, 0.1]]) + self.observed_data = numpy.array([[0, 1, 2], [1, 1, 0]]) + + def test_brier_score_calculation(self): + + # 1 bin + rate = 1 + prob = 1 - scipy.stats.poisson.pmf(0, rate) + brier_score_hit = -2 * (prob - 1)**2 + brier_score = brier._brier_score_ndarray(numpy.array([[rate]]), + numpy.array([[1]])) + numpy.testing.assert_allclose(brier_score_hit, brier_score) + + brier_score_nohit = -2 * prob**2 + brier_score = brier._brier_score_ndarray(numpy.array([[rate]]), + numpy.array([[0]])) + numpy.testing.assert_allclose(brier_score_nohit, brier_score) + # 2 bins + rate = np.array([1, 0.5]) + hits = np.array([0, 1]) + prob = 1 - scipy.stats.poisson.pmf(0, rate) + brier_score_2bins = (-2 * (prob - hits)**2).sum() / len(rate) + brier_score = brier._brier_score_ndarray(np.array(rate), + np.array([hits])) + numpy.testing.assert_allclose(brier_score_2bins, brier_score) + + hits = np.array([0, 0]) + brier_score_2bins = (-2 * (prob - hits)**2).sum() / len(rate) + brier_score = brier._brier_score_ndarray(np.array(rate), + np.array([hits])) + numpy.testing.assert_allclose(brier_score_2bins, brier_score) + + def test_brier_test(self): + + expected_sim = numpy.array([1, 1, 1, 1, 0, 0]) + qs, obs_brier, simulated_brier = brier._brier_score_test( + self.forecast_data, + self.observed_data, + num_simulations=1, + seed=self.seed, + verbose=True) + + probs = 1 - scipy.stats.poisson.pmf(0, self.forecast_data.ravel()) + sim_brier_analytic = ((-2 * (probs - expected_sim)**2).sum() / + len(probs)) + numpy.testing.assert_allclose(simulated_brier[0], sim_brier_analytic) + obs_brier_analytic = ( + (-2 * (probs - self.observed_data.ravel().astype(bool))**2).sum() / + len(probs)) + numpy.testing.assert_allclose(obs_brier, obs_brier_analytic) + + if __name__ == '__main__': unittest.main()