%matplotlib inline
import numpy as np
a = np.array([1, 2, 3])
a

c = np.array([4, 5, 6])
c

a*c

_

np.dot(a, c)

A = np.outer(c, a)
A

a * a

b = a ** 2
b

am = np.matrix(a)
cm = np.matrix(c)
am * cm.T

g = np.ones( (3, 3, 3) )
g

#h = np.matrix(g)
print "Error has a problem???"

np.exp(a)

np.log(_)

np.sqrt(a)

np.set_printoptions(precision=4)
print (np.sqrt(a))

help(np.set_printoptions)

np.set_printoptions(precision=8)
np.sqrt(a)

2 ** -24

b.sum()

c.mean()

dir(a)

help(a.argmax)

np.pi

_

y = np.tan(np.pi/6)
y

B = np.c_[[-3, 2, -1], [0, 5, 4], [-1, -7, 8]]
B

import numpy.linalg as nl
x = nl.solve(B, c)
x

nl.norm(np.dot(B,x) - c) / (nl.norm(B) * nl.norm(x))

np.set_printoptions(precision=5)  # make it the m*lab default
w, _ = nl.eig(B)
w

w, v = nl.eig(B)
v

v = np.arange(6)
v

w = np.arange(2, 10, 3)
w

y = np.arange(1, 0, -0.25)
y

y = np.arange(1, -0.001, -0.25)
y

C = np.c_[A, [8, 9, 10]]
C

D = np.r_[B, B, B]
D

E = np.vstack((B, B, a))
E

C[1, 2]

C[1:3, 0:2]

I3 = np.eye(3)
I3

Y = np.zeros((3,5))
Y

Z = np.ones((2,2))
Z

rn = np.random.RandomState()  # initialize a new RandomState object
rn.seed(20)
F = rn.rand(3, 3)
F

G = rn.randn(1, 5)
G

%who

%whos

g = 2.
for k in range(10):
    g = 1 + 1 / g
g

import scipy.constants as sconst
sconst.golden

#import mpld3
#mpld3.enable_notebook()
import matplotlib.pylab as plt
t = np.arange(0, 1.005, 0.005)
z = np.exp(10 * t * (t - 1)) * np.sin(12 * np.pi * t)
plt.plot(t, z)
plt.title("Figure 1.2. Basic 2D picture produced by plt.plot")

plt.hist(np.random.randn(1000))
plt.title("Figure 1.3. Histogram produced by plt.hist.")

%reset -f
import numpy as np
import pylab as plt
rd = np.random.RandomState()
rd.seed(100)
x = np.zeros(1000)
x[0:2] = 1, 2
for n in range(1, 999):
    x[n + 1] = x[n] + np.sign(rd.rand(1) - 0.5) * x[n - 1]
    
xx = np.arange(1, 1001)
plt.semilogy(xx, np.abs(x))
c = 1.13198824
plt.hold(True)
plt.semilogy(xx, c ** xx)
plt.title("Figure 1.4. Growth of a random Fibonacci sequence")
plt.hold(False)

x = 3
if x % 2 == 0:
    print ('x is even')
else:
    print ('x is odd')

#COLLATZ Collatz iteration. 

#n = int(raw_input(’Enter an integer bigger than 2: ’))
init_n = n = 27
narray = np.zeros(1000)  # only a maximum of 1000 iterations
narray[0] = n
count = 0
while n != 1:
    if n % 2 == 1:  # Remainder modulo 2. 
        n = 3 * n + 1
    else:
        n = n / 2
    count += 1
    narray[count] = n  # Store the current iterate

# Plot with * marker and solid line style.
# Only plot the non zero entries in narray
plt.plot(narray[narray != 0], '*-')
plt.title('Figure 1.5. Collatz iteration starting at %d' % init_n)

#COLLBAR Collatz iteration bar graph.
N = 29 # Use starting values 1,2,...,N.
niter = np.zeros(N);  # Preallocate array.
for i in range(N):
    count = 0
    n = i + 1
    while n != 1:
        if n % 2 == 1:
            n = 3 * n + 1
        else:
            n = n / 2

        count += 1
    niter[i] = count

