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

[2]:
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, task=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.
Checking folder existance...output exists.
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)
Seed: 1
Speciation rate: 0.1
Job Type: 1
Max time: 3600
Deme: 10
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
Getting geo transform...done.
Affine transform is -78.375, 0.00833333, 0, 0.858333, 0, -0.00833333
Importing ../../pycoalescence/tests/sample/SA_sample_fine.tif ...................done.
No data value is: 0
Getting geo transform...done.
Affine transform is -78.375, 0.00833333, 0, 0.858333, 0, -0.00833333
Using dispersal kernel.
Initial count is 380980
Setting up simulation...done.
Number of individuals simulating: 380980
*************************************************
Beginning simulations...done.
Finalising data...done.
Creating SQL database file...
        Checking for existing folders....
        Generating species list....
        Executing SQL commands....
No additional speciation rates to apply.
Speciation rate is: 0.1.
Time is: 0.
Applying speciation rate 0.1 at time 0...
        Generating biodiversity...
        Calculating coalescence tree...
        Assigning species IDs...
        Calculating species abundances...
        Number of species: 96296
        Generating SPECIES_ABUNDANCES table...
        Writing to output/data_1_1.db...
Total generations simulated (steps): 191.025 (963375)
Setup time was 0 minutes 0 seconds
Simulation time was 0 hours 0 minutes 1 seconds
File output and species calculation time was 0 minutes 1 seconds
SQL output time was 0 minutes 0 seconds
Total simulation and output time was 0 hours 0 minutes 2 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 COALESCENCE TREE CALCULATIONS
Input file is output/data_1_1.db
Speciation rates are: 0.1, 0.2, 0.4.
Beginning data import...
        Detected 665665 events in the coalescence tree.
Beginning data import...done.
Getting previous calculations...previous calculations detected.
Calculation already performed for speciation rate=0.1, time=0 and protracted parameters 0, 0
Applying speciation rate 0.2 at time 0...
        Generating biodiversity...
        Calculating coalescence tree...
        Assigning species IDs...
        Calculating species abundances...
        Number of species: 152176
        Generating SPECIES_ABUNDANCES table...
Applying speciation rate 0.4 at time 0...
        Generating biodiversity...
        Calculating coalescence tree...
        Assigning species IDs...
        Calculating species abundances...
        Number of species: 232205
        Generating SPECIES_ABUNDANCES table...
Writing out to output/data_1_1.db...

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_species_richness(reference)))
Species richness for speciation rate of 0.1 is 96296
Species richness for speciation rate of 0.2 is 152176
Species richness for speciation rate of 0.4 is 232205

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")
[10]:
sim2 = Simulation(logging_level=50)
sim2.set_simulation_parameters(seed=2, task=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_file=pristine_fine_file1, coarse_file=pristine_coarse_file1, time=100,
                                            rate=0.5)
sim2.add_historical_map(fine_file=pristine_fine_file2, coarse_file=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.

[11]:
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_species_richness(reference)))
Species richness at time 0.0 with speciation rate of 0.1 is 7386.
Species richness at time 50.0 with speciation rate of 0.1 is 7332.
Species richness at time 100.0 with speciation rate of 0.1 is 7420.
Species richness at time 0.0 with speciation rate of 0.2 is 8756.
Species richness at time 50.0 with speciation rate of 0.2 is 8742.
Species richness at time 100.0 with speciation rate of 0.2 is 8751.
Species richness at time 0.0 with speciation rate of 0.3 is 9568.
Species richness at time 50.0 with speciation rate of 0.3 is 9579.
Species richness at time 100.0 with speciation rate of 0.3 is 9607.

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.

[12]:
sim3 = Simulation(logging_level=50)
sim3.set_simulation_parameters(seed=3, task=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.

[13]:
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

[15]:
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_species_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: 61
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: 52
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: 52
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: 376
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: 405
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: 422
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: 259
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: 274
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: 280

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.

[16]:
sim1 = None
sim2 = None
sim3 = None
tree1 = None
tree2 = None
tree3 = None