文不可无观点,观点不可无论据。

转载请注明出处

MATLAB功能强大,编程方便,是国际广泛使用的计算软件。目前已有很多书籍介绍其在工程上的应用,但很少有从程序设计语言的角度写的书或文章。

在使用matlab时,总是会碰到这种情况:已经到了嗓子边,但就是想不出来、吐不出来的函数。这时候,matlab帮助里的"see also"可以帮助你。

譬如忘记了两个集合的交集的命令,但依稀记得并集函数为union,于是help union,结果发现其See also中intersect函数豁然在列。也许有时一次see also找不到函数,但点击后跳到另一个函数的帮助,此时的see also中可能就有需要的命令。

有个六度空间假说,世界上任意两个人之间可通过不超过六个人的相互介绍认识。也就是A和F之间可能不认识,但总能找到A认识的B、B认识的C,最后直至E认识的F。笔者突然想研究一下,see also是否有这种规律。自然,是通过MATLAB编写程序(32位windows下2010b完全安装版,其它操作系统或版本的结果数据不相同)。

see also直接引用统计

MATLAB函数统计

see also数量

函数总数

see also数量

函数总数

引用函数数目统计图

直接引用指本命令引用其它命令的数量,统计表明:

有169个函数不包含任何引用,譬如why、rank、waitfor、pcode、minus等;

大部分函数引用数目集中在1~5个;

引用最多的是specgraph,它是一个目录;

引用最多的函数是hdf,共引用其它19个函数,但看起来都是hdfdf24、hdfdfr8等不同界面的兄弟函数,其它的如ode15i(18次)、ode45(17次)、single(17次)也是类似;

plot(16次)引用了大量和它相关的绘图函数,是真正意义上的引用大户。

see also间接引用统计

从任意函数入手,通过引用的互相跳转,能找到多少函数呢。这包含两部分:一是从其他函数出发,多次跳转能找到A函数,这样的函数共有多少个;二是从A出发多次跳转,能搜寻到多少数量的其他函数。生成的统计图如下:

函数过渡关系(三轮跳转)

函数过渡关系(不计跳转数目)

直接看结论:

经过不超过3轮跳转(6轮跳转太多了),最容易找到的函数为function_handle,共有191个函数最多经过3轮跳转都能找到它;

iofun是最基础的文件夹,从它出发,不超过3次跳转可以发现1139个函数;

最基础的函数为string,从它出发不超过3次跳转可,可以发现411个函数。显然,后续函数和IO将是本公众号的核心;

如不计跳转数目,92%的函数可通过800多个函数找到;

通过59%的函数可最终找到1200多个函数。看起来see also没有六度空间理论那么神奇,但善用see also,绝对会有不一样的收获。

引用网络图

那些神奇的绘图软件

用画图板画图太low?用visio、smartdraw、coreldraw画就高级?都是浮云!什么趁手用什么!有那么几件小众的软件,在处理精确关系、海量数据时就很趁手。

MetaPost: 用过LaTeX,或学过算法的人应该都会知道Knuth。因为写《计算机程序设计艺术》,而且是还没有写完就可以得图灵奖;因为嫌出版社将此书排版的不好,中间花了10年时间开发了TeX和MetaFont。John Hobby的 MetaPost 深受 Metafont 的影响,继承了 Metafont 对直线、曲线、点和几何变换等图形优雅的定义语法的优点。想了解可以看图库“-graph/zoonekynd/metapost/metapost.html”。但建议不学,因为要安装LaTex,还要学习LaTeX一些语法。真要想了解,学Asymptote即可。

Asymptote: Asymptote灵感来自MetaPost,采用类似C++语法,原生的支持多段路径(因此不局限于单连通区域)、填充图案、Gouraud 着色(shading)以及PostScript 图像。感兴趣可以看图库“”。

使用程序及解读

统计引用的程序

1   % 通过See also列出MATLAB帮助间引用关系

2   clc;clear;close all

3   bdata_seealsodone=true; %此步耗时较长,如之前保存过则直接读文件

4   if(~bdata_seealsodone)