left = np.arange(29)
plt.bar(left, niter) # Bar graph.
plt.grid()   # Add horizontal and vertical grid lines.
plt.title('Figure 1.6. Col1atz iteration counts')
plt.xlabel('Starting value') # Label X-axis.
plt.ylabel('Number of iterations')  #Label y-axis.

#MANDEL Mandelbrot set.

x = np.linspace(-2.1, 0.6, 301)
y = np.linspace(-1.1, 1.1, 301)
[X,Y] = np.meshgrid(x, y)
C = X + 1j * Y

Z_max = 1e6
it_max = 50
Z = C

for k in range(it_max):
    Z = Z ** 2 + C

plt.contourf(x, y, np.float64(np.abs(Z) < Z_max))
plt.title("Figure 1.7  Mandelbrot Set")

import scipy.integrate as si

def lorenzde(y, t):
    '''LORENZDE Lorenz equations.
        YPRIME = LORENZDE(Y,T)
    '''
    yprime = np.array([10. * (y[1] - y[0]), 28. * y[0] - y[1] - y[0] * y[2], 
                       y[0] * y[1] - 8. * y[2] / 3.])
    return yprime

#lrun    ODE solving example:  Lortez.

t = np.arange(0, 50.01, .01)  # time points on which to solve
yzero = np.array([0., 1., 0.])
print (len(yzero))
y = si.odeint(lorenzde, yzero, t)

plt.plot(y[:, 0], y[:, 2])
plt.xlabel('y_0')
plt.ylabel('y_2')
plt.title('Figure 1.8  Lorenz equations')

def gasket(pa, pb, pc, level):
    '''
    GASKET Recursively generated Sierpinski gasket. 
    GASKET(PA, PB, PC, LEVEL) generates an approximation to 
    the Sierpinski gasket, where the 2-vectors PA, PB, and PC Z define the triangle vertices. 
    LEVEL is the level of recursion.
    '''
    if level == 0:
        # Fill the triangle with vertices Pa, Pb, Pc. 
        plt.fill([pa[0], pb[0], pc[0]], [pa[1], pb[1], pc[1]], 'g') 
        plt.hold(True)
    else:
        # Recursive calls for the three subtriangles. 
        gasket(pa, (pa + pb) / 2., (pa + pc) / 2., level - 1) 
        gasket(pb, (pb + pa) / 2., (pb + pc) / 2., level - 1) 
        gasket(pc, (pc + pa) / 2., (pc + pb) / 2., level - 1)
        
pa = np.array([0, 0])
pb = np.array([1, 0])
pc = np.array([0.5, np.sqrt(3)/2.])
level = 1
gasket(pa, pb, pc, level)
plt.hold(False)
plt.title("Figure 1.9  Gasket level = 1")
plt.axis('equal')

level = 5
gasket(pa, pb, pc, level)
plt.hold(False)
plt.title("Figure 1.10  Gasket level = 5")
plt.axis('equal')

#BARNSLEY Barnsley’s game to compute Sierpinski gasket. 

rd = np.random.RandomState()  # Initialize a new RandomState object
rd.seed(1)                    # and (re)set the seed

V = np.c_[[0, 0], [1, 0], [0.5, np.sqrt(3)/2]] # Columns give triangle vertices. 
point = V[:, 0]                    # Start at a Vertex. 

n = 1000     #  Change this number to experiment

for k in range(n):
    node = int(np.ceil(3 * np.random.rand()) - 1) # node is 0, 1, or 2 with equal prob. 
    point = (V[:, node] + point)/2; 
    plt.plot(point[0], point[1], ".", markersize=15) 
    plt.hold(True)

plt.axis('equal') 
plt.hold(False)

#  SWEEP   Generates a volume-swept 3D object.
import numpy as np
import matplotlib.pylab as pld
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

fig = plt.figure()
ax = fig.gca(projection='3d')

n = 10      # Number of increments - try increasing

zz = np.linspace(-5, 5, n).reshape(n, 1)
radius = np.sqrt(1 + zz ** 2)    # Try changing sqrt to cos, sin, log or abs
theta = 2 * np.pi * np.linspace(0, 1, n)
x = radius * np.cos(theta)
y = radius * np.sin(theta)
z = zz[:, n * [0]]   #  Tony's trick!   Who is Tony?  No idea.

surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=cm.coolwarm,
        linewidth=.2, antialiased=True)

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.title('Figure 1.12  3D picture produced by sweep')