数学专题-算法-矩阵-拟合

'''
@Author: Mikeliu@2020
@Date: 2020-03-08 10:32:13
@LastEditTime: 2020-03-12 09:26:47
@LastEditors: Mikeliu@2020
@Description: 
'''
#!/usr/bin/python
# -*- coding: utf-8 -*- 
import math
import datetime
'''
' Excel 表格操作'
'''
import pandas as pd
'''
' 多项式拟合画图'
'''
import matplotlib.pyplot  as plt
'''
' 斯密特正交化操作，矩阵相关'
'''
import numpy as np
from numpy import *
from scipy.linalg import *
'''
计算阶乘使用的库
'''
from functools import reduce   
from operator import mul
'''
计算圆周率精度使用的库
'''
from decimal import Decimal
from decimal import getcontext

# pip install -i https://pypi.tuna.tsinghua.edu.cn/simple --trusted-host pypi.tuna.tsinghua.edu.cn  模块名称
## 打印系统时间 
def print_LocalTime():
    μlocaltime = datetime.datetime.now()
    print("系统当前时间(YY-MM-DD Hour:min:sec.μs) %s",μlocaltime)
    return 
# 求向量的三种范式
A = np.array([[1,3],[-1,2]])
print ('\nA_array:\n',A)
ret = np.linalg.norm(A,ord=None) # 对象是矩阵，求矩阵所有元素平方和开方
print('所有元素平方和开方',ret)
# In[34]
ret = np.linalg.norm(A,ord=None,axis=1) # 对象是逐行，求每行元素平方和开方
print('每行元素平方和开方',ret)
# In[35]
ret = np.linalg.norm(A,ord=None,axis=0)  # 对象是逐列，求每列元素平方和开方
print('列元素平方和开方',ret)
# In[36]
# 对象是矩阵，求矩阵中逐行元素绝对值和的最大值
ret = np.linalg.norm(A,ord=np.Inf)
print('逐行元素绝对值和的最大值',ret)
# In[37]
ret = np.linalg.norm(A,ord=np.Inf,axis=1) # 对象是逐行，求每行元素绝对值的最大值
print('每行元素绝对值的最大值',ret)
# In[38]
ret = np.linalg.norm(A,ord=np.Inf,axis=0) # 对象是逐列，求每列元素绝对值的最大值
print('每列元素绝对值的最大值',ret)
# In[39]
# 对象是矩阵，求矩阵中逐行元素绝对值和的最小值
ret = np.linalg.norm(A,ord=-np.Inf)
print('逐行元素绝对值和的最小值',ret)

ret = np.linalg.norm(A,ord=-np.Inf,axis=1) # 对象是逐行，求行的元素绝对值的最小值
print('行的元素绝对值的最小值',ret)

ret = np.linalg.norm(A,ord=-np.Inf,axis=0) # 对象是逐列，求列的元素绝对值的最小值
print('列的元素绝对值的最小值',ret)
# In[40]
A = np.array([[3,-4,1],[-6,5,0]])
print ('\nArray:\n',A)
ret = np.linalg.norm(A,ord=0,axis=1) # 对象是逐行，求每行的非零元素个数
print('每行的非零元素个数',ret)
# In[41]
ret = np.linalg.norm(A,ord=0,axis=0) # 对象是逐列，求每列的非零元素个数
print('每列的非零元素个数',ret)
# In[42]
A = np.array([[100,4],[-6,5]])
print ('\nArray:\n',A)
ret = np.linalg.norm(A,ord=1)  # 对象是矩阵，求列的元素绝对值和的最大值
print('列的元素绝对值和的最大值',ret)

ret = np.linalg.norm(A,ord=1,axis=1)  # 对象是逐行，求行的元素绝对值和
print('行的元素绝对值和',ret)

ret = np.linalg.norm(A,ord=1,axis=0)  # 对象是逐列，求；列的元素绝对值和
print('列的元素绝对值和',ret)
# In[44]
A = np.array([[-4, 3, 4], [-6,4, 3]])
print ('\nArray:\n',A)
#10.0015300174, 对象是矩阵，求所有元素平方和开方，所有元素平方和102，开方=10.0015300174
ret = np.linalg.norm(A, ord=2)
print('元素平方和开方',ret)
#[ 6.40312424  7.81024968]，对象是每行，求行的元素平方和开方
ret = np.linalg.norm(A, ord=2, axis=1)
print('行的元素平方和开方',ret)
#[ 7.21110255  5. 5.]，对象是每列，求列的元素平方和开方
ret = np.linalg.norm(A, ord=2, axis=0)
print('列的元素平方和开方',ret)

