tskit-dev / pyslim

Tools for dealing with tree sequences coming to and from SLiM.
MIT License
27 stars 23 forks source link

add method to compute discretized population density #186

Closed giliapatterson closed 3 years ago

giliapatterson commented 3 years ago

Make a method to compute the population size in intervals of x, y, and time from a tree sequence. It should take a tree sequence, bins in the x and y directions, and bins in time and return an easy to plot 3-d array of the number of individuals alive in each x and y bin averaged across the times in each time bin.

def population_size(ts, x_bins = 1, y_bins = 1, time_bins = 1)

If x_bins, y_bins, or time_bins is a single number, the breakpoints will be np.linspace(0, max(x), xbins + 1). If any of the bins is an array it will be the breakpoints.

Returns a 3-d array with dimensions (number of x breakpoints - 1) by (number of y breakpoints - 1) by (number of time breakpoints - 1)

The goal is to make it easy to make plots like this:

import matplotlib.pyplot as plt
x_bins = [0, 10, 20]
time_bins = np.linspace(0, 1000, 101)
popsize = pyslim.population_size(ts, x_bins = x_bins, time_bins = time_bins)
plt.pcolormesh(time_bins, x_bins, popsize)

image