numpyarrays

numpy-arrays

Arrays

Arrays in numpy

• have elements all of the same type and occupy the same amount of storage, an element can be composed of other simple types,
• have a fixed size (cannot grow like lists),
• have a shape specified with a tuple, the tuple gives the size of each dimension,
• are indexed by integers.

A numpy array is very similar to arrays in static languages like C , Fortran, and Java. Arrays in numpy are multi-dimensional. Operation are arrays are, in general, faster than lists. Why?

Python meets APL.

Creating arrays

Arrays can be created with python sequences or initialized with constant values of 0 or 1, or uninitialized. Some of the array element types are byte, int, float, complex, uint8, uint16, uint64, int8, int16, int32, int64, float32, float64, float96, complex64, complex128, and complex192.

>>> from numpy import *
>>> zeros( (12), dtype=int8 )
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int8)
>>> zeros( (2,6), dtype=int16 )
array([[0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0]], dtype=int16)
>>> ones( (6,2), dtype=int32 )
array([[1, 1],
       [1, 1],
       [1, 1],
       [1, 1],
       [1, 1],
       [1, 1]])
>>> ones( (2,2,3), dtype=int64 )
array([[[1, 1, 1],
        [1, 1, 1]],

       [[1, 1, 1],
        [1, 1, 1]]], dtype=int64)
>>> a = ones( (2,3,2) )
>>> a.dtype # type of each element
dtype('float64')
>>> a.ndim  # number of dimension
3
>>> a.shape # tuple of dimension sizes
(2, 3, 2)
>>> a.size # total number of elements
12
>>> a.itemsize # number of bytes of storage per element
8
>>> array( [ [1,2,3], [4,5,6] ] )
array([[1, 2, 3],
       [4, 5, 6]])
>>> a = _
>>> a.dtype
dtype('int32')
>>> a.shape
(2, 3)
>>> array( range(7), float )
array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.])
>>> array( [ [ int(i != j) for j in xrange(3)] for i in xrange(4)] )
array([[0, 1, 1],
       [1, 0, 1],
       [1, 1, 0],
       [1, 1, 1]])
>>> _.shape
(4, 3)
>>> empty( (2,3), dtype=int16) # uninitialized, what advantage
array([[  344, 20318, -4920],
       [ 2644,    16,     0]], dtype=int16)
>>> 

Reshaping arrays

Arrays in numpy can be reshaped. The reshaped array must be the same size. Why?

>>> from numpy import *
>>> a = array( range(12) ) # equivalent to arange( 12, dtype=int)
>>> a
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
>>> a.shape 
(12,)
>>> b = a.reshape(2, 6) # b still refers to a's memory
>>> a.shape
(12,)
>>> b.shape
(2, 6)
>>> a[6] = 106 # access 7th element
>>> b # changes in a are reflected in b
array([[  0,   1,   2,   3,   4,   5],
       [106,   7,   8,   9,  10,  11]])
>>> a
array([  0,   1,   2,   3,   4,   5, 106,   7,   8,   9,  10,  11])
>>> b[1,0] = 6 # changes in b are reflected in a
>>> b
array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11]])
>>> a
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
>>> c = b.transpose() # swaps the indices order
>>> c.shape
(6, 2)
>>> b[1,4]
10
>>> c[4,1]
10
>>> 

Indexing and slicing arrays

Arrays are indexed with comma separated list of indices. Unlike list, slices do not copy the array, but provide another view into the same data.

>>> from numpy import *
>>> t = array( range(24), uint8 ) # unsigned 8 bit integer
>>> t 
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23], dtype=uint8)
>>> t = t.reshape(2,3,4)
>>> t
array([[[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11]],

       [[12, 13, 14, 15],
        [16, 17, 18, 19],
        [20, 21, 22, 23]]], dtype=uint8)
