Python

Python是一种解释型、面向对象、动态数据类型的高级程序设计语言。 由于其语言的简洁、易读以“胶水般”的可扩展性,越来越多的童鞋选择使用Python进行科学计算。 Numpy是Python下的一个数值计算库,提供了多维数组,线性代数、傅里叶变换等众多数值计算函数。 如果使用Python+Numpy+Scipy(提供众多科学计算函数,基于Numpy)+matplotlib(一个强大的绘图库,绘图效果秒杀matlab!),那么你几乎可以做在matlab下能做的任何事情。 网上有很多比较matlab与numpy的文章,他们的各自优缺点请自行谷歌。 在此提供一个matlab与numpy函数的对照表,从表中我们可以看出,numpy中的大多数函数的名称及使用方法与matlab中很相似,但是这并不意味这把程序从matlab移植到numpy中很容易, 其中的主要困难有:1、数组的引用方式(这是之后的重点)不同; 2、numpy中矩阵乘与数组乘的区别;3、数组的拼接方式不同,这些不同之处使得程序移植与把程序重新写一遍差不多。 matlab下有许多商业的工具箱,这些东西压根就没法移植。而如果你之前使用matlab积累了一些程序,那么把这些好不容易写成的程序都移植过去也不现实。 所以,保持两边各自的“画风”不变,仅仅交换数据,是matlab与numpy联合工作的最佳选择。 科学数据其实就是一大堆的数,其他的数据结构大都是用来辅助科学计算的,所以我们只需要实现两边多维数组数据的交换就可以了,交换其他的数据结构既费劲有没有太大的意义,好了,现在切入正题

numpy与matlab 的数组引用差别

首先我们来做一件事情,看看使用两种语言中同样功能的函数达到的不同效果:

创建一个一维数组1~120,120=2×3×4×5,然后reshape为2×3×4×5的多维数组

matlab中:

Make the height to be less than 400pixel

Matlab  


a = 1:120;
b = reshape(a, [2,3,4,5]);
b
b(:,:,1,1) =
     1     3     5
     2     4     6
b(:,:,2,1) =
     7     9    11
     8    10    12
.......
b(:,:,4,5) =
   115   117   119
   116   118   120

首先说明一下各个索引值的意义

b是一个(2,3,4,5)的四维数组,它可以这么理解

[每一列有2个元素,每一矩阵有3列,每一页有4个矩阵,一共有5页]

(其中行,列大家都懂,页就像是几个矩阵叠起来写满一页)

科学计算的时候,一定要搞清楚自己数据的数据结构:每一行,每一列都是什么物理或数学意义,一般我们习惯以列或者矩阵为一个具有物理或数学意义的单位实体。那么我这个4维数据可以这么理解。

我有一个四维数组,它存了4×5个矩阵,每一个矩阵的大小为2×3,我想引用第n页第m个矩阵,就使用b(:,:,m,n)

可以看出,让matlab自己输出这个矩阵,它就是按照这种理解输出的,或者

我所存的数据是二维坐标点,共存了3×4×5个,我可以使用b(:,m,n,k)的形式引用某一个坐标

matlab中比元素高一级的数组的单位是“列向量”,matlab储存数据是“按列”储存的。

我们拿到一个矩阵,通常使用“第几行第几列”来引用它的一个数据元素,matlab可不是这么理解的。

matlab认为:“第几行”=“第几个列向量”, “第几列”=“列向量中的第几个元素”

所以我们来看看reshape函数对a做了些神马:

reshape为(2,3,4,5),matlab先把数据分成5页,每章有24个数据,再把一页的数据分为4个矩阵,每个矩阵有6个数据,这时第一个矩阵的数据应为1~6,这时候再把每一个矩阵的数据分为3行,所以每一列数据就变成有两个元素的列向量了。

所以如果你看到一个matlab矩阵为

Make the height to be less than 400pixel

Matlab  


[1,2,3;
 4,5,6]

那么它在matlab内存中的存储顺序应该为[1,4,2,5,3,6],这一点非常的重要。

numpy中:

Make the height to be less than 400pixel

python  


a=np.array(range(1,121))
b=a.reshape((2,3,4,5))
b
array([[[[  1,   2,   3,   4,   5],
         [  6,   7,   8,   9,  10],
         [ 11,  12,  13,  14,  15],
         [ 16,  17,  18,  19,  20]],
     .........
        [[101, 102, 103, 104, 105],
         [106, 107, 108, 109, 110],
         [111, 112, 113, 114, 115],
         [116, 117, 118, 119, 120]]]])

numpy与C的数据储存方式相同,均是采用广义表的形式, 这一点上造成了numpy与C之间可以直接通过内存来传递数组(见ctypes)。