x = np.array([
    [0, 3, 4],
    [1, 6, 4]])

print('\nArray:\n',x)
#默认参数ord=None，axis=None，keepdims=False
print ("默认参数(矩阵整体元素平方和开根号，不保留矩阵二维特性)：\n",np.linalg.norm(x))
print ("矩阵整体元素平方和开根号，保留矩阵二维特性：\n",np.linalg.norm(x,keepdims=True))
 
print ("矩阵每个行向量求向量的2范数：\n",np.linalg.norm(x,axis=1,keepdims=True))
print ("矩阵每个列向量求向量的2范数：\n",np.linalg.norm(x,axis=0,keepdims=True))
 
print ("矩阵1范数：\n",np.linalg.norm(x,ord=1,keepdims=True))
print ("矩阵2范数：\n",np.linalg.norm(x,ord=2,keepdims=True))
print ("矩阵∞范数：\n",np.linalg.norm(x,ord=np.inf,keepdims=True))
print ("矩阵-∞范数：\n",np.linalg.norm(x,ord=-np.inf,keepdims=True))
 
print ("矩阵每个行向量求向量的1范数：\n",np.linalg.norm(x,ord=1,axis=1,keepdims=True))
print ("矩阵每个列向量求向量的1范数：\n",np.linalg.norm(x,ord=1,axis=0,keepdims=True))
print('\n\n')


'''
' 这部分是研究圆周率的计算收敛效率 '
'''
## 定义阶乘方法 (n)!
def fact(n):
    if (n <= 1):
        return 1
    return reduce(mul,range(1,n+1))
## 定义双阶乘方法 (n)!!
def fact_D(n):
    if (n <= 1):
        return 1
    elif (n&1==True):
        return reduce(mul,range(1,n+1,2))    # range(start,stop[,step]) []表示可选
    else:
        return reduce(mul,range(2,n+1,2))            
## BBP方法计算pi值 每次计算的精度是每算一项得到1.25个有效数字 
## ∑_(k=0)^∞▒〖1/〖16〗^k  (4/(8^k+1)+2/(8^k+4)+1/(8^k+5)+1/(8^k+6)) 〗
def cal_pi(precision):
    getcontext().prec=precision
    A = 0
    for k in range (precision):
        A += 1/Decimal(16)**k *(Decimal(4)/(8*k+1)-Decimal(2)/(8*k+4)-Decimal(1)/(8*k+5)-Decimal(1)/(8*k+6))
    return A
## 拉姆努金计算方法，每次计算的精度是每算一项得到 8个有效数字
## 1/((2√2)/〖99〗^2  ∑_(k=0)^∞▒(4k)!(26390k+1103)/((k!)^4 (〖396〗^4k ) ))
def lmr_pi(precision):
    getcontext().prec=precision
    A=0
    for k in range (precision):
        A += (fact(4*k)*(Decimal(26390)*k+1103))/((fact(k)**Decimal(4))*(Decimal(396)**(4*k)))
    return 1/(((2*(Decimal(2).sqrt()))/(Decimal(99)**2)) *A )
## 反正切法计算圆周率，是比BBP慢，每算一项得到0.3个有效数字
def atan_pi(precision):
    getcontext().prec=precision
    A=0
    for k in range (precision):
        A += (Decimal(fact_D(2*k))/Decimal(fact_D(2*k+1)))/(Decimal(2)**k)
    return 2*A
## 莱布尼茨法计算圆周率，100项误差能精确到小数点后2位。
def LBNC_pi(precision):
    getcontext().prec=precision
    A=0
    for k in range (precision):
        A += (Decimal(-1)**k)/(Decimal(2)*k+1)
    return 4*A

