From 0dc7c583f7d49957f3a614e1ee75343be8efd9f3 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 18 Dec 2023 13:50:12 -0500 Subject: [PATCH 1/8] Add probability field calc and plots --- ensembleperturbation/plotting/surrogate.py | 73 +++++++++ .../uncertainty_quantification/surrogate.py | 151 ++++++++++++++++++ 2 files changed, 224 insertions(+) diff --git a/ensembleperturbation/plotting/surrogate.py b/ensembleperturbation/plotting/surrogate.py index 861b2f3e..a41ee66f 100644 --- a/ensembleperturbation/plotting/surrogate.py +++ b/ensembleperturbation/plotting/surrogate.py @@ -564,3 +564,76 @@ def plot_kl_surrogate_fit( if output_filename is not None: figure.savefig(output_filename, dpi=200, bbox_inches='tight') + + +def plot_selected_probability_fields( + node_prob_field: xarray.Dataset, level_list: list, output_directory: PathLike +): + levels = node_prob_field.levels + + sources = node_prob_field['source'].values + if output_directory is not None: + if not isinstance(output_directory, Path): + output_directory = Path(output_directory) + + bounds = numpy.array( + [ + node_prob_field['x'].min(), + node_prob_field['y'].min(), + node_prob_field['x'].max(), + node_prob_field['y'].max(), + ] + ) + vmax = numpy.round_(levels.sel(source='model').quantile(0.98), decimals=1) + vmin = 0.0 + for lvl in level_list: + figure = pyplot.figure() + figure.set_size_inches(10, 10 / 1.61803398875) + figure.suptitle(f'comparison of probability fields at level: {lvl}') + for index, source in enumerate(sources): + map_axis = figure.add_subplot(2, len(sources), index + 1) + map_axis.title.set_text(f'{source}') + countries = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres')) + + map_axis.set_xlim((bounds[0], bounds[2])) + map_axis.set_ylim((bounds[1], bounds[3])) + + xlim = map_axis.get_xlim() + ylim = map_axis.get_ylim() + + countries.plot(color='lightgrey', ax=map_axis) + + points = numpy.vstack( + ( + node_prob_field['x'], + node_prob_field['y'], + levels.sel(level=lvl, source=source), + ) + ).T + if 'element' not in node_prob_field: + im = plot_points( + points=points, axis=map_axis, add_colorbar=False, vmax=vmax, vmin=vmin, + ) + else: + im = plot_surface( + points=points, + element_table=node_prob_field['element'].values, + axis=map_axis, + add_colorbar=False, + levels=numpy.linspace(vmin, vmax, 25 + 1), + extend='both', + ) + + map_axis.set_xlim(xlim) + map_axis.set_ylim(ylim) + + pyplot.subplots_adjust(wspace=0.02, right=0.96) + cax = pyplot.axes([0.95, 0.55, 0.015, 0.3]) + cbar = figure.colorbar(im, extend='both', cax=cax) + + if output_directory is not None: + figure.savefig( + output_directory / f'probability_exceeding_{lvl}m.png', + dpi=200, + bbox_inches='tight', + ) diff --git a/ensembleperturbation/uncertainty_quantification/surrogate.py b/ensembleperturbation/uncertainty_quantification/surrogate.py index db2c6163..49012c93 100644 --- a/ensembleperturbation/uncertainty_quantification/surrogate.py +++ b/ensembleperturbation/uncertainty_quantification/surrogate.py @@ -606,3 +606,154 @@ def compute_surrogate_percentiles( out = out.reshape(q.shape + shape) return out + + +def compute_surrogate_probability_field( + poly: numpoly.ndpoly, + levels: List[float], + dist: chaospy.Distribution, + sample: int = 10000, + minimum_allowable_value: float = None, + convert_from_log_scale: Union[bool, float] = False, + convert_from_depths: Union[bool, float] = False, + **kws, +): + + poly = chaospy.aspolynomial(poly) + shape = poly.shape + poly = poly.ravel() + + levels = numpy.asarray(levels).ravel() + dim = len(dist) + + # Interior + Z = dist.sample(sample, **kws).reshape(len(dist), sample) + poly1 = poly(*Z) + + # Min/max + ext = numpy.mgrid[(slice(0, 2, 1),) * dim].reshape(dim, 2 ** dim).T + ext = numpy.where(ext, dist.lower, dist.upper).T + poly2 = poly(*ext) + poly2 = numpy.array([_ for _ in poly2.T if not numpy.any(numpy.isnan(_))]).T + + # Finish + if poly2.shape: + poly1 = numpy.concatenate([poly1, poly2], -1) + if isinstance(convert_from_log_scale, float): + poly1 = convert_from_log_scale ** poly1 + elif convert_from_log_scale: + poly1 = numpy.exp(poly1) + samples = poly1.shape[1] + + # Adjustments and elev corrections + if isinstance(convert_from_depths, (float, numpy.ndarray)): + poly1 -= convert_from_depths + if minimum_allowable_value is not None: + too_small = poly1 < minimum_allowable_value + poly1[too_small] = numpy.nan + if isinstance(convert_from_depths, (float, numpy.ndarray)) or convert_from_depths: + # poly1 -= training_set['depth'] + poly1 -= training_set['depth'].values[:, None] + + out = (poly1[:, :, None] > (levels[None, None, :])).sum(axis=1).T / samples + + out = out.reshape(levels.shape + shape) + + return out + + +def probability_field_from_samples( + samples: xarray.DataArray, + levels: List[float], + surrogate_model: numpoly.ndpoly, + distribution: chaospy.Distribution, + minimum_allowable_value: float = None, + convert_from_log_scale: Union[bool, float] = False, + convert_from_depths: Union[bool, float] = False, +) -> xarray.DataArray: + + LOGGER.info(f'calculating {len(levels)} probability field(s): {levels}') + + surrogate_prob_field = compute_surrogate_probability_field( + poly=surrogate_model, + levels=levels, + dist=distribution, + minimum_allowable_value=minimum_allowable_value, + convert_from_log_scale=convert_from_log_scale, + convert_from_depths=convert_from_depths, + ) + + surrogate_prob_field = xarray.DataArray( + surrogate_prob_field, + coords={ + 'level': levels, + **{ + coord: values + for coord, values in samples.coords.items() + if coord not in ['run', 'type'] + }, + }, + dims=('level', *(dim for dim in samples.dims if dim not in ['run', 'type'])), + ) + + return surrogate_prob_field + + +def probability_field_from_surrogate( + levels: List[float], + surrogate_model: numpoly.ndpoly, + distribution: chaospy.Distribution, + training_set: xarray.Dataset, + minimum_allowable_value: float = None, + convert_from_log_scale: Union[bool, float] = False, + convert_from_depths: Union[bool, float] = False, + element_table: xarray.DataArray = None, + filename: PathLike = None, +) -> xarray.Dataset: + + if filename is not None and not isinstance(filename, Path): + filename = Path(filename) + + if filename is None or not filename.exists(): + surrogate_prob_field = probability_field_from_samples( + samples=training_set, + levels=levels, + surrogate_model=surrogate_model, + distribution=distribution, + minimum_allowable_value=minimum_allowable_value, + convert_from_log_scale=convert_from_log_scale, + convert_from_depths=convert_from_depths, + ) + + # before evaluating prob. field for model set null water elevation to the ground elevation + training_set = numpy.fmax(training_set, -training_set['depth']) + if minimum_allowable_value is not None: + too_small = (training_set + training_set['depth']).values < minimum_allowable_value + training_set.values[too_small] = numpy.nan + + ds1, ds2 = xarray.broadcast(training_set, surrogate_prob_field['level']) + modeled_prob_field = (ds1 >= ds2).sum(dim='run') / len(training_set.run) + + # modeled_prob_field.coords['level'] = surrogate_prob_field['level'] + + node_prob_field = xarray.combine_nested( + [surrogate_prob_field, modeled_prob_field], concat_dim='source' + ).assign_coords(source=['surrogate', 'model']) + + node_prob_field = node_prob_field.to_dataset(name='levels') + + node_prob_field = node_prob_field.assign( + differences=numpy.fabs(surrogate_prob_field - modeled_prob_field) + ) + + if element_table is not None: + node_prob_field = node_prob_field.assign_coords({'element': element_table}) + + if filename is not None: + LOGGER.info(f'saving prob_field to "{filename}"') + node_prob_field.to_netcdf(filename) + else: + LOGGER.info(f'loading prob_field from "{filename}"') + node_prob_field = xarray.open_dataset(filename) + + return node_prob_field From fad0e19877def5d7aba0bc81a1666a6da4ef309c Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Mon, 18 Dec 2023 14:12:15 -0500 Subject: [PATCH 2/8] Fix bugs --- .../uncertainty_quantification/surrogate.py | 116 +++++++++--------- 1 file changed, 58 insertions(+), 58 deletions(-) diff --git a/ensembleperturbation/uncertainty_quantification/surrogate.py b/ensembleperturbation/uncertainty_quantification/surrogate.py index 49012c93..ff58c055 100644 --- a/ensembleperturbation/uncertainty_quantification/surrogate.py +++ b/ensembleperturbation/uncertainty_quantification/surrogate.py @@ -608,62 +608,8 @@ def compute_surrogate_percentiles( return out -def compute_surrogate_probability_field( - poly: numpoly.ndpoly, - levels: List[float], - dist: chaospy.Distribution, - sample: int = 10000, - minimum_allowable_value: float = None, - convert_from_log_scale: Union[bool, float] = False, - convert_from_depths: Union[bool, float] = False, - **kws, -): - - poly = chaospy.aspolynomial(poly) - shape = poly.shape - poly = poly.ravel() - - levels = numpy.asarray(levels).ravel() - dim = len(dist) - - # Interior - Z = dist.sample(sample, **kws).reshape(len(dist), sample) - poly1 = poly(*Z) - - # Min/max - ext = numpy.mgrid[(slice(0, 2, 1),) * dim].reshape(dim, 2 ** dim).T - ext = numpy.where(ext, dist.lower, dist.upper).T - poly2 = poly(*ext) - poly2 = numpy.array([_ for _ in poly2.T if not numpy.any(numpy.isnan(_))]).T - - # Finish - if poly2.shape: - poly1 = numpy.concatenate([poly1, poly2], -1) - if isinstance(convert_from_log_scale, float): - poly1 = convert_from_log_scale ** poly1 - elif convert_from_log_scale: - poly1 = numpy.exp(poly1) - samples = poly1.shape[1] - - # Adjustments and elev corrections - if isinstance(convert_from_depths, (float, numpy.ndarray)): - poly1 -= convert_from_depths - if minimum_allowable_value is not None: - too_small = poly1 < minimum_allowable_value - poly1[too_small] = numpy.nan - if isinstance(convert_from_depths, (float, numpy.ndarray)) or convert_from_depths: - # poly1 -= training_set['depth'] - poly1 -= training_set['depth'].values[:, None] - - out = (poly1[:, :, None] > (levels[None, None, :])).sum(axis=1).T / samples - - out = out.reshape(levels.shape + shape) - - return out - - def probability_field_from_samples( - samples: xarray.DataArray, + samples: xarray.Dataset, levels: List[float], surrogate_model: numpoly.ndpoly, distribution: chaospy.Distribution, @@ -681,6 +627,7 @@ def probability_field_from_samples( minimum_allowable_value=minimum_allowable_value, convert_from_log_scale=convert_from_log_scale, convert_from_depths=convert_from_depths, + depths=samples['depth'], ) surrogate_prob_field = xarray.DataArray( @@ -726,16 +673,14 @@ def probability_field_from_surrogate( ) # before evaluating prob. field for model set null water elevation to the ground elevation - training_set = numpy.fmax(training_set, -training_set['depth']) if minimum_allowable_value is not None: too_small = (training_set + training_set['depth']).values < minimum_allowable_value training_set.values[too_small] = numpy.nan + training_set = numpy.fmax(training_set, -training_set['depth']) ds1, ds2 = xarray.broadcast(training_set, surrogate_prob_field['level']) modeled_prob_field = (ds1 >= ds2).sum(dim='run') / len(training_set.run) - # modeled_prob_field.coords['level'] = surrogate_prob_field['level'] - node_prob_field = xarray.combine_nested( [surrogate_prob_field, modeled_prob_field], concat_dim='source' ).assign_coords(source=['surrogate', 'model']) @@ -757,3 +702,58 @@ def probability_field_from_surrogate( node_prob_field = xarray.open_dataset(filename) return node_prob_field + + +def compute_surrogate_probability_field( + poly: numpoly.ndpoly, + levels: List[float], + dist: chaospy.Distribution, + sample: int = 10000, + minimum_allowable_value: float = None, + convert_from_log_scale: Union[bool, float] = False, + convert_from_depths: Union[bool, float] = False, + depths: xarray.DataArray = None, + **kws, +): + + poly = chaospy.aspolynomial(poly) + shape = poly.shape + poly = poly.ravel() + + levels = numpy.asarray(levels).ravel() + dim = len(dist) + + # Interior + Z = dist.sample(sample, **kws).reshape(len(dist), sample) + poly1 = poly(*Z) + + # Min/max + ext = numpy.mgrid[(slice(0, 2, 1),) * dim].reshape(dim, 2 ** dim).T + ext = numpy.where(ext, dist.lower, dist.upper).T + poly2 = poly(*ext) + poly2 = numpy.array([_ for _ in poly2.T if not numpy.any(numpy.isnan(_))]).T + + # Finish + if poly2.shape: + poly1 = numpy.concatenate([poly1, poly2], -1) + if isinstance(convert_from_log_scale, float): + poly1 = convert_from_log_scale ** poly1 + elif convert_from_log_scale: + poly1 = numpy.exp(poly1) + samples = poly1.shape[1] + + # adjustments and elev corrections + if isinstance(convert_from_depths, (float, numpy.ndarray)): + poly1 -= convert_from_depths + if minimum_allowable_value is not None: + too_small = poly1 < minimum_allowable_value + poly1[too_small] = numpy.nan + if isinstance(convert_from_depths, (float, numpy.ndarray)) or convert_from_depths: + # TODO: Sanity check for depth vs poly shapes + poly1 -= depths.values[:, None] + + out = (poly1[:, :, None] > (levels[None, None, :])).sum(axis=1).T / samples + + out = out.reshape(levels.shape + shape) + + return out From 4d82917f3e1fcce456fd7b5a587da07b7b8658e7 Mon Sep 17 00:00:00 2001 From: "Soroosh.Mani" Date: Tue, 19 Dec 2023 08:30:58 -0500 Subject: [PATCH 3/8] Update prob plot title --- ensembleperturbation/plotting/surrogate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ensembleperturbation/plotting/surrogate.py b/ensembleperturbation/plotting/surrogate.py index a41ee66f..7a337953 100644 --- a/ensembleperturbation/plotting/surrogate.py +++ b/ensembleperturbation/plotting/surrogate.py @@ -589,7 +589,7 @@ def plot_selected_probability_fields( for lvl in level_list: figure = pyplot.figure() figure.set_size_inches(10, 10 / 1.61803398875) - figure.suptitle(f'comparison of probability fields at level: {lvl}') + figure.suptitle(f'Probability of water level exceeding {lvl}-m') for index, source in enumerate(sources): map_axis = figure.add_subplot(2, len(sources), index + 1) map_axis.title.set_text(f'{source}') From 5364e2c878ae3144d279fd0ec1183e024864dd0c Mon Sep 17 00:00:00 2001 From: SorooshMani-NOAA Date: Tue, 19 Dec 2023 15:23:43 +0000 Subject: [PATCH 4/8] Swap order of elev model cleanup --- ensembleperturbation/uncertainty_quantification/surrogate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ensembleperturbation/uncertainty_quantification/surrogate.py b/ensembleperturbation/uncertainty_quantification/surrogate.py index ff58c055..5b04d7d0 100644 --- a/ensembleperturbation/uncertainty_quantification/surrogate.py +++ b/ensembleperturbation/uncertainty_quantification/surrogate.py @@ -673,10 +673,10 @@ def probability_field_from_surrogate( ) # before evaluating prob. field for model set null water elevation to the ground elevation + training_set = numpy.fmax(training_set, -training_set['depth']) if minimum_allowable_value is not None: too_small = (training_set + training_set['depth']).values < minimum_allowable_value training_set.values[too_small] = numpy.nan - training_set = numpy.fmax(training_set, -training_set['depth']) ds1, ds2 = xarray.broadcast(training_set, surrogate_prob_field['level']) modeled_prob_field = (ds1 >= ds2).sum(dim='run') / len(training_set.run) From e8006542dea76952f3b9452e8d94854d23141450 Mon Sep 17 00:00:00 2001 From: SorooshMani-NOAA Date: Tue, 19 Dec 2023 20:48:40 +0000 Subject: [PATCH 5/8] Update based on PR comments --- ensembleperturbation/plotting/surrogate.py | 6 +++--- .../uncertainty_quantification/surrogate.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ensembleperturbation/plotting/surrogate.py b/ensembleperturbation/plotting/surrogate.py index 7a337953..5d05904d 100644 --- a/ensembleperturbation/plotting/surrogate.py +++ b/ensembleperturbation/plotting/surrogate.py @@ -584,7 +584,7 @@ def plot_selected_probability_fields( node_prob_field['y'].max(), ] ) - vmax = numpy.round_(levels.sel(source='model').quantile(0.98), decimals=1) + vmax = 1.0 vmin = 0.0 for lvl in level_list: figure = pyplot.figure() @@ -621,7 +621,7 @@ def plot_selected_probability_fields( axis=map_axis, add_colorbar=False, levels=numpy.linspace(vmin, vmax, 25 + 1), - extend='both', + extend='neither', ) map_axis.set_xlim(xlim) @@ -629,7 +629,7 @@ def plot_selected_probability_fields( pyplot.subplots_adjust(wspace=0.02, right=0.96) cax = pyplot.axes([0.95, 0.55, 0.015, 0.3]) - cbar = figure.colorbar(im, extend='both', cax=cax) + cbar = figure.colorbar(im, extend='neither', cax=cax) if output_directory is not None: figure.savefig( diff --git a/ensembleperturbation/uncertainty_quantification/surrogate.py b/ensembleperturbation/uncertainty_quantification/surrogate.py index 5b04d7d0..91fd1844 100644 --- a/ensembleperturbation/uncertainty_quantification/surrogate.py +++ b/ensembleperturbation/uncertainty_quantification/surrogate.py @@ -673,7 +673,7 @@ def probability_field_from_surrogate( ) # before evaluating prob. field for model set null water elevation to the ground elevation - training_set = numpy.fmax(training_set, -training_set['depth']) + #training_set = numpy.fmax(training_set, -training_set['depth']) if minimum_allowable_value is not None: too_small = (training_set + training_set['depth']).values < minimum_allowable_value training_set.values[too_small] = numpy.nan From ef85570f89eb3056bfce510e021097fc2dfad840 Mon Sep 17 00:00:00 2001 From: SorooshMani-NOAA Date: Tue, 19 Dec 2023 20:50:47 +0000 Subject: [PATCH 6/8] lint fix --- ensembleperturbation/uncertainty_quantification/surrogate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ensembleperturbation/uncertainty_quantification/surrogate.py b/ensembleperturbation/uncertainty_quantification/surrogate.py index 91fd1844..45ce3145 100644 --- a/ensembleperturbation/uncertainty_quantification/surrogate.py +++ b/ensembleperturbation/uncertainty_quantification/surrogate.py @@ -673,7 +673,7 @@ def probability_field_from_surrogate( ) # before evaluating prob. field for model set null water elevation to the ground elevation - #training_set = numpy.fmax(training_set, -training_set['depth']) + # training_set = numpy.fmax(training_set, -training_set['depth']) if minimum_allowable_value is not None: too_small = (training_set + training_set['depth']).values < minimum_allowable_value training_set.values[too_small] = numpy.nan From 92b8e6d6f1c544ca30c9f0e1fa1afa9f0c50f180 Mon Sep 17 00:00:00 2001 From: SorooshMani-NOAA Date: Wed, 20 Dec 2023 02:33:00 +0000 Subject: [PATCH 7/8] Update dataarray name & vmax --- ensembleperturbation/plotting/surrogate.py | 10 +++++----- .../uncertainty_quantification/surrogate.py | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/ensembleperturbation/plotting/surrogate.py b/ensembleperturbation/plotting/surrogate.py index 5d05904d..446e2b0b 100644 --- a/ensembleperturbation/plotting/surrogate.py +++ b/ensembleperturbation/plotting/surrogate.py @@ -569,7 +569,7 @@ def plot_kl_surrogate_fit( def plot_selected_probability_fields( node_prob_field: xarray.Dataset, level_list: list, output_directory: PathLike ): - levels = node_prob_field.levels + probabilities = node_prob_field.probabilities sources = node_prob_field['source'].values if output_directory is not None: @@ -584,8 +584,8 @@ def plot_selected_probability_fields( node_prob_field['y'].max(), ] ) - vmax = 1.0 - vmin = 0.0 + vmax = 1 + vmin = 0 for lvl in level_list: figure = pyplot.figure() figure.set_size_inches(10, 10 / 1.61803398875) @@ -607,7 +607,7 @@ def plot_selected_probability_fields( ( node_prob_field['x'], node_prob_field['y'], - levels.sel(level=lvl, source=source), + probabilities.sel(level=lvl, source=source), ) ).T if 'element' not in node_prob_field: @@ -620,7 +620,7 @@ def plot_selected_probability_fields( element_table=node_prob_field['element'].values, axis=map_axis, add_colorbar=False, - levels=numpy.linspace(vmin, vmax, 25 + 1), + levels=numpy.linspace(vmin, vmax, 25), extend='neither', ) diff --git a/ensembleperturbation/uncertainty_quantification/surrogate.py b/ensembleperturbation/uncertainty_quantification/surrogate.py index 45ce3145..dc158094 100644 --- a/ensembleperturbation/uncertainty_quantification/surrogate.py +++ b/ensembleperturbation/uncertainty_quantification/surrogate.py @@ -685,7 +685,7 @@ def probability_field_from_surrogate( [surrogate_prob_field, modeled_prob_field], concat_dim='source' ).assign_coords(source=['surrogate', 'model']) - node_prob_field = node_prob_field.to_dataset(name='levels') + node_prob_field = node_prob_field.to_dataset(name='probabilities') node_prob_field = node_prob_field.assign( differences=numpy.fabs(surrogate_prob_field - modeled_prob_field) From 830ede3977bcce396aa49fa4ed9c23ecc7356cd2 Mon Sep 17 00:00:00 2001 From: SorooshMani-NOAA Date: Wed, 20 Dec 2023 02:41:41 +0000 Subject: [PATCH 8/8] Fix vmax cutoff issue --- ensembleperturbation/plotting/surrogate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ensembleperturbation/plotting/surrogate.py b/ensembleperturbation/plotting/surrogate.py index 446e2b0b..607ca144 100644 --- a/ensembleperturbation/plotting/surrogate.py +++ b/ensembleperturbation/plotting/surrogate.py @@ -584,7 +584,7 @@ def plot_selected_probability_fields( node_prob_field['y'].max(), ] ) - vmax = 1 + vmax = 1 + numpy.finfo(float).eps vmin = 0 for lvl in level_list: figure = pyplot.figure()