example_fragmented_landscape.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# pycoalescence fragmented landscape example\n",
    "\n",
    "The process outlined here simulates a community of 18010 individuals on a real landscape as it progresses towards a randomly-cleared landscape of 2495 individuals over 100 generations. The community is sampled every 10 generations. The simulation is repeated 10 times and the community is generated for three different speciation rates.\n",
    "\n",
    "### Inputs\n",
    "\n",
    "- The following map files in the \"maps\" folder\n",
    "\t- \"fine_present.tif\"\n",
    "\t- \"coarse_present.tif\"\n",
    "\t- \"fine_historical.tif\"\n",
    "\t- \"coarse_historical.tif\"\n",
    "- A speciation rates of 0.0001 (the minimum), 0.0005 and 0.001\n",
    "- A dispersal sigma of 16\n",
    "- A deme size of 10\n",
    "\n",
    "### Outputs\n",
    "\n",
    "- 10 databases (one per simulation) in the \"output\" folder, containing the following features in the respective SQL table:\n",
    "\t- The simulation parameters (SIMULATION_PARAMETERS)\n",
    "\t- The community parameters (COMMUNITY_PARAMETERS)\n",
    "\t- The full coalescence tree (SPECIES_LIST)\n",
    "\t- The species abundances for each community (SPECIES_ABUNDANCES)\n",
    "- A csv containing the species richness at each time\n",
    "- A csv containing the species abundances 0, 50 and 100 generations before present"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycoalescence import Simulation, CoalescenceTree\n",
    "import pandas as pd\n",
    "import os\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the simulation and generating the communities\n",
    "\n",
    "The `Simulation` class performs the simulation. The `CoalescenceTree` class calculates the coalescence tree for each parameter set and generates the community."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sim 1 complete.\n",
      "Sim 2 complete.\n",
      "Sim 3 complete.\n",
      "Sim 4 complete.\n",
      "Sim 5 complete.\n",
      "Sim 6 complete.\n",
      "Sim 7 complete.\n",
      "Sim 8 complete.\n",
      "Sim 9 complete.\n",
      "Sim 10 complete.\n",
      "Total time: 624.9451930522919 seconds\n"
     ]
    }
   ],
   "source": [
    "# Run 10 repeat simulations, each with a different random number seed\n",
    "start_time = time.time()\n",
    "for seed in range(1, 11, 1):\n",
    "    # Set up and run the simulation\n",
    "    sim = Simulation(logging_level=30)\n",
    "    sim.set_simulation_parameters(seed=seed, task=1, output_directory=\"output\",\n",
    "                                  min_speciation_rate=0.000001, sigma=16, deme=5)\n",
    "    if not os.path.exists(os.path.join(\"output\", \"data_{}_{}.db\".format(1, seed))):\n",
    "        sim.set_map_files(sample_file=\"null\", \n",
    "                          fine_file=os.path.join(\"maps\", \"fine_present.tif\"),\n",
    "                          coarse_file=os.path.join(\"maps\", \"coarse_present.tif\"))\n",
    "#         sim.set_map_files(sample_file=\"null\", fine_file=os.path.join(\"maps\", \"fine_historical.tif\"),\n",
    "#                                coarse_file=os.path.join(\"maps\", \"coarse_historical.tif\"))\n",
    "        sim.add_historical_map(fine_file=os.path.join(\"maps\", \"fine_historical.tif\"),\n",
    "                               coarse_file=os.path.join(\"maps\", \"coarse_historical.tif\"), time=99.9)\n",
    "        sim.add_sample_time([x for x in range(0, 160, 10)])\n",
    "        # This is the main time-consuming step\n",
    "        sim.run()\n",
    "    print(\"Sim {} complete.\".format(seed))\n",
    "    # Now generate the coalescence tree for each speciation rate and time\n",
    "    coal_tree = CoalescenceTree(sim)\n",
    "    coal_tree.set_speciation_parameters(speciation_rates=[0.000001, 0.000005, 0.00001],\n",
    "                                        times=[float(x) for x in range(0, 160, 10)])\n",
    "    coal_tree.apply()\n",
    "    # Calculate the species richness for each parameter and store it in a new table\n",
    "    coal_tree.calculate_richness()\n",
    "end_time = time.time()\n",
    "print(\"Total time: {} seconds\".format(end_time - start_time))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Store the species richness at each time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = []\n",
    "for seed in range(1, 11, 1):\n",
    "\tcoal_tree = CoalescenceTree(os.path.join(\"output\", \"data_{}_{}.db\".format(1, seed)))\n",
    "\tfor ref in coal_tree.get_community_references():\n",
    "\t\ttime = coal_tree.get_community_parameters(ref)[\"time\"]\n",
    "\t\tsr = coal_tree.get_species_richness(ref)\n",
    "\t\tspec_rate = coal_tree.get_community_parameters(ref)[\"speciation_rate\"]\n",
    "\t\toutput.append({\"seed\" : seed, \"time\" : time, \"species_richness\" : sr, \"speciation_rate\" : spec_rate})\n",
    "output_df = pd.DataFrame(output)\n",
    "output_df.to_csv(os.path.join(\"results\", \"richness.csv\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Store the species abundances at select times\n",
    "\n",
    "Just store the time for 0, 50 and 100 generations before the present. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = []\n",
    "for seed in range(1, 11, 1):\n",
    "\tcoal_tree = CoalescenceTree(os.path.join(\"output\", \"data_{}_{}.db\".format(1, seed)))\n",
    "\tfor ref in coal_tree.get_community_references():\n",
    "\t\ttime = coal_tree.get_community_parameters(ref)[\"time\"]\n",
    "\t\tif time in [0, 50, 100, 150]:\n",
    "\t\t\tspec_abun = coal_tree.get_species_abundances(reference=ref)\n",
    "\t\t\tspeciation_rate = coal_tree.get_community_parameters(ref)[\"speciation_rate\"]\n",
    "\t\t\tfor row in spec_abun:\n",
    "\t\t\t\toutput.append({\"time\" : time, \"species_id\" : row[0], \"abundance\" : row[1], \"seed\" : seed,\n",
    "\t\t\t\t\t\t\t   \"speciation_rate\" : speciation_rate})\n",
    "output_df = pd.DataFrame(output)\n",
    "output_df.to_csv(os.path.join(\"results\", \"species_abundances.csv\"))"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}