print_LocalTime()
X=cal_pi(300)
print_LocalTime()
print_LocalTime()
Y=lmr_pi(300)
print_LocalTime()
print_LocalTime()
Z=atan_pi(300)
print_LocalTime()
print_LocalTime()
T=LBNC_pi(300)
print_LocalTime()
print(' X :',X)
print(' Y :',Y)
print(' Z :',Z)
print(' T :',T)
print(' Delta_X_Y %i',abs(X-Y))
print(' Delta_Y_Z %i',abs(Z-Y))
print(' Delta_Y_T %i',abs(T-Y))

'''
' 这部分是计算矩阵操作，斯密特正交化，获取矩阵的某个元素，矩阵的特征值特征向量 '
'''
File = '.\Data.xlsx'
wb = pd.read_excel(File)    # 读取Excel表格 
print ('wb:\n',wb)

# print("wb.loc[m,:]:\n",wb.loc[0,:]) # 打印行

for m in range (wb.shape[1]):
    # print("wb.loc[:,m]:\n",wb.loc[:,m]) # 打印 i 列信息，并求这列的和
    SUM = np.sum(wb.loc[:,m])   
    # print ('SUM:',SUM)                    # 打印 i 列总和
    for n in range (wb.shape[0]):
        wb[m].at[n] = wb[m].at[n] / SUM
        # print('wb['+str(m)+'].at['+str(n)+']:',wb[m].at[n])   # 打印 （j行，i列) 元素信息
print ('wb_unit:\n',wb)

Ret = orth(wb)              # 求Excel表格的 斯密特 单位正交化
Ret = pd.DataFrame(Ret)
print('wb_ortn:\n',Ret)

# print('\n####  验证斯密特单位正交的计算结果 ### ')
for k in range (Ret.shape[1]):
    x = Ret.loc[:,k]
    for t in range (Ret.shape[1]):
        y = Ret.loc[:,t]
        z = np.dot(x,y)                      
        print('The result ('+str(k)+','+str(t)+') is: ',z)

# 计算特征值，特征向量 方法
x = random.randint(1,100,size=(4,4))
a,b=np.linalg.eig(x) ##特征值赋值给a，对应特征向量赋值给b 
print('\nStart Other Function, Array = \n',x)

for i in range(len(a)):
    print('\n特征值:',a[i],'\n对应特征向量为:',b[:,i])

print(np.trace(x))
print(np.sum(a))

'''
' 这部分是构建三角矩阵，对称矩阵，矩阵的迹 '
'''
# 构建三角矩阵 和对称矩阵
y = random.randint(-5,5,size=(5,5))
# print('y_random\n',y)
dig = y.diagonal()
# print('y_dig\n',dig)  ## 取矩阵的主对角线元素
# print('y_dig1\n',np.diag(dig))  ## 取矩阵的主对角线元素构建方阵
# print('y_trace\n',np.trace(y))  ## 取矩阵的迹
# print('y_trace1\n',np.sum(dig))  ## 取矩阵的迹
y = np.triu(y)                    ## 生产上三角矩阵
# print('y_tri\n',y)
y += y.T -np.diag(y.diagonal())     ## 生成对称矩阵
# print('y_Double\n',y)
a,b=np.linalg.eig(y) ##特征值赋值给a，对应特征向量赋值给b 
for i in range(len(a)):
    print('\n特征值:',a[i],'\n对应特征向量为:',b[:,i])
print(np.trace(y))
print(np.sum(a))

'''
' 这部分是多项式数据拟合图形化展示 '
'''
X_Data = range(0,wb.shape[0])
for kk in range (wb.shape[1]):
    Data_source = wb.loc[:,kk] 
    Data_source = np.array(Data_source)
    #print('Data_source is :\n',Data_source)
    #print('X_Data is:',X_Data)
    f1 = np.polyfit(X_Data,Data_source,3) ## 三次多项式拟合 系数表
    #print('f1 is :\n',f1)
    p1 = np.poly1d(f1)                   ## 三次多项式拟合 方程式
    print('p1 is :\n',p1)
    yvals = p1(X_Data)                        ## 三次多项式拟合 向量数值表
    plot1 = plt.plot(X_Data, Data_source, 's',label='original values')
    plot2 = plt.plot(X_Data, yvals, 'r',label='polyfit values')
    plt.xlabel('X_Data')
    plt.ylabel('Data_source')
    plt.legend(loc=4) #指定legend的位置右下角
    plt.title('Mikeliu@2020-'+str(kk)+'th Data')
    plt.show()