1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
|
import numeric as N
import types
import string as str_
# make translation table
_table = [None]*256
for k in range(256):
_table[k] = chr(k)
_table = ''.join(_table)
_numchars = str_.digits + ".-+jeEL"
del str_
_todelete = []
for k in _table:
if k not in _numchars:
_todelete.append(k)
_todelete = ''.join(_todelete)
def _eval(astr):
return eval(astr.translate(_table,_todelete))
def _convert_from_string(data):
rows = data.split(';')
newdata = []
count = 0
for row in rows:
trow = row.split(',')
newrow = []
for col in trow:
temp = col.split()
newrow.extend(map(_eval,temp))
if count == 0:
Ncols = len(newrow)
elif len(newrow) != Ncols:
raise ValueError, "Rows not the same size."
count += 1
newdata.append(newrow)
return newdata
_lkup = {'0':'000',
'1':'001',
'2':'010',
'3':'011',
'4':'100',
'5':'101',
'6':'110',
'7':'111'}
def _binary(num):
ostr = oct(num)
bin = ''
for ch in ostr[1:]:
bin += _lkup[ch]
ind = 0
while bin[ind] == '0':
ind += 1
return bin[ind:]
class matrix(N.ndarray):
__array_priority__ = 1.0
def __new__(self, data, dtype=None, copy=0):
if isinstance(data, matrix):
dtype2 = data.type
if (dtype is None):
dtype = dtype2
if (dtype2 is dtype) and (not copy):
return data
return data.astype(dtype)
if dtype is None:
dtype = N.intp
intype = N.totype(dtype)
if isinstance(data, types.StringType):
data = _convert_from_string(data)
# now convert data to an array
arr = N.array(data, dtype=intype, copy=copy)
ndim = arr.ndim
shape = arr.shape
if (ndim > 2):
raise ValueError, "matrix must be 2-dimensional"
elif ndim == 0:
shape = (1,1)
elif ndim == 1:
shape = (1,shape[0])
fortran = False;
if (ndim == 2) and arr.flags['FORTRAN']:
fortran = True
if not (fortran or arr.flags['CONTIGUOUS']):
arr = arr.copy()
ret = N.ndarray.__new__(matrix, shape, intype, buffer=arr, fortran=fortran,
swapped=(not arr.flags['NOTSWAPPED']))
return ret;
def __array_wrap__(self, obj):
try:
ret = matrix(obj,typecode=obj.type)
except:
ret = obj
return ret
def _update_meta(self, obj):
ndim = self.ndim
if ndim == 0:
self.shape = (1,1)
elif ndim == 1:
self.shape = (1, self.shape[0])
return
def __getitem__(self, index):
out = N.ndarray.__getitem__(self, index)
# Need to swap if slice is on first index
try:
n = len(index)
if (n > 1) and isinstance(index[0], types.SliceType) \
and (isinstance(index[1], types.IntType) or
isinstance(index[1], types.LongType)):
sh = out.shape
out.shape = (sh[1], sh[0])
except TypeError:
pass
return out
def __mul__(self, other):
return N.dot(self, other)
def __rmul__(self, other):
return N.dot(other, self)
def __pow__(self, other):
if len(shape)!=2 or shape[0]!=shape[1]:
raise TypeError, "matrix is not square"
if type(other) in (type(1), type(1L)):
if other==0:
return matrix(N.identity(shape[0]))
if other<0:
x = self.I
other=-other
else:
x=self
result = x
if other <= 3:
while(other>1):
result=result*x
other=other-1
return result
# binary decomposition to reduce the number of Matrix
# Multiplies for other > 3.
beta = _binary(other)
t = len(beta)
Z,q = x.copy(),0
while beta[t-q-1] == '0':
Z *= Z
q += 1
result = Z.copy()
for k in range(q+1,t):
Z *= Z
if beta[t-k-1] == '1':
result *= Z
return result
else:
raise TypeError, "exponent must be an integer"
def __rpow__(self, other):
raise NotImplementedError
def getA(self):
arr = self
fortran = False;
if (self.ndim == 2) and arr.flags['FORTRAN']:
fortran = True
if not (fortran or arr.flags['CONTIGUOUS']):
arr = arr.copy()
return N.ndarray.__new__(N.ndarray, self.shape, self.type, buffer=arr,
fortran=fortran,
swapped=(not arr.flags['NOTSWAPPED']))
def getT(self):
return N.transpose(self)
def getH(self):
return N.conjugate(N.transpose(self))
# inverse doesn't work yet....
def getI(self):
return self
A = property(getA, None, doc="Get base array")
T = property(getT, None, doc="transpose")
H = property(getH, None, doc="hermitian (conjugate) transpose")
|
