pinvred

来源：互联网发布：c语言web攻击脚本编辑：程序博客网时间：2024/06/06 22:45

#!/usr/bin/env python

import sys

import typedbytes

import numpy as np

import numpy.dual as nd

#x:2,3,5,10,12,15,16,17,18,26,27,28,32,35,36,38,39,40,42,43,44,45,46,47,48,49,51,52,53,55,57,58,59,61,63

#y:1,13,19,29,30,31,33,50

#xshape:8,4

#y_nums:2

#load cfg from test.cfg file

def load_config(d): #该函数以(key,value)形式读取配置文件。

f = open("test.cfg")

for line in f.readlines():

if line[0] == '#': continue

(key,val) = line.strip().split(':')

d[key] = [int(i) for i in val.split(',') if i.isdigit()]

if len(d[key]) == 0: d[key] = val.split(',')

def getA(X, Y, T, Ttype='percent' ): #Ttype='percent'表示Ttype的默认是percent。如果传进来的参数是其他的，将会被覆盖掉。

Adict = dict()

U,S,V = nd.svd(np.dot(X,np.conj(X.T))) #numpy.dot(a,b,out=None) Dot product of two arrays.两个数组的点积，numpy.conj()按元素返回

full_E = sum(S) #对矩阵S的各个元素求和 #共轭复数。X.T表示X的转置矩阵，所以np.conj(X.T)表示X的共轭转置 .T返回自身的转置 .H返回自身的共轭转置(师兄说这儿写成U,S,V = nd.svd(np.dot(X,X.H))会出错 .I返回自身的逆矩阵。

Z = np.dot(np.conj(U.T),X) #看论文P34

#for y in Y:#y is single channel signal

# Z_y_projection = np.dot(np.conj(Z),y.T)

# B = Z_y_projection/S

#Blist = [np.dot(np.conj(Z),y.T)/S for y in Y]

B = np.dot(np.conj(Z),Y.T).T/S #看P35的关于B的理论推导

for t in T:

A_list = list()

#print >>sys.stderr,U.shape,X.shape

#truncate

if Ttype == 'percent': #Ttype是指奇异值截断门限的类型，有两种类型可选，一种是num类型，此时T给出的值代表保留奇异值的个数，一种是percent

T_E = full_E * t #类型，此时T给出的值是指奇异值保留的能量百分比。

sumE = 0

for i in xrange(len(S)): #len(S)可以对矩阵用？？师兄说SVD分解里边的S矩阵都被抽象成了一个由对角线元素组成的向量

sumE += S[i]

if sumE > T_E: break

elif Ttype == 'num':

i = t-1 #这是干什么用的，是为了和40行的i保持一致，进而可以在54行进行统一处理。

print >>sys.stderr,'t:',t,'retain num:',i+1,'len(S):',len(S),'S',S

#for b in Blist:

# #PCA by Energy of S

# b[i+1:] = 0

# a = b.dot(np.conj(U.T))

# A_list.append(a)

#PCA by Energy of S

tB = B.copy() #复制B矩阵

tB[:,i+1:] = 0 #矩阵的第i+1列到最后一列为0，即进行奇异值截断

A = tB.dot(np.conj(U.T)) #论文p35

Adict[t] = A #every line in A is for every ch in Y (every line in Y)

return Adict

if __name__ == "__main__": #如果该包作为入口函数，那么这句话下边的程序将会得到执行，如果其他包为入口函数，那么下边的程序将不会被执行。

cfgdict = dict()

load_config(cfgdict)

x_ch_list = cfgdict['x']

y_ch_list = cfgdict['y']

train_cond_list = cfgdict['train_cond']

t_type = cfgdict['Ttype'][0]

T = cfgdict['T'] #T应该是一个列表。看看师兄论文的P44对于配置文件test.cfg的解释和P59对于截断奇异值的说明

if t_type == 'percent': T = [t/100.0 for t in cfgdict['T']]

Xshape = (len(x_ch_list),len(train_cond_list)) #5x7的矩阵 Xshape=(5,7)

Yshape = (len(y_ch_list),len(train_cond_list)) #8x7的矩阵 Yshape=（8，7）

input = typedbytes.PairedInput(sys.stdin)

output = typedbytes.PairedOutput(sys.stdout)

X = np.zeros(Xshape,dtype=np.complex) #X = 5x7的零矩阵

Y = np.zeros(Yshape,dtype=np.complex) #y is y_num * shape(1) matrix; y_num is number of y chs

#key:hz value:(mode,channel,fft_val)=(n,m,fft_value)

last_key = None

for (key, value) in input:

if last_key is not None and key != last_key:

Adict = getA(X,Y,T,Ttype=t_type) #每一个频点(key)都会对应T列表所对应的五个传递矩阵A

output.write((last_key,Adict)) #输出的key还是hz,value是T所对应的五种情况下求得的传递矩阵A

X = np.zeros(Xshape,dtype=np.complex)

Y = np.zeros(Yshape,dtype=np.complex) #y is y_num * shape(1) matrix; y_num is number of y chs

#jugg which this channel belongs to and index. etc. x[i] or y[i]

last_key = key

ch = int(value[1]) #ch=m

mod_i = train_cond_list.index(int(value[0]))#value:(mode,channel,fft_val)=(n,m,fft_value)，index为找到int(value[0])在train_cond_list中是第几个

#len(train_cond_list)=Xshape[1]=Yshape[1]

if mod_i >= Xshape[1]: continue #Xshape=(5,7)，所以Xshape[1]=7，又mod_i是从0开始的，所以当>=的时候就不满足了，得继续往下。

if ch in x_ch_list:

index_x = x_ch_list.index(ch)

X[index_x,mod_i] = value[2] #把某指定频点(key)下的fft赋给X矩阵中的对应元素。

elif ch in y_ch_list:

index_y = y_ch_list.index(ch)

Y[index_y,mod_i] = value[2] # encounter and record y[i] of this Hz 把某指定频点(key)下的fft赋给X矩阵中的对应元素。

else:

index_x,index_y = None,None

raise Exception('Unexcepted ch!')

Adict = getA(X,Y,T,Ttype=t_type) #这儿应该是为了求最后一个key所对应的传递矩阵A，因为上边的for循环不能求最后一个key所对应的A。

output.write((last_key,Adict))

0 0