[讨论]:cody基础数据操作题---Bullseye Matrix

bainhome

祁工发了一道题目，链接：cody基本数据操作练习：diff、正则应用及其他
我也发一道cody的题目，大家可以讨论一下做法：
Given n (always odd), return output a that has concentric rings of the numbers 1 through (n+1)/2 around the center point. Examples:

1. Input n = 3
   
2. 
   
3. Output a is
   
4. [ 2 2 2
   
5. 2 1 2
   
6. 2 2 2 ]
   
7. 
   
8. Input n = 5
   
9. 
   
10. Output a is
    
11. 
    
12. [ 3 3 3 3 3
    
13. 3 2 2 2 3
    
14. 3 2 1 2 3
    
15. 3 2 2 2 3
    
16. 3 3 3 3 3 ]

复制代码

题目相对简单，但想要做好也需要合理的构思。
ps，很久不上论坛回答问题，多说两句：
有些初学者对这种数据基础操作题目不以为然，认为没有必要，认为MATLAB是工具，简单会用两下即可，但是这样的思想往往在专业学习后期就带来恶劣后果。要不然就是索要程序→拿来不会用→运行出各类花式低级错误；或者索要程序未果，按自己想象胡编乱造一气；有的生怕别人看出自己什么都不懂，故意把MATLAB和专业问题揉在一起混淆视听，其实如同白痴的你，除了自己谁你也蒙不了(这段时间见了不少这种极品)。

突然想起鹿鼎记里韦小宝在少林寺跟师侄学武功的故事，其实很多人都是这样：一套罗汉拳都打得颠三倒四歪七扭八，这样的你也配学拈花指？
所以，后面会陆续出些这样的题目，大家一起讨论讨论。
做多了你会发现：很多人，甚至很多现在已经在出书的所谓“高手”，程序写出来也有许多值得商榷之处，不是流程布局不合理，就是技术已经落后。
做多了你还会发现：很多以前抽象枯燥、可望不可及的专业难题，随手花上几个小时就可以编出程序完成验算并惠及周围；
你最终发现：没必要”的事情做得多了，“有必要”的事情正慢慢变得简单、有序以及尽在掌握中。
啰嗦了半天，cody上做过见过的兄弟忍一忍，给其他伙计一个星期的缓冲时间，大家一起提高。也有好些伙计本身已经足金足两，不妨直接去cody的深水区游一游，看看世界级MATLAB高手的风采。不知道你们会有啥感受，反正我是很受促进。

nwcwww

先给个取巧的解，不是正途：
bwdistgeodesic(true(n, n), (n^2+1)/2)+1

这是通过图像处理工具箱做的，比较无脑，也是我看到题目后的第一反应。
用bwdist，选项chessboard也能达到相同效果。


提供一个矩阵操作的解，权当抛砖引玉：

1. function ans = bullseye(n)
   
2.      [x, y] = meshgrid(-(n-1)/2:(n-1)/2);
   
3.      max(abs(x), abs(y))+1
   
4. end

复制代码

liuyalong008

1. function a=concentric_rings(n)
   
2. %%
   
3. % 每个正方形上的值都相等，按照上下和左右分别赋值即可
   
4. for i=0:(1+n)/2-1
   
5.     a([(1+n)/2-i,(1+n)/2+i],[(1+n)/2-i:(1+n)/2+i])=i+1;
   
6.     a([(1+n)/2-i:(1+n)/2+i],[(1+n)/2-i,(1+n)/2+i])=i+1;
   
7. end
   

复制代码

抛砖！
等keys

lengyunfeng

楼上两位，你们以后出招的时候麻烦不要学q版，这种一两行精简到极点却不作注释或者说明的代码让人看了眼晕。除了阳春白雪以外，还是适当的加点柴米油盐嘛:L。

lin2009

3楼的代码可以进一步简化，见下面说明。

1. function a = concentric_rings(n)
   
2. %Given n (always odd), return output a that has concentric rings of the numbers 1 through (n+1)/2 around the center point.
   
3. %可以看作是大大小小正方形“重叠”起来。
   
4. % 先给最外围的最大的正方形赋值, 再向中心的正方形赋值。
   
5. 
   
