[frames] | no frames]

Source Code for Module gmisclib.matrix_arrange_entropy

``` 1
2
3  """This rotates a matrix by multiplying with a unitary matrix
4  so that the resulting elements are either nearly zero or relatively large.
5  It minimizes sum of |x_ij|*log(|x_ij|), i.e. the entropy (sortof).
6  """
7
8  import math
9  from gmisclib import Num
10  from gmisclib import mcmc
11  from gmisclib import mcmc_helper
12  from gmisclib import gpkmisc
13
14
15 -def make_min_entropy(x, extra_entropy=0.0):
16          """extra_entropy is a vector that helps choose what to optimize."""
17          x = Num.asarray(x, Num.Float)
18          print 'x.shape=', x.shape
19          ndim, nvec = x.shape
20
21          def unitary(p):
22                  assert p.shape[0] == nvec**2
23                  pm = Num.reshape(p, (nvec, nvec))
24                  psym = Num.dot(pm, Num.transpose(pm))
25                  evalues, umat = Num.LA.eigh(psym)
26                  return umat
27
28          def rotated(xt, p):
29                  frot = Num.dot(unitary(p), xt)
30                  assert frot.shape == xt.shape
31                  return frot
32
33          def fom(p, c):
34                  xt, nvec, extra_entropy, G = c
35                  frot = rotated(xt, p)
36                  assert frot.shape[0] == nvec
37                  # print 'frot=', frot
38                  frotfom = Num.absolute(frot)
39                  frotn = frotfom/Num.sum(frotfom, axis=1)[:,Num.NewAxis]
40                  assert abs(Num.sum(frotn[0])-1.0) < 0.01
41                  nege = Num.sum(frotn*(Num.log(frotn)-extra_entropy), axis=1)
42                  assert nege.shape == (nvec,)
43                  print 'fom=', Num.sum(nege)
44                  return G*Num.sum(nege)
45
46          def fixer(p, c):
47                  Num.divide(p, math.sqrt(Num.average(Num.square(p))), p)
48                  return p
49
50          prmvec = Num.RA.normal(0.0, 1.0, (nvec**2,))
51          stpr = mcmc.bootstepper(fom, [prmvec], v=gpkmisc.make_diag(Num.square(prmvec)),
52                                          c = (Num.transpose(x), nvec, extra_entropy, 20.0),
53                                          fixer=fixer
54                                          )
55          mcmch = mcmc_helper.stepper(stpr)
56          mcmch.run_to_bottom()
57          stpr.T = 0.3
58          mcmch.run_to_bottom()
59          stpr.T = 0.1
60          mcmch.run_to_bottom()
61          prms = stpr.current().prms()
62          return (Num.transpose(rotated(Num.transpose(x), prms)),
63                  unitary(prms)
64                  )
65
66
67 -def test2():
68          x = [[1.0, 1.0], [-1.0, 1.0]]
69          y,u = make_min_entropy(x)
70          assert abs(Num.sum(u**2) - u.shape[0]) < 0.001
71          print 'y=', y
72          ay = Num.absolute(y)
73          assert Num.sum(Num.less(ay, 0.001)) == 2
74          r2 = math.sqrt(2.0)
75          assert Num.sum(Num.less(Num.absolute(ay-r2), 0.001)) == 2
76
77 -def test3():
78          x = [[0, 0, 1], [1.0, -1.0, 0], [1.0, 1.0, 0]]
79          y,u = make_min_entropy(x)
80          assert abs(Num.sum(u**2) - u.shape[0]) < 0.001
81          print 'y=', y
82          ay = Num.absolute(y)
83          assert Num.sum(Num.less(ay, 0.01)) == 6
84          r2 = math.sqrt(2.0)
85          assert Num.sum(Num.less(Num.absolute(ay-r2), 0.01)) == 2
86          assert Num.sum(Num.less(Num.absolute(ay-1), 0.01)) == 1
87
88
89  if __name__ == '__main__':
90          test2()
91          test3()
92
<!--
expandto(location.href);
// -->

```

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