import sys
# Diese Zeile muss angepasst werden!
sys.path.append("/home/moser/MG_2016/pyMG/")
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import pymg
from project.helmholtz1d import Helmholtz1D
from project.helmholtz1d_periodic import Helmholtz1D_Periodic
from project.gauss_seidel import GaussSeidel
from project.weighted_jacobi import WeightedJacobi
from project.pfasst.plot_tools import eigvalue_plot_list, matrix_plot, matrix_row_plot
from project.pfasst.transfer_tools import to_dense
from project.pfasst.matrix_method_tools import matrix_power

def plot_3_eigvalueplots(A_p,A_z,A_m):
    eig_p.append(sp.linalg.eigvals(to_dense(A_p)))
    eig_z.append(sp.linalg.eigvals(to_dense(A_z)))
    eig_m.append(sp.linalg.eigvals(to_dense(A_m)))
    
    real_part_p = np.real(eig_p[-1])
    img_part_p = np.imag(eig_p[-1])
    
    real_part_z = np.real(eig_z[-1])
    img_part_z = np.imag(eig_z[-1])
    
    real_part_m = np.real(eig_m[-1])
    img_part_m = np.imag(eig_m[-1])
    fig1, (ax1, ax2, ax3) = plt.subplots(ncols=3,figsize=(15,3))
    
    ax1.plot(real_part_p,img_part_p,'ro')
    ax1.set_xlabel("real part")
    ax1.set_ylabel("img part")
    ax1.set_title('eigenvalues')
    
    ax2.plot(real_part_z,img_part_z,'bo')
    ax2.set_xlabel("real part")
    ax2.set_ylabel("img part")
    ax2.set_title('eigenvalues')
    
    ax3.plot(real_part_m,img_part_m,'go')
    ax3.set_xlabel("real part")
    ax3.set_ylabel("img part")
    ax3.set_title('eigenvalues')
    
    fig1.tight_layout()
    plt.show()

def plot_2_eigvalueplots(A_p,A_z):
    eig_p.append(sp.linalg.eigvals(to_dense(A_p)))
    eig_z.append(sp.linalg.eigvals(to_dense(A_z)))

    
    real_part_p = np.real(eig_p[-1])
    img_part_p = np.imag(eig_p[-1])
    
    real_part_z = np.real(eig_z[-1])
    img_part_z = np.imag(eig_z[-1])
    
    fig1, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,3))
    
    ax1.plot(real_part_p,img_part_p,'ro')
    ax1.set_xlabel("real part")
    ax1.set_ylabel("img part")
    ax1.set_title('eigenvalues')
    
    ax2.plot(real_part_z,img_part_z,'bo')
    ax2.set_xlabel("real part")
    ax2.set_ylabel("img part")
    ax2.set_title('eigenvalues')
    
    fig1.tight_layout()
    plt.show()

"""
Explanation: Interaktives Übungsblatt
Vorgeplänkel
End of explanation
"""


def system_matrix_hh1d(n,sig):
    hh1d = Helmholtz1D(n, sig)
    return hh1d.A

def system_matrix_hh1d_periodic(n,sig):
    hh1d = Helmholtz1D_Periodic(n, sig)
    return hh1d.A

def spec_rad(A):
    return np.max(np.abs(sp.linalg.eigvals(to_dense(A))))


"""
Explanation: Systemmatrizen
Wir werden im Weiteren pyMG nutzen um die Systemmatrix für gegebene Parameter $ n$ und $\sigma$ für das Helmholtz-Problem in 1D aufzustellen.
End of explanation
"""


matrix_plot(to_dense(system_matrix_hh1d(10,0)))
matrix_plot(to_dense(system_matrix_hh1d_periodic(10,0)))

"""
Explanation: Plotten Sie mithilfe von matrix_plot die Systemmatrizen für $\sigma = 0$ und $n=10$.
End of explanation
"""


n=30
for sigma in [1000,0,-1000]:
    plot_2_eigvalueplots(system_matrix_hh1d(n,sigma),system_matrix_hh1d_periodic(n,sigma))