>>> t[0,0,0] # index = i*(3*4) + j*4 + k
0
>>> t[1,1,1]
17
>>> t[1,2,3] # the last index advances serially
23
>>> # the transpose reverse the order of the indices
>>> s = t.transpose()
>>> s.shape
(4, 3, 2)
>>> s[3,2,1] # index i + j*4 + k*(3*4)
23
>>> s[0,1,0]
4
>>> s[0,0,1] # the first index advances serially
12
>>> t1 = t[1,2] # partial indices return slices
>>> t1
array([20, 21, 22, 23], dtype=uint8)
>>> t1.shape
(4,)
>>> t1[0] = 100 # in numpy a slice is similar to a reshape
>>> t           # t1 references the same array
array([[[  0,   1,   2,   3],
        [  4,   5,   6,   7],
        [  8,   9,  10,  11]],

       [[ 12,  13,  14,  15],
        [ 16,  17,  18,  19],
        [100,  21,  22,  23]]], dtype=uint8)
>>> t = array( range(24), uint8 ) ; t = t.reshape(2,3,4) # reset
>>> t[0,1]
array([4, 5, 6, 7], dtype=uint8)
>>> t[0,1] = 200 # assignment to slices changes all elements
>>> t
array([[[  0,   1,   2,   3],
        [200, 200, 200, 200],
        [  8,   9,  10,  11]],

       [[ 12,  13,  14,  15],
        [ 16,  17,  18,  19],
        [ 20,  21,  22,  23]]], dtype=uint8)
>>> 

Multi-dimensional indexing

More than one dimension can be indexed.

>>> from numpy import *
>>> z = zeros((5,5),int); z # 5 x 5 integers
array([[0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0]])
>>> z[0] = 1; z # first row
array([[1, 1, 1, 1, 1],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0]])
>>> z[:,0] = 2; z # first column
array([[2, 1, 1, 1, 1],
       [2, 0, 0, 0, 0],
       [2, 0, 0, 0, 0],
       [2, 0, 0, 0, 0],
       [2, 0, 0, 0, 0]])
>>> z[2:4,2:5] = 3; z # (2,3) block
array([[2, 1, 1, 1, 1],
       [2, 0, 0, 0, 0],
       [2, 0, 3, 3, 3],
       [2, 0, 3, 3, 3],
       [2, 0, 0, 0, 0]])
>>> z[:,-1] = 4; z # last column
array([[2, 1, 1, 1, 4],
       [2, 0, 0, 0, 4],
       [2, 0, 3, 3, 4],
       [2, 0, 3, 3, 4],
       [2, 0, 0, 0, 4]])
>>> z[-1,:] = 5; z # last row
array([[2, 1, 1, 1, 4],
       [2, 0, 0, 0, 4],
       [2, 0, 3, 3, 4],
       [2, 0, 3, 3, 4],
       [5, 5, 5, 5, 5]])
>>> 

arithmetic, vector and matrix operation

The standard math operators can be used with numpy arrays. Matrix operations are also defined.

>>> from numpy import *
>>> i = identity( 3, int16 )
>>> i
array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1]], dtype=int16)
>>> i + i # add element to element
array([[2, 0, 0],
       [0, 2, 0],
       [0, 0, 2]], dtype=int16)
>>> i + 4 # add a scalar to every entry
array([[5, 4, 4],
       [4, 5, 4],
       [4, 4, 5]], dtype=int16)
>>> a = array( range(1,10), int16 ) # equiv to arange(1,10,dtype=int16):w
>>> a = a.reshape(3,3)
>>> i * a # element to element
array([[1, 0, 0],
       [0, 5, 0],
       [0, 0, 9]], dtype=int16)
>>> x = array( [1,2,3], int32 )
>>> x
array([1, 2, 3])
>>> y = array( [ [4], [5], [6] ], int32 )
>>> y
array([[4],
       [5],
       [6]])
>>> x + y
array([[5, 6, 7],
       [6, 7, 8],
       [7, 8, 9]])
>>> array([4]) + array([1,2,3]) # one element with three
array([5, 6, 7])
>>> # 1 and 2d arrays can be turned into matrices with mat()
>>> # matrices redefine * to be dot
>>> mx = mat(x)
>>> mx.shape
(1, 3)
>>> my = mat(y)
>>> my.shape
(3, 1)
>>> mx * my
matrix([[32]])
>>> my * mx
matrix([[ 4,  8, 12],
        [ 5, 10, 15],
        [ 6, 12, 18]])