6. for i = (1+n)/2 - 1 : -1 : 0      % 最大正方形的中心点为基准进行排序编码。
   
7.      a([ (1+n)/2 - i : (1+n)/2 + i ], [ (1+n)/2 - i : (1+n)/2 + i ]) = i + 1;
   
8. end

复制代码

qibbxxt

简化一下2楼的代码
meshgrid从二维的角度来说，应该是bsxfun的一个特例

1. function ans = bullseye(n)
   
2. bsxfun(@(a,b)max(abs(a),abs(b)),-(n-1)/2:(n-1)/2,(-(n-1)/2:(n-1)/2)') + 1;
   
3. end

复制代码

bainhome

以上各位基础数据操作都用得相当精妙，其中两个求解基本是我个人心目中的leading solution，佩服。
出这道题有我自己的意图：这题解法非常之多，总结起来有两种大致思路：
1.利用MATLAB现有的特殊矩阵和测试矩阵，用built-in命令实现再构造；
2.老老实实循环，不管是几重，不管用不用各种arrayfun、bsxfun，或者就是for。
思路1中有如下几个比较好的思路：
1).“托普利茨”常对角阵打底，用其他基础矩阵操作旋转、相加等完成计算例如：

1. >> c = toeplitz(1:5)
   
2. c =
   
3.      1     2     3     4     5
   
4.      2     1     2     3     4
   
5.      3     2     1     2     3
   
6.      4     3     2     1     2
   
7.      5     4     3     2     1
   
8. >> c = (c + fliplr(c))/2
   
9. c =
   
10.      3     3     3     3     3
    
11.      3     2     2     2     3
    
12.      3     2     1     2     3
    
13.      3     2     2     2     3
    
14.      3     3     3     3     3

复制代码

2).托普利茨阵与rot90旋转阵的组合：

1. >> n=5;c=toeplitz(1:n)/2
   
2. c =
   
3.     0.5000    1.0000    1.5000    2.0000    2.5000
   
4.     1.0000    0.5000    1.0000    1.5000    2.0000
   
5.     1.5000    1.0000    0.5000    1.0000    1.5000
   
6.     2.0000    1.5000    1.0000    0.5000    1.0000
   
7.     2.5000    2.0000    1.5000    1.0000    0.5000
   
8. >> c + rot90(c)
   
9. ans =
   
10.      3     3     3     3     3
    
11.      3     2     2     2     3
    
12.      3     2     1     2     3
    
13.      3     2     2     2     3
    
14.      3     3     3     3     3

复制代码

3).使用螺旋阵spiral构造，这也是现在cody上真正的leading solution核心流程。

1. c=spiral(n)
   
2. c =
   
3.     21    22    23    24    25
   
4.     20     7     8     9    10
   
5.     19     6     1     2    11
   
6.     18     5     4     3    12
   
7.     17    16    15    14    13
   
8. c=sqrt(c)
   
9. c =
   
10.     4.5826    4.6904    4.7958    4.8990    5.0000
    
11.     4.4721    2.6458    2.8284    3.0000    3.1623
    
12.     4.3589    2.4495    1.0000    1.4142    3.3166
    
13.     4.2426    2.2361    2.0000    1.7321    3.4641
    
14.     4.1231    4.0000    3.8730    3.7417    3.6056
    
15. c=0.5+0.5*c
    
16. c =
    
17.     2.7913    2.8452    2.8979    2.9495    3.0000
    
18.     2.7361    1.8229    1.9142    2.0000    2.0811
    
19.     2.6794    1.7247    1.0000    1.2071    2.1583
    
20.     2.6213    1.6180    1.5000    1.3660    2.2321
    
21.     2.5616    2.5000    2.4365    2.3708    2.3028
    
22. c=ceil(c)
    
23. c =
    
24.      3     3     3     3     3
    
25.      3     2     2     2     3
    
26.      3     2     1     2     3
    
27.      3     2     2     2     3
    
28.      3     3     3     3     3

复制代码

4).mahatton权重距离阵做法：

1. 0.5*(mandist(1:n)+flipud(mandist(1:n)))+1;

复制代码

5).boxdist函数：

1. function a = bullseye(n)
   
2.   a = diag([(n+1)/2:-1:1 2:(n+1)/2]);
   
3. a = boxdist(a)+a;
   
