{ "cells": [ { "cell_type": "markdown", "id": "79e7c63b", "metadata": {}, "source": [ "# Sweet Butter 🧈\n", "\n", "Sweet Butter is a shortest paths problem by Greg Galperin from the [USACO training pages](https://train.usaco.org/):\n", "\n", "> *Farmer John has discovered the secret to making the sweetest butter in all of Wisconsin: sugar. By placing a sugar cube out in the pastures, he knows the cows will lick it and thus will produce super-sweet butter which can be marketed at better prices. Of course, he spends the extra money on luxuries for the cows.*\n", ">\n", "> *FJ is a sly farmer. Like Pavlov of old, he knows he can train the cows to go to a certain pasture when they hear a bell. He intends to put the sugar there and then ring the bell in the middle of the afternoon so that the evening's milking produces perfect milk.*\n", ">\n", "> *FJ knows each cow spends her time in a given pasture (not necessarily alone). Given the pasture location of the cows and a description of the paths that connect the pastures, find the pasture in which to place the sugar cube so that the total distance walked by the cows when FJ rings the bell is minimized. FJ knows the fields are connected well enough that some solution is always possible.*\n", "\n", "Rephrased:\n", "\n", "> Given a connected undirected weighted graph with some active nodes, find the node that minimizes the total distance taken from every active node to it.\n", "\n", "The key concepts for this problem are:\n", "\n", "- **Node**: pasture.\n", "- **Weighted Edge**: length of a direct connection between two pastures. \n", "- **Undirected**: edges can be traversed in both directions.\n", "- **Shortest Path**: minimal sum of edge weights to get from pasture A to pasture B. (Other problems may want this as a number of edges or sequence of nodes taken.)\n", "- **Adjacency List**: a data structure to represent a [graph](https://en.wikipedia.org/wiki/Graph_%28discrete_mathematics%29). Each node *x* in the graph is assigned an [adjacency list](https://en.wikipedia.org/wiki/Adjacency_list) that consists of nodes to where there is an edge from *x*.\n", "- **Dijkstra's Algorithm**: an algorithm for finding the shortest paths from a starting node to all other nodes. \n", " - **Greedy**: a type of algorithm (e.g. Dijkstra's) that constructs its solution by always making a choice that looks best at the moment.\n", " - **Heap**: a data structure that gives easy access and removal of the minimum or maximum element. Used to implement a priority queue.\n", "\n", "Related concepts are:\n", "\n", "- **Adjacency Matrix**: a data structure to represent a graph.\n", "- **Floyd-Warshall Algorithm**: a shortest paths algorithm.\n", "\n", "## Imports, Types, and Utility Functions" ] }, { "cell_type": "code", "execution_count": 1, "id": "4c465ed3", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import networkx as nx\n", "import heapq\n", "from typing import *\n", "from collections import defaultdict\n", "\n", "Cow = Dict[int, int] # A cow and its pasture\n", "Edge = Tuple[int, int] # A pasture and connection length\n", "Adj = DefaultDict[int, List[Edge]] # Adjacency list, e.g. {1: [(2, 1), (3, 5)], 2: [(1, 1), (3, 7), (4, 3)], ...}\n", "Data = Tuple[Cow, Adj] # Data read from input\n", "Dist = Dict[int, int] # Shortest paths" ] }, { "cell_type": "code", "execution_count": 2, "id": "a5e9d3c3", "metadata": {}, "outputs": [], "source": [ "def read_input(file: str='butter.in') -> Data:\n", " \"\"\"Store cow locations and adjacency list from input data.\"\"\"\n", " cow, adj = {}, defaultdict(list)\n", " with open(file) as f:\n", " C, V, E = map(int, f.readline().split()) # cows, vertices (nodes), edges\n", " for i in range(1, C + 1):\n", " cow[i] = int(f.readline())\n", " for _ in range(E):\n", " u, v, w = map(int, f.readline().split()) # undirected edge connecting u and v with weight w\n", " adj[u].append((v, w))\n", " adj[v].append((u, w))\n", " return cow, adj" ] }, { "cell_type": "code", "execution_count": 3, "id": "31be6d31", "metadata": {}, "outputs": [], "source": [ "def draw_graph(data: Data, is_large: bool=False):\n", " \"\"\"Visualize an adjacency list and show cows in each node.\"\"\"\n", " cow, adj = data\n", " G = nx.Graph()\n", " for u in adj:\n", " G.add_node(u, cow=0)\n", " for edge in adj[u]:\n", " v, w = edge[0], edge[1]\n", " G.add_edge(u, v, weight=w)\n", " for i in cow:\n", " G.nodes[cow[i]]['cow'] += 1\n", " \n", " if is_large: \n", " plt.figure(1, figsize=(75, 75))\n", " \n", " # draw nodes\n", " layout = nx.circular_layout(G) # positions for all nodes\n", " nx.draw(G, pos=layout, with_labels=True, node_color=\"lightgreen\")\n", " \n", " # draw cow attribute label above nodes\n", " pos_attrs = {} # positions for cow labels\n", " for u, xy in layout.items(): \n", " pos_attrs[u] = (xy[0], xy[1] + 0.12) \n", " cow_labels = {} \n", " for u, val in nx.get_node_attributes(G, 'cow').items():\n", " cow_labels[u] = '' if is_large else \"🐮 x\" + str(val) # no cow labels on large graphs; looks nicer\n", " nx.draw_networkx_labels(G, pos=pos_attrs, labels=cow_labels)\n", " \n", " # draw edge weights\n", " weight_labels = nx.get_edge_attributes(G, \"weight\")\n", " nx.draw_networkx_edge_labels(G, pos=layout, edge_labels=weight_labels) \n", " plt.show()" ] }, { "cell_type": "markdown", "id": "dc0f11c7", "metadata": {}, "source": [ "## Input\n", "\n", "The first line contains the number of cows, pastures, and edges respectively. The next cow number of lines is each cow and its pasture. Followed by lines containing connected pastures and the distance between them. \n", "\n", "> 3 4 5 \n", "> 2 \n", "> 3 \n", "> 4 \n", "> 1 2 1 \n", "> 1 3 5 \n", "> 2 3 7 \n", "> 2 4 3 \n", "> 3 4 5\n", "\n", "Cows and their pasture are stored in a dictionary and an adjacency list stores the graph of pastures and connections." ] }, { "cell_type": "code", "execution_count": 4, "id": "6ea6996a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "({1: 2, 2: 3, 3: 4},\n", " defaultdict(list,\n", " {1: [(2, 1), (3, 5)],\n", " 2: [(1, 1), (3, 7), (4, 3)],\n", " 3: [(1, 5), (2, 7), (4, 5)],\n", " 4: [(2, 3), (3, 5)]}))" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "read_input()" ] }, { "cell_type": "markdown", "id": "b20fe462", "metadata": {}, "source": [ "Here is the sample input visualized." ] }, { "cell_type": "code", "execution_count": 5, "id": "08de8c2b", "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgpUlEQVR4nO3df3BU9b3/8dcmS34AIUq4BgYFLj9qGzFpNIhoE5sbjJGGUr5CpePXKMhNQAp1RqIgVUavhIvVzO0UIVoGxNLKF7C3w28wFkqi00JkAkn5YWGa4YollaBJIGzM7p7vH14CIQnJJrt7zu55PmaYgezuyWcZOO88d8856zAMwxAAADYRYfYCAAAIJgYfAMBWGHwAAFth8AEAbIXBBwCwFQYfAMBWGHwAAFth8AEAbIXBBwCwFQYfAMBWGHxAL1RXV2vu3LkqLy/v9D5nz57VnDlztH379iCuDEBnGHxAD1VVVWn27Nm69957VVhYqJKSEuXk5Cg1NVWpqamaMGGCVq5cqalTpyolJUUbNmzQ+vXrff4++/btU2ZmpuLj4zVixAj/PxHAZhxcpBromWXLlikpKUnDhw/Xjh07VFZWptWrV2vUqFGSpPr6ei1atEgxMTEaP368srOzNWPGDO3du9en73Pw4EGdPHlSly9fVlFRkWpqagLwbAD7oPiAHsrOzlZxcbFGjhypW265RXv27GlTZAMGDNCbb76pjIwMZWZmavHixcrNzW23nQsXLujWW2/Vtm3bJEkXL17U6NGj9e6770qS7rnnHj3++OMaOXJkUJ4XEO6cZi8ACFXjxo3TkiVLlJ6erri4OMXFxamkpEQHDhyQJC1dulRut1s1NTVaunSpZsyYoQULFrTbzsCBA7V27Vrl5eXp6NGjWrJkib773e8qLy8v2E8JsAUGH9ALOTk5ys7O1vPPP6/jx49r7NixrbdFRUVp6NChOnDggCoqKhQVFdXpdrKzszV9+nRlZWWprq5OVVVVwVg+YEu81An0UkREhHJycrR161YtX75cX3/9tSRp/vz5qqysVGZm5g2H3hX5+fmqrq7WzJkzlZCQEOhlA7bF4AP8YMyYMRo2bJg2b96sjz/+WJI0b9483X777RozZkyXj/d4PCooKFBeXp5Wr16tU6dOBXrJgG0x+AA/GDJkiM6cOaOmpib169dPlZWVys3N1bp165SWltbl44uKiiRJa9eu1cKFC5WXlyePxyNJ8nq9crlcamlpkWEYcrlcrVUJwHeczgD4SXl5uaZMmaIpU6bI4XCovLxcjz32mF566aUbPu6TTz7RxIkTdejQIY0ePVoej0cZGRmaNGmSlixZov379yszM7PNYx544AHt378/gM8GCF8MPsCPVqxYoU2bNskwDKWkpGjVqlWKjY01e1kArsHgAwDYCqczAEHU5G3SseZjOu85r2ajWdGOaA2KHKSk6CT1jehr9vIAW6D4gCA45z6nCleFalpqJEkeeVpvi1SkJGlEnxFKi0nTYOdgM5YI2AaDDwiwo66jKrtcJrfcXd7XKafSY9OVHJMchJUB9sTpDEAA+TL0JMktt8oul+mo62iAVwbYF+/xAQFyzn2u3dBzN7u1eeFmffqnT9X0VZMG/esg/eDnP1DSg0lX7/O/wy/RmahEZ6IZSwfCGsUHBEiFq6Jd6XncHt009Cb9dPtPtbxmuSa9MEnrn1qvujN1be7nlluHXIeCuVzANhh8QAA0eZtaD2S5VnS/aD286GElDEtQRESE7njoDg0cNlCfVX7W7r41LTVq8jYFYbWAvTD4gAA41nysW/dr/Gejvjj9hQZ/u+MjObu7HQDdx+ADAuC853ybUxY64mnx6DcFv9G4GeOU+K327+V55FGdt66DRwLoDQYfEADNRvMNb/d6vdowZ4Mi+0Rq2mvTOt+O98bbAeA7Bh8QANGO6E5vMwxDG+dvVOMXjZq5fqYi+0R2vp2IzrcDoGcYfEAADIoc1HpFluttfnazaj+t1b//7t8VFdv5B9RGKlIJEXwgLeBvXLkFCIAmb5PW1q9t9z7fhf+5oFdSXpEz2qkI59WfO39c/GOlTW/7uX2RitSs+FlcwxPwM05gBwKg8YtGnas6p4TvJigi8uqAG3jbQP3Xhf/q8vFej1f96/ur780MPcDfeKkT8CPDMPTee+8pOTlZEcciFBXZ+UuZNxKpSL09/2099dRTqq+v9/MqAXtj8AF+Ultbq0ceeUSvvvqqtm/frlcWvKL02HQ5fXxhxSmnvt//+yr9f6Xq06eP7rzzTu3evTtAqwbsh8EH9NK1lfftb39bhw8f1rhx4yRJyTHJPg2/az+dIS4uTiUlJVq7dq3mzJlD/QF+wsEtQC/U1tZq7ty5OnnypN55553Wgdfufu5aHXId6vLz+MbFjOvwwtSNjY0qLCzUzp079fbbbysnJ8f/TwawCQYf0AOGYWjjxo165pln9NRTT2np0qWKju76nLsrn8Be563Tlq1bNO2H05QQkdDtT2AvLS3V7NmzlZWVpeLiYsXHx/vj6QC2wuADfNTdyuuKw+FQT/77UX9A7/AeH9BNN3ovL5h47w/oHQYf0A3XH7FZVFTUrZc2A2nixImqqqriyE/ARww+4AasUnmdof4A3zH4gE5YsfI6Q/0B3cfgA65j9crrDPUHdA+DD7hGKFVeZ6g/4MYYfIBCt/I6Q/0BnWPwwfbCofI6Q/0B7TH4YFvhVnmdof6Athh8sKVwrrzOUH/ANxh8sBW7VF5nqD+AwQcbsWPldYb6g50x+BD27F55naH+YFcMPoQ1Kq9r1B/shsGHsETl+Yb6g50w+BB2qLyeo/5gBww+hA0qzz+oP4Q7Bh/CApXnf9QfwhWDDyGNygss6g/hiMGHkEXlBQ/1h3DC4EPIofLMQf0hXDD4EFKoPPNRfwh1DD6EBCrPWqg/hDIGHyyPyrMu6g+hiMEHy6LyQgP1h1DD4IMlUXmhh/pDqGDwwVKovNBG/SEUMPhgGVRe+KD+YGUMPpiOygtP1B+sisEHU1F54Y/6g9Uw+GAKKs9eqD9YCYMPQUfl2Rf1Bytg8CFoqDxI1B/Mx+BDUFB5uB71B7Mw+BBQVB5uhPqDGRh8CBgqD91F/SGYGHzwOyoPPUH9IVgYfPArKg+9Rf0h0Bh88AsqD/5E/SGQGHzoNSoPgUL9IRAYfOgxKg/BQP3B3xh86BEqD8FG/cFfGHzwCZUHM/Wk/gzD0D/+8Y8grRChgMGHdv74xz/K7Xa3+zqVB6vwpf4++eQTTZo0ST//+c+DuEJYmcMwDMPsRcAaqqurVVBQoDNnzujYsWOKi4trvc3tdquoqEgul0tLly5l4PmBw+EQ//16r7S0VC+++KL27dunmJiYdrc3NjZq//79evTRR1VRUaGkpCQTVgkrYfBBkrR06VLt3r1bP/vZz1RcXKxFixZp2rRpbe7j8XgUGRlp0grDD4PPf668QuF0Otvd5vF4VFBQoAEDBqi4uDjYS4MFtf9XAtupra3VpUuXtHnzZg0bNkxffPGFPvroI02dOrXNoGPowao6GnhXrF+/Xn/9619VXl4uSfJ6vYqI4F0eO2PwQYmJiXr99ddb/9y3b19dunRJkZGR7CQQ0k6cOKGVK1fqtddeU2RkJK9aQBIHt+AaHo9HkpSbm6s//OEP+vTTTxUREcHLcQgphmHogw8+kCS98cYbSk9P18SJEyXxqgW+wXt8aOPKT8QLFizQ2LFjlZ+fb/aSwhbv8QVGY2Ojpk2bphMnTiguLk7V1dWSeIkTV/GvwOau3/Fe+Yk4NjZWdXV1Hd4HsLK4uDjt2bNH8+bN09///nd9+OGHksTQQyv+JdjYP//5T128eLHN164MudTUVH300UeSvikTINQ899xzOnz4sFwul7xeb+vXL126pCNHjpi4MpiNlzptyDAMbdy4Uc8884xee+015eXlMdxMwEud5qirq9Pdd9+trKwsFRcXKz4+3uwlIcgoPpu5/uorTzzxBEMPtpKQkMA1P22OwWcTXGMTuIpPfLA3Bp8NcI1NoGN84oM9MfjCGJUHdI36sx8GX5ii8gDfUH/2weALM1Qe0HPUnz0w+MIIlQf4B/UX3hh8YYDKA/yP+gtfDL4QR+UBgUX9hR8GX4ii8oDgof7CC4MvBFF5gDmov/DA4AshVB5gPuov9DH4QgSVB1gL9Re6GHwWR+UB1kX9hSYGn4VReUBooP5CC4PPgqg8IPRQf6GDwWcxVB4Q2qg/62PwWQSVB4QP6s/aGHwWQOUB4Yn6syYGn4moPCD8UX/Ww+AzCZUH2Av1Zx0MviCj8gD7ov6sgcEXRFQeAIn6MxuDLwioPADXo/7Mw+ALMCoPwI1Qf8HH4AsQKg9Ad1F/wcXgCwAqD0BPUH/BweDzIyoPQG9Rf4HH4PMTKg+AP1F/gcPg6yUqD0CgUH+BweDrBSoPQDBQf/7F4OsBKg9AsFF//sPg8xGVB8BM1F/vMfi6icoDYBXUX++YMviqq6s1d+5clZeXd3qfs2fPas6cOdq+fXsQV9YxKg+AFfmr/kJtn9xbQR98VVVVmj17tu69914VFhaqpKREOTk5Sk1NVWpqqiZMmKCVK1dq6tSpSklJ0YYNG7R+/Xqfv8++ffuUmZmp+Ph4jRgxokdrpfIAWF1v6y9Y+2TDMPT8888rISFBCQkJeu6552QYhs/b8QeHEeTvvGzZMiUlJWn48OHasWOHysrKtHr1ao0aNUqSVF9fr0WLFikmJkbjx49Xdna2ZsyYob179/r0fQ4ePKiTJ0/q8uXLKioqUk1NjU+Pr62t1dy5c3Xy5Em98847DDz4ncPhMO0/PsJTY2OjCgsLtXPnTr399tvKycnp8jHB2ie/9dZbKi4u1ocffiiHw6EHH3xQCxYs0Jw5c3r0XHvFCLKDBw8a3/ve94wvv/zSKCkpMbxer+F2u1tv93q9hsfjMX7/+98b586dM/Lz841f/vKX7bZTV1dnDB061Ni6dathGIbR2NhojBo1yli/fn2b+33wwQfG8OHDu70+r9dr/O53vzNuueUWY/HixYbL5erZEwW6YMJ/P9jElf3erFmzjK+++uqG9w3WPnnChAnGW2+91Xr/NWvWGOPHj/fH0/WZKf/zdu3aZYwdO9aYMGGC8dvf/tZIT09vve3FF180Fi9ebPzkJz8x7rzzTmPZsmWdbmfPnj1GYmKiUVtba8yePdt45JFH2t3Hl8F37tw5Y+rUqUZSUpJx8OBBn58X4AsGHwKpoaHBKCgoMG677TZj165dN7xvMPbJAwYMMP785z+3/vnQoUNG//79e/EMe84Z/MaUcnJylJ2dreeff17Hjx/X2LFjW2+LiorS0KFDdeDAAVVUVCgqKqrT7WRnZ2v69OnKyspSXV2dqqqqerQewzC0ceNGPfPMM3rqqaf03nvvcfAKgJB25b2/0tJSzZ49W1lZWSouLlZ8fHy7+wZjn3zx4sU23zs+Pl4XL16UYRhyOBx+etbdY9rpDBEREcrJydHWrVu1fPlyff3115Kk+fPnq7KyUpmZmTf8C74iPz9f1dXVmjlzphISEnxeB0dsAghn3T3yM9D75P79+6uhoaH1zw0NDerfv3/Qh55k8nl8Y8aM0bBhw7R582Z9/PHHkqR58+bp9ttv15gxY7p8vMfjUUFBgfLy8rR69WqdOnWq29/b4IhNADbR3SM/A7lPvuOOO3TkyJHWPx85ckR33HGHH56d70wdfEOGDNGZM2fU1NSkfv36qbKyUrm5uVq3bp3S0tK6fHxRUZEkae3atVq4cKHy8vLk8XgkSV6vVy6XSy0tLTIMQy6Xq/UnGCoPgB11VX+B3Cfn5eWpuLhYZ8+e1eeff6433nhDTz75pN+fY7eY8s7iNcrKyoyBAwcaM2fONGbNmmV861vfMl5++eUuH1dRUWHcdNNNxt/+9jfDMAzD7XYb9913n/Hqq68ahmEY+/btMyS1+fXAAw9wxCYswwL//WBjnR35Gah9stfrNQoLC42bb77ZuPnmm43CwkLD6/UG5sl1Iejn8XVkxYoV2rRpkwzDUEpKilatWqXY2Fi/fg/Oy4PVcB4fzNbZeX/B2CebyRKDL5CM647YXLp0KS9rwhIYfLCK7hz5GU5MOZ2hJ5q8TTrWfEznPefVbDQr2hGtQZGDlBSdpL4RfTt8zLWVt337dioPlvKnP/3J7CUAkq6+91dYWKg777yzW1d96ck+2SosX3zn3OdU4apQTUuNJMkjT+ttkYqUJI3oM0JpMWka7BwsicoDgJ7qqv56sk+2GksPvqOuoyq7XCa33F3e1ymn0mPTlVifyHt5ANALnb3315N9cnJMcqCX6zPLDj5f/oJbeaRd/7FLd/W7i8oDgF66tv7yX89XhSp82idbdfhZcvCdc5/T+43vt/sL/k3Bb/S3A39T86VmDUgcoH+b/2+akDehzX0ivBH6cfyPlehMDOaSASAsNTY26sVfvqhbZ92qPrF9Wr9e9usyHXzvoD4/9rnueuQuPfbmYx0+3imnpsVNs9Q+2ZKDb/vF7Trdcrrd1/9x/B/6l5H/Ime0U7Wf1mrlD1cqf2O+bvvubW3uN6rPKOX2zw3WcgEgrG2/uF2nvz4tXXN1sSPbjsgR4dCJP55Qi6ul08EnWW+fbOqVWzrS5G1qfdP0ekO+M0TO6P89ENXxzeHg5/9+vt39alpq1ORtCuAqAcAeWvfJ111SM2VyipJ/kKx+A/t1uQ2r7ZMtN/iONR+74e2bF25W4dBCLR+/XAMSByjpwaQebQcwU1VVlc6fb/9DmyQ1NzcHeTVA5/y1L7XSPtlyg++853ybw2OvN/316VpxZoUW7Fyg5NzkqwV4DY88qvPWBXKZQK888cQT6tPn6vslXq9XktTU1KSMjAyzlgW009U+uTustk+23OBrNrr+aTciMkIj7x2prz7/SuVryzu8z5atW+RwOPjFL1N+TZ48WV999VWn/4a9Xm+b86PuvvtuSVLfvn1bL6Z+4MAB058Hv/j1/rb3fdiDd67Za51XMix35ZZoR/dPQfC6var7e8c/RUz74TT92vi1v5YF+JXX69Xly5cVGxurhoYGnT59Wk1NTYqOjpbb/c3RzBkZGVzSDKbbfXG3Trac7PV2oiOsc3qZ5YpvUOSg1rP/r9X4RaMOv39YzReb5fV4dfzD4zr8+8Mak9H+M6IiFamECN8/lBYIlkcffVQPPfSQXnnlFf3oRz/SvHnzlJGRofvvv1/Tpk0ze3lAq872yR63Ry2uFnk9XhkeQy2uFnncHb8karV9suVOZ2jyNmlt/dp2rylfPH9R655cp7PVZ2V4DQ28baAy8jM04YkJ7bYRqUjNip9l+evFwd527typEydOKCcnR0lJSTp06JAMw9A999xj9tKAVp3tk3f95y7teW1Pm6899NxDenjRw+22YbV9suUGn9T5eXzdYXgNjYgYoR/d/CP/LgoAbKo3+2SJ8/i6JS0mTc4evv1ouA0tf2y5SktL/bwqALCn3uyTnXJqXIy1rplsycE32DlY6bHpPv9FO+VUVnyWXl7wsmbNmqWCggI1NDQEaJUAYA+92Senx6Zb6nJlkkUHnyQlxyT79Bd97cVQc3JyVFVVJY/Ho+TkZOoPlnSj0x0Aq+nNPtlqLPke37Vq3bU65DrU5Wc/jYsZ1+FPFbt371Z+fr4efvhh/eIXv9CAAQOCsm6gK5MnT9a2bdvMXgbgk97uk63A8oPviiuf9lvnrVOzt1nREdFKiEjo1qf91tfX69lnn1VpaanWrFmjiRMnBmnVQOccDgfn6SFk9WafbLaQGXz+QP3BShh8gDks+x5fIPDeHwDAVsV3LeoPZqP4AHPYqviuRf0BgD3ZtviuRf3BDBQfYA7bFt+1qD8AsA+K7zrUH4KF4gPMQfFdh/oDgPBG8d0A9YdAovgAc1B8N0D9AUD4ofi6ifqDv1F8gDkovm6i/gAgPFB8PUD9wR8oPsAcFF8PUH8AELoovl6i/tBTFB9gDoqvl6g/AAgtFJ8fUX/wBcUHmIPi8yPqDwCsj+ILEOoPXaH4AHNQfAFC/QGANVF8QUD9oSMUH2AOii8IqD8AsA6KL8ioP1xB8QHmoPiCjPoDAHNRfCai/uyN4gPMQfGZiPoDgOCj+CyC+rMfig8wB8VnEdQfAAQHxWdB1J89UHyAOSg+C6L+ACBwKD6Lo/7CF8UHmIPiszjqDwD8i+ILIdRfeKH4AHNQfCGE+gOA3qP4QhT1F/ooPsAcFF+Iov4AoGcovjBA/YUmig8wB8UXBqg/AOg+ii/MUH+hg+IDzEHxhRnqDwBujOILY9SftVF8gDkovjBG/QFAexSfTVB/1kPxAeag+GyC+gOAb1B8NtTd+vvVr36lxsZGvfDCC0FeoT1QfIA5KD4burb+tm7dKo/H0+4+1dXVWrFihfbu3avJkyersrIy+AsFgACg+GzO4/EoMjKy3dfvu+8+5efn68knn9SqVau0ZcsW7d69W1FRUSasMjxRfIA5KD6b62joVVdX67PPPtOJEyckSU8//bR27drF0AMQFhh8aGfs2LE6evSozp8/r+eee07Nzc0MPQBhg8GHVl9++WXr72+66SY9/fTTOnLkiKRvXpYDgHDA4EOrzZs3KysrS7W1tZKkv/zlL/rOd76jiAj+mQAIH06zFwDryM/PV0NDgzIzM5WSkqLGxkbNnTtXffr0MXtpAOA3HNWJdk6fPq3jx48rLS1NiYmJcjgc8nq9qqurU3R0NFd98ROO6gTMwWtYaGfUqFHKzc3V4MGDW9/bczgcKikp4aovAEIexQefcM1P/6H4AHNQfPAJ1/wEEOooPvQY9dc7FB9gDooPPUb9AQhFFB/8gvrzHcUHmIPig19QfwBCBcUHv6P+uofiA8xB8cHvqD8AVkbxIaCov85RfIA5KD4EFPUHwGooPgQN9dcWxQeYg+JD0FB/AKyA4oMpqD+KDzALxQdTUH8AzELxwXR2rT+KDzAHxQfTUX8Agonig6XYqf4oPsAcFB8shfoDEGgUHywr3OuP4gPMQfHBsqg/AIFA8SEkhGP9UXyAOSg+hATqD4C/UHwIOeFSfxQfYA6KDyGH+gPQGxQfQloo1x/FB5iD4kNIo/4A+IriQ9gItfqj+ABzUHwIG9QfgO6g+BCWQqH+KD7AHBQfwhL1B6AzFB/CnlXrj+IDzEHxIexRfwCuRfHBVqxUfxQfYA6KD7ZC/QFg8MF24uPjtWbNGpWUlGjWrFkqKChQQ0NDj7ZVXV2tuXPnqry8vNP7nD17VnPmzNH27dt7umQAfsTgg231tv6qqqo0e/Zs3XvvvSosLFRJSYlycnKUmpqq1NRUTZgwQStXrtTUqVOVkpKiDRs2aP369T6vc9++fcrMzFR8fLxGjBjh8+MBtMV7fIB69t7fsmXLlJSUpOHDh2vHjh0qKyvT6tWrNWrUKElSfX29Fi1apJiYGI0fP17Z2dmaMWOG9u7dK6n77/EdPHhQJ0+e1OXLl1VUVKSamppePVfA7ig+QD2rv+zsbBUXF2vkyJG65ZZbtGfPnjZFNmDAAL355pvKyMhQZmamFi9erNzc3HbbuXDhgm699VZt27ZNknTx4kWNHj1a7777riTpnnvu0eOPP66RI0f658kCNkfxAdfxpf52796twsJCxcXF6ac//alKSkp04MABSdJLL70kt9utmpoaVVdXa8aMGXrhhRdaH3tt8e3du1d5eXk6evSolixZoi+//FJbtmxp871KS0s1e/Zsig/oJYoPuI4v9ZeTk6MjR47o/vvv1/HjxzV27NjW26KiojR8+HCdOXNGFRUVbYbe9bKzszV9+nRlZWVpx44deuutt/z6nABcxeADOuDLkZ8RERHKycnR1q1btXz5cn399deSpPnz56uyslKZmZmKiorq8nvm5+erurpaM2fOVEJCgl+fD4CrGHzADXS3/saMGaNhw4Zp8+bN+vjjjyVJ8+bN0+23364xY8Z0+X08Ho8KCgqUl5en1atX69SpU359HgCuYvABXehO/Q0ZMkRnzpxRU1OT+vXrp8rKSuXm5mrdunVKS0vr8nsUFRVJktauXauFCxcqLy9PHo9HkuT1euVyudTS0iLDMORyuVqrEoDvOLgF8EF9fb2effZZlZaWas2aNZo4cWLrbeXl5ZoyZYqmTJkih8Oh8vJyPfbYY3rppZc63NaVg1s++eQTTZw4UYcOHdLo0aPl8XiUkZGhSZMmacmSJdq/f78yMzPbPPaBBx7Q/v37A/lUgbDF4AN6oLMjP1esWKFNmzbJMAylpKRo1apVio2N7XAbXKsTMAeDD+ihG9VfdzD4AHMw+IBe8uW8vyZvk441H9N5z3m9v+19PTL5EQ2KHKSk6CT1jegbxFUD9sXgA/ygq/o75z6nCleFalpqJEkeeVpvi1SkJGlEnxFKi0nTYOfgoK0bsCMGH+BHHdXfUddRlV0uk1vuLh/vlFPpselKjkkOwmoBe+J0BsCPrj/vb1Plpm4PPUlyy62yy2U66joa4JUC9kXxAQHy32X/rdOjTqtPbJ8Ob//i9Bda8b0VSvlhih5/6/E2tznl1LS4aUp0JgZjqYCtUHxAgPRJ7dPp0JOkLYVbNCx1WIe3ueXWIdehQC0NsDUGHxAATd6m1gNZOnL4/cOKjY/VmIzOL2dW01KjJm9TAFYH2BuDDwiAY83HOr3N1eDSrv/cpSn/MaVX2wHQMww+IADOe863OWXhWjuLdmr8/x2vm2+9+Ybb8MijOm9dIJYH2BqDDwiAZqO5w69/VvWZPv3Tp/r+3O93bzvejrcDoOecZi8ACEfRjugOv36q/JQu/M8FvZz8siSp+VKzDI+h10++roX7F7bfTkTH2wHQcww+IAAGRQ7SqZZT7V7uvO+J+3TX/7mr9c/7Vu7ThTMXNP2N6e22EalIJUTwgbSAvzH4gABIik7Sn11/bvf1qL5Riup79dPYo/pFyRnjVP9B/TvdDgD/4gR2IEC2X9yu0y2ne/z4UX1GKbd/rh9XBEDi4BYgYNJi0uTs4YsqTjk1Lmacn1cEQGLwAQEz2DlY6bHpPg+/Kxeq5nJlQGAw+IAASo5J9mn48ekMQODxHh8QBLXuWh1yHery8/jGxYyj9IAAY/ABQXTlE9jrvHVq9jYrOiJaCREJfAI7EEQMPgCArfAeHwDAVhh8AABbYfABAGyFwQcAsBUGHwDAVhh8AABbYfABAGyFwQcAsBUGHwDAVhh8AABbYfABAGyFwQcAsBUGHwDAVhh8AABbYfABAGyFwQcAsBUGHwDAVhh8AABbYfABAGyFwQcAsBUGHwDAVv4/TYTmPEozjg8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "draw_graph(read_input())" ] }, { "cell_type": "markdown", "id": "4cc0e474", "metadata": {}, "source": [ "Tougher inputs will be tested. [Test case 9](https://raw.githubusercontent.com/wilsjame/weblog/main/docs/source/_ipynb/butter/butter9.in) is the toughest with 500 cows, 800 pastures, and 1450 connections.\n", "\n", "Here is test case 9 visualized." ] }, { "cell_type": "code", "execution_count": 6, "id": "86e5440d", "metadata": {}, "outputs": [], "source": [ "# draw_graph(read_input('butter9.in'), is_large=True)" ] }, { "cell_type": "markdown", "id": "ef2743fe", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "id": "878eb2e4", "metadata": {}, "source": [ "Drawing the graph takes too long to load in some notebooks. Uncomment the line to try it or download the full-size image [butter.png](https://github.com/wilsjame/weblog/blob/main/docs/source/_ipynb/butter/butter9.png) (18 MB).\n", "\n", "## Output\n", "\n", "Find the minimum distance cows must walk to a pasture with a sugar cube.\n", "\n", "For the sample input shown above the answer is: **8**\n", "\n", "Putting the cube in pasture 4 means: cow 1 walks 3 units; cow 2 walks 5 units; cow 3 walks 0 units -- a total of 8.\n", "\n", "For test case 9 the answer is: **164290**\n", "\n", "## Shortest Paths: Dijkstra\n", "\n", "Dijkstra's algorithm finds the shortest paths from a source node to all nodes in a graph. It's greedy and facilitates this with a priority queue. The cool thing is every node needs only be processed once. Assuming we have the graph, three data structures are used:\n", "\n", "- An array `distance[]` stores the distance from the source to every node. Every node's initial distance is infinity.\n", "- An array `processed[]` keeps track of processed nodes. Once a node is processed it cannot be visited again and its distance is final. Initially false for every node.\n", "- A priority queue `min_pq` provides the node to process next. When a node is processed its neighbors are added. Nodes are ordered by their distance value.\n", "\n", "Set the source node distance `distance[source] = 0` and add the source and its distance to the priority queue `min_pq.push((0, source))`. Now, enter a loop until the queue is empty. The loop pops the minimum distance node from the queue `cur = min_pq.pop()` and processes it `processed[cur] = True`. Sweet one less node to go! Iterate the current node's neighbors. Skip neighbors that have been processed. For neighbors not yet processed we update their distances and add them to the queue. It's greedy `distance[neighbor] = min(distance[neighbor], distance[cur] + cur_to_neighbor_edge_weight)` and `min_pq.push(distance[neighbor], neighbor)`. When the queue is empty all nodes are processed and the loop terminates. `distance[]` holds the shortest paths for each node to the source. The worst-case performance is O(E\\*log(N)) where E is the number of edges and N is the number of nodes. We visit each node at most E times and for each visit, the priority queue has at most N nodes. " ] }, { "cell_type": "code", "execution_count": 7, "id": "be6dcee1", "metadata": {}, "outputs": [], "source": [ "def dijkstra(source: int, data: Data) -> Dist:\n", " \"\"\"Find shortest paths from source to each node in adjacency graph.\"\"\"\n", " _, adj = data\n", " nodes = len(adj)\n", " distance = dict.fromkeys([i for i in range(1, nodes + 1)], float(\"inf\"))\n", " processed = dict.fromkeys([i for i in range(1, nodes + 1)], False)\n", " min_pq = []\n", " distance[source] = 0\n", " heapq.heappush(min_pq, (0, source)) # distance before node for heap ordering\n", " while len(min_pq) > 0:\n", " _, cur = heapq.heappop(min_pq)\n", " processed[cur] = True\n", " for neighbor in adj[cur]:\n", " node, cur_to_node_dist = neighbor # node before distance in adjacency list\n", " if processed[node] is True:\n", " continue\n", " distance[node] = min(distance[node], distance[cur] + cur_to_node_dist)\n", " heapq.heappush(min_pq, (distance[node], node)) # distance before node for heap ordering\n", " return distance" ] }, { "cell_type": "markdown", "id": "b7b79a83", "metadata": {}, "source": [ "## Solve\n", "\n", "We try each node as the sugar cube source to find the minimum distance cows must walk. Running time is O(C\\*E\\*log(N)) where C is cows, E edges, and N nodes. " ] }, { "cell_type": "code", "execution_count": 8, "id": "484ff6b8", "metadata": {}, "outputs": [], "source": [ "def solve(data: Data) -> int:\n", " \"\"\"Find the minimum distance cows must walk.\"\"\"\n", " cow, adj = data\n", " res = float(\"inf\")\n", " nodes = len(adj)\n", " for source in range(1, nodes + 1):\n", " shortest_paths = dijkstra(source, data)\n", " total = 0\n", " for i in cow:\n", " total += shortest_paths[cow[i]]\n", " res = min(res, total)\n", " return res" ] }, { "cell_type": "code", "execution_count": 9, "id": "47d04a83", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solve(read_input()) == 8" ] }, { "cell_type": "code", "execution_count": 10, "id": "6be6a6d9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Wall time: 6.34 s\n" ] }, { "data": { "text/plain": [ "True" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time solve(read_input('butter9.in')) == 164290" ] }, { "cell_type": "markdown", "id": "9f9493d6", "metadata": {}, "source": [ "### Go Fast\n", "\n", "Test case 9 isn't *the* toughest input. There is a test case 10! The USACO grader allotted runtime is 1 second. The approach discussed here has the correct run time complexity O(C\\*E\\*log(N)), but the implementation falls short and fails the <1 second mark. \n", "\n", "![test case 9 py](usaco_grader_py.png)\n", "\n", "For starters, we can speed up our Dijkstra implementation. Sometimes the queue might pop an already processed node. When this happens there is no point processing it again because its shortest distance is already set. An if statement should do the trick." ] }, { "cell_type": "code", "execution_count": 11, "id": "e860aa20", "metadata": {}, "outputs": [], "source": [ "def dijkstra(source: int, data: Data) -> Dist:\n", " \"\"\"Find shortest paths from source to each node in adjacency graph.\"\"\"\n", " _, adj = data\n", " nodes = len(adj)\n", " distance = dict.fromkeys([i for i in range(1, nodes + 1)], float(\"inf\"))\n", " processed = dict.fromkeys([i for i in range(1, nodes + 1)], False)\n", " min_pq = []\n", " distance[source] = 0\n", " heapq.heappush(min_pq, (0, source)) # distance before node for heap ordering\n", " while len(min_pq) > 0:\n", " _, cur = heapq.heappop(min_pq)\n", " if processed[cur] == False: # <------------------------- speed up \n", " processed[cur] = True\n", " for neighbor in adj[cur]:\n", " node, cur_to_node_dist = neighbor # node before distance in adjacency list\n", " if processed[node] is True:\n", " continue\n", " distance[node] = min(distance[node], distance[cur] + cur_to_node_dist)\n", " heapq.heappush(min_pq, (distance[node], node)) # distance before node for heap ordering\n", " return distance" ] }, { "cell_type": "code", "execution_count": 12, "id": "8964c510", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Wall time: 1.39 s\n" ] }, { "data": { "text/plain": [ "True" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time solve(read_input('butter9.in')) == 164290" ] }, { "cell_type": "markdown", "id": "0fe90006", "metadata": {}, "source": [ "Sweet! A considerable speed up. Additional improvements could be writing your own Fibonacci heap discussed on the [Dijkstra\\'s algorithm](https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm) wiki page or storing the graph as an adjacency matrix. Trying to optimize look-ups by using list indices and values in place of a dict hashing keys to values. Implementing other search algorithms. And of course, experimenting with different languages. \n", "\n", "![test case 10 cpp](usaco_grader_cpp.png \"On average, PyPy is 4.2 times faster than CPython!\")\n", "\n", "Don't forget, on average, PyPy is 4.2 times faster than CPython! " ] } ], "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.8.8" } }, "nbformat": 4, "nbformat_minor": 5 }