>>> 

Relational operators

Relational operators are also defined for numpy arrays. An array of booleans is returned.

>>> from numpy import *
>>> a = array([3,6,8,9])
>>> a == 6
array([False,  True, False, False], dtype=bool)
>>> a >= 7
array([False, False,  True,  True], dtype=bool)
>>> a < 5
array([ True, False, False, False], dtype=bool)
>>> # count all the even numbers
>>> sum( (a%2) == 0 )
2
>>> b = array([2,6,7,10])
>>> a == b
array([False,  True, False, False], dtype=bool)
>>> a < b
array([False, False, False,  True], dtype=bool)
>>> # any is true if any value is true
>>> any( a == 3 )
True
>>> any( a == 10 )
False
>>> # all is true if all values are true
>>> a == a
array([ True,  True,  True,  True], dtype=bool)
>>> all( a == a )
True
>>> all( a != b )
False
>>> 

Indexing with boolean arrays

In addition to slices, a numpy array can be indexed with an equivalently shaped boolean array (i.e., their shape are the same).

>>> from numpy import *
>>> a = arange(1, 13).reshape(3,4); a
array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12]])
>>> i1 = (a % 2) == 1; i1 # odd entries
array([[ True, False,  True, False],
       [ True, False,  True, False],
       [ True, False,  True, False]], dtype=bool)
>>> a[ i1 ] # select all the odd entries
array([ 1,  3,  5,  7,  9, 11])
>>> i2 = (a >= 1) & (a < 4) ; i2 # and
array([[ True,  True,  True, False],
       [False, False, False, False],
       [False, False, False, False]], dtype=bool)
>>> a[ i2 ] # values greater than or equal to 1 and less than 4
array([1, 2, 3])
>>> 

Math functions

Most of the math functions can be used with numpy arrays. These functions are actually redefined. Some of the redefined math functions are: arccos, arcsin, arctan, arctan2, ceil, cos, cosh, degrees, exp, fabs, floor, fmod, frexp, hypot, ldexp, log, log10, modf, pow, radians, sin, sinh, sqrt, tan, and tanh.

>>> from numpy import *
>>> t = linspace(0, pi, 7) # 7 values linearly between 0 and pi
>>> sin(t)
array([  0.00000000e+00,   5.00000000e-01,   8.66025404e-01,
         1.00000000e+00,   8.66025404e-01,   5.00000000e-01,
         1.22460635e-16])
>>> arcsin( sin( t ) )
array([  0.00000000e+00,   5.23598776e-01,   1.04719755e+00,
         1.57079633e+00,   1.04719755e+00,   5.23598776e-01,
         1.22460635e-16])
>>> floor( t )
array([ 0.,  0.,  1.,  1.,  2.,  2.,  3.])
>>> log( t + 1 ) # adding 1 avoids log(0)
array([ 0.        ,  0.42107515,  0.71647181,  0.94421571,  1.12959244,
        1.28591969,  1.42108041])
>>> 

array to scalar operators

numpy also defines operators that reduce the number of elements. These operators are: max, mean, min, std, sum, var.

>>> from numpy import *
>>> s = random.standard_normal( 10 ) # mu = 0.0 sigma = 1.0
>>> s.shape
(10,)
>>> s
array([ 0.65711932, -1.19741536,  1.51470124,  0.60134355,  1.44631555,
        1.25936877, -1.347354  ,  0.33819449,  0.35765847,  0.84350667])
>>> s.max() # all elements
1.51470124183
>>> s.min()
-1.34735399976
>>> s.sum()  # sum of all elements
4.47343871645
>>> s.prod() # product of all elements
0.179456818521
>>> s.mean()
0.447343871645
>>> s.std() # square root of variance
0.946962388091
>>> s.var() # variance 
0.896737764459
>>> 

