diff --git a/CHANGELOG.md b/CHANGELOG.md index f3057973..f82f1d50 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,29 @@ +# v0.6.1 (12/12/2022) + +# Change-log +Added quadtree csv reader ([#186](https://github.com/SCECcode/pycsep/pull/186)) +Non-Poissonian tests +([#189](https://github.com/SCECcode/pycsep/pull/189), +[#205](https://github.com/SCECcode/pycsep/pull/205), +[#208](https://github.com/SCECcode/pycsep/pull/208), +[#209](https://github.com/SCECcode/pycsep/pull/209)) +Added plots for p-values, and confidence ranges for consistency tests ([#190](https://github.com/SCECcode/pycsep/pull/190)) +Added NZ testing and collection regions ([#198](https://github.com/SCECcode/pycsep/pull/198)) +Fixed region border plotting issue ([#199](https://github.com/SCECcode/pycsep/pull/199)) +Added documentation for non-Poissonian tests ([#202](https://github.com/SCECcode/pycsep/pull/202)) +Support for BSI catalog ([#201](https://github.com/SCECcode/pycsep/pull/201)) +Fixed compatibility with new version of matplotlib ([#206](https://github.com/SCECcode/pycsep/pull/206)) + +## Credits +Pablo Iturrieta (@pabloitu) +Jose Bayona (@bayonato89) +Khawaja Asim (@khawajasim) +William Savran (@wsavran) + # v0.6.0 (02/04/2022) ## Change-log -Adds support for quadtree regions [#184](https://github.com/SCECcode/pycsep/pull/184) +Adds support for quadtree regions ([#184])(https://github.com/SCECcode/pycsep/pull/184) ## Credits Khawaja Asim (@khawajasim) diff --git a/CITATION.cff b/CITATION.cff index 670b50d9..5b01f9cd 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -14,6 +14,6 @@ authors: given-names: Philip J. orcid: 0000-0002-9221-7068 title: "pyCSEP - Tools for Earthquake Forecast Developers" -version: 0.4.1 +version: 0.6.1 repository: https://github.com/SCECcode/pycsep date-released: 2021-04-20 diff --git a/codemeta.json b/codemeta.json index 31301e2c..c8f61d30 100644 --- a/codemeta.json +++ b/codemeta.json @@ -9,7 +9,7 @@ "downloadUrl": "https://github.com/SCECcode/pycsep", "issueTracker": "https://github.com/SCECcode/pycsep/issues", "name": "pyCSEP", - "version": "v0.6.0", + "version": "v0.6.1", "description": "The pyCSEP Toolkit helps earthquake forecast model developers evaluate their forecasts with the goal of understanding earthquake predictability.", "applicationCategory": "Seismology", "developmentStatus": "active", @@ -30,7 +30,7 @@ "NumPy 1.21.3 or later (https://numpy.org)", "SciPy 1.7.1 or later (https://scipy.org)", "pandas 1.3.4 or later (https://pandas.pydata.org)", - "cartopy 0.20.0 or later (https://scitools.org.uk/cartopy/docs/latest)", + "cartopy 0.21.5 or later (https://scitools.org.uk/cartopy/docs/latest)", "GEOS 3.7.2 or later (https://trac.osgeo.org/geos/)", "PROJ 8.0.0 or later (https://proj.org/)" ], diff --git a/csep/__init__.py b/csep/__init__.py index 4149d43c..de1a03f8 100644 --- a/csep/__init__.py +++ b/csep/__init__.py @@ -188,7 +188,9 @@ def load_catalog(filename, type='csep-csv', format='native', loader=None, apply_ def query_comcat(start_time, end_time, min_magnitude=2.50, min_latitude=31.50, max_latitude=43.00, - min_longitude=-125.40, max_longitude=-113.10, verbose=True, + min_longitude=-125.40, max_longitude=-113.10, + max_depth=1000, + verbose=True, apply_filters=False, **kwargs): """ Access Comcat catalog through web service @@ -201,11 +203,11 @@ def query_comcat(start_time, end_time, min_magnitude=2.50, max_latitude: max latitude of bounding box min_longitude: min latitude of bounding box max_longitude: max longitude of bounding box - region: :class:`csep.core.regions.CartesianGrid2D + max_depth: maximum depth of the bounding box verbose (bool): print catalog summary statistics Returns: - :class:`csep.core.catalogs.ComcatCatalog + :class:`csep.core.catalogs.CSEPCatalog """ # Timezone should be in UTC @@ -213,7 +215,8 @@ def query_comcat(start_time, end_time, min_magnitude=2.50, eventlist = readers._query_comcat(start_time=start_time, end_time=end_time, min_magnitude=min_magnitude, min_latitude=min_latitude, max_latitude=max_latitude, - min_longitude=min_longitude, max_longitude=max_longitude) + min_longitude=min_longitude, max_longitude=max_longitude, + max_depth=max_depth) t1 = time.time() comcat = catalogs.CSEPCatalog(data=eventlist, date_accessed=utc_now_datetime(), **kwargs) print("Fetched ComCat catalog in {} seconds.\n".format(t1 - t0)) @@ -234,6 +237,59 @@ def query_comcat(start_time, end_time, min_magnitude=2.50, return comcat + +def query_bsi(start_time, end_time, min_magnitude=2.50, + min_latitude=32.0, max_latitude=50.0, + min_longitude=2.0, max_longitude=21.0, + max_depth=1000, + verbose=True, + apply_filters=False, **kwargs): + """ + Access BSI catalog through web service + + Args: + start_time: datetime object of start of catalog + end_time: datetime object for end of catalog + min_magnitude: minimum magnitude to query + min_latitude: maximum magnitude to query + max_latitude: max latitude of bounding box + min_longitude: min latitude of bounding box + max_longitude: max longitude of bounding box + max_depth: maximum depth of the bounding box + verbose (bool): print catalog summary statistics + + Returns: + :class:`csep.core.catalogs.CSEPCatalog + """ + + # Timezone should be in UTC + t0 = time.time() + eventlist = readers._query_bsi(start_time=start_time, end_time=end_time, + min_magnitude=min_magnitude, + min_latitude=min_latitude, max_latitude=max_latitude, + min_longitude=min_longitude, max_longitude=max_longitude, + max_depth=max_depth) + t1 = time.time() + bsi = catalogs.CSEPCatalog(data=eventlist, date_accessed=utc_now_datetime(), **kwargs) + print("Fetched BSI catalog in {} seconds.\n".format(t1 - t0)) + + if apply_filters: + try: + bsi = bsi.filter().filter_spatial() + except CSEPCatalogException: + bsi = bsi.filter() + + if verbose: + print("Downloaded catalog from Bollettino Sismico Italiano (BSI) with following parameters") + print("Start Date: {}\nEnd Date: {}".format(str(bsi.start_time), str(bsi.end_time))) + print("Min Latitude: {} and Max Latitude: {}".format(bsi.min_latitude, bsi.max_latitude)) + print("Min Longitude: {} and Max Longitude: {}".format(bsi.min_longitude, bsi.max_longitude)) + print("Min Magnitude: {}".format(bsi.min_magnitude)) + print(f"Found {bsi.event_count} events in the BSI catalog.") + + return bsi + + def load_evaluation_result(fname): """ Load evaluation result stored as json file diff --git a/csep/_version.py b/csep/_version.py index d07e93fd..37416fa9 100644 --- a/csep/_version.py +++ b/csep/_version.py @@ -1,2 +1,2 @@ -__version__ = "0.6.0" +__version__ = "0.6.1" diff --git a/csep/core/binomial_evaluations.py b/csep/core/binomial_evaluations.py new file mode 100644 index 00000000..58fd7f50 --- /dev/null +++ b/csep/core/binomial_evaluations.py @@ -0,0 +1,395 @@ +import numpy +import scipy.stats +import scipy.spatial + +from csep.models import EvaluationResult +from csep.core.exceptions import CSEPCatalogException + + +def _nbd_number_test_ndarray(fore_cnt, obs_cnt, variance, epsilon=1e-6): + """ Computes delta1 and delta2 values from the Negative Binomial (NBD) number test. + + Args: + fore_cnt (float): parameter of negative binomial distribution coming from expected value of the forecast + obs_cnt (float): count of earthquakes observed during the testing period. + variance (float): variance parameter of negative binomial distribution coming from historical catalog. + A variance value of approximately 23541 has been calculated using M5.95+ earthquakes observed worldwide from 1982 to 2013. + epsilon (float): tolerance level to satisfy the requirements of two-sided p-value + + Returns + result (tuple): (delta1, delta2) + """ + var = variance + mean = fore_cnt + upsilon = 1.0 - ((var - mean) / var) + tau = (mean**2 /(var - mean)) + + delta1 = 1.0 - scipy.stats.nbinom.cdf(obs_cnt - epsilon, tau, upsilon, loc=0) + delta2 = scipy.stats.nbinom.cdf(obs_cnt + epsilon, tau, upsilon, loc=0) + + return delta1, delta2 + + +def negative_binomial_number_test(gridded_forecast, observed_catalog, variance): + """ Computes "negative binomial N-Test" on a gridded forecast. + + Computes Number (N) test for Observed and Forecasts. Both data sets are expected to be in terms of event counts. + We find the Total number of events in Observed Catalog and Forecasted Catalogs. Which are then employed to compute the + probablities of + (i) At least no. of events (delta 1) + (ii) At most no. of events (delta 2) assuming the negative binomial distribution. + + Args: + gridded_forecast: Forecast of a Model (Gridded) (Numpy Array) + A forecast has to be in terms of Average Number of Events in Each Bin + It can be anything greater than zero + observed_catalog: Observed (Gridded) seismicity (Numpy Array): + An Observation has to be Number of Events in Each Bin + It has to be a either zero or positive integer only (No Floating Point) + variance: Variance parameter of negative binomial distribution obtained from historical catalog. + + Returns: + out (tuple): (delta_1, delta_2) + """ + result = EvaluationResult() + + # observed count + obs_cnt = observed_catalog.event_count + + # forecasts provide the expeceted number of events during the time horizon of the forecast + fore_cnt = gridded_forecast.event_count + + epsilon = 1e-6 + + # stores the actual result of the number test + delta1, delta2 = _nbd_number_test_ndarray(fore_cnt, obs_cnt, variance, epsilon=epsilon) + + # store results + result.test_distribution = ('negative_binomial', fore_cnt) + result.name = 'NBD N-Test' + result.observed_statistic = obs_cnt + result.quantile = (delta1, delta2) + 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 + + +def binary_joint_log_likelihood_ndarray(forecast, catalog): + """ Computes Bernoulli log-likelihood scores, assuming that earthquakes follow a binomial distribution. + + Args: + forecast: Forecast of a Model (Gridded) (Numpy Array) + A forecast has to be in terms of Average Number of Events in Each Bin + It can be anything greater than zero + catalog: Observed (Gridded) seismicity (Numpy Array): + An Observation has to be Number of Events in Each Bin + It has to be a either zero or positive integer only (No Floating Point) + """ + # First, we mask the forecast in cells where we could find log=0.0 singularities: + forecast_masked = numpy.ma.masked_where(forecast.ravel() <= 0.0, forecast.ravel()) + # Then, we compute the log-likelihood of observing one or more events given a Poisson distribution, i.e., 1 - Pr(0) + target_idx = numpy.nonzero(catalog.ravel()) + y = numpy.zeros(forecast_masked.ravel().shape) + y[target_idx[0]] = 1 + first_term = y * (numpy.log(1.0 - numpy.exp(-forecast_masked.ravel()))) + # Also, we estimate the log-likelihood in cells no events are observed: + second_term = (1-y) * (-forecast_masked.ravel().data) + # Finally, we sum both terms to compute the joint log-likelihood score: + return sum(first_term.data + second_term.data) + + + +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 _binary_likelihood_test(forecast_data, observed_data, num_simulations=1000, random_numbers=None, + seed=None, use_observed_counts=True, verbose=True, normalize_likelihood=False): + """ Computes binary conditional-likelihood test from CSEP using an efficient simulation based approach. + + Args: + forecast_data (numpy.ndarray): nd array where [:, -1] are the magnitude bins. + observed_data (numpy.ndarray): same format as observation. + num_simulations: default number of simulations to use for likelihood based simulations + seed: used for reproducibility of the prng + random_numbers (numpy.ndarray): can supply an explicit list of random numbers, primarily used for software testing + use_observed_counts (bool): if true, will simulate catalogs using the observed events, if false will draw from poisson + distribution + """ + + # 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_ll = [] + n_active_cells = len(numpy.unique(numpy.nonzero(observed_data.ravel()))) + n_fore = numpy.sum(forecast_data) + expected_forecast_count = 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) + else: + sim_fore = _simulate_catalog(num_cells_to_simulate, sampling_weights, sim_fore, + random_numbers=random_numbers[idx,:]) + + # compute joint log-likelihood + current_ll = binary_joint_log_likelihood_ndarray(forecast_data.data, sim_fore) + + # append to list of simulated log-likelihoods + simulated_ll.append(current_ll) + + # just be verbose + if verbose: + if (idx + 1) % 100 == 0: + print(f'... {idx + 1} catalogs simulated.') + + # observed joint log-likelihood + obs_ll = binary_joint_log_likelihood_ndarray(forecast_data.data, observed_data) + + # quantile score + qs = numpy.sum(simulated_ll <= obs_ll) / num_simulations + + # float, float, list + return qs, obs_ll, simulated_ll + + +def binary_spatial_test(gridded_forecast, observed_catalog, num_simulations=1000, seed=None, random_numbers=None, verbose=False): + """ Performs the binary spatial test on the Forecast using the Observed Catalogs. + + Note: The forecast and the observations should be scaled to the same time period before calling this function. This increases + transparency as no assumptions are being made about the length of the forecasts. This is particularly important for + gridded forecasts that supply their forecasts as rates. + + Args: + gridded_forecast: csep.core.forecasts.GriddedForecast + observed_catalog: csep.core.catalogs.Catalog + num_simulations (int): number of simulations used to compute the quantile score + seed (int): used fore reproducibility, and testing + random_numbers (numpy.ndarray): random numbers used to override the random number generation. injection point for testing. + + Returns: + evaluation_result: csep.core.evaluations.EvaluationResult + """ + + # grid catalog onto spatial grid + gridded_catalog_data = observed_catalog.spatial_counts() + + # simply call likelihood test on catalog data and forecast + qs, obs_ll, simulated_ll = _binary_likelihood_test( + gridded_forecast.spatial_counts(), + gridded_catalog_data, + num_simulations=num_simulations, + seed=seed, + random_numbers=random_numbers, + use_observed_counts=True, + verbose=verbose, + normalize_likelihood=True + ) + + +# populate result data structure + result = EvaluationResult() + result.test_distribution = simulated_ll + result.name = 'Binary S-Test' + result.observed_statistic = obs_ll + result.quantile = qs + result.sim_name = gridded_forecast.name + result.obs_name = observed_catalog.name + result.status = 'normal' + try: + result.min_mw = numpy.min(gridded_forecast.magnitudes) + except AttributeError: + result.min_mw = -1 + return result + + +def binary_conditional_likelihood_test(gridded_forecast, observed_catalog, num_simulations=1000, seed=None, random_numbers=None, verbose=False): + """ Performs the binary conditional likelihood 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. + + Note: The forecast and the observations should be scaled to the same time period before calling this function. This increases + transparency as no assumptions are being made about the length of the forecasts. This is particularly important for + gridded forecasts that supply their forecasts as rates. + + Args: + gridded_forecast: csep.core.forecasts.GriddedForecast + observed_catalog: csep.core.catalogs.Catalog + num_simulations (int): number of simulations used to compute the quantile score + seed (int): used fore reproducibility, and testing + random_numbers (numpy.ndarray): random numbers used to override the random number generation. injection point for testing. + + Returns: + evaluation_result: csep.core.evaluations.EvaluationResult + """ + + # 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_ll, simulated_ll = _binary_likelihood_test( + gridded_forecast.data, + gridded_catalog_data, + num_simulations=num_simulations, + seed=seed, + random_numbers=random_numbers, + use_observed_counts=True, + verbose=verbose, + normalize_likelihood=False + ) + + # populate result data structure + result = EvaluationResult() + result.test_distribution = simulated_ll + result.name = 'Binary CL-Test' + result.observed_statistic = obs_ll + 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 + + +def matrix_binary_t_test(target_event_rates1, target_event_rates2, n_obs, n_f1, n_f2, catalog, alpha=0.05): + """ Computes binary T test statistic by comparing two target event rate distributions. + + We compare Forecast from Model 1 and with Forecast of Model 2. Information Gain per Active Bin (IGPA) is computed, which is then + employed to compute T statistic. Confidence interval of Information Gain can be computed using T_critical. For a complete + explanation see Rhoades, D. A., et al., (2011). Efficient testing of earthquake forecasting models. Acta Geophysica, 59(4), + 728-747. doi:10.2478/s11600-011-0013-5, and Bayona J.A. et al., (2022). Prospective evaluation of multiplicative hybrid earthquake + forecasting models in California. doi: 10.1093/gji/ggac018. + + Args: + target_event_rates1 (numpy.ndarray): nd-array storing target event rates + target_event_rates2 (numpy.ndarray): nd-array storing target event rates + n_obs (float, int, numpy.ndarray): number of observed earthquakes, should be whole number and >= zero. + n_f1 (float): Total number of forecasted earthquakes by Model 1 + n_f2 (float): Total number of forecasted earthquakes by Model 2 + catalog: csep.core.catalogs.Catalog + alpha (float): tolerance level for the type-i error rate of the statistical test + + Returns: + out (dict): relevant statistics from the t-test + """ + # Some Pre Calculations - Because they are being used repeatedly. + N_p = n_obs + N = len(np.unique(np.nonzero(catalog.spatial_magnitude_counts().ravel()))) # Number of active bins + N1 = n_f1 + N2 = n_f2 + X1 = numpy.log(target_event_rates1) # Log of every element of Forecast 1 + X2 = numpy.log(target_event_rates2) # Log of every element of Forecast 2 + + + # Information Gain, using Equation (17) of Rhoades et al. 2011 + information_gain = (numpy.sum(X1 - X2) - (N1 - N2)) / N + + # Compute variance of (X1-X2) using Equation (18) of Rhoades et al. 2011 + first_term = (numpy.sum(numpy.power((X1 - X2), 2))) / (N - 1) + second_term = numpy.power(numpy.sum(X1 - X2), 2) / (numpy.power(N, 2) - N) + forecast_variance = first_term - second_term + + forecast_std = numpy.sqrt(forecast_variance) + t_statistic = information_gain / (forecast_std / numpy.sqrt(N)) + + # Obtaining the Critical Value of T from T distribution. + df = N - 1 + t_critical = scipy.stats.t.ppf(1 - (alpha / 2), df) # Assuming 2-Tail Distribution for 2 tail, divide 0.05/2. + + # Computing Information Gain Interval. + ig_lower = information_gain - (t_critical * forecast_std / numpy.sqrt(N)) + ig_upper = information_gain + (t_critical * forecast_std / numpy.sqrt(N)) + + # If T value greater than T critical, Then both Lower and Upper Confidence Interval limits will be greater than Zero. + # If above Happens, Then It means that Forecasting Model 1 is better than Forecasting Model 2. + return {'t_statistic': t_statistic, + 't_critical': t_critical, + 'information_gain': information_gain, + 'ig_lower': ig_lower, + 'ig_upper': ig_upper} + + +def binary_paired_t_test(forecast, benchmark_forecast, observed_catalog, alpha=0.05, scale=False): + """ Computes the binary t-test for gridded earthquake forecasts. + + This score is positively oriented, meaning that positive values of the information gain indicate that the + forecast is performing better than the benchmark forecast + + Args: + forecast (csep.core.forecasts.GriddedForecast): nd-array storing gridded rates, axis=-1 should be the magnitude column + benchmark_forecast (csep.core.forecasts.GriddedForecast): nd-array storing gridded rates, axis=-1 should be the magnitude + column + observed_catalog (csep.core.catalogs.AbstractBaseCatalog): number of observed earthquakes, should be whole number and >= zero. + alpha (float): tolerance level for the type-i error rate of the statistical test + scale (bool): if true, scale forecasted rates down to a single day + + Returns: + evaluation_result: csep.core.evaluations.EvaluationResult + """ + + # needs some pre-processing to put the forecasts in the context that is required for the t-test. this is different + # for cumulative forecasts (eg, multiple time-horizons) and static file-based forecasts. + target_event_rate_forecast1p, n_fore1 = forecast.target_event_rates(observed_catalog, scale=scale) + target_event_rate_forecast2p, n_fore2 = benchmark_forecast.target_event_rates(observed_catalog, scale=scale) + + target_event_rate_forecast1 = forecast.data.ravel()[np.unique(np.nonzero(observed_catalog.spatial_magnitude_counts().ravel()))] + target_event_rate_forecast2 = benchmark_forecast.data.ravel()[np.unique(np.nonzero(observed_catalog.spatial_magnitude_counts(). + ravel()))] + + # call the primative version operating on ndarray + out = matrix_binary_t_test(target_event_rate_forecast1, target_event_rate_forecast2, observed_catalog.event_count, n_fore1, n_fore2, + observed_catalog, + alpha=alpha) + + # storing this for later + result = EvaluationResult() + result.name = 'binary paired T-Test' + result.test_distribution = (out['ig_lower'], out['ig_upper']) + result.observed_statistic = out['information_gain'] + result.quantile = (out['t_statistic'], out['t_critical']) + result.sim_name = (forecast.name, benchmark_forecast.name) + result.obs_name = observed_catalog.name + result.status = 'normal' + result.min_mw = np.min(forecast.magnitudes) + return result diff --git a/csep/core/catalogs.py b/csep/core/catalogs.py index 25603ad8..3b0fecfe 100644 --- a/csep/core/catalogs.py +++ b/csep/core/catalogs.py @@ -737,16 +737,16 @@ def spatial_magnitude_counts(self, mag_bins=None, tol=0.00001): """ Return counts of events in space-magnitude region. We figure out the index of the polygons and create a map that relates the spatial coordinate in the - Cartesian grid with with the polygon in region. + Cartesian grid with the polygon in region. Args: - mag_bins: magnitude bins (optional). tries to use magnitue bins associated with region + mag_bins (list, numpy.array): magnitude bins (optional), if empty tries to use magnitude bins associated with region + tol (float): tolerance for comparisons within magnitude bins Returns: output: unnormalized event count in each bin, 1d ndarray where index corresponds to midpoints """ - # make sure region is specified with catalog if self.region is None: raise CSEPCatalogException("Cannot create binned rates without region information.") @@ -784,8 +784,8 @@ def get_bvalue(self, mag_bins=None, return_error=True): If that fails, uses the default magnitude bins provided in constants. Args: - reterr (bool): returns errors mag_bins (list or array_like): monotonically increasing set of magnitude bin edges + return_error (bool): returns errors Returns: bval (float): b-value @@ -824,6 +824,10 @@ def p(): else: return bval + def b_positive(self): + """ Implements the b-positive indicator from Nicholas van der Elst """ + pass + def plot(self, ax=None, show=False, extent=None, set_global=False, plot_args=None): """ Plot catalog according to plate-carree projection @@ -1028,9 +1032,10 @@ def read_catalog_line(line): raise ValueError( "catalog_id should be monotonically increasing and events should be ordered by catalog_id") # yield final catalog, note: since this is just loading catalogs, it has no idea how many should be there - yield cls(data=events, catalog_id=prev_id, **kwargs) + cat = cls(data=events, catalog_id=prev_id, **kwargs) + yield cat - if os.path.isdir(filename): + elif os.path.isdir(filename): raise NotImplementedError("reading from directory or batched files not implemented yet!") @classmethod diff --git a/csep/core/forecasts.py b/csep/core/forecasts.py index 68cf27e7..4b2f863f 100644 --- a/csep/core/forecasts.py +++ b/csep/core/forecasts.py @@ -649,7 +649,7 @@ def magnitude_counts(self): self.get_expected_rates() return self.expected_rates.magnitude_counts() - def get_event_counts(self): + def get_event_counts(self, verbose=True): """ Returns a numpy array containing the number of event counts for each catalog. Note: This function can take a while to compute if called without already iterating through a forecast that @@ -661,7 +661,13 @@ def get_event_counts(self): """ if len(self._event_counts) == 0: # event counts is filled while iterating over the catalog - for _ in self: + t0 = time.time() + for i, _ in enumerate(self): + 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:.2f} seconds', flush=True) pass return numpy.array(self._event_counts) @@ -697,7 +703,7 @@ def get_expected_rates(self, verbose=False): 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) + print(f'Processed {i + 1} catalogs in {t1 - t0:.3f} seconds', flush=True) # after we iterate through the catalogs, we know self.n_cat data = data / self.n_cat self.expected_rates = GriddedForecast(self.start_time, self.end_time, data=data, region=self.region, diff --git a/csep/core/regions.py b/csep/core/regions.py index e9551f03..04795c01 100644 --- a/csep/core/regions.py +++ b/csep/core/regions.py @@ -631,14 +631,14 @@ def get_index_of(self, lons, lats): raise ValueError("at least one lon and lat pair contain values that are outside of the valid region.") if numpy.any(self.bbox_mask[idy, idx] == 1): raise ValueError("at least one lon and lat pair contain values that are outside of the valid region.") - return self.idx_map[idy,idx].astype(numpy.int) + return self.idx_map[idy, idx].astype(numpy.int64) def get_location_of(self, indices): """ Returns the polygon associated with the index idx. Args: - idx: index of polygon in region + indices: index of polygon in region Returns: Polygon diff --git a/csep/utils/basic_types.py b/csep/utils/basic_types.py index 5e1f5aa5..c8c5168e 100644 --- a/csep/utils/basic_types.py +++ b/csep/utils/basic_types.py @@ -51,7 +51,7 @@ def add(self, data): # need to know the range of the data to be inserted on discretized grid (min, max) # this is to determine the discretization of the data - eps=np.finfo(np.float).eps + eps = np.finfo(np.float64).eps disc_min = np.floor((data_min+eps-self.anchor)*self.rec_dh)/self.rec_dh+self.anchor disc_max = np.ceil((data_max+eps-self.anchor)*self.rec_dh)/self.rec_dh+self.anchor diff --git a/csep/utils/calc.py b/csep/utils/calc.py index e0071f47..dc208cda 100644 --- a/csep/utils/calc.py +++ b/csep/utils/calc.py @@ -79,9 +79,9 @@ def bin1d_vec(p, bins, tol=None, right_continuous=False): else: h = bins[1] - bins[0] - a0_tol = numpy.abs(a0) * numpy.finfo(numpy.float).eps - h_tol = numpy.abs(h) * numpy.finfo(numpy.float).eps - p_tol = numpy.abs(p) * numpy.finfo(numpy.float).eps + a0_tol = numpy.abs(a0) * numpy.finfo(numpy.float64).eps + h_tol = numpy.abs(h) * numpy.finfo(numpy.float64).eps + p_tol = numpy.abs(p) * numpy.finfo(numpy.float64).eps # absolute tolerance if tol is None: @@ -97,7 +97,7 @@ def bin1d_vec(p, bins, tol=None, right_continuous=False): try: idx[(idx < 0)] = -1 idx[(idx >= len(bins) - 1)] = len(bins) - 1 - except (TypeError): + except TypeError: if idx >= len(bins) - 1: idx = len(bins) - 1 if idx < 0: @@ -105,12 +105,12 @@ def bin1d_vec(p, bins, tol=None, right_continuous=False): else: try: idx[((idx < 0) | (idx >= len(bins)))] = -1 - except (TypeError): + except TypeError: if idx < 0 or idx >= len(bins): idx = -1 try: - idx = idx.astype(numpy.int) - except (AttributeError): + idx = idx.astype(numpy.int64) + except AttributeError: idx = int(idx) return idx diff --git a/csep/utils/comcat.py b/csep/utils/comcat.py index da0d2558..ecaea5f9 100644 --- a/csep/utils/comcat.py +++ b/csep/utils/comcat.py @@ -1,8 +1,9 @@ # python imports -from datetime import datetime, timedelta +from datetime import datetime, timedelta, timezone from urllib import request -from urllib.error import HTTPError +from urllib.error import HTTPError, URLError from urllib.parse import urlparse, urlencode +import ssl import json import time from collections import OrderedDict @@ -293,6 +294,39 @@ def _search(**newargs): except Exception as msg: raise Exception( 'Error downloading data from url %s. "%s".' % (url, msg)) + + except ssl.SSLCertVerificationError as SSLe: + # Fails to verify SSL certificate, when there is a hostname mismatch + if SSLe.verify_code == 62: + try: + context = ssl._create_unverified_context() + fh = request.urlopen(url, timeout=TIMEOUT, context=context) + data = fh.read().decode('utf8') + fh.close() + jdict = json.loads(data) + events = [] + for feature in jdict['features']: + events.append(SummaryEvent(feature)) + except Exception as msg: + raise Exception( + 'Error downloading data from url %s. "%s".' % (url, msg)) + + except URLError as URLe: + # Fails to verify SSL certificate, when there is a hostname mismatch + if isinstance(URLe.reason, ssl.SSLCertVerificationError) and URLe.reason.verify_code == 62: + try: + context = ssl._create_unverified_context() + fh = request.urlopen(url, timeout=TIMEOUT, context=context) + data = fh.read().decode('utf8') + fh.close() + jdict = json.loads(data) + events = [] + for feature in jdict['features']: + events.append(SummaryEvent(feature)) + except Exception as msg: + raise Exception( + 'Error downloading data from url %s. "%s".' % (url, msg)) + except Exception as msg: raise Exception( 'Error downloading data from url %s. "%s".' % (url, msg)) @@ -358,7 +392,11 @@ def id(self): Returns: str: Authoritative origin ID. """ - return self._jdict['id'] + ## comcat has an id key in each feature, whereas bsi has eventId within the properties dict + try: + return self._jdict['id'] + except: + return self._jdict['properties']['eventId'] @property def time(self): @@ -367,6 +405,10 @@ def time(self): datetime: Authoritative origin time. """ time_in_msec = self._jdict['properties']['time'] + # Comcat gives the event time in a ms timestamp, whereas bsi in datetime isoformat + if isinstance(time_in_msec, str): + event_dtime = datetime.fromisoformat(time_in_msec).replace(tzinfo=timezone.utc) + time_in_msec = event_dtime.timestamp() * 1000 time_in_sec = time_in_msec // 1000 msec = time_in_msec - (time_in_sec * 1000) dtime = datetime.utcfromtimestamp(time_in_sec) diff --git a/csep/utils/plots.py b/csep/utils/plots.py index ca75fa0d..f4369fb0 100644 --- a/csep/utils/plots.py +++ b/csep/utils/plots.py @@ -6,6 +6,7 @@ import scipy.stats import matplotlib +import matplotlib.lines from matplotlib import cm from matplotlib.collections import PatchCollection import matplotlib.pyplot as pyplot @@ -1512,6 +1513,7 @@ def plot_comparison_test(results_t, results_w=None, axes=None, plot_args=None): return ax + def plot_poisson_consistency_test(eval_results, normalize=False, one_sided_lower=False, axes=None, plot_args=None, show=False): """ Plots results from CSEP1 tests following the CSEP1 convention. @@ -1555,7 +1557,7 @@ def plot_poisson_consistency_test(eval_results, normalize=False, one_sided_lower figsize= plot_args.get('figsize', None) title = plot_args.get('title', results[0].name) title_fontsize = plot_args.get('title_fontsize', None) - xlabel = plot_args.get('xlabel', 'X') + xlabel = plot_args.get('xlabel', '') xlabel_fontsize = plot_args.get('xlabel_fontsize', None) xticks_fontsize = plot_args.get('xticks_fontsize', None) ylabel_fontsize = plot_args.get('ylabel_fontsize', None) @@ -1565,6 +1567,7 @@ def plot_poisson_consistency_test(eval_results, normalize=False, one_sided_lower hbars = plot_args.get('hbars', True) tight_layout = plot_args.get('tight_layout', True) percentile = plot_args.get('percentile', 95) + plot_mean = plot_args.get('mean', False) if axes is None: fig, ax = pyplot.subplots(figsize=figsize) @@ -1578,6 +1581,7 @@ def plot_poisson_consistency_test(eval_results, normalize=False, one_sided_lower if res.test_distribution[0] == 'poisson': plow = scipy.stats.poisson.ppf((1 - percentile/100.)/2., res.test_distribution[1]) phigh = scipy.stats.poisson.ppf(1 - (1 - percentile/100.)/2., res.test_distribution[1]) + mean = res.test_distribution[1] observed_statistic = res.observed_statistic # empirical distributions else: @@ -1594,12 +1598,14 @@ def plot_poisson_consistency_test(eval_results, normalize=False, one_sided_lower else: plow = numpy.percentile(test_distribution, (100 - percentile)/2.) phigh = numpy.percentile(test_distribution, 100 - (100 - percentile)/2.) + mean = numpy.mean(res.test_distribution) if not numpy.isinf(observed_statistic): # Check if test result does not diverges - low = observed_statistic - plow - high = phigh - observed_statistic - ax.errorbar(observed_statistic, index, xerr=numpy.array([[low, high]]).T, - fmt=_get_marker_style(observed_statistic, (plow, phigh), one_sided_lower), + percentile_lims = numpy.array([[mean - plow, phigh - mean]]).T + ax.plot(observed_statistic, index, + _get_marker_style(observed_statistic, (plow, phigh), one_sided_lower)) + ax.errorbar(mean, index, xerr=percentile_lims, + fmt='ko'*plot_mean, capsize=capsize, linewidth=linewidth, ecolor=color) # determine the limits to use xlims.append((plow, phigh, observed_statistic)) @@ -1883,15 +1889,16 @@ def add_labels_for_publication(figure, style='bssa', labelsize=16): ax.annotate(f'({annot})', (0.025, 1.025), xycoords='axes fraction', fontsize=labelsize) return - -def plot_consistency_test(eval_results, normalize=False, one_sided_lower=True, plot_args=None, variance=None): + + +def plot_consistency_test(eval_results, normalize=False, axes=None, one_sided_lower=False, variance=None, plot_args=None, show=False): """ Plots results from CSEP1 tests following the CSEP1 convention. Note: All of the evaluations should be from the same type of evaluation, otherwise the results will not be comparable on the same figure. Args: - results (list): Contains the tests results :class:`csep.core.evaluations.EvaluationResult` (see note above) + eval_results (list): Contains the tests results :class:`csep.core.evaluations.EvaluationResult` (see note above) normalize (bool): select this if the forecast likelihood should be normalized by the observed likelihood. useful for plotting simulation based simulation tests. one_sided_lower (bool): select this if the plot should be for a one sided test @@ -1924,8 +1931,10 @@ def plot_consistency_test(eval_results, normalize=False, one_sided_lower=True, p # Parse plot arguments. More can be added here if plot_args is None: plot_args = {} - figsize= plot_args.get('figsize', (7,8)) - xlabel = plot_args.get('xlabel', 'X') + figsize= plot_args.get('figsize', None) + title = plot_args.get('title', results[0].name) + title_fontsize = plot_args.get('title_fontsize', None) + xlabel = plot_args.get('xlabel', '') xlabel_fontsize = plot_args.get('xlabel_fontsize', None) xticks_fontsize = plot_args.get('xticks_fontsize', None) ylabel_fontsize = plot_args.get('ylabel_fontsize', None) @@ -1935,8 +1944,14 @@ def plot_consistency_test(eval_results, normalize=False, one_sided_lower=True, p hbars = plot_args.get('hbars', True) tight_layout = plot_args.get('tight_layout', True) percentile = plot_args.get('percentile', 95) + plot_mean = plot_args.get('mean', False) + + if axes is None: + fig, ax = pyplot.subplots(figsize=figsize) + else: + ax = axes + fig = ax.get_figure() - fig, ax = plt.subplots(figsize=figsize) xlims = [] for index, res in enumerate(results): @@ -1944,6 +1959,7 @@ def plot_consistency_test(eval_results, normalize=False, one_sided_lower=True, p if res.test_distribution[0] == 'poisson': plow = scipy.stats.poisson.ppf((1 - percentile/100.)/2., res.test_distribution[1]) phigh = scipy.stats.poisson.ppf(1 - (1 - percentile/100.)/2., res.test_distribution[1]) + mean = res.test_distribution[1] observed_statistic = res.observed_statistic elif res.test_distribution[0] == 'negative_binomial': @@ -1952,8 +1968,8 @@ def plot_consistency_test(eval_results, normalize=False, one_sided_lower=True, p mean = res.test_distribution[1] upsilon = 1.0 - ((var - mean) / var) tau = (mean**2 /(var - mean)) - phigh = nbinom.ppf((1 - percentile/100.)/2., tau, upsilon) - plow = nbinom.ppf(1 - (1 - percentile/100.)/2., tau, upsilon) + phigh = scipy.stats.nbinom.ppf((1 - percentile/100.)/2., tau, upsilon) + plow = scipy.stats.nbinom.ppf(1 - (1 - percentile/100.)/2., tau, upsilon) # empirical distributions else: @@ -1970,13 +1986,15 @@ def plot_consistency_test(eval_results, normalize=False, one_sided_lower=True, p else: plow = numpy.percentile(test_distribution, 2.5) phigh = numpy.percentile(test_distribution, 97.5) + mean = numpy.mean(res.test_distribution) if not numpy.isinf(observed_statistic): # Check if test result does not diverges - low = observed_statistic - plow - high = phigh - observed_statistic - ax.errorbar(observed_statistic, index, xerr=numpy.array([[low, high]]).T, - fmt=_get_marker_style(observed_statistic, (plow, phigh), one_sided_lower), - capsize=4, linewidth=linewidth, ecolor=color, markersize = 10, zorder=1) + percentile_lims = numpy.array([[mean - plow, phigh - mean]]).T + ax.plot(observed_statistic, index, + _get_marker_style(observed_statistic, (plow, phigh), one_sided_lower)) + ax.errorbar(mean, index, xerr=percentile_lims, + fmt='ko'*plot_mean, + capsize=capsize, linewidth=linewidth, ecolor=color) # determine the limits to use xlims.append((plow, phigh, observed_statistic)) # we want to only extent the distribution where it falls outside of it in the acceptable tail @@ -1998,133 +2016,137 @@ def plot_consistency_test(eval_results, normalize=False, one_sided_lower=True, p except ValueError: raise ValueError('All EvaluationResults have infinite observed_statistics') ax.set_yticks(numpy.arange(len(results))) - ax.set_yticklabels([res.sim_name for res in results], fontsize=14) + ax.set_yticklabels([res.sim_name for res in results], fontsize=ylabel_fontsize) ax.set_ylim([-0.5, len(results)-0.5]) if hbars: yTickPos = ax.get_yticks() if len(yTickPos) >= 2: ax.barh(yTickPos, numpy.array([99999] * len(yTickPos)), left=-10000, height=(yTickPos[1] - yTickPos[0]), color=['w', 'gray'], alpha=0.2, zorder=0) - ax.set_xlabel(xlabel, fontsize=14) - ax.tick_params(axis='x', labelsize=13) + ax.set_title(title, fontsize=title_fontsize) + ax.set_xlabel(xlabel, fontsize=xlabel_fontsize) + ax.tick_params(axis='x', labelsize=xticks_fontsize) if tight_layout: ax.figure.tight_layout() fig.tight_layout() + + if show: + pyplot.show() + return ax -def plot_pvalues_and_intervals(test_results, var=None): +def plot_pvalues_and_intervals(test_results, ax, var=None): + """ Plots p-values and intervals for a list of Poisson or NBD test results + + Args: + test_results (list): list of EvaluationResults for N-test. All tests should use the same distribution + (ie Poisson or NBD). + ax (matplotlib.axes.Axes.axis): axes to use for plot. create using matplotlib + var (float): variance of the NBD distribution. Must be used for NBD plots. + + Returns: + ax (matplotlib.axes.Axes.axis): axes handle containing this plot + + Raises: + ValueError: throws error if NBD tests are supplied without a variance + """ - variance= var + variance = var percentile = 97.5 p_values = [] - + + # Differentiate between N-tests and other consistency tests if test_results[0].name == 'NBD N-Test' or test_results[0].name == 'Poisson N-Test': - legend_elements = [Line2D([0], [0], marker='o', color='red', lw=0, label=r'p < 10e-5', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='#FF7F50', lw=0, label=r'10e-5 $\leq$ p < 10e-4', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='gold', lw=0, label=r'10e-4 $\leq$ p < 10e-3', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='white', lw=0, label=r'10e-3 $\leq$ p < 0.0125', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='skyblue', lw=0, label=r'0.0125 $\leq$ p < 0.025', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='blue', lw=0, label=r'p $\geq$ 0.025', markersize=10, markeredgecolor='k')] - + legend_elements = [matplotlib.lines.Line2D([0], [0], marker='o', color='red', lw=0, label=r'p < 10e-5', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='#FF7F50', lw=0, label=r'10e-5 $\leq$ p < 10e-4', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='gold', lw=0, label=r'10e-4 $\leq$ p < 10e-3', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='white', lw=0, label=r'10e-3 $\leq$ p < 0.0125', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='skyblue', lw=0, label=r'0.0125 $\leq$ p < 0.025', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='blue', lw=0, label=r'p $\geq$ 0.025', markersize=10, markeredgecolor='k')] ax.legend(handles=legend_elements, loc=4, fontsize=13, edgecolor='k') - + # Act on Negative binomial tests if test_results[0].name == 'NBD N-Test': - + if var is None: + raise ValueError("var must not be None if N-tests use the NBD distribution.") + for i in range(len(test_results)): mean = test_results[i].test_distribution[1] upsilon = 1.0 - ((variance - mean) / variance) tau = (mean**2 /(variance - mean)) - phigh97 = nbinom.ppf((1 - percentile/100.)/2., tau, upsilon) - plow97= nbinom.ppf(1 - (1 - percentile/100.)/2., tau, upsilon) + phigh97 = scipy.stats.nbinom.ppf((1 - percentile/100.)/2., tau, upsilon) + plow97 = scipy.stats.nbinom.ppf(1 - (1 - percentile/100.)/2., tau, upsilon) low97 = test_results[i].observed_statistic - plow97 high97 = phigh97 - test_results[i].observed_statistic ax.errorbar(test_results[i].observed_statistic, (len(test_results)-1) - i, xerr=numpy.array([[low97, high97]]).T, capsize=4, color='slategray', alpha=1.0, zorder=0) p_values.append(test_results[i].quantile[1] * 2.0) # Calculated p-values according to Meletti et al., (2021) - if p_values[i] < 10e-5: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'red', markersize= 8, zorder=2) - - if p_values[i] >= 10e-5 and p_values[i] < 10e-4: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = '#FF7F50', markersize= 8, zorder=2) - - if p_values[i] >= 10e-4 and p_values[i] < 10e-3: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'gold', markersize= 8, zorder=2) - - if p_values[i] >= 10e-3 and p_values[i] < 0.0125: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'white', markersize= 8, zorder=2) - - if p_values[i] >= 0.0125 and p_values[i] < 0.025: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'skyblue', markersize= 8, zorder=2) - - if p_values[i] >= 0.025: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'blue', markersize= 8, zorder=2) - + if p_values[i] < 10e-5: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='red', markersize=8, zorder=2) + if p_values[i] >= 10e-5 and p_values[i] < 10e-4: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='#FF7F50', markersize=8, zorder=2) + if p_values[i] >= 10e-4 and p_values[i] < 10e-3: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='gold', markersize=8, zorder=2) + if p_values[i] >= 10e-3 and p_values[i] < 0.0125: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='white', markersize=8, zorder=2) + if p_values[i] >= 0.0125 and p_values[i] < 0.025: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='skyblue', markersize=8, zorder=2) + if p_values[i] >= 0.025: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='blue', markersize=8, zorder=2) + # Act on Poisson N-test if test_results[0].name == 'Poisson N-Test': - for i in range(len(test_results)): plow97 = scipy.stats.poisson.ppf((1 - percentile/100.)/2., test_results[i].test_distribution[1]) - phigh97 = scipy.stats.poisson.ppf(1- (1 - percentile/100.)/2., test_results[i].test_distribution[1]) + phigh97 = scipy.stats.poisson.ppf(1 - (1 - percentile/100.)/2., test_results[i].test_distribution[1]) low97 = test_results[i].observed_statistic - plow97 high97 = phigh97 - test_results[i].observed_statistic ax.errorbar(test_results[i].observed_statistic, (len(test_results)-1) - i, xerr=numpy.array([[low97, high97]]).T, capsize=4, color='slategray', alpha=1.0, zorder=0) p_values.append(test_results[i].quantile[1] * 2.0) - - if p_values[i] < 10e-5: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'red', markersize= 8, zorder=2) - - if p_values[i] >= 10e-5 and p_values[i] < 10e-4: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = '#FF7F50', markersize= 8, zorder=2) - - if p_values[i] >= 10e-4 and p_values[i] < 10e-3: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'gold', markersize= 8, zorder=2) - - if p_values[i] >= 10e-3 and p_values[i] < 0.0125: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'white', markersize= 8, zorder=2) - - if p_values[i] >= 0.0125 and p_values[i] < 0.025: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'skyblue', markersize= 8, zorder=2) - - if p_values[i] >= 0.025: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'blue', markersize= 8, zorder=2) - - + if p_values[i] < 10e-5: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='red', markersize=8, zorder=2) + elif p_values[i] >= 10e-5 and p_values[i] < 10e-4: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='#FF7F50', markersize=8, zorder=2) + elif p_values[i] >= 10e-4 and p_values[i] < 10e-3: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='gold', markersize=8, zorder=2) + elif p_values[i] >= 10e-3 and p_values[i] < 0.0125: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='white', markersize=8, zorder=2) + elif p_values[i] >= 0.0125 and p_values[i] < 0.025: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='skyblue', markersize=8, zorder=2) + elif p_values[i] >= 0.025: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='blue', markersize=8, zorder=2) + # Operate on all other consistency tests else: for i in range(len(test_results)): plow97 = numpy.percentile(test_results[i].test_distribution, 2.5) - phigh97 = numpy.percentile(test_results[i].test_distribution, 100) + phigh97 = numpy.percentile(test_results[i].test_distribution, 97.5) low97 = test_results[i].observed_statistic - plow97 high97 = phigh97 - test_results[i].observed_statistic ax.errorbar(test_results[i].observed_statistic, (len(test_results)-1) -i, xerr=numpy.array([[low97, high97]]).T, capsize=4, color='slategray', alpha=1.0, zorder=0) p_values.append(test_results[i].quantile) - if p_values[i] < 10e-5: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'red', markersize= 8, zorder=2) - - if p_values[i] >= 10e-5 and p_values[i] < 10e-4: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = '#FF7F50', markersize= 8, zorder=2) - - if p_values[i] >= 10e-4 and p_values[i] < 10e-3: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'gold', markersize= 8, zorder=2) - - if p_values[i] >= 10e-3 and p_values[i] < 0.025: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'white', markersize= 8, zorder=2) - - if p_values[i] >= 0.025 and p_values[i] < 0.05: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'skyblue', markersize= 8, zorder=2) - - if p_values[i] >= 0.05: - ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color = 'blue', markersize= 8, zorder=2) + if p_values[i] < 10e-5: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='red', markersize=8, zorder=2) + elif p_values[i] >= 10e-5 and p_values[i] < 10e-4: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='#FF7F50', markersize=8, zorder=2) + elif p_values[i] >= 10e-4 and p_values[i] < 10e-3: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='gold', markersize=8, zorder=2) + elif p_values[i] >= 10e-3 and p_values[i] < 0.025: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='white', markersize=8, zorder=2) + elif p_values[i] >= 0.025 and p_values[i] < 0.05: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='skyblue', markersize=8, zorder=2) + elif p_values[i] >= 0.05: + ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='blue', markersize=8, zorder=2) - legend_elements = [Line2D([0], [0], marker='o', color='red', lw=0, label=r'p < 10e-5', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='#FF7F50', lw=0, label=r'10e-5 $\leq$ p < 10e-4', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='gold', lw=0, label=r'10e-4 $\leq$ p < 10e-3', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='white', lw=0, label=r'10e-3 $\leq$ p < 0.025', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='skyblue', lw=0, label=r'0.025 $\leq$ p < 0.05', markersize=10, markeredgecolor='k'), - Line2D([0], [0], marker='o', color='blue', lw=0, label=r'p $\geq$ 0.05', markersize=10, markeredgecolor='k')] + legend_elements = [ + matplotlib.lines.Line2D([0], [0], marker='o', color='red', lw=0, label=r'p < 10e-5', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='#FF7F50', lw=0, label=r'10e-5 $\leq$ p < 10e-4', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='gold', lw=0, label=r'10e-4 $\leq$ p < 10e-3', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='white', lw=0, label=r'10e-3 $\leq$ p < 0.025', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='skyblue', lw=0, label=r'0.025 $\leq$ p < 0.05', markersize=10, markeredgecolor='k'), + matplotlib.lines.Line2D([0], [0], marker='o', color='blue', lw=0, label=r'p $\geq$ 0.05', markersize=10, markeredgecolor='k')] ax.legend(handles=legend_elements, loc=4, fontsize=13, edgecolor='k') diff --git a/csep/utils/readers.py b/csep/utils/readers.py index b5352279..4569634b 100644 --- a/csep/utils/readers.py +++ b/csep/utils/readers.py @@ -655,8 +655,9 @@ def jma_csv(fname): return events def _query_comcat(start_time, end_time, min_magnitude=2.50, - min_latitude=31.50, max_latitude=43.00, - min_longitude=-125.40, max_longitude=-113.10, extra_comcat_params=None): + min_latitude=31.50, max_latitude=43.00, + min_longitude=-125.40, max_longitude=-113.10, + max_depth=1000, extra_comcat_params=None): """ Return eventlist from ComCat web service. @@ -668,6 +669,7 @@ def _query_comcat(start_time, end_time, min_magnitude=2.50, max_latitude (float): maximum latitude of query min_longitude (float): minimum longitude of query max_longitude (float): maximum longitude of query + max_depth (float): maximum depth of query extra_comcat_params (dict): additional parameters to pass to comcat search function Returns: @@ -679,10 +681,30 @@ def _query_comcat(start_time, end_time, min_magnitude=2.50, eventlist = search(minmagnitude=min_magnitude, minlatitude=min_latitude, maxlatitude=max_latitude, minlongitude=min_longitude, maxlongitude=max_longitude, - starttime=start_time, endtime=end_time, **extra_comcat_params) + starttime=start_time, endtime=end_time, + maxdepth=max_depth, **extra_comcat_params) return eventlist +def _query_bsi(start_time, end_time, min_magnitude=2.50, + min_latitude=32.0, max_latitude=50.0, + min_longitude=2.0, max_longitude=21.0, + max_depth=1000, extra_bsi_params=None): + """ + Queries INGV Bulletino Sismico Italiano, revised version. + :return: csep.core.Catalog object + """ + extra_bsi_params = extra_bsi_params or {} + bsi_host = 'webservices.rm.ingv.it' + extra_bsi_params.update({'host': bsi_host, 'limit': 15000, 'offset': 0}) + # get eventlist from Comcat + eventlist = search(minmagnitude=min_magnitude, + minlatitude=min_latitude, maxlatitude=max_latitude, + minlongitude=min_longitude, maxlongitude=max_longitude, + maxdepth=max_depth, + starttime=start_time, endtime=end_time, **extra_bsi_params) + + return eventlist def _parse_datetime_to_zmap(date, time): """ Helping function to return datetime in zmap format. diff --git a/docs/reference/api_reference.rst b/docs/reference/api_reference.rst index 9f8a6b36..6e3ce176 100644 --- a/docs/reference/api_reference.rst +++ b/docs/reference/api_reference.rst @@ -16,6 +16,7 @@ Loading catalogs and forecasts load_stochastic_event_sets load_catalog query_comcat + query_bsi load_gridded_forecast load_catalog_forecast diff --git a/examples/tutorials/catalog_filtering.py b/examples/tutorials/catalog_filtering.py index 406ee565..a21937ff 100644 --- a/examples/tutorials/catalog_filtering.py +++ b/examples/tutorials/catalog_filtering.py @@ -30,8 +30,9 @@ # Load catalog # ------------ # -# PyCSEP provides access to the ComCat web API using :func:`csep.query_comcat`. This function requires a -# :class:`datetime.datetime` to specify the start and end dates. +# PyCSEP provides access to the ComCat web API using :func:`csep.query_comcat` and to the Bollettino Sismico Italiano +# API using :func:`csep.query_bsi`. These functions require a :class:`datetime.datetime` to specify the start and end +# dates. start_time = csep.utils.time_utils.strptime_to_utc_datetime('2019-01-01 00:00:00.0') end_time = csep.utils.time_utils.utc_now_datetime() diff --git a/notebooks/CSEP_tests.ipynb b/notebooks/CSEP_tests.ipynb index da30f86c..b377641c 100644 --- a/notebooks/CSEP_tests.ipynb +++ b/notebooks/CSEP_tests.ipynb @@ -72,7 +72,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Fetched ComCat catalog in 5.9399449825286865 seconds.\n", + "Fetched ComCat catalog in 5.9891064167022705 seconds.\n", "\n", "Downloaded catalog from ComCat with following parameters\n", "Start Date: 2010-01-10 00:27:39.320000+00:00\n", @@ -134,6 +134,33 @@ "\n", "where $\\lambda(m_i, s_j)$ and $\\omega(m_i, s_j)$ are the expected counts from the forecast and observed counts in cell $m_i, s_j$ respectively. We can calculate the likelihood directly given the forecast and discretised observations.\n", "\n", + "#### Binary consistency tests\n", + "\n", + "It is known that the Poisson distribution insufficiently captures the variability of earthquakes, as it could be dispersed with respect to the true distribution of seismicity, especially in the presence of clustered earthquakes (Werner and Sornette, 2008; Lombardi and Marzocchi, 2010; Nandan et al., 2019). Therefore, CSEP has recently introduced a set of new consistency tests based on the negative binomial distribution (NBD) and a binary likelihood function, with the aim of reducing the sensitivity of traditional Poisson-based tests to clustering (Bayona et al., 2022).\n", + "\n", + "The NBD has been shown to better describe non-declustered seismicity, because it has a greater variance than the Poisson distribution to account for spatiotemporal earthquake clustering (Werner et al., 2011a,b). Hence, one can fit a NBD to the rates of observed earthquakes in each forecast's testing region to evaluate number consistencies between forecasts and observations. The probability mass function of an NBD is defined by:\n", + "\n", + "$P (\\omega\\hspace{0.1cm}|\\hspace{0.1cm}\\tau, \\nu) = \\frac{\\Gamma(\\tau + \\omega)}{\\Gamma(\\tau)\\omega!} \\nu^\\omega (1-\\nu)^{\\omega},$\n", + "\n", + "where $\\omega = 0,1,2,...$ is the number of events, $\\tau > 0$ and $0 \\leq \\nu \\leq 1$ are parameters, and $\\Gamma$ is the Gamma function. The mean $\\mu$ and variance $\\sigma^2$ of the NBD are given by\n", + "\n", + "$\\mu = \\tau\\frac{1 - \\nu}{\\nu}; \\hspace{0.1cm} \\sigma^2=\\tau \\frac{1 - \\nu}{\\nu^2}$.\n", + "\n", + "Following the approach by Werner et al. (2011b), one can use the expected number earthquakes by the forecast as the mean and estimate the variance from historical catalogs.\n", + "\n", + "The binary likelihood function can also make the assumption that earthquakes are Poisson distributed less critical by calculating the probability of observing one or more earthquakes in a forecast cell/bin, rather than the likelihood of observing $\\omega =$ 1,2, .., n events. The probability of observing $\\omega = $ 0 earthquakes given an expected number/rate $\\lambda$, assuming a Poisson distribution, is P$_\\mathrm{0}$ = exp(-$\\lambda)$, while the probability of observing one or more events is P$_\\mathrm{_1}$ = 1 - P$_\\mathrm{0}$ = 1 - exp($-\\lambda)$. Accordingly, one can compute the log-likelihood score in each cell/bin based on the binary (or Bernoulli) distribution as:\n", + "\n", + "$$ L(\\boldsymbol{\\Omega} | \\boldsymbol{\\Lambda}) = X(m_i, s_j) \\hspace{0.1cm} \\mathrm{log} \\hspace{0.1cm} \\mathrm{P_1} + (1 -X(m_i, s_j)) \\hspace{0.1cm} \\mathrm{log} \\hspace{0.1cm} \\mathrm{P_0} = X(m_i, s_j) \\hspace{0.1cm} \\mathrm{log} \\hspace{0.1cm} (1-\\mathrm{exp}(-\\lambda(m_i, s_j))) + (1 -X_i) \\hspace{0.1cm} \\mathrm{log} \\hspace{0.1cm} (\\mathrm{exp}(-\\lambda(m_i, s_j)))$$\n", + "\n", + "where the first term represents the contribution to the score if a bin $(m_i, s_j)$ contains one or more events, i.e. $X(m_i, s_j)$ = 1, and the second term is the contribution if that cell/bin contains no earthquakes, i.e. $X(m_i, s_j)$ = 0.\n", + "\n", + "Similar to the example above, the joint likelihood of observing events in individual bins, based on the binary likelihood function, is the sum of each individual log score, i.e.:\n", + "\n", + "\n", + "$$ L(\\boldsymbol{\\Omega} | \\boldsymbol{\\Lambda}) = \\sum_{m_i , s_j \\in \\boldsymbol{R}} \\hspace{0.05cm}[X(m_i, s_j) \\hspace{0.1cm} \\mathrm{log} \\hspace{0.1cm} (1-\\mathrm{exp}(-\\lambda(m_i, s_j))) + (1 - X(m_i, s_j)) \\hspace{0.1cm} \\mathrm{log} \\hspace{0.1cm} (\\mathrm{exp}(-\\lambda(m_i, s_j)))]$$\n", + "\n", + "\n", + "\n", " Forecast uncertainty \n", "\n", "A simulation based approach is used to account for uncertainty in the forecast. We simulate realizations of catalogs that are consistent with the forecast to obtain distributions of scores. In the pyCSEP package, as in the original CSEP tests, simulation is carried out using the cumulative probability density of the forecast obtained by ordering the rates in each bin. We shall call $F_{m_is_j}$ the cumulative probability density in cell $(m_i, s_j)$. The simulation approach then works as follows:\n", @@ -141,6 +168,7 @@ "* For each forecast bin, draw a random number $z$ from a uniform distribution between 0 and 1\n", "* Assign this event to a space-magnitude bin through the inverse cumulative density distribution at this point $F^{-1}_{m_i, s_j}(z)$\n", "* Iterate over all simulated events to generate a catalog containing $N_{sim}$ events consistent with the forecast\n", + "* In the case of binary tests, simulate catalogs that are consistent with the number of active cells/bins.\n", "\n", "For each of these tests, we can plot the distribution of likelihoods computed from theses simulated catalogs relative to the observations using the `plots.plot_poisson_consistency_test` function. We also calculate a quantile score to diagnose a particular forecast with repsect. The number of simulations can be supplied to the Poisson consistency test functions using the `num_simulations` argument: for best results we suggest 100,000 simulations to ensure convergence.\n", "\n", @@ -158,7 +186,9 @@ "\n", "Whether a forecast can be said to pass an evaluation depends on the significance level chosen for the testing process. The quantile score explicitly tells us something about the significance of the result: the observation is consistent with the forecast with $100(1-\\gamma)\\%$ confidence (Zechar, 2011). Low $\\gamma$ values demonstrate that the observed likelihood score is less than most of the simulated catalogs. The consistency tests, excluding the N-test, are considered to be one-sided tests: values which are too small are ruled inconsistent with the forecast, but very large values may not necessarily be inconsistent with the forecast and additional testing should be used to further clarify this (Schorlemmer et al, 2007). \n", "\n", - "Different CSEP experiments have used different sensitivity values. Schorlemmer et al (2010b) consider $\\gamma \\lt 0.05$ while the implementation in the Italian CSEP testing experiment uses $\\gamma$ < 0.01 (Taroni et al, 2018). However, the consistency tests are most useful as diagnostic tools where the quantile score assesses the level of consistency between observations and data. Temporal variations in seismicity make it difficult to formally reject a model from a consistency test over a single evaluation period." + "Different CSEP experiments have used different sensitivity values. Schorlemmer et al (2010b) consider $\\gamma \\lt 0.05$ while the implementation in the Italian CSEP testing experiment uses $\\gamma$ < 0.01 (Taroni et al, 2018). However, the consistency tests are most useful as diagnostic tools where the quantile score assesses the level of consistency between observations and data. Temporal variations in seismicity make it difficult to formally reject a model from a consistency test over a single evaluation period. \n", + "\n", + "Performing multiple tests among samples of the same distribution increases the probability of observing at least one rare events, and consequently, increases the probability of obtaining an apparently statistically significant result (Kato, 2019). Therefore, CSEP has recently included consistency test results at a Bonferroni-adjusted significance level $\\gamma_{B}$ = 0.05 /2 = 0.025 to control the overall type II error rate, that is, the false-positive rate or the \"false-inconsistency\" rate (see Bayona et al., 2022)." ] }, { @@ -187,7 +217,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQ8klEQVR4nO3de7BeVX3G8e9DQrkE5GIQxFuoBZEqpSZa8TahUrVYrResKFbRWiuiaGeoiGEg9VKlMLXeEMGZiFVq63hFxpqKAo4FQ4IhhFGpihRHUHBQBGwKYfWPd0VfD+eeQ87vnHw/M++c9a6919prnX3Oec7e7373m9YakiRVs8NsD0CSpNEYUJKkkgwoSVJJBpQkqSQDSpJUkgElSSrJgJIklWRASZJKMqAkTUmSPZO8blu31fbHgJK2M0lenOTcJG+eZhd7AtMNma1pq+2MASVtR5KcCDwKWAW8YJrdvBt4ZJL1Sc7s/b4syZpe9+EkC5IsSnJRkquTbEzy4tHaSmOJ9+KTtg9JdgX+BzgQWAjs1Vq7bhr9LAG+2Fp7TH/+aOAfgRe01u5OcjZwBXAn8KzW2l/39fYA9hpuK41n4WwPQNI2cwRwQ2vttv78luGFSb4C7DdKuxWttc+P0+/TgaXAlUkAdgF+ClwAnJXkDAah9PUke23lHLQdMaCk7cfTgNVjLWytHTnNfgOc31o75T4LkqXAUcC7kqwGPjbNbWg75GtQ0vbjScA3k+wDkOQBSQ6ZRj+/BHYfen4xcHSSB/V+907yiCT7A3e11j4OnAU8bpS20pg8gpK2A0mWA3cDK4BHJfkGcCFwzlT7aq39LMk3kmwEvtRa+7skpwKrk+zQt3MCsAdwZpJ7e93xo7WdiflpfvIiCUlSSZ7ikySVZEBJkkoyoCRJJRlQkqSSvIpPACxevLgtWbJktochaR5Zt27dra21fabb3oASAEuWLGHt2rWzPQxJ80iSG7amvaf4JEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgNG0rV66c7SFI2oa29e98WmvbdIOqadmyZW3t2rVTapMEf3400n5n7cdP7vzJfer3XbQvN5908yyMSDNlqr/zSda11pZNd3sTHkElWZJk4xQGtDLJSdMd0FA/b0qy61TXS/LWrd32JLZ5cJL1Sb6V5JFJXjqDff8wyeIZ6OeOmRiPNFWjhdN49dJYFs72AMbxJuDjwF1TXO+twD9MZUNJFrTWNk+hyfOAz7fWTk+yHHgpcMH9uL2yli9fPttDUDVHjL3InxdNxWRfg1qQ5Lwk1yZZnWSXfuTwH0nWJfl6koNHNkpySZL3JLksybeTPD7JZ5L8d5J39HUWJbkoydVJNiZ5cZITgf2BryX5Wl/vGUkuT3JVkk8l2W3kekneDezSj24+0du9LMmaXvfhJAt6/R1J3pbkm8Dho006yWlJruzjOjcDRzEIxVf3sb0beGrv/2+TLEhyZm+3Icnf9L6W9zFeAFwz2ryHNv2GPs9rtnxfk+yd5HO9zyuSHNrrd0uyqq+7IckLR8xhcf++PXuU+b0mydoka2+55ZZJ/ihI0jbSWhv3ASwB7gEO68//HXgZcDFwYK/7I+CrvbwSOKmXLwHO6OU3Aj8GHgzsBPwIeCDwQuC8oe3t0b/+EFjcy4uBy4BF/fnJwGkj1+vP7xgqPxq4ENixPz8beHkvN+AvJpj73kPlfwGeM8oclwNfHFrvNcCpvbwTsBY4oK93J3BAXzbevN/Qy68DPtLL7wdO7+U/Btb38hnAPw/1s9eW7wOwL/BN4E8m2s9Lly5tUzX48ZF+GysZ86G5baq/88DaNsHfnvEekz3Fd31rbX0vr2MQWk8CPpVkyzo7jdH2C/3rNcC1rbWbAJL8AHhYrz8ryRn9D/3XR+njicAhwDf69n4HuHwS4346sBS4srfbBfhpX7YZ+PQE7Y9I8mZgV2Bv4FoGgTeeZwCHJjm6P98DOBD4P2BNa+36Xj/evD/Tv64DXtDLT2EQarTWvprkgUn2AI4EjtnSsLV2Wy/uyOCfiBNaa5dOMGZJKmeyAbVpqLyZwX/mP2+tHTaFtveO6OdeYGFr7bokS4GjgHclWd1ae9uIPgL8Z2vtJZMc73C781trp4yy7H/bOK8DJdmZwRHXstbajUlWAjtPcptvaK19eUR/yxkcQQEwwby3fJ8285t99Ov/BIa0Xj/aZTX3MAi4ZwIGlLaZfRftO+ZVfNJUTPd9ULcD1yd5EUB/beYPptNRkv2Bu1prHwfOAh7XF/0S2L2XrwCenOT3eptdkxw0ynoAdyfZsZcvBo5O8qDebu8kj5jk0LaE0a1JdgOOHmO9kdv/MnD8ljEkOSjJopGNxpn3WC4Dju1tlwO3ttZuB1YDrx/qd69ebMCrgIOTvGWCvqfl9NNPvz+61Rx380k3005v93l4ifnct61/57fmKr5jgQ8lOZXB6aRPAldPo5/HAmcmuRe4Gzi+158LfCnJTa21I5IcB/xrki2nEk8Frhu5Xn++IclVrbVj+/hWJ9mh938CcMNEg2qt/TzJeQxOxf0QuHKMVTcA9yS5Gvgo8F4Gp0CvyuC84i0Mrvqb7LzHshJYlWQDgysWX9Hr3wF8MIO3AmwG/p5+irC1tjnJMcCFSW5vrZ09wTamxDfqStsX36irWTGdN+pK0nhyf79RV5Kk2VD5jbrbTJLPMrgUfNjJIy90kCRtOwYU0Fp7/myPQZL02zzFJ0kqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklpbU222NQAUluAW6Y7XEAi4FbZ3sQ25Dznd+29/k+orW2z3Q7M6BUSpK1rbVlsz2ObcX5zm/Od+t4ik+SVJIBJUkqyYBSNefO9gC2Mec7vznfreBrUJKkkjyCkiSVZEBJkkoyoDRrkrw9yYYk65OsTrL/0LJTknwvyXeTPHOofmmSa/qy9yXJ7Ix+apKcmeQ7fb6fTbJnr1+S5Ff9e7A+yTlDbebkXGHs+fZl82rfAiR5UZJrk9ybZNlQ/Xzdv6POty+buf3bWvPhY1YewAOGyicC5/TyIcDVwE7AAcD3gQV92RrgcCDAl4A/ne15THKuzwAW9vIZwBm9vATYOEabOTnXCeY77/ZtH/ujgUcBlwDLhurn6/4da74zun89gtKsaa3dPvR0EbDlip0/Bz7ZWtvUWrse+B7whCQPZhBql7fBT/zHgOdtyzFPV2ttdWvtnv70CuCh460/l+cK48533u1bgNbat1tr353s+vN4vjO6fw0ozaok70xyI3AscFqvfghw49BqP+p1D+nlkfVzzasY/Ae5xQFJvpXk0iRP7XXzZa7w2/Od7/t2NPN9/w6b0f27cEaHJo2Q5CvAfqMsWtFa+3xrbQWwIskpwOuB0xmcAhipjVNfwkRz7eusAO4BPtGX3QQ8vLX2syRLgc8l+X2KzxWmPd85uW9hcvMdxbzev6M1G6Vu2vvXgNL9qrV25CRXvQC4iEFA/Qh42NCyhwI/7vUPHaW+hInmmuQVwJ8BT++nOWitbQI29fK6JN8HDqL4XGF682WO7luY0s/ycJt5u3/HMKP711N8mjVJDhx6+lzgO738BeCYJDslOQA4EFjTWrsJ+GWSJ/YrgF4OjPWfXClJngWcDDy3tXbXUP0+SRb08u8ymOsP5vJcYez5Mg/37Xjm6/4dx8zu39m+GsTH9vsAPg1sBDYAFwIPGVq2gsEVQN9l6GofYFlv833gA/S7oVR/MHix+EZgfX9suWLxhcC1DK58ugp4zlyf63jznY/7to/9+QyOEjYBPwG+PM/376jznen9662OJEkleYpPklSSASVJKsmAkiSVZEBJkkoyoCRJJRlQ0hyS5I77oc+PJjm6lz+S5JCpbivJcUk+0MuvTfLyXr5k5N2uZ2jM90u/qsU7SUj6tdbaq2egj3MmXkuamEdQ0hyX5LAkVwx99tJevf7xve7yDD6faeMk+rrPkUmSxb2PZ/c7I3w6yZX98eRR+liZ5KShqhclWZPkui03S02yc5JV/fOBvpXkiAnqd0nyyT6ffwN2mf53THOFASXNfR8DTm6tHQpcw+B+hgCrgNe21g4HNk+n4yT7MrhH4mmttYuA9wLvaa09nsFdEj4yiW4WttaeALxpaGwnALTWHgu8BDg/yc7j1B8P3NXn+E5g6XTmo7nFU3zSHJZkD2DP1tqlvep84FMZfILt7q21/+r1FzC4cetU7AhcDJww1P+RwCFDH4b6gCS7T9DPZ/rXdQw+wA/gKcD7AVpr30lyA4ObqI5V/zTgfb1+Q5INU5yL5iADSpqfxvw47SSrgD8EftxaO2qcPu5hECrPBLYE1A7A4a21X43oc7yxbOpfN/ObvzljNRivI+/Ltp3xFJ80h7XWfgHcNvRBeH8JXNpau41+9+hef8xQm1e21g6bIJxgEAivAg5O8pZet5rB53YBg9e/pjn0yxh8SCVJDgIezuDmopOpfwxw6DS3qznEIyhpbtk1yfAnk/4T8ArgnCS7Aj8AXtmX/RVwXpI7gUuAX0x1Y621zUmOAS5McjtwIvDBfoptIYPgeO005nF2H/M1DI7UjmutbUoyVv2HgFV9u+uBNdPYpuYY72YuzVNJdmut3dHLbwEe3Fp74ywPS5o0j6Ck+evZSU5h8Ht+A3Dc7A5HmhqPoCRJJXmRhCSpJANKklSSASVJKsmAkiSVZEBJkkr6f1jYDgSAZ3c+AAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQ8klEQVR4nO3de7BeVX3G8e9DQrkE5GIQxFuoBZEqpSZa8TahUrVYrResKFbRWiuiaGeoiGEg9VKlMLXeEMGZiFVq63hFxpqKAo4FQ4IhhFGpihRHUHBQBGwKYfWPd0VfD+eeQ87vnHw/M++c9a6919prnX3Oec7e7373m9YakiRVs8NsD0CSpNEYUJKkkgwoSVJJBpQkqSQDSpJUkgElSSrJgJIklWRASZJKMqAkTUmSPZO8blu31fbHgJK2M0lenOTcJG+eZhd7AtMNma1pq+2MASVtR5KcCDwKWAW8YJrdvBt4ZJL1Sc7s/b4syZpe9+EkC5IsSnJRkquTbEzy4tHaSmOJ9+KTtg9JdgX+BzgQWAjs1Vq7bhr9LAG+2Fp7TH/+aOAfgRe01u5OcjZwBXAn8KzW2l/39fYA9hpuK41n4WwPQNI2cwRwQ2vttv78luGFSb4C7DdKuxWttc+P0+/TgaXAlUkAdgF+ClwAnJXkDAah9PUke23lHLQdMaCk7cfTgNVjLWytHTnNfgOc31o75T4LkqXAUcC7kqwGPjbNbWg75GtQ0vbjScA3k+wDkOQBSQ6ZRj+/BHYfen4xcHSSB/V+907yiCT7A3e11j4OnAU8bpS20pg8gpK2A0mWA3cDK4BHJfkGcCFwzlT7aq39LMk3kmwEvtRa+7skpwKrk+zQt3MCsAdwZpJ7e93xo7WdiflpfvIiCUlSSZ7ikySVZEBJkkoyoCRJJRlQkqSSvIpPACxevLgtWbJktochaR5Zt27dra21fabb3oASAEuWLGHt2rWzPQxJ80iSG7amvaf4JEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgNG0rV66c7SFI2oa29e98WmvbdIOqadmyZW3t2rVTapMEf3400n5n7cdP7vzJfer3XbQvN5908yyMSDNlqr/zSda11pZNd3sTHkElWZJk4xQGtDLJSdMd0FA/b0qy61TXS/LWrd32JLZ5cJL1Sb6V5JFJXjqDff8wyeIZ6OeOmRiPNFWjhdN49dJYFs72AMbxJuDjwF1TXO+twD9MZUNJFrTWNk+hyfOAz7fWTk+yHHgpcMH9uL2yli9fPttDUDVHjL3InxdNxWRfg1qQ5Lwk1yZZnWSXfuTwH0nWJfl6koNHNkpySZL3JLksybeTPD7JZ5L8d5J39HUWJbkoydVJNiZ5cZITgf2BryX5Wl/vGUkuT3JVkk8l2W3kekneDezSj24+0du9LMmaXvfhJAt6/R1J3pbkm8Dho006yWlJruzjOjcDRzEIxVf3sb0beGrv/2+TLEhyZm+3Icnf9L6W9zFeAFwz2ryHNv2GPs9rtnxfk+yd5HO9zyuSHNrrd0uyqq+7IckLR8xhcf++PXuU+b0mydoka2+55ZZJ/ihI0jbSWhv3ASwB7gEO68//HXgZcDFwYK/7I+CrvbwSOKmXLwHO6OU3Aj8GHgzsBPwIeCDwQuC8oe3t0b/+EFjcy4uBy4BF/fnJwGkj1+vP7xgqPxq4ENixPz8beHkvN+AvJpj73kPlfwGeM8oclwNfHFrvNcCpvbwTsBY4oK93J3BAXzbevN/Qy68DPtLL7wdO7+U/Btb38hnAPw/1s9eW7wOwL/BN4E8m2s9Lly5tUzX48ZF+GysZ86G5baq/88DaNsHfnvEekz3Fd31rbX0vr2MQWk8CPpVkyzo7jdH2C/3rNcC1rbWbAJL8AHhYrz8ryRn9D/3XR+njicAhwDf69n4HuHwS4346sBS4srfbBfhpX7YZ+PQE7Y9I8mZgV2Bv4FoGgTeeZwCHJjm6P98DOBD4P2BNa+36Xj/evD/Tv64DXtDLT2EQarTWvprkgUn2AI4EjtnSsLV2Wy/uyOCfiBNaa5dOMGZJKmeyAbVpqLyZwX/mP2+tHTaFtveO6OdeYGFr7bokS4GjgHclWd1ae9uIPgL8Z2vtJZMc73C781trp4yy7H/bOK8DJdmZwRHXstbajUlWAjtPcptvaK19eUR/yxkcQQEwwby3fJ8285t99Ov/BIa0Xj/aZTX3MAi4ZwIGlLaZfRftO+ZVfNJUTPd9ULcD1yd5EUB/beYPptNRkv2Bu1prHwfOAh7XF/0S2L2XrwCenOT3eptdkxw0ynoAdyfZsZcvBo5O8qDebu8kj5jk0LaE0a1JdgOOHmO9kdv/MnD8ljEkOSjJopGNxpn3WC4Dju1tlwO3ttZuB1YDrx/qd69ebMCrgIOTvGWCvqfl9NNPvz+61Rx380k3005v93l4ifnct61/57fmKr5jgQ8lOZXB6aRPAldPo5/HAmcmuRe4Gzi+158LfCnJTa21I5IcB/xrki2nEk8Frhu5Xn++IclVrbVj+/hWJ9mh938CcMNEg2qt/TzJeQxOxf0QuHKMVTcA9yS5Gvgo8F4Gp0CvyuC84i0Mrvqb7LzHshJYlWQDgysWX9Hr3wF8MIO3AmwG/p5+irC1tjnJMcCFSW5vrZ09wTamxDfqStsX36irWTGdN+pK0nhyf79RV5Kk2VD5jbrbTJLPMrgUfNjJIy90kCRtOwYU0Fp7/myPQZL02zzFJ0kqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklGVCSpJIMKElSSQaUJKkkA0qSVJIBJUkqyYCSJJVkQEmSSjKgJEklpbU222NQAUluAW6Y7XEAi4FbZ3sQ25Dznd+29/k+orW2z3Q7M6BUSpK1rbVlsz2ObcX5zm/Od+t4ik+SVJIBJUkqyYBSNefO9gC2Mec7vznfreBrUJKkkjyCkiSVZEBJkkoyoDRrkrw9yYYk65OsTrL/0LJTknwvyXeTPHOofmmSa/qy9yXJ7Ix+apKcmeQ7fb6fTbJnr1+S5Ff9e7A+yTlDbebkXGHs+fZl82rfAiR5UZJrk9ybZNlQ/Xzdv6POty+buf3bWvPhY1YewAOGyicC5/TyIcDVwE7AAcD3gQV92RrgcCDAl4A/ne15THKuzwAW9vIZwBm9vATYOEabOTnXCeY77/ZtH/ujgUcBlwDLhurn6/4da74zun89gtKsaa3dPvR0EbDlip0/Bz7ZWtvUWrse+B7whCQPZhBql7fBT/zHgOdtyzFPV2ttdWvtnv70CuCh460/l+cK48533u1bgNbat1tr353s+vN4vjO6fw0ozaok70xyI3AscFqvfghw49BqP+p1D+nlkfVzzasY/Ae5xQFJvpXk0iRP7XXzZa7w2/Od7/t2NPN9/w6b0f27cEaHJo2Q5CvAfqMsWtFa+3xrbQWwIskpwOuB0xmcAhipjVNfwkRz7eusAO4BPtGX3QQ8vLX2syRLgc8l+X2KzxWmPd85uW9hcvMdxbzev6M1G6Vu2vvXgNL9qrV25CRXvQC4iEFA/Qh42NCyhwI/7vUPHaW+hInmmuQVwJ8BT++nOWitbQI29fK6JN8HDqL4XGF682WO7luY0s/ycJt5u3/HMKP711N8mjVJDhx6+lzgO738BeCYJDslOQA4EFjTWrsJ+GWSJ/YrgF4OjPWfXClJngWcDDy3tXbXUP0+SRb08u8ymOsP5vJcYez5Mg/37Xjm6/4dx8zu39m+GsTH9vsAPg1sBDYAFwIPGVq2gsEVQN9l6GofYFlv833gA/S7oVR/MHix+EZgfX9suWLxhcC1DK58ugp4zlyf63jznY/7to/9+QyOEjYBPwG+PM/376jznen9662OJEkleYpPklSSASVJKsmAkiSVZEBJkkoyoCRJJRlQ0hyS5I77oc+PJjm6lz+S5JCpbivJcUk+0MuvTfLyXr5k5N2uZ2jM90u/qsU7SUj6tdbaq2egj3MmXkuamEdQ0hyX5LAkVwx99tJevf7xve7yDD6faeMk+rrPkUmSxb2PZ/c7I3w6yZX98eRR+liZ5KShqhclWZPkui03S02yc5JV/fOBvpXkiAnqd0nyyT6ffwN2mf53THOFASXNfR8DTm6tHQpcw+B+hgCrgNe21g4HNk+n4yT7MrhH4mmttYuA9wLvaa09nsFdEj4yiW4WttaeALxpaGwnALTWHgu8BDg/yc7j1B8P3NXn+E5g6XTmo7nFU3zSHJZkD2DP1tqlvep84FMZfILt7q21/+r1FzC4cetU7AhcDJww1P+RwCFDH4b6gCS7T9DPZ/rXdQw+wA/gKcD7AVpr30lyA4ObqI5V/zTgfb1+Q5INU5yL5iADSpqfxvw47SSrgD8EftxaO2qcPu5hECrPBLYE1A7A4a21X43oc7yxbOpfN/ObvzljNRivI+/Ltp3xFJ80h7XWfgHcNvRBeH8JXNpau41+9+hef8xQm1e21g6bIJxgEAivAg5O8pZet5rB53YBg9e/pjn0yxh8SCVJDgIezuDmopOpfwxw6DS3qznEIyhpbtk1yfAnk/4T8ArgnCS7Aj8AXtmX/RVwXpI7gUuAX0x1Y621zUmOAS5McjtwIvDBfoptIYPgeO005nF2H/M1DI7UjmutbUoyVv2HgFV9u+uBNdPYpuYY72YuzVNJdmut3dHLbwEe3Fp74ywPS5o0j6Ck+evZSU5h8Ht+A3Dc7A5HmhqPoCRJJXmRhCSpJANKklSSASVJKsmAkiSVZEBJkkr6f1jYDgSAZ3c+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -222,7 +252,7 @@ }, { "cell_type": "markdown", - "id": "34cf0e89", + "id": "missing-frequency", "metadata": { "tags": [] }, @@ -248,7 +278,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUIklEQVR4nO3de5BmdX3n8feHGZZrgsKg0SgMGBVRLnHGRCNmRyGaELPRVcQ4BGJKSakhMRVvMJYMVtzAYm12s4qCKdHIgLUaL4k3UAwXg1xmcG54Dxe1JBGSIMIkLMx894/z6+Wx7Znu6enp/rX9flU91ef5nfP7ne95+vJ5znlOn5OqQpKk3uwx1wVIkjQRA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSlqAkpyU5IYkG5N8O8nZSf4gybt3cdxHJHntXPXXTxcDSlpgkpwGvBl4SVUdDRwLbAGOBjbt4vCPAHYlYHa1v36KGFDSApLkZ4H/Abysqr4HUFX3VdX5wFHsekCdCzwhyfok57d1npLkxtZ2YZJFSfZL8ukkG5JsTnLy9vpr4Vo81wVImlUvBm6oqlsnmPc0YPMujv8W4GlVdSxAkqcAJwPPrqoHk1wArATuB75fVb/Zljtgov5a2AwoaWF5KrB+fGOSxwM/qqofbq9jki8APzfBrFVV9cntdDseWAbclARgH+AHwKXAO5OcB3yqqq7dmY3QwmBASQvL/QwhMd6knz9V1QnTWF+AD1bVmT8xI1kGnAj8eZIrqurt0xhfP8X8DEpaWD4DnJTk0QBJ9kryambm8yeAHwE/M/L8SuClSR7V1ndgkkOTPBbYUlWXAO8Enr6d/lrA4v2gpIUlye8CfwosYjiKcglwJPDrDAEBcGdVPWua41/KsEf22ap6YzsB4kyGN8QPAq8DDgDOB7a1ttdU1dqJ+k9rI/VTwYCSJHXJQ3ySpC4ZUJKkLhlQkqQuGVCSpC75f1ACYMmSJbV06dK5LkPSArRu3bq7q+rg8e0GlABYunQpa9eunesyJC1ASe6YqN1DfJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVDq0urVq+e6BEnTMJO/uwaUdt6aNbB0Keyxx/B1zZoZX8U555wz42NK2v1m8nd30oBKsjTJ5qkOmGR1kjfsWlmQ5PVJ9t3Z5ZKctavrnsI6j0iyPslXkjwhyStmcOzbkyyZgXHum4l6fsKaNXD66XDHHVA1fD399N0SUpIWtsVzXcAOvB64BNiyk8udBfy3nVlRkkVVtXUnurwI+GRVnZ1kBfAK4NLduL5+rFoFW8Z9S7ZsGdpXrpzRVa1YsWJGx5M0v0z1EN+iJO9LckuSK5Ls0/YcPpdkXZJrkxwxvlOSq5L8RZJrknwtyTOSfCzJt5L8WVtmvySfTrIhyeYkJyf5I+CxwN8n+fu23POTfDnJzUk+kmT/8cslORfYp+3drGn9TklyY2u7MMmi1n5fkrcnuQF41kQbneRtSW5qdV2UwYkMofiqVtu5wHPa+H+SZFGS81u/jUn+oI21otV4KbBpou0eWfUZbTs3jb2uSQ5M8ok25vVJjm7t+ye5uC27MclLxm3Dkva6/eYE23d6krVJ1t51111T+0n4znd2rl2SpquqdvgAlgIPAce25/8HOAW4Enhia/tl4IttejXwhjZ9FXBem/5j4PvAY4C9gO8BBwEvAd43sr4D2tfbgSVteglwDbBfe/5m4G3jl2vP7xuZfgrwd8Ce7fkFwKltuoCXTbLtB45Mfwj4rQm2cQXwqZHlTgfe2qb3AtYCh7Xl7gcOa/N2tN1ntOnXAn/Vpv83cHabfh6wvk2fB/zPkXEeOfY6AI8GbgB+bbLv87Jly2pKDj20aji49+OPQw+dWv8pGn40Jc030/ndBdbWBH+XproHdVtVrW/T6xhC61eAjyRZD1zIEDwT+dv2dRNwS1XdWVUPALcCj2/tJyQ5L8lzquqHE4zxTOBI4B/a+k4DDp1C3ccDy4CbWr/jgcPbvK3A30zS/7lJbkiyiSEUnjqFdT4fOLWt7waGEH5im3djVd3Wpne03R9rX8dea4DjGEKSqvoicFCSA4ATgHePdayqf2uTezK8iXhTVX1+CnVPzTveAfuO+2hw332HdkmaQVP9DOqBkemtDO/M76mqY3ei77Zx42wDFlfVN5MsA04E/jzJFVX19nFjBPh8Vf3OFOsd7ffBqjpzgnn/UTv4HCjJ3gx7XMur6rtJVgN7T3GdZ1TV5ePGW8GwBwXAJNs99jpt5eHvUSZYV7X2mmDeQwwB9wLg6inUPTVjnzOtWjUc1jvkkCGcZvjzJ0ma7mnm9wK3JTkJoH02c8x0BkryWGBLVV0CvBN4epv1I+Bn2vT1wLOT/ELrs2+SJ02wHMCDSfZs01cCL03yqNbvwCRT2fOCh8Po7iT7Ay/dznLj13858JqxGpI8Kcl+4zvtYLu35xpgZeu7Ari7qu4FrgD+cGTcR7bJAn4fOCLJWyYZe+esXAm33w7btg1fd0M4nX322TM+pqTdbyZ/d3flLL6VwHuSvJXhcNKHgQ3TGOco4Pwk24AHgde09ouAzya5s6qem+T3gMuS7NXmvxX45vjl2vONSW6uqpWtviuS7NHGfx1wx2RFVdU9Sd7HcCjuduCm7Sy6EXgoyQbgA8D/Yjgsd3OSAHcxnPU31e3entXAxUk2MpyxeFpr/zPg3Rn+FWArcA7tEGFVbU3ycuDvktxbVRdMso5u+I+60vw0k7+7GT6f0kK3fPnyWrt27VyXIWkBSrKuqpaPb/dKEpKkLvX8j7qzJsnHGU4FH/Xm8Sc6SJJmjwEFVNWL57oGSdKP8xCfJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLqaq5rkEdSHIXcMccrHoJcPccrHdXWffsm6+1W/fkDq2qg8c3GlCaU0nWVtXyua5jZ1n37JuvtVv39HmIT5LUJQNKktQlA0pz7aK5LmCarHv2zdfarXua/AxKktQl96AkSV0yoCRJXTKgNGuSnJ/k60k2Jvl4kkeMzDszybeTfCPJC0balyXZ1Ob9ZZLMQd0nJbklybYky0fa/1OSi1t9G5KsmCd175nkg62+ryU5c57UvTLJ+pHHtiTH9l53m3d0ki+3+ZuS7N173UmWJvn3kdf7vSPzZqfuqvLhY1YewPOBxW36POC8Nn0ksAHYCzgM+EdgUZt3I/AsIMBngd+Yg7qfAjwZuApYPtL+OuDiNv0oYB2wxzyo+xXAh9v0vsDtwNLe6x63zFHArSPPu60bWAxsBI5pzw+aJz/fS4HN2+kzK3W7B6VZU1VXVNVD7en1wOPa9G8z/MF8oKpuA74N/FKSxwA/W1VfruG34q+BF81B3V+rqm9MMOtI4Mq2zA+Ae4Dl86DuAvZLshjYB/i/wL3zoO5RvwNcBjAP6n4+sLGqNrTl/qWqts6Duic0m3UbUJorv8/wzgvg54Hvjsz7Xmv7+TY9vr0XG4DfTrI4yWHAMuDx9F/3R4H7gTuB7wDvrKp/pf+6R51MCyj6r/tJQCW5PMnNSd7U2nuvG+CwJF9JcnWS57S2Wat78e4YVAtXki8APzfBrFVV9cm2zCrgIWDNWLcJlq8dtM+4qdQ9gfczHB5Zy3Adw+sYtqv3un8J2Ao8FngkcG0bp/e6x/r+MrClqjaPNU2wWE91LwaOA54BbAGuTLIOuHeCZXuq+07gkKr6lyTLgE8keSqz+HobUJpRVXXCjuYnOQ14IXB8OzwAwzuwx48s9jjg+639cRO0z7jJ6t5On4eAPxl7nuQ64FvAv9Fx3QyfQX2uqh4EfpDkH4DlwLX0XfeYl/Pw3hN0/nPCUN/VVXU3QJLPAE8HLqHjuqvqAeCBNr0uyT8y7A3O2uvtIT7NmiS/DrwZ+C9VtWVk1t8CL0+yVztU9kTgxqq6E/hRkme2s4ROBXb47no2Jdk3yX5t+teAh6rqq73XzXBY73kZ7Ac8E/j6PKibJHsAJwEfHmubB3VfDhzdfl4WA/8Z6P7nJMnBSRa16cMZfi9vndW6Z/uMER8L98Fw8sN3gfXt8d6ReasYzt77BiNnBDG8s9/c5r2LdvWTWa77xQzvGh8A/hm4vLUvbfV+DfgCwy0D5kPd+wMfAW4Bvgq8cT7U3eatAK6foE/vdZ/SXu/NwH+fD3UDL2k1bwBuBn5rtuv2UkeSpC55iE+S1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNK6lySDyR5aZv+qyRHtumzxi133e5c91Tap7mOq8auop3kM0ke0a6kvXmyviNjrE7yhjb99iQntOnbkyyZiTrHrW+3jKsfZ0BJ80hVvaqqvtqenjVu3q/MQUkzqqpOrKp7dnGMt1XVF2aoJM0hA0qaQUlOzXC/qw1JPtTaDk1yZWu/Mskhrf0D7V461yW5dWQvKUneleSrST7NcCuPsfGvSrI8ybnAPu0+PWvavPtG+p+fZHOGe/ac3NpXtP4fzXBfrjXtSgAkeVuSm1qfi8bap7jNx2e4oOimJO9PsldrP7Gt50ttOz81hbF+Ys8kyeFt/GckeUKSzyVZl+TaJEdMMMb4vbszMlykddPY8kkOTPKJ9j25PsnRk7QflOSKVseFTHw9Os0wA0qaIRkupLkKeF5VHQP8cZv1LuCvq+pohgvk/uVIt8cwXEj0hcC5re3FDPfnOQp4NfATe0ZV9Rbg36vq2KpaOW72fwWOBY4BTgDOz3CLBIBfBF7PcKuQw4Fnj9VYVc+oqqcx3ILjhVPc5r2BDwAnV9VRDNf3fE1rv5DhqiDHAQdPZbwJxn8y8DfAK6vqJuAi4IyqWga8AbhgCsPcXVVPB97T+gCcA3ylfU/OYrhlxI7azwa+VFW/yHBprkOmsz3aOQaUNHOeB3y02kVBa7iFBQw3dru0TX+IIZDGfKKqtrXDdo9ubb8KXFZVW6vq+8AXd7KO40b6/zNwNcOVtGG4xuH3qmobw+Wmlrb25ya5Icmmth1PneK6ngzcVlXfbM8/2Oo/guG6bbe19ssm6jyJgxmu8XZKVa1Psj9DWH8kyXqGAHzMDvqP+Vj7uo6Ht/c4hu8FVfVF4KAkB+yg/VcZLu5KVX2a4YLA2s28mrk0c8LUbjswuswD4/pPtMx06tie0fVtBRa3vZ0LGO6m+t0kq4G9d3Fd260hyeUMYby2ql61g7F/yHDtxmczXBNuD+Ceqjp2irWNGdvmrTz8N286t3jxunCzzD0oaeZcCbwsyUEwfJ7R2q9juEUEwErgS5OMcw3D1d0XtUNzz93Ocg8m2XM7/U9u/Q9mePd/4w7WNxZGd7e9lJ05O+/rwNIkv9Ce/y7DHtvXgcOTLG3tJ491qKoXtEOTOwonGO70+yLg1CSvqKp7gduSnAT//7O2Y3ai1lHXMHwvSLKC4TDgvVNs/w2G+2hpN3MPSpohVXVLkncAVyfZCnwF+D3gj4D3J3kjcBfwykmG+jjDYbZNwDcZ/uBP5CJgY5Kbx30O9XGGw4obGN71v6mq/mmiEwpa3fckeV9b3+3ATZNt60jf/0jySobDbotb3/dW1QNJXgt8Lsnd7DggdzT+/UleCHw+yf0MIfGeJG8F9mS47caGaQy9Grg4yUaGmwieNkn7OcBlSW5m+H58Zzrbo53j1cwl7RZJ9q+q+9oZge8GvlVVfzHXdWn+8BCfpN3l1e1khluAAxhOapCmzD0oSVKX3IOSJHXJgJIkdcmAkiR1yYCSJHXJgJIkden/ASnV5vr0NT7jAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUIklEQVR4nO3de5BmdX3n8feHGZZrgsKg0SgMGBVRLnHGRCNmRyGaELPRVcQ4BGJKSakhMRVvMJYMVtzAYm12s4qCKdHIgLUaL4k3UAwXg1xmcG54Dxe1JBGSIMIkLMx894/z6+Wx7Znu6enp/rX9flU91ef5nfP7ne95+vJ5znlOn5OqQpKk3uwx1wVIkjQRA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSlqAkpyU5IYkG5N8O8nZSf4gybt3cdxHJHntXPXXTxcDSlpgkpwGvBl4SVUdDRwLbAGOBjbt4vCPAHYlYHa1v36KGFDSApLkZ4H/Abysqr4HUFX3VdX5wFHsekCdCzwhyfok57d1npLkxtZ2YZJFSfZL8ukkG5JsTnLy9vpr4Vo81wVImlUvBm6oqlsnmPc0YPMujv8W4GlVdSxAkqcAJwPPrqoHk1wArATuB75fVb/Zljtgov5a2AwoaWF5KrB+fGOSxwM/qqofbq9jki8APzfBrFVV9cntdDseWAbclARgH+AHwKXAO5OcB3yqqq7dmY3QwmBASQvL/QwhMd6knz9V1QnTWF+AD1bVmT8xI1kGnAj8eZIrqurt0xhfP8X8DEpaWD4DnJTk0QBJ9kryambm8yeAHwE/M/L8SuClSR7V1ndgkkOTPBbYUlWXAO8Enr6d/lrA4v2gpIUlye8CfwosYjiKcglwJPDrDAEBcGdVPWua41/KsEf22ap6YzsB4kyGN8QPAq8DDgDOB7a1ttdU1dqJ+k9rI/VTwYCSJHXJQ3ySpC4ZUJKkLhlQkqQuGVCSpC75f1ACYMmSJbV06dK5LkPSArRu3bq7q+rg8e0GlABYunQpa9eunesyJC1ASe6YqN1DfJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVDq0urVq+e6BEnTMJO/uwaUdt6aNbB0Keyxx/B1zZoZX8U555wz42NK2v1m8nd30oBKsjTJ5qkOmGR1kjfsWlmQ5PVJ9t3Z5ZKctavrnsI6j0iyPslXkjwhyStmcOzbkyyZgXHum4l6fsKaNXD66XDHHVA1fD399N0SUpIWtsVzXcAOvB64BNiyk8udBfy3nVlRkkVVtXUnurwI+GRVnZ1kBfAK4NLduL5+rFoFW8Z9S7ZsGdpXrpzRVa1YsWJGx5M0v0z1EN+iJO9LckuSK5Ls0/YcPpdkXZJrkxwxvlOSq5L8RZJrknwtyTOSfCzJt5L8WVtmvySfTrIhyeYkJyf5I+CxwN8n+fu23POTfDnJzUk+kmT/8cslORfYp+3drGn9TklyY2u7MMmi1n5fkrcnuQF41kQbneRtSW5qdV2UwYkMofiqVtu5wHPa+H+SZFGS81u/jUn+oI21otV4KbBpou0eWfUZbTs3jb2uSQ5M8ok25vVJjm7t+ye5uC27MclLxm3Dkva6/eYE23d6krVJ1t51111T+0n4znd2rl2SpquqdvgAlgIPAce25/8HOAW4Enhia/tl4IttejXwhjZ9FXBem/5j4PvAY4C9gO8BBwEvAd43sr4D2tfbgSVteglwDbBfe/5m4G3jl2vP7xuZfgrwd8Ce7fkFwKltuoCXTbLtB45Mfwj4rQm2cQXwqZHlTgfe2qb3AtYCh7Xl7gcOa/N2tN1ntOnXAn/Vpv83cHabfh6wvk2fB/zPkXEeOfY6AI8GbgB+bbLv87Jly2pKDj20aji49+OPQw+dWv8pGn40Jc030/ndBdbWBH+XproHdVtVrW/T6xhC61eAjyRZD1zIEDwT+dv2dRNwS1XdWVUPALcCj2/tJyQ5L8lzquqHE4zxTOBI4B/a+k4DDp1C3ccDy4CbWr/jgcPbvK3A30zS/7lJbkiyiSEUnjqFdT4fOLWt7waGEH5im3djVd3Wpne03R9rX8dea4DjGEKSqvoicFCSA4ATgHePdayqf2uTezK8iXhTVX1+CnVPzTveAfuO+2hw332HdkmaQVP9DOqBkemtDO/M76mqY3ei77Zx42wDFlfVN5MsA04E/jzJFVX19nFjBPh8Vf3OFOsd7ffBqjpzgnn/UTv4HCjJ3gx7XMur6rtJVgN7T3GdZ1TV5ePGW8GwBwXAJNs99jpt5eHvUSZYV7X2mmDeQwwB9wLg6inUPTVjnzOtWjUc1jvkkCGcZvjzJ0ma7mnm9wK3JTkJoH02c8x0BkryWGBLVV0CvBN4epv1I+Bn2vT1wLOT/ELrs2+SJ02wHMCDSfZs01cCL03yqNbvwCRT2fOCh8Po7iT7Ay/dznLj13858JqxGpI8Kcl+4zvtYLu35xpgZeu7Ari7qu4FrgD+cGTcR7bJAn4fOCLJWyYZe+esXAm33w7btg1fd0M4nX322TM+pqTdbyZ/d3flLL6VwHuSvJXhcNKHgQ3TGOco4Pwk24AHgde09ouAzya5s6qem+T3gMuS7NXmvxX45vjl2vONSW6uqpWtviuS7NHGfx1wx2RFVdU9Sd7HcCjuduCm7Sy6EXgoyQbgA8D/Yjgsd3OSAHcxnPU31e3entXAxUk2MpyxeFpr/zPg3Rn+FWArcA7tEGFVbU3ycuDvktxbVRdMso5u+I+60vw0k7+7GT6f0kK3fPnyWrt27VyXIWkBSrKuqpaPb/dKEpKkLvX8j7qzJsnHGU4FH/Xm8Sc6SJJmjwEFVNWL57oGSdKP8xCfJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLqaq5rkEdSHIXcMccrHoJcPccrHdXWffsm6+1W/fkDq2qg8c3GlCaU0nWVtXyua5jZ1n37JuvtVv39HmIT5LUJQNKktQlA0pz7aK5LmCarHv2zdfarXua/AxKktQl96AkSV0yoCRJXTKgNGuSnJ/k60k2Jvl4kkeMzDszybeTfCPJC0balyXZ1Ob9ZZLMQd0nJbklybYky0fa/1OSi1t9G5KsmCd175nkg62+ryU5c57UvTLJ+pHHtiTH9l53m3d0ki+3+ZuS7N173UmWJvn3kdf7vSPzZqfuqvLhY1YewPOBxW36POC8Nn0ksAHYCzgM+EdgUZt3I/AsIMBngd+Yg7qfAjwZuApYPtL+OuDiNv0oYB2wxzyo+xXAh9v0vsDtwNLe6x63zFHArSPPu60bWAxsBI5pzw+aJz/fS4HN2+kzK3W7B6VZU1VXVNVD7en1wOPa9G8z/MF8oKpuA74N/FKSxwA/W1VfruG34q+BF81B3V+rqm9MMOtI4Mq2zA+Ae4Dl86DuAvZLshjYB/i/wL3zoO5RvwNcBjAP6n4+sLGqNrTl/qWqts6Duic0m3UbUJorv8/wzgvg54Hvjsz7Xmv7+TY9vr0XG4DfTrI4yWHAMuDx9F/3R4H7gTuB7wDvrKp/pf+6R51MCyj6r/tJQCW5PMnNSd7U2nuvG+CwJF9JcnWS57S2Wat78e4YVAtXki8APzfBrFVV9cm2zCrgIWDNWLcJlq8dtM+4qdQ9gfczHB5Zy3Adw+sYtqv3un8J2Ao8FngkcG0bp/e6x/r+MrClqjaPNU2wWE91LwaOA54BbAGuTLIOuHeCZXuq+07gkKr6lyTLgE8keSqz+HobUJpRVXXCjuYnOQ14IXB8OzwAwzuwx48s9jjg+639cRO0z7jJ6t5On4eAPxl7nuQ64FvAv9Fx3QyfQX2uqh4EfpDkH4DlwLX0XfeYl/Pw3hN0/nPCUN/VVXU3QJLPAE8HLqHjuqvqAeCBNr0uyT8y7A3O2uvtIT7NmiS/DrwZ+C9VtWVk1t8CL0+yVztU9kTgxqq6E/hRkme2s4ROBXb47no2Jdk3yX5t+teAh6rqq73XzXBY73kZ7Ac8E/j6PKibJHsAJwEfHmubB3VfDhzdfl4WA/8Z6P7nJMnBSRa16cMZfi9vndW6Z/uMER8L98Fw8sN3gfXt8d6ReasYzt77BiNnBDG8s9/c5r2LdvWTWa77xQzvGh8A/hm4vLUvbfV+DfgCwy0D5kPd+wMfAW4Bvgq8cT7U3eatAK6foE/vdZ/SXu/NwH+fD3UDL2k1bwBuBn5rtuv2UkeSpC55iE+S1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNK6lySDyR5aZv+qyRHtumzxi133e5c91Tap7mOq8auop3kM0ke0a6kvXmyviNjrE7yhjb99iQntOnbkyyZiTrHrW+3jKsfZ0BJ80hVvaqqvtqenjVu3q/MQUkzqqpOrKp7dnGMt1XVF2aoJM0hA0qaQUlOzXC/qw1JPtTaDk1yZWu/Mskhrf0D7V461yW5dWQvKUneleSrST7NcCuPsfGvSrI8ybnAPu0+PWvavPtG+p+fZHOGe/ac3NpXtP4fzXBfrjXtSgAkeVuSm1qfi8bap7jNx2e4oOimJO9PsldrP7Gt50ttOz81hbF+Ys8kyeFt/GckeUKSzyVZl+TaJEdMMMb4vbszMlykddPY8kkOTPKJ9j25PsnRk7QflOSKVseFTHw9Os0wA0qaIRkupLkKeF5VHQP8cZv1LuCvq+pohgvk/uVIt8cwXEj0hcC5re3FDPfnOQp4NfATe0ZV9Rbg36vq2KpaOW72fwWOBY4BTgDOz3CLBIBfBF7PcKuQw4Fnj9VYVc+oqqcx3ILjhVPc5r2BDwAnV9VRDNf3fE1rv5DhqiDHAQdPZbwJxn8y8DfAK6vqJuAi4IyqWga8AbhgCsPcXVVPB97T+gCcA3ylfU/OYrhlxI7azwa+VFW/yHBprkOmsz3aOQaUNHOeB3y02kVBa7iFBQw3dru0TX+IIZDGfKKqtrXDdo9ubb8KXFZVW6vq+8AXd7KO40b6/zNwNcOVtGG4xuH3qmobw+Wmlrb25ya5Icmmth1PneK6ngzcVlXfbM8/2Oo/guG6bbe19ssm6jyJgxmu8XZKVa1Psj9DWH8kyXqGAHzMDvqP+Vj7uo6Ht/c4hu8FVfVF4KAkB+yg/VcZLu5KVX2a4YLA2s28mrk0c8LUbjswuswD4/pPtMx06tie0fVtBRa3vZ0LGO6m+t0kq4G9d3Fd260hyeUMYby2ql61g7F/yHDtxmczXBNuD+Ceqjp2irWNGdvmrTz8N286t3jxunCzzD0oaeZcCbwsyUEwfJ7R2q9juEUEwErgS5OMcw3D1d0XtUNzz93Ocg8m2XM7/U9u/Q9mePd/4w7WNxZGd7e9lJ05O+/rwNIkv9Ce/y7DHtvXgcOTLG3tJ491qKoXtEOTOwonGO70+yLg1CSvqKp7gduSnAT//7O2Y3ai1lHXMHwvSLKC4TDgvVNs/w2G+2hpN3MPSpohVXVLkncAVyfZCnwF+D3gj4D3J3kjcBfwykmG+jjDYbZNwDcZ/uBP5CJgY5Kbx30O9XGGw4obGN71v6mq/mmiEwpa3fckeV9b3+3ATZNt60jf/0jySobDbotb3/dW1QNJXgt8Lsnd7DggdzT+/UleCHw+yf0MIfGeJG8F9mS47caGaQy9Grg4yUaGmwieNkn7OcBlSW5m+H58Zzrbo53j1cwl7RZJ9q+q+9oZge8GvlVVfzHXdWn+8BCfpN3l1e1khluAAxhOapCmzD0oSVKX3IOSJHXJgJIkdcmAkiR1yYCSJHXJgJIkden/ASnV5vr0NT7jAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -290,7 +320,7 @@ "\n", "Aim: The number or N-test is the most conceptually simple test of a forecast: To test whether the number of observed events is consistent with that of the forecast.\n", "\n", - "Method: The originial N-test was introduced by Schorlemmer et al (2007) and modified by Zechar et al (2010). The observed number of events is given by,\n", + "Method: The original N-test was introduced by Schorlemmer et al (2007) and modified by Zechar et al (2010). The observed number of events is given by,\n", "\n", "$$N_{obs} = \\sum_{m_i, s_j \\in R} \\omega(m_i, s_j).$$\n", "\n", @@ -338,7 +368,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAATPklEQVR4nO3debRdZX3G8e9DoBAGUYwiqDUOuLRajCVO1bZxqNahjsGhUKFa5wGstCqLkuhyVRC1tNY6gAIq1IXz0IG4BMS6EEggJCBVVwWtCgrLsjCCVMKvf+z34uF6x3CT+97r97PWWXef97z73e97dnKes4ezd6oKSZJ6s9N8d0CSpIkYUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVDSApZka5KNSS5L8skku09R95lJ3ryD+7c2yY1J7j5StmVcnae0MWxMsiXJt9r0R2exnMOT7D+Xfdf8M6Ckhe2mqlpRVQ8F/g945WQVq+oLVXXcjuvaba4D3jjZi1V1VhvDCmA9cEh7/uJZLONwwIBaZAwoafH4GvCAJPsk+VySTUm+keRAuG0r45/a9MFtq+vSJOe1sockubBtvWxKckAr/6tW97IkR7ay5UmuSHJSksuTrEuydJJ+fQR4QZJ9ZjOYJIeO9OeDSZa0x6mtL5uTvCHJamAlcHqrO1k/tMAYUNIikGRn4KnAZuCtwCVVdSBwNDDRrrJjgadU1cOAZ7ayVwL/0LZkVgI/SHIQ8BfAo4BHAy9L8vBW/wDgfVX1EOB64HmTdG8LQ0gdMYvxPBh4AfDY1p+twCHACuCeVfXQqvpd4JSq+hS33/K6aabLUd8MKGlhW5pkI8MH9PeBDwOPAz4GUFVnA3dNsve4+b4OnJrkZcCSVnY+cHSSNwH3aR/0jwM+W1U/r6otwGeAP2j1r6yqjW16A7B8in7+I3BYkjvNcFxPBA4CLmrjeyJwP+C7wP2SvDfJnwA3zLA9LUA7z3cHJN0hN7UtjNskyQT1bnfRzap6ZZJHAU8HNiZZUVVnJLmglZ2V5C+Bidoac/PI9FZg0l1rVXV9kjOAV4/08zXAy9rTp1XVj0aHAZxWVW8Z31aShwFPAV4DPB94yRR91ALmFpS0+JzHsDuMJKuA66rqdlsaSe5fVRdU1bEMJzHcO8n9gO9W1T8CXwAObG09O8nuSfYAnsNwrGtbvAd4Be2LcVW9b+zkiHHhBPAVYPXY2X/tuNp9kiwDdqqqTwN/C/xeq/8zYK9t7Jc65RaUtPisBU5Jsgm4EThsgjontJMgwhAGlwJvBg5N8kvgGuBtVfXTJKcCF7b5Tq6qS5Isn22nquq6JJ8F3jCDut9McgywLslOwC8ZtphuamMb+3I9toV1KvCBJDcBj/E41OIQb7chSeqRu/gkSV0yoCRJXTKgJEldMqAkSV3yLD4BsGzZslq+fPl8d0PSb6ANGzZcV1V3G19uQAmA5cuXs379+vnuhqTfQEm+N1G5u/gkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqC0zdauXTvfXZDUmbn8XEhVzVljWrhWrlxZ69evn7bePd51D3788x//Wvm+e+zLNUddsz26JmkBScJscyXJhqpaOb582i2oJMuTXDaLBa1NctSsejdxO0cm2X229ZIcfUeXPYNlPijJxiSXJLl/kj+bw7avSrJsDtrZMhf9GW+icJqqXJK21c7z3YEpHAl8HLhxlvWOBv5uNgtKsqSqts5ilmcDn6+qNUlWAX8GnLEdl7cgrFq1ar67IGkRmekxqCVJTkpyeZJ1SZa2LYf/SLIhydeSPGj8TEnOTfL3Sc5LckWSRyT5TJLvJHl7q7NHkn9NcmmSy5K8IMnrgf2Bc5Kc0+o9Ocn5SS5O8skke46vl+Q4YGnbujm9zXdokgtb2QeTLGnlW5K8LckFwGMmGnSSY5Nc1Pr1oQyexhCKf9n6dhzwB639NyRZkuSENt+mJK9oba1qfTwD2DzRuEcW/bo2zs1j72uSfZJ8rrX5jSQHtvI9k5zS6m5K8rxxY1jW3renTzC+lydZn2T9tddeO8N/CpK0g1TVlA9gOXALsKI9PxM4FPgKcEArexRwdpteCxzVps8Fjm/TRwA/AvYDdgV+ANwVeB5w0sjy9m5/rwKWtellwHnAHu35m4Bjx9drz7eMTD8Y+CKwS3v+z8CL23QBz59m7PuMTH8M+NMJxrgK+NJIvZcDx7TpXYH1wH1bvZ8D922vTTXu17XpVwMnt+n3Amva9BOAjW36eODEkXbuMvY+APsCFwB/PN16Puigg2omWMukD0kaYmXW86yvCT6XZrqL78qq2timNzCE1u8Dn0wyVmfXSeb9Qvu7Gbi8qq4GSPJd4N6t/F1Jjm8f9F+boI1HA78DfL0t77eA82fQ7ycCBwEXtfmWAj9pr20FPj3N/I9P8jfA7sA+wOUMgTeVJwMHJlndnu8NHAD8H3BhVV3Zyqca92fa3w3Ac9v04xhCjao6O8ldk+wNPAl44diMVfW/bXIXhi8Rr6mqr07TZ0nqzkwD6uaR6a0M38yvr6oVs5j31nHt3ArsXFXfTnIQ8DTgHUnWVdXbxrUR4MtV9aIZ9nd0vtOq6i0TvPaLmuI4UJLdGLa4VlbV/yRZC+w2w2W+rqrOGtfeKoYtKACmGffY+7SVX62j274JjKhWPtEpM7cwBNxTgDkLqH332HfSs/gkaS5t6++gbgCuTHIwQDs287BtaSjJ/sCNVfVx4F3A77WXfgbs1aa/ATw2yQPaPLsneeAE9QB+mWSXNv0VYHWSu7f59klynxl2bSyMrkuyJ7B6knrjl38W8KqxPiR5YJI9xs80xbgncx5wSJt3FXBdVd0ArANeO9LuXdpkAS8BHpTkzdO0PWPXHHUNtaaoNcWaWnPbtKeYSwJYs2bNnLV1R87iOwR4f5JjGHYnfQK4dBva+V3ghCS3Ar8EXtXKPwT8e5Krq+rxSQ4H/iXJ2K7EY4Bvj6/Xnm9KcnFVHdL6ty7JTq391wDfm65TVXV9kpMYdsVdBVw0SdVNwC1JLgVOBf6BYRfoxRn2K17LcNbfTMc9mbXAKUk2MZyxeFgrfzvwvgw/BdgKvJW2i7CqtiZ5IfDFJDdU1T9Ps4xZ8Ye6ksbzh7qaczP9oa4kzbVs6w91JUmaDz3/UHeHSfJZhlPBR71p/IkOkqQdx4ACquo5890HSdLtuYtPktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1KVU1Xz3QR1Ici3wve24iGXAddux/fmyWMcFi3dsi3VcsHDHdp+qutv4QgNKO0SS9VW1cr77MdcW67hg8Y5tsY4LFt/Y3MUnSeqSASVJ6pIBpR3lQ/Pdge1ksY4LFu/YFuu4YJGNzWNQkqQuuQUlSeqSASVJ6pIBpTmV5N5JzklyRZLLkxzRytcm+WGSje3xtPnu62wl2S3JhUkubWN7ayvfJ8mXk3yn/b3LfPd1NqYY14JfZwBJliS5JMmX2vMFvb5GTTC2RbHOxngMSnMqyX7AflV1cZK9gA3As4HnA1uq6l3z2b87IkmAPapqS5JdgP8EjgCeC/y0qo5L8mbgLlX1pvns62xMMa4/YYGvM4AkfwWsBO5UVc9I8k4W8PoaNcHY1rII1tkYt6A0p6rq6qq6uE3/DLgCuOf89mpu1GBLe7pLexTwLOC0Vn4aQyAvGFOMa8FLci/g6cDJI8ULen2NmWRsi4oBpe0myXLg4cAFrei1STYl+chC3a3SdqlsBH4CfLmqLgD2raqrYQho4O7z2MVtMsm4YOGvsxOBvwFuHSlb8OurOZFfHxss/HV2GwNK20WSPYFPA0dW1Q3A+4H7AyuAq4F3z1/vtl1Vba2qFcC9gEcmeeg8d2lOTDKuBb3OkjwD+ElVbZjvvsy1Kca2oNfZeAaU5lw7jvFp4PSq+gxAVf24fQjeCpwEPHI++3hHVdX1wLkMx2l+3I69jR2D+8n89eyOGR3XIlhnjwWemeQq4BPAE5J8nMWxviYc2yJYZ7djQGlOtQPuHwauqKr3jJTvN1LtOcBlO7pvd1SSuyW5c5teCjwJ+C/gC8BhrdphwOfnpYPbaLJxLfR1VlVvqap7VdVy4IXA2VV1KAt8fcHkY1vo62y8nee7A1p0Hgv8ObC5HdMAOBp4UZIVDAffrwJeMR+du4P2A05LsoThy92ZVfWlJOcDZyZ5KfB94OD57OQ2mGxcH1sE62wix7Gw19dU3rmY1pmnmUuSuuQuPklSlwwoSVKXDChJUpcMKElSlwwoSVKXDCipQ0kqybtHnh/VLgQ6F22fmmT1XLQ1zXIOble1P2d7L2uCZR+eZP8dvVzNLQNK6tPNwHOTLJvvjoxqv5WaqZcCr66qx2+v/kzhcMCAWuAMKKlPtwAfAt4w/oXxW0BJtrS/q5J8NcmZSb6d5Lgkh7R7PW1Ocv+RZp6U5Gut3jPa/EuSnJDkonax0VeMtHtOkjOAzRP050Wt/cuSHN/KjgUeB3wgyQkTzPPXI8sZu//U8UlePVJnbZI3TlF/edtCOynDfazWJVna3puVwOkZ7om0tL0X32zzL4pbUfxGqCofPnx09gC2AHdiuBrA3sBRwNr22qnA6tG67e8q4HqGK0PsCvwQeGt77QjgxJH5/4PhC+oBwA+A3YCXA8e0OrsC64H7tnZ/Dtx3gn7uz3A1hrsxXJnmbODZ7bVzgZUTzPNkhvBN68OXgD9kuPL9V0fqfRP47SnqL2cI8hWt/pnAoeOXDewDfItfXZjgzvO9fn3M7OEWlNSpGq4C/1Hg9bOY7aIa7sl1M/DfwLpWvpnhA33MmVV1a1V9B/gu8CCGIHhxu0TVBcBdGQIM4MKqunKC5T0COLeqrq2qW4DTGcJjKk9uj0uAi9uyD6iqS4C7J9k/ycOA/62q709Wv7V1ZVVtbNMbxo1xzA3AL4CTkzwXuHGa/qkTXotP6tuJDB/Kp4yU3ULbPd8uzvtbI6/dPDJ968jzW7n9//fx1zgrhi2U11XVWaMvJFnFsAU1kUzT/8nmeUdVfXCC1z4FrAbuwXCV7knrt/uNjY53K7B0fINVdUuSRwJPZLiw6muBJ2xDv7WDuQUldayqfsqw6+qlI8VXAQe16Wcx3AF3tg5OslM7LnU/hl1gZwGvardLIckDk+wxTTsXAH+UZFk7geJFwFenmecs4CXtnmEkuWeSsZsGfoIhRFYzhNV09SfzM2CvVn9PYO+q+jfgSIZ7JWkBcAtK6t+7Gb71jzkJ+HySC4GvMPnWzVS+xRAk+wKvrKpfJDmZYRfZxW3L7FqmuR16VV2d5C3AOQxbOv9WVVPevqKq1iV5MHD+sBi2AIcy3IDv8iR7AT+sX931drL6W6dYzKkMJ2jcBDyV4f3arfXx1048UZ+8mrkkqUvu4pMkdcmAkiR1yYCSJHXJgJIkdcmAkiR1yYCSJHXJgJIkden/AfCEdh3W3Aw6AAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATPklEQVR4nO3debRdZX3G8e9DoBAGUYwiqDUOuLRajCVO1bZxqNahjsGhUKFa5wGstCqLkuhyVRC1tNY6gAIq1IXz0IG4BMS6EEggJCBVVwWtCgrLsjCCVMKvf+z34uF6x3CT+97r97PWWXef97z73e97dnKes4ezd6oKSZJ6s9N8d0CSpIkYUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVDSApZka5KNSS5L8skku09R95lJ3ryD+7c2yY1J7j5StmVcnae0MWxMsiXJt9r0R2exnMOT7D+Xfdf8M6Ckhe2mqlpRVQ8F/g945WQVq+oLVXXcjuvaba4D3jjZi1V1VhvDCmA9cEh7/uJZLONwwIBaZAwoafH4GvCAJPsk+VySTUm+keRAuG0r45/a9MFtq+vSJOe1sockubBtvWxKckAr/6tW97IkR7ay5UmuSHJSksuTrEuydJJ+fQR4QZJ9ZjOYJIeO9OeDSZa0x6mtL5uTvCHJamAlcHqrO1k/tMAYUNIikGRn4KnAZuCtwCVVdSBwNDDRrrJjgadU1cOAZ7ayVwL/0LZkVgI/SHIQ8BfAo4BHAy9L8vBW/wDgfVX1EOB64HmTdG8LQ0gdMYvxPBh4AfDY1p+twCHACuCeVfXQqvpd4JSq+hS33/K6aabLUd8MKGlhW5pkI8MH9PeBDwOPAz4GUFVnA3dNsve4+b4OnJrkZcCSVnY+cHSSNwH3aR/0jwM+W1U/r6otwGeAP2j1r6yqjW16A7B8in7+I3BYkjvNcFxPBA4CLmrjeyJwP+C7wP2SvDfJnwA3zLA9LUA7z3cHJN0hN7UtjNskyQT1bnfRzap6ZZJHAU8HNiZZUVVnJLmglZ2V5C+Bidoac/PI9FZg0l1rVXV9kjOAV4/08zXAy9rTp1XVj0aHAZxWVW8Z31aShwFPAV4DPB94yRR91ALmFpS0+JzHsDuMJKuA66rqdlsaSe5fVRdU1bEMJzHcO8n9gO9W1T8CXwAObG09O8nuSfYAnsNwrGtbvAd4Be2LcVW9b+zkiHHhBPAVYPXY2X/tuNp9kiwDdqqqTwN/C/xeq/8zYK9t7Jc65RaUtPisBU5Jsgm4EThsgjontJMgwhAGlwJvBg5N8kvgGuBtVfXTJKcCF7b5Tq6qS5Isn22nquq6JJ8F3jCDut9McgywLslOwC8ZtphuamMb+3I9toV1KvCBJDcBj/E41OIQb7chSeqRu/gkSV0yoCRJXTKgJEldMqAkSV3yLD4BsGzZslq+fPl8d0PSb6ANGzZcV1V3G19uQAmA5cuXs379+vnuhqTfQEm+N1G5u/gkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqC0zdauXTvfXZDUmbn8XEhVzVljWrhWrlxZ69evn7bePd51D3788x//Wvm+e+zLNUddsz26JmkBScJscyXJhqpaOb582i2oJMuTXDaLBa1NctSsejdxO0cm2X229ZIcfUeXPYNlPijJxiSXJLl/kj+bw7avSrJsDtrZMhf9GW+icJqqXJK21c7z3YEpHAl8HLhxlvWOBv5uNgtKsqSqts5ilmcDn6+qNUlWAX8GnLEdl7cgrFq1ar67IGkRmekxqCVJTkpyeZJ1SZa2LYf/SLIhydeSPGj8TEnOTfL3Sc5LckWSRyT5TJLvJHl7q7NHkn9NcmmSy5K8IMnrgf2Bc5Kc0+o9Ocn5SS5O8skke46vl+Q4YGnbujm9zXdokgtb2QeTLGnlW5K8LckFwGMmGnSSY5Nc1Pr1oQyexhCKf9n6dhzwB639NyRZkuSENt+mJK9oba1qfTwD2DzRuEcW/bo2zs1j72uSfZJ8rrX5jSQHtvI9k5zS6m5K8rxxY1jW3renTzC+lydZn2T9tddeO8N/CpK0g1TVlA9gOXALsKI9PxM4FPgKcEArexRwdpteCxzVps8Fjm/TRwA/AvYDdgV+ANwVeB5w0sjy9m5/rwKWtellwHnAHu35m4Bjx9drz7eMTD8Y+CKwS3v+z8CL23QBz59m7PuMTH8M+NMJxrgK+NJIvZcDx7TpXYH1wH1bvZ8D922vTTXu17XpVwMnt+n3Amva9BOAjW36eODEkXbuMvY+APsCFwB/PN16Puigg2omWMukD0kaYmXW86yvCT6XZrqL78qq2timNzCE1u8Dn0wyVmfXSeb9Qvu7Gbi8qq4GSPJd4N6t/F1Jjm8f9F+boI1HA78DfL0t77eA82fQ7ycCBwEXtfmWAj9pr20FPj3N/I9P8jfA7sA+wOUMgTeVJwMHJlndnu8NHAD8H3BhVV3Zyqca92fa3w3Ac9v04xhCjao6O8ldk+wNPAl44diMVfW/bXIXhi8Rr6mqr07TZ0nqzkwD6uaR6a0M38yvr6oVs5j31nHt3ArsXFXfTnIQ8DTgHUnWVdXbxrUR4MtV9aIZ9nd0vtOq6i0TvPaLmuI4UJLdGLa4VlbV/yRZC+w2w2W+rqrOGtfeKoYtKACmGffY+7SVX62j274JjKhWPtEpM7cwBNxTgDkLqH332HfSs/gkaS5t6++gbgCuTHIwQDs287BtaSjJ/sCNVfVx4F3A77WXfgbs1aa/ATw2yQPaPLsneeAE9QB+mWSXNv0VYHWSu7f59klynxl2bSyMrkuyJ7B6knrjl38W8KqxPiR5YJI9xs80xbgncx5wSJt3FXBdVd0ArANeO9LuXdpkAS8BHpTkzdO0PWPXHHUNtaaoNcWaWnPbtKeYSwJYs2bNnLV1R87iOwR4f5JjGHYnfQK4dBva+V3ghCS3Ar8EXtXKPwT8e5Krq+rxSQ4H/iXJ2K7EY4Bvj6/Xnm9KcnFVHdL6ty7JTq391wDfm65TVXV9kpMYdsVdBVw0SdVNwC1JLgVOBf6BYRfoxRn2K17LcNbfTMc9mbXAKUk2MZyxeFgrfzvwvgw/BdgKvJW2i7CqtiZ5IfDFJDdU1T9Ps4xZ8Ye6ksbzh7qaczP9oa4kzbVs6w91JUmaDz3/UHeHSfJZhlPBR71p/IkOkqQdx4ACquo5890HSdLtuYtPktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1KVU1Xz3QR1Ici3wve24iGXAddux/fmyWMcFi3dsi3VcsHDHdp+qutv4QgNKO0SS9VW1cr77MdcW67hg8Y5tsY4LFt/Y3MUnSeqSASVJ6pIBpR3lQ/Pdge1ksY4LFu/YFuu4YJGNzWNQkqQuuQUlSeqSASVJ6pIBpTmV5N5JzklyRZLLkxzRytcm+WGSje3xtPnu62wl2S3JhUkubWN7ayvfJ8mXk3yn/b3LfPd1NqYY14JfZwBJliS5JMmX2vMFvb5GTTC2RbHOxngMSnMqyX7AflV1cZK9gA3As4HnA1uq6l3z2b87IkmAPapqS5JdgP8EjgCeC/y0qo5L8mbgLlX1pvns62xMMa4/YYGvM4AkfwWsBO5UVc9I8k4W8PoaNcHY1rII1tkYt6A0p6rq6qq6uE3/DLgCuOf89mpu1GBLe7pLexTwLOC0Vn4aQyAvGFOMa8FLci/g6cDJI8ULen2NmWRsi4oBpe0myXLg4cAFrei1STYl+chC3a3SdqlsBH4CfLmqLgD2raqrYQho4O7z2MVtMsm4YOGvsxOBvwFuHSlb8OurOZFfHxss/HV2GwNK20WSPYFPA0dW1Q3A+4H7AyuAq4F3z1/vtl1Vba2qFcC9gEcmeeg8d2lOTDKuBb3OkjwD+ElVbZjvvsy1Kca2oNfZeAaU5lw7jvFp4PSq+gxAVf24fQjeCpwEPHI++3hHVdX1wLkMx2l+3I69jR2D+8n89eyOGR3XIlhnjwWemeQq4BPAE5J8nMWxviYc2yJYZ7djQGlOtQPuHwauqKr3jJTvN1LtOcBlO7pvd1SSuyW5c5teCjwJ+C/gC8BhrdphwOfnpYPbaLJxLfR1VlVvqap7VdVy4IXA2VV1KAt8fcHkY1vo62y8nee7A1p0Hgv8ObC5HdMAOBp4UZIVDAffrwJeMR+du4P2A05LsoThy92ZVfWlJOcDZyZ5KfB94OD57OQ2mGxcH1sE62wix7Gw19dU3rmY1pmnmUuSuuQuPklSlwwoSVKXDChJUpcMKElSlwwoSVKXDCipQ0kqybtHnh/VLgQ6F22fmmT1XLQ1zXIOble1P2d7L2uCZR+eZP8dvVzNLQNK6tPNwHOTLJvvjoxqv5WaqZcCr66qx2+v/kzhcMCAWuAMKKlPtwAfAt4w/oXxW0BJtrS/q5J8NcmZSb6d5Lgkh7R7PW1Ocv+RZp6U5Gut3jPa/EuSnJDkonax0VeMtHtOkjOAzRP050Wt/cuSHN/KjgUeB3wgyQkTzPPXI8sZu//U8UlePVJnbZI3TlF/edtCOynDfazWJVna3puVwOkZ7om0tL0X32zzL4pbUfxGqCofPnx09gC2AHdiuBrA3sBRwNr22qnA6tG67e8q4HqGK0PsCvwQeGt77QjgxJH5/4PhC+oBwA+A3YCXA8e0OrsC64H7tnZ/Dtx3gn7uz3A1hrsxXJnmbODZ7bVzgZUTzPNkhvBN68OXgD9kuPL9V0fqfRP47SnqL2cI8hWt/pnAoeOXDewDfItfXZjgzvO9fn3M7OEWlNSpGq4C/1Hg9bOY7aIa7sl1M/DfwLpWvpnhA33MmVV1a1V9B/gu8CCGIHhxu0TVBcBdGQIM4MKqunKC5T0COLeqrq2qW4DTGcJjKk9uj0uAi9uyD6iqS4C7J9k/ycOA/62q709Wv7V1ZVVtbNMbxo1xzA3AL4CTkzwXuHGa/qkTXotP6tuJDB/Kp4yU3ULbPd8uzvtbI6/dPDJ968jzW7n9//fx1zgrhi2U11XVWaMvJFnFsAU1kUzT/8nmeUdVfXCC1z4FrAbuwXCV7knrt/uNjY53K7B0fINVdUuSRwJPZLiw6muBJ2xDv7WDuQUldayqfsqw6+qlI8VXAQe16Wcx3AF3tg5OslM7LnU/hl1gZwGvardLIckDk+wxTTsXAH+UZFk7geJFwFenmecs4CXtnmEkuWeSsZsGfoIhRFYzhNV09SfzM2CvVn9PYO+q+jfgSIZ7JWkBcAtK6t+7Gb71jzkJ+HySC4GvMPnWzVS+xRAk+wKvrKpfJDmZYRfZxW3L7FqmuR16VV2d5C3AOQxbOv9WVVPevqKq1iV5MHD+sBi2AIcy3IDv8iR7AT+sX931drL6W6dYzKkMJ2jcBDyV4f3arfXx1048UZ+8mrkkqUvu4pMkdcmAkiR1yYCSJHXJgJIkdcmAkiR1yYCSJHXJgJIkden/AfCEdh3W3Aw6AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -416,7 +446,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU40lEQVR4nO3df7xldV3v8ddb4PLTi8IoRhJjCg8vyDQ5449CcyREs6uhYoJwA62USr3oxRDlyuj1FgSllj9KvEH+pChRSAIUUYhAmIFhAC008AcPRZkMaUD5MXzuH+t7YHs4P2fOzPkeej0fj/M4a3/Xd33XZ+199n6ftfbea6WqkCSpN4+Y7wIkSZqIASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElLUBJNiRZk+T6JGcl2WGKvi9O8pYtXN/KJJXkSSNtb2xty0fa9mvbsSbJD5Lc3KY/P4t1HZxkn7neBs0/A0pamH5UVUur6inAPcDRk3WsqnOq6qQtV9oDrgMOHbl9CPCV0Q5VdV3bjqXAOcCb2+0DZ7GegwED6mHIgJIWvkuBJyXZJcmnk6xNckWSJQBJjkryvjb98rbXdW2SS1rbvkmubHsua5Ps1drf1Ppen+SY1rY4yVeTnJbkhiQXJtl+kro+DfxaW+5ngR8Ct81kg5IclOTyJFe3PcSdWvtJSb7S6jw1yS8CLwZOafU/caPuQXXJgJIWsCRbA7/CsLfyDuCaqloCvBX4yASLvB14flX9HMMLOwx7X+9tezHLgVuSLANeBTwDeCbw20l+vvXfC3h/Ve0L3A68bJLy7gC+neQpwGHAX89wmxYBJwAHVtVTgVXAm5LsArwE2Ldt47uq6p/4yT2vf53JOrQwGFDSwrR9kjUML97fAv4f8CzgowBV9QVg1yQ7j1vuMuCMJL8NbNXaLgfemuQ4YM+q+lEb6+yqurOq1gOfAp7d+t9cVWva9Gpg8RR1nslwmO9g4OwZbtszGQ7ZXda28UhgT4bA+zHw4SQvBe6a4XhaoAwoaWEaew9qaVW9vqruATJBv5842WZVHc2wd7IHsCbJrlX1CYa9qR8BFyQ5YJKxxtw9Mr0B2HqKvucC/wP4VlXdMdaY5CUjH45YPm6ZAJ8b2b59quo3q+o+4OnA3zEE3vlTrFcPAwaU9PBxCXA4QJIVwLrRUGjtT6yqL1fV24F1wB7t/aGbqupPGQ6XLWljHZxkhyQ7Mhxau3S2BbW9seOA/zuu/eyRAFo1brErgP3HPgHYati7vQ+1c1WdBxwDLG39/wN45GxrU/+m+s9H0sKyEjg9yVqGw19HTtDnlPYhiAAXAdcCbwGOSHIvcCvwzqr6QZIzgCvbch+uqmuSLJ5tUVV15iz735bkKOCTSbZtzScwBNFnkmzX6n9jm3cmcFqSNwCH+D7Uw0e83IYkqUce4pMkdcmAkiR1yYCSJHXJgJIkdclP8QmARYsW1eLFi+e7DEn/iaxevXpdVT1msvkGlABYvHgxq1aN/zqKJG0+Sb451XwP8UmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQGmjrVy5cr5LkNSJzfF6kKqa80G18CxfvrxWrVo1q2WSMN3fz+NOfRzfu/N7D2nfbcfduPXYW2e1Pkn9msnrwQTLrK6q5ZPNn3YPKsniJNfPYoUrkxw70/5TjHNMkh1m2y/JWzd13TNY55OTrElyTZInJnnlHI79jSSL5mCc9XNRz6aaKJymapekMVvPdwFTOAb4GHDXLPu9FfiD2awoyVZVtWEWixwMfKaqTkyyAngl8InNuL5urVixYuoOz92EZSX9pzbT96C2SnJakhuSXJhk+7bncH6S1UkuTfLk8Qsl+WKSdye5JMlXkzwtyaeSfC3Ju1qfHZN8Nsm1Sa5P8ookbwB2By5OcnHrd1CSy5NcneSsJDuN75fkJGD7tnfz8bbcEUmubG1/kWSr1r4+yTuTfBn4hYk2Osnbk1zV6vpQBi9kCMXfarWdBDy7jf/GJFslOaUttzbJa9tYK1qNnwCum2i7R1b9+rad143dr0l2SfLpNuYVSZa09p2SnN76rk3ysnHbsKjdb786wfa9JsmqJKtuu+22Gf4pSNIWUlVT/gCLgfuApe323wBHABcBe7W2ZwBfaNMrgWPb9BeBk9v0/wS+A/wUsC1wC7Ar8DLgtJH17dx+fwNY1KYXAZcAO7bbxwFvH9+v3V4/Mv3fgHOBbdrtDwC/0aYL+PVptn2XkemPAi+aYBtXAH8/0u81wAlteltgFfCE1u9O4Alt3lTb/fo2/bvAh9v0nwEntukDgDVt+mTgPSPjPHrsfgB2A74MPG+6x3nZsmU1W8OfzzR9VjLpj6SHj5m8HkywzKqa4nVppof4bq6qNW16NUNo/SJwVpKxPttOsuw57fd1wA1V9V2AJDcBe7T2U5Oc3F7oL51gjGcC+wCXtfX9F+DyGdT9y8Ay4Kq23PbA99u8DcDfTbP8c5P8PrADsAtwA0PgTeUgYEmSQ9rtnYG9gHuAK6vq5tY+1XZ/qv1eDby0TT+LIdSoqi8k2TXJzsCBwKFjC1bVv7fJbRj+ifi9qvrSNDVLUndmGlB3j0xvYPjP/PaqWjqLZe8fN879wNZVdWOSZcALgT9McmFVvXPcGAE+V1WHzbDe0eX+qqqOn2Dej2uK94GSbMewx7W8qr6dZCWw3QzX+fqqumDceCsY9qAAmGa7x+6nDTz4GD3wn8CIau0TfXTmPoaAez4wbwG12467TfopPkmaysZ+D+oO4OYkLwdo78383MYMlGR34K6q+hhwKvDUNus/gEe26SuA/ZM8qS2zQ5K9J+gHcG+Sbdr0RcAhSR7bltslyZ4zLG0sjNYl2Qk4ZJJ+49d/AfA7YzUk2TvJjuMXmmK7J3MJcHhbdgWwrqruAC4EXjcy7qPbZAGvBp6c5C3TjL1RTjzxxGn73HrsrdSJ9ZAfP2IuPbzM5PVgtjblU3yHAx9McgLD4aQzgWs3Ypz9gFOS3A/cC/xOa/8Q8A9JvltVz01yFPDJJGOHEk8Abhzfr91em+Tqqjq81Xdhkke08X8P+OZ0RVXV7UlOYzgU9w3gqkm6rgXuS3ItcAbwXoZDoFdnOK54G8On/ma63ZNZCZyeZC3DJxaPbO3vAt6f4asAG4B30A4RVtWGJIcC5ya5o6o+MM06ZsUv6koa4xd1tdlszBd1JWlTZFO/qCtJ0nzo+Yu6W0ySsxk+Cj7quPEfdJAkbTkGFFBVL5nvGiRJP8lDfJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuparmuwZ1IMltwDfneNhFwLo5HnNzsda5t1DqhIVT60KpE2ZW655V9ZjJZhpQ2mySrKqq5fNdx0xY69xbKHXCwql1odQJc1Orh/gkSV0yoCRJXTKgtDl9aL4LmAVrnXsLpU5YOLUulDphDmr1PShJUpfcg5IkdcmAkiR1yYDSnEtySpJ/TrI2ydlJHtXat0nyV0muS/LVJMfPc6mT1trmLUlyeZIbWs3b9Vhnm/8zSdYnOXaeShytZbLH/3lJVrf7cnWSA3qss807PsnXk/xLkufPY5lj9by8/R3en2T5SHuPz6kJa23zZvWcMqC0OXwOeEpVLQFuBMaeNC8Htq2q/YBlwGuTLJ6fEh8wYa1JtgY+BhxdVfsCK4B756tIJr9Px7wb+IctXtXEJqt1HfCi9vgfCXx0nuobM9ljvw9wKLAv8ALgA0m2mrcqB9cDLwUuGdfe43Nqwlo35jllQGnOVdWFVXVfu3kF8PixWcCO7Q91e+Ae4I55KPEBU9R6ELC2qq5t/f6tqjbMR41t/ZPVSZKDgZuAG+ahtIeYrNaquqaqvtPabwC2S7LtfNTY6pnsPv014Myquruqbga+Djx9PmocU1Vfrap/mWgW/T2nJqt11s8pA0qb26t58D/7vwXuBL4LfAs4tap+MF+FTWC01r2BSnJBkquT/P481jXeA3Um2RE4DnjHvFY0udH7dNTLgGuq6u4tXM9kRuv8aeDbI/NuaW096v05NWrWz6mtt0BRehhK8nngcRPMeltVfab1eRtwH/DxNu/pwAZgd+DRwKVJPl9VN3VY69bAs4CnAXcBFyVZXVUXdVbnO4B3V9X6JJurtIfYyFrHlt0XOJnhP+oe65zojtzs38eZSa0T6PY5NYFZP6cMKG2UqjpwqvlJjgT+O/DL9eCX7V4JnF9V9wLfT3IZsJzh8FRvtd4CfKmq1rU+5wFPBTZbQG1knc8ADknyR8CjgPuT/Liq3re56tyEWknyeOBs4Deq6l83Z42wSY/9HiPdHg98Z/yyc226WifR5XNqErN+TnmIT3MuyQsYDju9uKruGpn1LeCADHYEngn883zUOGaKWi8AliTZoR3ffw7wlfmoESavs6qeXVWLq2ox8B7gDzZ3OE1nslrbp+Q+CxxfVZfNU3kPmOKxPwc4NMm2SZ4A7AVcOR81zkB3z6kpzPo55ZkkNOeSfB3YFvi31nRFVR2dZCfgdGAfhsMop1fVKfNUJjB5rW3eEQyf7CrgvKqat/ehpqpzpM9KYH1VnbqFy/sJUzz+JzDcn18b6X5QVX1/S9cI0z72b2N4X+o+4JiqmtdPSCZ5CfBnwGOA24E1VfX8Tp9TE9ba5s3qOWVASZK65CE+SVKXDChJUpcMKElSlwwoSVKXDChJUpcMKGkLSlJJ/njk9rHt4+FbsoYvjp1lOsl5GXdm9I0Yb0WSv5+qPcmLk7ylTZ+R5JBZjL++/d49yd+26aOSzPn3vTbXuNo4BpS0Zd0NvDTJoo1ZuH3Bcc5U1Qur6va5HHOS9ZxTVSdt4hjfqaoZB5sWPgNK2rLuAz4EvHH8jCR7Jrkow/WJLkryM639jCR/kuRi4OR2+4NJLk5yU5LnJPnLDNcDOmNkvA8mWZXh2jsTnkw2yTeSLEpydJI17efmti6SHJTh+j1XJzmrfTGUJC/IcC2lf2S4tMKUJtszSfJ/2vY8Ismbk1zVtv8h9SZZnOT6kabdk5yf5GvtVE9j/Q7LcK2h65OcPIP2VyW5McmXgP2n2xZtOQaUtOW9Hzg8yc7j2t8HfKRdn+jjwJ+OzNsbOLCq/le7/WjgAIagO5fhelD7AvslWdr6vK2qlgNLgOckWTJZQVX151W1lOFEnrcAf9L28k5o630qsAp4U4aLzJ0GvAh4NhOfNHRaLVQeC7wKOJDhlEJPB5YCy5L80jRDLAVeAewHvCLJHkl2ZzgR7QFt/tOSHDxF+08xnHB3f+B5DGdkUCcMKGkLq6o7gI8Abxg36xeAT7TpjzKc+XnMWeOunXNuO7npdcD3quq6qrqf4TpLi1ufX09yNXANQ3jN5MX3vcAXqupchvO67QNclmQNw0UG9wSeDNxcVV9rNXxsBuOO97+BR1XVa9sYB7Wfa4Cr2zr2mmaMi6rqh1X1Y4Zzuu3JELBfrKrb2rWePg780hTtzxhpvwf4643YFm0mns1cmh/vYXghPn2KPqPnIbtz3Lyx6yjdPzI9dnvrdpLTY4GnVdW/t0N/U15eO8lRDC/yrxtrAj5XVYeN67eUTb/8xFUMe0m7tOsXBfjDqvqLWYwxut0bGF7PJrvmyFTXIvF8b51yD0qaB+1F+W+A3xxp/ieGS40DHA784yas4r8yhNoPk+wG/MpUnZMsYwi0I9qeGAxXmd0/yZNanx2S7M1wtuwnJHli63fYQwac3vnAScBnkzyS4UzXrx55j+unkzx2I8b9MsPhzEUZLtN+GPCladpXJNk1yTYMl1BXJ9yDkubPH/Pg3goMh/z+MsmbgdsY3pvZKFV1bZJrGA753QRMd3mL1wG7ABdnuPDhqqr6rbZX9ck8eGn2E6rqxiSvYQiXdQxB+pSNqPGsFk7nAC9kOLx5eVv/euAIYFZnOq+q7yY5HriYYa/pvJELE07WvhK4nOGqtFcDW812W7R5eDZzSVKXPMQnSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSerS/weph5TRs1E/+QAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU40lEQVR4nO3df7xldV3v8ddb4PLTi8IoRhJjCg8vyDQ5449CcyREs6uhYoJwA62USr3oxRDlyuj1FgSllj9KvEH+pChRSAIUUYhAmIFhAC008AcPRZkMaUD5MXzuH+t7YHs4P2fOzPkeej0fj/M4a3/Xd33XZ+199n6ftfbea6WqkCSpN4+Y7wIkSZqIASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElLUBJNiRZk+T6JGcl2WGKvi9O8pYtXN/KJJXkSSNtb2xty0fa9mvbsSbJD5Lc3KY/P4t1HZxkn7neBs0/A0pamH5UVUur6inAPcDRk3WsqnOq6qQtV9oDrgMOHbl9CPCV0Q5VdV3bjqXAOcCb2+0DZ7GegwED6mHIgJIWvkuBJyXZJcmnk6xNckWSJQBJjkryvjb98rbXdW2SS1rbvkmubHsua5Ps1drf1Ppen+SY1rY4yVeTnJbkhiQXJtl+kro+DfxaW+5ngR8Ct81kg5IclOTyJFe3PcSdWvtJSb7S6jw1yS8CLwZOafU/caPuQXXJgJIWsCRbA7/CsLfyDuCaqloCvBX4yASLvB14flX9HMMLOwx7X+9tezHLgVuSLANeBTwDeCbw20l+vvXfC3h/Ve0L3A68bJLy7gC+neQpwGHAX89wmxYBJwAHVtVTgVXAm5LsArwE2Ldt47uq6p/4yT2vf53JOrQwGFDSwrR9kjUML97fAv4f8CzgowBV9QVg1yQ7j1vuMuCMJL8NbNXaLgfemuQ4YM+q+lEb6+yqurOq1gOfAp7d+t9cVWva9Gpg8RR1nslwmO9g4OwZbtszGQ7ZXda28UhgT4bA+zHw4SQvBe6a4XhaoAwoaWEaew9qaVW9vqruATJBv5842WZVHc2wd7IHsCbJrlX1CYa9qR8BFyQ5YJKxxtw9Mr0B2HqKvucC/wP4VlXdMdaY5CUjH45YPm6ZAJ8b2b59quo3q+o+4OnA3zEE3vlTrFcPAwaU9PBxCXA4QJIVwLrRUGjtT6yqL1fV24F1wB7t/aGbqupPGQ6XLWljHZxkhyQ7Mhxau3S2BbW9seOA/zuu/eyRAFo1brErgP3HPgHYati7vQ+1c1WdBxwDLG39/wN45GxrU/+m+s9H0sKyEjg9yVqGw19HTtDnlPYhiAAXAdcCbwGOSHIvcCvwzqr6QZIzgCvbch+uqmuSLJ5tUVV15iz735bkKOCTSbZtzScwBNFnkmzX6n9jm3cmcFqSNwCH+D7Uw0e83IYkqUce4pMkdcmAkiR1yYCSJHXJgJIkdclP8QmARYsW1eLFi+e7DEn/iaxevXpdVT1msvkGlABYvHgxq1aN/zqKJG0+Sb451XwP8UmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQEmSumRASZK6ZEBJkrpkQGmjrVy5cr5LkNSJzfF6kKqa80G18CxfvrxWrVo1q2WSMN3fz+NOfRzfu/N7D2nfbcfduPXYW2e1Pkn9msnrwQTLrK6q5ZPNn3YPKsniJNfPYoUrkxw70/5TjHNMkh1m2y/JWzd13TNY55OTrElyTZInJnnlHI79jSSL5mCc9XNRz6aaKJymapekMVvPdwFTOAb4GHDXLPu9FfiD2awoyVZVtWEWixwMfKaqTkyyAngl8InNuL5urVixYuoOz92EZSX9pzbT96C2SnJakhuSXJhk+7bncH6S1UkuTfLk8Qsl+WKSdye5JMlXkzwtyaeSfC3Ju1qfHZN8Nsm1Sa5P8ookbwB2By5OcnHrd1CSy5NcneSsJDuN75fkJGD7tnfz8bbcEUmubG1/kWSr1r4+yTuTfBn4hYk2Osnbk1zV6vpQBi9kCMXfarWdBDy7jf/GJFslOaUttzbJa9tYK1qNnwCum2i7R1b9+rad143dr0l2SfLpNuYVSZa09p2SnN76rk3ysnHbsKjdb786wfa9JsmqJKtuu+22Gf4pSNIWUlVT/gCLgfuApe323wBHABcBe7W2ZwBfaNMrgWPb9BeBk9v0/wS+A/wUsC1wC7Ar8DLgtJH17dx+fwNY1KYXAZcAO7bbxwFvH9+v3V4/Mv3fgHOBbdrtDwC/0aYL+PVptn2XkemPAi+aYBtXAH8/0u81wAlteltgFfCE1u9O4Alt3lTb/fo2/bvAh9v0nwEntukDgDVt+mTgPSPjPHrsfgB2A74MPG+6x3nZsmU1W8OfzzR9VjLpj6SHj5m8HkywzKqa4nVppof4bq6qNW16NUNo/SJwVpKxPttOsuw57fd1wA1V9V2AJDcBe7T2U5Oc3F7oL51gjGcC+wCXtfX9F+DyGdT9y8Ay4Kq23PbA99u8DcDfTbP8c5P8PrADsAtwA0PgTeUgYEmSQ9rtnYG9gHuAK6vq5tY+1XZ/qv1eDby0TT+LIdSoqi8k2TXJzsCBwKFjC1bVv7fJbRj+ifi9qvrSNDVLUndmGlB3j0xvYPjP/PaqWjqLZe8fN879wNZVdWOSZcALgT9McmFVvXPcGAE+V1WHzbDe0eX+qqqOn2Dej2uK94GSbMewx7W8qr6dZCWw3QzX+fqqumDceCsY9qAAmGa7x+6nDTz4GD3wn8CIau0TfXTmPoaAez4wbwG12467TfopPkmaysZ+D+oO4OYkLwdo78383MYMlGR34K6q+hhwKvDUNus/gEe26SuA/ZM8qS2zQ5K9J+gHcG+Sbdr0RcAhSR7bltslyZ4zLG0sjNYl2Qk4ZJJ+49d/AfA7YzUk2TvJjuMXmmK7J3MJcHhbdgWwrqruAC4EXjcy7qPbZAGvBp6c5C3TjL1RTjzxxGn73HrsrdSJ9ZAfP2IuPbzM5PVgtjblU3yHAx9McgLD4aQzgWs3Ypz9gFOS3A/cC/xOa/8Q8A9JvltVz01yFPDJJGOHEk8Abhzfr91em+Tqqjq81Xdhkke08X8P+OZ0RVXV7UlOYzgU9w3gqkm6rgXuS3ItcAbwXoZDoFdnOK54G8On/ma63ZNZCZyeZC3DJxaPbO3vAt6f4asAG4B30A4RVtWGJIcC5ya5o6o+MM06ZsUv6koa4xd1tdlszBd1JWlTZFO/qCtJ0nzo+Yu6W0ySsxk+Cj7quPEfdJAkbTkGFFBVL5nvGiRJP8lDfJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuparmuwZ1IMltwDfneNhFwLo5HnNzsda5t1DqhIVT60KpE2ZW655V9ZjJZhpQ2mySrKqq5fNdx0xY69xbKHXCwql1odQJc1Orh/gkSV0yoCRJXTKgtDl9aL4LmAVrnXsLpU5YOLUulDphDmr1PShJUpfcg5IkdcmAkiR1yYDSnEtySpJ/TrI2ydlJHtXat0nyV0muS/LVJMfPc6mT1trmLUlyeZIbWs3b9Vhnm/8zSdYnOXaeShytZbLH/3lJVrf7cnWSA3qss807PsnXk/xLkufPY5lj9by8/R3en2T5SHuPz6kJa23zZvWcMqC0OXwOeEpVLQFuBMaeNC8Htq2q/YBlwGuTLJ6fEh8wYa1JtgY+BhxdVfsCK4B756tIJr9Px7wb+IctXtXEJqt1HfCi9vgfCXx0nuobM9ljvw9wKLAv8ALgA0m2mrcqB9cDLwUuGdfe43Nqwlo35jllQGnOVdWFVXVfu3kF8PixWcCO7Q91e+Ae4I55KPEBU9R6ELC2qq5t/f6tqjbMR41t/ZPVSZKDgZuAG+ahtIeYrNaquqaqvtPabwC2S7LtfNTY6pnsPv014Myquruqbga+Djx9PmocU1Vfrap/mWgW/T2nJqt11s8pA0qb26t58D/7vwXuBL4LfAs4tap+MF+FTWC01r2BSnJBkquT/P481jXeA3Um2RE4DnjHvFY0udH7dNTLgGuq6u4tXM9kRuv8aeDbI/NuaW096v05NWrWz6mtt0BRehhK8nngcRPMeltVfab1eRtwH/DxNu/pwAZgd+DRwKVJPl9VN3VY69bAs4CnAXcBFyVZXVUXdVbnO4B3V9X6JJurtIfYyFrHlt0XOJnhP+oe65zojtzs38eZSa0T6PY5NYFZP6cMKG2UqjpwqvlJjgT+O/DL9eCX7V4JnF9V9wLfT3IZsJzh8FRvtd4CfKmq1rU+5wFPBTZbQG1knc8ADknyR8CjgPuT/Liq3re56tyEWknyeOBs4Deq6l83Z42wSY/9HiPdHg98Z/yyc226WifR5XNqErN+TnmIT3MuyQsYDju9uKruGpn1LeCADHYEngn883zUOGaKWi8AliTZoR3ffw7wlfmoESavs6qeXVWLq2ox8B7gDzZ3OE1nslrbp+Q+CxxfVZfNU3kPmOKxPwc4NMm2SZ4A7AVcOR81zkB3z6kpzPo55ZkkNOeSfB3YFvi31nRFVR2dZCfgdGAfhsMop1fVKfNUJjB5rW3eEQyf7CrgvKqat/ehpqpzpM9KYH1VnbqFy/sJUzz+JzDcn18b6X5QVX1/S9cI0z72b2N4X+o+4JiqmtdPSCZ5CfBnwGOA24E1VfX8Tp9TE9ba5s3qOWVASZK65CE+SVKXDChJUpcMKElSlwwoSVKXDChJUpcMKGkLSlJJ/njk9rHt4+FbsoYvjp1lOsl5GXdm9I0Yb0WSv5+qPcmLk7ylTZ+R5JBZjL++/d49yd+26aOSzPn3vTbXuNo4BpS0Zd0NvDTJoo1ZuH3Bcc5U1Qur6va5HHOS9ZxTVSdt4hjfqaoZB5sWPgNK2rLuAz4EvHH8jCR7Jrkow/WJLkryM639jCR/kuRi4OR2+4NJLk5yU5LnJPnLDNcDOmNkvA8mWZXh2jsTnkw2yTeSLEpydJI17efmti6SHJTh+j1XJzmrfTGUJC/IcC2lf2S4tMKUJtszSfJ/2vY8Ismbk1zVtv8h9SZZnOT6kabdk5yf5GvtVE9j/Q7LcK2h65OcPIP2VyW5McmXgP2n2xZtOQaUtOW9Hzg8yc7j2t8HfKRdn+jjwJ+OzNsbOLCq/le7/WjgAIagO5fhelD7AvslWdr6vK2qlgNLgOckWTJZQVX151W1lOFEnrcAf9L28k5o630qsAp4U4aLzJ0GvAh4NhOfNHRaLVQeC7wKOJDhlEJPB5YCy5L80jRDLAVeAewHvCLJHkl2ZzgR7QFt/tOSHDxF+08xnHB3f+B5DGdkUCcMKGkLq6o7gI8Abxg36xeAT7TpjzKc+XnMWeOunXNuO7npdcD3quq6qrqf4TpLi1ufX09yNXANQ3jN5MX3vcAXqupchvO67QNclmQNw0UG9wSeDNxcVV9rNXxsBuOO97+BR1XVa9sYB7Wfa4Cr2zr2mmaMi6rqh1X1Y4Zzuu3JELBfrKrb2rWePg780hTtzxhpvwf4643YFm0mns1cmh/vYXghPn2KPqPnIbtz3Lyx6yjdPzI9dnvrdpLTY4GnVdW/t0N/U15eO8lRDC/yrxtrAj5XVYeN67eUTb/8xFUMe0m7tOsXBfjDqvqLWYwxut0bGF7PJrvmyFTXIvF8b51yD0qaB+1F+W+A3xxp/ieGS40DHA784yas4r8yhNoPk+wG/MpUnZMsYwi0I9qeGAxXmd0/yZNanx2S7M1wtuwnJHli63fYQwac3vnAScBnkzyS4UzXrx55j+unkzx2I8b9MsPhzEUZLtN+GPCladpXJNk1yTYMl1BXJ9yDkubPH/Pg3goMh/z+MsmbgdsY3pvZKFV1bZJrGA753QRMd3mL1wG7ABdnuPDhqqr6rbZX9ck8eGn2E6rqxiSvYQiXdQxB+pSNqPGsFk7nAC9kOLx5eVv/euAIYFZnOq+q7yY5HriYYa/pvJELE07WvhK4nOGqtFcDW812W7R5eDZzSVKXPMQnSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSerS/weph5TRs1E/+QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -499,7 +529,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVQUlEQVR4nO3de7QlZX3m8e/DZZCLotCYYNRuMCgXRQwgopJ0ouKo0cEFJo6NIMYQHSOQGTMqsLDBy8AwK0ZHSWwcwRg0E9eIgg4BgrQkyK0bmm4ISxPlMghJYEXUBkUuv/mj3mNvjufa3fR5+/D9rLXXfvdbVW+9v304/VC196lKVSFJUm+2mOsJSJI0EQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDStrMJHkkyaokNyX5UpLtplj3DUnev4nn97wky9scb0mybNzyk9qyVSO1rEpy3Cz2ceLGn7l6E/8OStq8JFlbVTu09nnAyqr6kzme1s8luRg4q6q+2l6/oKrWTLLuz2uZ5T7WazttXjyCkjZvfwf8apKdknwlyeokVyfZFyDJ25J8srXf1I66bkxyRevbJ8m17QhmdZI9Wv9/buvelOSE1reoHRGdneTmJJck2XaCOe0K3Dn2YrJwGpVkyyRnJrmuzeMPWv+uSa4YOWI8JMnpwLat77wNevfUNQNK2kwl2Qp4DbAGOBW4oar2BU4E/mKCTU4BXl1VLwTe0PreCXy8qvYDDgDuTLI/cAxwEPAS4PeTvKitvwfwqaraB7gPOHyC/XwM+EaSi5L8UZKnzqCc3wN+WFUHAge2fe4GvAW4uM3vhcCqqno/8JOq2q+qlsxgbG2mDChp87NtklXACuAO4H8BLwc+D1BV3wB2TrLjuO2uBM5N8vvAlq3vKuDEJO8DFlbVT9pY51fV/VW1FvgycEhb/9aqWtXaK4FF4ydXVecAewFfAhYDVyfZZpqaDgWOanVdA+zMEIbXAcckWQq8oKp+PM04mke2musJSJq1n7Qjip9LkgnWe8wHzFX1ziQHAa8DViXZr6q+kOSa1ndxkncAE4015sGR9iPARKf4qKq7gM8Cn01yE/D8JH8IvAi4q6peO26TAO+pqovHj5Xk19v8Pp/kzKqa6OhQ85BHUNL8cAWwBCDJYuDeqvrR6ApJnlNV11TVKcC9wLOS7A58r6o+AVwA7NvGOizJdkm2B97I8FnXjCT590m2bu1fZjga+n5VHdNOy40PJ4CLgXeNbPfcJNsnWQj8a1WdzXCk+Gtt/YfG1tX85RGUND8sBc5Jshp4ADh6gnXObF+CCHAZcCPwfuDIJA8B/wycVlX/luRc4Nq23Weq6oYki2Y4l0OBjyf5aXv9x1X1z9Ns8xmG04XXt6PBe4DDGE4R/nGb31rgqLb+MmB1kuv9HGr+8mvmkqQueYpPktQlA0qS1CUDSpLUJQNKktQlv8UnABYsWFCLFi2a62lImodWrlx5b1XtMtvtDCgBsGjRIlasWDHX05A0DyW5fX228xSfJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUNoqlS5fO9RQkbSKb6vfdgNLsnXceLFoEW2wxPJ93Hqeeeupcz0rSJrKpft+nDagki5LcNNMBkyxN8t4NmxYkOSHJdrNdL8mJG7rvGexzzySrktyQ5DlJ3rIRx74tyYKNMM7ajTGfX3DeeXDssXD77VA1PB97LP/xcdmZpCeyreZ6AlM4AfhL4IFZrnci8NHZ7CjJllX1yCw2OQz4alV9MMli4C3AFx7H/fXjpJPggXE/kgce4KPA4sWL52JGkuapmZ7i2zLJ2UluTnJJkm3bkcPfJFmZ5O+S7Dl+oyTLk3wsyRVJbklyYJIvJ/nHJB9u62yf5OtJbkxyU5LfTXIc8Azg8iSXt/UOTXJVkuuTfCnJDuPXS3I6sG07ujmvbXdkkmtb36eTbNn61yY5Lck1wMETFZ3klCTXtXkty+C1DKH4jja304FD2vh/lGTLJGe27VYn+YM21uI2xy8Aayaqe2TX72l1rhl7X5PslOQrbcyrk+zb+ndIck5bd3WSw8fVsKC9b6+boL5jk6xIsuKee+6Z2X8Jd9wxYfezZ7a1JM1cVU35ABYBDwP7tdd/DRwJXAbs0foOAr7R2kuB97b2cuCM1j4euAvYFdgGuBPYGTgcOHtkfzu259uABa29ALgC2L69fh9wyvj12uu1I+29gAuBrdvrs4CjWruA35mm9p1G2p8HXj9BjYuBr42sdyxwcmtvA6wAdmvr3Q/s1pZNVfd7Wvs/AZ9p7f8JfLC1fwtY1dpnAH86Ms7Txt4H4JeAa4BXTfdz3n///WtGFi6sGk7uPeZxK8xse0mbPWb5+w6sqGn+DZroMdMjqFuralVrr2QIrZcCX0qyCvg0Q/BM5IL2vAa4uarurqoHge8Bz2r9r0xyRpJDquqHE4zxEmBv4Mq2v6OBhTOY9yuA/YHr2navAHZvyx4B/s802/9mkmuSrGEIhX1msM9DgaPa/q5hCOE92rJrq+rW1p6q7i+357H3GuDlDCFJVX0D2DnJjsArgU+NbVhVP2jNrRn+J+K/VtWlM5j3zHzkI7DduI8Gt9uOx/2DP0lPODMNqAdH2o8AOwH3VdV+I4+9ptn20XHjPApsVVXfYQiRNcB/S3LKBGMEuHRkX3tX1e/NYN4BPjey3fOqamlb9tOa4nOgJE9iOOI6oqpeAJwNPGmG+3zPyD53q6pL2rL7x1aapu6x9+kR1n1OmAn2Va2/Jlj2MEPAvXoGc565JUtg2TJYuBCS4XnZMr64UXciSev/NfMfAbcmeRNA+2zmheszUJJnAA9U1V8C/wP4tbbox8CTW/tq4GVJfrVts12S506wHsBDSbZu7cuAI5I8vW23U5KZHHnBujC6N8kOwBGTrDd+/xcD7xqbQ5LnJtl+/EZT1D2ZK4AlbdvFwL1V9SPgEuAPR8Z9WmsW8HZgzyTvn2bs2VmyBG67DR59dHhesoQPfvCDG3UXkvq1qX7fN+RbfEuAP0tyMsPppL8CblyPcV4AnJnkUeAh4F2tfxlwUZK7q+o3k7wN+GKSbdryk4HvjF+vvV6d5PqqWtLmd0mSLdr47wZun25SVXVfkrMZjnBuA66bZNXVwMNJbgTOBT7OcFru+iQB7mH41t9M657MUuCcJKsZvrF4dOv/MPCpDH8K8AhwKu0UYVU9kuTNwIVJflRVZ02zj/XmH+pKTxyb6vc9w+dXeqI74IADasWKFXM9DUnzUJKVVXXAbLfzShKSpC71/Ie6m0yS8xm+Cj7qfVV18VzMR5JkQAFQVW+c6zlIkh7LU3ySpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLqWq5noO6kCSe4Db53oeM7AAuHeuJ7GJWOv880SpEx5b68Kq2mW2AxhQ2qwkWVFVB8z1PDYFa51/nih1wsap1VN8kqQuGVCSpC4ZUNrcLJvrCWxC1jr/PFHqhI1Qq59BSZK65BGUJKlLBpQkqUsGlLqU5E1Jbk7yaJIDRvpfnGRVe9yY5I0jy/ZPsibJPyX5RJLMzexnZ4paX5VkZatpZZLfGlk232rdOcnlSdYm+eS4beZVrW3ZB1o9307y6pH+zbLWUUlemOSqVseFSZ4ysmzCuidVVT58dPcA9gKeBywHDhjp3w7YqrV3Bf515PW1wMFAgIuA18x1HRtY64uAZ7T284Hvjyybb7VuD7wceCfwyXHbzLda9wZuBLYBdgO+C2y5Odc6ru7rgN9o7bcDH5qu7skeHkGpS1V1S1V9e4L+B6rq4fbySUABJNkVeEpVXVXDb8NfAIdtqvluiClqvaGq7movbwaelGSbeVrr/VX198BPR/vnY63AfwD+qqoerKpbgX8CXrw51zrO84ArWvtS4PDWnrDuqQYyoLTZSXJQkpuBNcA7W2D9CnDnyGp3tr754nDghqp6kPlf66j5WOuvAP9v5PVYTfOl1puAN7T2m4BntfZkdU9qq40+NWmGkvwt8MsTLDqpqr462XZVdQ2wT5K9gM8luYjhlMgvrLpxZrrh1rfWtu0+wBnAoWNdE6w2L2qdaLgJ+jb3WierqetaR01VN8NpvU8kOQW4APjZ2GYTrD9lfQaU5kxVvXIDt78lyf0Mn8/cCTxzZPEzgbsm3HAOrG+tSZ4JnA8cVVXfbd3zstZJzMda72TdUQWsq6nrWkfNoO5DAZI8F3hd65us7kl5ik+blSS7JdmqtRcynO++raruBn6c5CXtm09HAbP9v/WuJHkq8HXgA1V15Vj/fKx1MvO01guAN7fPE3cD9gCunS+1Jnl6e94COBn487ZowrqnHGyuv/Hhw8dED+CNDP/H9SDwL8DFrf+tDF8YWAVcDxw2ss0BDOe/vwt8knallN4fU9R6MnB/q3Xs8fT5WGtbdhvwb8Dats7e87jWk1o932bkm3qba63j6j4e+E57nD5aw2R1T/bwUkeSpC55ik+S1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNK6liS25IsaO1vbYTx3jb+auEbU5ITkmw38vr/tr/nmmqbn9c4Wf9Y7UkWJ/naLOZzbpIjWvszSfZu7bUzHWM2Hq9xn6gMKOlxMvYHxRtLVb10Y473ODmB4YrzAFTVa6vqvg0ddGPUXlXvqKp/2NBxtOkYUNIkkixKckuSs9t9fS5Jsm1btl+Sq5OsTnJ+kqe1/uVJPprkm8Dx7fXHklzRxjowyZeT/GOSD4/s6ysZ7vl0c5JjJ5nP2vZ8WtbdE+v7Sc5p/Ucmubb1fzrJlq3/mCTfaXN62SRj/8bImDckeXI7Wrmi1fcPSf68XR2AJH+WZEWb76mt7zjgGcDlSS5vfaNHQdPWOMXP4heOTNp7eUOS3TPcR+mbbfyLM1wZfPz6y/PYe1B9JMM9xa5O8kutb2GSy9rP9bIkz56mf7cM9z66LsmHZlOTZmCu/+rYh49eH8Ai4GFgv/b6r4EjW3s16+55cxrwp629HDhrZIzlwBmtfTzDtcd2Zbgnzp3Azm3ZTu15W4YrCYz13wYsaO214+a3Y5vH/gz3HroQ2LotO4vhUjm7AncAuwD/DriScfdbautfCLystXdguE7nYobbX+wObMlw64Qjxs13y1bjvuPnO8H8p61x3Jx+ofY2p68BLwVWAs8Gtga+BezS1vld4LOtfe7InJfT7svEcJHS17f2fwdOHnkfjm7ttwNfmab/AobrJAK8e/zPyMeGPTyCkqZ2a1Wtau2VwKIkOwJPrapvtv7PAb8+ss3/HjfGBe15DXBzVd1dw20zvse6i2cel+RG4OrWt8dUk2rXajsP+FhVrQRewRBU1yVZ1V7vDhwELK+qe6rqZxPMbcyVwJ+0o6Cn1rp7bl1bVd+rqkeALzLcVBDgd5JcD9wA7MNwM7rpzKrGKewFLGMImDsYrsf4fODSVvvJPPaiqxP5GUPQQfu5tvbBwBda+/Osq3ey/pcxvC9j/dqIvJq5NLUHR9qPMPzf/3Tun2SMR8eN9yiwVZLFwCuBg6vqgSTLGW7GOJWlwJ1VdU57HeBzVfWB0ZWSHMYMbtlQVacn+TrwWuDqJGNXqx6/bWW40Od7gQOr6gdJzp1uvutZ42Tubtu+iOGINAzBf/Asxnio2mEPw891sn8LJ3vvagbraAN5BCXNUlX9EPhBkkNa11uBb06xyXR2BH7Q/uHeE3jJVCsn+W3gVcBxI92XAUdk3ZWkd8pwtfdrgMVJdk6yNcMN5CYa8zlVtaaqzgBWAHu2RS9un7NswXDq7O+BpzCE8A/bZzevGRnqx8CTN7TGadzHcAuHj7bg+zawS5KDWy1bZ7iH1vr4FvDm1l7CUO9U/VeO69dGZEBJ6+do4Mwkq4H9GD6HWl9/w3AktRr4EMMpsKn8F4YvI4x9IeK0Gr6ddjJwSRvnUmDXGm7hsBS4CvhbhivAT+SEJDe1U3A/AS5q/VcxXJH6JuBW4PyqupHh1N7NwGcZ/pEeswy4aOxLEhtQ45Sq6l+A1wOfYjiSOgI4o81/FcNnVOvjOOCYNs+3MnxuOFX/8cC7k1zHEMLaiLyauaQJtaOT91bVb8/xVPQE5RGUJKlLHkFJkrrkEZQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpS/8fL5YZrfxFHYUAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVQUlEQVR4nO3de7QlZX3m8e/DZZCLotCYYNRuMCgXRQwgopJ0ouKo0cEFJo6NIMYQHSOQGTMqsLDBy8AwK0ZHSWwcwRg0E9eIgg4BgrQkyK0bmm4ISxPlMghJYEXUBkUuv/mj3mNvjufa3fR5+/D9rLXXfvdbVW+9v304/VC196lKVSFJUm+2mOsJSJI0EQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDStrMJHkkyaokNyX5UpLtplj3DUnev4nn97wky9scb0mybNzyk9qyVSO1rEpy3Cz2ceLGn7l6E/8OStq8JFlbVTu09nnAyqr6kzme1s8luRg4q6q+2l6/oKrWTLLuz2uZ5T7WazttXjyCkjZvfwf8apKdknwlyeokVyfZFyDJ25J8srXf1I66bkxyRevbJ8m17QhmdZI9Wv9/buvelOSE1reoHRGdneTmJJck2XaCOe0K3Dn2YrJwGpVkyyRnJrmuzeMPWv+uSa4YOWI8JMnpwLat77wNevfUNQNK2kwl2Qp4DbAGOBW4oar2BU4E/mKCTU4BXl1VLwTe0PreCXy8qvYDDgDuTLI/cAxwEPAS4PeTvKitvwfwqaraB7gPOHyC/XwM+EaSi5L8UZKnzqCc3wN+WFUHAge2fe4GvAW4uM3vhcCqqno/8JOq2q+qlsxgbG2mDChp87NtklXACuAO4H8BLwc+D1BV3wB2TrLjuO2uBM5N8vvAlq3vKuDEJO8DFlbVT9pY51fV/VW1FvgycEhb/9aqWtXaK4FF4ydXVecAewFfAhYDVyfZZpqaDgWOanVdA+zMEIbXAcckWQq8oKp+PM04mke2musJSJq1n7Qjip9LkgnWe8wHzFX1ziQHAa8DViXZr6q+kOSa1ndxkncAE4015sGR9iPARKf4qKq7gM8Cn01yE/D8JH8IvAi4q6peO26TAO+pqovHj5Xk19v8Pp/kzKqa6OhQ85BHUNL8cAWwBCDJYuDeqvrR6ApJnlNV11TVKcC9wLOS7A58r6o+AVwA7NvGOizJdkm2B97I8FnXjCT590m2bu1fZjga+n5VHdNOy40PJ4CLgXeNbPfcJNsnWQj8a1WdzXCk+Gtt/YfG1tX85RGUND8sBc5Jshp4ADh6gnXObF+CCHAZcCPwfuDIJA8B/wycVlX/luRc4Nq23Weq6oYki2Y4l0OBjyf5aXv9x1X1z9Ns8xmG04XXt6PBe4DDGE4R/nGb31rgqLb+MmB1kuv9HGr+8mvmkqQueYpPktQlA0qS1CUDSpLUJQNKktQlv8UnABYsWFCLFi2a62lImodWrlx5b1XtMtvtDCgBsGjRIlasWDHX05A0DyW5fX228xSfJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUNoqlS5fO9RQkbSKb6vfdgNLsnXceLFoEW2wxPJ93Hqeeeupcz0rSJrKpft+nDagki5LcNNMBkyxN8t4NmxYkOSHJdrNdL8mJG7rvGexzzySrktyQ5DlJ3rIRx74tyYKNMM7ajTGfX3DeeXDssXD77VA1PB97LP/xcdmZpCeyreZ6AlM4AfhL4IFZrnci8NHZ7CjJllX1yCw2OQz4alV9MMli4C3AFx7H/fXjpJPggXE/kgce4KPA4sWL52JGkuapmZ7i2zLJ2UluTnJJkm3bkcPfJFmZ5O+S7Dl+oyTLk3wsyRVJbklyYJIvJ/nHJB9u62yf5OtJbkxyU5LfTXIc8Azg8iSXt/UOTXJVkuuTfCnJDuPXS3I6sG07ujmvbXdkkmtb36eTbNn61yY5Lck1wMETFZ3klCTXtXkty+C1DKH4jja304FD2vh/lGTLJGe27VYn+YM21uI2xy8Aayaqe2TX72l1rhl7X5PslOQrbcyrk+zb+ndIck5bd3WSw8fVsKC9b6+boL5jk6xIsuKee+6Z2X8Jd9wxYfezZ7a1JM1cVU35ABYBDwP7tdd/DRwJXAbs0foOAr7R2kuB97b2cuCM1j4euAvYFdgGuBPYGTgcOHtkfzu259uABa29ALgC2L69fh9wyvj12uu1I+29gAuBrdvrs4CjWruA35mm9p1G2p8HXj9BjYuBr42sdyxwcmtvA6wAdmvr3Q/s1pZNVfd7Wvs/AZ9p7f8JfLC1fwtY1dpnAH86Ms7Txt4H4JeAa4BXTfdz3n///WtGFi6sGk7uPeZxK8xse0mbPWb5+w6sqGn+DZroMdMjqFuralVrr2QIrZcCX0qyCvg0Q/BM5IL2vAa4uarurqoHge8Bz2r9r0xyRpJDquqHE4zxEmBv4Mq2v6OBhTOY9yuA/YHr2navAHZvyx4B/s802/9mkmuSrGEIhX1msM9DgaPa/q5hCOE92rJrq+rW1p6q7i+357H3GuDlDCFJVX0D2DnJjsArgU+NbVhVP2jNrRn+J+K/VtWlM5j3zHzkI7DduI8Gt9uOx/2DP0lPODMNqAdH2o8AOwH3VdV+I4+9ptn20XHjPApsVVXfYQiRNcB/S3LKBGMEuHRkX3tX1e/NYN4BPjey3fOqamlb9tOa4nOgJE9iOOI6oqpeAJwNPGmG+3zPyD53q6pL2rL7x1aapu6x9+kR1n1OmAn2Va2/Jlj2MEPAvXoGc565JUtg2TJYuBCS4XnZMr64UXciSev/NfMfAbcmeRNA+2zmheszUJJnAA9U1V8C/wP4tbbox8CTW/tq4GVJfrVts12S506wHsBDSbZu7cuAI5I8vW23U5KZHHnBujC6N8kOwBGTrDd+/xcD7xqbQ5LnJtl+/EZT1D2ZK4AlbdvFwL1V9SPgEuAPR8Z9WmsW8HZgzyTvn2bs2VmyBG67DR59dHhesoQPfvCDG3UXkvq1qX7fN+RbfEuAP0tyMsPppL8CblyPcV4AnJnkUeAh4F2tfxlwUZK7q+o3k7wN+GKSbdryk4HvjF+vvV6d5PqqWtLmd0mSLdr47wZun25SVXVfkrMZjnBuA66bZNXVwMNJbgTOBT7OcFru+iQB7mH41t9M657MUuCcJKsZvrF4dOv/MPCpDH8K8AhwKu0UYVU9kuTNwIVJflRVZ02zj/XmH+pKTxyb6vc9w+dXeqI74IADasWKFXM9DUnzUJKVVXXAbLfzShKSpC71/Ie6m0yS8xm+Cj7qfVV18VzMR5JkQAFQVW+c6zlIkh7LU3ySpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLqWq5noO6kCSe4Db53oeM7AAuHeuJ7GJWOv880SpEx5b68Kq2mW2AxhQ2qwkWVFVB8z1PDYFa51/nih1wsap1VN8kqQuGVCSpC4ZUNrcLJvrCWxC1jr/PFHqhI1Qq59BSZK65BGUJKlLBpQkqUsGlLqU5E1Jbk7yaJIDRvpfnGRVe9yY5I0jy/ZPsibJPyX5RJLMzexnZ4paX5VkZatpZZLfGlk232rdOcnlSdYm+eS4beZVrW3ZB1o9307y6pH+zbLWUUlemOSqVseFSZ4ysmzCuidVVT58dPcA9gKeBywHDhjp3w7YqrV3Bf515PW1wMFAgIuA18x1HRtY64uAZ7T284Hvjyybb7VuD7wceCfwyXHbzLda9wZuBLYBdgO+C2y5Odc6ru7rgN9o7bcDH5qu7skeHkGpS1V1S1V9e4L+B6rq4fbySUABJNkVeEpVXVXDb8NfAIdtqvluiClqvaGq7movbwaelGSbeVrr/VX198BPR/vnY63AfwD+qqoerKpbgX8CXrw51zrO84ArWvtS4PDWnrDuqQYyoLTZSXJQkpuBNcA7W2D9CnDnyGp3tr754nDghqp6kPlf66j5WOuvAP9v5PVYTfOl1puAN7T2m4BntfZkdU9qq40+NWmGkvwt8MsTLDqpqr462XZVdQ2wT5K9gM8luYjhlMgvrLpxZrrh1rfWtu0+wBnAoWNdE6w2L2qdaLgJ+jb3WierqetaR01VN8NpvU8kOQW4APjZ2GYTrD9lfQaU5kxVvXIDt78lyf0Mn8/cCTxzZPEzgbsm3HAOrG+tSZ4JnA8cVVXfbd3zstZJzMda72TdUQWsq6nrWkfNoO5DAZI8F3hd65us7kl5ik+blSS7JdmqtRcynO++raruBn6c5CXtm09HAbP9v/WuJHkq8HXgA1V15Vj/fKx1MvO01guAN7fPE3cD9gCunS+1Jnl6e94COBn487ZowrqnHGyuv/Hhw8dED+CNDP/H9SDwL8DFrf+tDF8YWAVcDxw2ss0BDOe/vwt8knallN4fU9R6MnB/q3Xs8fT5WGtbdhvwb8Dats7e87jWk1o932bkm3qba63j6j4e+E57nD5aw2R1T/bwUkeSpC55ik+S1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNK6liS25IsaO1vbYTx3jb+auEbU5ITkmw38vr/tr/nmmqbn9c4Wf9Y7UkWJ/naLOZzbpIjWvszSfZu7bUzHWM2Hq9xn6gMKOlxMvYHxRtLVb10Y473ODmB4YrzAFTVa6vqvg0ddGPUXlXvqKp/2NBxtOkYUNIkkixKckuSs9t9fS5Jsm1btl+Sq5OsTnJ+kqe1/uVJPprkm8Dx7fXHklzRxjowyZeT/GOSD4/s6ysZ7vl0c5JjJ5nP2vZ8WtbdE+v7Sc5p/Ucmubb1fzrJlq3/mCTfaXN62SRj/8bImDckeXI7Wrmi1fcPSf68XR2AJH+WZEWb76mt7zjgGcDlSS5vfaNHQdPWOMXP4heOTNp7eUOS3TPcR+mbbfyLM1wZfPz6y/PYe1B9JMM9xa5O8kutb2GSy9rP9bIkz56mf7cM9z66LsmHZlOTZmCu/+rYh49eH8Ai4GFgv/b6r4EjW3s16+55cxrwp629HDhrZIzlwBmtfTzDtcd2Zbgnzp3Azm3ZTu15W4YrCYz13wYsaO214+a3Y5vH/gz3HroQ2LotO4vhUjm7AncAuwD/DriScfdbautfCLystXdguE7nYobbX+wObMlw64Qjxs13y1bjvuPnO8H8p61x3Jx+ofY2p68BLwVWAs8Gtga+BezS1vld4LOtfe7InJfT7svEcJHS17f2fwdOHnkfjm7ttwNfmab/AobrJAK8e/zPyMeGPTyCkqZ2a1Wtau2VwKIkOwJPrapvtv7PAb8+ss3/HjfGBe15DXBzVd1dw20zvse6i2cel+RG4OrWt8dUk2rXajsP+FhVrQRewRBU1yVZ1V7vDhwELK+qe6rqZxPMbcyVwJ+0o6Cn1rp7bl1bVd+rqkeALzLcVBDgd5JcD9wA7MNwM7rpzKrGKewFLGMImDsYrsf4fODSVvvJPPaiqxP5GUPQQfu5tvbBwBda+/Osq3ey/pcxvC9j/dqIvJq5NLUHR9qPMPzf/3Tun2SMR8eN9yiwVZLFwCuBg6vqgSTLGW7GOJWlwJ1VdU57HeBzVfWB0ZWSHMYMbtlQVacn+TrwWuDqJGNXqx6/bWW40Od7gQOr6gdJzp1uvutZ42Tubtu+iOGINAzBf/Asxnio2mEPw891sn8LJ3vvagbraAN5BCXNUlX9EPhBkkNa11uBb06xyXR2BH7Q/uHeE3jJVCsn+W3gVcBxI92XAUdk3ZWkd8pwtfdrgMVJdk6yNcMN5CYa8zlVtaaqzgBWAHu2RS9un7NswXDq7O+BpzCE8A/bZzevGRnqx8CTN7TGadzHcAuHj7bg+zawS5KDWy1bZ7iH1vr4FvDm1l7CUO9U/VeO69dGZEBJ6+do4Mwkq4H9GD6HWl9/w3AktRr4EMMpsKn8F4YvI4x9IeK0Gr6ddjJwSRvnUmDXGm7hsBS4CvhbhivAT+SEJDe1U3A/AS5q/VcxXJH6JuBW4PyqupHh1N7NwGcZ/pEeswy4aOxLEhtQ45Sq6l+A1wOfYjiSOgI4o81/FcNnVOvjOOCYNs+3MnxuOFX/8cC7k1zHEMLaiLyauaQJtaOT91bVb8/xVPQE5RGUJKlLHkFJkrrkEZQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpS/8fL5YZrfxFHYUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] @@ -987,6 +1017,10 @@ "metadata": {}, "source": [ "### References\n", + "\n", + "Bayona, J.A., Savran, W.H., Rhoades, D.A. and Werner, M.J., 2022. Prospective evaluation of multiplicative hybrid earthquake forecasting models in California. Geophysical Journal International, 229(3), pp.1736-1753.\n", + "doi: https://doi.org/10.1093/gji/ggac018\n", + "\n", "Field, E. H., K. R. Milner, J. L. Hardebeck, M. T. Page, N. J. van der Elst, T. H. Jordan, A. J. Michael, B. E. Shaw, and M. J. Werner (2017). A spatiotemporal clustering model for the third Uniform California Earthquake Rupture Forecast (UCERF3-ETAS): Toward an operational earthquake forecast, Bull. Seismol. Soc. Am. 107, 1049–1081.\n", "\n", "Harte, D., and D. Vere-Jones (2005), The entropy score and its uses in earthquake forecasting, Pure Appl. Geophys. 162 , 6-7, 1229-1253, DOI: 10.1007/\n", @@ -1000,6 +1034,16 @@ "Imoto, M., and D.A. Rhoades (2010), Seismicity models of moderate earthquakes in Kanto, Japan utilizing multiple predictive parameters, Pure Appl. Geophys.\n", "167, 6-7, 831-843, DOI: 10.1007/s00024-010-0066-4.\n", "\n", + "Kato, M., 2019. On the apparently inappropriate use of multiple hypothesis testing in earthquake prediction studies. Seismological Research Letters, 90(3), pp.1330-1334. doi: \n", + "https://doi.org/10.1785/0220180378\n", + "\n", + "\n", + "Lombardi, A.M. and Marzocchi, W., 2010. The assumption of Poisson seismic-rate variability in CSEP/RELM experiments. Bulletin of the Seismological Society of America, 100(5A), pp.2293-2300, doi: \n", + "https://doi.org/10.1785/0120100012\n", + "\n", + "Nandan, S., Ouillon, G., Sornette, D. and Wiemer, S., 2019. Forecasting the full distribution of earthquake numbers is fair, robust, and better. Seismological Research Letters, 90(4), pp.1650-1659. doi: \n", + "https://doi.org/10.1785/0220180374\n", + "\n", "Rhoades, D.A, D., Schorlemmer, M.C.Gerstenberger, A. Christophersen, J. D. Zechar & M. Imoto (2011) Efficient testing of earthquake forecasting models, Acta Geophysica 59\n", "\n", "Savran, W., M. J. Werner, W. Marzocchi, D. Rhoades, D. D. Jackson, K. R. Milner, E. H. Field, and A. J. Michael (2020). Pseudoprospective evaluation of UCERF3-ETAS forecasts during the 2019 Ridgecrest Sequence, Bulletin of the Seismological Society of America.\n", @@ -1014,8 +1058,12 @@ "\n", "M. Taroni, W. Marzocchi, D. Schorlemmer, M. J. Werner, S. Wiemer, J. D. Zechar, L. Heiniger, F. Euchner; Prospective CSEP Evaluation of 1‐Day, 3‐Month, and 5‐Yr Earthquake Forecasts for Italy. Seismological Research Letters 2018;; 89 (4): 1251–1261. doi: https://doi.org/10.1785/0220180031\n", "\n", + "\n", "Werner, M. J., A. Helmstetter, D. D. Jackson, and Y. Y. Kagan (2011a). High-Resolution Long-Term and Short-Term Earthquake Forecasts for California, Bulletin of the Seismological Society of America 101 1630-1648\n", "\n", + "Werner, M.J. and Sornette, D., 2008. Magnitude uncertainties impact seismic rate estimates, forecasts, and predictability experiments. Journal of Geophysical Research: Solid Earth, 113(B8). doi: https://doi.org/10.1029/2007JB005427\n", + "\n", + "\n", "Werner, M.J. J.D. Zechar, W. Marzocchi, and S. Wiemer (2011b), Retrospective evaluation of the five-year and ten-year CSEP-Italy earthquake forecasts, Annals of Geophysics 53, no. 3, 11–30, doi:10.4401/ag-4840.\n", "\n", "Zechar, 2011: Evaluating earthquake predictions and earthquake forecasts: a guide for students and new researchers, CORSSA (http://www.corssa.org/en/articles/theme_6/)\n", @@ -1025,6 +1073,14 @@ "Zechar, J.D., D. Schorlemmer, M. Liukis, J. Yu, F. Euchner, P.J. Maechling, and T.H. Jordan (2010b), The Collaboratory for the Study of Earthquake Predictability perspective on computational earthquake science, Concurr. Comp-Pract. E., doi:10.1002/cpe.1519.\n", "\n" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fifty-wright", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -1043,7 +1099,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.8" + "version": "3.8.6" } }, "nbformat": 4, diff --git a/requirements.txt b/requirements.txt index 78dfc73b..6a95af5b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy<=1.21.5 +numpy scipy pandas matplotlib diff --git a/requirements.yml b/requirements.yml index 85779ee0..14caee7e 100644 --- a/requirements.yml +++ b/requirements.yml @@ -4,7 +4,7 @@ channels: - defaults dependencies: - python>=3.7 - - numpy<=1.21.5 + - numpy - pandas - scipy - matplotlib diff --git a/setup.py b/setup.py index f11037a0..5390de9f 100644 --- a/setup.py +++ b/setup.py @@ -12,6 +12,7 @@ def get_version(): raise RuntimeError("Unable to find version string in %s." % (VERSIONFILE,)) return verstr + with open("README.md",'r') as fh: long_description = fh.read() @@ -25,7 +26,7 @@ def get_version(): description='Python tools from the Collaboratory for the Study of Earthquake Predictability', long_description=long_description, long_description_content_type='text/markdown', - install_requires = [ + install_requires=[ 'numpy', 'scipy', 'pandas', diff --git a/tests/artifacts/BSI/vcr_search.yaml b/tests/artifacts/BSI/vcr_search.yaml new file mode 100644 index 00000000..76053c90 --- /dev/null +++ b/tests/artifacts/BSI/vcr_search.yaml @@ -0,0 +1,46 @@ +interactions: +- request: + body: null + headers: + Connection: + - close + Host: + - webservices.rm.ingv.it + User-Agent: + - Python-urllib/3.10 + method: GET + uri: https://webservices.rm.ingv.it/fdsnws/event/1/query?format=geojson&starttime=2009-04-06T00%3A00%3A00&endtime=2009-04-07T00%3A00%3A00&limit=15000&maxdepth=1000&maxmagnitude=10.0&mindepth=-100&minmagnitude=5.5&offset=0&orderby=time-asc&eventtype=earthquake + response: + body: + string: '{"type":"FeatureCollection","features":[{"type":"Feature","properties":{"eventId":1895389,"originId":755599,"time":"2009-04-06T01:32:40.400000","author":"BULLETIN-SISPICK","magType":"Mw","mag":6.1,"magAuthor":"--","type":"earthquake","place":"2 + km SW L''Aquila (AQ)","version":1000,"geojson_creationTime":"2022-09-21T12:18:03"},"geometry":{"type":"Point","coordinates":[13.38,42.342,8.3]}}]}' + headers: + Access-Control-Allow-Origin: + - '*' + Cache-Control: + - public, max-age=60 + Connection: + - close + Content-Type: + - application/json + Date: + - Wed, 21 Sep 2022 12:18:03 GMT + Server: + - nginx + Transfer-Encoding: + - chunked + Vary: + - Accept-Encoding + - Accept-Encoding + X-Cache-Status: + - EXPIRED + X-RateLimit-Limit: + - '10' + X-RateLimit-Reset: + - '1' + X-UA-Compatible: + - IE=Edge + status: + code: 200 + message: OK +version: 1 diff --git a/tests/artifacts/BSI/vcr_summary.yaml b/tests/artifacts/BSI/vcr_summary.yaml new file mode 100644 index 00000000..9349aa3d --- /dev/null +++ b/tests/artifacts/BSI/vcr_summary.yaml @@ -0,0 +1,46 @@ +interactions: +- request: + body: null + headers: + Connection: + - close + Host: + - webservices.rm.ingv.it + User-Agent: + - Python-urllib/3.10 + method: GET + uri: https://webservices.rm.ingv.it/fdsnws/event/1/query?format=geojson&starttime=2009-04-06T00%3A00%3A00&endtime=2009-04-07T00%3A00%3A00&limit=15000&maxdepth=1000&maxmagnitude=10.0&mindepth=-100&minmagnitude=5.5&offset=0&orderby=time-asc&eventtype=earthquake + response: + body: + string: '{"type":"FeatureCollection","features":[{"type":"Feature","properties":{"eventId":1895389,"originId":755599,"time":"2009-04-06T01:32:40.400000","author":"BULLETIN-SISPICK","magType":"Mw","mag":6.1,"magAuthor":"--","type":"earthquake","place":"2 + km SW L''Aquila (AQ)","version":1000,"geojson_creationTime":"2022-09-21T12:14:55"},"geometry":{"type":"Point","coordinates":[13.38,42.342,8.3]}}]}' + headers: + Access-Control-Allow-Origin: + - '*' + Cache-Control: + - public, max-age=60 + Connection: + - close + Content-Type: + - application/json + Date: + - Wed, 21 Sep 2022 12:14:55 GMT + Server: + - nginx + Transfer-Encoding: + - chunked + Vary: + - Accept-Encoding + - Accept-Encoding + X-Cache-Status: + - MISS + X-RateLimit-Limit: + - '10' + X-RateLimit-Reset: + - '1' + X-UA-Compatible: + - IE=Edge + status: + code: 200 + message: OK +version: 1 diff --git a/tests/test_bsi.py b/tests/test_bsi.py new file mode 100644 index 00000000..4673bfc8 --- /dev/null +++ b/tests/test_bsi.py @@ -0,0 +1,42 @@ +from datetime import datetime +import os.path +import vcr +from csep.utils.comcat import search + +HOST = 'webservices.rm.ingv.it' + + +def get_datadir(): + root_dir = os.path.dirname(os.path.abspath(__file__)) + data_dir = os.path.join(root_dir, 'artifacts', 'BSI') + return data_dir + + +def test_search(): + datadir = get_datadir() + tape_file = os.path.join(datadir, 'vcr_search.yaml') + with vcr.use_cassette(tape_file): + # L'Aquila + eventlist = search(starttime=datetime(2009, 4, 6, 0, 0, 0), + endtime=datetime(2009, 4, 7, 0, 0, 0), + minmagnitude=5.5, host=HOST, limit=15000, offset=0) + event = eventlist[0] + assert event.id == 1895389 + + +def test_summary(): + datadir = get_datadir() + tape_file = os.path.join(datadir, 'vcr_summary.yaml') + with vcr.use_cassette(tape_file): + eventlist = search(starttime=datetime(2009, 4, 6, 0, 0, 0), + endtime=datetime(2009, 4, 7, 0, 0, 0), + minmagnitude=5.5, host=HOST, limit=15000, offset=0) + event = eventlist[0] + cmp = '1895389 2009-04-06 01:32:40.400000 (42.342,13.380) 8.3 km M6.1' + assert str(event) == cmp + assert event.id == 1895389 + assert event.time == datetime(2009, 4, 6, 1, 32, 40, 400000) + assert event.latitude == 42.342 + assert event.longitude == 13.380 + assert event.depth == 8.3 + assert event.magnitude == 6.1 diff --git a/tests/test_csep1_evaluations.py b/tests/test_csep1_evaluations.py index a733003f..d095e677 100644 --- a/tests/test_csep1_evaluations.py +++ b/tests/test_csep1_evaluations.py @@ -21,6 +21,7 @@ def get_datadir(): data_dir = os.path.join(root_dir, 'artifacts', 'example_csep1_forecasts') return data_dir + class TestCSEP1NTestThreeMonthsEEPAS(unittest.TestCase): def __init__(self, *args, **kwargs): @@ -61,8 +62,8 @@ def test_ntest_three_months_eepas_model(self): test_evaluation_dict['event_count_forecast'] = fore.event_count test_evaluation_dict['event_count'] = cata.event_count # comparing floats, so we will just map to ndarray and use allclose - expected = numpy.array([v for k,v in result_dict.items()]) - computed = numpy.array([v for k,v in test_evaluation_dict.items()]) + expected = numpy.array([v for k, v in result_dict.items()]) + computed = numpy.array([v for k, v in test_evaluation_dict.items()]) numpy.testing.assert_allclose(expected, computed, rtol=1e-5) @@ -82,13 +83,13 @@ def _parse_xml_result(self): xml_result['event_count'] = float(child.find('ns0:eventCount', ns).text) return xml_result -class TestGriddedForecastTests(unittest.TestCase): +class TestGriddedForecastTests(unittest.TestCase): def test_n_test(self): - forecast = numpy.zeros((10,10))+0.0015 + forecast = numpy.zeros((10, 10))+0.0015 forecast = forecast / forecast.size - observation = numpy.zeros((10,10)) + observation = numpy.zeros((10, 10)) expected_output = (1.0, 0.9985011244377109) print('N Test: Running Unit Test') numpy.testing.assert_allclose(_number_test_ndarray(forecast.sum(), observation.sum()), expected_output) @@ -116,18 +117,18 @@ def test_t_test(self): the equations from """ - forecast_A = numpy.array([[8, 2], [3, 5]]) - forecast_B = numpy.array([[6, 4], [2, 8]]) + forecast_a = numpy.array([[8, 2], [3, 5]]) + forecast_b = numpy.array([[6, 4], [2, 8]]) obs = numpy.array([[5, 8], [4, 2]]) - t_test_expected = {'t_critical': 2.10092204024096, - 't_statistic': 1.5385261717159382, + t_test_expected = {'t_statistic': 1.5385261717159382, + 't_critical': 2.10092204024096, 'information_gain': 0.08052612477654024, 'ig_lower': -0.029435677283374914, 'ig_upper': 0.19048792683645538} - - print('T Test: Running Unit Test') - self.assertEqual(_t_test_ndarray(forecast_A, forecast_B, numpy.sum(obs), forecast_A.sum(), forecast_B.sum()), t_test_expected, 'Failed T Test') + numpy.testing.assert_allclose( + [v for k, v in _t_test_ndarray(forecast_a, forecast_b, numpy.sum(obs), forecast_a.sum(), forecast_b.sum()).items()], + [v for k, v in t_test_expected.items()]) if __name__ == '__main__': diff --git a/tests/test_evaluations.py b/tests/test_evaluations.py index df6e22e6..c7468528 100644 --- a/tests/test_evaluations.py +++ b/tests/test_evaluations.py @@ -2,13 +2,15 @@ import numpy import unittest -from csep.core.poisson_evaluations import _simulate_catalog, _poisson_likelihood_test +import csep.core.poisson_evaluations as poisson +import csep.core.binomial_evaluations as binary def get_datadir(): root_dir = os.path.dirname(os.path.abspath(__file__)) data_dir = os.path.join(root_dir, 'artifacts', 'Comcat') return data_dir + class TestPoissonLikelihood(unittest.TestCase): def __init__(self, *args, **kwargs): @@ -22,7 +24,7 @@ def __init__(self, *args, **kwargs): def test_simulate_catalog(self): # expecting the sampling weights to be [0.25, 0.5, 0.75, 1.0] - # assuming the random numbers are equal to thhe following: + # assuming the random numbers are equal to the following: random_numbers = numpy.array([[0.5488135, 0.71518937, 0.60276338, 0.54488318]]) num_events = 4 @@ -45,21 +47,21 @@ def test_simulate_catalog(self): # this is taken from the test likelihood function sim_fore = numpy.empty(sampling_weights.shape) - sim_fore = _simulate_catalog(num_events, sampling_weights, sim_fore, + sim_fore = poisson._simulate_catalog(num_events, sampling_weights, sim_fore, random_numbers=self.random_matrix) # final statement numpy.testing.assert_allclose(expected_catalog, sim_fore) # test again to ensure that fill works properply - sim_fore = _simulate_catalog(num_events, sampling_weights, sim_fore, + sim_fore = poisson._simulate_catalog(num_events, sampling_weights, sim_fore, random_numbers=self.random_matrix) # final statement numpy.testing.assert_allclose(expected_catalog, sim_fore) def test_likelihood(self): - qs, obs_ll, simulated_ll = _poisson_likelihood_test(self.forecast_data, self.observed_data, num_simulations=1, + qs, obs_ll, simulated_ll = poisson._poisson_likelihood_test(self.forecast_data, self.observed_data, num_simulations=1, random_numbers=self.random_matrix, use_observed_counts=True) # very basic result to pass "laugh" test @@ -71,3 +73,47 @@ def test_likelihood(self): # calculated by hand given the expected data, see explanation in zechar et al., 2010. numpy.testing.assert_allclose(simulated_ll[0], -7.178053830347945) + +class TestBinomialLikelihood(unittest.TestCase): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.seed = 0 + numpy.random.seed(self.seed) + self.forecast_data = numpy.array([[0.1, 0.3, 0.4], [0.2, 0.1, 0.1]]) + self.observed_data = numpy.array([[0, 1, 2], [1, 1, 0]]) + self.random_matrix = numpy.random.rand(1, 9) + + def test_joint_likelihood_calculation(self): + bill = binary.binary_joint_log_likelihood_ndarray(self.forecast_data, self.observed_data) + numpy.testing.assert_allclose(bill, -6.7197988064) + + def test_simulate_active_cells(self): + #With fixed seed we get the same random numbers if we get all the number at once or one by one. + #Making sure random number generated by seed 0 match. + expected_random_numbers = numpy.array([[0.5488135, 0.71518937, 0.60276338, 0.54488318, 0.4236548, 0.64589411, + 0.4375872112626925, 0.8917730007820798, 0.9636627605010293]]) + + numpy.testing.assert_allclose(expected_random_numbers, self.random_matrix) + + #We can expect the following catalog, if we get the above random numbers. + #We get 4 active cells after 9th random sample. + expected_catalog = [0, 0, 1, 1, 1, 1] + + sampling_weights = numpy.cumsum(self.forecast_data.ravel()) / numpy.sum(self.forecast_data) + sim_fore = numpy.zeros(sampling_weights.shape) + obs_active_cells = len(numpy.unique(numpy.nonzero(self.observed_data.ravel()))) + #resetting seed again to 0, to make sure _simulate_catalog uses this. + seed = 0 + numpy.random.seed(seed) + sim_fore = binary._simulate_catalog(obs_active_cells, sampling_weights, sim_fore) + numpy.testing.assert_allclose(expected_catalog, sim_fore) + + def test_binomial_likelihood(self): + qs, bill, simulated_ll = binary._binary_likelihood_test(self.forecast_data,self.observed_data, num_simulations=1,seed=0, verbose=True) + numpy.testing.assert_allclose(bill, -6.7197988064) + numpy.testing.assert_allclose(qs, 1) + numpy.testing.assert_allclose(simulated_ll[0], -7.921741654647629) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_forecast.py b/tests/test_forecast.py index f63d8a3e..7b427586 100644 --- a/tests/test_forecast.py +++ b/tests/test_forecast.py @@ -1,14 +1,24 @@ -import os, unittest +import os +import unittest import numpy from csep import load_catalog_forecast + def get_test_catalog_root(): root_dir = os.path.dirname(os.path.abspath(__file__)) data_dir = os.path.join(root_dir, 'artifacts', 'test_ascii_catalogs') return data_dir + class TestCatalogForecastCreation(unittest.TestCase): + def test_all_present(self): + fname = os.path.join(get_test_catalog_root(), 'all_present.csv') + test_fore = load_catalog_forecast(fname) + total_event_count = numpy.array([cat.event_count for cat in test_fore]).sum() + self.assertEqual(10, test_fore.n_cat) + self.assertEqual(10, total_event_count) + def test_ascii_load_all_empty(self): fname = os.path.join(get_test_catalog_root(), 'all_empty.csv') test_fore = load_catalog_forecast(fname) @@ -56,5 +66,12 @@ def test_get_event_counts(self): test_fore = load_catalog_forecast(fname) numpy.testing.assert_array_equal(numpy.ones(10), test_fore.get_event_counts()) + def test_multiple_iterations(self): + fname = os.path.join(get_test_catalog_root(), 'all_present.csv') + test_fore = load_catalog_forecast(fname) + ec1 = [cat.event_count for cat in test_fore] + ec2 = [cat.event_count for cat in test_fore] + numpy.testing.assert_array_equal(ec1, ec2) + if __name__ == '__main__': unittest.main()