"""
Explanation: Aufgabe: Plotten Sie mithilfe von plot_3_eigvalueplots die Eigenwerte der Systemmatrix für $n \in [5,10,20]$ und $\sigma = 100$,$\sigma = -100$ und $\sigma = 0$.
Frage: Wie unterscheiden sich die Spektren der verschiedenen Systemmatrizen?
End of explanation
"""


def iteration_matrix_wjac(n, sigma, periodic=True):
    if periodic:
        A = system_matrix_hh1d_periodic(n,sigma)
    else:
        A = system_matrix_hh1d(n,sigma)
        
    wjac = WeightedJacobi(A, 2.0/3.0)
    P_inv = wjac.Pinv
    return np.eye(n) - P_inv.dot(A)

"""
Explanation: Iterationsmatrizen des Glätters
Frage
Weitaus spannender sind die Spektralradiien der Iterationsmatrizen eines Glätters. Warum?
End of explanation
"""


matrix_plot(iteration_matrix_wjac(10,-100))

"""
Explanation: Aufgabe
Plotten Sie die Iterationsmatrix des gewichteten Jacobi für verschiedene $\sigma$. Zu welcher Klasse gehört diese Iterationsmatrix.
End of explanation
"""


n = 10
sigma_range = np.linspace(-100,100,100)
sr_wjac_periodic = map(lambda sig : spec_rad(iteration_matrix_wjac(n, sig,periodic=True)), sigma_range)
sr_wjac = map(lambda sig : spec_rad(iteration_matrix_wjac(n, sig,periodic=False)), sigma_range)

# Achsen festhalten

fig1, (ax1, ax2, ax3) = plt.subplots(ncols=3,figsize=(15,4))

ax1.plot(sigma_range, sr_wjac_periodic,'k-')
ax1.set_xlabel('$\sigma$')
ax1.set_ylabel("spectral radius")
ax1.set_title('periodic')

ax2.plot(sigma_range, sr_wjac,'k-')
ax2.set_xlabel('$\sigma$')
ax2.set_ylabel("spectral radius")
ax2.set_title('non-periodic')

ax3.plot(sigma_range, np.abs(np.asarray(sr_wjac) - np.asarray(sr_wjac_periodic)),'k-')
ax3.set_xlabel('$\sigma$')
ax3.set_ylabel("spectral radius")
ax3.set_title('difference')

fig1.tight_layout()
plt.show()

"""
Explanation: Gehen wir das ganze mal systematischer an!
End of explanation
"""


def iteration_matrix_gs(n, sigma, periodic=True):
    
    if periodic:
        A = system_matrix_hh1d_periodic(n,sigma)
    else:
        A = system_matrix_hh1d(n,sigma)
        
    gs = GaussSeidel(A)
    P_inv = gs.Pinv
    return np.eye(n) - P_inv.dot(A)

"""
Explanation: Frage : Wie verhalten sich die Spektren für das periodische Problem zu den Problemen mit Dirichletrandbedingungen für verschiedene $\sigma$ und $n$? Erkenntnis durch ausprobieren!
Aufgabe
Nutzen Sie die folgende Funktion, um die Iterationsmatrix für Gauß-Seidel abhängig von $\sigma$ und $n$ zu berechne. Finden Sie heraus wie sich der Spektralradius für verschiedene $\sigma$ und den periodischen, sowie nicht periodischen Fall verhält.
End of explanation
"""


matrix_plot(iteration_matrix_gs(10,0,True))

"""
Explanation: Frage
Ist die Iterationsmatrix zyklisch?
End of explanation
"""


sr_gs_periodic = map(lambda sig : spec_rad(iteration_matrix_gs(n, sig,periodic=True)), sigma_range)
sr_gs = map(lambda sig : spec_rad(iteration_matrix_gs(n, sig,periodic=False)), sigma_range)


fig1, (ax1, ax2, ax3) = plt.subplots(ncols=3,figsize=(15,4))

