From 2129ea0454033e68f945e172819c79cc4fe71f3c Mon Sep 17 00:00:00 2001 From: Hauke Schulz <43613877+observingClouds@users.noreply.github.com> Date: Tue, 16 Apr 2024 15:05:29 -0700 Subject: [PATCH] WIP --- .github/workflows/ci.yaml | 4 +- docs/ArtificialInformation_Filter.ipynb | 335 ++++++++++++++++++++++++ docs/environment.yml | 6 +- docs/index.rst | 10 +- tests/test_get_keepbits.py | 140 ++++++++++ xbitinfo/xbitinfo.py | 136 +++++++++- 6 files changed, 624 insertions(+), 7 deletions(-) create mode 100644 docs/ArtificialInformation_Filter.ipynb diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 5669aaf0..1a877e3e 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -96,11 +96,11 @@ jobs: auto-update-conda: false channels: conda-forge miniforge-variant: Mambaforge - activate-environment: bitinfo + activate-environment: bitinfo-docs python-version: '3.11' - name: Set up conda environment run: | - mamba env update -f environment.yml + mamba env update -f docs/environment.yml - name: Remove julia (issue #212) run: | conda remove julia diff --git a/docs/ArtificialInformation_Filter.ipynb b/docs/ArtificialInformation_Filter.ipynb new file mode 100644 index 00000000..c9f8d4cc --- /dev/null +++ b/docs/ArtificialInformation_Filter.ipynb @@ -0,0 +1,335 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "38ac6d1a", + "metadata": {}, + "source": [ + "**<<<<<<< local**" + ] + }, + { + "cell_type": "markdown", + "id": "895e1d5a", + "metadata": {}, + "source": [ + "# Artificial information filtering\n", + "\n", + "In simple terms the bitinformation is retrieved by checking how variable a bit pattern is. However, this approach cannot distinguish between actual information content and artifical information content. By studying the distribution of the information content the user can often identify clear cut-offs of real information content and artificial information content.\n", + "\n", + "The following example shows how such a separation of real information and artificial information can look like. To do so, artificial information is artificially added to an example dataset by applying linear quantization. Linear quantization is often applied to climate datasets (e.g. ERA5) and needs to be accounted for in order to retrieve meaningful bitinformation content. An algorithm that aims at detecting this artificial information itself is introduced." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3c37dd36", + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import xbitinfo as xb\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "e8e1424f", + "metadata": {}, + "source": [ + "## Loading example dataset\n", + "We use here the openly accessible CONUS dataset. The dataset is available at full precision." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b18b9e24", + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_zarr(\n", + " \"s3://hytest/conus404/conus404_hourly.zarr\",\n", + " storage_options={\n", + " \"anon\": True,\n", + " \"requester_pays\": False,\n", + " \"client_kwargs\": {\"endpoint_url\": \"https://usgs.osn.mghpcc.org\"},\n", + " },\n", + ")\n", + "# selecting water vapor mixing ratio at 2 meters\n", + "data = ds[\"ACSWUPB\"]\n", + "# select subset of data for demonstration purposes\n", + "chunk = data.isel(time=slice(0, 9), y=slice(0, 525), x=slice(0, 525))\n", + "chunk" + ] + }, + { + "cell_type": "markdown", + "id": "535ce421", + "metadata": {}, + "source": [ + "## Creating dataset copy with artificial information\n", + "### Functions to encode and decode" + ] + }, + { + "cell_type": "markdown", + "id": "69543b4c", + "metadata": {}, + "source": [ + "**=======**" + ] + }, + { + "cell_type": "markdown", + "id": "1842f792", + "metadata": {}, + "source": [ + "# Artificial information filtering\n", + "\n", + "In simple terms the bitinformation is retrieved by checking how variable a bit pattern is. However, this approach cannot distinguish between actual information content and artifical information content. By studying the distribution of the information content the user can often identify clear cut-offs of real information content and artificial information content.\n", + "\n", + "The following example shows how such a separation of real information and artificial information can look like. To do so, artificial information is artificially added to an example dataset by applying linear quantization. Linear quantization is often applied to climate datasets (e.g. ERA5) and needs to be accounted for in order to retrieve meaningful bitinformation content. An algorithm that aims at detecting this artificial information itself is introduced." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb998fbb", + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import xbitinfo as xb\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "32ac97e0", + "metadata": {}, + "source": [ + "## Loading example dataset\n", + "We use here the openly accessible CONUS dataset. The dataset is available at full precision." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9639a618", + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_zarr(\n", + " \"s3://hytest/conus404/conus404_monthly.zarr\",\n", + " storage_options={\n", + " \"anon\": True,\n", + " \"requester_pays\": False,\n", + " \"client_kwargs\": {\"endpoint_url\": \"https://usgs.osn.mghpcc.org\"},\n", + " },\n", + ")\n", + "# selecting water vapor mixing ratio at 2 meters\n", + "data = ds[\"ACSWDNT\"]\n", + "# select subset of data for demonstration purposes\n", + "chunk = data.isel(time=slice(0, 2), y=slice(0, 1015), x=slice(0, 1050))\n", + "chunk" + ] + }, + { + "cell_type": "markdown", + "id": "3d735e4b", + "metadata": {}, + "source": [ + "## Creating dataset copy with artificial information\n", + "### Functions to encode and decode" + ] + }, + { + "cell_type": "markdown", + "id": "0d30feaa", + "metadata": {}, + "source": [ + "**>>>>>>> remote**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3a7c7ae", + "metadata": {}, + "outputs": [], + "source": [ + "# Encoding function to compress data\n", + "def encode(chunk, scale, offset, dtype, astype):\n", + " enc = (chunk - offset) * scale\n", + " enc = np.around(enc)\n", + " enc = enc.astype(astype, copy=False)\n", + " return enc\n", + "\n", + "\n", + "# Decoding function to decompress data\n", + "def decode(enc, scale, offset, dtype, astype):\n", + " dec = (enc / scale) + offset\n", + " dec = dec.astype(dtype, copy=False)\n", + " return dec" + ] + }, + { + "cell_type": "markdown", + "id": "fa6f26c7", + "metadata": {}, + "source": [ + "### Transform dataset to introduce artificial information" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c09e3cf3", + "metadata": {}, + "outputs": [], + "source": [ + "xmin = np.min(chunk)\n", + "xmax = np.max(chunk)\n", + "scale = (2**16 - 1) / (xmax - xmin)\n", + "offset = xmin\n", + "enc = encode(chunk, scale, offset, \"f4\", \"u2\")\n", + "dec = decode(enc, scale, offset, \"f4\", \"u2\")" + ] + }, + { + "cell_type": "markdown", + "id": "7126810d", + "metadata": {}, + "source": [ + "## Comparison of bitinformation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05ef8a94", + "metadata": {}, + "outputs": [], + "source": [ + "# original dataset without artificial information\n", + "orig_info = xb.get_bitinformation(\n", + " xr.Dataset({\"w/o artif. info\": chunk}),\n", + " dim=\"x\",\n", + " implementation=\"python\",\n", + ")\n", + "\n", + "# dataset with artificial information\n", + "arti_info = xb.get_bitinformation(\n", + " xr.Dataset({\"w artif. info\": dec}),\n", + " dim=\"x\",\n", + " implementation=\"python\",\n", + ")\n", + "\n", + "# plotting distribution of bitwise information content\n", + "info = xr.merge([orig_info, arti_info])\n", + "plot = xb.plot_bitinformation(info)" + ] + }, + { + "cell_type": "markdown", + "id": "de1ecb7e", + "metadata": {}, + "source": [ + "The figure reveals that artificial information is introduced by applying linear quantization. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8600d4b8", + "metadata": {}, + "outputs": [], + "source": [ + "keepbits = xb.get_keepbits(info, inflevel=[0.99])\n", + "print(\n", + " f\"The number of keepbits increased from {keepbits['w/o artif. info'].item(0)} bits in the original dataset to {keepbits['w artif. info'].item(0)} bits in the dataset with artificial information.\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "fa80f988", + "metadata": {}, + "source": [ + "In the following, a gradient based filter is introduced to remove this artificial information again so that even in case artificial information is present in a dataset the number of keepbits remains similar." + ] + }, + { + "cell_type": "markdown", + "id": "3f7a7c2e", + "metadata": {}, + "source": [ + "## Artificial information filter\n", + "The filter `gradient` works as follows:\n", + "\n", + "1. It determines the Cumulative Distribution Function(CDF) of the bitwise information content\n", + "2. It computes the gradient of the CDF to identify points where the gradient becomes close to a given tolerance indicating a drop in information.\n", + "3. Simultaneously, it keeps track of the minimum cumulative sum of information content which is threshold here, which signifies at least this much fraction of total information needs to be passed.\n", + "4. So the bit where the intersection of the gradient reaching the tolerance and the cumulative sum exceeding the threshold is our TrueKeepbits. All bits beyond this index are assumed to contain artificial information and are set to zero in order to cut them off.\n", + "5. You can see the above concept implemented in the function get_cdf_without_artificial_information in xbitinfo.py\n", + "\n", + "Please note that this filter relies on a clear separation between real and artificial information content and might not work in all cases." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0ab6633", + "metadata": {}, + "outputs": [], + "source": [ + "xb.get_keepbits(\n", + " arti_info,\n", + " inflevel=[0.99],\n", + " information_filter=\"Gradient\",\n", + " **{\"threshold\": 0.7, \"tolerance\": 0.001}\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "21c6369d", + "metadata": {}, + "source": [ + "With the application of the filter the keepbits are closer/identical to their original value in the dataset without artificial information. The plot of the bitinformation visualizes this:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8e9183b2", + "metadata": {}, + "outputs": [], + "source": [ + "plot = xb.plot_bitinformation(arti_info, information_filter=\"Gradient\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/environment.yml b/docs/environment.yml index 23e08aa9..10d2f483 100644 --- a/docs/environment.yml +++ b/docs/environment.yml @@ -2,7 +2,9 @@ name: bitinfo-docs channels: - conda-forge dependencies: - - python=3.9 + - python + - julia<1.9.0 + - pyjulia - matplotlib-base - numpy - pooch @@ -21,6 +23,8 @@ dependencies: - sphinx-book-theme>=0.1.7 - myst-nb - numcodecs>=0.10.0 + - intake-xarray + - s3fs - pip - pip: - -e ../. diff --git a/docs/index.rst b/docs/index.rst index 4b38309f..de27d775 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -96,7 +96,6 @@ Credits quick-start.ipynb - **User Guide** * :doc:`chunking` @@ -108,6 +107,15 @@ Credits chunking.ipynb +* :doc:`artificialinformation` + +.. toctree:: + :maxdepth: 1 + :hidden: + :caption: User Guide + + ArtificialInformation_Filter.ipynb + **Help & Reference** * :doc:`api` diff --git a/tests/test_get_keepbits.py b/tests/test_get_keepbits.py index cb24cd13..9525354a 100644 --- a/tests/test_get_keepbits.py +++ b/tests/test_get_keepbits.py @@ -30,3 +30,143 @@ def test_get_keepbits_inflevel_dim(rasm_info_per_bit, inflevel): if isinstance(inflevel, (int, float)): inflevel = [inflevel] assert (keepbits.inflevel == inflevel).all() + + +def test_get_keepbits_informationFilter(): + """ + Test the `get_keepbits` function with different information filters. + + This test function checks the behavior of the `get_keepbits` function when applying gradient information filter. + The dataset contains artificial information and thus applying the filter should result in lesser number of bits + than what should be when filter is None. + + + Raises: + AssertionError: If the test conditions are not met. + + """ + + bit32_values = [ + "±", + "e1", + "e2", + "e3", + "e4", + "e5", + "e6", + "e7", + "e8", + "m1", + "m2", + "m3", + "m4", + "m5", + "m6", + "m7", + "m8", + "m9", + "m10", + "m11", + "m12", + "m13", + "m14", + "m15", + "m16", + "m17", + "m18", + "m19", + "m20", + "m21", + "m22", + "m23", + ] + data_variable = xr.DataArray( + data=[ + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 0.00000000e00, + 1.11799129e-01, + 8.19114977e-01, + 4.41578500e-01, + 3.25470303e-01, + 4.35195738e-01, + 2.81462993e-01, + 2.10719742e-01, + 1.46638224e-01, + 9.24534031e-02, + 4.41090879e-02, + 1.13842504e-02, + 8.20088050e-04, + 2.62239097e-06, + 7.11284508e-07, + 1.18183485e-06, + 9.49338973e-09, + 1.80859255e-07, + 7.72662891e-07, + 1.37865391e-05, + 2.11117224e-06, + 2.01353088e-07, + 3.20755770e-02, + 6.06721012e-03, + 2.25987148e-04, + 1.71530452e-06, + 5.13067595e-03, + ], + coords={"bitfloat32": bit32_values, "dim": "x"}, + dims=["bitfloat32"] + ) + info_ds = xr.Dataset({"RH2": data_variable}) + Keepbits_FilterNone = xb.get_keepbits( + info_ds, + inflevel=[0.99], + information_filter=None, + **{"threshold": 0.7, "tolerance": 0.001} + ) + Keepbits_FilterNone_Value = Keepbits_FilterNone["RH2"].values + assert Keepbits_FilterNone_Value == 19 + + Keepbits_FilterGradient = xb.get_keepbits( + info_ds, + inflevel=[0.99], + information_filter="Gradient", + **{"threshold": 0.7, "tolerance": 0.001} + ) + Keepbits_FilterGradient_Value = Keepbits_FilterGradient["RH2"].values + assert Keepbits_FilterGradient_Value == 7 + + +def test_get_keepbits_informationFilter_1(): + """ + Test the `get_keepbits` function with different information filters. + + This test function checks the behavior of the `get_keepbits` function when applying gradient information filter. + The dataset does not contain artificial information and thus the number of keepbits when gradient filter is applied + should be equal to when filter is None. + + Raises: + AssertionError: If the test conditions are not met. + + """ + + ds = xr.tutorial.load_dataset("air_temperature") + info = xb.get_bitinformation(ds, dim="lat") + Keepbits_FilterNone = xb.get_keepbits( + info, + inflevel=[0.99], + information_filter=None, + **{"threshold": 0.7, "tolerance": 0.001} + ) + Keepbits_FilterNone_Value = Keepbits_FilterNone["air"].values + + Keepbits_FilterGradient = xb.get_keepbits( + info, + inflevel=[0.99], + information_filter="Gradient", + **{"threshold": 0.7, "tolerance": 0.001} + ) + + Keepbits_FilterGradient_Value = Keepbits_FilterGradient["air"].values + assert Keepbits_FilterNone_Value == Keepbits_FilterGradient_Value diff --git a/xbitinfo/xbitinfo.py b/xbitinfo/xbitinfo.py index 9a98f408..10cf1d42 100644 --- a/xbitinfo/xbitinfo.py +++ b/xbitinfo/xbitinfo.py @@ -385,7 +385,121 @@ def load_bitinformation(label): raise FileNotFoundError(f"No bitinformation could be found at {label+'.json'}") -def get_keepbits(info_per_bit, inflevel=0.99): +def get_cdf_without_artificial_information( + info_per_bit, bitdim, threshold, tolerance, bit_vars +): + """ + Calculate a Cumulative Distribution Function (CDF) with artificial information removal. + + This function calculates a modified CDF for a given set of bit information and variable dimensions, + removing artificial information while preserving the desired threshold of information content. + + 1.)The function's aim is to return the cdf in a way that artificial information gets removed. + 2.)This function calculates the CDF using the provided information content per bit dataset. + 3.)It then computes the gradient of the CDF values to identify points where the gradient becomes close to the given tolerance, + indicating a drop in information. + 4.)Simultaneously, it keeps track of the minimum cumulative sum of information content which is threshold here, which signifies atleast + this much fraction of total information needs to be passed. + 5.)So the bit where the intersection of the gradient reaching the tolerance and the cumulative sum exceeding the threshold. All bits beyond this + index are assumed to contain artificial information and are set to zero in the resulting CDF. + + + Parameters: + ----------- + info_per_bit : :py:class: 'xarray.Dataset' + Information content of each bit. This is the output from :py:func:`xbitinfo.xbitinfo.get_bitinformation`. + bitdim : str + The dimension representing the bit information. + threshold : float + Minimum cumulative sum of information content before artificial information filter is applied. + tolerance : float + The tolerance is the value below which gradient starts becoming constant + bit_vars : list + List of variable names of the dataset. + + Returns: + -------- + xarray.Dataset + A modified CDF dataset with artificial information removed. + + Example: + -------- + >>> ds = xr.tutorial.load_dataset("air_temperature") + >>> info = xb.get_bitinformation(ds) + >>> get_keepbits( + ... info, + ... inflevel=[0.99], + ... information_filter="Gradient", + ... **{"threshold": 0.7, "tolerance": 0.001} + ... ) + Size: 80B + Dimensions: (dim: 3, inflevel: 1) + Coordinates: + * dim (dim) = threshold * infSum: + infbits = i + break + + for i in range(0, infbits + 1): + # Normalize CDF values for elements up to 'infbits'. + cdf_array[i] = cdf_array[i] / cdf_array[infbits] + + cdf_array[(infbits + 1) :] = 1 + return cdf + + +def get_keepbits(info_per_bit, inflevel=0.99, information_filter=None, **kwargs): """Get the number of mantissa bits to keep. To be used in :py:func:`xbitinfo.bitround.xr_bitround` and :py:func:`xbitinfo.bitround.jl_bitround`. Parameters @@ -395,6 +509,13 @@ def get_keepbits(info_per_bit, inflevel=0.99): inflevel : float or list Level of information that shall be preserved. + Kwargs + threshold(` `float ``) : defaults to ``0.7`` + Minimum cumulative sum of information content before artificial information filter is applied. + tolerance(` `float ``) : defaults to ``0.001`` + The tolerance is the value below which gradient starts becoming constant + + Returns ------- keepbits : dict @@ -458,7 +579,16 @@ def get_keepbits(info_per_bit, inflevel=0.99): # get only variables of bitdim bit_vars = [v for v in info_per_bit.data_vars if bitdim in info_per_bit[v].dims] if bit_vars != []: - cdf = _cdf_from_info_per_bit(info_per_bit[bit_vars], bitdim) + if information_filter == "Gradient": + cdf = get_cdf_without_artificial_information( + info_per_bit[bit_vars], + bitdim, + kwargs["threshold"], + kwargs["tolerance"], + bit_vars, + ) + else: + cdf = _cdf_from_info_per_bit(info_per_bit[bit_vars], bitdim) data_type = np.dtype(bitdim.replace("bit", "")) n_bits, _, _, n_mant = bit_partitioning(data_type) bitdim_non_mantissa_bits = n_bits - n_mant @@ -701,7 +831,7 @@ class JsonCustomEncoder(json.JSONEncoder): def default(self, obj): if isinstance(obj, (np.ndarray, np.number)): return obj.tolist() - elif isinstance(obj, complex): + elif isinstance(obj, (complex, np.complex)): return [obj.real, obj.imag] elif isinstance(obj, set): return list(obj)