from IPython.display import YouTubeVideo
YouTubeVideo("y2R3FvS4xr4")

# A silly function implemented in Python

def func_python(N):
    d = 0.0
    for i in range(N):
        d += (i % 3 - 1) * i
    return d

# Use IPython timeit magic to time the execution
%timeit func_python(10000)

%%file func_fortran.f

      subroutine func_fort(n, d)
           integer, intent(in) :: n
           double precision, intent(out) :: d
           integer :: i
           d = 0
           do i = 0, n - 1
                d = d + (mod(i, 3) - 1) * i
           end do
      end subroutine func_fort


!f2py -c func_fortran.f -m func_fortran > /dev/null

from func_fortran import func_fort
%timeit func_fort(10000)

import numpy as np
x = np.random.random(4)
print x
print x + 1  # add 1 to each element of x

print x * 2  # multiply each element of x by 2

print x * x  # multiply each element of x by itself

print x[1:] - x[:-1]

-x

np.sin(x)

x = np.random.random(10000)

%%timeit
# compute element-wise x + 1 via a ufunc 
y = np.zeros_like(x)
y = x + 1

%%timeit
# compute element-wise x + 1 via a loop
y = np.zeros_like(x)
for i in range(len(x)):
    y[i] = x[i] + 1

x = np.random.random(5)
print x
print np.minimum(x, 0.5)
print np.maximum(x, 0.5)

# contrast this behavior with that of min() and max()
print np.min(x)
print np.max(x)

%matplotlib inline
# On older IPython versions, use %pylab inline
import matplotlib.pyplot as plt

x = np.linspace(0, 10, 1000)
plt.plot(x, np.sin(x))

from scipy.special import gammaln
plt.plot(x, gammaln(x))

x = np.arange(5)
y = np.arange(1, 6)
np.add(x, y)

np.add.accumulate(x)

np.multiply.accumulate(x)

np.multiply.accumulate(y)

np.add.identity

np.multiply.identity

np.add.outer(x, y)

# make a times-table
x = np.arange(1, 13)
print np.multiply.outer(x, x)

# 1. computing the element-wise sine + cosine

from math import sin, cos
def slow_sincos(x):
    """x is a 1-dimensional array"""
    y = np.zeros_like(x)
    for i in range(len(x)):
        y[i] = sin(x[i]) + cos(x[i])
    return y

x = np.random.random(5)
print slow_sincos(x)

# write a fast_sincos function


# 2. computing the difference between adjacent squares

def slow_sqdiff(x):
    """x is a 1-dimensional array"""
    y = np.zeros(len(x) - 1)
    for i in range(len(y)):
        y[i] = x[i + 1] ** 2 - x[i] ** 2
    return y

x = np.random.random(5)
print slow_sqdiff(x)

# write a fast_sqdiff function


# 3. computing the outer-product of each consecutive pair

def slow_pairprod(x):
    """x is a 1-dimensional array"""
    if len(x) % 2 != 0:
        raise ValueError("length of x must be even")
    N = len(x) / 2
    y = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            y[i, j] = x[2 * i] * x[2 * j + 1]
    return y

x = np.arange(1, 9)
print slow_pairprod(x)

# write a fast_pairprod function


# 10 x 10 array drawn from a standard normal
x = np.random.randn(10, 10)

print x.mean()

print np.mean(x)

print x.std()

print x.var()

print x.sum()

print x.prod()

print np.median(x)

print np.percentile(x, 50)

print np.percentile(x, (25, 75))

x = np.random.rand(3, 5)
print x

print x.sum(0)  # sum along rows

print x.sum(1)  # sum along columns

print np.median(x, 1)

print np.mean(x, 1)

np.sum(x, 1)

np.add.reduce(x, 1)

np.prod(x, 1)

np.multiply.reduce(x, 1)

np.divide.reduce(x, 1)

np.add.reduce(x)

np.sum(x)

x = np.random.random(10000)

%timeit np.sum(x)
%timeit sum(x)

def slow_cubesum(x):
    """x is a 1D array"""
    result = 0
    for i in range(len(x)):
        result += x[i] ** 3
    return result

x = np.random.random(100)
print slow_cubesum(x)

# implement fast_cubesum



def slow_rms(x):
    """x is a 1D array"""
    m = np.mean(x)
    rms = 0
    for i in range(len(x)):
        rms += (x[i] - m) ** 2
    rms /= len(x)
    return np.sqrt(rms)

x = np.random.random(100)
print slow_rms(x)

# implement fast_rms


def slow_sillyfunc(N):
    """N is an integer"""
    d = 0.0
    for i in range(N):
        d += (i % 3 - 1) * i
    return d

print slow_sillyfunc(100)

# Implement fast_sillyfunc using ufuncs & aggragates



x = np.arange(10)
print x ** 2

Y = x * x[:, np.newaxis]
print Y

print Y + 10 * x

print Y + 10 * x[:, np.newaxis]

Y = np.random.random((2, 3, 4))
x = 10 * np.arange(3)

print Y + x[:, np.newaxis]

X = np.random.random((1000, 5))  # 1000 points in 5 dimensions

# 1. Compute the mean of the 1000 points in X



# 2. Compute the mean using np.add



# 3. Compute the standard deviation across the 1000 points



# 4. Compute the standard deviation using np.add only



# 5. Compute the whitened version of the array



x = np.array([1, 2, 3, -999, 2, 4, -999])

for i in range(len(x)):
    if x[i] < 0:
        x[i] = 0
print x

x = np.array([1, 2, 3, -999, 2, 4, -999])

mask = (x < 0)
print mask

x[mask] = 0
print x

x = np.array([1, 2, 3, -999, 2, 4, -999])
x[x < 0] = 0
print x

x = np.random.random(5)
print x

x[x > 0.5] = np.nan
print x

x[np.isnan(x)] = np.inf
print x

np.nan == np.nan

x[np.isinf(x)] = 0
print x

x = np.array([1, 0, -np.inf, np.inf, np.nan])
print "input   ", x
print "x < 0   ", (x < 0)
print "x > 0   ", (x > 0)
print "isinf   ", np.isinf(x)
print "isnan   ", np.isnan(x)
print "isposinf", np.isposinf(x)
print "isneginf", np.isneginf(x)

x = np.arange(16).reshape((4, 4))
print x

print (x < 5)

print ~(x < 5)

print (x < 10) & (x % 2 == 0)

print (x > 3) & (x < 8)

x = np.random.random(100)
print "array is length", len(x), "and has"
print (x > 0.5).sum(), "elements are greater than 0.5"

# clip is a useful function:
x = np.clip(x, 0.3, 0.6)

print np.sum(x < 0.3)
print np.sum(x > 0.6)

# works for 2D arrays as well
X = np.random.random((10, 10))
print (X < 0.1).sum()

x = np.random.random((3, 3))
print x

print np.where(x < 0.3)

print x[x < 0.3]

print x[np.where(x < 0.3)]

X = np.arange(16).reshape((4, 4))
print X

X[(0, 1), (1, 0)]

X[range(4), range(4)]

X.diagonal()

X.diagonal() = 100

X[range(4), range(4)] = 100

print X

X = np.arange(24).reshape((6, 4))
print X

i = np.arange(6)
np.random.shuffle(i)
print i

print X[i]  # X[i, :] is identical

i2 = i.reshape(3, 2)
print X[i2]

print X[i2].shape