ax1.plot(sigma_range, sr_gs_periodic,'k-')
ax1.set_xlabel('$\sigma$')
ax1.set_ylabel("spectral radius")
ax1.set_title('periodic')

ax2.plot(sigma_range, sr_gs,'k-')
ax2.set_xlabel('$\sigma$')
ax2.set_ylabel("spectral radius")
ax2.set_title('non-periodic')

ax3.plot(sigma_range, np.abs(np.asarray(sr_gs) - np.asarray(sr_gs_periodic)),'k-')
ax3.set_xlabel('$\sigma$')
ax3.set_ylabel("spectral radius")
ax3.set_title('difference')

fig1.tight_layout()
plt.show()

"""
Explanation: Nochmal das Ganze mit mehr System!
End of explanation
"""


def transformation_matrix_fourier_basis(N):
    psi = np.zeros((N,N),dtype=np.complex128)
    for i in range(N):
        for j in range(N):
            psi[i,j] = np.exp(2*np.pi*1.0j*j*i/N)
    return psi/np.sqrt(N)

def plot_fourier_transformed(A):
    A = to_dense(A)
    n = A.shape[0]
    PSI_trafo = transformation_matrix_fourier_basis(n)
    PSI_trafo_inv =  sp.linalg.inv(PSI_trafo)
    A_traf = np.dot(PSI_trafo_inv, np.dot(A,PSI_trafo))
    matrix_row_plot([A,np.abs(A_traf)])

plot_fourier_transformed(iteration_matrix_wjac(16,0))
plot_fourier_transformed(iteration_matrix_gs(16,0))

"""
Explanation: Frage : Wie verhalten sich die Spektren für das periodische Problem zu den Problemen mit Dirichletrandbedingungen für verschiedene $\sigma$ und $n$? Erkenntnis durch ausprobieren!
Das Leben im Fourier-Raum
Wir können mithilfe der Moden
$$v^{(m)} = \frac{1}{\sqrt{n}} \begin{pmatrix}1 \ e^{-2\pi i m/n} \ \vdots \  e^{-2\pi i m(n-1)/n} \end{pmatrix} $$
eine Transformation definieren. Die uns den Operatoren/Matrizen in den Fourier-Raum übersetzt.
End of explanation
"""


def get_theta_eigvals(A, plot=False,which='all'):
    A = to_dense(A)
    n = A.shape[0]
    PSI_trafo = transformation_matrix_fourier_basis(n)
    PSI_trafo_inv =  sp.linalg.inv(PSI_trafo)
    A_traf = np.dot(PSI_trafo_inv, np.dot(A,PSI_trafo))
    if plot:
        matrix_plot(np.abs(A_traf))
    eigvals =  np.asarray(map(lambda k : A_traf[k,k],range(n)))
    if which is 'high':
        return eigvals[np.ceil(n/4):np.floor(3.0*n/4)]
    elif which is 'low':
        return np.hstack([eigvals[:np.floor(n/4)],eigvals[np.ceil(3.0*n/4):]])
    else:
        return eigvals

"""
Explanation: Frage: Was ist hier passiert? Und was passiert für unterschiedliche $\sigma$?
Die hohen Eigenwerte extrahiert man nun durch auspicken der richtigen Diagonalwerte nach der Transformation, falls die Matrix zyklisch ist.
End of explanation
"""


It_gs = iteration_matrix_gs(16,0)
eigvals = sp.linalg.eigvals(It_gs)
diagonals = get_theta_eigvals(It_gs)

"""
Explanation: Frage
Verhält sich der gewichtete Jacobi Glätter wie in der Vorlesung vorraussgesagt?
Frage: 
Wie gut passen eigentlich die Diagonalwerte der Fourier-transformierten Iterationsmatrix mit den Eigenwerten der Matrix für Gauss-Seidel zusammen? Was könnte man machen um Sie zu vergleichen.
End of explanation
"""


sum_eig = np.sum(np.abs(eigvals))
sum_diag = np.sum(np.abs(diagonals))
print sum_eig
print sum_diag

