# Scientific computing

- Numpy: advanced array operations, multidimensional arrays
- Scipy: scientific computing by examples
    - singular value decomposition, with scipy.linalg
    - scipy.signal and scipy.fftpack: Signal theory
    - scipy.optimize: Local and global optimization, fitting and root finding
    - scipy.interpolate: Cubic interpolation
    - scipy.integrate: Integration and ODE solvers
    - scipy.ndimage - Image processing
- Simpy: symbolic math

## Numpy

- using numpy improves RAM space and speed
- numpy enforces strong typing, while Python is a dynamic typed language
- translates in Numpy using less heap space for representing data
- multidimensional array operations are the core of scientific computing

Further reading:
- https://docs.scipy.org/doc/numpy/user/basics.html


```python
import numpy as np
L = range(1000)
%timeit [i**2 for i in L]
a = np.arange(1000)
%timeit a**2
print(type(L[1]))
print(a.dtype)
```

    293 µs ± 1.46 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
    1.46 µs ± 19.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
    <class 'int'>
    int32



```python
import numpy as np
##Get help!
#np.lookfor('create array')
#np.array?
#np.arr*?
```


```python
a = np.array([0, 1, 2, 3])
print("a:\n", a)
b = np.array([[0, 1, 2], [3, 4, 5]])
print("b:\n", b)

print("b.shape:\n", b.shape)
```

    a:
     [0 1 2 3]
    b:
     [[0 1 2]
     [3 4 5]]
    b.shape:
     (2, 3)


#### Data types

There are 5 basic numerical types representing booleans (bool), integers (int), unsigned integers (uint) floating point (float) and complex. Those with numbers in their name indicate the bitsize of the type (i.e. how many bits are needed to represent a single value in memory). 


```python
print("np.float32(1.0) :", np.float32(1.0))
print("np.arange(3, dtype=np.uint8) :", np.arange(3, dtype=np.uint8))
```

    np.float32(1.0) : 1.0
    np.arange(3, dtype=np.uint8) : [0 1 2]



```python
z = np.array([1, 2, 3], dtype='f')
print(z)
```

    [1. 2. 3.]



```python
z = np.arange(3, dtype=np.uint8)
print(z)
print(z.astype(float))
print(z.dtype)
```

    [0 1 2]
    [0. 1. 2.]
    uint8


#### Array creation


```python
# extrinsic
x = np.array([2,3,1,0])
print(x)
print()
x = np.array([[ 1., 2.], [ 0., 0.], [ 1., 3.]])
print(x)
```

    [2 3 1 0]
    
    [[1. 2.]
     [0. 0.]
     [1. 3.]]



```python
#intrinsic
b = np.arange(1, 9, 2)
print(b)
c = np.linspace(0, 1, 6)
print(c)
```

    [1 3 5 7]
    [0.  0.2 0.4 0.6 0.8 1. ]



```python
print(np.arange(35).reshape(5,7))
```

    [[ 0  1  2  3  4  5  6]
     [ 7  8  9 10 11 12 13]
     [14 15 16 17 18 19 20]
     [21 22 23 24 25 26 27]
     [28 29 30 31 32 33 34]]



```python
x = np.random.rand(35).reshape(5,7)
print(x)
```

    [[0.9670321  0.41898259 0.20363724 0.35449615 0.16035759 0.71796865
      0.91560078]
     [0.7941912  0.11327121 0.39895543 0.1938332  0.25015112 0.7301424
      0.8902718 ]
     [0.38559724 0.80936697 0.87522364 0.2348649  0.50799691 0.96275272
      0.25653361]
     [0.11539546 0.29497229 0.82328139 0.47741048 0.08609102 0.39704512
      0.5578069 ]
     [0.72175903 0.77159655 0.4940947  0.16196267 0.405733   0.1140635
      0.94903288]]



```python
%pylab inline
import matplotlib.pyplot as plt

image = np.random.rand(30, 30)
plt.imshow(image, cmap=plt.cm.hot)    
plt.colorbar()    
```

#### Indexing, slicing and selection


```python
a = np.arange(10)
print(a)
print(a[0], a[2], a[-1], a[-3])
print(a[2:5], a[2:], a[:-2], a[::2], a[2::2])
```

    [0 1 2 3 4 5 6 7 8 9]
    0 2 9 7
    [2 3 4] [2 3 4 5 6 7 8 9] [0 1 2 3 4 5 6 7] [0 2 4 6 8] [2 4 6 8]



