Untitled

import numpy as np
from numpy import array , eye , outer
import scipy
from scipy.linalg import norm
import math
from matplotlib.image import imread
from matplotlib import pyplot as plt



def proj(v,U):
    vp = np.zeros(np.size(v))
    for u in U:
            vp += (np.dot(u,v)/np.dot(u,u))*u
    return vp

def gramschmidt(A):
    n,m = np.shape(A)
    u1 = A[:,0]
    e1 = u1/math.sqrt(np.dot(u1,u1))
    U = [e1]
        # A[:,1:] har alla kolumner förutom första
    for v in A[:,1:].T:
        u = v - proj(v,U)
        e = u/math.sqrt(np.dot(u,u))
        U.append(e)

    Q = np.zeros(np.shape(A))
    for i,e in enumerate(U):
        Q[:,i] = e
    return Q

def tests(Q):
    m,n = np.shape(Q)
    print(f'2Norm of (Q^TQ) = {norm(Q,ord=2)}')
    print(f'2Norm of (Q^TQ) - I = {norm(Q.T@Q-np.identity(n),ord=2)}')
    # eigenvärden också av Q^TQ
    print(f'Eigenvalues of (Q^TQ) = {np.linalg.eigvals(Q.T@Q)}')

def Householder_step(A):
    m,n = np.shape(A)
    a = A[:,0]
    ahat=array([norm(a)]+(m-1)*[0.])
    v = a - ahat
    v=v/norm(v)
    H = eye(m)-2*outer(v,v)
    return H

def Givens_step(A):
    m,n = np.shape(A)
    a = A[:,0]
    # kanske behöver börja nerifrån om fel
    a0 = a[0]
    G = np.eye(m)
    for i in range(m-1,0,-1):
        ahat = a[i]
        theta = math.atan(ahat/a0)
        Gi = np.block([[math.cos(theta), np.zeros( (1,i-1) ), -math.sin(theta),  np.zeros(   ( 1, m-(i+1) ) )],
        [np.zeros((i-1,1)), np.eye(i-1), np.zeros((i-1,m-i))   ],
        [math.sin(theta),  np.zeros((1,i-1)), math.cos(theta),  np.zeros( (1,m-(i+1)) )],
        [np.zeros((m-(i+1),i+1)), np.eye(m-(i+1))] ])
        G = Gi@G

    return G


def QRH(A):
    m,n = np.shape(A)
    tempA = A
    H = np.eye(m)
    for i in range(n):
        # byt ut denna för reflection vs rotation
        Hi = Householder_step(tempA)
        tempA = Hi@tempA
        # ta bort första rad första kolumn
        tempA = np.delete(tempA,0,0)
        tempA = np.delete(tempA,0,1)
        H = H@np.block([[np.eye(i),  np.zeros((i,m-i))],[np.zeros((m-i,i)), Hi]])
    R = H.T@A
    Q = H
    return Q,R

def QRG(A):
    m,n = np.shape(A)
    tempA = A
    H = np.eye(m)
    for i in range(n):
        # byt ut denna för reflection vs rotation
        Hi = Givens_step(tempA)
        tempA = Hi@tempA
        # ta bort första rad första kolumn
        tempA = np.delete(tempA,0,0)
        tempA = np.delete(tempA,0,1)
        H = H@np.block([[np.eye(i),  np.zeros((i,m-i))],[np.zeros((m-i,i)), Hi]])
    R = H.T@A
    Q = H
    return Q,R

n = 4
m = n + 2
A = np.random.rand(m,n)
Q1,_ = scipy.linalg.qr(A)
Q2,R2 = QRG(A)
#print(f'Q1Q1^T = {Q1@Q1.T}')
print(f'Det (Q1) = {np.linalg.det(Q1)}')
#print(f'Q2Q2^T = {Q2@Q2.T}')
print(f'Det (Q2) = {np.linalg.det(Q2)}')

print(f'A = {A}')
print(f'Q2R2 = {Q2@R2}')
R2[R2<10**(-10)] = 0
print(f'R2 = {R2}')


# B = imread('kvinnagr.jpeg')
# print(B)
# U,s,VT = scipy.linalg.svd(B)
# m,n = np.shape(B)
# S = scipy.linalg.diagsvd(s,m,n)

#S = np.zeros((m,n))
#for i in range(0,min(m,n)):
#    S[i,i] = s[i]
#
# b = np.round(U@S@VT)
# b = b.astype(int)
#
# print(f'Difference of B and USVT = {norm(B-b,ord=2)}')
# plt.imshow(b,cmap='gray', vmin=0, vmax=255.)

#Q2 = gramschmidt(A)
#tests(Q1)
#tests(Q2)