"""
Explanation: Eine einfache Methode wäre:
End of explanation
"""


from project.linear_transfer import LinearTransfer
from project.linear_transfer_periodic import LinearTransferPeriodic

"""
Explanation: Aufgabe
Überlegt euch eigene Vergleichsmethoden und variiert $n$,$\sigma$ und die Periodizität, um herauszufinden, wann die 
get_theta_eigvals Methode die Eigenwerte gut schätzen kann.
Zweigitter-Iterationsmatrix
End of explanation
"""


def coarse_grid_correction(n,nc, sigma):
    A_fine = to_dense(system_matrix_hh1d(n,sigma))
    A_coarse = to_dense(system_matrix_hh1d(nc,sigma))
    A_coarse_inv = sp.linalg.inv(A_coarse)
    lin_trans = LinearTransfer(n, nc)
    prolong = to_dense(lin_trans.I_2htoh)
    restrict = to_dense(lin_trans.I_hto2h)
    return np.eye(n)- np.dot(prolong.dot(A_coarse_inv.dot(restrict)), A_fine)

"""
Explanation: Im Folgenden werden wir nun mithilfe des pymg frameworks die Zweigitter-Iterationsmatrix für ein einfaches Multigrid 
aufstellen. Wir beginnen mit der Grobgitterkorrektur.
End of explanation
"""


plot_fourier_transformed(coarse_grid_correction(31,15,0))

"""
Explanation: Buntebilderaufgabe: Nutze plot_fourier_transformed um für $n=31$, $n_c=15$ und verschiedene $\sigma\in[-1000,1000]$ um die Grobgitterkorrekturiterationsmatrizen und deren Fourier-transformierten zu plotten.
End of explanation
"""


def coarse_grid_correction_periodic(n,nc, sigma):
    A_fine = to_dense(system_matrix_hh1d_periodic(n,sigma))
    A_coarse = to_dense(system_matrix_hh1d_periodic(nc,sigma))
    A_coarse_inv = sp.linalg.inv(A_coarse)
    lin_trans = LinearTransferPeriodic(n, nc)
    prolong = to_dense(lin_trans.I_2htoh)
    restrict = to_dense(lin_trans.I_hto2h)
    return np.eye(n)- np.dot(prolong.dot(A_coarse_inv.dot(restrict)), A_fine)

"""
Explanation: Und nun einmal der periodische Fall.
End of explanation
"""


def two_grid_it_matrix(n,nc, sigma, nu1=3,nu2=3,typ='wjac'):
    cg = coarse_grid_correction(n,nc,sigma)
    if typ is 'wjac':
        smoother = iteration_matrix_wjac(n,sigma, periodic=False)
    if typ is 'gs':
        smoother = iteration_matrix_gs(n,sigma, periodic=False)
    
    pre_sm = matrix_power(smoother, nu1)
    post_sm = matrix_power(smoother, nu2)
    
    return pre_sm.dot(cg.dot(post_sm))    

"""
Explanation: Aufgabe: 
Nutzen Sie coarse_grid_correction_periodic für die Grobgitterkorrektur für das periodische Problem und plotten Sie nochmal für verschiedene $\sigma$ die Matrizen und ihre Fourier-transformierten.Für welche $n_f$ und $n_c$ ist die Grobgitterkorrektur für ein periodisches Problem sinnvoll?
Frage: 
Was genau passiert bei $\sigma = 0$ und in der Nähe davon?  Was fällt sonst auf? Und was hat das Ganze mit der Vorlesung zu tun?
Nun verwenden wir die Grobgitterkorrektur und die Iterationsmatrizen der Glätter um die Zweigitteriterationsmatrix zu berrechnen.
End of explanation
"""


plot_fourier_transformed(two_grid_it_matrix(15,7,0,typ='wjac'))

"""
Explanation: Buntebilderaufgabe:
Nutzen Sie plot_fourier_transformed um für $n=15$, $n_c=7$ und verschiedene $\sigma\in[-1000,1000]$ um die Zweigittermatrizen und deren Fourier-transformierten zu plotten.
End of explanation
"""


