diff --git a/docs/guide/experiment_config.rst b/docs/guide/experiment_config.rst index ba0b7c2..c1a535d 100644 --- a/docs/guide/experiment_config.rst +++ b/docs/guide/experiment_config.rst @@ -14,9 +14,10 @@ Configuration files are written in ``YAML`` format and are divided into differen 2. **Temporal Configuration** (``time_config``): Temporal characteristics of the experiment, such as the start and end dates, experiment class (time-independent or time-dependent), the testing intervals, etc. 3. **Spatial and Magnitude Configuration** (``region_config``): Describes the testing region, such as its geographic bounds, magnitude ranges, and depth ranges. 4. **Seismic Catalog** (``catalog``): Defines the seismicity data source to test the models. It can either link to a seismic network API, or an existing file in the system. -5. **Models** (``models``): Configuration of forecasting models. It can direct to an additional configuration file with ``model_config`` for readability. See :ref:`model_config`. -6. **Evaluation Tests** (``tests``): Configuration of the statistical tests to evaluate the models. It can direct to an additional configuration file with ``test_config`` for readability. See :ref:`evaluation_config`. -7. **Postprocessing** (``postprocess``): Instructions on how to process and visualize the experiment's results, such as plotting forecasts or generating reports. See :ref:`postprocess`. +5. **Input Data** (``input_config``): Describes the data handed to the models (e.g., the input catalog of time-dependent models), when it should differ from the testing settings. +6. **Models** (``models``): Configuration of forecasting models. It can direct to an additional configuration file with ``model_config`` for readability. See :ref:`model_config`. +7. **Evaluation Tests** (``tests``): Configuration of the statistical tests to evaluate the models. It can direct to an additional configuration file with ``test_config`` for readability. See :ref:`evaluation_config`. +8. **Postprocessing** (``postprocess``): Instructions on how to process and visualize the experiment's results, such as plotting forecasts or generating reports. See :ref:`postprocess`. .. note:: @@ -228,7 +229,7 @@ Configuring the spatial and magnitude definitions is done through the ``region_c resulting in the ``[4.0, 4.1)``, ``[4.1, 4.2)`` and ``[4.2, 4.3)``. -3. Depths: The minimum and maximum depths are just required to filter out seismicity outside those ranges. +3. Depths: The minimum and maximum depths are just required to filter out seismicity outside those ranges. The testing catalogs contain the events with ``depth_min <= depth < depth_max``. Some example of region configurations would be: @@ -282,6 +283,47 @@ The seismicity catalog can be defined with the ``catalog`` parameter. It represe The main catalog will be stored, and consecutively filtered to the extent of each testing time-window, as well as to the experiment's spatial domain, and magnitude- and depth- ranges. +Input Data +---------- + +Time-dependent models receive, at each time window, an **input catalog** with the seismicity that occurred *before* the forecast start (see :ref:`model_config`). By default, this catalog is filtered from the main catalog using the testing magnitude range only (``mag_min``, ``mag_max``), without spatial, depth or lower time bounds. The ``input_config`` option allows to set the input data independently from the testing settings, as is usual in forecasting experiments (e.g., models are calibrated with small events in a wider region than the one they are tested in). + +.. code-block:: yaml + + input_config: + catalog: + region: italy_csep_collection_region + mag_min: 2.5 + mag_max: 8.0 + depth_min: 0 + depth_max: 30 + start_date: 2005-01-01T00:00:00 + +The ``catalog`` block admits the following parameters, all optional: + +.. list-table:: + :widths: 20 80 + + * - **region** + - Spatial domain of the input catalog, with the same options as the testing ``region`` (a **CSEP** region or a file). Default: no spatial filter. + * - **mag_min** + - Minimum magnitude of the input catalog. Default: the testing ``mag_min``. + * - **mag_max** + - Maximum magnitude of the input catalog (exclusive). Default: the testing ``mag_max``. + * - **depth_min** + - Minimum depth (in `km`) of the input catalog. Default: no filter. + * - **depth_max** + - Maximum depth (in `km`, exclusive) of the input catalog. Default: no filter. + * - **start_date** + - Earliest event (``datetime``) of the input catalog. Default: no lower bound. The upper bound is always the start of the forecast time window. + +.. note:: + + When the main catalog is queried from a network API, the query bounds are extended to cover both the testing and the input settings (i.e., the lowest magnitude, the largest depth, the earliest date and the union of the regions' bounding boxes). + +The input data configuration is written in the reproducibility configuration file, so an experiment reproduced with ``floatcsep reproduce`` hands the same input catalogs to the models. The ``input_config`` block is meant to be extended with other types of input data in the future. + + Run Configuration ----------------- diff --git a/docs/guide/model_config.rst b/docs/guide/model_config.rst index 9cad164..5d1a54c 100644 --- a/docs/guide/model_config.rst +++ b/docs/guide/model_config.rst @@ -344,17 +344,18 @@ We distinguish **input data** versus **input arguments**. The input data is give Therefore, the model source-code should be at least able to dynamically read the obligatory arguments (simply the time window of the issued forecast) -2. **Input Data**: Correspond to any data source outside the control of the modeler (e.g., authoritative input catalog, testing region). For now, **floatCSEP** just handles an input **catalog**, which are all the events within the **main catalog** until the forecast **start_date**. The catalog is written by default in ``model_path/input/catalog.csv`` in the CSEP ascii format (see :doc:`pycsep:concepts/catalogs`) as: +2. **Input Data**: Correspond to any data source outside the control of the modeler (e.g., authoritative input catalog, testing region). For now, **floatCSEP** just handles an input **catalog**, which are the events of the **main catalog** that occurred before the forecast **start_date**. By default, the input catalog contains the events within the testing magnitude range, without spatial or depth filters. The ``input_config`` option of the experiment configuration (see :ref:`experiment_config`) sets the region, magnitude range, depth range and earliest date of the input catalog, independently from the testing settings. The catalog is written by default in ``model_path/input/catalog.csv`` in the CSEP ascii format (see :doc:`pycsep:concepts/catalogs`), with a header line, as: .. code-block:: none - longitude, latitude, magnitude, time_string, depth, event_id + lon,lat,mag,time_string,depth,catalog_id,event_id - - **longitude**: Decimal degrees of the forecasted event location. - - **latitude**: Decimal degrees of the forecasted event location. - - **magnitude**: Magnitude of the forecasted event. - - **time_string**: Timestamp in UTC following the ISO8601 format (`%Y-%m-%dT%H:%M:%S`). + - **lon**: Decimal degrees of the event longitude. + - **lat**: Decimal degrees of the event latitude. + - **mag**: Magnitude of the event. + - **time_string**: Timestamp in UTC following the ISO8601 format (`%Y-%m-%dT%H:%M:%S`, with fractional seconds `.%f` when the event has them). - **depth**: Depth of the event in kilometers. + - **catalog_id**: Always ``-1`` for an observed catalog (as opposed to the index of a synthetic catalog in a forecast). - **event_id**: The event ID in case is necessary to map the event to an additional table. diff --git a/docs/reference/api_reference.rst b/docs/reference/api_reference.rst index e85de58..2c44e54 100644 --- a/docs/reference/api_reference.rst +++ b/docs/reference/api_reference.rst @@ -125,6 +125,8 @@ These are the helper functions of ``floatCSEP`` parse_timedelta_string read_time_cfg read_region_cfg + read_input_cfg + parse_region time_windows_ti time_windows_td @@ -241,8 +243,11 @@ catalogs, forecasts), abstracting the experiment logic from the pyCSEP io functi CatalogRepository CatalogRepository.set_main_catalog CatalogRepository.catalog + CatalogRepository.input_cat_config + CatalogRepository.query_bounds CatalogRepository.get_test_cat CatalogRepository.set_test_cats + CatalogRepository.get_input_cat CatalogRepository.set_input_cats CatalogRepository.filter_catalog diff --git a/docs/tutorials/case_h.rst b/docs/tutorials/case_h.rst index e2fb02a..c74de91 100644 --- a/docs/tutorials/case_h.rst +++ b/docs/tutorials/case_h.rst @@ -81,7 +81,7 @@ As in :ref:`Tutorial G`, each **Model** requires to build and execute a .. literalinclude:: ../../tutorials/case_h/models.yml :caption: tutorials/case_h/models.yml :language: yaml - :lines: 1-3, 12-13, 22-23 + :lines: 1-3, 12-13, 23-24 2. A ``path`` needs to be indicated for each model, to both download the repository contents therein and from where the source code will be executed. @@ -132,7 +132,7 @@ As in :ref:`Tutorial G`, each **Model** requires to build and execute a .. literalinclude:: ../../tutorials/case_h/models.yml :caption: tutorials/case_h/models.yml :language: yaml - :lines: 1,7,12,15,22,25 + :lines: 1,7,12,15,23,26 .. important:: Please refer to :ref:`Tutorial G` for example of how to set up ``func`` for the model and interface it to **floatCSEP**. @@ -142,7 +142,7 @@ As in :ref:`Tutorial G`, each **Model** requires to build and execute a .. literalinclude:: ../../tutorials/case_h/models.yml :caption: tutorials/case_h/models.yml :language: yaml - :lines: 12,16,22, 26 + :lines: 12,16,23,27 The experiment will read the forecasts as: @@ -157,7 +157,7 @@ As in :ref:`Tutorial G`, each **Model** requires to build and execute a .. literalinclude:: ../../tutorials/case_h/models.yml :caption: tutorials/case_h/models.yml :language: yaml - :lines: 12,18-21,22,28-33 + :lines: 12,18-22,23,29-35 .. note:: @@ -178,6 +178,18 @@ Catalog The catalog was obtained *prior* to the experiment using :func:`query bsi `, but it was filtered from 2006 onwards, so it has enough data for the model calibration. +Input Catalog +~~~~~~~~~~~~~ + + The models are time-dependent, so they receive an **input catalog** at each time window, with the seismicity before the forecast start. The ``input_config`` block sets how this catalog is filtered, independently from the testing settings of ``region_config``: here the models are calibrated with events down to magnitude 2.5 within the Italy *collection* region, whereas the forecasts are tested for magnitudes above 3.5 within the Italy *testing* region. + + .. literalinclude:: ../../tutorials/case_h/config.yml + :caption: tutorials/case_h/config.yml + :language: yaml + :lines: 17-22 + + Without ``input_config``, the input catalog would contain the events of the testing magnitude range only, with no spatial or depth filter. See :ref:`experiment_config` for the admissible parameters. + Tests ~~~~~ @@ -201,7 +213,7 @@ Custom Post-Process .. literalinclude:: ../../tutorials/case_h/config.yml :caption: tutorials/case_h/config.yml :language: yaml - :lines: 22-23 + :lines: 29-30 This option provides `hook` for a Python script and a function within as: diff --git a/floatcsep/experiment.py b/floatcsep/experiment.py index ce852a5..929ce73 100644 --- a/floatcsep/experiment.py +++ b/floatcsep/experiment.py @@ -22,6 +22,7 @@ NoAliasLoader, read_time_cfg, read_region_cfg, + read_input_cfg, timewindow2str, parse_nested_dicts, ) @@ -75,6 +76,21 @@ class Experiment: - depth_min (:class:`float`): Minimum depth. Defaults to -2 - depth_max (:class:`float`): Maximum depth. Defaults to 6000 + input_config (dict): Settings of the data handed to the models, as opposed to the + testing settings above. Contains a ``catalog`` block with the filters applied to + the main catalog before it is written as a model's input catalog: + + - region: Spatial domain of the input catalog (same options as the testing + ``region``). Defaults to no spatial filter. + - mag_min, mag_max (:class:`float`): Magnitude range of the input catalog. + Default to the testing magnitude range. + - depth_min, depth_max (:class:`float`): Depth range of the input catalog. + Default to no depth filter. + - start_date (:class:`datetime.datetime`): Earliest event of the input catalog. + Defaults to no lower bound. The upper bound is always the forecast start. + + See :func:`~floatcsep.utils.helpers.read_input_cfg`. + model_config (str): Path to the models' configuration file test_config (str): Path to the evaluations' configuration file run_mode (str): 'sequential' or 'parallel' @@ -106,6 +122,7 @@ def __init__( LICENSE: str = None, time_config: dict = None, region_config: dict = None, + input_config: dict = None, catalog: str = None, models: str = None, tests: str = None, @@ -156,6 +173,7 @@ def __init__( self.seed = kwargs.get("seed", None) self.time_config = read_time_cfg(time_config, **kwargs) self.region_config = read_region_cfg(region_config, **kwargs) + self.input_config = read_input_cfg(input_config, self.region_config, **kwargs) self.model_config = models if isinstance(models, str) else None self.test_config = tests if isinstance(tests, str) else None @@ -182,6 +200,15 @@ def __init__( else "Time-Independent" ) log.info(f"\tExperiment class: {exp_class_str}") + if self.input_config: + inp = self.input_config["catalog"] + inp_region = inp["region"].name if inp["region"] is not None else None + log.info( + f"\tInput catalog: region {inp_region}, " + f"magnitudes [{inp['mag_min']}, {inp['mag_max']}), " + f"depths [{inp.get('depth_min')}, {inp.get('depth_max')}), " + f"from {inp.get('start_date')}" + ) self.catalog = None self.models = [] @@ -189,7 +216,9 @@ def __init__( self.postprocess = postprocess if postprocess else {} self.default_test_kwargs = default_test_kwargs - self.catalog_repo.set_main_catalog(catalog, self.time_config, self.region_config) + self.catalog_repo.set_main_catalog( + catalog, self.time_config, self.region_config, self.input_config + ) self.models = self.set_models( models or kwargs.get("model_config"), kwargs.get("order", None) @@ -582,13 +611,17 @@ def make_repr(self) -> None: repr_config = self.registry.get_attr("repr_config") # Dropping region to results folder if it is a file - region_path = self.region_config.get("path", False) - if isinstance(region_path, str): - if isfile(region_path) and region_path: - new_path = join(self.registry.run_dir, self.region_config["path"]) - shutil.copy2(region_path, new_path) - self.region_config.pop("path") - self.region_config["region"] = self.registry.rel(new_path) + configs = [self.region_config] + if self.input_config: + configs.append(self.input_config["catalog"]) + for cfg in configs: + region_path = cfg.get("path", False) + if isinstance(region_path, str): + if isfile(region_path) and region_path: + new_path = join(self.registry.run_dir, cfg["path"]) + shutil.copy2(region_path, new_path) + cfg.pop("path") + cfg["region"] = self.registry.rel(new_path) # Dropping catalog to results folder target_cat = join( @@ -632,10 +665,20 @@ def as_dict(self, extra: Sequence = (), extended=False) -> dict: for i, j in self.region_config.items() if (i not in ("magnitudes", "depths") or extended) }, - "catalog": self.registry.rel(self.catalog_repo.cat_path).as_posix(), - "models": [i.as_dict() for i in self.models], - "tests": [i.as_dict() for i in self.tests], } + if self.input_config: + dict_walk["input_config"] = { + "catalog": { + i: j for i, j in self.input_config["catalog"].items() if j is not None + } + } + dict_walk.update( + { + "catalog": self.registry.rel(self.catalog_repo.cat_path).as_posix(), + "models": [i.as_dict() for i in self.models], + "tests": [i.as_dict() for i in self.tests], + } + ) dict_walk.update(extra) return parse_nested_dicts(dict_walk) diff --git a/floatcsep/infrastructure/repositories.py b/floatcsep/infrastructure/repositories.py index 2f84832..a5db9ae 100644 --- a/floatcsep/infrastructure/repositories.py +++ b/floatcsep/infrastructure/repositories.py @@ -55,6 +55,7 @@ def __init__(self, registry: ExperimentFileRegistry): self.registry = registry self.time_config = {} self.region_config = {} + self.input_config = {} def __dir__(self): """Adds time and region configs keys to instance scope.""" @@ -87,20 +88,45 @@ def as_dict(self): return def set_main_catalog( - self, catalog: Union[str, Callable, CSEPCatalog], time_config: dict, region_config: dict + self, + catalog: Union[str, Callable, CSEPCatalog], + time_config: dict, + region_config: dict, + input_config: dict = None, ): """ Sets the catalog to be used for the experiment. Args: catalog: Experiment's main catalog. - region_config: Experiment instantiation - time_config: + time_config: Experiment temporal configuration + region_config: Experiment region configuration (testing domain) + input_config: Input data configuration (see + :func:`~floatcsep.utils.helpers.read_input_cfg`). Its ``catalog`` block + sets how the models' input catalogs are filtered. """ self.time_config = time_config self.region_config = region_config + self.input_config = input_config or {} self.catalog = catalog + @property + def input_cat_config(self) -> dict: + """ + Returns: + The settings used to filter the models' input catalogs. Defaults to the testing + magnitude range, without spatial, depth or lower time bounds. + """ + cfg = self.input_config.get("catalog") or {} + return { + "region": cfg.get("region", None), + "mag_min": cfg.get("mag_min", self.region_config.get("mag_min")), + "mag_max": cfg.get("mag_max", self.region_config.get("mag_max")), + "depth_min": cfg.get("depth_min", None), + "depth_max": cfg.get("depth_max", None), + "start_date": cfg.get("start_date", None), + } + @property def catalog(self) -> CSEPCatalog: """ @@ -126,22 +152,7 @@ def catalog(self, cat: Union[Callable, CSEPCatalog, str]) -> None: self.name = cat else: query_function = parse_csep_func(cat) - bounds = { - "start_time": min([item for sublist in self.time_windows for item in sublist]), - "end_time": max([item for sublist in self.time_windows for item in sublist]), - "min_magnitude": self.magnitudes.min(), - "max_depth": self.depths.max(), - } - if self.region: - bounds.update( - { - i: j - for i, j in zip( - ["min_longitude", "max_longitude", "min_latitude", "max_latitude"], - self.region.get_bbox(), - ) - } - ) + bounds = self.query_bounds() self._catalog = query_function(catalog_id="catalog", **bounds) self.cat_path = self.registry.rel("catalog.json") @@ -154,6 +165,44 @@ def catalog(self, cat: Union[Callable, CSEPCatalog, str]) -> None: else: log.info(f"\tCatalog: '{cat}'") + def query_bounds(self) -> dict: + """ + Bounds to query the main catalog from a network API, wide enough to cover both the + testing settings (``region_config``, ``time_config``) and the input catalog settings + (``input_config``). + + Returns: + Keyword arguments for the catalog query function. + """ + inp = self.input_cat_config + start = min([item for sublist in self.time_windows for item in sublist]) + if inp["start_date"]: + start = min(start, inp["start_date"]) + mag_min = self.magnitudes.min() + if inp["mag_min"] is not None: + mag_min = min(mag_min, inp["mag_min"]) + depth_max = self.depths.max() + if inp["depth_max"] is not None: + depth_max = max(depth_max, inp["depth_max"]) + bounds = { + "start_time": start, + "end_time": max([item for sublist in self.time_windows for item in sublist]), + "min_magnitude": mag_min, + "max_depth": depth_max, + } + regions = [r for r in (self.region, inp["region"]) if r is not None] + if regions: + bboxes = numpy.array([r.get_bbox() for r in regions]) + bounds.update( + { + "min_longitude": bboxes[:, 0].min(), + "max_longitude": bboxes[:, 1].max(), + "min_latitude": bboxes[:, 2].min(), + "max_latitude": bboxes[:, 3].max(), + } + ) + return bounds + def get_test_cat(self, tstring: str = None, fmt: str = "json") -> CSEPCatalog: """ Filters the complete experiment catalog to a test sub-catalog bounded by the test @@ -170,31 +219,51 @@ def get_test_cat(self, tstring: str = None, fmt: str = "json") -> CSEPCatalog: return test_catalog + def get_input_cat(self, tstring: str) -> CSEPCatalog: + """ + Filters the main catalog to the input catalog of a time window: all the events before + the window start, within the settings of ``input_config.catalog`` (region, magnitude, + depth and lower time bound). See :func:`~floatcsep.utils.helpers.read_input_cfg`. + + Args: + tstring (str): Time window string + + Returns: + The input catalog as a :class:`csep.core.catalogs.CSEPCatalog` + """ + start, end = str2timewindow(tstring) + cfg = self.input_cat_config + filters = [f"origin_time < {start.timestamp() * 1000}"] + if cfg["start_date"]: + filters.append(f"origin_time >= {cfg['start_date'].timestamp() * 1000}") + if cfg["mag_min"] is not None: + filters.append(f"magnitude >= {cfg['mag_min']}") + if cfg["mag_max"] is not None: + filters.append(f"magnitude < {cfg['mag_max']}") + if cfg["depth_min"] is not None: + filters.append(f"depth >= {cfg['depth_min']}") + if cfg["depth_max"] is not None: + filters.append(f"depth < {cfg['depth_max']}") + sub_cat = self.catalog.filter(filters, in_place=False) + if cfg["region"] is not None and sub_cat.get_number_of_events() > 0: + sub_cat.filter_spatial(region=cfg["region"], in_place=True) + return sub_cat + def set_input_cats(self, tstring: str, models: List["Model"], fmt: str = "ascii") -> None: """ - Filters the complete experiment catalog to input sub-catalog filtered to the beginning - of the test time-window. + Filters the complete experiment catalog to the input sub-catalog of a time window (see + :meth:`get_input_cat`) and writes it to the input directory of each model. Args: tstring (str): Time window string - model (:class:`~floatcsep.model.Model`): Model to give the input - catalog + models (list of :class:`~floatcsep.model.Model`): Models to give the input catalog fmt (str): Output catalog format """ - start, end = str2timewindow(tstring) - log.debug(f"[Catalogs] Filtering input catalog and saving to models' input directory") + log.debug("[Catalogs] Filtering input catalog and saving to models' input directory") + sub_cat = self.get_input_cat(tstring) + writer = getattr(CatalogSerializer, fmt) for model in models: - input_cat_name = model.registry.get_input_catalog_key(tstring) - sub_cat = self.catalog.filter( - [ - f"origin_time < {start.timestamp() * 1000}", - f"magnitude >= {self.mag_min}", - f"magnitude < {self.mag_max}", - ], - in_place=False, - ) - writer = getattr(CatalogSerializer, fmt) - writer(catalog=sub_cat, filename=input_cat_name) + writer(catalog=sub_cat, filename=model.registry.get_input_catalog_key(tstring)) def set_test_cats(self, tstring: str, fmt: str = "json") -> None: """ @@ -213,16 +282,20 @@ def set_test_cats(self, tstring: str, fmt: str = "json") -> None: f"{self.registry.rel(test_cat_name)}" ) start, end = str2timewindow(tstring) - sub_cat = self.catalog.filter( - [ - f"origin_time < {end.timestamp() * 1000}", - f"origin_time >= {start.timestamp() * 1000}", - f"magnitude >= {self.mag_min}", - f"magnitude < {self.mag_max}", - ], - in_place=False, - ) - if self.region: + filters = [ + f"origin_time < {end.timestamp() * 1000}", + f"origin_time >= {start.timestamp() * 1000}", + f"magnitude >= {self.mag_min}", + f"magnitude < {self.mag_max}", + ] + depth_min = self.region_config.get("depth_min", None) + depth_max = self.region_config.get("depth_max", None) + if depth_min is not None: + filters.append(f"depth >= {depth_min}") + if depth_max is not None: + filters.append(f"depth < {depth_max}") + sub_cat = self.catalog.filter(filters, in_place=False) + if self.region and sub_cat.get_number_of_events() > 0: sub_cat.filter_spatial(region=self.region, in_place=True) writer = getattr(CatalogSerializer, fmt) diff --git a/floatcsep/utils/helpers.py b/floatcsep/utils/helpers.py index 5d446a7..3093984 100644 --- a/floatcsep/utils/helpers.py +++ b/floatcsep/utils/helpers.py @@ -172,6 +172,55 @@ def read_time_cfg(time_config, **kwargs): return time_config +def parse_region(region_data, magnitudes=None, path=""): + """ + Resolves a region specification into a :class:`csep.core.regions.CartesianGrid2D`. + + Args: + region_data (str, dict, CartesianGrid2D): Name of a pyCSEP/floatCSEP region function, + a path to a text file with the cells' origins, or a region dictionary. + magnitudes (array-like): Magnitude bins bound to the region. Can be None if the region + is only used for spatial filtering. + path (str): Directory from which to resolve a region file. + + Returns: + A tuple ``(region, region_file)``. ``region_file`` is the path of the region file as + given, or None if the region was not read from a file. + """ + if region_data is None: + return None, None + if isinstance(region_data, CartesianGrid2D): + return region_data, None + try: + return parse_csep_func(region_data)(name=region_data, magnitudes=magnitudes), None + except AttributeError: + pass + + if isinstance(region_data, str): + filename = os.path.join(path, region_data) + with open(filename, "r") as file_: + parsed_region = file_.readlines() + try: + data = numpy.array( + [re.split(r"\s+|,", i.strip()) for i in parsed_region], dtype=float + ) + except ValueError: + data = numpy.array( + [re.split(r"\s+|,", i.strip()) for i in parsed_region[1:]], dtype=float + ) + dh1 = scipy.stats.mode(numpy.diff(numpy.unique(data[:, 0]))).mode + dh2 = scipy.stats.mode(numpy.diff(numpy.unique(data[:, 1]))).mode + dh = numpy.nanmin([dh1, dh2]) + region = CartesianGrid2D.from_origins( + data, name=region_data, magnitudes=magnitudes, dh=dh + ) + return region, region_data + + region_data = dict(region_data) + region_data["magnitudes"] = magnitudes + return CartesianGrid2D.from_dict(region_data), None + + def read_region_cfg(region_config, **kwargs): """ Builds the region configuration of an experiment. @@ -202,42 +251,74 @@ def read_region_cfg(region_config, **kwargs): magbin = region_config["mag_bin"] magnitudes = cleaner_range(magmin, magmax, magbin) - region_data = region_config.get("region", None) - try: - region = ( - parse_csep_func(region_data)(name=region_data, magnitudes=magnitudes) - if region_data - else None - ) - except AttributeError: - if isinstance(region_data, str): - filename = os.path.join(kwargs.get("path", ""), region_data) - with open(filename, "r") as file_: - parsed_region = file_.readlines() - try: - data = numpy.array( - [re.split(r"\s+|,", i.strip()) for i in parsed_region], dtype=float - ) - except ValueError: - data = numpy.array( - [re.split(r"\s+|,", i.strip()) for i in parsed_region[1:]], dtype=float - ) - dh1 = scipy.stats.mode(numpy.diff(numpy.unique(data[:, 0]))).mode - dh2 = scipy.stats.mode(numpy.diff(numpy.unique(data[:, 1]))).mode - dh = numpy.nanmin([dh1, dh2]) - region = CartesianGrid2D.from_origins( - data, name=region_data, magnitudes=magnitudes, dh=dh - ) - region_config.update({"path": region_data}) - else: - region_data["magnitudes"] = magnitudes - region = CartesianGrid2D.from_dict(region_data) + region, region_file = parse_region( + region_config.get("region", None), magnitudes, kwargs.get("path", "") + ) + if region_file: + region_config["path"] = region_file region_config.update({"depths": depths, "magnitudes": magnitudes, "region": region}) return region_config +def read_input_cfg(input_config, region_config=None, **kwargs): + """ + Builds the configuration of the input data handed to the models (e.g., the input catalog of + a time-dependent model), as opposed to the testing data defined by ``region_config`` and + ``time_config``. + + The ``catalog`` block defines how the main catalog is filtered before being written as a + model's input catalog, for each time window. Only the window end (the forecast start date) + is set by the experiment. Admissible keys are ``region``, ``mag_min``, ``mag_max``, + ``depth_min``, ``depth_max`` and ``start_date``. Missing magnitudes default to those of + ``region_config``. A missing region, depths or start date means no filter, which is also + what an absent ``input_config`` does. + + Args: + input_config (dict): Dictionary with an (optional) ``catalog`` block. + region_config (dict): Parsed experiment region configuration, used for defaults. + **kwargs: ``path`` is used to resolve a region file relative to the config file. + + Returns: + A dictionary with a ``catalog`` block holding the resolved input catalog settings, or + an empty dictionary if no input configuration was given. + """ + if not input_config: + return {} + input_config = copy.deepcopy(input_config) + region_config = region_config or {} + cat_cfg = input_config.get("catalog") or {} + _attrs = ["region", "mag_min", "mag_max", "depth_min", "depth_max", "start_date"] + unknown = set(cat_cfg) - set(_attrs) - {"path"} + if unknown: + raise ValueError( + f"Unknown keys in input_config.catalog: {sorted(unknown)}. " + f"Admissible keys are {_attrs}" + ) + + cat_cfg.setdefault("mag_min", region_config.get("mag_min")) + cat_cfg.setdefault("mag_max", region_config.get("mag_max")) + + start = cat_cfg.get("start_date") + if isinstance(start, str): + cat_cfg["start_date"] = datetime.fromisoformat(start) + elif isinstance(start, date) and not isinstance(start, datetime): + cat_cfg["start_date"] = datetime(start.year, start.month, start.day) + + region, region_file = parse_region( + cat_cfg.get("region", None), magnitudes=None, path=kwargs.get("path", "") + ) + cat_cfg["region"] = region + if region_file: + cat_cfg["path"] = region_file + else: + cat_cfg.pop("path", None) + + input_config["catalog"] = cat_cfg + return input_config + + def timewindow2str(datetimes: Sequence) -> Union[str, list[str]]: """ Converts a time window (list/tuple of datetimes) to a string that represents it. Can be a diff --git a/tests/unit/test_experiment.py b/tests/unit/test_experiment.py index 237fe8a..c8eee96 100644 --- a/tests/unit/test_experiment.py +++ b/tests/unit/test_experiment.py @@ -98,6 +98,55 @@ def test_to_dict(self): } self.assertEqual(dict_, exp_a.as_dict()) + def test_input_config(self): + input_config = { + "catalog": { + "region": "italy_csep_collection_region", + "mag_min": 2.5, + "depth_max": 30, + "start_date": datetime(2005, 1, 1), + } + } + exp = Experiment( + name="test", + **_time_config, + **_region_config, + input_config=input_config, + catalog=_cat, + ) + self.assertEqual(2.5, exp.catalog_repo.input_cat_config["mag_min"]) + self.assertEqual(10.0, exp.catalog_repo.input_cat_config["mag_max"]) + self.assertEqual( + "italy_csep_collection_region", exp.input_config["catalog"]["region"].name + ) + # the testing settings are untouched + self.assertEqual(1.0, exp.mag_min) + + dict_ = exp.as_dict() + self.assertEqual( + { + "catalog": { + "region": "italy_csep_collection_region", + "mag_min": 2.5, + "mag_max": 10.0, + "depth_max": 30, + "start_date": datetime(2005, 1, 1), + } + }, + dict_["input_config"], + ) + keys = list(dict_) + self.assertLess(keys.index("region_config"), keys.index("input_config")) + exp_b = Experiment(**{**dict_, "path": os.getcwd()}) + self.assertEqual( + exp.catalog_repo.input_cat_config["mag_min"], + exp_b.catalog_repo.input_cat_config["mag_min"], + ) + + exp_c = Experiment(name="test", **_time_config, **_region_config, catalog=_cat) + self.assertEqual({}, exp_c.input_config) + self.assertNotIn("input_config", exp_c.as_dict()) + def test_to_yml(self): time_config = { "start_date": datetime(2021, 1, 1), diff --git a/tests/unit/test_repositories.py b/tests/unit/test_repositories.py index 5b5ce58..38fbf69 100644 --- a/tests/unit/test_repositories.py +++ b/tests/unit/test_repositories.py @@ -1,8 +1,16 @@ import datetime +import json +import os +import tempfile import unittest from unittest.mock import MagicMock, patch, PropertyMock, mock_open +import numpy +from csep.core.catalogs import CSEPCatalog from csep.core.forecasts import GriddedForecast +from csep.utils.time_utils import datetime_to_utc_epoch + +from floatcsep.utils.helpers import read_time_cfg, read_region_cfg, read_input_cfg from floatcsep.utils.file_io import GriddedForecastParsers from floatcsep.infrastructure.registries import ModelFileRegistry @@ -214,5 +222,166 @@ def test_set_catalog(self, mock_reader, mock_isfile): self.assertEqual(self.catalog_repo._catalog, "csep catalog") +class TestCatalogFiltering(unittest.TestCase): + """Input and test catalog filtering, with a synthetic catalog around the Italy region.""" + + window = "2016-08-24_2016-08-25" + + @staticmethod + def _epoch(*args): + return datetime_to_utc_epoch(datetime.datetime(*args, tzinfo=datetime.timezone.utc)) + + @classmethod + def setUpClass(cls): + # id, time, lat, lon, depth, mag + rows = [ + (1, cls._epoch(2010, 1, 1), 42.5, 13.0, 10, 2.6), # small, before window + (2, cls._epoch(2012, 1, 1), 42.5, 13.0, 10, 4.5), # before window + (3, cls._epoch(2014, 1, 1), 42.5, 13.0, 40, 4.5), # deep, before window + (4, cls._epoch(2015, 1, 1), 48.0, 13.0, 10, 4.5), # outside region, before window + (5, cls._epoch(2016, 8, 24, 5), 42.5, 13.0, 10, 4.5), # in window + (6, cls._epoch(2016, 8, 24, 6), 42.5, 13.0, 50, 4.5), # in window, deep + (7, cls._epoch(2008, 1, 1), 42.5, 13.0, 10, 4.5), # before input start_date + (8, cls._epoch(2016, 8, 24, 7), 48.0, 13.0, 10, 4.5), # in window, outside region + ] + dtype = numpy.dtype( + [ + ("id", "S256"), + ("origin_time", "