sym_gs_eigs.py

from __future__ import division
from sympy import MatrixSymbol, Matrix, latex, simplify
import baker

def SymmetricMatrix(name,n):
    J = Matrix(MatrixSymbol(name,n,n))
    for i in range(n):
        for j in range(i):
            J[i,j] = J[j,i]
    return J

def split(J):
    n = J.shape[0]
    return Matrix([[(J[i,j] if j <= i else 0) for j in range(n)] for i in range(n)]), \
            Matrix([[(-J[i,j] if j > i else 0) for j in range(n)] for i in range(n)])

def update(J):
    L, U = split(J)
    return L.inv() * U

def update_evals(J):
    return update(J).eigenvals()

@baker.command(default=True)
def main(n=3):
    dct = update_evals(SymmetricMatrix('J',n))
    for valexpr, count in dct.iteritems():
        print count
        print latex(simplify(valexpr))
        print

@baker.command
def somezeros():
    J = SymmetricMatrix('J',3)
    J[0,2] = J[2,0] = 0
    dct = update_evals(J)
    for valexpr, count in dct.iteritems():
        print count
        print latex(simplify(valexpr))
        print

@baker.command
def somezeros2():
    J = SymmetricMatrix('J',4)
    J[1,1] = J[2,2] = J[3,3] = J[0,0]
    J[0,1] = J[1,0] = J[1,2] = J[2,1] = J[2,3] = J[3,2] = J[0,3] = J[3,0]
    J[0,2] = J[2,0] = J[3,1] = J[1,3] = 0
    dct = update_evals(J)
    for valexpr, count in dct.iteritems():
        print count
        print latex(simplify(valexpr))
        print

if __name__ == '__main__':
    baker.run()

contour_plot.py

from __future__ import division
from pylab import *
from blah import spectrum

np.random.seed(0)

def f(s):
    return 2*s/(1+s)

t = np.linspace(-1,1,2000,endpoint=True)
X,Y = meshgrid(t,t)
pts = (X + 1j*Y)
pts[np.abs(pts) > 1] = 0.

figure(figsize=(8,8))
subplot(111,aspect='equal')
imshow(np.log(np.abs(f(pts))),extent=(-1,1,-1,1))
CS = contour(X,Y,np.log(np.abs(f(pts))),np.log((0.5,1.,2.,3.,4.,5.)),colors='k')
manual_locations = [(-0.1,0.2),(-0.2,0.38),(-0.5,0.),(-0.7,-0.53),(-0.8,-0.38),(-0.85,-0.34)]
clabel(CS,fontsize=12,inline=1,fmt=lambda x: str(np.round(np.exp(x),3)),manual=manual_locations)
grid('off')
axis('off')

pairs = spectrum(n=50,plot=False)
plt.plot(np.real(pairs.flat),np.imag(pairs.flat),'kx')

# pairs = spectrum(n=50,k=2,eps=100.,plot=False)
# plt.plot(np.real(pairs.flat),np.imag(pairs.flat),'rx')

show()

gs_eigs.py

from __future__ import division
import baker
import numpy as np
from numpy import *
from numpy.random import *
from matplotlib import pyplot as plt

def rand_psd(n,k=None,eps=0.):
    k = k if k is not None else n
    A = randn(n,k)
    return A.dot(A.T) + np.eye(n)*eps

### gs utils

def gs_split(A):
    return np.tril(A), np.triu(A,k=1)

def gs_update(A):
    L, U = gs_split(A)
    return np.linalg.solve(L,-U)

def gs_spectrum(A):
    G = gs_update(A)
    return np.linalg.eig(G)

### experiments

@baker.command(default=True)
def spectrum(n=50,k=-1,eps=0.0):
    k = k if k != -1 else n
    ws, Vs = zip(*[gs_spectrum(rand_psd(n,k)) for i in range(2)])
    pairs = np.outer(*ws)
    plt.plot(np.real(pairs.flat),np.imag(pairs.flat),'b.')

if __name__ == '__main__':
    baker.run()
    plt.show()