The Basics¶
An example¶
>>> import numpy as np
>>> a = np.arange(15).reshape(3, 5)
>>> a
array([[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9],
[10, 11, 12, 13, 14]])
>>> a.shape
(3, 5)
>>> a.ndim
2
>>> a.dtype.name
'int64'
>>> a.itemsize
8
>>> a.size
15
>>> type(a)
<type 'numpy.ndarray'>
>>> b = np.array([6, 7, 8])
>>> b
array([6, 7, 8])
>>> type(b)
<type 'numpy.ndarray'>
>>> import Numeric.LinearAlgebra
>>> import Foreign.Storable (sizeOf)
>>> let a = (3><5) [0::I ..] -- or ((3><5) . toList . range) 15
>>> a -- for vectors that are just indices, using @range 15@ is fine
(3><5)
[ 0, 1, 2, 3, 4
, 5, 6, 7, 8, 9
, 10, 11, 12, 13, 14 ]
>>> size a
(3,5)
-- there is no @ndim@ for hmatrix, since there are only scalar, vectors, and (2D) matrices
>>> :t (a ! 0 ! 0)
(a ! 0 ! 0) :: I
>>> sizeOf (a ! 0 ! 0)
4
>>> size (flatten a)
15
>>> :t a
a :: Matrix I
>>> let b = fromList [6,7,8]
>>> b
[6.0,7.0,8.0]
>>> :t b
b :: (Foreign.Storable.Storable a, Num a) => Vector a
Array creation¶
>>> import numpy as np
>>> a = np.array([2,3,4])
>>> a
array([2, 3, 4])
>>> a.dtype
dtype('int64')
>>> b = np.array([1.2, 3.5, 5.1])
>>> b.dtype
dtype('float64')
>>> import Numeric.LinearAlgebra
>>> let a = fromList [2,3,4::Z]
>>> a
[2,3,4]
>>> :t a
a :: Vector Z
>>> let b = fromList [1.2, 3.5, 5.1::R]
>>> :t b
b :: Vector R
-- If you do not cast the numerical literal, you don't get a fully-resolved-type vector
>>> b = np.array([(1.5,2,3), (4,5,6)])
>>> b
array([[ 1.5, 2. , 3. ],
[ 4. , 5. , 6. ]])
>>> let b = fromLists [[1.5,2,3], [4,5,6::R]]
>>> b
(2><3)
[ 1.5, 2.0, 3.0
, 4.0, 5.0, 6.0 ]
>>> c = np.array( [ [1,2], [3,4] ], dtype=complex )
>>> c
array([[ 1.+0.j, 2.+0.j],
[ 3.+0.j, 4.+0.j]])
>>> let c = (fromLists [ [1,2], [3,4] ]) :: Matrix C
>>> c
(2><2)
[ 1.0 :+ 0.0, 2.0 :+ 0.0
, 3.0 :+ 0.0, 4.0 :+ 0.0 ]
>>> np.zeros( (3,4) )
array([[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.]])
>>> np.ones( (2,3,4), dtype=np.int16 ) # dtype can also be specified
array([[[ 1, 1, 1, 1],
[ 1, 1, 1, 1],
[ 1, 1, 1, 1]],
[[ 1, 1, 1, 1],
[ 1, 1, 1, 1],
[ 1, 1, 1, 1]]], dtype=int16)
>>> np.empty( (2,3) ) # uninitialized, output may vary
array([[ 3.73603959e-262, 6.02658058e-154, 6.55490914e-260],
[ 5.30498948e-313, 3.14673309e-307, 1.00000000e+000]])
>>> konst (0::R) (3::Int,4::Int)
(3><4)
[ 0.0, 0.0, 0.0, 0.0
, 0.0, 0.0, 0.0, 0.0
, 0.0, 0.0, 0.0, 0.0 ]
-- or
>>> ((3><4) . repeat) 0
(3><4)
[ 0.0, 0.0, 0.0, 0.0
, 0.0, 0.0, 0.0, 0.0
, 0.0, 0.0, 0.0, 0.0 ]
-- there are no 3D tensors in hmatrix
-- there are no "empty" matrices in hmatrix
>>> np.arange( 10, 30, 5 )
array([10, 15, 20, 25])
>>> np.arange( 0, 2, 0.3 ) # it accepts float arguments
array([ 0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8])
>>> vector (init [10, 15.. 30])
[10.0,15.0,20.0,25.0]
>>> vector (init [0, 0.3.. 2])
[0.0,0.3,0.6,0.8999999999999999,1.1999999999999997,1.4999999999999996,1.7999999999999994]
>>> from numpy import pi
>>> np.linspace( 0, 2, 9 ) # 9 numbers from 0 to 2
array([ 0. , 0.25, 0.5 , 0.75, 1. , 1.25, 1.5 , 1.75, 2. ])
>>> x = np.linspace( 0, 2*pi, 100 ) # useful to evaluate function at lots of points
>>> f = np.sin(x)
>>> linspace 9 (0, 2::Double)
[0.0,0.25,0.5,0.75,1.0,1.25,1.5,1.75,2.0]
>>> let x = linspace 100 (0, 2*pi) :: Vector Double
>>> let f = cmap sin x
Printing arrays¶
>>> a = np.arange(6) # 1d array
>>> print(a)
[0 1 2 3 4 5]
>>>
>>> b = np.arange(12).reshape(4,3) # 2d array
>>> print(b)
[[ 0 1 2]
[ 3 4 5]
[ 6 7 8]
[ 9 10 11]]
>>>
>>> c = np.arange(24).reshape(2,3,4) # 3d array
>>> print(c)
[[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
[[12 13 14 15]
[16 17 18 19]
[20 21 22 23]]]
>>> let a = 6 |> [0..]
>>> disp 2 (asRow a)
1x6
0 1 2 3 4 5
>>> let b = (4><3) [0..]
>>> disp 2 b
4x3
0 1 2
3 4 5
6 7 8
9 10 11
-- there are no 3D matrices in hmatrix
>>> print(np.arange(10000))
[ 0 1 2 ..., 9997 9998 9999]
>>>
>>> print(np.arange(10000).reshape(100,100))
[[ 0 1 2 ..., 97 98 99]
[ 100 101 102 ..., 197 198 199]
[ 200 201 202 ..., 297 298 299]
...,
[9700 9701 9702 ..., 9797 9798 9799]
[9800 9801 9802 ..., 9897 9898 9899]
[9900 9901 9902 ..., 9997 9998 9999]]
>>> dispShort 6 6 0 (asRow (10000 |> [0..]))
1x10000
0 1 2 .. 9998 9999
>>> dispShort 6 6 0 ((100><100) [0..])
100x100
0 1 2 .. 98 99
100 101 102 .. 198 199
200 201 202 .. 298 299
::
9800 9801 9802 .. 9898 9899
9900 9901 9902 .. 9998 9999
Basic Operations¶
>>> a = np.array( [20,30,40,50] )
>>> b = np.arange( 4 )
>>> b
array([0, 1, 2, 3])
>>> c = a-b
>>> c
array([20, 29, 38, 47])
>>> b**2
array([0, 1, 4, 9])
>>> 10*np.sin(a)
array([ 9.12945251, -9.88031624, 7.4511316 , -2.62374854])
>>> a<35
array([ True, True, False, False], dtype=bool)
>>> let a = fromList [20,30,40,50::R]
>>> let b = 4 |> [0::R ..]
>>> b
[0.0,1.0,2.0,3.0]
>>> let c = a-b
>>> c
[20.0,29.0,38.0,47.0]
>>> cmap (^ 2) b
[0.0,1.0,4.0,9.0]
>>> disp 8 (asRow (cmap ((*10) . sin) a))
1x4
9.12945251 -9.88031624 7.45113160 -2.62374854
>>> (fromList . (fmap (<35)) . toList) a -- hmatrix does not have an Element Bool definition
[True,True,False,False]
>>> A = np.array( [[1,1],
... [0,1]] )
>>> B = np.array( [[2,0],
... [3,4]] )
>>> A*B # elementwise product
array([[2, 0],
[0, 4]])
>>> A.dot(B) # matrix product
array([[5, 4],
[3, 4]])
>>> np.dot(A, B) # another matrix product
array([[5, 4],
[3, 4]])
>>> let a = fromLists [[1,1],[0,1]] :: Matrix R
>>> let b = fromLists [[2,0],[3,4]] :: Matrix R
>>> a * b -- capitals are reserved for data constructors
(2><2)
[ 2.0, 0.0
, 0.0, 4.0 ]
>>> a <> b
(2><2)
[ 5.0, 4.0
, 3.0, 4.0 ]
>>> import Numeric.LinearAlgebra.HMatrix (mul)
>>> mul a b -- dot is reserved for vector <.> vector
(2><2)
[ 5.0, 4.0
, 3.0, 4.0 ]
>>> a = np.ones((2,3), dtype=int)
>>> b = np.random.random((2,3))
>>> a *= 3
>>> a
array([[3, 3, 3],
[3, 3, 3]])
>>> b += a
>>> b
array([[ 3.417022 , 3.72032449, 3.00011437],
[ 3.30233257, 3.14675589, 3.09233859]])
>>> a += b # b is not automatically converted to integer type
Traceback (most recent call last):
...
TypeError: Cannot cast ufunc add output from dtype('float64') to dtype('int64') with casting rule 'same_kind'
>>> import Numeric.LinearAlgebra.HMatrix (Seed, RandDist(..), randomVector, rand)
>>> let a' = ((2><3) . repeat) 1 -- inplace mutation not useful here
>>> let rand' seed r c = reshape c (randomVector seed Uniform (r*c))
>>> let b' = rand' 1 2 3 -- or unseeded: rand 2 3; rand' not actually deterministic (?)
>>> let a = cmap (*3) a' :: Matrix Z -- haskell cannot reassign variables using itself
>>> a
(2><3)
[ 3, 3, 3
, 3, 3, 3 ]
>>> let b = (+ (fromZ a)) b' -- need to cast away from integer (Z)
>>> disp 7 b -- unseeded would need fmap due to being IO (...)
2x3
3.4583777 3.1595528 3.4802844
3.6986673 3.8688027 3.5846251
>>> a = np.ones(3, dtype=np.int32)
>>> b = np.linspace(0,pi,3)
>>> b.dtype.name
'float64'
>>> c = a+b
>>> c
array([ 1. , 2.57079633, 4.14159265])
>>> c.dtype.name
'float64'
>>> d = np.exp(c*1j)
>>> d
array([ 0.54030231+0.84147098j, -0.84147098+0.54030231j,
-0.54030231-0.84147098j])
>>> d.dtype.name
'complex128'
>>> :set -XFlexibleContexts
>>> let a = konst 1 3 -- type not resolved here
>>> let b = linspace 3 (0, pi) -- set FlexibleContexts to get type inference
>>> :t b
b :: (Container Vector e, Floating e) => Vector e
>>> let c = a + b :: Vector R -- we can force resolving type of only c
>>> disp 8 (asRow (c :: Vector R)) -- a and b will remain unresolved
1x3
1.00000000 2.57079633 4.14159265
>>> :t c
c :: Vector R
>>> let d = cmap (exp . (*iC)) (complex c)
>>> putStrLn (dispcf 8 (asRow d))
1x3
0.54030231+0.84147098i -0.84147098+0.54030231i -0.54030231-0.84147098i
>>> :t d
d :: Vector C
>>> a = np.random.random((2,3))
>>> a
array([[ 0.18626021, 0.34556073, 0.39676747],
[ 0.53881673, 0.41919451, 0.6852195 ]])
>>> a.sum()
2.5718191614547998
>>> a.min()
0.1862602113776709
>>> a.max()
0.6852195003967595
>>> let a = rand' 1 2 3
>>> a
(2><3)
[ 0.18626021, 0.34556073, 0.39676747
, 0.53881673, 0.41919451, 0.6852195 ]
>>> sumElements a :: R
2.57181915
>>> minElement a :: R
0.18626021
>>> maxElement a :: R
0.6852195
>>> b = np.arange(12).reshape(3,4)
>>> b
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>>
>>> b.sum(axis=0) # sum of each column
array([12, 15, 18, 21])
>>>
>>> b.min(axis=1) # min of each row
array([0, 4, 8])
>>>
>>> b.cumsum(axis=1) # cumulative sum along each row
array([[ 0, 1, 3, 6],
[ 4, 9, 15, 22],
[ 8, 17, 27, 38]])
>>> let b = (3><4) [0::I ..]
>>> b
(3><4)
[ 0, 1, 2, 3
, 4, 5, 6, 7
, 8, 9, 10, 11 ]
>>> (fromList . fmap sumElements . toColumns) b
[12,15,18,21]
>>> (fromList . fmap minElement . toRows) b
[0,4,8]
>>> (fromLists . (fmap (scanl1 (+))) . toLists) b
(3><4)
[ 0, 1, 3, 6
, 4, 9, 15, 22
, 8, 17, 27, 38 ]
Universal Functions¶
>>> B = np.arange(3)
>>> B
array([0, 1, 2])
>>> np.exp(B)
array([ 1. , 2.71828183, 7.3890561 ])
>>> np.sqrt(B)
array([ 0. , 1. , 1.41421356])
>>> C = np.array([2., -1., 4.])
>>> np.add(B, C)
array([ 2., 0., 6.])
>>> let b = 3 |> [0::R ..] -- if we used @I@, we'd need @fromInt b@ below
>>> b
[0.0,1.0,2.0]
>>> cmap exp b
[1.0,2.718281828459045,7.38905609893065]
>>> cmap sqrt b
[0.0,1.0,1.4142135623730951]
>>> let c = vector [2, -1, 4]
>>> b + c
[2.0,0.0,6.0]
Indexing, Slicing, and Iterating¶
>>> a = np.arange(10)**3
>>> a
array([ 0, 1, 8, 27, 64, 125, 216, 343, 512, 729])
>>> a[2]
8
>>> a[2:5]
array([ 8, 27, 64])
>>> a[:6:2] = -1000 # equivalent to a[0:6:2] = -1000; from start to position 6, exclusive, set every 2nd element to -1000
>>> a
array([-1000, 1, -1000, 27, -1000, 125, 216, 343, 512, 729])
>>> a[ : :-1] # reversed a
array([ 729, 512, 343, 216, 125, -1000, 27, -1000, 1, -1000])
>>> for i in a:
... print(i**(1/3.))
...
nan
1.0
nan
3.0
nan
5.0
6.0
7.0
8.0
9.0
>>> let a = build 10 (^3) :: Vector I -- helper function with implicit counter
>>> a
[0,1,8,27,64,125,216,343,512,729]
>>> a ! 2 -- also can do @a `atIndex` 2@
8
>>> subVector 2 (5-2) a
[8,27,64]
-- that uses @Data.Vector.slice@ underneath, without copying
-- we could convert to a matrix for nice indexing
-- using @(flatten . ((flip (??)) (All, Range 2 1 4)) . asRow) a@
-- we could also do @(flatten . remap (asColumn (scalar 0)) (asRow (fromList [2..4])) . asRow) b@
-- or even @fromList ((sequence . fmap (flip (!))) [2..4] b)@
>>> let a' = accum a const (zip [0,2..4] (repeat (negate 1000)))
>>> a'
[-1000,1,-1000,27,-1000,125,216,343,512,729]
>>> import qualified Data.Vector.Storable as V
>>> V.reverse a' -- unlike in Python, this is not a view, and is O(n)
[729,512,343,216,125,-1000,27,-1000,1,-1000]
>>> import Numeric.LinearAlgebra.Devel (mapVectorM_)
>>> mapVectorM_ print (cmap (**(1/3.0)) (fromInt a') :: Vector R)
NaN
1.0
NaN
3.0
NaN
4.999999999999999
6.0
6.999999999999998
8.0
8.999999999999998
>>> def f(x,y):
... return 10*x+y
...
>>> b = np.fromfunction(f,(5,4),dtype=int)
>>> b
array([[ 0, 1, 2, 3],
[10, 11, 12, 13],
[20, 21, 22, 23],
[30, 31, 32, 33],
[40, 41, 42, 43]])
>>> b[2,3]
23
>>> b[0:5, 1] # each row in the second column of b
array([ 1, 11, 21, 31, 41])
>>> b[ : ,1] # equivalent to the previous example
array([ 1, 11, 21, 31, 41])
>>> b[1:3, : ] # each column in the second and third row of b
array([[10, 11, 12, 13],
[20, 21, 22, 23]])
>>> b[-1] # the last row. Equivalent to b[-1,:]
array([40, 41, 42, 43])
>>> c = np.array( [[[ 0, 1, 2], # a 3D array (two stacked 2D arrays)
... [ 10, 12, 13]],
... [[100,101,102],
... [110,112,113]]])
>>> c.shape
(2, 2, 3)
>>> c[1,...] # same as c[1,:,:] or c[1]
array([[100, 101, 102],
[110, 112, 113]])
>>> c[...,2] # same as c[:,:,2]
array([[ 2, 13],
[102, 113]])
>>> let f x y = 10*x+y
>>> let b = build (5,4) f :: Matrix Z
>>> b
(5><4)
[ 0, 1, 2, 3
, 10, 11, 12, 13
, 20, 21, 22, 23
, 30, 31, 32, 33
, 40, 41, 42, 43 ]
>>> b ! 2 ! 3
23
>>> flatten (subMatrix (0,1) (5,1) b) -- @subMatrix@ always return a Matrix
[1,11,21,31,41]
>>> flatten (b ?? (All, Pos (scalar 1)))
[1,11,21,31,41]
>>> b ?? (Range 1 1 2, All)
(2><4)
[ 10, 11, 12, 13
, 20, 21, 22, 23 ]
>>> flatten (b ? [rows b - 1])
[40,41,42,43]
-- hmatrix does not have 3D matrices; there are no ellipsis equivalents
>>> for row in b:
... print(row)
...
[0 1 2 3]
[10 11 12 13]
[20 21 22 23]
[30 31 32 33]
[40 41 42 43]
>>> for element in b.flat:
... print(element)
...
0
1
2
3
10
11
12
13
20
21
22
23
30
31
32
33
40
41
42
43
>>> let f x y = 10*x+y
>>> let b = build (5,4) f :: Matrix Z
>>> b
(5><4)
[ 0, 1, 2, 3
, 10, 11, 12, 13
, 20, 21, 22, 23
, 30, 31, 32, 33
, 40, 41, 42, 43 ]
>>> b ! 2 ! 3
23
>>> flatten (subMatrix (0,1) (5,1) b) -- @subMatrix@ always return a Matrix
[1,11,21,31,41]
>>> flatten (b ?? (All, Pos (scalar 1)))
[1,11,21,31,41]
>>> b ?? (Range 1 1 2, All)
(2><4)
[ 10, 11, 12, 13
, 20, 21, 22, 23 ]
>>> flatten (b ? [rows b - 1])
[40,41,42,43]
-- hmatrix does not have 3D matrices; there are no ellipsis equivalents