sr_2grid_var_sigma = map(lambda sig : spec_rad(two_grid_it_matrix(15,7,sig)), sigma_range)

plt.semilogy(sigma_range, sr_2grid_var_sigma,'k-')
plt.title('$n_f = 15, n_c = 7$')
plt.xlabel('$\sigma$')
plt.ylabel("spectral radius")

nf_range = map(lambda k: 2**k-1,range(3,10))
nc_range = map(lambda k: 2**k-1,range(2,9))
sr_2grid_m1000 = map(lambda nf,nc : spec_rad(two_grid_it_matrix(nf,nc,-1000)), nf_range, nc_range)
sr_2grid_0 = map(lambda nf,nc : spec_rad(two_grid_it_matrix(nf,nc,0)), nf_range, nc_range)
sr_2grid_p1000 = map(lambda nf,nc : spec_rad(two_grid_it_matrix(nf,nc,1000)), nf_range, nc_range)

plt.semilogy(nf_range, sr_2grid_m1000,'k-',nf_range, sr_2grid_0,'k--',nf_range, sr_2grid_p1000,'k:')
plt.xlabel('$n_f$')
plt.ylabel("spectral radius")
plt.legend(("$\sigma = -1000$","$\sigma = 0$","$\sigma = 1000$"),'upper right',shadow = True)

"""
Explanation: Nun etwas genauer!
End of explanation
"""


def two_grid_it_matrix_periodic(n,nc, sigma, nu1=3,nu2=3,typ='wjac'):
    cg = coarse_grid_correction_periodic(n,nc,sigma)
    if typ is 'wjac':
        smoother = iteration_matrix_wjac(n,sigma, periodic=True)
    if typ is 'gs':
        smoother = iteration_matrix_gs(n,sigma, periodic=True)
    
    pre_sm = matrix_power(smoother, nu1)
    post_sm = matrix_power(smoother, nu2)

    return pre_sm.dot(cg.dot(post_sm))    

"""
Explanation: Frage 
Was sieht man in den beiden vorhergehenden Plots?
Als nächstes behandeln wir den periodischen Fall.
End of explanation
"""


plot_fourier_transformed(two_grid_it_matrix_periodic(16,8,-100,typ='wjac'))

"""
Explanation: Buntebilderaufgabe:
Nutzen Sie plot_fourier_transformed um für $n=16$, $n_c=8$ und verschiedene $\sigma\in[-1000,1000]$ um die Zweigittermatrizen und deren Fourier-transformierten zu plotten.
End of explanation
"""


n_range = np.arange(10,100)
hs_sysmat_m1000 = map(lambda n: hs_norm(to_dense(system_matrix_hh1d(n,-1000))-to_dense(system_matrix_hh1d_periodic(n,-1000))),n_range)
hs_sysmat_0 = map(lambda n: hs_norm(to_dense(system_matrix_hh1d(n,0.001))-to_dense(system_matrix_hh1d_periodic(n,0.001))),n_range)
hs_sysmat_p1000 = map(lambda n: hs_norm(to_dense(system_matrix_hh1d(n,1000))-to_dense(system_matrix_hh1d_periodic(n,1000))),n_range)


plt.plot(hs_sysmat_m1000)
plt.plot(hs_sysmat_0)
plt.plot(hs_sysmat_p1000)

