Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Dive into Scientific Python

Yoav Ram
March 15, 2017

Dive into Scientific Python

I introduce the Python programming language and demonstrate how Scientific Python can be used to study evolutionary theory using mathematical and computational models. You'll see how to run fast evolutionary simulations with NumPy and Cython, analyze and visualize simulation results with Pandas and Seaborn, and find solutions to evolutionary models using SciPy. This presentation is a wonderful opportunity to learn about Scientific Python through actual research-based examples, as well as an occasion to to discover how theoretical evolutionary biologists approach their research.

See full repo with source code at: https://github.com/yoavram/DataTalks2017

Yoav Ram

March 15, 2017
Tweet

More Decks by Yoav Ram

Other Decks in Research

Transcript

  1. Who I am • Yoav Ram @yoavram • Postdoc at

    Stanford University • PhD in BioMath from Tel-Aviv University • Using Python since 2002 • Using & teaching Scientific Python since 2011 • Python training for engineers & data scientists Presentation & source code on GitHub: https://github.com/yoavram/DataTalks2017 License: CC-BY-SA 4.0 2
  2. Why Python? I used Matlab. Now I use Python. by

    Steve Tjoa Why use Python for scientific computing? for scientific computing…
  3. Gratis: Free as in Beer • MATLAB is expensive (as

    of Feb 2017) – Individuals: $2,350 – Academia: $550 – Personal: $85 – Student: $29-55 – Batteries (toolboxes…) not included • Python is totally free – Batteries included (NumPy, SciPy…) • R is also free MathWork Pricing
  4. Libre: Free as in Speech • MATLAB source code is

    closed and proprietary – You cannot see the code – You cannot change the code – You can participate in the discussion as a client • Python source code is open – You can see, you can change, you can contribute code and documentation (python, numpy) – You can participate in the discussion as a peer (python, numpy) • R is also open
  5. Python is used for: • Scientific computing • Enterprise software

    • Web design • Back-end • Front-end • Everything in between
  6. Python is used at Google, Rackspace, Microsoft, Intel, Walt Disney,

    MailChimp, twilio, Bank of America, Facebook, Instagram, HP, Linkedin, Elastic, Mozilla, YouTube, ILM, Thawte, CERN, Yahoo!, NASA, Trac, Civilization IV, reddit, LucasFilms, D-Link, Phillips, AstraZeneca, KLA-Tencor, Nerua https://us.pycon.org/2016/sponsors/ https://www.python.org/about/quotes/ https://en.wikipedia.org/wiki/Python_%28programming_language%29#Use https://en.wikipedia.org/wiki/List_of_Python_software https://www.python.org/about/success/
  7. Python is portable More or less same code runs on

    Windows, Linux, macOS, and any platform with a Python interpreter Python for "other" platforms
  8. Python is high level, easy to learn, and fast to

    develop So is MATLAB, Ruby, R…
  9. Python is fast enough Written in C Easy to wrap

    more C Easy to parallelize Benchmark Game | NumFocus Benchmarks
  10. Popularity • Python is 7th most popular tag on StackOverflow

    • R is 24th most popular tag • MATLAB is 58th most popular tag stackoverflow.com/tags, Feb 2017
  11. Active community • 3rd most active repositories on GitHub after

    Java (incl. Android) and JavaScript (incl. node.js) • ~4.8-fold more than R (12th) • ~27-fold more than MATLAB (24th) • As of Feb 2017 • See breakdown at githut
  12. Many new libraries released every month During 2016 >2,000 new

    packages released every month. See more stats at PyGarden/stats.
  13. Python can do nearly everything MATLAB and R can do

    With libraries like NumPy, SciPy, Matplotlib, IPython/Jupyter, Scikit-image, Scikit-learn, and more
  14. Deep learning with tensorflow.Session() as s: readout = s.graph.get_tensor_by_name('softmax:0') predictions

    = s.run(readout, {’Image': image_data}) pred_id = predictions.argmax() Label = node_lookup.id_to_string(pred_id) score = predictions[pred_id] basketball (score = 0.98201)
  15. First language at Israeli universities • TAU: CS & Engineering

    use Python • Technion: CS use C, some courses in Python • HUJI: CS & Humanities, use Python • BGU: CS use Java, Engineering use C
  16. History of Python • Developed in 1989-91 by Guido van

    Rossum in the Netherlands • Python 2.0 released Oct 2000 (support ends 2020) • Python 3.0 released Dec 2008 • Python 3.6 released Dec 2016 • Python 3 is widely used Guido
  17. 2 vs 3 If you use Python 2.x: • 2012

    called, they want their print back • Seriously, consider moving to 3.x ASAP • But at least 3.4 • See www.python3statement.org http://pygarden.com/stats
  18. • Formally this field is Population Genetics • Study changes

    in frequency of gene variants within populations • Focus on two forces: –Natural selection –Random genetic drift • Methods from applied math, statistics, CS, theoretical physics Population genetics 35 Scientific Python in action: Theoretical Evolution
  19. Wright-Fisher Model Standard model for change in frequency of gene

    variants 39 R.A. Fisher 1890-1962 UK & Australia Sewall Wright 1889-1988 USA
  20. Wright-Fisher Model Standard model for change in frequency of gene

    variants Two gene variants: 0 and 1. Number of individuals with each variant is n0 and n1 . Total population size is N = n0 + n1 . Frequency of each variant is p0 =n0 /N and p1 =n1 /N. 40
  21. Wright-Fisher Model Assume that variant 1 is favored by selection

    due to better survival or reproduction. The frequency of variant 1 after the effect of selection natural (p1 ) is: p" = n" ⋅ 1 + s n) + n" ⋅ 1 + s s is a selection coefficient, representing how much variant 1 is favored over variant 0. 41
  22. Wright-Fisher Model Random genetic drift accounts for the effect of

    random sampling. Due to genetic drift, the number of individuals with variant 1 in the next generation (n’1 ) is: n" * ∼ Binomial(N, p" ) The Binomial distribution is the distribution of the number of successes in N independent trials with probability of success p1 . 42
  23. Fixation Probability Assume a single copy variant 1 in a

    population of size N. What is the probability that variant 1 will take over the population rather than go extinct? 43
  24. NumPy The fundamental package for scientific computing with Python: •

    N-dimensional arrays • Random number generators • Array functions • Broadcasting • Tools for integrating C/C++ and Fortran code • Linear algebra • Fourier transform numpy.org 44 Loosing your loops
  25. from numpy.random import binomial n1 = 1 while 0 <

    n1 < N: n0 = N - n1 p1 = n1*(1+s) / (n0 + n1*(1+s)) n1 = binomial(N, p1) fixation = n1 == N Random drift Natural Selection n" * ∼ Binomial(N, p" ) p" = n" ⋅ 1 + s n) + n" ⋅ 1 + s Import a binomial random number generator from NumPy
  26. from numpy.random import binomial n1 = 1 while 0 <

    n1 < N: n0 = N - n1 p1 = n1*(1+s) / (n0 + n1*(1+s)) n1 = binomial(N, p1) fixation = n1 == N Random drift Natural Selection n" * ∼ Binomial(N, p" ) p" = n" ⋅ 1 + s n) + n" ⋅ 1 + s Start with a single copy of variant 1
  27. from numpy.random import binomial n1 = 1 while 0 <

    n1 < N: n0 = N - n1 p1 = n1*(1+s) / (n0 + n1*(1+s)) n1 = binomial(N, p1) fixation = n1 == N Random drift Natural Selection n" * ∼ Binomial(N, p" ) p" = n" ⋅ 1 + s n) + n" ⋅ 1 + s Until number of individuals with variant 1 is 0 or N: extinction or fixation
  28. from numpy.random import binomial n1 = 1 while 0 <

    n1 < N: n0 = N - n1 p1 = n1*(1+s) / (n0 + n1*(1+s)) n1 = binomial(N, p1) fixation = n1 == N Random drift Natural Selection n" * ∼ Binomial(N, p" ) p" = n" ⋅ 1 + s n) + n" ⋅ 1 + s The frequency of variant 1 after selection is p1
  29. from numpy.random import binomial n1 = 1 while 0 <

    n1 < N: n0 = N - n1 p1 = n1*(1+s) / (n0 + n1*(1+s)) n1 = binomial(N, p1) fixation = n1 == N Random drift Natural Selection n" * ∼ Binomial(N, p" ) p" = n" ⋅ 1 + s n) + n" ⋅ 1 + s Due to genetic drift, the number of individuals with variant 1 in the next generation is n1
  30. from numpy.random import binomial n1 = 1 while 0 <

    n1 < N: n0 = N - n1 p1 = n1*(1+s) / (n0 + n1*(1+s)) n1 = binomial(N, p1) fixation = n1 == N Random drift Natural Selection n" * ∼ Binomial(N, p" ) p" = n" ⋅ 1 + s n) + n" ⋅ 1 + s Fixation: n1 equals N Extinction: n1 equals 0
  31. NumPy vs. Pure Python NumPy is useful for random number

    generation: n1 = binomial(N, p1) Pure Python version would replace this with: from random import random rands = (random() for _ in range(N)) n1 = sum(1 for r in rands if r < p1) random is a standard library module 52
  32. NumPy vs. Pure Python %timeit simulation(N=1000, s=0.1) %timeit simulation(N=1000000, s=0.01)

    Pure Python version: 100 loops, best of 3: 6.42 ms per loop 1 loop, best of 3: 528 ms per loop NumPy version: 10000 loops, best of 3: 150 µs per loop x42 faster 1000 loops, best of 3: 313 µs per loop x1680 faster! 53
  33. • Optimizing compiler • Declare the static type of variables

    • Makes writing C extensions for Python as easy as Python itself • Foreign function interface for invoking C/C++ routines http://cython.org 55
  34. def simulation(np.uint64_t N, np.float64_t s): cdef np.uint64_t n1 = 1

    cdef np.uint64_t n0 cdef np.float64_t p while 0 < n1 < N: n0 = N - n1 p1 = n1 * (1 + s) / (n0 + n1 * (1 + s)) n1 = np.random.binomial(N, p1) return n1 == N 56
  35. Cython %timeit simulation(N=1000, s=0.1) %timeit simulation(N=1000000, s=0.01) Cython vs. NumPy:

    10000 loops, best of 3: 87.8 µs per loop x2 faster 10000 loops, best of 3: 177 µs per loop x1.75 faster 57
  36. To approximate the fixation probability we need to run many

    simulations. Thousands. 58 In principle, the standard error of our approximation decreases with the square root of the number of simulations: SEM ∼ 1/ Death to the Stock Photo
  37. Multiple simulations: for loop fixations = [ simulation(N, s) for

    _ in range(1000) ] fixations [False, True, False, ..., False, False] sum(fixations) / len(fixations) 0.195 60
  38. Multiple simulations: for loop %%timeit fixations = [ simulation(N, s)

    for _ in range(1000) ] 1 loop, best of 3: 8.05 s per loop 61
  39. Multiple simulations: NumPy def simulation(N, s, repetitions): n1 = np.ones(repetitions)

    update = np.array([True] * repetitions) while update.any(): p1 = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p1[update]) update = (n1 > 0) & (n1 < N) return n1 == N 62 Initialize multiple simulations
  40. Multiple simulations: NumPy def simulation(N, s, repetitions): n1 = np.ones(repetitions)

    update = np.array([True] * repetitions) while update.any(): p1 = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p1[update]) update = (n1 > 0) & (n1 < N) return n1 == N 63 Natural selection: n1 is an array so operations are element-wise
  41. Multiple simulations: NumPy def simulation(N, s, repetitions): n1 = np.ones(repetitions)

    update = np.array([True] * repetitions) while update.any(): p1 = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p1[update]) update = (n1 > 0) & (n1 < N) return n1 == N 64 Genetic drift: p1 is an array so binomial(N, p1) draws from multiple distributions
  42. Multiple simulations: NumPy def simulation(N, s, repetitions): n1 = np.ones(repetitions)

    update = np.array([True] * repetitions) while update.any(): p1 = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p1[update]) update = (n1 > 0) & (n1 < N) return n1 == N 65 update follows the simulations that didn’t finish yet
  43. Multiple simulations: NumPy def simulation(N, s, repetitions): n1 = np.ones(repetitions)

    update = np.array([True] * repetitions) while update.any(): p1 = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p1[update]) update = (n1 > 0) & (n1 < N) return n1 == N 66 update follows the simulations that didn’t finish yet
  44. Multiple simulations: NumPy def simulation(N, s, repetitions): n1 = np.ones(repetitions)

    update = np.array([True] * repetitions) while update.any(): p1 = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p1[update]) update = (n1 > 0) & (n1 < N) return n1 == N 67 result is array of Booleans: for each simulation, did variant 1 fix?
  45. Fixation probability as a function of N Nrange = np.logspace(1,

    6, 20, dtype=np.uint64) N must be an integer for this to evaluate to True: (n1 > 0) & (n1 < N) 69
  46. Fixation probability as a function of N fixations = [

    simulation( N, s, repetitions ) for N in Nrange ]
  47. Fixation probability as a function of N fixations = np.array(fixations)

    fixations array([[False, False, ..., False, False], [False, True, ..., False, False], , ..., [False, False, ..., True, False], [False, False, ..., False, False]], dtype=bool) 71
  48. Fixation probability as a function of N fixations = np.array(fixations)

    mean = fixations.mean(axis=1) sem = fixations.std( axis=1, ddof=1 ) / np.sqrt(repetitions) 72
  49. Approximation Kimura’s equation: e=>? − 1 e=>A? − 1 def

    kimura(N, s): return np.expm1(-2 * s) / np.expm1(-2 * N * s) expm1(x) is ex-1 with better precision for small values of x 73 Motoo Kimura 1994 - 1924 Japan & USA
  50. kimura works on arrays out-of-the-box %timeit [kimura(N=N, s=s) for N

    in Nrange] %timeit kimura(N=Nrange, s=s) 1 loop, best of 3: 752 ms per loop 1000 loops, best of 3: 3.91 ms per loop X200 faster! 74
  51. Numexpr Fast evaluation of element-wise array expressions using a vector-based

    virtual machine def kimura(N, s): return numexpr.evaluate( "expm1(-2 * s) / expm1(-2 * N * s)“) %timeit kimura(N=Nrange, s=s) 1000 loops, best of 3: 803 µs per loop x5 faster 75 https://github.com/pydata/numexpr
  52. Fixation time How much time does it take variant 1

    to go extinct or to fix? We want to keep track of time*. 77 time is measured in number of generations https://commons.wikimedia.org/wiki/File:Prim_clockwork.jpg
  53. def simulation(N, s, repetitions): n1 = np.ones(repetitions) T = np.empty_like(n1)

    update = (n1 > 0) & (n1 < N) t = 0 while update.any(): t += 1 p = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p[update]) T[update] = t update = (n1 > 0) & (n1 < N) return n1 == N, T 78 t keeps track of time
  54. def simulation(N, s, repetitions): n1 = np.ones(repetitions) T = np.empty_like(n1)

    update = (n1 > 0) & (n1 < N) t = 0 while update.any(): t += 1 p = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p[update]) T[update] = t update = (n1 > 0) & (n1 < N) return n1 == N, T 79 T holds time for extinction/fixation
  55. def simulation(N, s, repetitions): n1 = np.ones(repetitions) T = np.empty_like(n1)

    update = (n1 > 0) & (n1 < N) t = 0 while update.any(): t += 1 p = n1 * (1 + s) / (N + n1 * s) n1[update] = binomial(N, p[update]) T[update] = t update = (n1 > 0) & (n1 < N) return n1 == N, T 80 Return both Booleans and times (T)
  56. • Visualization library based on matplotlib and Pandas • High-level

    interface for attractive statistical graphics By Michael Waskom, postdoc at NYU http://seaborn.pydata.org 81 Statistical data visualization with Seaborn
  57. Plot with Seaborn from seaborn import distplot fixations, times =

    simulation(…) distplot(times[fixations]) distplot(times[~fixations]) 83
  58. Plot with Seaborn from seaborn import distplot fixations, times =

    simulation(…) distplot(times[fixations]) distplot(times[~fixations]) 84
  59. Plot with Seaborn from seaborn import distplot fixations, times =

    simulation(…) distplot(times[fixations]) distplot(times[~fixations]) 85
  60. from functools import partial from scipy.integrate import quad def integral(f,

    N, s, a, b): f = partial(f, N, s) return quad(f, a, b)[0] 89 integral will calculate ∫ , ,
  61. from functools import partial from scipy.integrate import quad def integral(f,

    N, s, a, b): f = partial(f, N, s) return quad(f, a, b)[0] 90 partial freezes N and s in f(N, s, x) to create f(x)
  62. from functools import partial from scipy.integrate import quad def integral(f,

    N, s, a, b): f = partial(f, N, s) return quad(f, a, b)[0] 91 SciPy’s quad computes a definite integral ∫ (using a technique from the Fortran library QUADPACK)
  63. def I1(N, s, x): … def I2(N, s, x): …

    92 I1 and I2 are defined according to the equations
  64. T_kimura is the fixation time given a single copy of

    variant 1: frequency x=1/N @np.vectorize def T_kimura(N, s): x = 1.0 / N J1 = -1.0 / (s * expm1(-2 * N * s)) * integral(I1, N, s, x, 1) J2 = -1.0 / (s * expm1(-2 * N *s)) * integral(I2, N, s, 0, x) u = expm1(-2 * N * s * x) / expm1(-2 * N * s) return J1 + ((1 - u) / u) * J2 93
  65. @np.vectorize def T_kimura(N, s): x = 1.0 / N J1

    = -1.0 / (s * expm1(-2 * N * s)) * integral(I1, N, s, x, 1) J2 = -1.0 / (s * expm1(-2 * N *s)) * integral(I2, N, s, 0, x) u = expm1(-2 * N * s * x) / expm1(-2 * N * s) return J1 + ((1 - u) / u) * J2 94 J1 and J2 are calculated using integrals of I1 and I2
  66. @np.vectorize def T_kimura(N, s): x = 1.0 / N J1

    = -1.0 / (s * expm1(-2 * N * s)) * integral(I1, N, s, x, 1) J2 = -1.0 / (s * expm1(-2 * N *s)) * integral(I2, N, s, 0, x) u = expm1(-2 * N * s * x) / expm1(-2 * N * s) return J1 + ((1 - u) / u) * J2 95 Tfix is the return value
  67. np.vectorize creates a function that takes a sequence and returns

    an array - x2 faster @np.vectorize def T_kimura(N, s): x = 1.0 / N J1 = -1.0 / (s * expm1(-2 * N * s)) * integral(I1, N, s, x, 1) J2 = -1.0 / (s * expm1(-2 * N *s)) * integral(I2, N, s, 0, x) u = expm1(-2 * N * s * x) / expm1(-2 * N * s) return J1 + ((1 - u) / u) * J2 96
  68. 97

  69. Dig Deeper Online Jupyter notebook: github.com/yoavram/PyConIL2016 Multi-type simulation: Includes L

    variants, with mutation. Follow n0 , n1 , …, nL until nL =N. 98
  70. • Numba: JIT compiler, array-oriented and math-heavy Python syntax to

    machine code • IPyParallel: IPython’s sophisticated and powerful architecture for parallel and distributed computing. • IPyWidgets: Interactive HTML Widgets for Jupyter notebooks and the IPython kernel 99 Dig Deeper Online Jupyter notebook: github.com/yoavram/PyConIL2016
  71. Thank You! Presentation & Jupyter notebook: https://github.com/yoavram/DataTalks2017 100 [email protected] @yoavram

    github.com/yoavram www.yoavram.com python.yoavram.com Death to the Stock Photo