```python
a = np.diag(np.arange(3))
a[2, 1] = 10 # !third line, !second column

print(a)
print(a[1, 1])
#print(a[1])
#print(a[:,1], a[1,:])
#print(a[1:,1:])
```

    [[ 0  0  0]
     [ 0  1  0]
     [ 0 10  2]]
    1



```python
# array indexes
x = np.arange(10,1,-1)
print(x)
print()
print(x[np.array([3,3,-3,8])])
print()
print(x[np.array([[1,1],[2,3]])])
```

    [10  9  8  7  6  5  4  3  2]
    
    [7 7 4 2]
    
    [[9 9]
     [8 7]]



```python
# 10 random numbers 0 - 20
a = np.random.randint(0, 20, 10)
print(a)
print(a%3==0)
print(a[a%3==0])

a[a % 3 == 0] = -1
print(a)
```

    [18 10  3 15 12  1  6  2  3  4]
    [ True False  True  True  True False  True False  True False]
    [18  3 15 12  6  3]
    [-1 10 -1 -1 -1  1 -1  2 -1  4]


Task:
- What does this do:


```python
# How does it work?
# Print the primes!
def get_primes():
    primes = np.ones((100,), dtype=bool)
    primes[:2] = 0
    N_max = int(np.sqrt(len(primes)))
    for j in range(2, N_max):
        primes[2*j::j] = 0
    return primes
print(get_primes())

```

    [False False  True  True False  True False  True False False False  True
     False  True False False False  True False  True False False False  True
     False False False False False  True False  True False False False False
     False  True False False False  True False  True False False False  True
     False False False False False  True False False False False False  True
     False  True False False False False False  True False False False  True
     False  True False False False False False  True False False False  True
     False False False False False  True False False False False False False
     False  True False False]


#### Broadcasting, assignment, structured arrays


```python
a = np.arange(10)
b = a
print(np.may_share_memory(a, b))

a = np.arange(10)
c = a.copy()  # force a copy
print(np.may_share_memory(a, c))
```

    True
    False



```python
#Array operations
a = np.array([1, 2, 3, 4])
print("a: ", a)
print("a + 1, 2**a: ", a + 1, 2**a)
print ("2**(a + 1) - a: ", 2**(a + 1) - a)
```

    a:  [1 2 3 4]
    a + 1, 2**a:  [2 3 4 5] [ 2  4  8 16]
    2**(a + 1) - a:  [ 3  6 13 28]



```python
a = np.array([1, 2, 3, 4])
b = np.ones(4)+1
print("a: ",a)
print("b: ",b)
print("a - b, a * b: ", a - b, a * b)
```

    a:  [1 2 3 4]
    b:  [2. 2. 2. 2.]
    a - b, a * b:  [-1.  0.  1.  2.] [2. 4. 6. 8.]



```python
c = np.ones((3, 3))
print(c)
print(2*c + 1)
```

    [[1. 1. 1.]
     [1. 1. 1.]
     [1. 1. 1.]]
    [[3. 3. 3.]
     [3. 3. 3.]
     [3. 3. 3.]]



```python
# matrix multiplication
a = np.ones((3, 2)) + 1
b = np.ones((2, 3)) + 1
c = a.dot(b)
print(a, b, c, sep="\n\n")
```

    [[2. 2.]
     [2. 2.]
     [2. 2.]]
    
    [[2. 2. 2.]
     [2. 2. 2.]]
    
    [[8. 8. 8.]
     [8. 8. 8.]
     [8. 8. 8.]]



```python
a = np.arange(5)
print(np.sin(a))
print(np.log(a))
print(np.exp(a))
```

    [ 0.          0.84147098  0.90929743  0.14112001 -0.7568025 ]
    [      -inf 0.         0.69314718 1.09861229 1.38629436]
    [ 1.          2.71828183  7.3890561  20.08553692 54.59815003]


    D:\programs\anaconda\envs\biopy37\lib\site-packages\ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in log
      This is separate from the ipykernel package so we can avoid doing imports until



```python
# shape manipulation
x = np.array([1, 2, 3])
y = x[:, np.newaxis]
z = x[np.newaxis, :]
print(x, y, z, sep='\n\n')
```

    [1 2 3]
    
    [[1]
     [2]
     [3]]
    
    [[1 2 3]]



```python
# flatten
a = np.array([[1, 2, 3], [4, 5, 6]])
print(a)
print(a.ravel())
```

    [[1 2 3]
     [4 5 6]]
    [1 2 3 4 5 6]



