import sys
import math
import profile
import numarray as na
import numarray.linear_algebra as la

def pscircle(xy,r):
  x = xy[0]
  y = xy[1]
  print x,y+r,'moveto'
  print x+0.55*r, y+r,    x+r,y+0.45*r,    x+r,y,'curveto'
  print x+r,y-0.45*r,   x+0.55*r, y-r,   x,y-r,'curveto'
  print x-0.55*r,y-r,   x-r, y-0.45*r,   x-r,y,'curveto'
  print x-r,y+0.55*r,   x-0.55*r, y+r,   x,y+r,'curveto fill '

def psline(start, end):
  print start[0],start[1],'moveto'
  print end[0],end[1],'lineto'
  print '2 setlinewidth stroke'

def pscol(n):
  if (n < 0):
    print 0.1, 0.9+0.8*n, 0.1-0.8*n, 'setrgbcolor'
  else:
    print 0.1+0.8*n, 0.9-0.8*n, 0.1, 'setrgbcolor'

centre = (0,0)
ringrad = 1
ringdots = 6
rds = math.pi * 2 / ringdots
nrings = 12
sc = 0.8/nrings

def circlecent(ring, posn):
  if (ring != 0):
    return (centre[0] + ring*ringrad*math.cos(posn*rds/ring), centre[1] + ring*ringrad*math.sin(posn*rds/ring))
  else:
    return centre

# OK, this is a not-entirely-trivial sparse-matrix problem

# we start by setting up a labelling for where the circles are

circle_posn = [centre]
circle_labels = {}
circle_labels[(0,0)] = 0
circle_labelling = [(0,0)]
circle_number = 1
circle_isboundary = [False]
circles_nonboundary = 1

for rr in range(0, 1+nrings):
  for n in range(ringdots):
    for s in range(rr):
      dotnum = n*rr+s
      lab = (rr, dotnum)
      circle_labelling.append(lab)
      circle_labels[lab] = circle_number
      circle_posn.append(circlecent(rr,dotnum))
      circle_number = 1+circle_number
      if (rr == nrings):
        circle_isboundary.append(True)
      else:
        circle_isboundary.append(False)
        circles_nonboundary = 1+circles_nonboundary

sys.stderr.write(str(circle_number)+" circles handled\n")
sys.stderr.write(str(circles_nonboundary)+" of them are non-boundary\n")

neighbours = [[] for i in range(circle_number)]

def addlink(l1, l2):
  neighbours[l1].append(l2)
  neighbours[l2].append(l1)

for rr in range(0,1+nrings):
  for n in range(ringdots):
    if (rr != nrings):
      if (n != 0):
        addlink(circle_labels[(rr,n*rr)], circle_labels[(rr+1, n*(rr+1)-1)])
      else:
        addlink(circle_labels[(rr,n*rr)], circle_labels[(rr+1, ringdots*(rr+1)-1)])
    for s in range(rr):
      dotnum = n*rr+s
      if (dotnum != ringdots*rr-1):
        addlink(circle_labels[(rr,dotnum)], circle_labels[(rr, 1+dotnum)])
      else:
        addlink(circle_labels[(rr,dotnum)], circle_labels[(rr, 0)])
      if (rr != nrings):
        addlink(circle_labels[(rr,dotnum)], circle_labels[(rr+1, n*(rr+1)+s)])
        addlink(circle_labels[(rr,dotnum)], circle_labels[(rr+1, n*(rr+1)+s+1)])

# OK, that's set up the position set (the nodes, I think it's called)
# now for the matrix

nbc = [i for i in range(circle_number) if circle_isboundary[i] == False]
mat = na.zeros(shape=(circles_nonboundary, circles_nonboundary), type="Float64")

for node in nbc:
  nodep = circle_posn[node]
  # 1 x y x^2 xy y^2
  m0 = [[1,0,0,0,0,0]]
  for node2 in neighbours[node]:
    nodep2 = circle_posn[node2]
    dx = nodep2[0]-nodep[0]
    dy = nodep2[1]-nodep[1]
    line = [1, dx, dy, dx*dx, dx*dy, dy*dy]
    m0.append(line)
  # we want to get d/dy^2 and d/dx^2 out, and add them together
  m0 = na.array(m0)
  m0alg = na.matrixmultiply(m0,la.inverse(na.matrixmultiply(na.transpose(m0),m0)))
  # OK, if we multiply m0alg by a vector [val[node], val[neighbours[node]]]
  # we should get a six-vector estimating 1, d/dx, d/dy, d/dx^2, d/dxdy, d/dy^2
  m0alg = na.transpose(m0alg)
  
  coeffs = m0alg[3] + m0alg[5]

  places = [node] + neighbours[node]

  # recall that the boundary is fixed at zero
  for a in range(len(places)):
    if (circle_isboundary[places[a]] == False):
      mat[node, places[a]] = coeffs[a]

# that's produced the matrix, and truly it's a hideous thing
# do the eigenfunction computation

sys.stderr.write("Doing eigenvector calculation on "+str(len(nbc))+"^2 matrix\n\n")
profile.run("ev = la.eigenvectors(mat)", filename="lagrange.profile")
sys.stderr.write("done")

# extract and plot the solutions

sz = zip(ev[0],ev[1])
sz = [[abs(z[0]), z[1]] for z in sz if 
(z[1].imaginary * z[1].imaginary).sum() < 1.0e-10]


szi = zip([z[0] for z in sz],range(len(sz)))
szi.sort()
sz = [sz[z[1]] for z in szi]

solns = []

for res in sz:
  rr = res[1].real
  rm = max([abs(r) for r in rr.tolist()])
  rr = rr/rm
  solns.append(rr)

# paper is 210x297mm = 590 x 840 points

sys.stderr.write(str(len(solns))+" solutions found\n");

# with margins, we have a region of roughly
# 520 by 770, let's call it 500 x 750
# so we want `size` across and `3*size/2` down
# so we want size to be even, and 3*size**2/2 to be
# at least len(solns)
size = 2
while (3 * size * size / 2 < len(solns)):
  size = 1+size

e=0;

blocksize = 500.0/size
scalsize = (0.8 * blocksize) / (2*nrings+1)
dotsize = scalsize / 2.5

for s in solns:
  x0 = 45 + blocksize/2 + (e%size) * blocksize
  y0 = 45 + 750 - blocksize/2 - (e/size) * blocksize
  for cel in range(len(s)):
    pscol(s[cel])
    x,y = circle_posn[nbc[cel]]
    pscircle((x0+x*scalsize, y0+y*scalsize),dotsize)
  e = 1+e