The operator can also select an axis (direction) for grouping the elements to be operated on.

>>> from numpy import *
>>> a = array( [[4, -9, 5, 7], [100, -100, 50, 21], [20, 40, 0, 11]] )
>>> a
array([[   4,   -9,    5,    7],
       [ 100, -100,   50,   21],
       [  20,   40,    0,   11]])
>>> a.max()
100
>>> a.min()
-100
>>> a.sum()
149
>>> a.min(0) # along the columns, varies last index
array([   4, -100,    0,    7])
>>> a.max(0) # along the columns 
array([100,  40,  50,  21])
>>> a.sum(0) # along the columns
array([124, -69,  55,  39])
>>> a.min(1) # along the rows, varies first index
array([  -9, -100,    0])
>>> a.max(1) # along the rows
array([  7, 100,  40])
>>> a.sum(1) # along the rows
array([ 7, 71, 71])
>>> 

polynomials

Polynomials are supported by numpy. Some examples are:

>>> from numpy import *
>>> p1 = poly1d([2.0,1.0]) # 2*x + 1, using coefficients
>>> print p1
 
2 x + 1
>>> # using the polynomials roots
>>> p1a = poly1d([-0.5],r=1) # p1(-0.5) = 0
>>> print p1a
 
1 x + 0.5
>>> # quadratic with roots (x-5)(x-3)
>>> p2 = poly1d( [5.0, 3.0], r=1 )
>>> print p2
   2
1 x - 8 x + 15
>>> # polynomials can be multiplied
>>> p3 = poly1d([1.0,-5]) * poly1d([1.0,3.0])
>>> print p3
   2
1 x - 2 x - 15
>>> # or added
>>> poly1d([1.0,-5]) + poly1d([1.0,3.0])
poly1d([ 2., -2.])
>>> # or subtracted
>>> poly1d([3.0,1.0,-5]) - poly1d([1.0,3.0])
poly1d([ 3.,  0., -8.])
>>> # or division with remainder
>>> quot, rem =  poly1d([3.0,1.0,-5]) /  poly1d([1.0,3.0])
>>> print quot
 
3 x - 8
>>> print rem
 
19
>>> # polynomials can be evaluated
>>> p3(-3)
0.0
>>> p3(5)
0.0
>>> p3(1)
-16.0
>>> # the polynomials roots are given by
>>> p3.r
array([ 5., -3.])
>>> # its coefficients
>>> p3.c
array([  1.,  -2., -15.])
>>> # its order
>>> p3.o
2
>>> # the first derivative
>>> print p3.deriv()
 
2 x - 2
>>> print p3.deriv(1)
 
2 x - 2
>>> # the second derivative
>>> print p3.deriv(2)
 
2
>>> 

Course record

The students records for a course is stored in the following csv file.

course.csv

Name, Assignments, Midterm1, Midterm2, Final
A Student, 80, 65, 72, 74 
B Student, 75, 52, 62, 64
C Student, 100, 70, 80, 80
D Student, 94, 45, 65, 55

The weights for the term work, midterm 1, midterm 2, and the final is 10, 20, 20, and 50.

Course report

A report giving statistics on each work product, and the student's course grade is generated by:

course_report.py

#!/usr/bin/env python
from numpy import *
import csv

# get head and names
f = file( 'course.csv' )
reader = csv.reader( f )
headings = reader.next()
# get student names
students = []
for row in reader :
    students.append( row[0] )
f.close()

# get workproducts, skipping name column
f = file( 'course.csv' )
wps = loadtxt( f, delimiter=",",
                usecols=(1,2,3,4), skiprows=1 )
f.close()

# value of each workproduct item
value = array( [10., 20., 20., 50.] )
# report final mark
print 'Course mark'
marks = [] # save the marks
for i, name in enumerate( students ) :
    mark = dot( value, wps[i]/100.0 )
    marks.append( mark )
    print "%-20s: %5.1f" % (name, mark)

marks = array( marks ) # convert to array

