{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "49deba2e",
   "metadata": {},
   "source": [
    "# GRLP random networks: generate and analyze at scale\n",
    "\n",
    "The hand-built [5-segment network](example_network.ipynb) showed the mechanics of confluences. This capstone does two things that only make sense at scale:\n",
    "\n",
    "1. **Generate** a whole river network with no DEM required — a random **Shreve (1966, 1974)** binary tree, wired up with discharges and widths automatically.\n",
    "2. **Analyze its structure** — Strahler stream orders, Horton ratios, and Hack-type scaling — the descriptors that characterize real drainage networks.\n",
    "\n",
    "This is the quickest way to get a runnable GRLP network for experiments. Reference: [McNab et al. (2025, ESurf)](https://doi.org/10.5194/esurf-13-1059-2025)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8aec77f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import grlp\n",
    "\n",
    "# Seed both RNGs for a reproducible network; drop these to draw a new one.\n",
    "random.seed(7)\n",
    "np.random.seed(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bba9688a",
   "metadata": {},
   "source": [
    "## 1. Generate a Shreve random network\n",
    "\n",
    "`generate_random_network` builds the topology and populates it with node spacing, discharges (accumulating downstream), and valley widths. The main knobs:\n",
    "\n",
    "* `magnitude` — the number of channel heads (the network's \"size\").\n",
    "* `max_length` — the length of the longest source-to-outlet path [m].\n",
    "* `mean_discharge` — sets segment discharges from drainage area.\n",
    "\n",
    "It returns the `Network` and its topology object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65f30029",
   "metadata": {},
   "outputs": [],
   "source": [
    "net, topo = grlp.generate_random_network(\n",
    "    magnitude=8,\n",
    "    max_length=2.0e4,\n",
    "    mean_discharge=10.,\n",
    ")\n",
    "net.set_niter(3)\n",
    "net.get_z_lengths()\n",
    "print('%d segments, %d channel heads'\n",
    "      % (len(net.segments),\n",
    "         len(net.list_of_channel_head_segment_IDs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee7c8662",
   "metadata": {},
   "source": [
    "## 2. Evolve to steady state and plot the long profiles\n",
    "\n",
    "Exactly as for the hand-built network — one call evolves the whole graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8347443",
   "metadata": {},
   "outputs": [],
   "source": [
    "net.evolve_threshold_width_river_network(nt=200, dt=3.15e11)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "for lp in net.segments:\n",
    "    if lp.downstream_segment_IDs:\n",
    "        ds = net.segments[lp.downstream_segment_IDs[0]]\n",
    "        plt.plot([lp.x[-1] / 1000., ds.x[0] / 1000.],\n",
    "                 [lp.z[-1], ds.z[0]], 'k-', lw=1, alpha=0.5)\n",
    "    else:\n",
    "        plt.plot([lp.x[-1] / 1000., lp.x_ghost_downstream / 1000.],\n",
    "                 [lp.z[-1], lp.z_bl], 'k-', lw=1, alpha=0.5)\n",
    "    plt.plot(lp.x / 1000., lp.z, '-', lw=2)\n",
    "plt.xlabel('Downstream distance [km]', fontsize=14)\n",
    "plt.ylabel('Elevation [m]', fontsize=14)\n",
    "plt.title('Shreve random network at steady state')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3607db86",
   "metadata": {},
   "source": [
    "## 3. The network planform\n",
    "\n",
    "`Network.plot()` draws the branching map view — the dendritic pattern of the random tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1faf1c15",
   "metadata": {},
   "outputs": [],
   "source": [
    "net.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16f298e1",
   "metadata": {},
   "source": [
    "## 4. Analyze the network structure\n",
    "\n",
    "Random Shreve networks reproduce the statistical laws of real drainage networks. GRLP computes the classic descriptors.\n",
    "\n",
    "**Strahler stream order** — headwater streams are order 1; two streams of order *n* join to make order *n+1*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bea9a4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "net.compute_strahler_orders()\n",
    "print('Strahler order of each segment:',\n",
    "      [int(o) for o in net.segment_orders])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "820da59a",
   "metadata": {},
   "source": [
    "**Horton ratios** — the near-constant ratios between successive stream orders (Horton 1945; Schumm 1956): how quickly stream *number*, *length*, and *discharge* change from one order to the next. (A *stream* of a given order is a maximal chain of same-order segments, so there are fewer streams than segments at the higher orders.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c40735c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "net.compute_horton_ratios()\n",
    "print('Number of streams per order:',\n",
    "      {int(k): int(v) for k, v in net.order_counts.items()})\n",
    "print('Bifurcation ratio (R_B): %.2f' % net.bifurcation_ratio)\n",
    "print('Length ratio      (R_L): %.2f' % net.length_ratio)\n",
    "print('Discharge ratio   (R_Q): %.2f' % net.discharge_ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c040f90",
   "metadata": {},
   "source": [
    "**Hack-type scaling** — Hack's law relates downstream distance to upstream drainage (here discharge) as a power law. `find_hack_parameters` fits the coefficient `k` and exponent `p`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bfac4ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "hack = net.find_hack_parameters()\n",
    "print('Hack coefficient k = %.4g' % hack['k'])\n",
    "print('Hack exponent    p = %.3f' % hack['p'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58805e22",
   "metadata": {},
   "source": [
    "## Wrap-up\n",
    "\n",
    "That's the three-step path through GRLP: a single [long profile](example_1d.ipynb), a hand-built [network](example_network.ipynb), and a generated-and-analyzed random network. For larger network studies (including the McNab et al., 2025 experiments) and the full API, see [grlp.readthedocs.io](https://grlp.readthedocs.io)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}