geonomics_biogeo_2020.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biogeography Notebook 3\n",
    "\n",
    "The goal of this notebook is to get you to think critically about speciation and extinction. This notebook will also introduce you to the usefulness of modeling and simulations for understanding evolutionary scenarios. We are using a Python package call **\"Geonomics\"** that was conceived and created by Drew Hart, a PhD candidate at UC Berkeley in Ian Wang's lab in the College of Natural Resources, Environmental Science, Policy and Management Department.\n",
    "\n",
    "* **Geographic variation**: There is often a geographic component to genetic divergence, because both genetic drift and natural selection are facilitated, and gene flow is impeded, by geographic isolation. Thus, it makes sense to simulate evolution in a programatic way that allows for explicit spatial landscape scenarios. \n",
    "\n",
    "*Geonomics* allows the user to test hypotheses and run simulations about populations of organisms while incorporting any realistic or modeled geospatial **layers** of a landscape. The so-called forces of evolution can all be simulated: mutation, natural selection, genetic drift, gene-flow. The result is a very flexible platform to allow prediction of evolutionary outcomes during population divergence leading to speciation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models and simulations ##\n",
    "*A model is an abstract representation of a system or process.*\n",
    "* Models are useful because they allow us to precisely define the problem, articulate the relevant concepts, and then provide a means of analyzing data and communicating results.\n",
    "* Models allow us to predict the logical outcomes of how we think a system works and then explore the suite of conditions that vary in time and space.\n",
    "* Because knowledge is always incomplete, and all data needed to build a model are never available, all models require assumptions to “fill in the blanks.”\n",
    "* Most models are used to explore the consequences of our assumptions and hypotheses rather than to represent system structure and dynamics definitively.\n",
    "* When researchers are faced with answering questions in a large and complex system, it is difficult—sometimes impossible—to collect specimen data or conduct experiments at the ideal spatial and temporal scales. \n",
    "* The cost of experiments in time and money is often prohibitive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flowchart illustrating the major steps in building a model "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"img/model.PNG\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Glossary\n",
    "\n",
    "* **Speciation**: The process whereby two or more species evolve from a single ancestral population.\n",
    "* **Extinction**: The end of a species, although the processes that begin a species on the path to extinction typically occur at the scale of individual populations (i.e. local extinction). A geographic event, marking the final stage of the processes influencing geographic range collapse.\n",
    "* **Mutation**: Changes in DNA sequences of one or more nucleotides.\n",
    "* **Genetic drift**: The occurance of random changes in allele frequency between generations; a relatively weak force because it involves changes in the genetic constitution of a population caused solely by chance, genetic drift can be an important force in small, isolated populations.\n",
    "* **Neutral evolution**: Evolution that is due to neutral processes and *not due* to selection on traits \n",
    "* **Selection**: The change in a population that occurs because individuals express genetic traits that alter their interactions with their environment in such a way that their survival and reproductive success are enhanced relative to other individuals in the population.\n",
    "* **Gene flow**: The movement of alleles within or between populations because of the dispersal of gametes or offspring -- often tends to counter genetic drift and natural selection and impedes genetic divergence. However, importantly, some amount of gene flow, when populations are diverging (admixture) or after two populations have speciated (introgression) is thought to add novel genetic variation on which selection in the new environment can act, and thus actually can lead to increased genetic divergence. \n",
    "* **Carrying capacity**: The maximum population size of the species that the environment can sustain indefinitely, given the food, habitat, water, and other necessities available in the environment.\n",
    "* **Cost distance analysis**: The analysis of movement over continuous space, in which the cost of moving through any location is variable, so that some paths of movement have higher costs than others.\n",
    "* **Locus (plural loci)**: The position on a chromosome where a particular gene is located.\n",
    "* **Gene**: A specific sequence of nucleotides in DNA that is the functional unit of inheritance controlling the transmission and expression of one or more traits.\n",
    "* **Allele**: One of two or more alternative forms of a gene.\n",
    "* **Neutral landscape models (NLM)**: Random maps against which effects of processes that structure actual landscapes may be tested.\n",
    "* **Null model**: The hypothesis of no effect. A proper null model provides the required reference point against which alternatives may be contrasted. A null model may be a pattern in the absence of any mechanism (geological, ecological or otherwise) for which to compare alternative hypothesis patterns.\n",
    "* **Population**: A group of individuals belonging to the same species that live in the same region at the same time. \n",
    "* **Population density**: A measure of the number of organisms that make up a population in a defined area.\n",
    "* **Phenotype:** The physical appearance of an organism according to the individual's genetic makeup and environment. The expression of a particular trait, for example, skin color, height, behavior, biochemical characteristic, physiology, and so on.\n",
    "* **Genotype:** The genetic constitution of an individual organism. The set of genes in our DNA which is responsible for a particular trait.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Important notes: \n",
    "\n",
    "*Please don't be overwhelmed by the length of the code blocks in this notebook.* We intentionally left the code for you to see -- and exercises for you to modify the code -- so that you can think about the different parameters than need to be defined to construct a somewhat realistic evolutionary simulation.\n",
    "\n",
    "You need to read the **text cells** carefully and follow the instructions. Some **code cells** you will *run without changing*. Other times you will need to *edit the cells as instructed* to adjust the parameters of the model so we can simulate different scenarios. Reading the green `# comments` in Geonomics code will help to understand what you are actually doing during model set-up,  running of the model, and model results evaluation, and for **creating your own scenario** at the end. It is not necessary for you to read all of the `# comments` or code during this Notebook for biogeography, this is what a user of the Geonomics package will refer to while running their simulations.\n",
    "\n",
    "Recall that `def` means that we are defining a function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helpful Reminders: \n",
    "\n",
    "### Jupyter notebooks\n",
    "When a cell is running the circle in the top right of the notebook next to Python 3 will be solid grey. Some cells take several seconds to run. Just be patient.\n",
    "\n",
    "\n",
    "### Text cells\n",
    "In a notebook, each rectangle containing text or code is called a *cell*.\n",
    "\n",
    "Text cells (like this one) can be edited by double-clicking on them. They're written in a simple format called [Markdown](http://daringfireball.net/projects/markdown/syntax) to add formatting and section headings.  You don't need to learn Markdown, but you might want to.\n",
    "\n",
    "After you edit a text cell, click the \"run cell\" button at the top to confirm any changes. (Try not to delete the instructions of the lab.)\n",
    "\n",
    "**The only text cells that need to be modified are labeled \"YOUR RESPONSE HERE\" and are right below yellow question boxes. To edit a response, double click on YOUR RESPONSE HERE and type in your answer. Afterwards, run the cell with Shift-Enter.**\n",
    "\n",
    "### Code cells\n",
    "Other cells contain code in the Python 3 language. Running a code cell will execute all of the code it contains.\n",
    "\n",
    "To run the code in a code cell, first click on that cell to activate it.  It'll be highlighted with a little green or blue rectangle.  Next, either press the Run button or hold down the `shift` key and press `return` or `enter`.\n",
    "\n",
    "The only code cells that need to be modified are right below a blue exercise box.\n",
    "\n",
    "### Comments\n",
    "Comments are statements in English that the computer ignores. We use comments to explain what the surrounding code does. Comments appear in green after the `#` symbol."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h1><u>Evolutionary biogeography</u>:</h1></b>\n",
    "\n",
    "   <h3>using Geonomics</h3>\n",
    "\n",
    "\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "**Agenda**:\n",
    "<ul>\n",
    "    <li>1. short intro about Geonomics package</li>\n",
    "    <li>2. see briefly how Geonomics works (create parameters file, edit, make model, run model, visualize)</li>\n",
    "    <li>3. run a two-stage simulation of neutral evolution</li>\n",
    "    <li>4. run a simulation of natural selection to a heterogeneous environment</li>\n",
    "    <li>5. create and run your own model scenario</li>\n",
    "</ul>\n",
    "\n",
    "<p>\n",
    "\n",
    "</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>What does Geonomics do?</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geonomics is a new Python package that simulates the genetic evolution of species distributed across landscapes. The landscapes can be arbitrarily complex, and can undergo environmental change to over time. Species are composed of individuals, each of which has its own genome. The genes in those genomes can be used to determine the phenotypes of individuals, which can then undergo differential survival because of how poorly or well they are adapted to their local environments - that is, natural selection.\n",
    "\n",
    "The model operates by iterating through time steps (which you can think of being similar to generations). Each time step, a series of basic biological operations are carried out for all individuals. Over time, patterns of genetic diversity can build up on the landscape. These are what we are interested in studying, because these patterns in real-world species can help us make inferences about important things that have happened or are happening in nature, such as migration, adaptation, extinction, and the like.\n",
    "\n",
    "Below are a conceptual diagram and its caption, taken from [the Geonomics homepage](https://github.com/drewhart/geonomics). They give a simple explanation of how each time step works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"img/conceptual_diagram.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the center is a species on a multi-layer landscape that includes a selection layer (above) and a movement and carrying capacity layer (below). Surrounding the landscape is a flow-diagram of the major operations during a time step.\n",
    "\n",
    "- **movement**: During the movement stages (top-left), individuals move along movement vectors drawn from various distribution options.\n",
    "- **mating**: During the mating stage (top-right), an individual (purple outline) randomly chooses a mate (green outline) from all potential mates within its mating radius (dashed circle). The resulting offspring (dashed outline) disperses from its parents' midpoint along a randomly drawn dispersal vector.\n",
    "- **mortality**: During the mortality stage (bottom-right), deaths are modeled as a Bernoulli process, with the probability of mortality being a product of density-dependence and selection on all traits.\n",
    "- **changes**: During this stage (bottom-left), demographic change events (not pictured) and environmental change events (represented as a series of change rasters corresponding to scheduled time steps, t1, t2, …, tn), take place."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Okay! Now onto the code!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Import needed packages</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the following cell to bring in the necessary Python packages for our exercise today. Recall that you need to run this cell and any cells above the one you are working on each time you relaunch the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geonomics as gnx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# these are for submissions\n",
    "import otter\n",
    "from otter.export import export_notebook \n",
    "\n",
    "# and also set Matplotlib's default plotting style and plot size\n",
    "%matplotlib inline\n",
    "plt.rcParams[\"figure.figsize\"] = (9,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Define some helper functions</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will define some functions to create simple, geometric arrays, which we'll use as landscapes for our simulations.\n",
    "\n",
    "__No need to alter the code here.__ Just run it. We'll use the functions below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "# RUN THIS CODE WITHOUT CHANGING\n",
    "################################\n",
    "\n",
    "def make_unif_array(n):\n",
    "    \"\"\"Makes a square array of ones, size n x n cells.\"\"\"\n",
    "    array = np.ones((n,n))\n",
    "    return array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "# RUN THIS CODE WITHOUT CHANGING\n",
    "################################\n",
    "\n",
    "def add_vert_barrier(array, width = 5):\n",
    "    \"\"\"Returns a new array, which is the input array but with an added vertical barrier of the defined width.\"\"\"\n",
    "    if array.shape[0] % 2 != width % 2:\n",
    "        width += 1\n",
    "    assert width <= array.shape[0] - 2, (\"The width of the barrier should be \"\n",
    "                                       \"at least 2 less than the width of the landscape.\")\n",
    "    left_extent = int((array.shape[0] - width)/2)\n",
    "    right_extent = left_extent + width\n",
    "    new_array = np.copy(array)\n",
    "    new_array[:, left_extent:right_extent] = 0\n",
    "    return new_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "# RUN THIS CODE WITHOUT CHANGING\n",
    "################################\n",
    "\n",
    "def make_horz_grad_array(n, grad_width=8):\n",
    "    \"\"\"Makes a square array with a horizontal gradient\n",
    "    running from 0 to 1, size n x n cells.\"\"\"\n",
    "    assert grad_width <= n, \"Argument 'grad_width' must be <= argument 'n'.\"\n",
    "    if n % 2 != grad_width % 2:\n",
    "        grad_width += 1\n",
    "    grad_vals = np.linspace(0, 1, grad_width)\n",
    "    append_len = int((n - grad_width) / 2)\n",
    "    grad_vals = np.hstack([np.zeros((append_len)), grad_vals, np.ones((append_len))])\n",
    "    array = np.vstack([grad_vals.T]*n)\n",
    "    return array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also load some functions that will use the multivariate statistical method called Principal Component Analysis (PCA) to assign genetic-relatedness values to all the individuals in a species. \n",
    "\n",
    "The first function will plot all individuals in \"genetic relatedness space\", such that distances between individuals represent how genetically different they are. \n",
    "\n",
    "The second function will map individuals on the landscape, colored by their genetic relatedness, such that more similar-colored individuals are more related.\n",
    "\n",
    "We can use these visualizations to get a sense for how the genetic relatedness of the individuals in our species changes over time, because of distance, landscape features, or other influences.\n",
    "\n",
    "__No need to alter the code here.__ Just run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "# RUN THIS CODE WITHOUT CHANGING\n",
    "################################\n",
    "\n",
    "# function for running and plotting genetic PCA\n",
    "def plot_PCA(mod):\n",
    "    from copy import deepcopy\n",
    "    from sklearn.decomposition import PCA\n",
    "    figsize = 6\n",
    "    species = mod.comm[0]\n",
    "    # get array of resulting genomic data (i.e. 'speciome'),\n",
    "    # genotypes averaged by individual\n",
    "    speciome = np.mean(species._get_genotypes(), axis=2)\n",
    "    # run PCA on speciome\n",
    "    pca = PCA(n_components=2)\n",
    "    PCs = pca.fit_transform(speciome)\n",
    "    # normalize the PC results\n",
    "    norm_PCs = (PCs - np.min(PCs,\n",
    "                             axis=0)) / (np.max(PCs,\n",
    "                                                axis=0) - np.min(PCs,\n",
    "                                                                 axis=0))\n",
    "    # assign a value to each species, 0 or 1, indicating whether they're located on\n",
    "    # the left or right vertical half of the landscape\n",
    "    ind_colors = ['#00ffff' if ind.x < mod.land.dim[0]/2 else '#ff00ff' for ind in species.values()]\n",
    "    # plot individuals on PCs 1 and 2, colored by their landscape half\n",
    "    fig = plt.figure(figsize=(figsize, figsize), dpi= 80, facecolor='w', edgecolor='k')\n",
    "    plt.scatter(norm_PCs[:,0], norm_PCs[:, 1], color = ind_colors)\n",
    "    plt.xlabel('genetic PC 1')\n",
    "    plt.ylabel('genetic PC 2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "# RUN THIS CODE WITHOUT CHANGING\n",
    "################################\n",
    "\n",
    "# function for running and plotting genetic PCA\n",
    "def map_PCA(mod, lyr_num=0, mask=True):\n",
    "    from copy import deepcopy\n",
    "    from sklearn.decomposition import PCA\n",
    "    cmaps = {0: plt.cm.RdBu, 1: plt.cm.BrBG_r}\n",
    "    mark_size = 60\n",
    "    figsize = 8\n",
    "    species = mod.comm[0]\n",
    "    land = mod.land\n",
    "    # get array of resulting genomic data (i.e. 'speciome'),\n",
    "    # genotypes meaned by individual\n",
    "    speciome = np.mean(species._get_genotypes(), axis=2)\n",
    "    # run PCA on speciome\n",
    "    pca = PCA(n_components=3)\n",
    "    PCs = pca.fit_transform(speciome)\n",
    "    # normalize the PC results\n",
    "    norm_PCs = (PCs - np.min(PCs,\n",
    "                             axis=0)) / (np.max(PCs,\n",
    "                                                axis=0) - np.min(PCs,\n",
    "                                                                 axis=0))\n",
    "    # use first 3 PCs to get normalized values for R, G, & B colors\n",
    "    PC_colors = norm_PCs * 255\n",
    "    # scatter all individuals on top of landscape, colored by the\n",
    "    # RBG colors developed from the first 3 geonmic PCs\n",
    "    xs = mod.comm[0]._get_x()\n",
    "    ys = mod.comm[0]._get_y()\n",
    "    # get environmental raster, with barrier masked out\n",
    "    env = deepcopy(mod.land[lyr_num].rast)\n",
    "    if mask:\n",
    "        env[mod.land[lyr_num].rast == 0] = np.nan\n",
    "    # create light colormap for plotting landscape\n",
    "    # bot = plt.cm.get_cmap('Blues', 256)(np.linspace(0.4, 0.45, 2))[0]\n",
    "    # top = plt.cm.get_cmap('Reds', 256)(np.linspace(0.4, 0.45, 2))[0]\n",
    "    # cols = np.vstack((top, bot))\n",
    "    # cmap = mpl.colors.ListedColormap(cols, name='OrangeBlue')\n",
    "    cmap = cmaps[lyr_num]\n",
    "    cmap.set_bad(color='#8C8C8C')\n",
    "    # plot landscape\n",
    "    fig = plt.figure(figsize=(figsize, figsize), dpi= 80, facecolor='w', edgecolor='k')\n",
    "    # plt.imshow(masked_env, cmap=cmap, alpha=0.8)\n",
    "    plt.pcolormesh(land._x_cell_bds, land._y_cell_bds, env, cmap=cmap)\n",
    "    # scatter plot of individuals, colored by composite PC score\n",
    "    plt.scatter(xs, ys, c=PC_colors/255.0, s=mark_size, edgecolors='black')\n",
    "    # fix x and y limits\n",
    "    [f([dim - 0.5 for dim in (0, mod.land.dim[0])]) for f in (plt.xlim,\n",
    "                                                              plt.ylim)]\n",
    "    # get rid of x and y ticks\n",
    "    [f([]) for f in (plt.xticks, plt.yticks)]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>PART 1: Neutral evolution</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following section will walk through a neutral evolutionary simulation, on a very simple landscape.\n",
    "\n",
    "This simulation will highlight the effects of geographic and cost distance on genetic relatedness.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Set the parameter values </h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following long block of code creates a Geonomics parameters object.\n",
    "\n",
    "That object in turn will be used to instantiate (i.e. make a new, independent example of) a Geonomics model, which we can use to run simulations.\n",
    "\n",
    "The model will be composed of:\n",
    "\n",
    "- **a landscape**, with one **layer** that we'll use as both the environment and the carrying-capacity map; and\n",
    "- **a species**, which moves around on the landscape and has genomes containing only neutral loci."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# neutral_demo.py\n",
    "\n",
    "# This is a parameters file generated by Geonomics\n",
    "# (by the gnx.make_parameters_file() function).\n",
    "\n",
    "\n",
    "                   #   ::::::          :::    :: :::::::::::#\n",
    "             #::::::    ::::   :::      ::    :: :: ::::::::::: ::#\n",
    "          #:::::::::     ::            ::   :::::::::::::::::::::::::#\n",
    "        #::::::::::                      :::::::::: :::::: ::::::::  ::#\n",
    "      #  : ::::  ::                    ::::  : ::    :::::::: : ::  :    #\n",
    "     # GGGGG :EEEE: OOOOO   NN   NN   OOOOO   MM   MM IIIIII  CCCCC SSSSS #\n",
    "    # GG     EE    OO   OO  NNN  NN  OO   OO  MM   MM   II   CC     SS     #\n",
    "    # GG     EE   OO     OO NN N NN OO     OO MMM MMM   II   CC     SSSSSS #\n",
    "    # GG GGG EEEE OO     OO NN  NNN OO     OO MM M MM   II   CC         SS #\n",
    "    # GG   G EE    OO   OO  NN   NN  OO   OO  MM   MM   II   CC        SSS #\n",
    "     # GGGGG :EEEE: OOOOO   NN   NN   OOOOO   MM   MM IIIIII  CCCCC SSSSS #\n",
    "      #    : ::::::::               :::::::::: ::              ::  :   : #\n",
    "        #:    :::::                    :::::: :::             :::::::  #\n",
    "          #    :::                      :::::  ::              ::::: #\n",
    "             #  ::                      ::::                      #\n",
    "                   #                                        #\n",
    "                      #  :: ::    :::             #\n",
    "\n",
    "\n",
    "params = {\n",
    "###############################################################################\n",
    "\n",
    "###################\n",
    "#### LANDSCAPE ####\n",
    "###################\n",
    "    'landscape': {\n",
    "\n",
    "    ##############\n",
    "    #### main ####\n",
    "    ##############\n",
    "        'main': {\n",
    "            #x,y (a.k.a. j,i) dimensions of the Landscape\n",
    "            'dim':                      (20,20),\n",
    "            #x,y resolution of the Landscape\n",
    "            'res':                      (1,1),\n",
    "            #x,y coords of upper-left corner of the Landscape\n",
    "            'ulc':                      (0,0),\n",
    "            #projection of the Landscape\n",
    "            'prj':                      None,\n",
    "            }, # <END> 'main'\n",
    "\n",
    "    ################\n",
    "    #### layers ####\n",
    "    ################\n",
    "        'layers': {\n",
    "\n",
    "            #layer name (LAYER NAMES MUST BE UNIQUE!)\n",
    "            'lyr_0': {\n",
    "\n",
    "        #######################################\n",
    "        #### layer num. 0: init parameters ####\n",
    "        #######################################\n",
    "\n",
    "                #initiating parameters for this layer\n",
    "                'init': {\n",
    "\n",
    "                    #parameters for a 'defined'-type Layer\n",
    "                    'defined': {\n",
    "                        #raster to use for the Layer\n",
    "                        'rast':                   np.ones((20,20)),\n",
    "                        #point coordinates\n",
    "                        'pts':                    None,\n",
    "                        #point values\n",
    "                        'vals':                   None,\n",
    "                        #interpolation method {None, 'linear', 'cubic',\n",
    "                        #'nearest'}\n",
    "                        'interp_method':          None,\n",
    "\n",
    "                        }, # <END> 'defined'\n",
    "\n",
    "                    }, # <END> 'init'\n",
    "\n",
    "                }, # <END> layer num. 0\n",
    "\n",
    "\n",
    "\n",
    "    #### NOTE: Individual Layers' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Layers.\n",
    "\n",
    "\n",
    "            } # <END> 'layers'\n",
    "\n",
    "        }, # <END> 'landscape'\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "\n",
    "###################\n",
    "#### COMMUNITY ####\n",
    "###################\n",
    "    'comm': {\n",
    "\n",
    "        'species': {\n",
    "\n",
    "            #species name (SPECIES NAMES MUST BE UNIQUE!)\n",
    "            'spp_0': {\n",
    "\n",
    "            #####################################\n",
    "            #### spp num. 0: init parameters ####\n",
    "            #####################################\n",
    "\n",
    "                'init': {\n",
    "                    #starting number of individs\n",
    "                    'N':                250,\n",
    "                    #carrying-capacity Layer name\n",
    "                    'K_layer':          'lyr_0',\n",
    "                    #multiplicative factor for carrying-capacity layer\n",
    "                    'K_factor':         1,\n",
    "                    }, # <END> 'init'\n",
    "\n",
    "            #######################################\n",
    "            #### spp num. 0: mating parameters ####\n",
    "            #######################################\n",
    "\n",
    "                'mating'    : {\n",
    "                    #age(s) at sexual maturity (if tuple, female first)\n",
    "                    'repro_age':                0,\n",
    "                    #whether to assign sexes\n",
    "                    'sex':                      False,\n",
    "                    #ratio of males to females\n",
    "                    'sex_ratio':                1/1,\n",
    "                    #whether P(birth) should be weighted by parental dist\n",
    "                    'dist_weighted_birth':       False,\n",
    "                    #intrinsic growth rate\n",
    "                    'R':                        0.5,\n",
    "                    #intrinsic birth rate (MUST BE 0<=b<=1)\n",
    "                    'b':                        0.2,\n",
    "                    #expectation of distr of n offspring per mating pair\n",
    "                    'n_births_distr_lambda':    1,\n",
    "                    #whether n births should be fixed at n_births_dist_lambda\n",
    "                    'n_births_fixed':           True,\n",
    "                    #radius of mate-search area\n",
    "                    'mating_radius':            10,\n",
    "                    }, # <END> 'mating'\n",
    "\n",
    "            ##########################################\n",
    "            #### spp num. 0: mortality parameters ####\n",
    "            ##########################################\n",
    "\n",
    "                'mortality'     : {\n",
    "                    #maximum age\n",
    "                    'max_age':                      None,\n",
    "                    #min P(death) (MUST BE 0<=d_min<=1)\n",
    "                    'd_min':                        0,\n",
    "                    #max P(death) (MUST BE 0<=d_max<=1)\n",
    "                    'd_max':                        1,\n",
    "                    #width of window used to estimate local pop density\n",
    "                    'density_grid_window_width':    None,\n",
    "                    }, # <END> 'mortality'\n",
    "\n",
    "            #########################################\n",
    "            #### spp num. 0: movement parameters ####\n",
    "            #########################################\n",
    "\n",
    "                'movement': {\n",
    "                    #whether or not the species is mobile\n",
    "                    'move':                                 True,\n",
    "                    #mode of distr of movement direction\n",
    "                    'direction_distr_mu':                   0,\n",
    "                    #concentration of distr of movement direction\n",
    "                    'direction_distr_kappa':                0,\n",
    "                    #1st param of distr of movement distance\n",
    "                    'movement_distance_distr_param1':       0.5,\n",
    "                    #2nd param of distr of movement distance\n",
    "                    'movement_distance_distr_param2':       5e-8,\n",
    "                    #movement distance distr to use ('levy' or 'wald')\n",
    "                    'movement_distance_distr':              'levy',\n",
    "                    #1st param of distr of dispersal distance\n",
    "                    'dispersal_distance_distr_param1':      0,\n",
    "                    #2nd param of distr of dispersal distance\n",
    "                    'dispersal_distance_distr_param2':      5e-14,\n",
    "                    #dispersal distance distr to use ('levy' or 'wald')\n",
    "                    'dispersal_distance_distr':             'levy',\n",
    "                    'move_surf'     : {\n",
    "                        #move-surf Layer name\n",
    "                        'layer':                'lyr_0',\n",
    "                        #whether to use mixture distrs\n",
    "                        'mixture':              True,\n",
    "                        #concentration of distrs\n",
    "                        'vm_distr_kappa':       12,\n",
    "                        #length of approximation vectors for distrs\n",
    "                        'approx_len':           5000,\n",
    "                        }, # <END> 'move_surf'\n",
    "\n",
    "                    },    # <END> 'movement'\n",
    "\n",
    "\n",
    "            #####################################################\n",
    "            #### spp num. 0: genomic architecture parameters ####\n",
    "            #####################################################\n",
    "\n",
    "                'gen_arch': {\n",
    "                    #file defining custom genomic arch\n",
    "                    'gen_arch_file':            None,\n",
    "                    #num of loci\n",
    "                    'L':                        100,\n",
    "                    #starting allele frequency (None to draw freqs randomly)\n",
    "                    'start_p_fixed':            0.5,\n",
    "                    #whether to start neutral locus freqs at 0\n",
    "                    'start_neut_zero':          False,\n",
    "                    #genome-wide per-base neutral mut rate (0 to disable)\n",
    "                    'mu_neut':                  0,\n",
    "                    #genome-wide per-base deleterious mut rate (0 to disable)\n",
    "                    'mu_delet':                 0,\n",
    "                    #shape of distr of deleterious effect sizes\n",
    "                    'delet_alpha_distr_shape':  0.2,\n",
    "                    #scale of distr of deleterious effect sizes\n",
    "                    'delet_alpha_distr_scale':  0.2,\n",
    "                    #alpha of distr of recomb rates\n",
    "                    'r_distr_alpha':            None,\n",
    "                    #beta of distr of recomb rates\n",
    "                    'r_distr_beta':             None,\n",
    "                    #whether loci should be dominant (for allele '1')\n",
    "                    'dom':                      False,\n",
    "                    #whether to allow pleiotropy\n",
    "                    'pleiotropy':               False,\n",
    "                    #custom fn for drawing recomb rates\n",
    "                    'recomb_rate_custom_fn':    None,\n",
    "                    #number of recomb paths to hold in memory\n",
    "                    'n_recomb_paths_mem':       int(1e4),\n",
    "                    #total number of recomb paths to simulate\n",
    "                    'n_recomb_paths_tot':       int(1e5),\n",
    "                    #num of crossing-over events (i.e. recombs) to simulate\n",
    "                    'n_recomb_sims':            10_000,\n",
    "                    #whether to generate recombination paths at each timestep\n",
    "                    'allow_ad_hoc_recomb':      False,\n",
    "                    #whether to save mutation logs\n",
    "                    'mut_log':                  False,\n",
    "\n",
    "                    }, # <END> 'gen_arch'\n",
    "\n",
    "\n",
    "                }, # <END> spp num. 0\n",
    "\n",
    "\n",
    "\n",
    "    #### NOTE: individual Species' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Species.\n",
    "\n",
    "\n",
    "            }, # <END> 'species'\n",
    "\n",
    "        }, # <END> 'comm'\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "\n",
    "###############\n",
    "#### MODEL ####\n",
    "###############\n",
    "    'model': {\n",
    "        #total Model runtime (in timesteps)\n",
    "        'T':            100,\n",
    "        #min burn-in runtime (in timesteps)\n",
    "        'burn_T':       30,\n",
    "        #seed number\n",
    "        'num':          None,\n",
    "        #time step interval for simplication of tskit tables\n",
    "        'tskit_simp_interval':      100,\n",
    "\n",
    "\n",
    "        ###############################\n",
    "        #### iterations parameters ####\n",
    "        ###############################\n",
    "        'its': {\n",
    "            #num iterations\n",
    "            'n_its':            1,\n",
    "            #whether to randomize Landscape each iteration\n",
    "            'rand_landscape':   False,\n",
    "            #whether to randomize Community each iteration\n",
    "            'rand_comm':        False,\n",
    "            #whether to burn in each iteration\n",
    "            'repeat_burn':      False,\n",
    "            }, # <END> 'iterations'\n",
    "\n",
    "\n",
    "\n",
    "        } # <END> 'model'\n",
    "\n",
    "    } # <END> params\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Create a model from those parameters</h2></b> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make our params dict into a proper Geonomics ParamsDict object\n",
    "params = gnx.make_params_dict(params, 'neutral_demo')\n",
    "#then use it to make a model\n",
    "mod = gnx.make_model(parameters=params, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#take a look at the resulting object\n",
    "mod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll plot the model, to see how many individuals we have, and where on the landscape they're located.\n",
    "\n",
    "------------------------\n",
    "\n",
    "\n",
    "### NOTE:\n",
    "**This and all following figures will be interactive! You can use the buttons to the bottom left of the figure to zoom and move the figure as desired. The 'home' button will return the figure to its default view. This may later prove helpful when you are inspecting your figures and thinking about the associated questions!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the resulting object\n",
    "mod.plot(spp=0, lyr=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Concept check**: *What are the black dots on the gray background?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Run the model</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geonomics allows us to use our model in either of two ways:\n",
    "\n",
    "- **'run'**: Run the model to completion (i.e. for the number of iterations and the number of timesteps per iteration that are stipulated in the model's parameters). This is more convenient for leaving a model running and generating data.\n",
    "\n",
    "- **'walk'**: Run the model for any number of timesteps, then stop it. This is more convenient for running a model interactively and inspecting it. This is what we'll use.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#First we have to run the model burn-in\n",
    "\n",
    "#T=10000 allows the model run up to 10000 steps, but it will burn in much faster than that.\n",
    "\n",
    "#After it runs, scroll to the bottom to see how the time steps progressed.\n",
    "\n",
    "mod.walk(T=10000, mode='burn')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Concept check** *How many individuals did you start and end with? How many births and deaths were there at the beginning / and births and deaths during the last time step?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we can plot it again to see the burned-in stable state.\n",
    "mod.plot(0, 0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we can run the main phase of the model for any desired number of timesteps, then plot it again.\n",
    "mod.walk(10, 'main')\n",
    "mod.plot(0, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Concept check** *As you run the main phase of the model, and plot the result, what is changing?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Make the landscape heterogeneous</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What if we make the same model, but use a more interesting landscape?**\n",
    "\n",
    "Here we'll repeat the same process, but make our landscape using a **neutral landscape model (NLM)** algorithm (by changing some of the parameters in the 'lyr_0' section; compare that to the parameters code above). The NLM will produce a spatially autocorrelated pattern - i.e. one in which the landscape variable varies across space, but in a way that nearby values are more similar than distant values.\n",
    "\n",
    "Landscape ecologists often use NLMs as null models against which real-world patterns can be compared, to draw inferences about what might have generated real-world patterns. NLMs do not represent actual landscapes, but serve as the standard against which actual landscapes may be compared. This is referred to as a **null model**.\n",
    "\n",
    "*Here we won't use the NLM as a null model for hypothesis testing. Instead, we'll just use it as a quick way to get a more interesting landscape with a realistic pattern.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# neutral_demo_NLM.py\n",
    "\n",
    "# This is a parameters file generated by Geonomics\n",
    "# (by the gnx.make_parameters_file() function).\n",
    "\n",
    "\n",
    "                   #   ::::::          :::    :: :::::::::::#\n",
    "             #::::::    ::::   :::      ::    :: :: ::::::::::: ::#\n",
    "          #:::::::::     ::            ::   :::::::::::::::::::::::::#\n",
    "        #::::::::::                      :::::::::: :::::: ::::::::  ::#\n",
    "      #  : ::::  ::                    ::::  : ::    :::::::: : ::  :    #\n",
    "     # GGGGG :EEEE: OOOOO   NN   NN   OOOOO   MM   MM IIIIII  CCCCC SSSSS #\n",
    "    # GG     EE    OO   OO  NNN  NN  OO   OO  MM   MM   II   CC     SS     #\n",
    "    # GG     EE   OO     OO NN N NN OO     OO MMM MMM   II   CC     SSSSSS #\n",
    "    # GG GGG EEEE OO     OO NN  NNN OO     OO MM M MM   II   CC         SS #\n",
    "    # GG   G EE    OO   OO  NN   NN  OO   OO  MM   MM   II   CC        SSS #\n",
    "     # GGGGG :EEEE: OOOOO   NN   NN   OOOOO   MM   MM IIIIII  CCCCC SSSSS #\n",
    "      #    : ::::::::               :::::::::: ::              ::  :   : #\n",
    "        #:    :::::                    :::::: :::             :::::::  #\n",
    "          #    :::                      :::::  ::              ::::: #\n",
    "             #  ::                      ::::                      #\n",
    "                   #                                        #\n",
    "                      #  :: ::    :::             #\n",
    "\n",
    "\n",
    "params = {\n",
    "###############################################################################\n",
    "\n",
    "###################\n",
    "#### LANDSCAPE ####\n",
    "###################\n",
    "    'landscape': {\n",
    "\n",
    "    ##############\n",
    "    #### main ####\n",
    "    ##############\n",
    "        'main': {\n",
    "            #x,y (a.k.a. j,i) dimensions of the Landscape\n",
    "            'dim':                      (20,20),\n",
    "            #x,y resolution of the Landscape\n",
    "            'res':                      (1,1),\n",
    "            #x,y coords of upper-left corner of the Landscape\n",
    "            'ulc':                      (0,0),\n",
    "            #projection of the Landscape\n",
    "            'prj':                      None,\n",
    "            }, # <END> 'main'\n",
    "\n",
    "    ################\n",
    "    #### layers ####\n",
    "    ################\n",
    "        'layers': {\n",
    "\n",
    "            #layer name (LAYER NAMES MUST BE UNIQUE!)\n",
    "            'lyr_0': {\n",
    "\n",
    "        #######################################\n",
    "        #### layer num. 0: init parameters ####\n",
    "        #######################################\n",
    "\n",
    "                #initiating parameters for this layer\n",
    "                'init': {\n",
    "\n",
    "                    #parameters for an 'nlmpy'-type Layer\n",
    "                    'nlmpy': {\n",
    "                        #nlmpy function to use the create this Layer\n",
    "                        'function':                 'mpd',\n",
    "                        #number of rows (MUST EQUAL LANDSCAPE DIMENSION y!)\n",
    "                        'nRow':                     20,\n",
    "                        #number of cols (MUST EQUAL LANDSCAPE DIMENSION x!)\n",
    "                        'nCol':                     20,\n",
    "                        #level of spatial autocorrelation in element values\n",
    "                        'h':                        1,\n",
    "\n",
    "                        }, # <END> 'nlmpy'\n",
    "\n",
    "                    }, # <END> 'init'\n",
    "\n",
    "                }, # <END> layer num. 0\n",
    "\n",
    "\n",
    "\n",
    "    #### NOTE: Individual Layers' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Layers.\n",
    "\n",
    "\n",
    "            } # <END> 'layers'\n",
    "\n",
    "        }, # <END> 'landscape'\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "\n",
    "###################\n",
    "#### COMMUNITY ####\n",
    "###################\n",
    "    'comm': {\n",
    "\n",
    "        'species': {\n",
    "\n",
    "            #species name (SPECIES NAMES MUST BE UNIQUE!)\n",
    "            'spp_0': {\n",
    "\n",
    "            #####################################\n",
    "            #### spp num. 0: init parameters ####\n",
    "            #####################################\n",
    "\n",
    "                'init': {\n",
    "                    #starting number of individs\n",
    "                    'N':                250,\n",
    "                    #carrying-capacity Layer name\n",
    "                    'K_layer':          'lyr_0',\n",
    "                    #multiplicative factor for carrying-capacity layer\n",
    "                    'K_factor':         1,\n",
    "                    }, # <END> 'init'\n",
    "\n",
    "            #######################################\n",
    "            #### spp num. 0: mating parameters ####\n",
    "            #######################################\n",
    "\n",
    "                'mating'    : {\n",
    "                    #age(s) at sexual maturity (if tuple, female first)\n",
    "                    'repro_age':                0,\n",
    "                    #whether to assign sexes\n",
    "                    'sex':                      False,\n",
    "                    #ratio of males to females\n",
    "                    'sex_ratio':                1/1,\n",
    "                    #whether P(birth) should be weighted by parental dist\n",
    "                    'dist_weighted_birth':       False,\n",
    "                    #intrinsic growth rate\n",
    "                    'R':                        0.5,\n",
    "                    #intrinsic birth rate (MUST BE 0<=b<=1)\n",
    "                    'b':                        0.2,\n",
    "                    #expectation of distr of n offspring per mating pair\n",
    "                    'n_births_distr_lambda':    1,\n",
    "                    #whether n births should be fixed at n_births_dist_lambda\n",
    "                    'n_births_fixed':           True,\n",
    "                    #radius of mate-search area\n",
    "                    'mating_radius':            10,\n",
    "                    }, # <END> 'mating'\n",
    "\n",
    "            ##########################################\n",
    "            #### spp num. 0: mortality parameters ####\n",
    "            ##########################################\n",
    "\n",
    "                'mortality'     : {\n",
    "                    #maximum age\n",
    "                    'max_age':                      None,\n",
    "                    #min P(death) (MUST BE 0<=d_min<=1)\n",
    "                    'd_min':                        0,\n",
    "                    #max P(death) (MUST BE 0<=d_max<=1)\n",
    "                    'd_max':                        1,\n",
    "                    #width of window used to estimate local pop density\n",
    "                    'density_grid_window_width':    None,\n",
    "                    }, # <END> 'mortality'\n",
    "\n",
    "            #########################################\n",
    "            #### spp num. 0: movement parameters ####\n",
    "            #########################################\n",
    "\n",
    "                'movement': {\n",
    "                    #whether or not the species is mobile\n",
    "                    'move':                                 True,\n",
    "                    #mode of distr of movement direction\n",
    "                    'direction_distr_mu':                   0,\n",
    "                    #concentration of distr of movement direction\n",
    "                    'direction_distr_kappa':                0,\n",
    "                    #1st param of distr of movement distance\n",
    "                    'movement_distance_distr_param1':       0.5,\n",
    "                    #2nd param of distr of movement distance\n",
    "                    'movement_distance_distr_param2':       5e-8,\n",
    "                    #movement distance distr to use ('levy' or 'wald')\n",
    "                    'movement_distance_distr':              'levy',\n",
    "                    #1st param of distr of dispersal distance\n",
    "                    'dispersal_distance_distr_param1':      0,\n",
    "                    #2nd param of distr of dispersal distance\n",
    "                    'dispersal_distance_distr_param2':      5e-14,\n",
    "                    #dispersal distance distr to use ('levy' or 'wald')\n",
    "                    'dispersal_distance_distr':             'levy',\n",
    "                    'move_surf'     : {\n",
    "                        #move-surf Layer name\n",
    "                        'layer':                'lyr_0',\n",
    "                        #whether to use mixture distrs\n",
    "                        'mixture':              True,\n",
    "                        #concentration of distrs\n",
    "                        'vm_distr_kappa':       12,\n",
    "                        #length of approximation vectors for distrs\n",
    "                        'approx_len':           5000,\n",
    "                        }, # <END> 'move_surf'\n",
    "\n",
    "                    },    # <END> 'movement'\n",
    "\n",
    "\n",
    "            #####################################################\n",
    "            #### spp num. 0: genomic architecture parameters ####\n",
    "            #####################################################\n",
    "\n",
    "                'gen_arch': {\n",
    "                    #file defining custom genomic arch\n",
    "                    'gen_arch_file':            None,\n",
    "                    #num of loci\n",
    "                    'L':                        100,\n",
    "                    #starting allele frequency (None to draw freqs randomly)\n",
    "                    'start_p_fixed':            0.5,\n",
    "                    #whether to start neutral locus freqs at 0\n",
    "                    'start_neut_zero':          False,\n",
    "                    #genome-wide per-base neutral mut rate (0 to disable)\n",
    "                    'mu_neut':                  0,\n",
    "                    #genome-wide per-base deleterious mut rate (0 to disable)\n",
    "                    'mu_delet':                 0,\n",
    "                    #shape of distr of deleterious effect sizes\n",
    "                    'delet_alpha_distr_shape':  0.2,\n",
    "                    #scale of distr of deleterious effect sizes\n",
    "                    'delet_alpha_distr_scale':  0.2,\n",
    "                    #alpha of distr of recomb rates\n",
    "                    'r_distr_alpha':            None,\n",
    "                    #beta of distr of recomb rates\n",
    "                    'r_distr_beta':             None,\n",
    "                    #whether loci should be dominant (for allele '1')\n",
    "                    'dom':                      False,\n",
    "                    #whether to allow pleiotropy\n",
    "                    'pleiotropy':               False,\n",
    "                    #custom fn for drawing recomb rates\n",
    "                    'recomb_rate_custom_fn':    None,\n",
    "                    #number of recomb paths to hold in memory\n",
    "                    'n_recomb_paths_mem':       int(1e4),\n",
    "                    #total number of recomb paths to simulate\n",
    "                    'n_recomb_paths_tot':       int(1e5),\n",
    "                    #num of crossing-over events (i.e. recombs) to simulate\n",
    "                    'n_recomb_sims':            10_000,\n",
    "                    #whether to generate recombination paths at each timestep\n",
    "                    'allow_ad_hoc_recomb':      False,\n",
    "                    #whether to save mutation logs\n",
    "                    'mut_log':                  False,\n",
    "\n",
    "                    }, # <END> 'gen_arch'\n",
    "\n",
    "\n",
    "                }, # <END> spp num. 0\n",
    "\n",
    "\n",
    "\n",
    "    #### NOTE: individual Species' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Species.\n",
    "\n",
    "\n",
    "            }, # <END> 'species'\n",
    "\n",
    "        }, # <END> 'comm'\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "\n",
    "###############\n",
    "#### MODEL ####\n",
    "###############\n",
    "    'model': {\n",
    "        #total Model runtime (in timesteps)\n",
    "        'T':            100,\n",
    "        #min burn-in runtime (in timesteps)\n",
    "        'burn_T':       30,\n",
    "        #seed number\n",
    "        'num':          None,\n",
    "        #time step interval for simplication of tskit tables\n",
    "        'tskit_simp_interval':      100,\n",
    "\n",
    "\n",
    "        ###############################\n",
    "        #### iterations parameters ####\n",
    "        ###############################\n",
    "        'its': {\n",
    "            #num iterations\n",
    "            'n_its':            1,\n",
    "            #whether to randomize Landscape each iteration\n",
    "            'rand_landscape':   False,\n",
    "            #whether to randomize Community each iteration\n",
    "            'rand_comm':        False,\n",
    "            #whether to burn in each iteration\n",
    "            'repeat_burn':      False,\n",
    "            }, # <END> 'iterations'\n",
    "\n",
    "\n",
    "\n",
    "        } # <END> 'model'\n",
    "\n",
    "    } # <END> params\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll do the same steps as before:\n",
    "- make the model\n",
    "- run it\n",
    "- plot the results\n",
    "\n",
    "To do that, go ahead and **run the next few cells**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = gnx.make_params_dict(params, 'heterogeneous_demo')\n",
    "mod = gnx.make_model(parameters=params, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.walk(T=10000, mode='burn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.walk(T=10, mode='main')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.plot(lyr=0, spp=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- BEGIN QUESTION -->\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>QUESTION 1: Briefly explain how Geonomics works and compare the models we built and ran (one without spatial heterogeneity and one with spatial heterogenetity using the NLM as our layer).</b>\n",
    "    <br />    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Q1: YOUR RESPONSE HERE\n",
    "\n",
    "<!-- END QUESTION -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  **Hint**: Recall that in the second model we used an NLM as the habitat map and the carrying-capacity map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>PART 2: Environmental Change</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've seen a quick example of how Geonomics works, we'll make a more complicated model, then use the results to explore questions about gene flow, isolation, and genetic divergence in species distributed across space.\n",
    "\n",
    "For example, **what would happen if a barrier were to arise on the landscape?** Let's simulate one to find out!\n",
    "\n",
    "This time, for the first 500 time steps, our species will evolve on a homogeneous landscape (just like our first model). We will make some plots to observe the model, we'll create another model with a homogeneous landscape layer, then run that model for 500 timesteps and observe how genetic relatedness changes over time. *This will be equivalent to watching genetic evolution in the first model we made.* \n",
    "\n",
    "But then, a big central barrier will arise. We will see this barrier arise right after it arises, simulate another 500 time steps, then use the same plots to see how the barrier influenced genetic relatedness.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start off, we need to create and save an array for our landscape layer. To do that, we'll use one of the helper functions that we loaded at the beginning. To do that, **run the next cell**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save a new file\n",
    "barr_lyr = add_vert_barrier(make_unif_array(40))\n",
    "plt.imshow(barr_lyr, cmap = 'RdBu_r')\n",
    "np.savetxt('barrier.txt', barr_lyr, fmt='%0.5f')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll create a new model, using extra parameters that set up a landscape change event that uses that file. \n",
    "\n",
    "Here are the key changes we've made in order to characterize that landscape change event:\n",
    "\n",
    "- **dim** parameter: *We set the landscape dimensions to __(40, 40)__.*\n",
    "\n",
    "- **rast** parameter: *We set the landscape raster to a 40-cell by 40-cell uniform array of ones. We do this using a function we loaded at the beginning: type __`make_unif_array(40)`__.*\n",
    "\n",
    "- **change_rast** parameter: *We set the value of this parameter, in the landscape-change section, to __'./barrier.txt'__ (the barrier layer that we just created)*.\n",
    "\n",
    "- **start_t** parameter: *We set the landscape-change event to start at timestep __500__*.\n",
    "\n",
    "- **end_t** parameter: *We set the landscape-change event to end at timestep __600__*.\n",
    "\n",
    "- **n_steps** parameter: *To keep things simple, we set the number of stepwise-changes in the landscape-change event to just __2__*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# barrier_demo.py\n",
    "\n",
    "# This is a parameters file generated by Geonomics\n",
    "# (by the gnx.make_parameters_file() function).\n",
    "\n",
    "\n",
    "                   #   ::::::          :::    :: :::::::::::#\n",
    "             #::::::    ::::   :::      ::    :: :: ::::::::::: ::#\n",
    "          #:::::::::     ::            ::   :::::::::::::::::::::::::#\n",
    "        #::::::::::                      :::::::::: :::::: ::::::::  ::#\n",
    "      #  : ::::  ::                    ::::  : ::    :::::::: : ::  :    #\n",
    "     # GGGGG :EEEE: OOOOO   NN   NN   OOOOO   MM   MM IIIIII  CCCCC SSSSS #\n",
    "    # GG     EE    OO   OO  NNN  NN  OO   OO  MM   MM   II   CC     SS     #\n",
    "    # GG     EE   OO     OO NN N NN OO     OO MMM MMM   II   CC     SSSSSS #\n",
    "    # GG GGG EEEE OO     OO NN  NNN OO     OO MM M MM   II   CC         SS #\n",
    "    # GG   G EE    OO   OO  NN   NN  OO   OO  MM   MM   II   CC        SSS #\n",
    "     # GGGGG :EEEE: OOOOO   NN   NN   OOOOO   MM   MM IIIIII  CCCCC SSSSS #\n",
    "      #    : ::::::::               :::::::::: ::              ::  :   : #\n",
    "        #:    :::::                    :::::: :::             :::::::  #\n",
    "          #    :::                      :::::  ::              ::::: #\n",
    "             #  ::                      ::::                      #\n",
    "                   #                                        #\n",
    "                      #  :: ::    :::             #\n",
    "\n",
    "\n",
    "params = {\n",
    "###############################################################################\n",
    "\n",
    "###################\n",
    "#### LANDSCAPE ####\n",
    "###################\n",
    "    'landscape': {\n",
    "\n",
    "    ##############\n",
    "    #### main ####\n",
    "    ##############\n",
    "        'main': {\n",
    "            #x,y (a.k.a. j,i) dimensions of the Landscape\n",
    "            'dim':                      (40,40),\n",
    "            #x,y resolution of the Landscape\n",
    "            'res':                      (1,1),\n",
    "            #x,y coords of upper-left corner of the Landscape\n",
    "            'ulc':                      (0,0),\n",
    "            #projection of the Landscape\n",
    "            'prj':                      None,\n",
    "            }, # <END> 'main'\n",
    "\n",
    "    ################\n",
    "    #### layers ####\n",
    "    ################\n",
    "        'layers': {\n",
    "\n",
    "            #layer name (LAYER NAMES MUST BE UNIQUE!)\n",
    "            'lyr_0': {\n",
    "\n",
    "        #######################################\n",
    "        #### layer num. 0: init parameters ####\n",
    "        #######################################\n",
    "\n",
    "                #initiating parameters for this layer\n",
    "                'init': {\n",
    "\n",
    "                    #parameters for a 'defined'-type Layer\n",
    "                    'defined': {\n",
    "                        #raster to use for the Layer\n",
    "                        'rast':                   make_unif_array(40),\n",
    "                        #point coordinates\n",
    "                        'pts':                    None,\n",
    "                        #point values\n",
    "                        'vals':                   None,\n",
    "                        #interpolation method {None, 'linear', 'cubic',\n",
    "                        #'nearest'}\n",
    "                        'interp_method':          None,\n",
    "\n",
    "                        }, # <END> 'defined'\n",
    "\n",
    "                    }, # <END> 'init'\n",
    "\n",
    "            #########################################\n",
    "            #### layer num. 0: change parameters ####\n",
    "            #########################################\n",
    "\n",
    "                #landscape-change events for this Layer\n",
    "                'change': {\n",
    "\n",
    "                    0: {\n",
    "                        #array or file for final raster of event, or directory\n",
    "                        #of files for each stepwise change in event\n",
    "                        'change_rast':      './barrier.txt',\n",
    "                        #starting timestep of event\n",
    "                        'start_t':          500,\n",
    "                        #ending timestep of event\n",
    "                        'end_t':            600,\n",
    "                        #number of stepwise changes in event\n",
    "                        'n_steps':          1,\n",
    "                        }, # <END> event 0\n",
    "\n",
    "                    }, # <END> 'change'\n",
    "\n",
    "                }, # <END> layer num. 0\n",
    "\n",
    "\n",
    "\n",
    "    #### NOTE: Individual Layers' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Layers.\n",
    "\n",
    "\n",
    "            } # <END> 'layers'\n",
    "\n",
    "        }, # <END> 'landscape'\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "\n",
    "###################\n",
    "#### COMMUNITY ####\n",
    "###################\n",
    "    'comm': {\n",
    "\n",
    "        'species': {\n",
    "\n",
    "            #species name (SPECIES NAMES MUST BE UNIQUE!)\n",
    "            'spp_0': {\n",
    "\n",
    "            #####################################\n",
    "            #### spp num. 0: init parameters ####\n",
    "            #####################################\n",
    "\n",
    "                'init': {\n",
    "                    #starting number of individs\n",
    "                    'N':                250,\n",
    "                    #carrying-capacity Layer name\n",
    "                    'K_layer':          'lyr_0',\n",
    "                    #multiplicative factor for carrying-capacity layer\n",
    "                    'K_factor':         1,\n",
    "                    }, # <END> 'init'\n",
    "\n",
    "            #######################################\n",
    "            #### spp num. 0: mating parameters ####\n",
    "            #######################################\n",
    "\n",
    "                'mating'    : {\n",
    "                    #age(s) at sexual maturity (if tuple, female first)\n",
    "                    'repro_age':                0,\n",
    "                    #whether to assign sexes\n",
    "                    'sex':                      False,\n",
    "                    #ratio of males to females\n",
    "                    'sex_ratio':                1/1,\n",
    "                    #whether P(birth) should be weighted by parental dist\n",
    "                    'dist_weighted_birth':       False,\n",
    "                    #intrinsic growth rate\n",
    "                    'R':                        0.5,\n",
    "                    #intrinsic birth rate (MUST BE 0<=b<=1)\n",
    "                    'b':                        0.2,\n",
    "                    #expectation of distr of n offspring per mating pair\n",
    "                    'n_births_distr_lambda':    1,\n",
    "                    #whether n births should be fixed at n_births_dist_lambda\n",
    "                    'n_births_fixed':           True,\n",
    "                    #radius of mate-search area\n",
    "                    'mating_radius':            10,\n",
    "                    }, # <END> 'mating'\n",
    "\n",
    "            ##########################################\n",
    "            #### spp num. 0: mortality parameters ####\n",
    "            ##########################################\n",
    "\n",
    "                'mortality'     : {\n",
    "                    #maximum age\n",
    "                    'max_age':                      None,\n",
    "                    #min P(death) (MUST BE 0<=d_min<=1)\n",
    "                    'd_min':                        0,\n",
    "                    #max P(death) (MUST BE 0<=d_max<=1)\n",
    "                    'd_max':                        1,\n",
    "                    #width of window used to estimate local pop density\n",
    "                    'density_grid_window_width':    None,\n",
    "                    }, # <END> 'mortality'\n",
    "\n",
    "            #########################################\n",
    "            #### spp num. 0: movement parameters ####\n",
    "            #########################################\n",
    "\n",
    "                'movement': {\n",
    "                    #whether or not the species is mobile\n",
    "                    'move':                                 True,\n",
    "                    #mode of distr of movement direction\n",
    "                    'direction_distr_mu':                   0,\n",
    "                    #concentration of distr of movement direction\n",
    "                    'direction_distr_kappa':                0,\n",
    "                    #1st param of distr of movement distance\n",
    "                    'movement_distance_distr_param1':       0.5,\n",
    "                    #2nd param of distr of movement distance\n",
    "                    'movement_distance_distr_param2':       5e-8,\n",
    "                    #movement distance distr to use ('levy' or 'wald')\n",
    "                    'movement_distance_distr':              'levy',\n",
    "                    #1st param of distr of dispersal distance\n",
    "                    'dispersal_distance_distr_param1':      0,\n",
    "                    #2nd param of distr of dispersal distance\n",
    "                    'dispersal_distance_distr_param2':      5e-14,\n",
    "                    #dispersal distance distr to use ('levy' or 'wald')\n",
    "                    'dispersal_distance_distr':             'levy',\n",
    "                    'move_surf'     : {\n",
    "                        #move-surf Layer name\n",
    "                        'layer':                'lyr_0',\n",
    "                        #whether to use mixture distrs\n",
    "                        'mixture':              True,\n",
    "                        #concentration of distrs\n",
    "                        'vm_distr_kappa':       12,\n",
    "                        #length of approximation vectors for distrs\n",
    "                        'approx_len':           5000,\n",
    "                        }, # <END> 'move_surf'\n",
    "\n",
    "                    },    # <END> 'movement'\n",
    "\n",
    "\n",
    "            #####################################################\n",
    "            #### spp num. 0: genomic architecture parameters ####\n",
    "            #####################################################\n",
    "\n",
    "                'gen_arch': {\n",
    "                    #file defining custom genomic arch\n",
    "                    'gen_arch_file':            None,\n",
    "                    #num of loci\n",
    "                    'L':                        100,\n",
    "                    #starting allele frequency (None to draw freqs randomly)\n",
    "                    'start_p_fixed':            0.5,\n",
    "                    #whether to start neutral locus freqs at 0\n",
    "                    'start_neut_zero':          False,\n",
    "                    #genome-wide per-base neutral mut rate (0 to disable)\n",
    "                    'mu_neut':                  0,\n",
    "                    #genome-wide per-base deleterious mut rate (0 to disable)\n",
    "                    'mu_delet':                 0,\n",
    "                    #shape of distr of deleterious effect sizes\n",
    "                    'delet_alpha_distr_shape':  0.2,\n",
    "                    #scale of distr of deleterious effect sizes\n",
    "                    'delet_alpha_distr_scale':  0.2,\n",
    "                    #alpha of distr of recomb rates\n",
    "                    'r_distr_alpha':            None,\n",
    "                    #beta of distr of recomb rates\n",
    "                    'r_distr_beta':             None,\n",
    "                    #whether loci should be dominant (for allele '1')\n",
    "                    'dom':                      False,\n",
    "                    #whether to allow pleiotropy\n",
    "                    'pleiotropy':               False,\n",
    "                    #custom fn for drawing recomb rates\n",
    "                    'recomb_rate_custom_fn':    None,\n",
    "                    #number of recomb paths to hold in memory\n",
    "                    'n_recomb_paths_mem':       int(1e4),\n",
    "                    #total number of recomb paths to simulate\n",
    "                    'n_recomb_paths_tot':       int(1e5),\n",
    "                    #num of crossing-over events (i.e. recombs) to simulate\n",
    "                    'n_recomb_sims':            10_000,\n",
    "                    #whether to generate recombination paths at each timestep\n",
    "                    'allow_ad_hoc_recomb':      False,\n",
    "                    #whether to save mutation logs\n",
    "                    'mut_log':                  False,\n",
    "\n",
    "                    }, # <END> 'gen_arch'\n",
    "\n",
    "\n",
    "                }, # <END> spp num. 0\n",
    "\n",
    "\n",
    "\n",
    "    #### NOTE: individual Species' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Species.\n",
    "\n",
    "\n",
    "            }, # <END> 'species'\n",
    "\n",
    "        }, # <END> 'comm'\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "\n",
    "###############\n",
    "#### MODEL ####\n",
    "###############\n",
    "    'model': {\n",
    "        #total Model runtime (in timesteps)\n",
    "        'T':            100,\n",
    "        #min burn-in runtime (in timesteps)\n",
    "        'burn_T':       30,\n",
    "        #seed number\n",
    "        'num':          None,\n",
    "        #time step interval for simplication of tskit tables\n",
    "        'tskit_simp_interval':      100,\n",
    "\n",
    "\n",
    "        ###############################\n",
    "        #### iterations parameters ####\n",
    "        ###############################\n",
    "        'its': {\n",
    "            #num iterations\n",
    "            'n_its':            1,\n",
    "            #whether to randomize Landscape each iteration\n",
    "            'rand_landscape':   False,\n",
    "            #whether to randomize Community each iteration\n",
    "            'rand_comm':        False,\n",
    "            #whether to burn in each iteration\n",
    "            'repeat_burn':      False,\n",
    "            }, # <END> 'iterations'\n",
    "\n",
    "\n",
    "\n",
    "        } # <END> 'model'\n",
    "\n",
    "    } # <END> params\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, just like before, we'll create the model, burn it in, then run it and observe.\n",
    "\n",
    "We'll start by just running the portion of the model **before** the landscape change event, to observe what happens during that time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make our params dict into a proper Geonomics ParamsDict object\n",
    "params = gnx.make_params_dict(params, 'barrier_demo')\n",
    "#then use it to make a model\n",
    "mod = gnx.make_model(parameters=params, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#First we have to run the model burn-in, as always\n",
    "\n",
    "#SCROLL DOWN TO OBSERVE BURN-IN AND SEE PLOT\n",
    "\n",
    "mod.walk(T=10000, mode='burn')\n",
    "#plot the resulting object\n",
    "mod.plot(spp=0, lyr=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Now we can run the main phase of the model for the first 500 timesteps.\n",
    "mod.walk(500, 'main')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then let's plot the individuals in \"genetic-relatedness space\",\n",
    "and also map them in geographic space but colored by their genetic relatedness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We'll use the genetic PCA functions we defined at the beginning.\n",
    "plot_PCA(mod)\n",
    "map_PCA(mod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- BEGIN QUESTION -->\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>QUESTION 2: What can you say about gene flow among individuals of this species from the plot of genetic relatedness? Is there any evidence of genetic drift happening?</b>\n",
    "    <br />    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Q2: YOUR RESPONSE HERE\n",
    "\n",
    "<!-- END QUESTION -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>The effect of the barrier</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will continue running this model. Over the next 50 timesteps, our landscape-change event will occur. The barrier will \"arise\" on the landscape, and will gradually reduce individuals' ability to move across it.\n",
    "\n",
    "We will plot our same genetic relatedness plots a few more times, so that we can see if the barrier has any effect on relatedness, and if so, how that effect unfolds over time.\n",
    "\n",
    "First, we'll run the model 10 more timesteps, the see the landscape right after the barrier has arisen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.walk(10, 'main')\n",
    "mod.plot(spp=0, lyr=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_PCA(mod)\n",
    "map_PCA(mod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we'll run the model for another 500 time steps, to give any effects ample time to play out. Then we'll plot genetic relatedness one last time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then let's plot genetic relatedness again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mod.walk(500, 'main')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_PCA(mod)\n",
    "map_PCA(mod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Geonomics* also has some built in functionality for visualizing the spatial pedigree of the current population. In other words, for any individual currently 'alive' in the model, *Geonomics* can serially map the birth locations of that individual's ancestors, going backward through generational time, all the way to the birth location of the ancestor who was alive when the model began. (These 'ancestry tracks' are plotted as confetti-colored lines, a different line for each currently living individual who was plotted. The tracks get thinner as they go further back in time toward the location of the original ancestor in the simulation.)\n",
    "\n",
    "This may be a bit complicated, but if we inspect these visualizations carefully then we can use them to help us undertsand gene flow in our model -- that is, how genes have moved around our landscape during the simulation.\n",
    "\n",
    "Here is a pair of such visualizations, one for each of two arbitrary 'loci' (i.e. genetic positions; 'locus' is singular, 'loci' is plural) in our simulated genomes. (To avoid making a crazy-looking plot, only a random subset of 10 individuals has their gene flow plotted.) **Observe these plots closely. What general patterns stand out to you? How would you summarize what those patterns suggest about the dynamics of gene flow and their relationship to landscape structure?**\n",
    "\n",
    "(**NOTE: We've made these plots interactive, so that you can zoom in closely if need be, to explore individual ancestry tracks.**)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax1.set_title('Locus 5')\n",
    "spp = mod.comm[0]\n",
    "spp._plot_gene_flow(5, 'lineage', mod.land, individs=spp._get_random_individuals(10))\n",
    "ax2 = fig.add_subplot(122)\n",
    "ax2.set_title('Locus 17')\n",
    "spp._plot_gene_flow(17, 'lineage', mod.land, individs=spp._get_random_individuals(10))\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- BEGIN QUESTION -->\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>QUESTION 3: What was the influence of the rise of this barrier (a vicariance event) on gene flow and genetic drift? Be specific.</b>\n",
    "    <br />    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Q3: YOUR RESPONSE HERE\n",
    "\n",
    "<!-- END QUESTION -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Part 3: Natural Selection</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the remainder of today's lab, we'll explore a different simulation scenario, one in which natural selection is operating.\n",
    "\n",
    "Remember, **natural selection occurs because**:\n",
    "1. Individuals vary phenotypically.\n",
    "2. They vary phenotypically because they vary genotypically.\n",
    "3. Genotypic variations are heritable.\n",
    "4. Individuals' variation makes some more likely to survive and reproduce than others (i.e. more fit).\n",
    "5. More fit individuals pass their genes on more often.\n",
    "6. Over time, this changes the genetic and phenotypic constitution of a population, and thus of a species.\n",
    "\n",
    "Geonomics will simulate that process - in space! As you're watching these simulations, keep in mind\n",
    "the concepts you just explored in the previous simulation about gene flow and isolation.\n",
    "Those things still happen in this model, but natural selection adds an extra level of complexity. It can change the probability that certain genes move and persist on the landscape, and thus the nature of isolation between\n",
    "different neighborhoods on the lanscape.\n",
    "\n",
    "The following cells' code are very similar to what we saw in the first section. As before, run the code, changing any values as directed.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Set the parameter values</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, we'll create a model with <u>__two layers__</u>:\n",
    "\n",
    "1. A layer to describe the __distribution of habitat__, which will determine both the species' local population density and its movement throughout the landscape. \n",
    "    - We will use a uniform array from the neutral landscape model for this layer, as we saw in the previous simulation.\n",
    "    \n",
    "2. A layer to describe an __environmental gradient__ across the landscape, which will exert spatially divergent selection on the species.\n",
    "    - We will define this layer's array using the function that we defined at the beginning of this notebook.\n",
    "    \n",
    "We again want a model with a <u>__single species__</u>, but this time that species will have a __one trait__. Each individual's phenotype for this trait will determined by its genotype at the genetic loci that underlie the trait, and its fitness (i.e. probability of survival) will be determined by the difference between that phenotype and the optimal phenotype at the individual's location, according to the environmental gradient layer. Thus, the environmental gradient will exert spatially divergent natural selection on this trait."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To accomplish all of that, we have changed the following parameter values:\n",
    "\n",
    "- **dim** parameter: *We set the landscape dimensions to __(40, 40)__*.\n",
    "\n",
    "- **rast** parameter in the __'lyr_0'__ section: *We set this to a horizontal gradient, using __`make_unif_array(40)`__*.\n",
    "\n",
    "- **rast** parameter in the __'lyr_1'__ section: *We set this to a uniform array, using __`make_horz_grad_array(40)`__*\n",
    "\n",
    "- **layer** parameter in the __'trait_0'__ section: *We set this to __`'lyr_1'`__, to indicate that layer 1 should serve as trait 0's selective force*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# selection_demo.py\n",
    "\n",
    "# This is a parameters file generated by Geonomics\n",
    "# (by the gnx.make_parameters_file() function).\n",
    "\n",
    "\n",
    "                   #   ::::::          :::    :: :::::::::::#\n",
    "             #::::::    ::::   :::      ::    :: :: ::::::::::: ::#\n",
    "          #:::::::::     ::            ::   :::::::::::::::::::::::::#\n",
    "        #::::::::::                      :::::::::: :::::: ::::::::  ::#\n",
    "      #  : ::::  ::                    ::::  : ::    :::::::: : ::  :    #\n",
    "     # GGGGG :EEEE: OOOOO   NN   NN   OOOOO   MM   MM IIIIII  CCCCC SSSSS #\n",
    "    # GG     EE    OO   OO  NNN  NN  OO   OO  MM   MM   II   CC     SS     #\n",
    "    # GG     EE   OO     OO NN N NN OO     OO MMM MMM   II   CC     SSSSSS #\n",
    "    # GG GGG EEEE OO     OO NN  NNN OO     OO MM M MM   II   CC         SS #\n",
    "    # GG   G EE    OO   OO  NN   NN  OO   OO  MM   MM   II   CC        SSS #\n",
    "     # GGGGG :EEEE: OOOOO   NN   NN   OOOOO   MM   MM IIIIII  CCCCC SSSSS #\n",
    "      #    : ::::::::               :::::::::: ::              ::  :   : #\n",
    "        #:    :::::                    :::::: :::             :::::::  #\n",
    "          #    :::                      :::::  ::              ::::: #\n",
    "             #  ::                      ::::                      #\n",
    "                   #                                        #\n",
    "                      #  :: ::    :::             #\n",
    "\n",
    "\n",
    "params = {\n",
    "#--------------------------------------------------------------------------#\n",
    "\n",
    "#-----------------#\n",
    "#--- LANDSCAPE ---#\n",
    "#-----------------#\n",
    "    'landscape': {\n",
    "\n",
    "    #------------#\n",
    "    #--- main ---#\n",
    "    #------------#\n",
    "        'main': {\n",
    "            #x,y (a.k.a. j,i) dimensions of the Landscape\n",
    "            'dim':                      (40,40),\n",
    "            #x,y resolution of the Landscape\n",
    "            'res':                      (1,1),\n",
    "            #x,y coords of upper-left corner of the Landscape\n",
    "            'ulc':                      (0,0),\n",
    "            #projection of the Landscape\n",
    "            'prj':                      None,\n",
    "            }, # <END> 'main'\n",
    "\n",
    "    #--------------#\n",
    "    #--- layers ---#\n",
    "    #--------------#\n",
    "        'layers': {\n",
    "\n",
    "            #layer name (LAYER NAMES MUST BE UNIQUE!)\n",
    "            'lyr_0': {\n",
    "\n",
    "        #-------------------------------------#\n",
    "        #--- layer num. 0: init parameters ---#\n",
    "        #-------------------------------------#\n",
    "\n",
    "                #initiating parameters for this layer\n",
    "                'init': {\n",
    "\n",
    "                    #parameters for a 'defined'-type Layer\n",
    "                    'defined': {\n",
    "                        #raster to use for the Layer\n",
    "                        'rast':                   make_unif_array(40),\n",
    "                        #point coordinates\n",
    "                        'pts':                    None,\n",
    "                        #point values\n",
    "                        'vals':                   None,\n",
    "                        #interpolation method {None, 'linear', 'cubic',\n",
    "                        #'nearest'}\n",
    "                        'interp_method':          None,\n",
    "\n",
    "                        }, # <END> 'defined'\n",
    "\n",
    "                    }, # <END> 'init'\n",
    "\n",
    "                }, # <END> layer num. 0\n",
    "\n",
    "\n",
    "            #layer name (LAYER NAMES MUST BE UNIQUE!)\n",
    "            'lyr_1': {\n",
    "\n",
    "        #-------------------------------------#\n",
    "        #--- layer num. 1: init parameters ---#\n",
    "        #-------------------------------------#\n",
    "\n",
    "                #initiating parameters for this layer\n",
    "                'init': {\n",
    "\n",
    "                    #parameters for a 'defined'-type Layer\n",
    "                    'defined': {\n",
    "                        #raster to use for the Layer\n",
    "                        'rast':                   make_horz_grad_array(40),\n",
    "                        #point coordinates\n",
    "                        'pts':                    None,\n",
    "                        #point values\n",
    "                        'vals':                   None,\n",
    "                        #interpolation method {None, 'linear', 'cubic',\n",
    "                        #'nearest'}\n",
    "                        'interp_method':          None,\n",
    "\n",
    "                        }, # <END> 'defined'\n",
    "\n",
    "                    }, # <END> 'init'\n",
    "\n",
    "                }, # <END> layer num. 1\n",
    "\n",
    "\n",
    "\n",
    "    #### NOTE: Individual Layers' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Layers.\n",
    "\n",
    "\n",
    "            } # <END> 'layers'\n",
    "\n",
    "        }, # <END> 'landscape'\n",
    "\n",
    "\n",
    "#-------------------------------------------------------------------------#\n",
    "    \n",
    "#-----------------#\n",
    "#--- COMMUNITY ---#\n",
    "#-----------------#\n",
    "    'comm': {\n",
    "\n",
    "        'species': {\n",
    "\n",
    "            #species name (SPECIES NAMES MUST BE UNIQUE!)\n",
    "            'spp_0': {\n",
    "\n",
    "            #-----------------------------------#\n",
    "            #--- spp num. 0: init parameters ---#\n",
    "            #-----------------------------------#\n",
    "\n",
    "                'init': {\n",
    "                    #starting number of individs\n",
    "                    'N':                250,\n",
    "                    #carrying-capacity Layer name\n",
    "                    'K_layer':          'lyr_0',\n",
    "                    \n",
    "###########################################################################\n",
    "# CONSIDER CHANGING THIS PARAMETER                                        #\n",
    "                    #multiplicative factor for carrying-capacity layer    #\n",
    "                    'K_factor':         1,                                #\n",
    "###########################################################################\n",
    "                    \n",
    "                    }, # <END> 'init'\n",
    "\n",
    "            #-------------------------------------#\n",
    "            #--- spp num. 0: mating parameters ---#\n",
    "            #-------------------------------------#\n",
    "\n",
    "                'mating'    : {\n",
    "                    #age(s) at sexual maturity (if tuple, female first)\n",
    "                    'repro_age':                0,\n",
    "                    #whether to assign sexes\n",
    "                    'sex':                      False,\n",
    "                    #ratio of males to females\n",
    "                    'sex_ratio':                1/1,\n",
    "                    #whether P(birth) should be weighted by parental dist\n",
    "                    'dist_weighted_birth':       False,\n",
    "                    #intrinsic growth rate\n",
    "                    'R':                        0.5,\n",
    "                    #intrinsic birth rate (MUST BE 0<=b<=1)\n",
    "                    'b':                        0.2,\n",
    "                    #expectation of distr of n offspring per mating pair\n",
    "                    'n_births_distr_lambda':    1,\n",
    "                    #whether n births should be fixed at n_births_dist_lambda\n",
    "                    'n_births_fixed':           True,\n",
    "                    \n",
    "########################################################\n",
    "# CONSIDER CHANGING THIS PARAMETER                     #\n",
    "                    #radius of mate-search area        #   \n",
    "                    'mating_radius':            10,    #\n",
    "########################################################            \n",
    "                    \n",
    "                    }, # <END> 'mating'\n",
    "\n",
    "            #----------------------------------------#\n",
    "            #--- spp num. 0: mortality parameters ---#\n",
    "            #----------------------------------------#\n",
    "\n",
    "                'mortality'     : {\n",
    "                    #maximum age\n",
    "                    'max_age':                      None,\n",
    "                    #min P(death) (MUST BE 0<=d_min<=1)\n",
    "                    'd_min':                        0,\n",
    "                    #max P(death) (MUST BE 0<=d_max<=1)\n",
    "                    'd_max':                        1,\n",
    "                    #width of window used to estimate local pop density\n",
    "                    'density_grid_window_width':    None,\n",
    "                    }, # <END> 'mortality'\n",
    "\n",
    "            #---------------------------------------#\n",
    "            #--- spp num. 0: movement parameters ---#\n",
    "            #---------------------------------------#\n",
    "\n",
    "                'movement': {\n",
    "                    #whether or not the species is mobile\n",
    "                    'move':                                 True,\n",
    "                    #mode of distr of movement direction\n",
    "                    'direction_distr_mu':                   0,\n",
    "                    #concentration of distr of movement direction\n",
    "                    'direction_distr_kappa':                0,\n",
    "##################################################################\n",
    "# CONSIDER CHANGING THIS PARAMETER                               #\n",
    "                    #1st param of distr of movement distance     #   \n",
    "                    'movement_distance_distr_param1':       0.5, #\n",
    "##################################################################\n",
    "                    #2nd param of distr of movement distance\n",
    "                    'movement_distance_distr_param2':       5e-8,\n",
    "                    #movement distance distr to use ('levy' or 'wald')\n",
    "                    'movement_distance_distr':              'levy',\n",
    "                    #1st param of distr of dispersal distance\n",
    "                    'dispersal_distance_distr_param1':      0,\n",
    "                    #2nd param of distr of dispersal distance\n",
    "                    'dispersal_distance_distr_param2':      5e-14,\n",
    "                    #dispersal distance distr to use ('levy' or 'wald')\n",
    "                    'dispersal_distance_distr':             'levy',\n",
    "                    'move_surf'     : {\n",
    "                        #move-surf Layer name\n",
    "                        'layer':                'lyr_0',\n",
    "                        #whether to use mixture distrs\n",
    "                        'mixture':              True,\n",
    "                        #concentration of distrs\n",
    "                        'vm_distr_kappa':       12,\n",
    "                        #length of approximation vectors for distrs\n",
    "                        'approx_len':           5000,\n",
    "                        }, # <END> 'move_surf'\n",
    "\n",
    "                    },    # <END> 'movement'\n",
    "\n",
    "\n",
    "            #---------------------------------------------------#\n",
    "            #--- spp num. 0: genomic architecture parameters ---#\n",
    "            #---------------------------------------------------#\n",
    "\n",
    "                'gen_arch': {\n",
    "                    #file defining custom genomic arch\n",
    "                    'gen_arch_file':            None,\n",
    "                    #num of loci\n",
    "                    'L':                        100,\n",
    "                    #starting allele frequency (None to draw freqs randomly)\n",
    "                    'start_p_fixed':            0.5,\n",
    "                    #whether to start neutral locus freqs at 0\n",
    "                    'start_neut_zero':          False,\n",
    "                    #genome-wide per-base neutral mut rate (0 to disable)\n",
    "                    'mu_neut':                  0,\n",
    "                    #genome-wide per-base deleterious mut rate (0 to disable)\n",
    "                    'mu_delet':                 0,\n",
    "                    #shape of distr of deleterious effect sizes\n",
    "                    'delet_alpha_distr_shape':  0.2,\n",
    "                    #scale of distr of deleterious effect sizes\n",
    "                    'delet_alpha_distr_scale':  0.2,\n",
    "                    #alpha of distr of recomb rates\n",
    "                    'r_distr_alpha':            None,\n",
    "                    #beta of distr of recomb rates\n",
    "                    'r_distr_beta':             None,\n",
    "                    #whether loci should be dominant (for allele '1')\n",
    "                    'dom':                      False,\n",
    "                    #whether to allow pleiotropy\n",
    "                    'pleiotropy':               False,\n",
    "                    #custom fn for drawing recomb rates\n",
    "                    'recomb_rate_custom_fn':    None,\n",
    "                    #number of recomb paths to hold in memory\n",
    "                    'n_recomb_paths_mem':       int(1e4),\n",
    "                    #total number of recomb paths to simulate\n",
    "                    'n_recomb_paths_tot':       int(1e5),\n",
    "                    #num of crossing-over events (i.e. recombs) to simulate\n",
    "                    'n_recomb_sims':            10_000,\n",
    "                    #whether to generate recombination paths at each timestep\n",
    "                    'allow_ad_hoc_recomb':      False,\n",
    "                    #whether to save mutation logs\n",
    "                    'mut_log':                  False,\n",
    "\n",
    "                    'traits': {\n",
    "\n",
    "                        #--------------------------#\n",
    "                        #--- trait 0 parameters ---#\n",
    "                        #--------------------------#\n",
    "                        #trait name (TRAIT NAMES MUST BE UNIQUE!)\n",
    "                        'trait_0': {\n",
    "                            #trait-selection Layer name\n",
    "                            'layer':                'lyr_1',\n",
    "#################################################################\n",
    "# CONSIDER CHANGING THESE PARAMETERS                            #\n",
    "                            #polygenic selection coefficient    #\n",
    "                            'phi':                  0.05,       # \n",
    "                            #number of loci underlying trait    #\n",
    "                            'n_loci':               1,          #\n",
    "#################################################################\n",
    "                            #mutation rate at loci underlying trait\n",
    "                            'mu':                   0,\n",
    "                            #mean of distr of effect sizes\n",
    "                            'alpha_distr_mu' :      0.1,\n",
    "                            #variance of distr of effect size\n",
    "                            'alpha_distr_sigma':    0,\n",
    "                            #max allowed magnitude for an alpha value\n",
    "                            'max_alpha_mag':        None,\n",
    "                            #curvature of fitness function\n",
    "                            'gamma':                1,\n",
    "                            #whether the trait is universally advantageous\n",
    "                            'univ_adv':             False\n",
    "                            }, # <END> trait 0\n",
    "\n",
    "\n",
    "    #### NOTE: Individual Traits' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Traits.\n",
    "\n",
    "\n",
    "                        }, # <END> 'traits'\n",
    "\n",
    "                    }, # <END> 'gen_arch'\n",
    "\n",
    "\n",
    "                }, # <END> spp num. 0\n",
    "\n",
    "\n",
    "\n",
    "    #### NOTE: individual Species' sections can be copy-and-pasted (and\n",
    "    #### assigned distinct keys and names), to create additional Species.\n",
    "\n",
    "\n",
    "            }, # <END> 'species'\n",
    "\n",
    "        }, # <END> 'comm'\n",
    "\n",
    "\n",
    "#------------------------------------------------------------------------#\n",
    "\n",
    "#-------------#\n",
    "#--- MODEL ---#\n",
    "#-------------#\n",
    "    'model': {\n",
    "        #total Model runtime (in timesteps)\n",
    "        'T':            100,\n",
    "        #min burn-in runtime (in timesteps)\n",
    "        'burn_T':       30,\n",
    "        #seed number\n",
    "        'num':          None,\n",
    "        #time step interval for simplication of tskit tables\n",
    "        'tskit_simp_interval':      100,\n",
    "\n",
    "\n",
    "        #-----------------------------#\n",
    "        #--- iterations parameters ---#\n",
    "        #-----------------------------#\n",
    "        'its': {\n",
    "            #num iterations\n",
    "            'n_its':            1,\n",
    "            #whether to randomize Landscape each iteration\n",
    "            'rand_landscape':   False,\n",
    "            #whether to randomize Community each iteration\n",
    "            'rand_comm':        False,\n",
    "            #whether to burn in each iteration\n",
    "            'repeat_burn':      False,\n",
    "            }, # <END> 'iterations'\n",
    "\n",
    "\n",
    "\n",
    "        } # <END> 'model'\n",
    "\n",
    "    } # <END> params\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Make the model</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the following code again to make the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make our params dict into a proper Geonomics ParamsDict object\n",
    "params = gnx.make_params_dict(params, 'selection_demo')\n",
    "#then use it to make a model\n",
    "mod = gnx.make_model(parameters=params, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then take a look at the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mod)\n",
    "%matplotlib inline\n",
    "mod.plot(spp=0, lyr=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\" role=\"alert\" style=\"font-size:120%\">\n",
    "    <b><h2>Run the model</h2></b> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll burn the model in, then see what we've got:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SCROLL DOWN TO FOLLOW BURN-IN AND SEE PLOT\n",
    "\n",
    "mod.walk(T=10000, mode='burn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.plot(spp=0, lyr=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, before we start running the main section of the model, we can use a built-in Geonomics function to make a map of the individuals, colored by their phenotypes for the trait.\n",
    "\n",
    "Note that the closer an individual's color is to its environmental background, the closer that individual is to the optimum phenotype for its location, and thus the higher the individual's fitness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.plot_phenotype(spp=0, trt=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's also create our genetic relatedness plots again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_PCA(mod)\n",
    "map_PCA(mod, lyr_num=1, mask=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's run our model for 50 timesteps, then create our plots again, up until time step 150. Keep a close eye on all three plots, and think about why you see what you do."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.walk(50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.plot_phenotype(spp=0, trt=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_PCA(mod)\n",
    "map_PCA(mod, lyr_num=1, mask=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.walk(50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.plot_phenotype(spp=0, trt=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_PCA(mod)\n",
    "map_PCA(mod, lyr_num=1, mask=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mod.walk(50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod.plot_phenotype(spp=0, trt=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_PCA(mod)\n",
    "map_PCA(mod, lyr_num=1, mask=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- BEGIN QUESTION -->\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>QUESTION 4: How can you tell natural selection is acting on the population? Use evidence from the plots. Be specific.</b>\n",
    "    <br />    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Q4: YOUR RESPONSE HERE\n",
    "\n",
    "<!-- END QUESTION -->\n",
    "\n",
    "<!-- BEGIN QUESTION -->\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>QUESTION 5: Give an example of a condition that needs to be met for the genetic divergence of these populations to lead to speciation. Discuss. </b>\n",
    "    <br />    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q5: YOUR RESPONSE HERE\n",
    "\n",
    "<!-- END QUESTION -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>EXERCISE:</b>\n",
    "    <br />\n",
    "    Tweak the model and rerun.\n",
    "    In the remaining time, go back and <u>create your own scenario, by changing one or more of the following parameters</u> as you wish. Note that <b>each of these parameters is outlined by a comment box within the code above</b>.\n",
    "\n",
    "    \n",
    "<b>Parameters to tweak</b>:\n",
    "- ***phi***: *This controls the strength of selection on the trait. It can be set to 0 <= value <= 1.*\n",
    "\n",
    "- ***n_loci***: *This controls the number of genetic loci underlying the trait's phenotype. For this simulation, it can be set to 1 <= value <= 100.*\n",
    "\n",
    "- ***K_factor***: *This controls the population density, by setting the number of individuals that can inhabit each cell. It can be set to any value > 0, technically, but values larger than 3 or 4 will probably require too much memory to run on this server, so be careful!*\n",
    "\n",
    "- ***movement_distance_distr_param_1***: *This controls the average distance that an individual moves each time step (in units of cell-widths), and thus how mobile the species is. For this simulation, it can be set to any value > 0 (but very large values might cause problems!).*\n",
    "\n",
    "- ***mating_radius***: *This controls the radius (expressed in cell-width units) within which an indivudal can choose a mate. For this simulation, it can be set to any value > 0.*\n",
    "\n",
    "If you wish, you should feel free to also toy with any other parameters. But please know that by doing so you may break the code and be unable to run your simulation.\n",
    "\n",
    "Once you've made your desired changes, <u>rerun the code above</u>, and make observations about the effect of the change you made.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<!-- BEGIN QUESTION -->\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>QUESTION 6: Which parameter did you change? What effect does this have and what change did you expect this would have on the evolutionary scenario?</b>\n",
    "    <br />    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Q6: YOUR RESPONSE HERE\n",
    "\n",
    "<!-- END QUESTION -->\n",
    "\n",
    "<!-- BEGIN QUESTION -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>QUESTION 7: What do the results suggest? Use evidence from the plots of the simulation to formulate your answer. </b>\n",
    "    <br />    \n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Q7: YOUR RESPONSE HERE\n",
    "\n",
    "<!-- END QUESTION -->\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- BEGIN QUESTION -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# run this cell to generate your final plots for grading.\n",
    "mod.plot_phenotype(spp=0, trt=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# run this cell to generate your final plots for grading.\n",
    "plot_PCA(mod)\n",
    "map_PCA(mod, lyr_num=1, mask=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- END QUESTION -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run this cell to generate pdf\n",
    "# right click on the link and open in a new window/tab\n",
    "from IPython.display import display, HTML\n",
    "export_notebook(\"geonomics_biogeo_2020.ipynb\", filtering = True, pagebreaks = False) \n",
    "display(HTML(\"Download your PDF <a href='geonomics_biogeo_2020.pdf' download>here</a>.\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Turning in your answers to bCourses\n",
    "\n",
    "**Make sure you have completed the exercise of creating your own scenario and answered questions 1-7.**\n",
    "\n",
    "You are finished with this notebook! Please run the above cell to generate a download link for your submission file.\n",
    "\n",
    "Open a new tab and go to https://datahub.berkeley.edu, navigate to the \"geonomics_biogeo_demo\" folder, click the `geonomics_biogeo_2020.pdf` to view your submission and download it.\n",
    "\n",
    "**If you would like to resubmit, delete the pdf from datahub and run the cell again.** (You can delete it by checking the box and clicking the red button that appears)\n",
    "\n",
    "Then upload this answer sheet to bCourses for grading."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook was developed by Drew Hart, Natalie Graham, Monica Wilkinson and Keeley Takimoto\n",
    "\n",
    "October 2019"
   ]
  }
 ],
 "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": 2
}