[frames] | no frames]

# Source Code for Module gmisclib.matrix_arrange

```  1  """Take a covariance-like matrix, and re-order it to be close to a diagonal matrix.
2  We calculate a map that renames the indices (the same map for both left and right),
3  to bring large off-diagonal elements close to the diagonal.
4  """
5
6  import math
7  import random
8  import numpy
9
10
11 -def _eval_map(a, m, wt):
12          """This is a measure of how far the matrix is from diagonal,
13          after application of the mapping.
14          Matrix "a" is a numeric python 2-d array.
15          Map is a numeric python 1-d integer array."""
16
17          n = a.shape[0]
18          assert a.shape == (n, n)
19          # print "mapped a=", map_array(a, m)
20          # print "m=", m
21          return numpy.sum(map_array(a, m)*wt)
22
23
24 -def _offdiag(n):
25          index = numpy.arange(n)
26          offdiag = numpy.absolute(index[:,numpy.newaxis]-index[numpy.newaxis,:])
27          assert offdiag.shape == (n, n)
28          assert offdiag[0, 0] == 0
29          assert offdiag[1, 0] == 1
30          # print "WT=", offdiag
31          return offdiag
32
33
34
35 -def _optimistic(n):
36          return numpy.ones((n, n), numpy.int) - numpy.identity(n)
37
38
39
40 -def diagonalize(a0):
41          EPS = 1e-10
42          a = numpy.array(a0, copy=True)
43          n = a.shape[0]
44          assert a.shape == (n, n)
45          wt = numpy.square(_offdiag(n))
46          m = numpy.arange(n)
47          opval = _eval_map(a, m, _optimistic(n))
48          # print "OPVAL=", opval
49          bad = _eval_map(a, m, wt)
51          f = 1.0
52          fr = math.pow(0.5, 1.0/n)
53          nobetter = 0
54          while 1:
55                  # i = random.choice(_worst(a, m, wt))
56                  i = random.randrange(n)
57                  j = random.randrange(n)
58                  if j == i:
59                          continue
60                  # print "TESTING", i, j
61                  tmp = m[i]
62                  m[i] = m[j]
63                  m[j] = tmp
64                  bt = _eval_map(a, m, wt)
65                  if bt <= opval + bad*EPS:
67                          break
68                  if bt <= bad * (1 + f/float(n)):
69                          # print "ACCEPT", bad, bt
71                                  f *= fr
72                                  nobetter = 0
74                          nobetter = 0
75                  else:
76                          # print "NG", bad, bt
77                          tmp = m[i]
78                          m[i] = m[j]
79                          m[j] = tmp
80                          nobetter += 1
81
82                  if nobetter >= n*n:
83                          break
84          if m[0] > m[-1]:        # Pick one of two equivalent orderings
85                  m = m[::-1]
87
88
89
90 -def map_array(a, m):
91          return numpy.take(numpy.take(a, m, axis=0), m, axis=1)
92
93
94 -def test():
95          a = numpy.array([[1, 0, 3, 2], [0, 1, 0, 0], [3, 0, 1, 0], [2, 0, 0, 1]])
96          m, b0, b1 = diagonalize(a)
97          assert b1<b0
98          print m
99          print map_array(a, m)
100
101
102  if __name__ == '__main__':
103          test()
104
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:14 2011 http://epydoc.sourceforge.net