diff --git a/users/awf/linterp_and_monotonic.ipynb b/users/awf/linterp_and_monotonic.ipynb new file mode 100644 index 000000000..611328c51 --- /dev/null +++ b/users/awf/linterp_and_monotonic.ipynb @@ -0,0 +1,369 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "#!conda install -y ipympl\n", + "\n", + "import torch\n", + "import numpy as np\n", + "\n", + "%matplotlib ipympl\n", + "\n", + "import matplotlib.pyplot as plt\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …", + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "bb8d6f0e7c5b447c81fca376986d44af" + } + }, + "metadata": {} + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "[]" + ] + }, + "metadata": {}, + "execution_count": 2 + } + ], + "source": [ + "\n", + "x_dense = torch.arange(0,3.01,0.1)\n", + "f = lambda x: x ** 2 - 0.2 * x**3\n", + "plt.plot(x_dense, f(x_dense))\n", + "\n", + "n_samples = 10\n", + "samples_x = torch.rand(n_samples) * 2.8 + 0.1\n", + "\n", + "# Sample, but lower\n", + "samples_y = f(samples_x)\n", + "\n", + "plt.plot(samples_x, samples_y, '.')\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …", + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "35375136dd694f1db1dff11bffd5f974" + } + }, + "metadata": {} + } + ], + "source": [ + "# Grab from \n", + "def interp1d(x, y, xnew):\n", + " \"\"\"\n", + " PyTorch interp1d, as in scipy.interpolate(x,y,assume_sorted=True)(xnew)\n", + " \n", + " Limited to 1D-1D for now, should not be hard to enhance if needed\n", + " \n", + " github.com/awf\n", + " \"\"\"\n", + " # For each point in x, we want to find index of the knot above it, and hence that below it\n", + " # so knots_x[ind-1] <= x < knots_x[ind]\n", + " inds = torch.bucketize(xnew, x, right=True)\n", + "\n", + " # Call those points xlo, xhi\n", + " xlo = x[inds-1]\n", + " xhi = x[inds]\n", + " ylo = y[inds-1]\n", + " yhi = y[inds]\n", + " \n", + " dx = xhi - xlo\n", + " dy = yhi - ylo\n", + "\n", + " # Then t = (xnew - xlo)/dx\n", + " # ynew = ylo + t * dy\n", + " t = (xnew - xlo) / dx\n", + " return ylo + t * dy\n", + "\n", + "def test_interp1d():\n", + " # scipy example from https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html\n", + " import scipy.interpolate\n", + " import numpy as np\n", + " x = np.arange(0, 10)\n", + " y = np.exp(-x/3.0)\n", + " f = scipy.interpolate.interp1d(x, y)\n", + " xnew = np.arange(0, 9, 0.1)\n", + " ynew = f(xnew) # use interpolation function returned by `interp1d`\n", + "\n", + " # Now do it with torch\n", + " tx = torch.arange(0, 10)\n", + " ty = torch.exp(-tx/3.0)\n", + " txnew = torch.arange(0, 9, 0.1)\n", + " tynew = interp1d(tx,ty,txnew)\n", + " plt.figure()\n", + " plt.plot(x, y, 'ro', xnew, ynew, 'r-',\n", + " tx, ty, 'kx', txnew, tynew, 'k:')\n", + " plt.show()\n", + " \n", + "test_interp1d()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …", + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "d7833134ba9a4b5f82281af887dc1700" + } + }, + "metadata": {} + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "SGD (\n", + "Parameter Group 0\n", + " dampening: 0\n", + " lr: 0.01\n", + " momentum: 0.9\n", + " nesterov: False\n", + " weight_decay: 0\n", + ")" + ] + }, + "metadata": {}, + "execution_count": 4 + } + ], + "source": [ + "\n", + "\n", + "# Fitting. Define approximator.\n", + "knots_x = torch.linspace(0.0,3.01,6)\n", + "\n", + "# Model is piecewise linear interpolation between knots (\"linear spline\")\n", + "model = lambda weights, x: interp1d(knots_x, weights, x)\n", + "\n", + "# Define losses\n", + "\n", + "# 1. Data term: function should agree with samples\n", + "loss_data = lambda weights: (model(weights, samples_x) - samples_y).abs().sum()\n", + "\n", + "# 2. Monotonicity: pay a penalty wherever weights[i+1] < weights[i]\n", + "loss_monotonic = lambda weights: torch.relu(weights[:-1] - weights[1:]).sum()\n", + "\n", + "loss = lambda weights: loss_data(weights) + loss_monotonic(weights)\n", + "\n", + "\n", + "# Initial estimate of parameters\n", + "torch.manual_seed(41)\n", + "weights = torch.cumsum(torch.randn(knots_x.shape),0)\n", + "\n", + "weights.requires_grad = True\n", + "\n", + "fig = plt.figure()\n", + "with torch.no_grad():\n", + " samples_pred_y = model(weights, samples_x)\n", + " plt.plot(samples_x, samples_y, '.', samples_x, samples_pred_y, 'x')\n", + " y_dense = model(weights, x_dense)\n", + " plt.plot(x_dense, y_dense, '-', knots_x, weights, 'o')\n", + "\n", + "\n", + "# Initialize optimizer\n", + "\n", + "#optimizer = torch.optim.Adam((weights,), lr=0.01, betas=(0.9, 0.999), eps=1e-08, weight_decay=0, amsgrad=False)\n", + "optimizer = torch.optim.SGD((weights,), lr=0.01,momentum=0.9)\n", + "iter_count = 0\n", + "optimizer" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "for k in range(150):\n", + " # Run one iteration\n", + " iter_count += 1\n", + " optimizer.zero_grad()\n", + " loss_value = loss(weights)\n", + " loss_value.backward()\n", + " optimizer.step()\n", + "\n", + " with torch.no_grad():\n", + " fig.clear()\n", + " plt.plot(samples_x, samples_y, 'r.')\n", + " plt.plot(x_dense, model(weights, x_dense), 'b-', knots_x, weights, 'bo')\n", + " plt.title(f'iter {iter_count} loss {loss_value:.2f} {optimizer.__class__}')\n", + " plt.show()\n", + " fig.canvas.draw()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …", + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "7fa7b95ef085468f944fa70fbc214a1c" + } + }, + "metadata": {} + } + ], + "source": [ + "# One-sided Fitting. Same model as above, but redefine here anyway\n", + "n_knots = 16\n", + "\n", + "# Fitting. Define approximator.\n", + "knots_x = torch.linspace(0.0,3.01,n_knots)\n", + "\n", + "# Model is piecewise linear interpolation between knots (\"linear spline\")\n", + "model = lambda weights, x: interp1d(knots_x, weights, x)\n", + "\n", + "# Define losses\n", + "\n", + "# 1. Data term: function should not exceed samples: relu(data - model)\n", + "p = 2 # p-norm loss on one-sided term. In principle, we should pay infinity for violating the bound, but quadratic should be enough\n", + "loss_data = lambda weights: (torch.relu(samples_y - model(weights, samples_x)) ** p).sum()\n", + "\n", + "# 2. Monotonicity: pay a penalty wherever weights[i] > weights[i+1]\n", + "loss_monotonic = lambda weights: torch.relu(weights[:-1] - weights[1:]).sum()\n", + "\n", + "weight_decay = 0.1\n", + "\n", + "# 3. Weight decay: \"gravity\" to pull model down\n", + "loss_weight_decay = lambda weights: 0.1 * weights.square().sum()\n", + "\n", + "loss = lambda weights: loss_data(weights) + loss_monotonic(weights) + loss_weight_decay(weights)\n", + "\n", + "# Initial estimate of parameters\n", + "torch.manual_seed(41)\n", + "weights = torch.cumsum(torch.randn(knots_x.shape),0)\n", + "\n", + "# Parameters => requires gradient\n", + "weights.requires_grad = True\n", + "\n", + "# Plot the samples, their model predictions, and the model parameters (i.e. spline knots)\n", + "fig = plt.figure()\n", + "with torch.no_grad():\n", + " samples_pred_y = model(weights, samples_x)\n", + " plt.plot(samples_x, samples_y, '.', samples_x, samples_pred_y, 'x')\n", + " y_dense = model(weights, x_dense)\n", + " plt.plot(x_dense, y_dense, '-', knots_x, weights, 'o')\n", + "fig.canvas.draw()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Initialize optimizer: SGD/LBFGS are the only ones with decent behaviour\n", + "optimizer = torch.optim.LBFGS\n", + "max_iter = 150\n", + "if optimizer == torch.optim.Adam:\n", + " optimizer = optimizer((weights,), lr=0.01, betas=(0.9, 0.999), eps=1e-08, weight_decay=0, amsgrad=False)\n", + "if optimizer == torch.optim.Adagrad:\n", + " optimizer = optimizer((weights,), lr=1.0, lr_decay=0, weight_decay=0, initial_accumulator_value=0, eps=1e-10)\n", + "if optimizer == torch.optim.SGD:\n", + " optimizer = optimizer((weights,), lr=0.01,momentum=0.9)\n", + "if optimizer == torch.optim.LBFGS:\n", + " max_iter = 150\n", + " optimizer = optimizer((weights,), lr=.1, max_iter=20, max_eval=None, \n", + " tolerance_grad=1e-07, tolerance_change=1e-09, \n", + " history_size=4, line_search_fn=None)\n", + "iter_count = 0\n", + "\n", + "weights_best = weights\n", + "loss_best = 1e20\n", + "best_iter = iter_count\n", + "for k in range(max_iter):\n", + " # Run one iteration\n", + " iter_count += 1\n", + " def closure():\n", + " optimizer.zero_grad()\n", + " loss_value = loss(weights)\n", + " loss_value.backward()\n", + " return loss_value\n", + " \n", + " optimizer.step(closure)\n", + "\n", + " with torch.no_grad():\n", + " loss_value = loss(weights)\n", + " if loss_value < loss_best:\n", + " weights_best = weights.clone().detach()\n", + " loss_best = loss_value\n", + " best_iter = iter_count\n", + " fig.clear()\n", + " plt.plot(samples_x, samples_y, 'r.', samples_x, model(weights_best, samples_x), 'rx')\n", + " plt.plot(x_dense, model(weights_best, x_dense), 'b-', knots_x, weights_best, 'bo')\n", + " plt.title(f'iter {iter_count}, loss {loss_best:.3f} at iter {best_iter}, weight_decay={weight_decay}')\n", + " plt.show()\n", + " fig.canvas.draw()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "name": "python388jvsc74a57bd021ad436eab6efc5a79d7ba41ee8a80a7e2ac2e5a53e8d2a30f44497713f6797f", + "display_name": "Python 3.8.8 64-bit ('ptnightly': conda)" + }, + "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" + }, + "metadata": { + "interpreter": { + "hash": "21ad436eab6efc5a79d7ba41ee8a80a7e2ac2e5a53e8d2a30f44497713f6797f" + } + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file