# This is a block of code, below you'll see its output
print "Welcome to the world of scientific computing with Python!"

def f(x):
    return (x-3)*(x-5)*(x-7)+85

import numpy as np
x = np.linspace(0, 10, 200)
y = f(x)

a, b = 1, 9
xint = x[logical_and(x>=a, x<=b)][::30]
yint = y[logical_and(x>=a, x<=b)][::30]

import matplotlib.pyplot as plt
plt.plot(x, y, lw=2)
plt.axis([0, 10, 0, 140])
plt.fill_between(xint, 0, yint, facecolor='gray', alpha=0.4)
plt.text(0.5 * (a + b), 30,r"$\int_a^b f(x)dx$", horizontalalignment='center', fontsize=20);

from scipy.integrate import quad, trapz
integral, error = quad(f, 1, 9)
trap_integral = trapz(yint, xint)
print "The integral is: %g +/- %.1e" % (integral, error)
print "The trapezoid approximation with", len(xint), "points is:", trap_integral
print "The absolute error is:", abs(integral - trap_integral)

import numpy

import numpy as np

lst = [10, 20, 30, 40]
arr = np.array([10, 20, 30, 40])

lst[0]

arr[0]

arr[-1]

arr[2:]

lst[-1] = 'a string inside a list'
lst

arr[-1] = 'a string inside an array'

arr.dtype

arr[-1] = 1.234
arr

np.zeros(5, float)

np.zeros(3, int)

np.zeros(3, complex)

print '5 ones:', np.ones(5)

a = empty(4)
a.fill(5.5)
a

np.arange(5)

print "A linear grid between 0 and 1:", np.linspace(0, 1, 5)
print "A logarithmic grid between 10**1 and 10**4: ", np.logspace(1, 4, 4)

np.random.randn(5)

norm10 = np.random.normal(10, 3, 5)
norm10

mask = norm10 > 9
mask

print 'Values above 9:', norm10[mask]

print 'Resetting all values above 9 to 0...'
norm10[mask] = 0
print norm10

lst2 = [[1, 2], [3, 4]]
arr2 = np.array([[1, 2], [3, 4]])
arr2

print lst2[0][1]
print arr2[0,1]

np.zeros((2,3))

np.random.normal(10, 3, (2, 4))

arr = np.arange(8).reshape(2,4)
print arr

print 'Slicing in the second row:', arr[1, 2:4]
print 'All rows, third column   :', arr[:, 2]

print 'First row:  ', arr[0]
print 'Second row: ', arr[1]

print 'Data type                :', arr.dtype
print 'Total number of elements :', arr.size
print 'Number of dimensions     :', arr.ndim
print 'Shape (dimensionality)   :', arr.shape
print 'Memory used (in bytes)   :', arr.nbytes

print 'Minimum and maximum             :', arr.min(), arr.max()
print 'Sum and product of all elements :', arr.sum(), arr.prod()
print 'Mean and standard deviation     :', arr.mean(), arr.std()

print 'For the following array:\n', arr
print 'The sum of elements along the rows is    :', arr.sum(axis=1)
print 'The sum of elements along the columns is :', arr.sum(axis=0)

np.zeros((3,4,5,6)).sum(2).shape

print 'Array:\n', arr
print 'Transpose:\n', arr.T

arr1 = np.arange(4)
arr2 = np.arange(10, 14)
print arr1, '+', arr2, '=', arr1+arr2

print arr1, '*', arr2, '=', arr1*arr2

arr1 + 1.5*np.ones(4)

arr1 + 1.5

b = np.array([2, 3, 4, 5])
print arr, '\n\n+', b , '\n----------------\n', arr + b

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

(c[:, np.newaxis]).shape

arr + c[:, np.newaxis]

x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x)

v1 = np.array([2, 3, 4])
v2 = np.array([1, 0, 1])
print v1, '.', v2, '=', v1.dot(v2)

A = np.arange(6).reshape(2, 3)
print A, 'x', v1, '=', A.dot(v1)

print A.dot(A.T)

print A.T.dot(A)

arr = np.arange(10).reshape(2, 5)
np.savetxt('test.out', arr, fmt='%.2e', header="My dataset")
!cat test.out

arr2 = np.loadtxt('test.out')
print arr2

np.save('test.npy', arr2)
# Now we read this back
arr2n = np.load('test.npy')
# Let's see if any element is non-zero in the difference.
# A value of True would be a problem.
print 'Any differences?', np.any(arr2-arr2n)

np.savez('test.npz', arr, arr2)
arrays = np.load('test.npz')
arrays.files

np.savez('test.npz', array1=arr, array2=arr2)
arrays = np.load('test.npz')
arrays.files

print 'First row of first array:', arrays['array1'][0]
# This is an equivalent way to get the same field
print 'First row of first array:', arrays.f.array1[0]

import matplotlib.pyplot as plt

x = np.linspace(0, 2*np.pi)
y = np.sin(x)
plt.plot(x,y, label='sin(x)')
plt.legend()
plt.grid()
plt.title('Harmonic')
plt.xlabel('x')
plt.ylabel('y');

plt.plot(x, y, linewidth=2);

plt.plot(x, y, 'o', markersize=5, color='r');

# example data
x = np.arange(0.1, 4, 0.5)
y = np.exp(-x)

# example variable error bar values
yerr = 0.1 + 0.2*np.sqrt(x)
xerr = 0.1 + yerr

# First illustrate basic pyplot interface, using defaults where possible.
plt.figure()
plt.errorbar(x, y, xerr=0.2, yerr=0.4)
plt.title("Simplest errorbars, 0.2 in x, 0.4 in y");

x = np.linspace(-5, 5)
y = np.exp(-x**2)
plt.semilogy(x, y);

mu, sigma = 100, 15
x = mu + sigma * np.random.randn(10000)

# the histogram of the data
n, bins, patches = plt.hist(x, 50, normed=1, facecolor='g', alpha=0.75)

plt.xlabel('Smarts')
plt.ylabel('Probability')
plt.title('Histogram of IQ')
# This will put a text fragment at the position given:
plt.text(55, .027, r'$\mu=100,\ \sigma=15$', fontsize=14)
plt.axis([40, 160, 0, 0.03])
plt.grid(True)

from matplotlib import cm
plt.imshow(np.random.rand(5, 10), cmap=cm.gray, interpolation='nearest');

img = plt.imread('stinkbug.png')
print 'Dimensions of the array img:', img.shape
plt.imshow(img);

from mpl_toolkits.mplot3d import Axes3D

from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,
        linewidth=0, antialiased=False)
ax.set_zlim3d(-1.01, 1.01);
% This cell is for the pdflatex output only
\begin{figure}[htbp]
\centering
\includegraphics[width=3in]{ipython_qtconsole2.png}
\caption{The IPython Qt console: a lightweight terminal for scientific exploration, with code, results and graphics in a soingle environment.}
\end{figure}% This cell is for the pdflatex output only
\begin{figure}[htbp]
\centering
\includegraphics[width=3in]{ipython-notebook-specgram-2.png}
\caption{The IPython Notebook: text, equations, code, results, graphics and other multimedia in an open format for scientific exploration and collaboration}
\end{figure}