『科学计算』线性代数部分作业

from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D


A=np.array([[1],[2],[3]])
B=np.array([[1],[1],[1]])

x=np.linspace(-0.5,1,10)

x.shape=(1,10)
xx=A.dot(x)

C=A.T.dot(B)
AA=np.linalg.inv(A.T.dot(A))

P=A.dot(AA).dot(C)
E=B-P


fig = plt.figure() 
ax = fig.add_subplot(111,projection=‘3d‘)

ax.plot(xx[0,:],xx[1,:],xx[2,:],label="lineA")
ax.plot(A[0],A[1],A[2],‘ko‘,label="A")
ax.plot([0,B[0]],[0,B[1]],[0,B[2]],‘m-o‘,label="0B")
ax.plot([B[0][0],P[0][0]],[B[1][0],P[1][0]],[B[2][0],P[2][0]],‘r-o‘,label="BP")
ax.plot([0,E[0]],[0,E[1]],[0,E[2]],‘y-o‘,label="0E")

ax.legend()
ax.axis(‘equal‘)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import copy

# 产生一个方波(x,y)
x = np.linspace(-10,10,300) 
y=[]
for i in np.cos(x):
    if i>0:
        y.append(0)
    else:
        y.append(2)
y=np.array(y)

def projection(A,b):
    ####
    # return A*inv(AT*A)*AT*b
    ####
    AA = A.T.dot(A)
    w=np.linalg.inv(AA).dot(A.T).dot(b)
    print(w)
    return A.dot(w)

def fourier(x,y,N):
    A = []
    for i in x:
        A.append(copy.deepcopy([]))
        for j in range(N):
            A[-1].append(np.sin((j + 1) * i))
            A[-1].append(np.cos((j + 1) * i))
        A[-1].append(1)
    b = y.T
    yw = projection(np.array(A),b)
    return yw

# write Your code, Fourier function    
plt.plot(x,y,color=‘g‘,label=‘origin‘) 
plt.plot(x,fourier(x,y,3),color=‘r‘,label=‘3‘)
plt.plot(x,fourier(x,y,8),color=‘b‘,label=‘8‘)
plt.plot(x,fourier(x,y,23),color=‘k‘,label=‘23‘)

plt.legend()
plt.axis(‘equal‘)
plt.show()

import numpy as np
import copy
import matplotlib.pyplot as plt
from matplotlib import animation

A=np.array([[3,1],[2,4]])/4.0

# write Your code
fig = plt.figure()
ax = fig.add_subplot(111)
line1, = ax.plot([],[],c=‘r‘)
line2, = ax.plot([],[],c=‘b‘)
ax.axis(‘equal‘)
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)
ax.set_xticks(np.linspace(-2,2,5))
ax.set_yticks(np.linspace(-2,2,5))
ax.grid(True)

def data():
    point_num = 100
    x = np.cos(np.linspace(0,2 * np.pi,point_num))
    y = np.sin(np.linspace(0,2 * np.pi,point_num))
    rot_x,rot_y = copy.deepcopy([]),copy.deepcopy([])
    int_x,int_y = copy.deepcopy([]),copy.deepcopy([])
    for i in range(point_num):
        rot_x.append(A.dot([x[i],y[i]])[0])
        rot_y.append(A.dot([x[i],y[i]])[1])
        int_x.append(x[i])
        int_y.append(y[i])
        yield (rot_x,rot_y,int_x,int_y)

def update(data):
    line1.set_xdata(data[0])
    line1.set_ydata(data[1])
    line2.set_xdata(data[2])
    line2.set_ydata(data[3])
    return line1,line2

def init():
    line1.set_data([],[])
    line2.set_data([],[])
    return line1,line2

ani = animation.FuncAnimation(fig,update,data,init_func=init,interval=100,repeat=False)
plt.show()

import numpy as np
import matplotlib.pyplot as plt  


import time
def time_cost(f):
    def _f(*arg, **kwarg):
        start = time.clock()
        a=f(*arg,**kwarg)
        end = time.clock()
        print(f.__name__,"run cost time is ",end-start)
        return a
    return _f

   
@time_cost
def fib_opt_seq(seq):
    return [fib_opt(i) for i in seq]
    
def fib_opt(n):
    a,b,i=0,1,0
    
    while i<n:
        a,b=b,a+b
        i+=1
    else:
        #print(b)
        return b    

        
import random
#seq = [random.randint(800,1000) for i in xrange(1000)]
seq = range(1000)
a=fib_opt_seq(seq)

# write Your code fib_eig_seq function
A = np.array([[1,1],[1,0]])
eigval,eigvect = np.linalg.eig(A)
S_inv = np.linalg.inv(eigvect)

@time_cost
def fib_eig_seq(seq):
    return [fib_eig(i) for i in seq]


def fib_eig(n):
    af = (np.diag(eigval)**(n+1)).dot(S_inv)
    #print((eigvect.dot(af)).dot(np.array([[1],[0]]))[1])
    return (eigvect.dot(af)).dot(np.array([[1],[0]]))[1]

b=fib_eig_seq(seq)

import numpy as np

A = np.array([[ 0.8 ,  0.1],
              [ 0.2 ,  0.9]])

N_year = 1000

x=np.array([3200,
            4000])

def hw_3_5(a,x,n):
    eigval, eigvact = np.linalg.eig(a)
    res = (eigvact.dot((np.diag(eigval)**n).dot(np.linalg.inv(eigvact)).dot(x.T)))
    print(res)
    return res

hw_3_5(A,x,N_year)

import os
import numpy as np
import matplotlib.pyplot as plt  
from tensorflow.examples.tutorials.mnist import input_data

print(os.getcwd())
os.environ[‘TF_CPP_MIN_LOG_LEVEL‘] = ‘3‘

batch_size = 10
mnist = input_data.read_data_sets(‘../../Mnist_data/‘, one_hot=True)
batch_xs,batch_ys = mnist.train.next_batch(batch_size)
print(batch_xs.shape)

# write Your code
def pca(data,k=0.5):
    # data 原始的图片
    # k是保留特征值的百分比
    # return 返回降维后重建的图片
    res = np.empty_like(data).reshape(batch_size, 28, 28)
    rec = np.empty_like(res)
    data = data.reshape(batch_size, 28, 28)
    for i in range(batch_size):
        cov = np.cov(data[i])
        eigVal, eigVact = np.linalg.eig(cov)
        print(cov.shape, eigVal.shape, eigVact.shape)
        for j in range(len(eigVal) - 1):
            # print(sum(eigVal[:j]),sum(eigVal[:]),sum(eigVal[:j])/sum(eigVal[:]))
            if sum(eigVal[:j])/sum(eigVal[:]) > k:
                print(‘提取前面‘,j, ‘个特征‘)
                break
        # res[i] = eigVact[:,:j].T.dot(data[i].T)# np.dot(eigVact[:j],data[i].T)
        a = eigVact#[:,:j]
        b = data[i]
        print(a.shape,b.shape)
        res[i] = np.dot(b, a)
        rec[i] = np.dot(res[i],a.T)
    f,a = plt.subplots(2,batch_size,figsize=(10,2))
    for i in range(batch_size):
        a[0][i].imshow(data[i],‘gray‘)
        a[0][i].set_xticks([])
        a[0][i].set_yticks([])
        a[1][i].imshow(rec[i].reshape(28,28),‘gray‘)
        a[1][i].set_xticks([])
        a[1][i].set_yticks([])



# N = 20
recon_data = pca(batch_xs)
# show_pic(data[:N,:],recon_data[:N,:])