```python
# sorting matrices
a = np.array([[4, 3, 5], [1, 2, 1]])
b = np.sort(a, axis=1) #sorting per row
print(a)
print(b)
```

    [[4 3 5]
     [1 2 1]]
    [[3 4 5]
     [1 1 2]]



```python
# sorting arguments
a = np.array([4, 3, 1, 2])
j = np.argsort(a)
print(j, a[j])
```

    [2 3 1 0] [1 2 3 4]


#### Reductions


```python
# unidimensional
x = np.array([1, 2, 3, 4])
print(np.sum(x), x.sum(), x.sum(axis=0))
```

    10 10 10



```python
# multidimensional
x = np.array([[1, 1], [2, 2]])
print(x)
print(x.sum(axis=0))   # rows (first dimension)
print(x.sum(axis=1))   # columns (second dimension)
```

    [[1 1]
     [2 2]]
    [3 3]
    [2 4]



```python
x = np.array([1, 3, 2])
print(x)
print(x.min(), x.max(), x.argmin(), x.argmax())

print(np.all([True, True, False]), np.any([True, True, False]))

x = np.array([1, 2, 3, 1])
y = np.array([[1, 2, 3], [5, 6, 1]])
print(x.mean(), x.std(), np.median(x), np.median(y, axis=-1))
```

    [1 3 2]
    1 3 0 1
    False True
    1.75 0.82915619758885 1.5 [2. 5.]



```python
a = np.zeros((100, 100))
print(np.any(a != 0))

a = np.array([1, 2, 3, 2])
b = np.array([2, 2, 3, 2])
c = np.array([6, 4, 4, 5])
print(((a <= b) & (b <= c)).all())
```

    False
    True


**Tricksy task:**
- Replace all values greater than 25 with 9 and all values smaller than 10 with 29.


```python
a = np.random.randint(0, 50, 10)
print(a)
```

    [ 3  9 33 19 46 34 39 39 32 27]


### Sympy

Symbolic math is sometimes important, especially if we are weak at calculus or if we need to perform automated calculus on long formulas. We are briefly going through a few test cases, to get the feel of it. Symbolic math is especially developed for [Mathematica](http://www.wolfram.com/mathematica/?source=nav), or [Sage](http://www.sagemath.org/) which is an open-source equivalent.


```python
import sympy
print sympy.sqrt(8)
import math
print math.sqrt(8)
```

    2*sqrt(2)
    2.82842712475



```python
from sympy import symbols
x, y, z, t = symbols('x y z t')
expr = x + 2*y
print expr
print x * expr
from sympy import expand, factor, simplify
expanded_expr = expand(x*expr)
print expanded_expr
print factor(expanded_expr)
exp = expanded_expr.subs(x, z**t)
print exp
print simplify(exp)
```

    x + 2*y
    x*(x + 2*y)
    x**2 + 2*x*y
    x*(x + 2*y)
    2*y*z**t + z**(2*t)
    z**t*(2*y + z**t)


In the scipy.optimize paragraph we needed the Hessian matrix for a function f. Here is how you can obtain it in sympy:


```python
import sympy
x, y = sympy.symbols('x y')
f = .5*(1 - x)**2 + (y - x**2)**2
h = sympy.hessian(f, [x,y])
print(h)
from IPython.display import Latex
Latex(sympy.latex(h))

```

    Matrix([[12*x**2 - 4*y + 1.0, -4*x], [-4*x, 2]])





\left[\begin{matrix}12 x^{2} - 4 y + 1.0 & - 4 x\\- 4 x & 2\end{matrix}\right]




```python
from IPython.display import HTML
HTML('')
```









```python

```

bc9ca0637d778699d0254fb57df3ed74526d8507
ipynb
Jupyter Notebook
day2/scicomp_numpy.ipynb
grokkaine/biopycourse
cb8b554abb987e6f657c5e522c7e28ecbc9fb4d5
CC0-1.0
2017-05-16T06:07:22.000Z
2021-08-06T14:58:28.000Z
day2/scicomp_numpy.ipynb
grokkaine/biopycourse
cb8b554abb987e6f657c5e522c7e28ecbc9fb4d5
CC0-1.0
day2/scicomp_numpy.ipynb
grokkaine/biopycourse
cb8b554abb987e6f657c5e522c7e28ecbc9fb4d5
CC0-1.0
2017-05-16T07:25:08.000Z
2021-04-22T19:22:53.000Z
Qwen/Qwen-72B
1. YES
2. YES