diff --git a/auxiliary_tools/cdat_regression_testing/880-eamxx/run_e3sm_diags_1996.py b/auxiliary_tools/cdat_regression_testing/880-eamxx/run_e3sm_diags_1996.py new file mode 100644 index 000000000..af15ebae9 --- /dev/null +++ b/auxiliary_tools/cdat_regression_testing/880-eamxx/run_e3sm_diags_1996.py @@ -0,0 +1,37 @@ +import os +from e3sm_diags.parameter.core_parameter import CoreParameter +from e3sm_diags.run import runner + +param = CoreParameter() + +# param.reference_data_path = '/global/cfs/cdirs/e3sm/diagnostics/observations/Atm/climatology' +# param.test_data_path = '/global/cfs/cdirs/e3sm/zhang40/e3sm_diags_for_EAMxx/data/Cess' +# param.reference_data_path = '/global/cfs/cdirs/e3sm/diagnostics/observations/Atm/climatology' +param.reference_data_path = ( + "/global/cfs/cdirs/e3sm/diagnostics/observations/Atm/time-series" +) +param.test_data_path = "/global/cfs/cdirs/e3sm/chengzhu/eamxx/post/data/rgr" +param.test_name = "eamxx_decadal" +param.seasons = ["ANN"] +# param.save_netcdf = True + +param.ref_timeseries_input = True +# Years to slice the ref data, base this off the years in the filenames. +param.ref_start_yr = "1996" +param.ref_end_yr = "1996" + +prefix = "/global/cfs/cdirs/e3sm/www/vo13/tests/eamxx" +param.results_dir = os.path.join(prefix, "eamxx_decadal_1996_1212_edv3") + +runner.sets_to_run = [ + "lat_lon", + "zonal_mean_xy", + "zonal_mean_2d", + "zonal_mean_2d_stratosphere", + "polar", + "cosp_histogram", + "meridional_mean_2d", + "annual_cycle_zonal_mean", +] + +runner.run_diags([param]) diff --git a/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/regression_test.ipynb b/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/regression_test.ipynb new file mode 100644 index 000000000..a1bce678d --- /dev/null +++ b/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/regression_test.ipynb @@ -0,0 +1,1714 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# CDAT Migration Regression Testing Notebook (`.nc` files)\n", + "\n", + "This notebook is used to perform regression testing between the development and\n", + "production versions of a diagnostic set.\n", + "\n", + "## How it works\n", + "\n", + "It compares the relative differences (%) between ref and test variables between\n", + "the dev and `main` branches.\n", + "\n", + "## How to use\n", + "\n", + "PREREQUISITE: The diagnostic set's netCDF stored in `.json` files in two directories\n", + "(dev and `main` branches).\n", + "\n", + "1. Make a copy of this notebook under `auxiliary_tools/cdat_regression_testing/`.\n", + "2. Run `mamba create -n cdat_regression_test -y -c conda-forge \"python<3.12\" xarray netcdf4 dask pandas matplotlib-base ipykernel`\n", + "3. Run `mamba activate cdat_regression_test`\n", + "4. Update `SET_DIR` and `SET_NAME` in the copy of your notebook.\n", + "5. Run all cells IN ORDER.\n", + "6. Review results for any outstanding differences (>=1e-5 relative tolerance).\n", + " - Debug these differences (e.g., bug in metrics functions, incorrect variable references, etc.)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup Code\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import glob\n", + "\n", + "import numpy as np\n", + "import xarray as xr\n", + "from e3sm_diags.derivations.derivations import DERIVED_VARIABLES\n", + "\n", + "SET_NAME = \"qbo\"\n", + "SET_DIR = \"912-qbo-cwt\"\n", + "\n", + "DEV_PATH = f\"/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/{SET_DIR}/{SET_NAME}/**\"\n", + "DEV_GLOB = sorted(glob.glob(DEV_PATH + \"/*.nc\"))\n", + "DEV_NUM_FILES = len(DEV_GLOB)\n", + "\n", + "MAIN_DIR = \"912-qbo-cwt-main\"\n", + "MAIN_PATH = f\"/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/{MAIN_DIR}/{SET_NAME}/**\"\n", + "MAIN_GLOB = sorted(glob.glob(MAIN_PATH + \"/*.nc\"))\n", + "MAIN_NUM_FILES = len(MAIN_GLOB)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def _check_if_files_found():\n", + " if DEV_NUM_FILES == 0 or MAIN_NUM_FILES == 0:\n", + " raise IOError(\n", + " \"No files found at DEV_PATH and/or MAIN_PATH. \"\n", + " f\"Please check {DEV_PATH} and {MAIN_PATH}.\"\n", + " )\n", + "\n", + "\n", + "def _check_if_matching_filecount():\n", + " if DEV_NUM_FILES != MAIN_NUM_FILES:\n", + " raise IOError(\n", + " \"Number of files do not match at DEV_PATH and MAIN_PATH \"\n", + " f\"({DEV_NUM_FILES} vs. {MAIN_NUM_FILES}).\"\n", + " )\n", + "\n", + " print(f\"Matching file count ({DEV_NUM_FILES} and {MAIN_NUM_FILES}).\")\n", + "\n", + "\n", + "def _check_if_missing_files():\n", + " missing_count = 0\n", + "\n", + " for fp_main in MAIN_GLOB:\n", + " fp_dev = fp_main.replace(SET_DIR, \"main\")\n", + "\n", + " if fp_dev not in MAIN_GLOB:\n", + " print(f\"No production file found to compare with {fp_dev}!\")\n", + " missing_count += 1\n", + "\n", + " for fp_dev in DEV_GLOB:\n", + " fp_main = fp_dev.replace(\"main\", SET_DIR)\n", + "\n", + " if fp_main not in DEV_GLOB:\n", + " print(f\"No development file found to compare with {fp_main}!\")\n", + " missing_count += 1\n", + "\n", + " print(f\"Number of files missing: {missing_count}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def _get_relative_diffs():\n", + " # We are mainly focusing on relative tolerance here (in percentage terms).\n", + " atol = 0\n", + " rtol = 1e-5\n", + "\n", + " for fp_main in MAIN_GLOB:\n", + " if \"test.nc\" in fp_main or \"ref.nc\" in fp_main:\n", + " fp_dev = fp_main.replace(\"main-qbo-wavelet\", SET_DIR)\n", + "\n", + " print(\"Comparing:\")\n", + " print(f\" * {fp_dev}\")\n", + " print(f\" * {fp_main}\")\n", + "\n", + " ds1 = xr.open_dataset(fp_dev)\n", + " ds2 = xr.open_dataset(fp_main)\n", + "\n", + " var_keys = [\"U\"]\n", + " for key in var_keys:\n", + " print(f\" * var_key: {key}\")\n", + "\n", + " dev_data = ds1[key].values\n", + " main_data = ds2[key].values\n", + "\n", + " if dev_data is None or main_data is None:\n", + " print(\" * Could not find variable key in the dataset(s)\")\n", + " continue\n", + "\n", + " try:\n", + " np.testing.assert_allclose(\n", + " dev_data,\n", + " main_data,\n", + " atol=atol,\n", + " rtol=rtol,\n", + " )\n", + " except (KeyError, AssertionError) as e:\n", + " print(f\" {e}\")\n", + " else:\n", + " print(f\" * All close and within relative tolerance ({rtol})\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Check for matching and equal number of files\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "_check_if_files_found()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt/qbo/QBO-ERA-Interim/qbo_diags_qbo_ref.nc',\n", + " '/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt/qbo/QBO-ERA-Interim/qbo_diags_qbo_test.nc']" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DEV_GLOB" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt-main/qbo/QBO-ERA-Interim/qbo_diags_qbo_ref.nc',\n", + " '/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt-main/qbo/QBO-ERA-Interim/qbo_diags_qbo_test.nc']" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "MAIN_GLOB" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "No production file found to compare with /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/main-main/qbo/QBO-ERA-Interim/qbo_diags_qbo_ref.nc!\n", + "No production file found to compare with /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/main-main/qbo/QBO-ERA-Interim/qbo_diags_qbo_test.nc!\n", + "Number of files missing: 2\n" + ] + } + ], + "source": [ + "_check_if_missing_files()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Matching file count (2 and 2).\n" + ] + } + ], + "source": [ + "_check_if_matching_filecount()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Let's ignore `qbo_diags_level_ref.nc` and `qbo_diags_level_test.nc`.\n", + "\n", + "- Those files are just the Z dimension of the variable found in the `qbo_diags_qbo_ref.nc` and `qbo_diags_qbo_test.nc`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "MAIN_GLOB = [filename for filename in MAIN_GLOB if \"_level\" not in filename]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2 Compare the netCDF files between branches\n", + "\n", + "- Compare \"ref\" and \"test\" files\n", + "- \"diff\" files are ignored because getting relative diffs for these does not make sense (relative diff will be above tolerance)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Comparing:\n", + " * /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt-main/qbo/QBO-ERA-Interim/qbo_diags_qbo_ref.nc\n", + " * /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt-main/qbo/QBO-ERA-Interim/qbo_diags_qbo_ref.nc\n", + " * var_key: U\n", + " * All close and within relative tolerance (1e-05)\n", + "Comparing:\n", + " * /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt-main/qbo/QBO-ERA-Interim/qbo_diags_qbo_test.nc\n", + " * /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt-main/qbo/QBO-ERA-Interim/qbo_diags_qbo_test.nc\n", + " * var_key: U\n", + " * All close and within relative tolerance (1e-05)\n" + ] + } + ], + "source": [ + "_get_relative_diffs()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Results\n", + "\n", + "- Reference file diffs are massive because the CDAT codebase does not correctly sort the data by the Z axis (`plev`). I opened an issue to address this on `main` here: https://github.com/E3SM-Project/e3sm_diags/issues/825\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Validation: Sorting the CDAT produced reference file by the Z axis in ascending fixes the issue. We can move forward with the changes in this PR.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/global/homes/v/vo13/miniforge3/envs/e3sm_diags_dev_912/lib/python3.12/site-packages/esmpy/interface/loadESMF.py:94: VersionWarning: ESMF installation version 8.8.0, ESMPy version 8.8.0b0\n", + " warnings.warn(\"ESMF installation version {}, ESMPy version {}\".format(\n" + ] + } + ], + "source": [ + "import xcdat as xc\n", + "import xarray as xr\n", + "\n", + "ds_xc = xc.open_dataset(\n", + " \"/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/664-qbo/qbo/QBO-ERA-Interim/qbo_diags_qbo_ref.nc\"\n", + ")\n", + "ds_cdat = xc.open_dataset(\n", + " \"/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/main/qbo/QBO-ERA-Interim/qbo_diags_qbo_ref.nc\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'plev' (plev: 37)> Size: 296B\n",
+       "array([   1.,    2.,    3.,    5.,    7.,   10.,   20.,   30.,   50.,   70.,\n",
+       "        100.,  125.,  150.,  175.,  200.,  225.,  250.,  300.,  350.,  400.,\n",
+       "        450.,  500.,  550.,  600.,  650.,  700.,  750.,  775.,  800.,  825.,\n",
+       "        850.,  875.,  900.,  925.,  950.,  975., 1000.])\n",
+       "Coordinates:\n",
+       "  * plev     (plev) float64 296B 1.0 2.0 3.0 5.0 7.0 ... 925.0 950.0 975.0 1e+03
" + ], + "text/plain": [ + " Size: 296B\n", + "array([ 1., 2., 3., 5., 7., 10., 20., 30., 50., 70.,\n", + " 100., 125., 150., 175., 200., 225., 250., 300., 350., 400.,\n", + " 450., 500., 550., 600., 650., 700., 750., 775., 800., 825.,\n", + " 850., 875., 900., 925., 950., 975., 1000.])\n", + "Coordinates:\n", + " * plev (plev) float64 296B 1.0 2.0 3.0 5.0 7.0 ... 925.0 950.0 975.0 1e+03" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds_xc[\"plev\"]" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'plev' (plev: 37)> Size: 296B\n",
+       "array([1000.,  975.,  950.,  925.,  900.,  875.,  850.,  825.,  800.,  775.,\n",
+       "        750.,  700.,  650.,  600.,  550.,  500.,  450.,  400.,  350.,  300.,\n",
+       "        250.,  225.,  200.,  175.,  150.,  125.,  100.,   70.,   50.,   30.,\n",
+       "         20.,   10.,    7.,    5.,    3.,    2.,    1.])\n",
+       "Coordinates:\n",
+       "  * plev     (plev) float64 296B 1e+03 975.0 950.0 925.0 ... 5.0 3.0 2.0 1.0\n",
+       "Attributes:\n",
+       "    axis:           Z\n",
+       "    units:          hPa\n",
+       "    standard_name:  air_pressure\n",
+       "    long_name:      pressure\n",
+       "    positive:       down\n",
+       "    realtopology:   linear
" + ], + "text/plain": [ + " Size: 296B\n", + "array([1000., 975., 950., 925., 900., 875., 850., 825., 800., 775.,\n", + " 750., 700., 650., 600., 550., 500., 450., 400., 350., 300.,\n", + " 250., 225., 200., 175., 150., 125., 100., 70., 50., 30.,\n", + " 20., 10., 7., 5., 3., 2., 1.])\n", + "Coordinates:\n", + " * plev (plev) float64 296B 1e+03 975.0 950.0 925.0 ... 5.0 3.0 2.0 1.0\n", + "Attributes:\n", + " axis: Z\n", + " units: hPa\n", + " standard_name: air_pressure\n", + " long_name: pressure\n", + " positive: down\n", + " realtopology: linear" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds_cdat[\"plev\"]" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "ds_cdat = ds_cdat.sortby(\"plev\", ascending=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'plev' (plev: 37)> Size: 296B\n",
+       "array([   1.,    2.,    3.,    5.,    7.,   10.,   20.,   30.,   50.,   70.,\n",
+       "        100.,  125.,  150.,  175.,  200.,  225.,  250.,  300.,  350.,  400.,\n",
+       "        450.,  500.,  550.,  600.,  650.,  700.,  750.,  775.,  800.,  825.,\n",
+       "        850.,  875.,  900.,  925.,  950.,  975., 1000.])\n",
+       "Coordinates:\n",
+       "  * plev     (plev) float64 296B 1.0 2.0 3.0 5.0 7.0 ... 925.0 950.0 975.0 1e+03
" + ], + "text/plain": [ + " Size: 296B\n", + "array([ 1., 2., 3., 5., 7., 10., 20., 30., 50., 70.,\n", + " 100., 125., 150., 175., 200., 225., 250., 300., 350., 400.,\n", + " 450., 500., 550., 600., 650., 700., 750., 775., 800., 825.,\n", + " 850., 875., 900., 925., 950., 975., 1000.])\n", + "Coordinates:\n", + " * plev (plev) float64 296B 1.0 2.0 3.0 5.0 7.0 ... 925.0 950.0 975.0 1e+03" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds_xc.plev" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "ename": "AssertionError", + "evalue": "\nNot equal to tolerance rtol=1e-07, atol=0\n\nMismatched elements: 4440 / 4440 (100%)\nMax absolute difference among violations: 0.18704352\nMax relative difference among violations: 7.69507964\n ACTUAL: array([[-16.930712, -42.569729, -25.370665, ..., -2.711329, -2.530502,\n -2.283684],\n [ -2.24533 , -41.558418, -27.585657, ..., -2.62069 , -2.403724,...\n DESIRED: array([[-16.94262 , -42.5402 , -25.36748 , ..., -2.710924, -2.53099 ,\n -2.285837],\n [ -2.284392, -41.54002 , -27.588021, ..., -2.624998, -2.409103,...", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[16], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01mnumpy\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtesting\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43massert_allclose\u001b[49m\u001b[43m(\u001b[49m\u001b[43mds_xc\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mU\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mds_cdat\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mU\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n", + " \u001b[0;31m[... skipping hidden 1 frame]\u001b[0m\n", + "File \u001b[0;32m~/miniforge3/envs/e3sm_diags_dev_912/lib/python3.12/contextlib.py:81\u001b[0m, in \u001b[0;36mContextDecorator.__call__..inner\u001b[0;34m(*args, **kwds)\u001b[0m\n\u001b[1;32m 78\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(func)\n\u001b[1;32m 79\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21minner\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwds):\n\u001b[1;32m 80\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_recreate_cm():\n\u001b[0;32m---> 81\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwds\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/miniforge3/envs/e3sm_diags_dev_912/lib/python3.12/site-packages/numpy/testing/_private/utils.py:885\u001b[0m, in \u001b[0;36massert_array_compare\u001b[0;34m(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf, strict, names)\u001b[0m\n\u001b[1;32m 880\u001b[0m err_msg \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m'\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;241m.\u001b[39mjoin(remarks)\n\u001b[1;32m 881\u001b[0m msg \u001b[38;5;241m=\u001b[39m build_err_msg([ox, oy], err_msg,\n\u001b[1;32m 882\u001b[0m verbose\u001b[38;5;241m=\u001b[39mverbose, header\u001b[38;5;241m=\u001b[39mheader,\n\u001b[1;32m 883\u001b[0m names\u001b[38;5;241m=\u001b[39mnames,\n\u001b[1;32m 884\u001b[0m precision\u001b[38;5;241m=\u001b[39mprecision)\n\u001b[0;32m--> 885\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mAssertionError\u001b[39;00m(msg)\n\u001b[1;32m 886\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m:\n\u001b[1;32m 887\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01mtraceback\u001b[39;00m\n", + "\u001b[0;31mAssertionError\u001b[0m: \nNot equal to tolerance rtol=1e-07, atol=0\n\nMismatched elements: 4440 / 4440 (100%)\nMax absolute difference among violations: 0.18704352\nMax relative difference among violations: 7.69507964\n ACTUAL: array([[-16.930712, -42.569729, -25.370665, ..., -2.711329, -2.530502,\n -2.283684],\n [ -2.24533 , -41.558418, -27.585657, ..., -2.62069 , -2.403724,...\n DESIRED: array([[-16.94262 , -42.5402 , -25.36748 , ..., -2.710924, -2.53099 ,\n -2.285837],\n [ -2.284392, -41.54002 , -27.588021, ..., -2.624998, -2.409103,..." + ] + } + ], + "source": [ + "import numpy as np\n", + "\n", + "np.testing.assert_allclose(ds_xc[\"U\"], ds_cdat[\"U\"])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compare Maxes and Mins -- Really close\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "61.54721945135814 61.36017592984254\n", + "-66.54760399615296 -66.52449748057968\n" + ] + } + ], + "source": [ + "print(ds_xc[\"U\"].max().item(), ds_cdat[\"U\"].max().item())\n", + "print(ds_xc[\"U\"].min().item(), ds_cdat[\"U\"].min().item())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compare Sum and Mean -- Really close\n" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-3.7398468780963863 -3.7455298743231142\n", + "-16604.92013874794 -16630.15264199463\n" + ] + } + ], + "source": [ + "print(ds_xc[\"U\"].mean().item(), ds_cdat[\"U\"].mean().item())\n", + "print(ds_xc[\"U\"].sum().item(), ds_cdat[\"U\"].sum().item())" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "e3sm_diags_dev_912", + "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.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/regression_test_png.ipynb b/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/regression_test_png.ipynb new file mode 100644 index 000000000..7b8f57abb --- /dev/null +++ b/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/regression_test_png.ipynb @@ -0,0 +1,198 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# CDAT Migration Regression Testing Notebook (`.png` files)\n", + "\n", + "This notebook is used to perform regression testing between the development and\n", + "production versions of a diagnostic set.\n", + "\n", + "## How to use\n", + "\n", + "PREREQUISITE: The diagnostic set's netCDF stored in `.json` files in two directories\n", + "(dev and `main` branches).\n", + "\n", + "1. Make a copy of this notebook under `auxiliary_tools/cdat_regression_testing/`.\n", + "2. Run `mamba create -n cdat_regression_test -y -c conda-forge \"python<3.12\" xarray netcdf4 dask pandas matplotlib-base ipykernel`\n", + "3. Run `mamba activate cdat_regression_test`\n", + "4. Update `SET_DIR` and `SET_NAME` in the copy of your notebook.\n", + "5. Run all cells IN ORDER.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup Code\n" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "import glob\n", + "\n", + "from auxiliary_tools.cdat_regression_testing.utils import get_image_diffs\n", + "\n", + "SET_NAME = \"qbo\"\n", + "SET_DIR = \"912-qbo-cwt\"\n", + "\n", + "DEV_PATH = f\"/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/{SET_DIR}/{SET_NAME}/**\"\n", + "DEV_GLOB = sorted(glob.glob(DEV_PATH + \"/*.png\"))\n", + "DEV_NUM_FILES = len(DEV_GLOB)\n", + "\n", + "MAIN_DIR = \"912-qbo-cwt-main\"\n", + "MAIN_PATH = f\"/global/cfs/cdirs/e3sm/www/cdat-migration-fy24/{MAIN_DIR}/{SET_NAME}/**\"\n", + "MAIN_GLOB = sorted(glob.glob(MAIN_PATH + \"/*.png\"))\n", + "MAIN_NUM_FILES = len(MAIN_GLOB)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "def _check_if_files_found():\n", + " if DEV_NUM_FILES == 0 or MAIN_NUM_FILES == 0:\n", + " raise IOError(\n", + " \"No files found at DEV_PATH and/or MAIN_PATH. \"\n", + " f\"Please check {DEV_PATH} and {MAIN_PATH}.\"\n", + " )\n", + "\n", + "\n", + "def _check_if_matching_filecount():\n", + " if DEV_NUM_FILES != MAIN_NUM_FILES:\n", + " raise IOError(\n", + " \"Number of files do not match at DEV_PATH and MAIN_PATH \"\n", + " f\"({DEV_NUM_FILES} vs. {MAIN_NUM_FILES}).\"\n", + " )\n", + "\n", + " print(f\"Matching file count ({DEV_NUM_FILES} and {MAIN_NUM_FILES}).\")\n", + "\n", + "\n", + "def _check_if_missing_files():\n", + " missing_count = 0\n", + "\n", + " for fp_main in MAIN_GLOB:\n", + " fp_dev = fp_main.replace(MAIN_DIR, SET_DIR)\n", + "\n", + " if fp_dev not in DEV_GLOB:\n", + " print(f\"No development file found to compare with {fp_dev}!\")\n", + " missing_count += 1\n", + "\n", + " for fp_dev in DEV_GLOB:\n", + " if \"diff\" not in fp_dev:\n", + " fp_main = fp_dev.replace(SET_DIR, MAIN_DIR)\n", + "\n", + " if fp_main not in MAIN_GLOB:\n", + " print(f\"No development file found to compare with {fp_main}!\")\n", + " missing_count += 1\n", + "\n", + " print(f\"Number of files missing: {missing_count}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Check for matching and equal number of files\n" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "_check_if_files_found()" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of files missing: 0\n" + ] + } + ], + "source": [ + "_check_if_missing_files()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2 Compare the plots between branches\n", + "\n", + "- Compare \"ref\" and \"test\" files\n", + "- \"diff\" files are ignored because getting relative diffs for these does not make sense (relative diff will be above tolerance)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Comparing:\n", + " * /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt-main/qbo/QBO-ERA-Interim/qbo_diags.png\n", + " * /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt/qbo/QBO-ERA-Interim/qbo_diags.png\n", + " * Difference path /global/cfs/cdirs/e3sm/www/cdat-migration-fy24/912-qbo-cwt/qbo/QBO-ERA-Interim_diff/qbo_diags.png\n" + ] + } + ], + "source": [ + "for main_path, dev_path in zip(MAIN_GLOB, DEV_GLOB):\n", + " print(\"Comparing:\")\n", + " print(f\" * {main_path}\")\n", + " print(f\" * {dev_path}\")\n", + "\n", + " get_image_diffs(dev_path, main_path)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Results\n", + "\n", + "- There is a large difference with the QBO wavelet plot.\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "e3sm_diags_dev_912", + "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.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/run_script.py b/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/run_script.py new file mode 100644 index 000000000..e584d2caf --- /dev/null +++ b/auxiliary_tools/cdat_regression_testing/912-qbo-wavelet-cwt-fix/run_script.py @@ -0,0 +1,10 @@ +# %% +from auxiliary_tools.cdat_regression_testing.base_run_script import run_set + +SET_NAME = "qbo" +SET_DIR = "912-qbo-cwt" +CFG_PATH: str | None = None +MULTIPROCESSING = False + +# %% +run_set(SET_NAME, SET_DIR, CFG_PATH, MULTIPROCESSING) diff --git a/conda-env/ci.yml b/conda-env/ci.yml index 511df6df5..66f1d4250 100644 --- a/conda-env/ci.yml +++ b/conda-env/ci.yml @@ -23,7 +23,8 @@ dependencies: - netcdf4 - output_viewer >=1.3.0 - numpy >=2.0.0,<3.0.0 - - scipy <1.15 + - pywavelets + - scipy - shapely >=2.0.0,<3.0.0 - xarray >=2024.3.0 - xcdat >=0.7.3,<1.0.0 diff --git a/conda-env/dev.yml b/conda-env/dev.yml index 9c75f377a..8ebceabb9 100644 --- a/conda-env/dev.yml +++ b/conda-env/dev.yml @@ -21,7 +21,8 @@ dependencies: - netcdf4 - output_viewer >=1.3.0 - numpy >=2.0.0,<3.0.0 - - scipy <1.15 + - pywavelets + - scipy - shapely >=2.0.0,<3.0.0 - xarray >=2024.3.0 - xcdat >=0.7.3,<1.0.0 diff --git a/e3sm_diags/driver/qbo_driver.py b/e3sm_diags/driver/qbo_driver.py index 3379f4c46..cbf59d253 100644 --- a/e3sm_diags/driver/qbo_driver.py +++ b/e3sm_diags/driver/qbo_driver.py @@ -5,6 +5,7 @@ from typing import TYPE_CHECKING, Dict, Literal, Tuple, TypedDict import numpy as np +import pywt import scipy.fftpack import xarray as xr import xcdat as xc @@ -412,10 +413,10 @@ def _get_psd_from_wavelet(data: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: """ deg = 6 period = np.arange(1, 55 + 1) - freq = 1 / period + widths = deg / (2 * np.pi / period) - widths = deg / (2 * np.pi * freq) - cwtmatr = scipy.signal.cwt(data, scipy.signal.morlet2, widths=widths, w=deg) - psd = np.mean(np.square(np.abs(cwtmatr)), axis=1) + [cfs, freq] = pywt.cwt(data, scales=widths, wavelet="cmor1.5-1.0") + psd = np.mean(np.square(np.abs(cfs)), axis=1) + period = 1 / freq return (period, psd) diff --git a/pyproject.toml b/pyproject.toml index df177d55d..e8a80e39a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,7 +41,8 @@ dependencies = [ "netcdf4", "output_viewer >=1.3.0", "numpy >=2.0.0,<3.0.0", - "scipy <1.15", + "pywavelets", + "scipy", "shapely >=2.0.0,<3.0.0", "xarray >=2024.03.0", "xcdat >=0.7.3,<1.0.0",