print "\n"*4
print 'Statistics'
# report on statistics of each workproduct
wp_names = headings[1:]
for i, w in enumerate( wp_names ) :
    print w.strip()
    print '    minimum %6.2f' %  wps[:,i].min()
    print '    maximum %6.2f' %  wps[:,i].max()
    print '       mean %6.2f' %  wps[:,i].mean()
    print '        std %6.2f' %  wps[:,i].std()

# generated histogram of As, Bs, Cs, Ds, Fs
bins = array([0.0, 49.5, 54.5, 64.5, 79.5, 100.0] ) # F,D,C,B,A
N,bins = histogram(marks,bins)

N = list(N) # convert to list
N.reverse()
N.pop( 0 ) # discard overflow, XXX could check
for letter,n in zip( "ABCDF", N) :
    print letter, n

The produced report is:

Command:

python course_report.py

Output:

Course mark
A Student           :  72.4
B Student           :  62.3
C Student           :  80.0
D Student           :  58.9





Statistics
Assignments
    minimum  75.00
    maximum 100.00
       mean  87.25
        std  10.13
Midterm1
    minimum  45.00
    maximum  70.00
       mean  58.00
        std   9.97
Midterm2
    minimum  62.00
    maximum  80.00
       mean  69.75
        std   6.94
Final
    minimum  55.00
    maximum  80.00
       mean  68.25
        std   9.55
A 1
B 1
C 2
D 0
F 0

Conway's game of life

The games consists of an infinite board composed of a grid of squares. Each square can contain a zero or a one. Each square has eight neighbours, the next square connected either horizontally, vertically, or on a diagonal. Give a board configuration, the next configuration is determined by:

• A square with a 1 and fewer than two neighbours with a 1 is set to 0.
• A square with a 1 and more than three neighbours with a 1 is set to 0.
• A square with a 1 and two or three neighbours with a 1, remains set to 0.
• A square with a 0, and exactly three neighbours with a 1 is set to a 0.

Conway's game of life in numpy

Conway's game can be easily implemented with numpy's array operations. A bounded implementation (i.e., the board is not infinite) where the borders are always zero is realized by:

simple_life.py

#!/usr/bin/env python
from numpy import *

# compute next generation of conway's game of life
def next_generation( current, next ) :
    next[:,:] = 0 # zero out next board
    # bound examination area
    bend0 = (current.shape[0]-3)+1
    bend1 = (current.shape[1]-3)+1
    for i in xrange( bend0 ) :
        for j in xrange( bend1 ) :
            neighbours = sum( current[i:i+3, j:j+3] )
            if current[i+1, j+1] == 1 :
                neighbours -= 1 # do not count yourself
                if 2 <= neighbours <= 3 :
                    next[i+1, j+1] = 1
            else:
                if neighbours == 3 :
                    next[i+1,j+1] = 1

board = zeros( (5,5), uint8)
next_board = zeros_like( board ) # same shape

# create inital configuration
board[2,1:4] = 1 # 3 squares 

for i in xrange(4) :
    next_generation( board, next_board )
    board, next_board = next_board, board # swap boards
    print board[1:-1,1:-1]

Four generation of this configuration yeilds:

Command:

python simple_life.py

Output:

[[0 1 0]
 [0 1 0]
 [0 1 0]]
[[0 0 0]
 [1 1 1]
 [0 0 0]]
[[0 1 0]
 [0 1 0]
 [0 1 0]]
[[0 0 0]
 [1 1 1]
 [0 0 0]]

Game of life with user input

An input string of ones and zeros can be used to initialize the board. This version is identical to the previous version with the exception of the command line arguments.

life.py

#!/usr/bin/env python
from numpy import *
import sys
import math

# compute next generation of conway's game of life
def next_generation( current, next ) :
    next[:,:] = 0 # zero out next board
    # bound examination area
    bend0 = (current.shape[0]-3)+1
    bend1 = (current.shape[1]-3)+1
    for i in xrange( bend0 ) :
        for j in xrange( bend1 ) :
            neighbours = sum( current[i:i+3, j:j+3] )
            if current[i+1, j+1] == 1 :
                neighbours -= 1 # do not count yourself
                if 2 <= neighbours <= 3 :
                    next[i+1, j+1] = 1
            else:
                if neighbours == 3 :
                    next[i+1,j+1] = 1

