def f(x):
    return x**3

n = 5                 # no of points
dx = 1.0/(n-1)        # x spacing in [0,1]
xlist = [i*dx for i in range(n)]
ylist = [f(x) for x in xlist]

n = 5                     # number of points
x = np.linspace(0, 1, n) # n points in [0, 1]
y = np.zeros(n)          # n zeros (float data type)
for i in xrange(n):
    y[i] = f(x[i])

from math import sin

for i in xrange(len(x)):
    y[i] = sin(x[i])

y = np.sin(x)                  # x: array, y: array

n = 100000

import numpy as np

x = np.linspace(0, 2*np.pi, n+1)

y = np.zeros(len(x))

%timeit for i in xrange(len(x)): \

%timeit y = np.sin(x)*np.exp(-x)

247/4.77

from numpy import sin, exp, linspace

def f(x):
    return x**3 + sin(x)*exp(-3*x)

x = 1.2                     # float object
y = f(x)                    # y is float

x = linspace(0, 3, 10001)   # 10000 intervals in [0,3]
y = f(x)                    # y is array

import math, numpy
x = numpy.linspace(0, 1, 11)
math.sin(x[3])

math.sin(x)

numpy.sin(x)

from numpy import *

n = 100
x = linspace(0, 4*pi, n+1)
y = 2.5 + x**2*exp(-0.5*x)*sin(x-pi/3)

%matplotlib inline

from scitools.std import *  # import numpy and plotting

# Make points along the curve
t = linspace(0, 3, 51)      # 50 intervals in [0, 3]
y = t**2*exp(-t**2)         # vectorized expression

plot(t, y)           # make plot on the screen
savefig('fig.pdf')   # make PDF image for reports
savefig('fig.png')   # make PNG image for web pages

from scitools.std import *  # import numpy and plotting

def f(t):
    return t**2*exp(-t**2)

t = linspace(0, 3, 51)    # t coordinates
y = f(t)                  # corresponding y values

plot(t, y)

xlabel('t')               # label on the x axis
ylabel('y')               # label on the y axix
legend('t^2*exp(-t^2)')   # mark the curve
axis([0, 3, -0.05, 0.6])  # [tmin, tmax, ymin, ymax]
title('My First Easyviz Demo')

plot(t, y,
     xlabel='t',
     ylabel='y',
     legend='t^2*exp(-t^2)',
     axis=[0, 3, -0.05, 0.6],
     title='My First Easyviz Demo',
     savefig='tmp1.png',
     show=True)   # display on the screen (default)

from scitools.std import *   # curve plotting + array computing

def f1(t):
    return t**2*exp(-t**2)

def f2(t):
    return t**2*f1(t)

t = linspace(0, 3, 51)
y1 = f1(t)
y2 = f2(t)

plot(t, y1)
hold('on')   # continue plotting in the same plot
plot(t, y2)

xlabel('t')
ylabel('y')
legend('t^2*exp(-t^2)', 't^4*exp(-t^2)')
title('Plotting two curves in the same plot')
savefig('tmp2.png')

plot(t, y1, t, y2,
     xlabel='t', ylabel='y',
     legend=('t^2*exp(-t^2)', 't^4*exp(-t^2)'),
     title='Plotting two curves in the same plot',
     savefig='tmp2.pdf')

# equivalent to
plot(t, y1)
hold('on')
plot(t, y2)

xlabel('t')
ylabel('y')
legend('t^2*exp(-t^2)', 't^4*exp(-t^2)')
title('Plotting two curves in the same plot')
savefig('tmp2.pdf')

plot(t, y1, 'r-')   # red (r) line (-)
hold('on')
plot(t, y2, 'bo')   # blue (b) circles (o)

# or
plot(t, y1, 'r-', t, y2, 'bo')

t = linspace(0, 3, 51)
plot(t, t**2*exp(-t**2), t, t**4*exp(-t**2))

from scitools.std import *
import time