4. end

复制代码

当然还有人用nwcwww所给的图像处理函数的做法，他们基本都是属于第1中思路，即利用特殊矩阵做底阵，然后通过一些矩阵基础操作，或凑或算，达到计算目的
(3)是size最小的leading solution的思想，但这种做法有些取巧，其实那个size最小的用了动态正则替换，把原本18的size降到了11，我认为这种无厘头做法，纯粹属于偏执，不值得提倡，可是人家的正则替换用得真心是牛A<X<牛C，这也不得不承认：

1. regexp('','(?@y=ceil(sqrt(spiral(n))/2+0.5);)');

复制代码

另一种思路最实在，也被很多MATLAB learner很看不起——循环：
1).

1. for i = 0 : n / 2
   
2.     j = i + 1 : n - i;
   
3.     a( j, j ) = ( n + 1 ) / 2 - i;
   
4. end
   

复制代码

2).当然这种做法也有很多变体，比如liuyalong给的：

1. for i=0:(1+n)/2-1
   
2. 
   
3.     a([(1+n)/2-i,(1+n)/2+i],[(1+n)/2-i:(1+n)/2+i])=i+1;
   
4. 
   
5.     a([(1+n)/2-i:(1+n)/2+i],[(1+n)/2-i,(1+n)/2+i])=i+1;
   
6. 
   
7. end
   

复制代码

3).

1. n1=(n+1)/2;
   
2. a=n1*ones(n);
   
3. x=[n1,1,length(a)];
   
4. for i=n1-1:-1:1
   
5.     a(x(2)+1:x(3)-1,x(2)+1:x(3)-1)=(x(1)-1)*ones(2*i-1);
   
6.     x=[x(1)-1,x(2)+1,x(3)-1];
   
7. end

复制代码

上面这个思路和lin2009提供的思想有些许类似。
4).递归计算：

1. function a = bullseye(n)
   
2. try
   
3.     a = ones(n)*ceil(n/2);
   
4.     a(2:end-1,2:end-1) = bullseye(n-2);
   
5. end
   
6. end

复制代码

还有一种思路介于二者之间，无论bsxfun这种逐元素操作还是meshgrid（ndgrid）布点再用max做逻辑判断都属于这种，动作分解如下：

1. c=meshgrid(1:n)
   
2. c =
   
3.      1     2     3     4     5
   
4.      1     2     3     4     5
   
5.      1     2     3     4     5
   
6.      1     2     3     4     5
   
7.      1     2     3     4     5

复制代码

1. c=c-(n+1)/2
   
2. c =
   
3.     -2    -1     0     1     2
   
4.     -2    -1     0     1     2
   
5.     -2    -1     0     1     2
   
6.     -2    -1     0     1     2
   
7.     -2    -1     0     1     2

复制代码

1. c=abs(c)
   
2. c =
   
3.      2     1     0     1     2
   
4.      2     1     0     1     2
   
5.      2     1     0     1     2
   
6.      2     1     0     1     2
   
7.      2     1     0     1     2

复制代码