5       allfuns=getfunsbytoolboxdir('\toolbox\matlab\');   %获取matlab目录下的列出的函数

6       disp('--------getfunsbytoolboxdir done---------------');

7       allseealsos=cellfun(@getseealso, allfuns, 'UniformOutput',false);    %获取上述函数指向的函数

8       disp('--------getseealso done---------------');

9       save data_seealso  allfunsallseealsos;     % 存盘后直接载入

10       else

11          load data_seealso;

12       end

13       % 得到函数及其指向的邻接矩阵

14       conjmat=genconjmatrix(allfuns, allseealsos);

15       disp('--------genconjmatrix done---------------');

16       % 从邻接矩阵得到每个函数被其它函数指向的总数

17       s=sum(conjmat,2);n=max(s);beref=zeros(n+1,2);

18       fori=0:n, beref(i+1, :)=[i, length(find(s==i))];end

19       beref(find(beref(:,2)==0), :)=[];

20       disp(sprintf('%d %5d %5.2f\n', [beref(:,1:2) beref(:,2)/length(allfuns)*100]'));

21       figure;bar(beref(:,1), beref(:,2));grid on;xlabel('引用函数数目(个)');ylabel('函数个数(个)');xlim([-1 20]);

23       arrimat=graph_arrivemat(conjmat);

24       disp('--------graph_arrivemat done---------------');

25       figure;plot(1:length(arrimat), sum(arrimat,1),'b+');grid on;xlabel('函数编号');ylabel('函数个数(个)');title('多少个其它函数能过渡到本函数');  %每个函数被它引次数

26       figure;plot(1:length(arrimat), sum(arrimat,2),'b+');grid on;xlabel('函数编号');title('本函数能过渡到多少个其它函数');  %每个函数引用其它

28       writedot('graph.dot', allfuns, allseealsos);

第2行,清屏,清除变量,关闭所有图形;

第5行,获取MATLAB目录列出的函数;

第7行,获取每一个函数的引用。其中使用了cellfun,其含义为对于元胞数组每一个数组,依次返回函数操作结果。本函数本身的操作对象为函数,是函数式编程的一种简单形式,关于此种编程方法,将在后续介绍;

第14行,生成函数及其引用的邻接矩阵;

第17行,统计函数指向的总数,其中sum第二个参数为2,表示按行求和;

第21行,画图;

第23行,求出邻接矩阵对应的可达矩阵;

第28行,将函数及其引用的关系按graphviz格式输出,用于绘图。

getfunsbytoolboxdir.m

程序设计思路:通过搜索MATLAB文件夹,同时调用MATLAB帮助,获取文件夹内所有函数。

1   functionallfuns=getfunsbytoolboxdir(toolboxname)

2   % 列出安装目录下toolbox/matlab所有分类名称,再从分类中找到其中所有函数

3   excludes={'.', '..'};   %排除'.'和'..'文件夹

4   allfuns={};

5   dirs=dir([matlabroot,toolboxname]); % matlabroot为安装目录

6   for i=1:length(dirs)

7       d=dirs(i).name     %文件夹名

8       if(strmatch(d, excludes, 'exact')) continue;end  %滤去'.'和'..'文件夹

9       seealso=getseealso(d);      % 从Contents.m文件得到函数参考列表

10           allfuns=[allfuns, seealso, {d}];  % 将列表和文件夹本身保存到函数列表中

11       end

12       allfuns=unique(allfuns);   % 去除相同的项

第5行,matlabroot为MATLAB内置函数,输出MATLAB安装文件夹。Dir命令列出制定目录下所有文件夹,并存储在元胞数组内;

第6行,对文件夹循环;

第7行,获取文件夹名;

第8行,从文件夹名中滤去’.’(当前目录)和’..’(上层目录)。strmatch表示按行访问数组或元胞数组,看需要匹配字符与数组开头是否匹配,如果要精确匹配则使用’exact’关键字;

第9行,得到本目录下的引用。根据MATLAB帮助所述,显示MATLAB路径DIR下每个函数的简要描述,如果目录下存在Contents.m文件,则会列出此文件的帮助;

第10行,将新发现的目录和所有引用增加到所有函数列表中;

第12行,由于函数列表中可能存在相同项,采用unique命令去除相同项。

getseealso.m

1   functionseealso=getseealso(comm)

2   % 列出相应函数的See Also列表,如果没有,则为空

3   if(nargin==0) comm='help';end

4   str=help('-help', comm);      % 获取带超链接的帮助文本

5   n=regexp(str, '([\n]\s+Overloaded methods:|[\n]\s+Reference page in Help browser)');

6   if(isempty(n)) n=length(str);end;str=str(1:min(n));    % 滤掉Overloaded methods等后文本(此文本含超链接)

7   seealso=regexp(str, '<a href="matlab:help.*?>(.*?)</a>',  'tokens');   % 正则表达式匹配超链接文本,超链接文本即为函数

8   seealso=cellfun(@(s) s{:}, seealso,'UniformOutput',false);  %将匹配项中的1×1元胞数组展开

9   end

在设计此程序时,遇到了很多困难。最初的想法是直接过滤See also文本,然后找到其中引用的函数。后来发现,还存在不少例外,如有些文本带括号,括号内有文字说明;对于文件夹,其帮助格式完全不一样等等。在打了很多补丁后,突然灵机一动,MATLAB的help命令使用了.m文件的注释,注释中并没有加超链接,但help中增加了超链接,因此MATLAB的help命令里必然含有一个See also分析器。跟踪后发现help未公布于文档的使用方法,即第4行的调用。实际上,如果help不带有这种方法,就只能对help函数进行简单修改以完成本程序了。

第3行,如果调用时不附加任何参数,则使用默认值,类似于c语言的默认型参;

第4行,调用help,此用法MATLAB文档中没有给出,但跟踪help程序可发现;

第5行,正则表达式,表示发现所有以回车开头,包括1~n个空格字符,然后连接Overloaded methods或Reference page的字符。关于正则表达式,将在后续详细阐述;

第6行,从最早发现字符处滤至结尾;

第7行,再次使用正则表达式,发现超链接文本;

第8行,由于正则表达式搜索结果被压缩到了1×1元胞数组内,此处展开,因此输出结果为n×1元胞数组,数组每项为字符型,表示命令指向的引用。

genconjmatrix.m

程序设计思路:有了引用关系后,很容易填写邻接矩阵。在邻接矩阵中需要用到字符串->数字的对应关系,当然可以在每个地方采用strmatch函数搜索到字符串对应的数目。但还有什么比字典更简便快捷的方法呢?

1   functionconjmat=genconjmatrix(allfuns, allseealso)

2   % 写出所有函数指向函数的邻接矩阵

3   n=length(allfuns);conjmat=zeros(n);

4   dict=containers.Map(allfuns,num2cell(1:n));      % 将数据压缩到字典,方便访问

5   fori=1:length(allfuns)

6       root=allfuns{i};iroot=dict(root);

7       forj=1:length(allseealso{i})

8           leaf=(allseealso{i}{j});

9           if(isKey(dict, leaf)) conjmat(iroot, dict(leaf))=1;end       % 填写邻接矩阵

10           end

11       end

第3行,如果调用时不附加任何参数,则使用默认值,类似于c语言的默认型参;

第4行,将函数->编号压缩为字典,以提高后续访问效率,其中num2cell生成了1~n的元胞数组;

第9行,将引用关系压缩到邻接矩阵。语句中采用了字典以提高效率,MATLAB自带了很多高级的数据类型,因此大大方便了其操作。

graph_arrivemat.m

1   functionQr=graph_arrivemat(A, n)

2   % 计算可达矩阵

3   % 可达矩阵=A+A^2+...+A^n

4   if(nargin==1) n=length(A);end

5   A=single(A);

6   sqrtn=nextpow2(n+1);

7   An{1}=A;

8   fori=2:sqrtn-1

9       i;

10           buf=An{i-1}^2;buf(buf>0)=1;An{i}=buf;

11       end

12       I=eye(size(A));

13       Qr=sumA1ton(n);

14           functionQ=sumA1ton(n)

15           % 计算A+A^2+...+A^n的递归程序,其中An和I变量继承自父程序

16               n;

17               if(n==1)

18                   Q=A;

19               elseif(rem(n,2)==0)

20                   Q=(I+A_power_n(n/2))*sumA1ton(n/2);

21               else

22                   Q=A*(I+  (I+A_power_n((n-1)/2))*sumA1ton((n-1)/2) );

23               end

24               Q(Q>0)=1;

25           end

26           functionQ=A_power_n(n)

27           % 计算A^n的快捷算法

28               bin=fliplr(dec2bin(n));

29               Q=I;

30               fori=1:length(bin)

31                   if(bin(i)=='1') Q=Q*An{i};end

32               end

33               Q(Q>0)=1;

34           end

35       end

程序设计思路:图论业已证明,图的n次可达矩阵可以表示为A+A2+…+An。剩下的就是通过MATLAB计算它。最简单的方法为按循环实现,结果效率非常低。开始笔者的思路是A+A2+…+An=A(I-An)/(I-A),并用一个10×10矩阵测试,但由于截断误差,无法得出正确结果。因此最终使用了分治法,如在n为偶数时,将A+A2+…+An拆分为A+A2+…+An/2 +An/2+1+…+An=(I+An/2)×(A+A2+…+An/2),从而将问题规模减半,依次类推,再对后者应用上述公式,直至规模为1。在计算时,还发现,计算多个An/2时存在重复计算的情形,因此事先将A1、A2、A4、A8…等计算好后存储,在计算An/2时直接使用存储数值组合。在1.80GHz的双核CPU上,使用分治法后,程序运行时间由小时级降低到62s,使用存储后,运行时间进一步降低到46s,使用单精度存储后,运行时间降低到18s。

第5行,使用单精度储存数据。毋庸置疑,由于邻接矩阵为整型的、存在较多0元素的矩阵,采用更精简(如整型、单精度、稀疏矩阵等)的格式存储数据可提高计算效率。不幸的是,MATLAB不支持整型乘法、也不支持单精度稀疏矩阵等,经测试使用单精度效果最好;

第6行,获取2p>n的第一个p,用于计算A1、A2、A4、A8…;

第12行,保存一个单位矩阵待用;

第13行,调用sumA1ton子程序;

第14行,嵌套子程序,嵌套子程序可继承父程序内的变量,特别是An和I,不用作为函数传入,方便了程序的编写。与普通程序相比,嵌套程序要求在相关函数后加end符,形成function-end对。当然,命名时An和I这种简短的符号是不推荐的,除非是特别短小的函数;

第17~23行,分治法递归程序,分为递归终止(n为1),n为偶数和n为奇数三种情况;

第24行,将不为0的Q全部置为1,此种方法称为逻辑切片,Q内的输入A为与Q同维数,一般由逻辑操作生成的矩阵,则指定为Q中对应A不为0位置的元素,与下标切片Q(find(Q>0))相比,此种语法更简捷,也运行的快一点;

第28行,dec2bin将十进制数字转换为二进制字符串;fliplr将字符串按从右至做顺序重写;

第30到32行,计算An的快捷程序,实现思路为任意n可写为2m之和,譬如9=8+1(即二进制1001),13=8+4+1(即二进制1101),将n展开后,用之前计算的An乘即可得到。

1   functionwritedot(fn, allfuns, allseealso)

2   % 将函数引用关系写为dot文件,由graphviz/dot生成图形

3   fp=fopen(fn,'w');

4   fprintf(fp,sprintf('#!"c:/Program Files/Graphviz2.26.3/bin/dot.exe" -Temf %s -o %s.emf\r\n', fn, fn));

5   fprintf(fp,'digraph g\r\n{\r\n');

6   fori=1:length(allfuns)

7       if(~isempty(strfind(allfuns{i}, '.')) | ~isempty(strfind(allfuns{i},'/'))) continue;end   % dot中不识别.和/,简单滤去

8       forj=1:length(allseealso{i})

9           if(~isempty(strfind(allseealso{i}{j},'.')) | ~isempty(strfind(allseealso{i}{j}, '/'))) continue;end

10               fprintf(fp,'%s->%s;\r\n',(allfuns{i}), (allseealso{i}{j}));

11           end

12       end

13       fprintf(fp,'}\r\n');

14       fclose(fp);

上面的程序中,涉及了矩阵生成、数据结构、编程方法、计算效率、正则表达式、IO操作、图形绘制等多方面,有一定的准入门槛,初学者可以先行跳过。后续分析其中的方法,同时笔者也提醒自己,以后不能再出现这么长的程序了,不然就没有人看了!

关于帮助的几个小tip

当输入命令时,在输入字符的前面一部分后,按TAB键具有补全功能。此时,如果以此字符开头的函数或变量唯一,则直接补全函数或变量,否则弹出列表以供选择。此功能除提高输入效率外,也可以用来作为命令提示。譬如想看看MATLAB有无进制转换函数,大致感觉十六进制为HEX,则输入hex,然后按两次TAB键,hex2dec就赫然在列了。再执行help hex2dec,其它进制转换均在see also中列出了。具体详见MATLAB / User’s Guide / Programming Fundamentals / Programming Tips / Command and Function的帮助。

代码编辑(edit)

MATLAB是收费软件吗?no no no,MATLAB只收了1%,其它的都是开源的!EDIT FUN将弹出函数对应的代码。譬如在写这篇时笔者发现,使用命令help help可以出现超链接,但t=help(‘help’)则不出现。于是edit help设置断点后跟踪,发现helpProcess这个类的初始化函数中,如果无返回值,则加超链接,否则不加。但同时,如果help函数的第一个参数是’-help’,’-helpwin’或’-doc’,则可以在返回值中显示超链接。还有,笔者对MATLAB如何将self-documenting文档的see also加上超链接很感兴趣,于是继续一通跟踪,发现在MATLAB安装目录/matlab/helptools/+helpUtils/@helpProcess/hotlinkHelp.m中,它将帮助文档的最后一个See Also作为开头,以overloaded或者note或者两行连续回车或者文本结束作为结束,搜索到的内容作为see also字符,然后识别出关键词,再增加超链接标签。当然,很多时候大家不会这么无聊,但当你的程序出错时,有些时候,MATLAB仍会将你引导到相应的M文件中。

编写自己的帮助

在命令窗口执行edit help,会发现help.m文件中出现的一堆注释恰好与help help输出结果一致。这不是巧合,而是MATLAB的设计,并赋予它很眩的名字,如MATLAB称为self-documenting,不要小看前4个字母,现代程序语言越来越注重这种self设计,它代表了一种内省或自省的威力。一个东西如果能看懂自己,一定很强大。除了self-documenting之外,帮助本身也是自省的,因为程序的帮助本身也是程序。试想,有多少种语言能将在程序中将帮助输出呢?

每一个函数下连续的以%开头的注释,都将作为文档出现,包括MATLAB自带,他人编写,或者自己编写的程序。一个典型的程序为:

1   function f = fact(n)

2   % Compute a factorial value.

3   % FACT(N) returns the factorial of N,

4   % usually denoted by N!

5   % See also factorial, prod.

7   % Put simply, FACT(N) is PROD(1:N).

8   f = prod(1:n);

第1行,函数定义行,定义函数名称,输入和返回值;

第2行,H1行,对于程序的一句话总结,当在整个文件夹使用help,或使用lookfor时显示;

第3、4行,连续的以注释字符%开头的行,将作为帮助文档在help命令中出现;

第5行,See also行(see和also用一个空格隔开,两个单词可任意大小写),其后词汇均会被help识别为超链接;

第7行,由于第5行未以%开头,它仅为单纯的注释,而不是帮助文本;

第8行,函数体。

将以上文件保存到fact.m后,输入help fact命令将显示

Compute a factorial value.

FACT(N) returns the factorial of N,

usually denoted by N!

See also factorial, prod.

自己编写的程序怎么在帮助浏览器中显示公式和图形呢?这需要单独编写文件,已与self-documenting无关了。

往期文章:

微信扫一扫

关注“理念世界的影子”

版权声明:

本文是"洞穴之外"作者原创文章,欢迎转载,须署名并注明来自“理念世界的影子”公众号。