pinvred
2015-06-28 11:39
330 查看
#!/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))
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))
相关文章推荐
- MongoDB基本语法
- 分页函数优化注意点(一):关注业务数据
- Yii查询生成器(Query Builder)用法实例教程
- pinvmap
- 关于linux中交叉编译器的配置
- 项目报错:Cannot find class file for javax/servlet/ServletException
- Java读书笔记一(异常处理)
- LTE中的数据传输(2)——下行数据传输
- fft
- UVA 10791
- 移动端 meta 标签笔记
- 用JavaScript做浏览器对象事件的兼容性
- 黑马程序员——Java基础——内部类
- 【HDU 4609】3-idiots(FFT)
- 解决Tomcat中设置默认项目时只能访问静态页面的问题
- C/C++ sort函数的用法
- 2015 股市小牛 记载
- Win32SDK编辑框控件的简单操作
- js中bind、call、apply函数的用法
- C#基于SQLiteHelper类似SqlHelper类实现存取Sqlite数据库的方法