RESULT05_performance_with_number_of_members.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4795a486-c7a9-4ccf-a3c3-937b93553e2f",
   "metadata": {},
   "source": [
    "# Brier Scores with varying ensemble member sizes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2ebf2c51-59f0-4726-8c35-5c8c3d3aa4b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# general tools\n",
    "import os\n",
    "import sys\n",
    "import time\n",
    "import h5py\n",
    "# import zarr\n",
    "import numba as nb\n",
    "from glob import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b9290788-e3d9-4a03-9902-b6392707d62b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "sys.path.insert(0, '/glade/u/home/ksha/GAN_proj/')\n",
    "sys.path.insert(0, '/glade/u/home/ksha/GAN_proj/libs/')\n",
    "\n",
    "from namelist import *\n",
    "import data_utils as du\n",
    "import verif_utils as vu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d0472ee8-fe68-4224-b21a-d8ef0ef30de9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "19890738-dfa8-4036-8edf-a6e2118c407f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bs(y_true, y_pred):\n",
    "    return np.nanmean((y_true - y_pred)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a3ae0362-790b-48bf-8b29-48959c334954",
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(save_dir+'CCPA_domain.hdf', 'r') as h5io:\n",
    "    lon_CCPA = h5io['lon_CCPA'][...]\n",
    "    lat_CCPA = h5io['lat_CCPA'][...]\n",
    "    land_mask_CCPA = h5io['land_mask_CCPA'][...]\n",
    "\n",
    "land_mask = land_mask_CCPA > 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "03c7bc8c-7273-48b6-b9cd-959d55f77f9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'prob_by_ens_size3.hdf'\n",
    "with h5py.File(result_dir+filename, 'r') as h5io:\n",
    "    LDM_prob_p99 = h5io['LDM_prob_p99'][...]\n",
    "    LDM_prob_40 = h5io['LDM_prob_40'][...]\n",
    "\n",
    "LDM_prob_p99 = LDM_prob_p99[:, :, :, land_mask]\n",
    "LDM_prob_40 = LDM_prob_40[:, :, :, land_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fed3f537-f106-4478-9a99-7d52502365c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(camp_dir+'CCPA/CCPA_lead_y2021.hdf', 'r') as h5io:\n",
    "    CCPA_true = h5io['CCPA_lead'][...]\n",
    "\n",
    "with h5py.File(camp_dir+'CCPA/CCPA_CDFs_lead_2002_2019.hdf', 'r') as h5io:\n",
    "    CCPA_CDFs = h5io['CCPA_CDFs_lead'][...]\n",
    "q99 = CCPA_CDFs[99, ...]\n",
    "\n",
    "CCPA_true_40 = CCPA_true > 40\n",
    "CCPA_true_40 = CCPA_true_40[:, :, land_mask]\n",
    "\n",
    "CCPA_true_99 = CCPA_true > q99[None, ...]\n",
    "CCPA_true_99 = CCPA_true_99[:, :, land_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41150f2b-01b6-4230-be3f-964ff405bb3d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "463ad865-af78-4695-ad5e-2aea0f4af0b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# hist_bins = np.linspace(0, 1, 19)\n",
    "\n",
    "# y_true = CCPA_true_40[:, 0:1, :].ravel()\n",
    "# y_pred = LDM_prob_40[1, :, 0:1, ...].ravel()\n",
    "\n",
    "# output_bundle = vu.reliability_diagram_bootstrap(y_true, y_pred, N_boost=1, hist_bins=hist_bins)\n",
    "# prob_true, prob_pred, hist_bins_, use_, o_bar, prob_pred_mean, prob_true_mean = output_bundle\n",
    "\n",
    "# fake_x = np.linspace(0, 1, 10)\n",
    "# plt.figure(figsize=(6, 6))\n",
    "\n",
    "# plt.plot(prob_pred, prob_true, 'r-', lw=1)\n",
    "# plt.plot(prob_pred_mean, prob_true_mean, 'k-', lw=3)\n",
    "\n",
    "# plt.plot(fake_x, fake_x, 'k--')\n",
    "# plt.xlim([0, 1])\n",
    "# plt.ylim([0, 1])\n",
    "# plt.xlabel('Predicted probability', fontsize=14)\n",
    "# plt.ylabel('Observed relative frequency', fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "64cefdf4-8f38-4835-a456-466c90954f79",
   "metadata": {},
   "outputs": [],
   "source": [
    "BS_40_0_8 = np.empty(8)\n",
    "BS_99_0_8 = np.empty(8)\n",
    "\n",
    "BS_40_8_16 = np.empty(8)\n",
    "BS_99_8_16 = np.empty(8)\n",
    "\n",
    "BS_40_16_24 = np.empty(8)\n",
    "BS_99_16_24 = np.empty(8)\n",
    "\n",
    "for i in range(8):\n",
    "    LDM_prob_40_ = LDM_prob_40[i, ...]\n",
    "    LDM_prob_p99_ = LDM_prob_p99[i, ...]\n",
    "    \n",
    "    BS_40_0_8[i] = bs(CCPA_true_40[:, 0:8, :].ravel(), LDM_prob_40_[:, 0:8, :].ravel())\n",
    "    BS_99_0_8[i] = bs(CCPA_true_99[:, 0:8, :].ravel(), LDM_prob_p99_[:, 0:8, :].ravel())\n",
    "\n",
    "    BS_40_8_16[i] = bs(CCPA_true_40[:, 8:16, :].ravel(), LDM_prob_40_[:, 8:16, :].ravel())\n",
    "    BS_99_8_16[i] = bs(CCPA_true_99[:, 8:16, :].ravel(), LDM_prob_p99_[:, 8:16, :].ravel())\n",
    "\n",
    "    BS_40_16_24[i] = bs(CCPA_true_40[:, 16:24, :].ravel(), LDM_prob_40_[:, 16:24, :].ravel())\n",
    "    BS_99_16_24[i] = bs(CCPA_true_99[:, 16:24, :].ravel(), LDM_prob_p99_[:, 16:24, :].ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "33c83dfc-c25c-40a8-b1fb-49917c19d342",
   "metadata": {},
   "outputs": [],
   "source": [
    "BS_40 = np.concatenate((BS_40_0_8[:, None], \n",
    "                        BS_40_8_16[:, None], \n",
    "                        BS_40_16_24[:, None]), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "882de61f-34fc-4b12-9512-abc128fa601b",
   "metadata": {},
   "outputs": [],
   "source": [
    "BS_99 = np.concatenate((BS_99_0_8[:, None], \n",
    "                        BS_99_8_16[:, None], \n",
    "                        BS_99_16_24[:, None]), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "07ae4cc9-75b0-4163-96e7-870dceb2901a",
   "metadata": {},
   "outputs": [],
   "source": [
    "save_dict = {}\n",
    "save_dict['BS_40'] = BS_40\n",
    "save_dict['BS_99'] = BS_99\n",
    "#np.save(result_dir+'ENS_num_verif.npy', save_dict, allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "4123a31e-2b36-4156-b4c3-60c64fb29158",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.plot(BS_40_0_8_)\n",
    "#plt.plot(BS_40_8_16_)\n",
    "#plt.plot(BS_40_16_24_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "03830534-1ea3-4fb5-a35b-e7ff7d9ae504",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.plot(BS_99_0_8_)\n",
    "#plt.plot(BS_99_8_16_)\n",
    "#plt.plot(BS_99_16_24_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "436c7340-0a17-422b-aa7b-65d777d73c98",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}