对于这里的b的索引值,我们必须这样来理解:

[一共有两章,每一章有3页,每一页有4行,每一行有5个元素]

所以在这里我们作为单位实体的矩阵应该是4×5的。为什么会这样呢?下面有请栗子君:

如果你在numpy中看到这样一个二维数组(或者你想称之为“矩阵”,但在numpy中矩阵和二维数组是截然不同的两种东西)

Make the height to be less than 400pixel

python  


np.array([[1, 2, 3],
          [4, 5, 6]])

一定不要觉得numpy和你想的一样,因为它对于电脑本身更自然一点的表达方式应该是

Make the height to be less than 400pixel

python  


np.array([[1, 2, 3], [4, 5, 6]])

哈?这有神马不同?

首先python自带一种数据结构叫“列表”,numpy只是把一个列表内全为数的列表变成了相同结构的“数组”而已,所以我们只看列表的数据结构

Make the height to be less than 400pixel

python  


[[1, 2, 3], [4, 5, 6]]

计算机是这么理解这个列表的:

恩….这个列表里有两个元素,每一个元素都是一个列表,啊…..这两个列表内分别有三个元素,每一个元素是个整数。

对于更高维度的列表有着类似的理解,广义表有种类似于“一层套一层”的结构,我们引用b[1] [2][3][4]代表的意思是

我要看b这个广义表中第二个元素中的第三个元素中的第四个元素中的第五个元素(numpy中下标0代表第一个元素),啊…这也就是这个列表中的最后一个元素了

然后我们看看reshape对a做了些什么

reshape为(2,3,4,5),首先把120个数据分为两章,每章60个数据,分为3页,每页20个数据,这个时候第一页的数据应为[1,2,3,4…20],再把每一页分为4行,每一行有5个数据。

这能说明神马呢?好像有点乱…我们还是请栗子君吧..

我有4个2×3的矩阵,我要把他们分别存在matlab和numpy中,并且想输出第三个矩阵的第二列的所对应的那个列向量,我的矩阵为

matlab格式表述[1,2,3;4,5,6],[7,8,9;10,11,12],[13,14,15;16,17,18],[19,20,21;22,23,24]

所以我其实是想让你输出[14,17]^T,恩…

matlab中:

Make the height to be less than 400pixel

Matlab  


%正确的做法
a = 1:24;
b = reshape(a,[3,2,4]);
b(2,:,3)'
%输出
14
17
%错误的做法
a = 1:24;
b = a.reshape(a,[2,3,4])
b(:,2,3)
%输出
15
16

我们其实并没有把以上4个矩阵按照那种形式存到matlab中,

b其实是4个3行2列的矩阵,不是2行3列的,我们这样做以后在其转置后才是所要求的2行3列的矩阵。为什么要这么做呢?是因为reshape的时候matlab是按列储存数据的。

numpy中:

Make the height to be less than 400pixel

python  


a=np.array(range(1,25))
b=a.reshape(4,2,3)
b[2,:,1]
#output
array([14, 17]) #这是一个(2,)的一维数组,而不是一个(2,1)的二位数组

总结起来就是:

matlab中,对于一个(m,n,a,b,c,d….)的多维数组,我们认为它是a×b×c×d…个m×n的矩阵

而在numpy中[….,c,b,a,m,n]的多维数组,我们认为他是….×c×b×a个m×n的二维数组(注意还不是矩阵)

numpy与matlab通过文件进行数据交互

前面讲了这么一大堆其实是为了给这一章做铺垫

在matlab中写两个函数getNpBin和toMatBin

getNpBin.m

Make the height to be less than 400pixel

Matlab  