"""
Explanation: Frage
Was fällt auf? (Das sind die besten Fragen... zumindest für den Übungsleiter.)
Aufgabe: 
Nutzen Sie die Funktion two_grid_it_matrix_periodic für den periodischen Fall und plotten Sie den Spektralradius über $\sigma$ und den Spektralradius über $n$ für 3 verschiedene $\sigma$.
Aufgabe: Plotten sie die Differenzen zwischen dem periodischem und nicht-periodischem Fall.
Bonusbuntebilderaufgabe
Vergleicht analog zu den Eigenwertplots der Systemmatrizen, die Eigenwertplots der Zweigitteriterationsmatrizen.
Asymptotische Äquivalenz  zwischen periodisch und nicht-periodisch
Wir sehen, dass die Spektralradiien auf den ersten Blick gut übereinstimmen. Wir wollen nun empirisch ergründen ob die Matrizenklassen der periodischen und nicht periodischen Fällen möglicherweise zueinander asymptotisch äquivalent sind.
Zur Erinnerung:
Hilbert-Schmidt Norm
Wir definieren die Hilbert-Schmidt Norm einer Matrx $A \in K^{n \times n}$ als
$$ |A| = \left( \frac{1}{n}\sum_{i = 0}^{n-1}\sum_{i = 0}^{n-1} |a_{i,j}|^2 \right)^{1/2}.$$
Es gilt
1. $|A| = \left( \frac{1}{n}\mbox{Spur}(A^A) \right)^{1/2}$
1. $|A| = \left( \frac{1}{n}\sum_{k=0}^{n-1}\lambda_k\right)^{1/2}$, wobei $\lambda_k$ die Eigenwerte von $A^A$ sind
1. $|A| \leq \|A\|$
Asymptotisch äquivalente Folgen von Matrizen
Seien ${A_n}$ und ${B_n}$ Folgen von $n\times n$ Matrizen, welche 
beschränkt bzgl. der starken Norm sind:
$$ \|A_n\|,\|B_n\| \leq M \le \infty, n=1,2,\ldots $$ 
und bzgl. der schwachen Norm konvergieren
$$\lim_{n \to \infty} |A_n -B_n| = 0.$$
Wir nennen diese Folgen asymptotisch äquivalent und notieren dies als $A_n \sim B_n$. 
Für ${A_n}$ , ${B_n}$ und ${C_n}$, 
welche jeweils die Eigenwerte ${\alpha_{n,i}}$,${\beta_{n,i}}$ und ${\zeta_{n,i}}$ haben gelten
folgende Zusammenhänge.

Wenn $A_n \sim B_n$, dann $\lim_{n \to \infty} |A_n| = \lim_{n \to \infty} |B_n| $
Wenn $A_n \sim B_n$ und $B_n \sim C_n$, dann  $A_n \sim C_n$
Wenn $A_nB_n \sim C_n$ und $\|A_n^{-1}\|\leq K \le \infty$, dann gilt $B_n \sim A_n^{-1}C_n$
Wenn $A_n \sim B_n$, dann $\exists -\infty \le m,M\le \infty$, s.d. $m\leq \alpha_{n,i}, \beta_{n,i}\leq M \; \forall n\geq 1 \mbox{und}\; k\geq 0$  
Wenn $A_n \sim B_n$, dann gilt $ \lim_{n \to \infty} \frac{1}{n} \sum_{k=0}^{n-1} (\alpha_{n,k}^s - \beta_{n,k}^s) = 0$

Aufgabe: 
Schreiben sie eine Funktion hs_norm, welche die Hilbert-Schmidt Norm berechnet in maximal 3 Zeilen. 
Tipp: 'numpy norm frobenius' -> Suchmaschine
Aufgabe: 
Überprüfen Sie empirisch ob die

Systemmatrizenklassen 
Glättungsiterationsmatrizenklassen 
Grobgitterkorrekturmatrizenklassen 
Zweigitteriterationsmatrizenklassen 

asymptotisch äquivalent sind für $\sigma = { -1000, 0.001, 1000 }$.
Systemmatrizen:
End of explanation
"""


n_range = 2**np.arange(1,11)
hs_wjac_m1000 = map(lambda n: hs_norm(to_dense(iteration_matrix_wjac(n,-1000))-to_dense(iteration_matrix_wjac(n,-1000,False))),n_range)
hs_wjac_0 = map(lambda n: hs_norm(to_dense(iteration_matrix_wjac(n,0))-to_dense(iteration_matrix_wjac(n,0,False))),n_range)
hs_wjac_p1000 = map(lambda n: hs_norm(to_dense(iteration_matrix_wjac(n,1000))-to_dense(iteration_matrix_wjac(n,1000,False))),n_range)