def f(x, m, s):
    return (1.0/(sqrt(2*pi)*s))*exp(-0.5*((x-m)/s)**2)

m = 0;  s_start = 2;  s_stop = 0.2
s_values = linspace(s_start, s_stop, 30)

x = linspace(m -3*s_start, m + 3*s_start, 1000)
# f is max for x=m (smaller s gives larger max value)
max_f = f(m, m, s_stop)

# Show the movie on the screen
# and make hardcopies of frames simultaneously
import time
frame_counter = 0

for s in s_values:
    y = f(x, m, s)
    plot(x, y, axis=[x[0], x[-1], -0.1, max_f],
         xlabel='x', ylabel='f', legend='s=%4.2f' % s,
         savefig='tmp_%04d.png' % frame_counter)
    frame_counter += 1
    #time.sleep(0.2)  # pause to control movie speed

def H(x):
    return (0 if x < 0 else 1)

x = linspace(-10, 10, 5)  # few points (simple curve)
y = H(x)
plot(x, y)

x = linspace(-10,10,5)
x

b = x < 0
b

bool(b)  # evaluate b in a boolean context

def H_loop(x):
    r = zeros(len(x))  # or r = x.copy()
    for i in xrange(len(x)):
        r[i] = H(x[i])
    return r

n = 5
x = linspace(-5, 5, n+1)
y = H_loop(x)

from numpy import vectorize

# Automatic vectorization of function H
Hv = vectorize(H)
# Hv(x) works with array x

def Hv(x):
    return where(x < 0, 0.0, 1.0)

# Make next plot command ready
from scitools.std import *

x = linspace(-10, 10, 5)   # linspace(-10, 10, 50)
y = Hv(x)
plot(x, y, axis=[x[0], x[-1], -0.1, 1.1])

plot([-10, 0, 0, 10], [0, 0, 1, 1],
     axis=[x[0], x[-1], -0.1, 1.1])

def f(x):
    return sin(1.0/x)

x1 = linspace(-1, 1, 10)    # use 10 points
x2 = linspace(-1, 1, 1000)  # use 1000 points
plot(x1, f(x1), label='%d points' % len(x))
plot(x2, f(x2), label='%d points' % len(x))

a = x
a[-1] = 1000

a = x.copy()

a = x[r:]      # a refers to a part of the x array
a[-1] = 1000   # changes x[-1]!
a = x[r:].copy()
a[-1] = 1000   # does not change x[-1]

a = (3*x**4 + 2*x + 4)/(x + 1)

r1 = x**4;  r2 = 3*r1;  r3 = 2*x;  r4 = r1 + r3
r5 = r4 + 4;  r6 = x + 1; r7 = r5/r6; a = r7

a = x.copy()
a **= 4
a *= 3
a += 2*x
a += 4
a /= x + 1

def expression(x):
    return (3*x**4 + 2*x + 4)/(x + 1)

def expression_inplace(x):
    a = x.copy()
    a **= 4
    a *= 3
    a += 2*x
    a += 4
    a /= x + 1
    return a

import numpy as np

x = np.linspace(0, 1, 10000000)

%timeit  expression(x)

%timeit  expression_inplace(x)

from numpy import *

# x is numpy array
a = x.copy()

# or
a = zeros(x.shape, x.dtype)

# or
a = zeros_like(x)  # zeros and same size as x

type(a)

isinstance(a, ndarray)

def f(x):
    return 2

def fv(x):
    return zeros(x.shape, x.dtype) + 2

def f(x):
    if isinstance(x, (float, int)):
        return 2
    elif isinstance(x, ndarray):
        return zeros(x.shape, x.dtype) + 2
    else:
        raise TypeError(
           'x must be int/float/ndarray, not %s' % type(x))

a = linspace(1, 8, 8)
a

a[[1,6,7]] = 10
a

a[range(2,8,3)] = -2   # same as a[2:8:3] = -2
a

a < 0

a[a < 0]           # pick out all negative elements

a[a < 0] = a.max() # if a[i]<10, set a[i]=10
a