function np=getNpBin()
%author:wujinnnnn@qq.com
%本函数自动读取./npbin文件夹下的*.npbin文件
%返回一个结构体np,np.XX(XX为变量名)包含了./npbin下的所有变量文件
    if not (exist('./npbin')==7)
        error('no data dir!');
    end
    np=struct();
    lsStr = ls('npbin');
    lsStr = regexprep(lsStr,'\s{2,}','+');
    fileNameCell = regexp(lsStr,'+','split');
    fileNum = length(fileNameCell);
    for ii=1:fileNum
        eachFileName = fileNameCell{ii};
        eachFileName = deblank(eachFileName);
        fprintf('loding file %s\n', eachFileName)
        parts = regexp(eachFileName,'_','split');
        if length(parts)~=4
            fprintf('file name error of %s\n',eachFileName);
            continue
        end

        varName = parts{1};
        type = parts{3};
        shapeStr = parts{2};
        shapeCell = regexp(shapeStr, '-', 'split');
        n = length(shapeCell);
        switchArray = 1:n;
        temp = switchArray(1);
        switchArray(1) = switchArray(2);
        switchArray(2) = temp;
        switchStr = '[';
        for i=1:n
            switchStr = [switchStr,num2str(switchArray(i)),','];
        end
        switchStr = [switchStr,']'];

        shapeArray = zeros(1,n);
        for i=1:n
            shapeArray(i) = str2num(shapeCell{n-i+1});
        end
        fid=fopen(['./npbin/',eachFileName],'r');
        eval(['tempRead','=fread(fid,inf,''',type,''');']);
        fclose(fid);
        tempRead = reshape(tempRead,shapeArray);
        eval(['np.',varName,'=','permute(tempRead,',switchStr,');']);
        fprintf([varName,' loaded!!\n']);
    end
end

toMatBin.m

Make the height to be less than 400pixel

Matlab  


function toMatBin(var,varStr)
%author:wujinnnnn@qq.com
%将变量var存为./matbin下的数据文件
%当它被numpy的读入后的变量名称叫varStr
   info = whos('var');
   type = info.class;
   if not (exist('./matbin')==7)
       mkdir matbin;
   end
   shape = info.size;
   type = info.class;
   toSwitch = 1:length(shape);
   toSwitch(1)=2;
   toSwitch(2)=1;
   shapeStr = regexprep(num2str(shape),'\s{2,}','.');
   filename = ['./matbin/',varStr,'_',shapeStr,'_',type,'_','.matbin'];
   fid = fopen(filename,'w');
   fwrite(fid,permute(var,toSwitch),type);
   fclose(fid);
end

对于python写一个文件matlabIO.py

matlabIO.py

Make the height to be less than 400pixel

python  


# author:wujinnnnn@qq.com
import numpy as np
import os
import subprocess

class getMatBin():
	def __init__(self):
		if not os.path.exists('matbin'):
			raise(Exception("no matbin folder!!!"))

		lsInfo = subprocess.Popen('ls matbin', shell=True, stdout = subprocess.PIPE).stdout.readlines()

		if len(lsInfo)==0:
			raise(Exception("no data files!"))

		print '\n\n\n'
		loadedVars = []
		for eachStr in lsInfo:
			tempStrArray = eachStr.split('_')
			if len(tempStrArray)!=4:
				print "error file format of %s" %(eachStr)
				continue

			varName = tempStrArray[0]
			varSize = tempStrArray[1].split('.')
			varSize.reverse()
			shapeArray = []
			for eachTerm in varSize:
				shapeArray.append(int(eachTerm))

			typeStr = tempStrArray[2]
			shapeArray[-2], shapeArray[-1] = shapeArray[-1], shapeArray[-2]
			exec('self.'+varName + '=' + 'np.fromfile(' +
					'"./matbin/' + eachStr[0:-1] + '", dtype=np.' +
					typeStr + ')')
			#pdb.set_trace()
			exec('self.'+varName + '=' + 'self.'+varName + '.reshape(' + 
					str(shapeArray) + ')')
			print 'variable ' + varName + " loaded!!"

def toNpBin(var,varName):
	if not os.path.exists('npbin'):
		os.mkdir('npbin')
	#else:
		#os.system('rm ./npbin/*.npbin')

	typeStr = str(var.dtype)
	if typeStr=='float64':
		typeStr='double'

	shape = list(var.shape)
	#pdb.set_trace()
	if len(shape)==1:
		shapeArray=[shape[0],1]
	else:
		#shape.reverse()
		shapeArray=shape
		#shapeArray[0], shapeArray[1] = shapeArray[1], shapeArray[0]
		filename = './npbin/'+varName+'_'+str(shapeArray)[1:-1].replace(', ','-')+ \
					'_' + typeStr + '_' + '.npbin'
		var.tofile(filename)

	print 'write ' + varName +' done!'

好了,现在进入实战:

在matlab中有一个2×3×4×5的多维数组A,我的一个逻辑单元是一个2×3的矩阵,用A(:,:,m,n)表示,我要把他弄到numpy中

在matlab中执行toMatBin(A,’varA’)

在python中执行mat = getMatBin(),则mat.varA.shape=(5,4,2,3),这时的一个逻辑单元是2×3的二维数组,用mat.varA[n-1,m-1,:,:]表示

在numpy中有一个B.shape=(2,3,4,5)的多维矩阵,我的一个逻辑单元是一个4×5的二维数组,用B[n-1,m-1,:,:]表示,我要把它弄到matlab中去

在python中执行toNpBin(B,’varB’)

在matlab中执行np = getNpBin(),则np.varB是一个(4×5×3×2)的多维数组,这时一个逻辑单元是4×5的矩阵,用np.varB(:,:,m,n)表示