%matplotlib inline
import numpy as np
import scipy as sp
import scipy.linalg as la
import matplotlib.pyplot as plt

import scipy.integrate as sint

A = np.array([[1, -20],
              [3, 4]], dtype=np.double)

def yprime(y, t):
    return sp.dot(A, y)

t = np.linspace(0, 1, 51)
yzero = np.array([1, 0])

# Solve equation using odeint
y = sint.odeint(yprime, yzero, t) # Calls the Fortran library odeint
plt.plot(y[:, 0], y[:, 1]) # Blue line is y(t) for t in [0, 1]

plt.xlim([-4, 9])
plt.ylim([-3, 5])

# Solve for 10 t values using the exponential
tvals = np.linspace(0, 1, 11)

for tval in tvals:
    yval = sp.dot(la.expm(A*tval), yzero)
    plt.plot(yval[0], yval[1], 'rs') # Red squares are evaluated via the matrix exponential


import networkx as nx

Adj = np.array([[0, 1, 1, 0, 0, 0], 
                [1, 0, 1, 1, 1, 0], 
                [1, 1, 0, 0, 0, 0], 
                [0, 1, 0, 0, 0, 1], 
                [0, 1, 0, 0, 0, 0],
                [0, 0, 0, 1, 0, 0]])

G = nx.from_numpy_matrix(Adj)
nx.draw(G, node_color='y', node_size=1000)

centralities = np.diag(la.expm(np.array(Adj, dtype=np.double)))
nodeorder = np.argsort(centralities)[::-1]

print np.array([nodeorder, centralities[nodeorder]])

# Note: This is already built into networkx using the following command
# print nx.communicability_centrality_exp(G)

sconst = 1

def TaylorSS(A):
    taylordegree = 15 # Use order 15 Taylor approximation
    s = np.ceil(sp.log2(la.norm(A))) + sconst # Find s such that norm(A/2^s) is small.
    X = A/(2**s)
    eX = np.eye(np.shape(A)[0])
    for k in range(taylordegree): # Compute the Taylor series
        eX = eX + X/sp.misc.factorial(k+1) 
        X = sp.dot(X, X)
        
    for k in range(np.int64(s)):
        eX = sp.dot(eX, eX) # Do the squaring phase of the algorithm
        
    return eX

#Let's test it against la.expm
A = np.random.randn(4, 4)
E1 = TaylorSS(A)
Eexact = la.expm(A)
print la.norm(E1 - Eexact)/la.norm(Eexact) # Relative error

relerrs = np.zeros(15)
for sconst in range(15):
    E = TaylorSS(A)
    relerrs[sconst] = la.norm(E - Eexact)/la.norm(Eexact)
    
plt.plot(relerrs, color='blue', lw=2)
plt.yscale('log')
plt.xlabel('sconst', fontsize=12)
plt.ylabel('rel. err.', fontsize=12)