Example simulations and analysis using pycoalescence¶
This jupyter notebook shows a variety of full example simulations that can be run using the files provided in the sample folder in pycoalescence. They are intended to be used as templates for your own simulations.
Note that this notebook requires Python 3.x and a full install of pycoalescence.
First import the necessary modules
[1]:
import os
from pycoalescence import Simulation, CoalescenceTree
Define our input and output directories - change these as required
[3]:
input_dir = "pycoalescence/tests/sample"
output_dir = "output"
Basic spatial simulation and analysis¶
Basic spatial simulation on just one map file, with a closed landscape. For this example we run with full logging (at ‘debug’ level). All other examples will be run using ‘critical’ level, which shouldn’t output any information unless there is a problem.
Use a speciation rate of 0.1, a dispersal of 2.0 (sigma) and 10 individuals per cell (deme).
The output will be stored in output/data_1_1.db
[4]:
sim1 = Simulation(logging_level=10)
sim1.set_simulation_parameters(seed=1, job_type=1, output_directory=output_dir, min_speciation_rate=0.1,
sigma=2.0, deme=10)
# Use automatic detection of the map file dimensions
sim1.set_map(os.path.join(input_dir, "SA_sample_fine.tif"))
sim1.run()
Checking folder existance...pycoalescence/tests/sample/SA_sample_fine.tif exists!
Output folder does not exist... creating...done!
Checking folder existance...done!
Checking for unfinished simulations...done!
No files found containing unfinished simulations.
*************************************************
Setting up simulation...
Dispersal (tau, sigma): 1, 2
Dispersal method: normal
Fine map
-file: pycoalescence/tests/sample/SA_sample_fine.tif
-dimensions: (13, 13)
-offset: (0, 0)
Coarse map
-file: none
-dimensions: (13, 13)
-offset: (0, 0)
-scale: 1
Sample grid
-dimensions: (13, 13)
-optimised area: (13, 13)
-optimised offsets: (0, 0)
-deme: 10
Seed: 1
Speciation rate: 0.1
Job Type: 1
Max time: 3600
-deme sample: 1
Output directory: output
Disp Rel Cost: 1
Times: 0.0
Checking folder existance...pycoalescence/tests/sample/SA_sample_fine.tif exists!
Importing pycoalescence/tests/sample/SA_sample_fine.tif No data value is: 0
Importing pycoalescence/tests/sample/SA_sample_fine.tif ...................done!
No data value is: 0
Using dispersal kernel.
Initial count is 380980
Setting up simulation...done!
Number of individuals simulating: 380980
*************************************************
Beginning simulations...done!
No additional speciation rates to apply.
Speciation rate is: 0.1.
Finalising data...done!
Creating SQL database file...
Executing SQL commands....done!
0.000000,Calculating generation 0.000000
Calculating tree structure...done!
Number of species: 96371
Generating new SQL table for speciation rate 0.1...done!
Writing to output/data_1_1.db .... done!
Total generations simulated (steps): 223.285 (964919)
Setup time was 0 minutes 0 seconds
Simulation time was 0 hours 0 minutes 2 seconds
File output and species calculation time was 0 minutes 5 seconds
SQL output time was 0 minutes 0 seconds
Total simulation and output time was 0 hours 0 minutes 7 seconds
Re-create the coalescence tree for different speciation rates
[5]:
tree1 = CoalescenceTree(sim1, logging_level=10)
tree1.set_speciation_parameters(speciation_rates=[0.1, 0.2, 0.4])
tree1.apply()
No sample file provided, defaulting to null.No times provided, defaulting to 0.0.***************************
STARTING CALCULATIONS
Input file is output/data_1_1.db
Speciation rates are: 0.1, 0.2, 0.4.
Protracted speciation parameters (min, max) are:
0, 0
Beginning data import...done
Calculating speciation rate 0.1
Calculating generation 0
calculation already performed for 0.1 at time 0
calculation already performed for 0.1 at time 0
Calculating speciation rate 0.2
Calculating generation 0
Calculating tree structure...done!
Number of species: 152584
Generating new SQL table for speciation rate 0.2...done!
Calculating generation 0
Calculating speciation rate 0.4
Calculating generation 0
Calculating tree structure...done!
Number of species: 232829
Generating new SQL table for speciation rate 0.4...done!
Writing out to output/data_1_1.db...done!
Print out the species richness for each speciation rate by looping over our community parameters.
[6]:
for reference in tree1.get_community_references():
spec_rate = tree1.get_community_parameters(reference=reference)["speciation_rate"]
print("Species richness for speciation rate of {} is {}".format(spec_rate,
tree1.get_richness(reference)))
Species richness for speciation rate of 0.1 is 96371
Species richness for speciation rate of 0.2 is 152584
Species richness for speciation rate of 0.4 is 232829
More complex spatial example¶
Here we use a sample area, a fine map file (defining high-resolution density around the area) and a coarse map file (defining the low-resolution density over a larger area).
Additionally, we provide a historical fine and coarse maps to define the density at 100 generations in the past (with a rate of 0.5) and 200 generations in the past (with a rate of 0.0).
Within the simulation, we sample the community at times 0, 50 and 100 generations in the past.
The output will be stored in output/data_1_2.db
Define our file paths
[7]:
sample_file = os.path.join(input_dir, "SA_samplemaskINT.tif")
fine_file = os.path.join(input_dir, "SA_sample_fine.tif")
coarse_file = os.path.join(input_dir, "SA_sample_coarse.tif")
pristine_fine_file1 = os.path.join(input_dir, "SA_sample_fine_pristine1.tif")
pristine_coarse_file1 = os.path.join(input_dir, "SA_sample_coarse_pristine1.tif")
pristine_fine_file2 = os.path.join(input_dir, "SA_sample_fine_pristine2.tif")
pristine_coarse_file2 = os.path.join(input_dir, "SA_sample_coarse_pristine2.tif")
[8]:
sim2 = Simulation(logging_level=50)
sim2.set_simulation_parameters(seed=2, job_type=1, output_directory=output_dir, min_speciation_rate=0.1,
sigma=2.0, deme=1)
sim2.set_speciation_rates([0.1, 0.2, 0.3])
# Use automatic detection of the map file dimensions
sim2.set_map_files(sample_file=sample_file, fine_file=fine_file, coarse_file=coarse_file)
sim2.add_historical_map(fine_map=pristine_fine_file1, coarse_map=pristine_coarse_file1, time=100,
rate=0.5)
sim2.add_historical_map(fine_map=pristine_fine_file2, coarse_map=pristine_coarse_file2, time=100,
rate=0.5)
sim2.add_sample_time(50)
sim2.add_sample_time(100)
sim2.run()
Print the species richness for each time and speciation rate.
[9]:
tree2 = CoalescenceTree(sim2)
for reference in tree2.get_community_references():
# Contains a dictionary of the parameters
community_parameters = tree2.get_community_parameters(reference)
spec_rate = community_parameters["speciation_rate"]
time = community_parameters["time"]
print("Species richness at time {} with speciation rate of {} is {}.".format(time,
spec_rate,
tree2.get_richness(reference)))
Species richness at time 0.0 with speciation rate of 0.1 is 7392.
Species richness at time 50.0 with speciation rate of 0.1 is 7312.
Species richness at time 100.0 with speciation rate of 0.1 is 7456.
Species richness at time 0.0 with speciation rate of 0.2 is 8808.
Species richness at time 50.0 with speciation rate of 0.2 is 8691.
Species richness at time 100.0 with speciation rate of 0.2 is 8764.
Species richness at time 0.0 with speciation rate of 0.3 is 9656.
Species richness at time 50.0 with speciation rate of 0.3 is 9515.
Species richness at time 100.0 with speciation rate of 0.3 is 9647.
Example using protracted speciation and a metacommunity¶
This is the same as the first example, but with protracted speciation preventing speciation from occuring before 1000 generations, and forcing speciation to occur at 10000 generations.
The coalescence tree is generated with and without a metacommunity as well.
[10]:
sim3 = Simulation(logging_level=50)
sim3.set_simulation_parameters(seed=3, job_type=1, output_directory=output_dir,
min_speciation_rate=0.1, sigma=2.0, deme=1,
protracted=True, min_speciation_gen=1000, max_speciation_gen=10000)
sim3.set_speciation_rates([0.1, 0.2, 0.3])
# Use automatic detection of the map file dimensions
sim3.set_map(os.path.join(input_dir, "SA_sample_fine.tif"))
sim3.run()
Re-create the coalescence tree using sampling from a metacommunity instead of a speciation rate.
[11]:
tree3 = CoalescenceTree(sim3)
tree3.set_speciation_parameters(speciation_rates=[0.1, 0.2, 0.3], protracted_speciation_min=1000,
protracted_speciation_max=10000, metacommunity_speciation_rate=0.001,
metacommunity_size=100000)
tree3.add_protracted_parameters(10, 5000)
tree3.add_protracted_parameters(50, 800)
tree3.apply()
Print out the species richness for each parameter set
[11]:
for reference in tree3.get_community_references():
params = tree3.get_community_parameters(reference)
spec_rate = params["speciation_rate"]
proc_min = params["min_speciation_gen"]
proc_max = params["max_speciation_gen"]
print("Speciation rate of {}, (min gen={}, max gen={})".format(spec_rate, proc_min, proc_max))
if params["metacommunity_reference"] != 0:
meta_params = tree3.get_metacommunity_parameters(params["metacommunity_reference"])
meta_spec = meta_params["speciation_rate"]
meta_size = meta_params["metacommunity_size"]
print("Metacommunity used of size {} with speciation rate {}.".format(meta_size, meta_spec))
print("Species richness: {}".format(tree3.get_richness(reference)))
Speciation rate of 0.1, (min gen=1000.0, max gen=10000.0)
Species richness: 71
Speciation rate of 0.2, (min gen=1000.0, max gen=10000.0)
Species richness: 71
Speciation rate of 0.3, (min gen=1000.0, max gen=10000.0)
Species richness: 71
Speciation rate of 0.1, (min gen=1000.0, max gen=10000.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 71
Speciation rate of 0.2, (min gen=1000.0, max gen=10000.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 71
Speciation rate of 0.3, (min gen=1000.0, max gen=10000.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 71
Speciation rate of 0.1, (min gen=10.0, max gen=5000.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 657
Speciation rate of 0.2, (min gen=10.0, max gen=5000.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 657
Speciation rate of 0.3, (min gen=10.0, max gen=5000.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 657
Speciation rate of 0.1, (min gen=50.0, max gen=800.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 656
Speciation rate of 0.2, (min gen=50.0, max gen=800.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 657
Speciation rate of 0.3, (min gen=50.0, max gen=800.0)
Metacommunity used of size 100000.0 with speciation rate 0.001.
Species richness: 657
Clean up¶
Delete the objects from memory to clear from RAM - this is not usually required as objects will be deleted when they fall out of scope.
[12]:
sim1 = None
sim2 = None
sim3 = None
tree1 = None
tree2 = None
tree3 = None