{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Numeric Optimization of cost functions of 1 parameter" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# do the usual imports\n", "import numpy as np\n", "import scipy.signal\n", "import scipy.integrate\n", "import scipy.optimize\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def hill_descent_1D( fn, p0 ):\n", " ''' Find the parameter that minimizes the function fn\n", " where fn is a function of just one parameter.\n", " p0 is the initial guess for the parameter.\n", " The approach taken to find the minimal value of fn is \n", " a hill-descending approach.\n", " Returns the best parameter value and its cost (as a tuple)'''\n", " # The current value of the parameter\n", " currP = p0\n", " # How much we will adjust the parameter by to take a setp\n", " delP = 0.05*p0\n", " # how much we will adjust the size of our next parameter step\n", " accel = 1.2\n", " # how close the costs need to be before we stop taking steps.\n", " tol = 1e-4\n", " # candidate step sizes (small and big, negative and positive)\n", " candidates = np.array([-accel,-1/accel,1/accel,accel])\n", " while True:\n", " currCost = fn(currP) # evaluate the cost at the current parameter value\n", " # then evaluate at 2 steps to the \"left\" and 2 steps to the \"right\"\n", " costs = np.array(list(map(fn,currP+candidates*delP)))\n", " # find the minium cost and the index of which step took us too it\n", " mCost = costs.min()\n", " mIdx = np.argmin(costs)\n", " if currCost < mCost:\n", " # no step made the cost any smaller. we should take a smaller step and try again\n", " delP = delP/accel\n", " else:\n", " # update the parameter to the best step taken\n", " currP = currP + candidates[mIdx]*delP\n", " # and adjust the step size\n", " delP = delP*candidates[mIdx]\n", " if np.abs(currCost-mCost) < tol:\n", " # if the last step and this step have costs within the tolerance,\n", " # then it is time to stop looping.\n", " return (currP, mCost)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# hill-climbing on quadratic function.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Does the hill climbing algorithm work for abs value function?\n", "# result: sort of. it does take more steps, presumably because the values don't flatten out at the bottom.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Van der Pol Oscillator\n", "\n", "Define the model function for the van der Pol oscillator. It is defines as\n", "\n", "\\begin{align*}\n", "\\frac{dy_1}{dt} &= y_2\\\\\n", "\\frac{dy_2}{dt} &= \\mu (1-y_1^2)y_2 - y_1.\\\\\n", "\\end{align*}\n", "\n", "We will assume the unit of time is $h$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def vdp( t, y, mu ):\n", " y1 = y[0]\n", " y2 = y[1]\n", " dydt = np.zeros((2,))\n", " dydt[0] = y2\n", " dydt[1] = mu*(1-y1**2)*y2 - y1\n", " return dydt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the model\n", "\n", "Use $\\mu$=0.1. Run it for 500 h at a tight tolerance. Then run it again for 100 h and plot the result." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### define get_period\n", "Stephanie copied her code from HW2." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def get_period(t,x):\n", " \"\"\" Approximate the period of a 1-D x, given the time-steps t.\n", " Returns a tuple with the period and the standard deviation of the period over time.\n", " if the value of the standard deviation is not smaller than 0.1, then\n", " it means the period estimate is dodgy and you shouldn't use it. Instead,\n", " plot your simulation and figure out why it isn't periodic - maybe it just\n", " hasn't reached the limit cycle yet.\"\"\"\n", " idxs = scipy.signal.find_peaks(x)\n", " idxs = idxs[0]\n", " times = t[idxs]\n", " period = np.mean(np.diff(times));\n", " sdperiod = np.sqrt(np.var(np.diff(times)));\n", " return (period,sdperiod)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define the van der Pol error function\n", "\n", "It should take the parameter $\\mu$ as input and return the squared difference between the simulated period and 10 h. \n", "\n", "To find the simulated period, first run the model for 200 hours and ignore the result. The simulate for another 50 h and call get_period." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the cost function landscape\n", "\n", "Plot the cost over 100 values of $\\mu$ ranging from 0.01 to 5.\n", "\n", "Then plot the result from the Hill descent method beginning at $\\mu=0.2$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do built-in local methods do any better? Not really." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What if we have multiple parameters?\n", "\n", "Let's code up the Goodwin Oscillator" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def goodwin_oscillator( t, statevars, params ):\n", " X = statevars[0];\n", " Y = statevars[1];\n", " Z = statevars[2];\n", " alpha_1 = params[0]\n", " K = params[1]\n", " n = params[2]\n", " gamma_1 = params[3]\n", " alpha_2 = params[4]\n", " gamma_2 = params[5]\n", " alpha_3 = params[6]\n", " gamma_3 = params[7]\n", " \n", " dydt = np.zeros( (3,) )\n", " dydt[0] = alpha_1*K**n/(K**n+Z**n) - gamma_1*X\n", " dydt[1] = alpha_2*X - gamma_2*Y;\n", " dydt[2] = alpha_3*Y - gamma_3*Z;\n", "\n", " return dydt\n", "\n", "def goodwin_circadian_error( params ):\n", " mf = lambda t,y : goodwin_oscillator(t,y, params)\n", "\n", " sol1 = scipy.integrate.solve_ivp( fun=mf, t_span=(0,500), y0=(.1,.1,.1) )\n", " sol2 = scipy.integrate.solve_ivp( fun=mf, t_span=(0,100), y0=sol1.y[:,-1], \n", " t_eval=np.arange(0,100,.01) )\n", " (period,sdper) = get_period( sol2.t, sol2.y[0,:] )\n", " return (period-24)**2\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "alpha_1 = 1\n", "K = 0.61\n", "n = 11\n", "gamma_1 = 0.15\n", "alpha_2 = 1\n", "gamma_2 = 0.15\n", "alpha_3 = 1\n", "gamma_3 = 0.25\n", "\n", "params = (alpha_1,K,n,gamma_1,alpha_2,gamma_2,alpha_3,gamma_3)\n", "\n", "mf = lambda t,y : goodwin_oscillator(t,y, params)\n", "\n", "sol1 = scipy.integrate.solve_ivp( fun=mf, t_span=(0,200), y0=(.1,.1,.1), \n", " t_eval=np.arange(0,200,.1) )\n", "sol2 = scipy.integrate.solve_ivp( fun=mf, t_span=(0,100), y0=sol1.y[:,-1], \n", " t_eval=np.arange(0,100,.1) )\n", "\n", "plt.plot( sol2.t, sol2.y.T )\n", "(period,sdperiod) = get_period( sol2.t, sol2.y[0,:])\n", "plt.title( 'Period = %f' % period );\n", "\n", "err = goodwin_circadian_error( params )\n", "print( err )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Run a built-in optimizer and print the results\n", "res = scipy.optimize.minimize( goodwin_circadian_error, method='Nelder-Mead', x0=params )\n", "print(res)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Use the results from the optimization to run the model\n", "print( goodwin_circadian_error( res.x ))\n", "\n", "mf = lambda t,y : goodwin_oscillator(t,y, res.x)\n", "\n", "sol1 = scipy.integrate.solve_ivp( fun=mf, t_span=(0,200), y0=(.1,.1,.1), \n", " t_eval=np.arange(0,200,.1) )\n", "sol2 = scipy.integrate.solve_ivp( fun=mf, t_span=(0,100), y0=sol1.y[:,-1], \n", " t_eval=np.arange(0,100,.1) )\n", "\n", "plt.plot( sol2.t, sol2.y.T )\n", "(period,sdperiod) = get_period( sol2.t, sol2.y[0,:])\n", "plt.title( 'Period = %f' % period );\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.5" } }, "nbformat": 4, "nbformat_minor": 4 }