# +----------------------------------------------------+
# |                  Lekce 9 - demo 2                  |
# +----------------------------------------------------+
#  Linearni algebra - numpy
#   vlastni implementace LU rozkladu s pivotingem
#   a jeji testik

import numpy as np
import matplotlib.pyplot as plt

#S PIVOTACI, ale bez ulozeni P[], implementace LU rozkladu pro testy
#  matice a je znicena a nahrazena lu faktory
#  matice musi byt ctvercova n x n (netestuje se)
def LUfactor(a):
  n=a.shape[0]
  for k in range(n-1):
      i=k+np.argmax(abs(a[k:,k]))
      v=a.copy()[i,:]
      a[i,:]= a[k,:]
      a[k,:]= v
      for j in range(k+1,n):
          a[j,k]=a[j,k]/a[k,k]
          a[j,k+1:n]=a[j,k+1:n]-a[j,k]*a[k,k+1:n]

#Testik na matici jejiz rozklad zname:
Lfactor=np.array([[1.0,0.0,0.0,0.0],[2.0,1.0,0.0,0.0],[3.0,2.0,1.0,0.0],[4.0,3.0,2.0,1.0]])
Ufactor=np.array([[1.0,2.0,1.0,2.0],[0.0,1.0,1.0,2.0],[0.0,0.0,1.0,3.0],[0.0,0.0,0.0,1.0]])

aa=np.dot(Lfactor,Ufactor)
a=aa.copy()
LUfactor(a)

#Testik podminenosti (v norme)
d=np.dot(np.tril(a,-1)+np.eye(len(a)),np.triu(a))

#Statistika
Nstat=100
Nmax=50
an=[]
ak=[]
akl=[]
aku=[]
for n in range(5,Nmax+1):
  for i in range(Nstat):
    an.append(n)
    a=np.random.rand(n,n)
    ak.append(np.linalg.cond(a))
    LUfactor(a)
    akl.append(np.linalg.cond(np.tril(a,-1)+np.eye(len(a))))
    aku.append(np.linalg.cond(np.triu(a)))

#Rust prumerneho kappa(a)
#plt.plot(np.log10(an),np.log10(ak),'ro')
#plt.plot(np.log10(an),np.log10(np.array(an)**2),'b')
#plt.show()

#Rust prumerneho kappa(l)
#plt.plot(np.log10(an),np.log10(akl),'ro')
#plt.plot(np.log10(an),np.log10(np.array(an)**2),'b')
#plt.plot(np.log10(an),np.log10(np.array(an)**1.5),'g')
#plt.show()

#Rust prumerneho kappa(u)
#plt.plot(np.log10(an),np.log10(aku),'ro')
#plt.plot(np.log10(an),np.log10(np.array(an)**2),'b')
#plt.plot(np.log10(an),np.log10(np.array(an)**2.5),'g')
#plt.show()

#Zhorseni cisla podminenosti:
plt.plot(np.log10(an),np.log10(np.array(aku)*np.array(akl)),'ro')
plt.plot(np.log10(an),np.log10(ak),'bo')
plt.plot(np.log10(an),np.log10(np.array(an)**2),'b')
plt.plot(np.log10(an),np.log10(np.array(an)**2.5),'g')
plt.show()