#!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
 import numpy, random
 
 ## Hilfsprozeduren
 
 # Konstruiert aus dem Hauptdiagonalenvektor alphas und dem unteren
 # Nebendiagonalvektor betas die zugehörige hermitische n*n Tridiagonalmatrix.
 # Als obere Nebendiagonale wird die komplexe Konjugation der übergebenen
 # unteren Nebendiagonale verwendet.
 # Es wird in jedem Fall eine Matrix mit komplexen Einträgen zurückgeliefert,
 # auch, wenn die Eingabevektoren reell waren.
 # len(alphas) = n, len(betas) = n - 1
 def tridiag(alphas, betas):
     if len(alphas) != len(betas) + 1:
 	raise "..."
     n = len(alphas)
 
     matrix = numpy.zeros((n,n), dtype=complex)
 
     for i in range(0,n):
 	matrix[i,i] = alphas[i]
 	if i < n - 1:
 	    matrix[i,i+1] = complex(betas[i]).conjugate()
 	    matrix[i+1,i] = complex(betas[i])
     return matrix
 
 # Wie numpy.linalg.eigvalsh, nur werden die Eigenwerte als reelle float-Objekte
 # zurückgegeben, statt unnötigerweise (und Sortierung verhindernde) komplexe
 # Zahlen.
 def eigvalsh(matrix):
     return map(lambda x: x.real, numpy.linalg.eigvalsh(matrix))
 
 # (Reelle) Multiplikation einer spärlich besetzten Matrix (im
 # Dictionary-Format) mit einem (numpy-)Vektor.
 def sparse_dot(A, v, N):
     assert len(v) == N
 
     w = numpy.array([0.0] * N)
     for (i,j) in A:
 	w[i] += A[(i,j)] * v[j]
 
     return w
 
 # Länge, ab der wir einen Vektor als Nullvektor ansehen und das Verfahren
 # abbrechen
 epsilon = 0.0001
 
 # Hilfsfunktion: Nimmt eine numpy-Matrix und gibt sie als spärlich besetzte
 # Matrix im Dictionary-Format zurück.
 def dense_to_sparse(A, eps=epsilon):
     mat = {}
 
     for i in range(A.shape[0]):
 	for j in range(A.shape[1]):
 	    if(abs(A[i,j]) >= eps):
 		mat[(i,j)] = A[i,j]
     
     return mat
 
 def extremal_eigenvalues(H, max_iterations):
     assert H.shape[0] == H.shape[1]
     N = H.shape[0]
 
     return extremal_eigenvalues_sparse(dense_to_sparse(H), N, max_iterations)
 
 def extremal_eigenvalues_sparse(H, N, max_iterations):
     # Zufälliger normierter Startvektor
     v0  = numpy.array([ random.random() + 1 for i in range(0,N) ])
     v0 /= numpy.linalg.norm(v0)
 
     # Beginn des Lanczos-Algorithmus.
     # Initialisierung der Variablen:
     # vs:     Liste der erzeugten Orthonormalbasis des Krylowraums
     #         span { v0, H v0, H^2 v0, H^3 v0, ... }
     # alphas: Hauptdiagonalelemente in der Darstellung von H bzgl. vs
     # betas:  Untere = obere Nebendiagonalelemente (rein reell)
     vs     = [v0]
     alphas = []
     betas  = [0]
 
     # Iterationsschleife
     for j in range(N):
 	w  = sparse_dot(H, vs[j], N)
 	if j > 0:
 	    w -= betas[j] * vs[j-1]
 	alphas.append(numpy.dot(vs[j], w))
 	w -= alphas[j] * vs[j]
 
 	# Mehr als N Basisvektoren kann es nicht geben, wird sind fertig.
 	if j == N - 1 or j == max_iterations - 1: break
 
 	# Haben wir es geschafft, einen vom Nullvektor echt verschiedenen
 	# Basisvektor zu finden?
 	beta = numpy.linalg.norm(w)
 	if beta < epsilon:
 	    print "[LANCZOS] Norm < kleiner epsilon: %.2f, j=%d" % (beta,j)
 	    break
 	else:
 	    betas.append(beta)
 	    vs.append(w / beta)
 
     H_ = tridiag(alphas,betas[1:len(betas)])
     eigenvalues = eigvalsh(H_)
 
     return eigenvalues[0], eigenvalues[-1]