{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Genetic Algorithm \n", "\n", "In this notebook, we implement a new selection operator (tournament selection) and a new cost function (for the Goodwin oscillator). At the end, there is also code to plot the parameters as they change over the generations." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import os.path\n", "import scipy.linalg\n", "import scipy.integrate\n", "import scipy.signal\n", "import matplotlib.cm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Support function for GA\n", "Copy-paste your copy of generateIndividual, generateChild, and GA below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Update code to include linear ranking selection\n", "\n", "Copy-paste your GASettings definition and add to it a field named linearRankingRatio. Give it a default value of 4." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Copy-paste your select_pool below. Add to it support for linear ranking selection (call it \"linearRanking\").\n", "\n", "In linear rankng selection, we compute a discrete probability distribution with the probability of each individual in the previous generation being chosen for the breeding pool. The probabilities are constructed so that the probability of choosing the lowest-cost individual is settings.linearRankingRatio times the probability of choosing the highest cost individual. For each individual we want to add to the pool we use this probability distribution by choosing a number from a uniform distribution between 0 and 1 and then determining which index in the cumulative distribution is the smallest such that the cumulative distribution is larger than the random number.\n", "\n", "The algorithm is this. If $n$ is the number of individuals in the previous generation and $r$ is the linearRankingRatio, then we construct the probability distribution according to\n", "\n", "\n", "\n", "And selection individuals in the breeding pool\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test the selection operators\n", "\n", "This includes code that demonstrates the distribution of costs in the pool when selecting with linear ranking selection (ratios of 4 and 8) and tournament selection (tournament sizes of 2 and 5)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Test the selection operators\n", "# Make fake parameters, using the index as the values (so we can differentiate them)\n", "lam = 10000 # number of individuals in generation\n", "mu = 1000 # number of individuals in breeding pool\n", "NP = 5 # num of parameters per individual\n", "G = np.zeros( (lam,NP) )\n", "for i in range(lam):\n", " G[i,:] = i\n", "Gcost = np.linspace(0.1,500,lam)\n", "settings = GASettings(numParents=mu, numChildren=lam)\n", "\n", "settings.selection = 'linearRanking'\n", "settings.linearRankingRatio = 4\n", "(P, Pcost) = select_pool(G, Gcost, settings)\n", "plt.figure( figsize=(8,4) )\n", "plt.hist( Pcost );\n", "plt.title( \"Breeding Pool Chosen With Linear Ranking Selection (r=4)\")\n", "plt.xlabel( 'Cost')\n", "plt.ylabel('Frequency')\n", "\n", "settings.selection = 'linearRanking'\n", "settings.linearRankingRatio = 8\n", "(P, Pcost) = select_pool(G, Gcost, settings)\n", "plt.figure( figsize=(8,4) )\n", "plt.hist( Pcost );\n", "plt.title( \"Breeding Pool Chosen With Linear Ranking Selection (r=8)\")\n", "plt.xlabel( 'Cost')\n", "plt.ylabel('Frequency')\n", "\n", "settings.selection = 'tournament'\n", "settings.tournamentSize = 2\n", "(P, Pcost) = select_pool(G, Gcost, settings)\n", "plt.figure( figsize=(8,4) )\n", "plt.hist( Pcost );\n", "plt.title( \"Breeding Pool Chosen with Tournament Selection (k=2)\")\n", "plt.xlabel( 'Cost')\n", "plt.ylabel('Frequency')\n", "\n", "settings.selection = 'tournament'\n", "settings.tournamentSize = 5\n", "(P, Pcost) = select_pool(G, Gcost, settings)\n", "plt.figure( figsize=(8,4) )\n", "plt.hist( Pcost );\n", "plt.title( \"Breeding Pool Chosen with Tournament Selection (k=6)\")\n", "plt.xlabel( 'Cost')\n", "plt.ylabel('Frequency');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test the GA on the pattern cost function\n", "\n", "Test it with the pattern cost function. Plot the costs for all individuals in all the generations." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def pattern_cost_function( params ):\n", " return scipy.linalg.norm( params-np.array([1,5,2.5,2,4,3,0]) )\n", "\n", "\n", "PN = 7\n", "lb = np.zeros( (PN,) );\n", "ub = 5 + np.zeros( (PN,) );\n", "cost_fcn = pattern_cost_function;\n", "settings = GASettings(numParents=10, numChildren=50)\n", "settings.numGenerations = 5;\n", "settings.eliteCount = 1;\n", "settings.selection = 'linearRanking';\n", "settings.linearRankingRatio = 8;\n", "\n", "(params, cost, all_params, all_costs) = GA(pattern_cost_function, lb, ub, settings);\n", "\n", "plt.figure( figsize=(8,4) )\n", "legs = []\n", "for i in range(len(all_costs)):\n", " plt.plot( np.arange(all_costs[i].size), all_costs[i], '*' )\n", " legs.append( \"Generation\" + str(i) )\n", "plt.legend( legs );" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cost function for 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 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)\n", "\n", "def goodwin_circadian_error( params, debug=False ):\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", " # Compute the period of the first state. It is a scalar.\n", " (period,sdper) = get_period( sol2.t, sol2.y[0,:] )\n", " # compute the peak-to-trough amp for each state\n", " # by computing the min and max over time and subtracting.\n", " # amps will be an ndarray of length 3.\n", " amps = sol2.y.max(axis=1)-sol2.y.min(axis=1)\n", " # an amplitude larger than 0.1 is going to have low cost\n", " # scale the amplitudes and then compute the exponential for\n", " # each amp. Since log(0.001)<0, the curve will be decreasing\n", " # as the amplitude is increasing. \n", " # amp_errvals will be an ndarray of length 3.\n", " rate = np.log(0.001)/0.1;\n", " amp_errvals = np.exp(np.multiply(rate,amps));\n", " # The amplitude cost is the sum of the 3 amplitude error values.\n", " acost = amp_errvals.sum()\n", " # amplitude should not be shrinking or growing\n", " # compute the 3 amplitudes for the first half of the simulation\n", " # and then for the second\n", " half = sol2.y.shape[1]//2\n", " amps1 = sol2.y[:,:half].max(axis=1)-sol2.y[:,:half].min(axis=1)\n", " amps2 = sol2.y[:,half:].max(axis=1)-sol2.y[:,half:].min(axis=1)\n", " if debug:\n", " print( amps1 )\n", " print(amps2)\n", " # sum the amplitudes and take their difference. Then multiply\n", " # the difference by 100 (to amplify it) and square the result.\n", " bcost = (100*(amps1.sum()-amps2.sum()))**2\n", " if debug:\n", " print((period-24)**2, acost, bcost, amps1.sum(), amps2.sum() )\n", "\n", " # The total cost is the squared difference between the period and 24 hours \n", " # plue the two amplitude costs. We weight the acost with 5 because it values\n", " # tend to be small.\n", " return (period-24)**2 + 5*acost + bcost" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the GA for the Goodwin Oscillator\n", "\n", "Run the GA once for a given selection operator. Then plot the costs for each generation.\n", "\n", "With 50 parents and 250 kids per generation, this takes Stephanie about 1 minute to run." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# For reference, here is a set of working parameters.\n", "# I am going to construct the lower and upper bounds of parameter\n", "# space so that they capture these values.\n", "# 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", "param_names = ['alpha_1','K','n','gamma_1','alpha_2','gamma_2','alpha _3','gamma_3']\n", "lb = np.array([0.001,0.001,4,0.001,0.05,0.001,0.05,0.001])\n", "ub = np.array( [2,2,15,1,2,1,2,1])\n", "\n", "# We will use more parents and children than with the toy problem.\n", "settings = GASettings(numParents=50, numChildren=250)\n", "settings.numGenerations = 5;\n", "settings.eliteCount = 1;\n", "settings.selection = 'linearRanking';\n", "settings.linearRankingRatio = 8;\n", "\n", "# Run algorithm\n", "(params, cost, all_params, all_costs) = GA(goodwin_circadian_error, lb, ub, settings);\n", "\n", "# Plot costs\n", "legs = []\n", "for i in range(len(all_costs)):\n", " plt.plot( np.arange(all_costs[i].size), all_costs[i], '*' )\n", " legs.append( \"Generation\" + str(i) )\n", "\n", "plt.legend( legs ); " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Re-Plot costs using a log-scale for the y axis.\n", "# This allows us to see the differences a little bit better.\n", "legs = []\n", "for i in range(len(all_costs)):\n", " plt.semilogy( np.arange(all_costs[i].size), all_costs[i], '*' )\n", " legs.append( \"Generation\" + str(i) )\n", "plt.legend( legs ); " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example analysis of GA with Goodwin Oscillator\n", "\n", "Run the GA with linear ranking selection with a ratio of 8. Run the algorithm 5 times, saving the results in .npz files, because this will take awhile.\n", "\n", "Subsequent cells read in the files and plot the results." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# For reference, here is a set of working parameters.\n", "# I am going to construct the lower and upper bounds of parameter\n", "# space so that they capture these values.\n", "# 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", "param_names = ['alpha_1','K','n','gamma_1','alpha_2','gamma_2','alpha _3','gamma_3']\n", "lb = np.array([0.001,0.001,4,0.001,0.05,0.001,0.05,0.001])\n", "ub = np.array( [2,2,15,1,2,1,2,1])\n", "\n", "# We will use more parents and children than with the toy problem.\n", "settings = GASettings(numParents=100, numChildren=500)\n", "settings.numGenerations = 5;\n", "settings.eliteCount = 1;\n", "settings.selection = 'linearRanking';\n", "settings.linearRankingRatio = 8;\n", "\n", "for runNum in range(1,6):\n", " fn = \"Goodwin_linearRanking%d_run%d.npz\" % (settings.linearRankingRatio,runNum)\n", " if os.path.isfile( fn ):\n", " continue\n", "\n", " (params, cost, all_params, all_costs) = GA(goodwin_circadian_error, lb, ub, settings);\n", " np.savez( fn, params=params, cost=cost, all_params=all_params, all_costs=all_costs )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the costs for each generation for all 5 runs of a GA\n", "\n", "For a given tournament size, read in the files with GA results and plot the costs for each generation in a subplot." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lr_Ratio = 8\n", "\n", "plt.figure( figsize=(8,8) )\n", "for runNum in range(1,6):\n", " plt.subplot( 5, 2, (runNum-1)*2+1 )\n", " fn = \"Goodwin_linearRanking%d_run%d.npz\" % (lr_Ratio,runNum)\n", " f = np.load( fn, allow_pickle=True )\n", " all_costs = f['all_costs']\n", " all_params = f['all_params']\n", " legs = []\n", " for i in range(len(all_costs)):\n", " plt.plot( np.arange(all_costs[i].size), all_costs[i], '*' )\n", " legs.append( \"Generation\" + str(i) )\n", " plt.ylabel( \"Cost\")\n", " plt.xlabel( 'Individual')\n", "\n", " if runNum == 1:\n", " plt.legend( legs ); \n", " plt.subplot( 5, 2, (runNum-1)*2+2 )\n", " for i in range(len(all_costs)):\n", " plt.semilogy( np.arange(all_costs[i].size), all_costs[i], '*' )\n", " plt.xlabel( 'Individual')\n", "\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot results from a single GA run\n", "\n", "Load in the run, then plot the costs for each generation in one figure. Then, plot the values of the parameters in a jitter plot. Construct the jitter plot so that we can adjust the a threshold so that we can include parameters from all children or just those from children with costs below the threshold." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "runNum = 2\n", "tournamentSize = 5\n", "fn = \"Goodwin_tournament%d_run%d.npz\" %(tournamentSize,runNum)\n", "f = np.load( fn, allow_pickle=True )\n", "all_costs = f['all_costs']\n", "all_params = f['all_params']\n", "params = f['params']\n", "cost = f['cost']\n", "\n", "plt.figure( figsize=(12,4) )\n", "plt.subplot( 1, 2, 1 )\n", "legs = []\n", "for i in range(len(all_costs)):\n", " plt.plot( np.arange(all_costs[i].size), all_costs[i], '*' )\n", " legs.append( \"Generation\" + str(i) )\n", "plt.legend( legs );\n", "plt.ylabel(\"Cost\")\n", "plt.xlabel(\"Individual\")\n", "plt.subplot( 1, 2, 2 )\n", "for i in range(len(all_costs)):\n", " plt.semilogy( np.arange(all_costs[i].size), all_costs[i], '*' )\n", "plt.xlabel(\"Individual\")\n", "\n", "plt.figure(figsize=(12,4));\n", "PN = lb.size\n", "dp = 1/(2*settings.numGenerations)\n", "cost_threshold = np.inf #5\n", "\n", "for g in range(settings.numGenerations+1):\n", " idx = all_costs[g] < cost_threshold\n", " for i in range(lb.size):\n", " plt.plot( i+g*dp + np.random.rand(np.count_nonzero(idx))*dp, all_params[g][idx,i], '.C'+ str(g))\n", "plt.xticks( ticks=np.arange(0,PN)+0.25, labels=param_names);\n", "plt.ylabel( \"Parameter Value\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the simulation from the best individual in the GA run we have been studying." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "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", "\n", "plt.figure()\n", "plt.plot( sol2.t, sol2.y.T );\n", "plt.xlabel( \"Time (h)\")\n", "plt.ylabel( \"Concentration\")\n", "plt.legend( ['X','Y','Z'])\n", "\n", "for i in range(len(params)):\n", " print(param_names[i], ':', params[i])\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 }