1. c=max(c,c')
   
2. c =
   
3.      2     2     2     2     2
   
4.      2     1     1     1     2
   
5.      2     1     0     1     2
   
6.      2     1     1     1     2
   
7.      2     2     2     2     2

复制代码

1. c=c+1
   
2. c =
   
3.      3     3     3     3     3
   
4.      3     2     2     2     3
   
5.      3     2     1     2     3
   
6.      3     2     2     2     3
   
7.      3     3     3     3     3

复制代码

bxcfun变体，实际和前面一个道理：

1. b=[0.5*n+0.5:-1:1 2:0.5*n+0.5]
   
2. a=bsxfun(@max,b,b')

复制代码

点评：这两种思路(最后meshgrid+max的做法其实是属于第1种)难说优劣，多数MATLAB Fans偏重于前一种：想必长泡论坛很多用循环写代码的伙计都不好意思在论坛里公开用循环解决问题，怕遭人笑，实在没法避免也要用arrayfun包装一下，以前我也有这种想法，而且很严重...
前两天某论坛里又有个关于90%循环在MATLAB中是否可以避免的争执(也算争执)，我感觉MATLAB代码写得成矢量化固然是功力的体现，可循环这种东西现在一来效率不算低，当然我承认也有较极端的例子反向说明这个观点的绝对和武断，可实际问题求解中有几个问题要循环到10w次以上？吴鹏这种整天跟卫星高清图片过不去的变态鸟人除外。所以我觉得没事for两下也没啥大不了；二来还有一个比较实际的好处：写代码写得老实，该for就for，这样的代码容易看懂、移植性好，前面那个循环的办法size=38，高于leading solution整整27个，但是如果改成C#代码，用5min稍微调一下语法，搞定！至于那个飘逸潇洒的leading solution估计就头疼了。毕竟有几个人喜欢整天带着MCRinstaller到处跑的？现在好多人吭哧吭哧改MATLAB算法为其他语言，都被爽得痛不欲生，矢量化是其中一原因。所以前天我跟祁彬彬qq开玩笑说那个用for写bullseye的兄弟是个实在人，呵呵。不过事先也说好：别较劲，跟我说“你既然这么强调底层，干脆把啥东西都往二进制上写！”
悄悄告诉你：我真想写，可惜水平没你奇葩...
个人观点+一点闲话，大家听听就好，我是答应lengyunfeng写点儿东西捧场。各位水准都很高，simwe虽如今阵容不齐，可如此残阵依旧猛男云集，水底下藏的都是大鱼。这样的问题放在其他论坛估计参与的不会多(祁彬彬试过，压根儿没人搭理这厮，那天结贴的时候写那么长，我估计是看到大家突然搭理，热泪盈眶一时之间没想开，就难得地总结了那多玩意儿，哈哈)。建议兄弟们cody上转转吧！除了lin2009是多软件混成部队，nwcwww和liuyalong你俩明明能玩儿大江大海，还赖在游泳池浅水区里带着救生圈猛练狗刨？

nwcwww

其实这个问题用for loop确实挺好的。简单易懂，而且效率挺高。
矢量化更多是由于早期matlab在循环上惨不忍睹的表现吧，结果现在大家有点走火入魔的感觉，程序里出现循环都羞于见人。

实际上随着JIT的逐渐完善，非嵌套的循环在不调用自定义函数的情况下表现往往已经足够好了。这时又何苦去牺牲程序的可读性和可维护性呢？而且还能少死一些脑细胞。

矢量化代码的一大问题是不好找接盘侠~ 老老实实写出来的程序，别人本来花几分钟略作修改就可以投入使用的。矢量化无疑加大了移植的难度。
有师兄抱怨过，拿到矢量化的程序时，心情就和要破解没脱壳的软件差不多。

qibbxxt

nwcwww 发表于 2013-5-14 00:32
其实这个问题用for loop确实挺好的。简单易懂，而且效率挺高。
矢量化更多是由于早期matlab在循环上惨不忍 ...

从写程序的效率和可读性来说，老老实实写程序，固然是最好的

但是我觉得矢量化是Matlab的特色和魅力所在，如果离开这一点，似乎和别的语言区别不是很大，Matlab有工具箱，人家有各种库可以调用

另外一方面，在cody上面，大家以size最小为目标，所以时间长了，就习惯于代码简短的写法了

一种折中的做法是二者兼顾，比如，有时候，用一个矩阵的每一行和一个行向量操作是，我就习惯了用bsxfun,

这样其实也不用费太多脑筋，反而比写循环效率要高一些

这个问题大家的看法不同，以上的只是我的个人见解

bainhome

呵呵，祁工太谨慎了。
其实点评目的一方面提醒初学MATLAB的learner，向量化应循序渐进，开始还是把程序写实在，各种程序共通之处先写好写懂，先通解后特解，不要一上来数据类型没整明白就二百五地搞向量化，这一点吴鹏和老谢估计有感受，整天猫在那儿回答些四六不沾的弱智问题，诸如非要用arrayfun结合符号函数之类,不是不能这么实现，而是不该这么实现：数值、符号两颗树上的果实，能碰到一起是MATLAB能，不是你能。因此整天琢磨着MATLAB为啥没帮我把所有事都干完的人，建议写好基本内容再思考其他更妥当，尤对写稍大一点论文规模的学生而言，这样做后期有好处；另一方面关于向量化似乎有矫枉过正的趋势，尤其一些学了两天的半吊子更是整天“矢量化”不离口。我的主要意思还是别妖魔化循环，当然你几个伙计不会这样想这样做，尤其祁工你自己不能做例子，因为你对MATLAB中high-level函数的理解在我知道的User中几乎无出其右(应该有几个同水平线)，你现在对高级函数的运用理解境界只有少数人可以达到(我肯定做不到)。但你们好赖有点影响力，你们不写，很多人模仿着都不写...结果最后，写论文导师要求更换语言，一脑门汗全傻眼，才知道自己语言基本功有多龌龊。
另外，不知道各位发现没有：向量化也是Mathworks的战略之一，用习惯MATLAB之后很难转其他语言，至少比其他语言互转要困难，总有“一夜回到解放前”的郁闷。
最后，没有一丝一毫否定向量化的意思，MATLAB与数学语言的天然契合性是它的最大特征。我写东西也都按向量化思路尽量写(祁工最清楚不过，呵呵)，只是提醒初学MATLAB的各位同学，基础打好，向量化思想的转变是个自然而然水到渠成的过程，不需刻意强求。

等有一天你感觉有了，那你就有了。

其实在那时之前，不知不觉你早就有了。

lengyunfeng

我觉得矢量化操作对大部分人而言还是比较困难的（我就是其中一个），因为这涉及到个人的数学功底或者怎么看这个世界的问题，就像一维数据我好理解、二维也将就、三维有点勉强，一到4维就发晕了，再高维的话脑袋就成浆糊了，因为现实世界中可以想象的就只有空间和时间，最高维只能到4维。5维？想像一下，另一个平行的世界还有一个自已在瞎搞闯祸？干，我闲得蛋疼啊，自己的屁事都搞不定，哪有时间去管他！所以个人觉得这全是惯性知识造的孽，搞得现在即使接写过一个向量程序，隔段时间让我再写，我还是会自然而然的写回循环去。我觉得正如bainhome说的那样，对于我们大部分人而言，向量化编程是个渐进学习、自然而然的过程，没必要刻意的去为之。

rocwoods

大家都冒泡了，不吐不快！
其实cody上面偏执狂很多的，好多leading代码别看size小，效率未必高。这个就跟优化一样，看你的目标是啥，以size为最优目标，得到的就是cody上的leading代码，以方便快捷为目标得到的程序肯定可读性强，如果以运行效率为目标得到的和上面两种都不同。
我个人倾向于兼顾效率和易读性。作为铁杆MATLAB fans，以前很痛恨那些不了解MATLAB就贬低MATLAB，停留在过去认识的人。经常在其他论坛，尤其是C++论坛，看到有人说MATLAB就是个小玩具（靠，口气不小），搞不了大的应用，C++爆MATLAB n条街。每每看到这样的言论，气儿就不打一处来。于是，就想方设法证明，使用正确情况下，MATLAB在多种场合，效率完全不输C++。
后来工作多了，写的程序多了，真正认识到每种语言都有其特长后，看到类似哪个语言好坏的争论后就一笑置之了，也懒得和不懂MATLAB的争论了。还是那句话，for循环不可怕，向量化未必效率高，可怕的是“勤拿少取，本来叫一个箱货就可以一趟搞定的搬家，非得打出租车分10趟搞定！”

liuyalong008

以后进行算法比较的时候，可以把运行时间(以及内存开销等)也作为考虑因素进行探讨，代码简洁和高效应该兼具

bainhome

把bullseye matrix(红心阵)玩儿到底，这回在原有基础上扩展一下：
cody原帖链接见这里
为给大家省点儿时间，简单解释一下：
要求给定一个最大的正整数，然后形成先从1始顺序增加，后倒序减小回1的红心矩阵
例1：n=3

1. ans =
   
2. 1 1 1 1 1 1 1 1 1
   
3. 1 2 2 2 2 2 2 2 1
   
4. 1 2 3 3 3 3 3 2 1
   
5. 1 2 3 2 2 2 3 2 1
   
6. 1 2 3 2 1 2 3 2 1
   
7. 1 2 3 2 2 2 3 2 1
   
8. 1 2 3 3 3 3 3 2 1
   
9. 1 2 2 2 2 2 2 2 1
   
10. 1 1 1 1 1 1 1 1 1

复制代码

例2：n=4

1. ans =
   
2. 1 1 1 1 1 1 1 1 1 1 1 1 1
   
3. 1 2 2 2 2 2 2 2 2 2 2 2 1
   
4. 1 2 3 3 3 3 3 3 3 3 3 2 1
   
5. 1 2 3 4 4 4 4 4 4 4 3 2 1
   
6. 1 2 3 4 3 3 3 3 3 4 3 2 1
   
7. 1 2 3 4 3 2 2 2 3 4 3 2 1
   
8. 1 2 3 4 3 2 1 2 3 4 3 2 1
   
9. 1 2 3 4 3 2 2 2 3 4 3 2 1
   
10. 1 2 3 4 3 3 3 3 3 4 3 2 1
    
11. 1 2 3 4 4 4 4 4 4 4 3 2 1
    
12. 1 2 3 3 3 3 3 3 3 3 3 2 1
    
13. 1 2 2 2 2 2 2 2 2 2 2 2 1
    
14. 1 1 1 1 1 1 1 1 1 1 1 1 1

复制代码

函数大致出入口:

1. bullseyeExtended=yourfun(n)

复制代码

大家继续high一把。

liuyalong008

bainhome 发表于 2013-5-22 22:21
把bullseye matrix(红心阵)玩儿到底，这回在原有基础上扩展一下：
cody原帖链接见这里
为给大家省点儿时间 ...

我首先贴一个长的，借鉴nwcwww兄的部分代码

1. function t=yourfun(a)
   
2. n=4*a-3;
   
3. [x, y] = meshgrid(-(n-1)/2:(n-1)/2);
   
4. % 生成矩阵从中心向边上递减
   
5. t=-(max(abs(x), abs(y))+1)+a*2;
   
6. n=a*2-1;
   
7. [x, y] = meshgrid(-(n-1)/2:(n-1)/2);
   
8. %把中心的再替换掉
   
9. t(a:end-a+1,a:end-a+1)=max(abs(x), abs(y))+1;

复制代码

bainhome

发现没有？现在贴代码形容词很多：长的、短的、省力的、容易理解的、丑陋的、暴力的、抛砖引玉的、shamless的、偷懒的、啰嗦的、古典的、矢量化的、实在的、取巧的...
人活着不容易啊...
我贴个偷懒+shamless+抛砖引玉+比较短的：

1. function ans=bullseyeExtend(n)
   
2. (toeplitz(4*n-3:-1:1)+fliplr(toeplitz(4*n-3:-1:1)))/2-2*n+2;
   
3. ans(n:(3*n-2),n:(3*n-2))=ceil(sqrt(spiral(2*n-1))/2+0.5);

复制代码

坐等nwcwww和祁工的暴力美学版。

liuyalong008

贴一个循环，

1. function a=bullseyeExtend(n)
   
2. t=[1:n n-1:-1:1];
   
3. n=4*n-3;
   
4. for i = (1+n)/2 - 1 : -1 : 0      
   
5.      a([ (1+n)/2 - i : (1+n)/2 + i ], [ (1+n)/2 - i : (1+n)/2 + i ]) =t(i+1);
   
6. end

复制代码

lengyunfeng

我也来个循环的金字塔赋值

1. function bullseyeExtended=yourfun(n)
   
2. bullseyeExtended=ones(4*n-3);
   
3. k=[2:n,n-1:-1:1];
   
4. for m=2:2*n-1
   
5.     bullseyeExtended(m:4*n-m-2,m:4*n-m-2)=k(m-1);
   
6. end

复制代码

qibbxxt

1. 老兄都点名了，不得不出来啊

复制代码

1. function ans = yourfun(n)
   
2. [2*n-1:-1:1, 2:2*n-1];
   
3. min(bsxfun(@min,2*n-ans,2*n-ans'),bsxfun(@max,ans,ans'));

复制代码

nwcwww

1. function ans = your_fcn_name(n)
   
2.   2*n-ceil(0.5*sqrt(spiral(4*n-3))+0.5);
   
3.   ans(ans>n) = 2*n-ans(ans>n);
   
4. end

复制代码

借鉴了上面的spiral。 凑出来这个矩阵的那位真是神人。