plt.plot(n_range, hs_wjac_m1000)
plt.plot(n_range, hs_wjac_0)
plt.plot(n_range, hs_wjac_p1000)

"""
Explanation: Glättung:
Jacobi
End of explanation
"""


n_range = 2**np.arange(1,11)
hs_gs_m1000 = map(lambda n: hs_norm(to_dense(iteration_matrix_gs(n,-1000))-to_dense(iteration_matrix_gs(n,-1000,False))),n_range)
hs_gs_0 = map(lambda n: hs_norm(to_dense(iteration_matrix_gs(n,0))-to_dense(iteration_matrix_gs(n,0,False))),n_range)
hs_gs_p1000 = map(lambda n: hs_norm(to_dense(iteration_matrix_gs(n,1000))-to_dense(iteration_matrix_gs(n,1000,False))),n_range)

plt.plot(n_range, hs_gs_m1000)
plt.plot(n_range, hs_gs_0)
plt.plot(n_range, hs_gs_p1000)

"""
Explanation: Gauss-Seidel
End of explanation
"""


def einmal_einpacken(A):
    return np.r_[[np.zeros(A.shape[0]+1)],np.c_[np.zeros(A.shape[0]),A]]

"""
Explanation: Grobgitterkorrektur
Hier trifft man mal wieder auf das Problem, dass die Freiheitsgrade im periodischen und nicht periodischen Fall unterschiedlich verteilt sind.
End of explanation
"""


n_f_range = 2**np.arange(3,10)
n_c_range = 2**np.arange(2,9)
hs_cgc_m1000 =  map(lambda nf,nc: hs_norm(einmal_einpacken(coarse_grid_correction(nf-1,nc-1,-1000))-coarse_grid_correction_periodic(nf,nc,-1000)),n_f_range ,n_c_range)
hs_cgc_0 =  map(lambda nf,nc: hs_norm(einmal_einpacken(coarse_grid_correction(nf-1,nc-1,0))-coarse_grid_correction_periodic(nf,nc,0.001)),n_f_range ,n_c_range)
hs_cgc_p1000 =  map(lambda nf,nc: hs_norm(einmal_einpacken(coarse_grid_correction(nf-1,nc-1,1000))-coarse_grid_correction_periodic(nf,nc,1000)),n_f_range ,n_c_range)

plt.semilogy(n_f_range, hs_cgc_m1000)
plt.semilogy(n_f_range, hs_cgc_0)
plt.semilogy(n_f_range, hs_cgc_p1000)
# plt.semilogy(n_f_range, 1/np.sqrt(n_f_range))

"""
Explanation: Frage
Welcher Trick wird hier verwendet um mit den unterschiedlichen Dimensionen der Matrizen umzugehen?
End of explanation
"""


n_f_range = 2**np.arange(3,12)
n_c_range = 2**np.arange(2,11)
hs_2grid_m1000 =  map(lambda nf,nc: hs_norm(
        einmal_einpacken(two_grid_it_matrix(nf-1,nc-1,-1000))-two_grid_it_matrix_periodic(nf,nc,-1000))
                    ,n_f_range ,n_c_range)
hs_2grid_0 =  map(lambda nf,nc: hs_norm(
        einmal_einpacken(two_grid_it_matrix(nf-1,nc-1,0.001))-two_grid_it_matrix_periodic(nf,nc,0.001))
                ,n_f_range ,n_c_range)
hs_2grid_p1000 =  map(lambda nf,nc: hs_norm(
        einmal_einpacken(two_grid_it_matrix(nf-1,nc-1,1000))-two_grid_it_matrix_periodic(nf,nc,1000))
                    ,n_f_range ,n_c_range)

plt.semilogy(n_f_range, hs_2grid_m1000)
plt.semilogy(n_f_range, hs_2grid_0)
plt.semilogy(n_f_range, hs_2grid_p1000)
plt.semilogy(n_f_range, 1/np.sqrt(n_f_range)*30)

"""
Explanation: Zweigitter
End of explanation
"""