A = zeros((3,4))   # 3x4 table of numbers
A[0,0] = -1
A[1,0] =  1
A[2,0] = 10
A[0,1] = -5
A[2,3] = -100

# can also write (as for nested lists)
A[2][3] = -100

Cdegrees = [-30 + i*10 for i in range(3)]
Fdegrees = [9./5*C + 32 for C in Cdegrees]
table = [[C, F] for C, F in zip(Cdegrees, Fdegrees)]
print table

table2 = array(table)
print table2

table2.shape   # see the number of elements in each dir.

for i in range(table2.shape[0]):
    for j in range(table2.shape[1]):
        print 'table2[%d,%d] = %g' % (i, j, table2[i,j])

for index_tuple, value in np.ndenumerate(table2):
    print 'index %s has value %g' % \
          (index_tuple, table2[index_tuple])

type(index_tuple)

table2[0:, 1]             # same

table2[:, 1]              # same

t = linspace(1, 30, 30).reshape(5, 6)
t[1:-1:2, 2:]

t

t

[  7.,   8.,   9.,  10.,  11.,  12.]
[ 19.,  20.,  21.,  22.,  23.,  24.]

[  9.,  10.,  11.,  12.]
[ 21.,  22.,  23.,  24.]

from matplotlib.pyplot import *
#from scitools.std import *

# Make some data first
from numpy import *
x = linspace(0, 4*pi, 501)
y = exp(-0.1*x)*sin(x)
t1 = x
y1 = t1/t1[-1]
t2 = x[0::20]  # every 20th element in x
y2 = exp(-0.1*t2)*cos(t2)
t3 = t2
y3 = sin(0.1*t3)

# Plot

plot(x, y)         # simplest command

plot(t1, y1, 'r',  # curve 1, red line
     t2, y2, 'b',  # curve 2, blue line
     t3, y3, 'o')  # curve 3, circles at data points
axis([t1[0], t1[-1], -1.1, 1.1])
legend(['model 1', 'model 2', 'measurements'])
xlabel('time');  ylabel('force')
savefig('myframe_%04d.png' % plot_counter)

from scitools.std import *

plot(t1, y1, 'r',  # curve 1, red line
     t2, y2, 'b',  # curve 2, blue line
     t3, y3, 'o',  # curve 3, circles at data points
     axis=[t1[0], t1[-1], -1.1, 1.1],
     legend=('model 1', 'model 2', 'measurements'),
     xlabel='time', ylabel='force',
     savefig='myframe_%04d.png' % plot_counter)

from scitools.std import plot  # convenient for animations

def animate(tmax, dt, x, function, ymin, ymax, t0=0,
            xlabel='x', ylabel='y', filename='tmp_'):
    t = t0
    counter = 0
    while t <= tmax:
        y = function(x, t)
        plot(x, y,
             axis=[x[0], x[-1], ymin, ymax],
             title='time=%g' % t,
             xlabel=xlabel, ylabel=ylabel,
             savefig=filename + '%04d.png' % counter)
        t += dt
        counter += 1

# remove old plot files:
import glob, os
for filename in glob.glob('tmp_*.png'): os.remove(filename)

def T(z, t):
    # T0, A, k, and omega are global variables
    a = sqrt(omega/(2*k))
    return T0 + A*exp(-a*z)*cos(omega*t - a*z)

k = 1E-6     # heat conduction coefficient (in m*m/s)
P = 24*60*60.# oscillation period of 24 h (in seconds)
omega = 2*pi/P
dt = P/24    # time lag: 1 h
tmax = 3*P   # 3 day/night simulation
T0 = 10      # mean surface temperature in Celsius
A = 10       # amplitude of the temperature variations (in C)
a = sqrt(omega/(2*k))
D = -(1/a)*log(0.001) # max depth
n = 501      # no of points in the z direction

z = linspace(0, D, n)
animate(tmax, dt, z, T, T0-A, T0+A, 0, 'z', 'T')