graphcast_demo.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-jAYlxeKxvAJ"
      },
      "source": [
        "# GraphCast\n",
        "\n",
        "This colab lets you run several versions of GraphCast.\n",
        "\n",
        "The model weights, normalization statistics, and example inputs are available on [Google Cloud Bucket](https://console.cloud.google.com/storage/browser/dm_graphcast).\n",
        "\n",
        "A Colab runtime with TPU/GPU acceleration will substantially speed up generating predictions and computing the loss/gradients. If you're using a CPU-only runtime, you can switch using the menu \"Runtime \u003e Change runtime type\"."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IIWlNRupdI2i"
      },
      "source": [
        "\u003e \u003cp\u003e\u003csmall\u003e\u003csmall\u003eCopyright 2023 DeepMind Technologies Limited.\u003c/small\u003e\u003c/p\u003e\n",
        "\u003e \u003cp\u003e\u003csmall\u003e\u003csmall\u003eLicensed under the Apache License, Version 2.0 (the \"License\"); you may not use this file except in compliance with the License. You may obtain a copy of the License at \u003ca href=\"http://www.apache.org/licenses/LICENSE-2.0\"\u003ehttp://www.apache.org/licenses/LICENSE-2.0\u003c/a\u003e.\u003c/small\u003e\u003c/small\u003e\u003c/p\u003e\n",
        "\u003e \u003cp\u003e\u003csmall\u003e\u003csmall\u003eUnless required by applicable law or agreed to in writing, software distributed under the License is distributed on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.\u003c/small\u003e\u003c/small\u003e\u003c/p\u003e"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yMbbXFl4msJw"
      },
      "source": [
        "# Installation and Initialization\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "-W4K9skv9vh-"
      },
      "outputs": [],
      "source": [
        "# @title Pip install graphcast and dependencies\n",
        "\n",
        "%pip install --upgrade https://github.com/deepmind/graphcast/archive/master.zip"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "MA5087Vb29z2"
      },
      "outputs": [],
      "source": [
        "# @title Workaround for cartopy crashes\n",
        "\n",
        "# Workaround for cartopy crashes due to the shapely installed by default in\n",
        "# google colab kernel (https://github.com/anitagraser/movingpandas/issues/81):\n",
        "!pip uninstall -y shapely\n",
        "!pip install shapely --no-binary shapely"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "Z_j8ej4Pyg1L"
      },
      "outputs": [],
      "source": [
        "# @title Imports\n",
        "\n",
        "import dataclasses\n",
        "import datetime\n",
        "import functools\n",
        "import math\n",
        "import re\n",
        "from typing import Optional\n",
        "\n",
        "import cartopy.crs as ccrs\n",
        "from google.cloud import storage\n",
        "from graphcast import autoregressive\n",
        "from graphcast import casting\n",
        "from graphcast import checkpoint\n",
        "from graphcast import data_utils\n",
        "from graphcast import graphcast\n",
        "from graphcast import normalization\n",
        "from graphcast import rollout\n",
        "from graphcast import xarray_jax\n",
        "from graphcast import xarray_tree\n",
        "from IPython.display import HTML\n",
        "import ipywidgets as widgets\n",
        "import haiku as hk\n",
        "import jax\n",
        "import matplotlib\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib import animation\n",
        "import numpy as np\n",
        "import xarray\n",
        "\n",
        "\n",
        "def parse_file_parts(file_name):\n",
        "  return dict(part.split(\"-\", 1) for part in file_name.split(\"_\"))\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4wagX1TL_f15"
      },
      "outputs": [],
      "source": [
        "# @title Authenticate with Google Cloud Storage\n",
        "\n",
        "gcs_client = storage.Client.create_anonymous_client()\n",
        "gcs_bucket = gcs_client.get_bucket(\"dm_graphcast\")\n",
        "dir_prefix = \"graphcast/\""
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5JUymx84dI2m"
      },
      "outputs": [],
      "source": [
        "# @title Plotting functions\n",
        "\n",
        "def select(\n",
        "    data: xarray.Dataset,\n",
        "    variable: str,\n",
        "    level: Optional[int] = None,\n",
        "    max_steps: Optional[int] = None\n",
        "    ) -\u003e xarray.Dataset:\n",
        "  data = data[variable]\n",
        "  if \"batch\" in data.dims:\n",
        "    data = data.isel(batch=0)\n",
        "  if max_steps is not None and \"time\" in data.sizes and max_steps \u003c data.sizes[\"time\"]:\n",
        "    data = data.isel(time=range(0, max_steps))\n",
        "  if level is not None and \"level\" in data.coords:\n",
        "    data = data.sel(level=level)\n",
        "  return data\n",
        "\n",
        "def scale(\n",
        "    data: xarray.Dataset,\n",
        "    center: Optional[float] = None,\n",
        "    robust: bool = False,\n",
        "    ) -\u003e tuple[xarray.Dataset, matplotlib.colors.Normalize, str]:\n",
        "  vmin = np.nanpercentile(data, (2 if robust else 0))\n",
        "  vmax = np.nanpercentile(data, (98 if robust else 100))\n",
        "  if center is not None:\n",
        "    diff = max(vmax - center, center - vmin)\n",
        "    vmin = center - diff\n",
        "    vmax = center + diff\n",
        "  return (data, matplotlib.colors.Normalize(vmin, vmax),\n",
        "          (\"RdBu_r\" if center is not None else \"viridis\"))\n",
        "\n",
        "def plot_data(\n",
        "    data: dict[str, xarray.Dataset],\n",
        "    fig_title: str,\n",
        "    plot_size: float = 5,\n",
        "    robust: bool = False,\n",
        "    cols: int = 4\n",
        "    ) -\u003e tuple[xarray.Dataset, matplotlib.colors.Normalize, str]:\n",
        "\n",
        "  first_data = next(iter(data.values()))[0]\n",
        "  max_steps = first_data.sizes.get(\"time\", 1)\n",
        "  assert all(max_steps == d.sizes.get(\"time\", 1) for d, _, _ in data.values())\n",
        "\n",
        "  cols = min(cols, len(data))\n",
        "  rows = math.ceil(len(data) / cols)\n",
        "  figure = plt.figure(figsize=(plot_size * 2 * cols,\n",
        "                               plot_size * rows))\n",
        "  figure.suptitle(fig_title, fontsize=16)\n",
        "  figure.subplots_adjust(wspace=0, hspace=0)\n",
        "  figure.tight_layout()\n",
        "\n",
        "  images = []\n",
        "  for i, (title, (plot_data, norm, cmap)) in enumerate(data.items()):\n",
        "    ax = figure.add_subplot(rows, cols, i+1)\n",
        "    ax.set_xticks([])\n",
        "    ax.set_yticks([])\n",
        "    ax.set_title(title)\n",
        "    im = ax.imshow(\n",
        "        plot_data.isel(time=0, missing_dims=\"ignore\"), norm=norm,\n",
        "        origin=\"lower\", cmap=cmap)\n",
        "    plt.colorbar(\n",
        "        mappable=im,\n",
        "        ax=ax,\n",
        "        orientation=\"vertical\",\n",
        "        pad=0.02,\n",
        "        aspect=16,\n",
        "        shrink=0.75,\n",
        "        cmap=cmap,\n",
        "        extend=(\"both\" if robust else \"neither\"))\n",
        "    images.append(im)\n",
        "\n",
        "  def update(frame):\n",
        "    if \"time\" in first_data.dims:\n",
        "      td = datetime.timedelta(microseconds=first_data[\"time\"][frame].item() / 1000)\n",
        "      figure.suptitle(f\"{fig_title}, {td}\", fontsize=16)\n",
        "    else:\n",
        "      figure.suptitle(fig_title, fontsize=16)\n",
        "    for im, (plot_data, norm, cmap) in zip(images, data.values()):\n",
        "      im.set_data(plot_data.isel(time=frame, missing_dims=\"ignore\"))\n",
        "\n",
        "  ani = animation.FuncAnimation(\n",
        "      fig=figure, func=update, frames=max_steps, interval=250)\n",
        "  plt.close(figure.number)\n",
        "  return HTML(ani.to_jshtml())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WEtSV8HEkHtf"
      },
      "source": [
        "# Load the Data and initialize the model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G50ORsY_dI2n"
      },
      "source": [
        "## Load the model params\n",
        "\n",
        "Choose one of the two ways of getting model params:\n",
        "- **random**: You'll get random predictions, but you can change the model architecture, which may run faster or fit on your device.\n",
        "- **checkpoint**: You'll get sensible predictions, but are limited to the model architecture that it was trained with, which may not fit on your device. In particular generating gradients uses a lot of memory, so you'll need at least 25GB of ram (TPUv4 or A100).\n",
        "\n",
        "Checkpoints vary across a few axes:\n",
        "- The mesh size specifies the internal graph representation of the earth. Smaller meshes will run faster but will have worse outputs. The mesh size does not affect the number of parameters of the model.\n",
        "- The resolution and number of pressure levels must match the data. Lower resolution and fewer levels will run a bit faster. Data resolution only affects the encoder/decoder.\n",
        "- All our models predict precipitation. However, ERA5 includes precipitation, while HRES does not. Our models marked as \"ERA5\" take precipitation as input and expect ERA5 data as input, while model marked \"ERA5-HRES\" do not take precipitation as input and are specifically trained to take HRES-fc0 as input (see the data section below).\n",
        "\n",
        "We provide three pre-trained models.\n",
        "1. `GraphCast`, the high-resolution model used in the GraphCast paper (0.25 degree resolution, 37 pressure levels), trained on ERA5 data from 1979 to 2017,\n",
        "\n",
        "2. `GraphCast_small`, a smaller, low-resolution version of GraphCast (1 degree resolution, 13 pressure levels, and a smaller mesh), trained on ERA5 data from 1979 to 2015, useful to run a model with lower memory and compute constraints,\n",
        "\n",
        "3. `GraphCast_operational`, a high-resolution model (0.25 degree resolution, 13 pressure levels) pre-trained on ERA5 data from 1979 to 2017 and fine-tuned on HRES data from 2016 to 2021. This model can be initialized from HRES data (does not require precipitation inputs).\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KGaJ6V9MdI2n"
      },
      "outputs": [],
      "source": [
        "# @title Choose the model\n",
        "\n",
        "params_file_options = [\n",
        "    name for blob in gcs_bucket.list_blobs(prefix=dir_prefix+\"params/\")\n",
        "    if (name := blob.name.removeprefix(dir_prefix+\"params/\"))]  # Drop empty string.\n",
        "\n",
        "random_mesh_size = widgets.IntSlider(\n",
        "    value=4, min=4, max=6, description=\"Mesh size:\")\n",
        "random_gnn_msg_steps = widgets.IntSlider(\n",
        "    value=4, min=1, max=32, description=\"GNN message steps:\")\n",
        "random_latent_size = widgets.Dropdown(\n",
        "    options=[int(2**i) for i in range(4, 10)], value=32,description=\"Latent size:\")\n",
        "random_levels = widgets.Dropdown(\n",
        "    options=[13, 37], value=13, description=\"Pressure levels:\")\n",
        "\n",
        "\n",
        "params_file = widgets.Dropdown(\n",
        "    options=params_file_options,\n",
        "    description=\"Params file:\",\n",
        "    layout={\"width\": \"max-content\"})\n",
        "\n",
        "source_tab = widgets.Tab([\n",
        "    widgets.VBox([\n",
        "        random_mesh_size,\n",
        "        random_gnn_msg_steps,\n",
        "        random_latent_size,\n",
        "        random_levels,\n",
        "    ]),\n",
        "    params_file,\n",
        "])\n",
        "source_tab.set_title(0, \"Random\")\n",
        "source_tab.set_title(1, \"Checkpoint\")\n",
        "widgets.VBox([\n",
        "    source_tab,\n",
        "    widgets.Label(value=\"Run the next cell to load the model. Rerunning this cell clears your selection.\")\n",
        "])\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lYQgrPgPdI2n"
      },
      "outputs": [],
      "source": [
        "# @title Load the model\n",
        "\n",
        "source = source_tab.get_title(source_tab.selected_index)\n",
        "\n",
        "if source == \"Random\":\n",
        "  params = None  # Filled in below\n",
        "  state = {}\n",
        "  model_config = graphcast.ModelConfig(\n",
        "      resolution=0,\n",
        "      mesh_size=random_mesh_size.value,\n",
        "      latent_size=random_latent_size.value,\n",
        "      gnn_msg_steps=random_gnn_msg_steps.value,\n",
        "      hidden_layers=1,\n",
        "      radius_query_fraction_edge_length=0.6)\n",
        "  task_config = graphcast.TaskConfig(\n",
        "      input_variables=graphcast.TASK.input_variables,\n",
        "      target_variables=graphcast.TASK.target_variables,\n",
        "      forcing_variables=graphcast.TASK.forcing_variables,\n",
        "      pressure_levels=graphcast.PRESSURE_LEVELS[random_levels.value],\n",
        "      input_duration=graphcast.TASK.input_duration,\n",
        "  )\n",
        "else:\n",
        "  assert source == \"Checkpoint\"\n",
        "  with gcs_bucket.blob(f\"{dir_prefix}params/{params_file.value}\").open(\"rb\") as f:\n",
        "    ckpt = checkpoint.load(f, graphcast.CheckPoint)\n",
        "  params = ckpt.params\n",
        "  state = {}\n",
        "\n",
        "  model_config = ckpt.model_config\n",
        "  task_config = ckpt.task_config\n",
        "  print(\"Model description:\\n\", ckpt.description, \"\\n\")\n",
        "  print(\"Model license:\\n\", ckpt.license, \"\\n\")\n",
        "\n",
        "model_config"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rQWk0RRuCjDN"
      },
      "source": [
        "## Load the example data\n",
        "\n",
        "Several example datasets are available, varying across a few axes:\n",
        "- **Source**: fake, era5, hres\n",
        "- **Resolution**: 0.25deg, 1deg, 6deg\n",
        "- **Levels**: 13, 37\n",
        "- **Steps**: How many timesteps are included\n",
        "\n",
        "Not all combinations are available.\n",
        "- Higher resolution is only available for fewer steps due to the memory requirements of loading them.\n",
        "- HRES is only available in 0.25 deg, with 13 pressure levels.\n",
        "\n",
        "The data resolution must match the model that is loaded.\n",
        "\n",
        "Some transformations were done from the base datasets:\n",
        "- We accumulated precipitation over 6 hours instead of the default 1 hour.\n",
        "- For HRES data, each time step corresponds to the HRES forecast at leadtime 0, essentially providing an \"initialisation\" from HRES. See HRES-fc0 in the GraphCast paper for further description. Note that a 6h accumulation of precipitation is not available from HRES, so our model taking HRES inputs does not depend on precipitation. However, because our models predict precipitation, we include the ERA5 precipitation in the example data so it can serve as an illustrative example of ground truth.\n",
        "- We include ERA5 `toa_incident_solar_radiation` in the data. Our model uses the radiation at -6h, 0h and +6h as a forcing term for each 1-step prediction. If the radiation is missing from the data (e.g. in an operational setting), it will be computed using a custom implementation that produces values similar to those in ERA5."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-DJzie5me2-H"
      },
      "outputs": [],
      "source": [
        "# @title Get and filter the list of available example datasets\n",
        "\n",
        "dataset_file_options = [\n",
        "    name for blob in gcs_bucket.list_blobs(prefix=dir_prefix+\"dataset/\")\n",
        "    if (name := blob.name.removeprefix(dir_prefix+\"dataset/\"))]  # Drop empty string.\n",
        "\n",
        "def data_valid_for_model(\n",
        "    file_name: str, model_config: graphcast.ModelConfig, task_config: graphcast.TaskConfig):\n",
        "  file_parts = parse_file_parts(file_name.removesuffix(\".nc\"))\n",
        "  return (\n",
        "      model_config.resolution in (0, float(file_parts[\"res\"])) and\n",
        "      len(task_config.pressure_levels) == int(file_parts[\"levels\"]) and\n",
        "      (\n",
        "          (\"total_precipitation_6hr\" in task_config.input_variables and\n",
        "           file_parts[\"source\"] in (\"era5\", \"fake\")) or\n",
        "          (\"total_precipitation_6hr\" not in task_config.input_variables and\n",
        "           file_parts[\"source\"] in (\"hres\", \"fake\"))\n",
        "      )\n",
        "  )\n",
        "\n",
        "\n",
        "dataset_file = widgets.Dropdown(\n",
        "    options=[\n",
        "        (\", \".join([f\"{k}: {v}\" for k, v in parse_file_parts(option.removesuffix(\".nc\")).items()]), option)\n",
        "        for option in dataset_file_options\n",
        "        if data_valid_for_model(option, model_config, task_config)\n",
        "    ],\n",
        "    description=\"Dataset file:\",\n",
        "    layout={\"width\": \"max-content\"})\n",
        "widgets.VBox([\n",
        "    dataset_file,\n",
        "    widgets.Label(value=\"Run the next cell to load the dataset. Rerunning this cell clears your selection and refilters the datasets that match your model.\")\n",
        "])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Yz-ekISoJxeZ"
      },
      "outputs": [],
      "source": [
        "# @title Load weather data\n",
        "\n",
        "if not data_valid_for_model(dataset_file.value, model_config, task_config):\n",
        "  raise ValueError(\n",
        "      \"Invalid dataset file, rerun the cell above and choose a valid dataset file.\")\n",
        "\n",
        "with gcs_bucket.blob(f\"{dir_prefix}dataset/{dataset_file.value}\").open(\"rb\") as f:\n",
        "  example_batch = xarray.load_dataset(f).compute()\n",
        "\n",
        "assert example_batch.dims[\"time\"] \u003e= 3  # 2 for input, \u003e=1 for targets\n",
        "\n",
        "print(\", \".join([f\"{k}: {v}\" for k, v in parse_file_parts(dataset_file.value.removesuffix(\".nc\")).items()]))\n",
        "\n",
        "example_batch"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "lXjFvdE6qStr"
      },
      "outputs": [],
      "source": [
        "# @title Choose data to plot\n",
        "\n",
        "plot_example_variable = widgets.Dropdown(\n",
        "    options=example_batch.data_vars.keys(),\n",
        "    value=\"2m_temperature\",\n",
        "    description=\"Variable\")\n",
        "plot_example_level = widgets.Dropdown(\n",
        "    options=example_batch.coords[\"level\"].values,\n",
        "    value=500,\n",
        "    description=\"Level\")\n",
        "plot_example_robust = widgets.Checkbox(value=True, description=\"Robust\")\n",
        "plot_example_max_steps = widgets.IntSlider(\n",
        "    min=1, max=example_batch.dims[\"time\"], value=example_batch.dims[\"time\"],\n",
        "    description=\"Max steps\")\n",
        "\n",
        "widgets.VBox([\n",
        "    plot_example_variable,\n",
        "    plot_example_level,\n",
        "    plot_example_robust,\n",
        "    plot_example_max_steps,\n",
        "    widgets.Label(value=\"Run the next cell to plot the data. Rerunning this cell clears your selection.\")\n",
        "])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "kIK-EgMdkHtk"
      },
      "outputs": [],
      "source": [
        "# @title Plot example data\n",
        "\n",
        "plot_size = 7\n",
        "\n",
        "data = {\n",
        "    \" \": scale(select(example_batch, plot_example_variable.value, plot_example_level.value, plot_example_max_steps.value),\n",
        "              robust=plot_example_robust.value),\n",
        "}\n",
        "fig_title = plot_example_variable.value\n",
        "if \"level\" in example_batch[plot_example_variable.value].coords:\n",
        "  fig_title += f\" at {plot_example_level.value} hPa\"\n",
        "\n",
        "plot_data(data, fig_title, plot_size, plot_example_robust.value)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "tPVy1GHokHtk"
      },
      "outputs": [],
      "source": [
        "# @title Choose training and eval data to extract\n",
        "train_steps = widgets.IntSlider(\n",
        "    value=1, min=1, max=example_batch.sizes[\"time\"]-2, description=\"Train steps\")\n",
        "eval_steps = widgets.IntSlider(\n",
        "    value=example_batch.sizes[\"time\"]-2, min=1, max=example_batch.sizes[\"time\"]-2, description=\"Eval steps\")\n",
        "\n",
        "widgets.VBox([\n",
        "    train_steps,\n",
        "    eval_steps,\n",
        "    widgets.Label(value=\"Run the next cell to extract the data. Rerunning this cell clears your selection.\")\n",
        "])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "Ogp4vTBvsgSt"
      },
      "outputs": [],
      "source": [
        "# @title Extract training and eval data\n",
        "\n",
        "train_inputs, train_targets, train_forcings = data_utils.extract_inputs_targets_forcings(\n",
        "    example_batch, target_lead_times=slice(\"6h\", f\"{train_steps.value*6}h\"),\n",
        "    **dataclasses.asdict(task_config))\n",
        "\n",
        "eval_inputs, eval_targets, eval_forcings = data_utils.extract_inputs_targets_forcings(\n",
        "    example_batch, target_lead_times=slice(\"6h\", f\"{eval_steps.value*6}h\"),\n",
        "    **dataclasses.asdict(task_config))\n",
        "\n",
        "print(\"All Examples:  \", example_batch.dims.mapping)\n",
        "print(\"Train Inputs:  \", train_inputs.dims.mapping)\n",
        "print(\"Train Targets: \", train_targets.dims.mapping)\n",
        "print(\"Train Forcings:\", train_forcings.dims.mapping)\n",
        "print(\"Eval Inputs:   \", eval_inputs.dims.mapping)\n",
        "print(\"Eval Targets:  \", eval_targets.dims.mapping)\n",
        "print(\"Eval Forcings: \", eval_forcings.dims.mapping)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Q--ZRhpTdI2o"
      },
      "outputs": [],
      "source": [
        "# @title Load normalization data\n",
        "\n",
        "with gcs_bucket.blob(dir_prefix+\"stats/diffs_stddev_by_level.nc\").open(\"rb\") as f:\n",
        "  diffs_stddev_by_level = xarray.load_dataset(f).compute()\n",
        "with gcs_bucket.blob(dir_prefix+\"stats/mean_by_level.nc\").open(\"rb\") as f:\n",
        "  mean_by_level = xarray.load_dataset(f).compute()\n",
        "with gcs_bucket.blob(dir_prefix+\"stats/stddev_by_level.nc\").open(\"rb\") as f:\n",
        "  stddev_by_level = xarray.load_dataset(f).compute()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "ke2zQyuT_sMA"
      },
      "outputs": [],
      "source": [
        "# @title Build jitted functions, and possibly initialize random weights\n",
        "\n",
        "def construct_wrapped_graphcast(\n",
        "    model_config: graphcast.ModelConfig,\n",
        "    task_config: graphcast.TaskConfig):\n",
        "  \"\"\"Constructs and wraps the GraphCast Predictor.\"\"\"\n",
        "  # Deeper one-step predictor.\n",
        "  predictor = graphcast.GraphCast(model_config, task_config)\n",
        "\n",
        "  # Modify inputs/outputs to `graphcast.GraphCast` to handle conversion to\n",
        "  # from/to float32 to/from BFloat16.\n",
        "  predictor = casting.Bfloat16Cast(predictor)\n",
        "\n",
        "  # Modify inputs/outputs to `casting.Bfloat16Cast` so the casting to/from\n",
        "  # BFloat16 happens after applying normalization to the inputs/targets.\n",
        "  predictor = normalization.InputsAndResiduals(\n",
        "      predictor,\n",
        "      diffs_stddev_by_level=diffs_stddev_by_level,\n",
        "      mean_by_level=mean_by_level,\n",
        "      stddev_by_level=stddev_by_level)\n",
        "\n",
        "  # Wraps everything so the one-step model can produce trajectories.\n",
        "  predictor = autoregressive.Predictor(predictor, gradient_checkpointing=True)\n",
        "  return predictor\n",
        "\n",
        "\n",
        "@hk.transform_with_state\n",
        "def run_forward(model_config, task_config, inputs, targets_template, forcings):\n",
        "  predictor = construct_wrapped_graphcast(model_config, task_config)\n",
        "  return predictor(inputs, targets_template=targets_template, forcings=forcings)\n",
        "\n",
        "\n",
        "@hk.transform_with_state\n",
        "def loss_fn(model_config, task_config, inputs, targets, forcings):\n",
        "  predictor = construct_wrapped_graphcast(model_config, task_config)\n",
        "  loss, diagnostics = predictor.loss(inputs, targets, forcings)\n",
        "  return xarray_tree.map_structure(\n",
        "      lambda x: xarray_jax.unwrap_data(x.mean(), require_jax=True),\n",
        "      (loss, diagnostics))\n",
        "\n",
        "def grads_fn(params, state, model_config, task_config, inputs, targets, forcings):\n",
        "  def _aux(params, state, i, t, f):\n",
        "    (loss, diagnostics), next_state = loss_fn.apply(\n",
        "        params, state, jax.random.PRNGKey(0), model_config, task_config,\n",
        "        i, t, f)\n",
        "    return loss, (diagnostics, next_state)\n",
        "  (loss, (diagnostics, next_state)), grads = jax.value_and_grad(\n",
        "      _aux, has_aux=True)(params, state, inputs, targets, forcings)\n",
        "  return loss, diagnostics, next_state, grads\n",
        "\n",
        "# Jax doesn't seem to like passing configs as args through the jit. Passing it\n",
        "# in via partial (instead of capture by closure) forces jax to invalidate the\n",
        "# jit cache if you change configs.\n",
        "def with_configs(fn):\n",
        "  return functools.partial(\n",
        "      fn, model_config=model_config, task_config=task_config)\n",
        "\n",
        "# Always pass params and state, so the usage below are simpler\n",
        "def with_params(fn):\n",
        "  return functools.partial(fn, params=params, state=state)\n",
        "\n",
        "# Our models aren't stateful, so the state is always empty, so just return the\n",
        "# predictions. This is requiredy by our rollout code, and generally simpler.\n",
        "def drop_state(fn):\n",
        "  return lambda **kw: fn(**kw)[0]\n",
        "\n",
        "init_jitted = jax.jit(with_configs(run_forward.init))\n",
        "\n",
        "if params is None:\n",
        "  params, state = init_jitted(\n",
        "      rng=jax.random.PRNGKey(0),\n",
        "      inputs=train_inputs,\n",
        "      targets_template=train_targets,\n",
        "      forcings=train_forcings)\n",
        "\n",
        "loss_fn_jitted = drop_state(with_params(jax.jit(with_configs(loss_fn.apply))))\n",
        "grads_fn_jitted = with_params(jax.jit(with_configs(grads_fn)))\n",
        "run_forward_jitted = drop_state(with_params(jax.jit(with_configs(\n",
        "    run_forward.apply))))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VBNutliiCyqA"
      },
      "source": [
        "# Run the model\n",
        "\n",
        "Note that the cell below may take a while (possibly minutes) to run the first time you execute them, because this will include the time it takes for the code to compile. The second time running will be significantly faster.\n",
        "\n",
        "This use the python loop to iterate over prediction steps, where the 1-step prediction is jitted. This has lower memory requirements than the training steps below, and should enable making prediction with the small GraphCast model on 1 deg resolution data for 4 steps."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "7obeY9i9oTtD"
      },
      "outputs": [],
      "source": [
        "# @title Autoregressive rollout (loop in python)\n",
        "\n",
        "assert model_config.resolution in (0, 360. / eval_inputs.sizes[\"lon\"]), (\n",
        "  \"Model resolution doesn't match the data resolution. You likely want to \"\n",
        "  \"re-filter the dataset list, and download the correct data.\")\n",
        "\n",
        "print(\"Inputs:  \", eval_inputs.dims.mapping)\n",
        "print(\"Targets: \", eval_targets.dims.mapping)\n",
        "print(\"Forcings:\", eval_forcings.dims.mapping)\n",
        "\n",
        "predictions = rollout.chunked_prediction(\n",
        "    run_forward_jitted,\n",
        "    rng=jax.random.PRNGKey(0),\n",
        "    inputs=eval_inputs,\n",
        "    targets_template=eval_targets * np.nan,\n",
        "    forcings=eval_forcings)\n",
        "predictions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "ft298eZskHtn"
      },
      "outputs": [],
      "source": [
        "# @title Choose predictions to plot\n",
        "\n",
        "plot_pred_variable = widgets.Dropdown(\n",
        "    options=predictions.data_vars.keys(),\n",
        "    value=\"2m_temperature\",\n",
        "    description=\"Variable\")\n",
        "plot_pred_level = widgets.Dropdown(\n",
        "    options=predictions.coords[\"level\"].values,\n",
        "    value=500,\n",
        "    description=\"Level\")\n",
        "plot_pred_robust = widgets.Checkbox(value=True, description=\"Robust\")\n",
        "plot_pred_max_steps = widgets.IntSlider(\n",
        "    min=1,\n",
        "    max=predictions.dims[\"time\"],\n",
        "    value=predictions.dims[\"time\"],\n",
        "    description=\"Max steps\")\n",
        "\n",
        "widgets.VBox([\n",
        "    plot_pred_variable,\n",
        "    plot_pred_level,\n",
        "    plot_pred_robust,\n",
        "    plot_pred_max_steps,\n",
        "    widgets.Label(value=\"Run the next cell to plot the predictions. Rerunning this cell clears your selection.\")\n",
        "])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "_tTdx6fmmj1I"
      },
      "outputs": [],
      "source": [
        "# @title Plot predictions\n",
        "\n",
        "plot_size = 5\n",
        "plot_max_steps = min(predictions.dims[\"time\"], plot_pred_max_steps.value)\n",
        "\n",
        "data = {\n",
        "    \"Targets\": scale(select(eval_targets, plot_pred_variable.value, plot_pred_level.value, plot_max_steps), robust=plot_pred_robust.value),\n",
        "    \"Predictions\": scale(select(predictions, plot_pred_variable.value, plot_pred_level.value, plot_max_steps), robust=plot_pred_robust.value),\n",
        "    \"Diff\": scale((select(eval_targets, plot_pred_variable.value, plot_pred_level.value, plot_max_steps) -\n",
        "                        select(predictions, plot_pred_variable.value, plot_pred_level.value, plot_max_steps)),\n",
        "                       robust=plot_pred_robust.value, center=0),\n",
        "}\n",
        "fig_title = plot_pred_variable.value\n",
        "if \"level\" in predictions[plot_pred_variable.value].coords:\n",
        "  fig_title += f\" at {plot_pred_level.value} hPa\"\n",
        "\n",
        "plot_data(data, fig_title, plot_size, plot_pred_robust.value)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Pa78b64bLYe1"
      },
      "source": [
        "# Train the model\n",
        "\n",
        "The following operations require a large amount of memory and, depending on the accelerator being used, will only fit the very small \"random\" model on low resolution data. It uses the number of training steps selected above.\n",
        "\n",
        "The first time executing the cell takes more time, as it include the time to jit the function."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "Nv-u3dAP7IRZ"
      },
      "outputs": [],
      "source": [
        "# @title Loss computation (autoregressive loss over multiple steps)\n",
        "loss, diagnostics = loss_fn_jitted(\n",
        "    rng=jax.random.PRNGKey(0),\n",
        "    inputs=train_inputs,\n",
        "    targets=train_targets,\n",
        "    forcings=train_forcings)\n",
        "print(\"Loss:\", float(loss))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "mBNFq1IGZNLz"
      },
      "outputs": [],
      "source": [
        "# @title Gradient computation (backprop through time)\n",
        "loss, diagnostics, next_state, grads = grads_fn_jitted(\n",
        "    inputs=train_inputs,\n",
        "    targets=train_targets,\n",
        "    forcings=train_forcings)\n",
        "mean_grad = np.mean(jax.tree_util.tree_flatten(jax.tree_util.tree_map(lambda x: np.abs(x).mean(), grads))[0])\n",
        "print(f\"Loss: {loss:.4f}, Mean |grad|: {mean_grad:.6f}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "J4FJFKWD8Loz"
      },
      "outputs": [],
      "source": [
        "# @title Autoregressive rollout (keep the loop in JAX)\n",
        "print(\"Inputs:  \", train_inputs.dims.mapping)\n",
        "print(\"Targets: \", train_targets.dims.mapping)\n",
        "print(\"Forcings:\", train_forcings.dims.mapping)\n",
        "\n",
        "predictions = run_forward_jitted(\n",
        "    rng=jax.random.PRNGKey(0),\n",
        "    inputs=train_inputs,\n",
        "    targets_template=train_targets * np.nan,\n",
        "    forcings=train_forcings)\n",
        "predictions"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "name": "GraphCast",
      "private_outputs": true,
      "provenance": [],
      "toc_visible": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}