if len(sys.argv) != 3 :
    print "usage:", sys.argv[0], "init-board generations"
    sys.exit( 1 )

init = sys.argv[1]
generations = int( sys.argv[2] )

board_sz = math.ceil( math.sqrt(len(init)) ) 

# board_sz+2 includes the border of zeros
board = zeros( (board_sz+2,board_sz+2), uint8)
next_board = zeros_like( board ) # same shape

# fill the board
i = 0
j = 0
for index,ch in enumerate(init) :
    if ch == '1' :
        board[i+1,j+1] = 1
    j = (j+1) % board_sz
    if ((index+1) % board_sz) == 0 : i += 1

for gen in xrange(generations) :
    next_generation( board, next_board )
    board, next_board = next_board, board # swap boards
    print board[1:-1,1:-1] # do not print the border

Solving 2 unknowns with 2 linear equations

>>> from numpy import *
>>> from numpy.linalg import solve
>>> 
>>> # 1 * x0 - 2 * x1 = 1
>>> # 2 * x0 - 6 * x1 = -2
>>> a = array( [[1,-2], [2,-6]] )
>>> b = array( [1, -2 ] )
>>> x = solve( a, b ); x
array([ 5.,  2.])
>>> dot( a, x ) # check solution
array([ 1., -2.])
>>> 

Solving 2 unknowns with 2 linear equations (2)

>>> from numpy import *
>>> from numpy.linalg import solve
>>> 
>>> # 1 * x0 - 2 * x1 = 1
>>> # 2 * x0 - 4 * x1 = 2
>>> # any problems with the above equations
>>> a = array( [[1,-2], [2,-4]] )
>>> b = array( [1, 2 ] )
>>> x = solve( a, b )
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "/usr/lib/python2.5/site-packages/numpy/linalg/linalg.py", line 184, in solve
    raise LinAlgError, 'Singular matrix'
LinAlgError: Singular matrix
>>> 

Solving 3 unknowns with 3 linear equations

>>> from numpy import *
>>> from numpy.linalg import solve
>>> from numpy.linalg import inv
>>> 
>>> # 1 * x0 - 2 * x1 + 5 * x2 = 1
>>> # 2 * x0 - 6 * x1 - 10 *x2 = -2
>>> #-4 * x0 + 1 * x1 + 3 * x2 = 3
>>> a = array( [[1, -2, 5], [2,-6,-10], [-4,1,3] ] )
>>> b = array( [1, -2, 3 ] )
>>> x = solve( a, b ); x
array([-0.64516129, -0.25806452,  0.22580645])
>>> dot( a, x ) # check solution
array([ 1., -2.,  3.])
>>> 
>>> # a * x = b
>>> # inv(a) * a * x = inv(a) * b
>>> # I * x = inv(a) * b
>>> # x = inv(a) * b
>>> x = dot( inv(a), b ); x
array([-0.64516129, -0.25806452,  0.22580645])
>>> # the identity
>>> dot( inv(a), a)
array([[  1.00000000e+00,   5.55111512e-17,   2.77555756e-17],
       [ -2.77555756e-17,   1.00000000e+00,   1.24900090e-16],
       [  6.93889390e-18,  -1.90819582e-17,   1.00000000e+00]])
>>> dot( a, inv(a))
array([[  1.00000000e+00,   5.55111512e-17,   1.56125113e-17],
       [ -2.77555756e-17,   1.00000000e+00,   5.20417043e-17],
       [  1.38777878e-17,   1.38777878e-17,   1.00000000e+00]])
>>> 

2D line segment intersection

Many geometric problems are solved using vector notation. Intesection of lines can be solved with:

lines.py

#
# line segment intersection using vectors
# see Computer Graphics by F.S. Hill
#
from numpy import *
def perp( a ) :
    b = empty_like(a)
    b[0] = -a[1]
    b[1] = a[0]
    return b

# line segment a given by endpoints a1, a2
# line segment b given by endpoints b1, b2
# return 
def seg_intersect(a1,a2, b1,b2) :
    da = a2-a1
    db = b2-b1
    dp = a1-b1
    dap = perp(da)
    denom = dot( dap, db)
    num = dot( dap, dp )
    return (num / denom)*db + b1

p1 = array( [0.0, 0.0] )
p2 = array( [1.0, 0.0] )

p3 = array( [4.0, -5.0] )
p4 = array( [4.0, 2.0] )

print seg_intersect( p1,p2, p3,p4)

p1 = array( [2.0, 2.0] )
p2 = array( [4.0, 3.0] )

p3 = array( [6.0, 0.0] )
p4 = array( [6.0, 3.0] )

print seg_intersect( p1,p2, p3,p4)

Its sample execution is:

Command:

python lines.py

Output:

[ 4.  0.]
[ 6.  4.]

Points on circle "without" sin/cos

All the points around a circle can be generated with matrix operations. Multiplication is usually cheaper than calculationg sin or cos.

circle.py

from numpy import *
import math

# sin/cos just once
a = math.cos( math.radians(-30.0) )
b = math.sin( math.radians(-30.0) )

# rotation matrix, 30 degrees clockwise
R = mat( [[a, -b], [b, a]] )
v = mat( [1,0] )   # unit vector along x-axis

print "Repeated rotations by 30"
# rotate (1,0) by 30 degree increments
for i in xrange(12) :
    nv = v * R
    vx,vy = v.T  # transpose
    x = math.cos( math.radians(30.0 * i ) )
    y = math.sin( math.radians(30.0 * i ) )
    print "%12.9f %12.9f -- %12.9f %12.9f" % (vx,vy, x,y)
    v = nv
print

print "Rotation by multiples of 30"
# rotate (1,0) by 30, 60, 90, ...
for i in xrange(12) :
    r = v * pow(R,i)
    vx,vy = r.T  # transpose
    x = math.cos( math.radians(30.0 * i ) )
    y = math.sin( math.radians(30.0 * i ) )
    print "%12.9f %12.9f -- %12.9f %12.9f" % (vx,vy, x,y)

The points are:

Command:

python circle.py

Output:

Repeated rotations by 30
 1.000000000  0.000000000 --  1.000000000  0.000000000
 0.866025404  0.500000000 --  0.866025404  0.500000000
 0.500000000  0.866025404 --  0.500000000  0.866025404
 0.000000000  1.000000000 --  0.000000000  1.000000000
-0.500000000  0.866025404 -- -0.500000000  0.866025404
-0.866025404  0.500000000 -- -0.866025404  0.500000000
-1.000000000  0.000000000 -- -1.000000000  0.000000000
-0.866025404 -0.500000000 -- -0.866025404 -0.500000000
-0.500000000 -0.866025404 -- -0.500000000 -0.866025404
-0.000000000 -1.000000000 -- -0.000000000 -1.000000000
 0.500000000 -0.866025404 --  0.500000000 -0.866025404
 0.866025404 -0.500000000 --  0.866025404 -0.500000000

Rotation by multiples of 30
 1.000000000 -0.000000000 --  1.000000000  0.000000000
 0.866025404  0.500000000 --  0.866025404  0.500000000
 0.500000000  0.866025404 --  0.500000000  0.866025404
 0.000000000  1.000000000 --  0.000000000  1.000000000
-0.500000000  0.866025404 -- -0.500000000  0.866025404
-0.866025404  0.500000000 -- -0.866025404  0.500000000
-1.000000000  0.000000000 -- -1.000000000  0.000000000
-0.866025404 -0.500000000 -- -0.866025404 -0.500000000
-0.500000000 -0.866025404 -- -0.500000000 -0.866025404
-0.000000000 -1.000000000 -- -0.000000000 -1.000000000
 0.500000000 -0.866025404 --  0.500000000 -0.866025404
 0.866025404 -0.500000000 --  0.866025404 -0.500000000