建模编程擂台系列之一
mathematica吧
全部回复
仅看楼主
level 12
eyofdu 楼主
@waynebuaa 同学发了一个编程擂台系列,但是后面有评论说题目不太看快速解决的能力。
正好最近在做Mathematica教程(大概开学了之后会陆陆续续发出来偶~),有一些题目,和大家分享一下,顺便看看有没有更好的解法吧。
题目1:阈渗模型与果树种植。
果农在种植果树的时候需要关注许多问题,问题之一就是果树的传染病,本题旨在用Mathematica辅助果农找到最适合的栽培方法。
我们知道如果果树种植的较为密集,那么如果一颗果树染病,很快就会不可遏制的扩散到整个群落导致颗粒无收。若种植稍微稀疏一些则可以防止整体灭族,因为果树染病后传染给别的果树的几率较小,而在一段时间内果树因病会死亡,所以也就不会再把病传染下去了。但是假如过于稀疏,果农收入又会降低。所以需要来寻找一个种树的最适间隔,使得在免遭灭族的危险的同时可以让果农收入极大化。
模型的基本假设如下:
果树被种在一个很大的区域内,成正方格点形的培植。
初始条件下有0.1%的树(数量可有小波动)很不幸地感染上了某种植物传染病
传染病的传播是有特点的,在某棵病树旁边,距离该树x的另一颗树每天受到感染的几率为0.1 E^(-x^2),如果该树受到两颗病树的影响,几率相加。
一颗病树会在50天内死亡并不再传染,也不会再生。
请写出可以模拟以上过程的程序,并且在不同的种树间隔下分别多次模拟以分析间隔对最终死亡树的比例的影响并作出(死亡比例~间距)图,你发现了什么?尝试用它来指导果树的栽培。
2016年01月28日 06点01分 1
吧务
level 12
楼主的做法也发出来交流下吧
2016年01月31日 12点01分 3
level 12
eyofdu 楼主
dim = 50;
initp = 0.01;
p = 0.1;
lifespan = 50;
trees = Array[
RandomChoice[{initp, 1 - initp} -> {2 lifespan + 1, 0}] &, {dim,
dim}];
ArrayPlot@trees
pos[trees_, dis_] :=
Total[p E^(-(
#[[1]]^2 + #
[[2]]^2)*dis^2) RotateRight[
1 - UnitStep[-trees], #] & /@
Select[Tuples[Range[Floor[-2 dis], Ceiling[2 dis]], 2],
EuclideanDistance[{0, 0}, #] <= Ceiling[2 dis] &]];
f[-1, _] := -1;
f[x_ /; x >= 1, _] := x - 2;
f[0, x_] :=
RandomChoice[{If[x >= 1, 1, x],
If[x >= 1, 0, 1 - x]} -> {2 lifespan + 1, 0}]
SetAttributes[f, Listable]
t = SessionTime[];
res = Table[
Total[Flatten@
NestWhile[f[
#, pos[#
, dis]] &,
trees, (Total@Flatten@UnitStep[# - 1] > 0) &]]/-dim^2, {dis,
0.3, 2.1, 0.2}]
这是一个简单而慢速的版本,同样是在1.5处跃变。
2016年01月31日 15点01分 6
level 6
用 256*256 的矩阵表示树的状态。健康 (-1);感染 (0 ~ 49);死亡 (-2)。当矩阵的最大值小于零,则终止迭代,计算存活的比例。
用 ListConvolve 可以加速感染概率的计算。
ParallelTable[{l,Total[NestWhile[
#-52UnitStep[#
-50]&[#+UnitStep[ListConvolve[.1Exp[-l^2Array[Norm[{#
#}-5.]^2&,{9,9}]],UnitStep@#
,{{5,5}}]-1~RandomReal~{256,256}](1-Unitize[
#+1])+UnitStep@#
]&,UnitStep[.001-1~RandomReal~{256,256}]-1,Max@#>=0&]+2,2]/256.^2},{l,1.2,1.7,.02}];
ListPlot[%]
2016年02月01日 09点02分 7
这个做法太巧妙了!用一个矩阵表示出了所有树的状态,那几个Unitize和UnitStep看了好久才看懂,学习了!
2016年02月02日 02点02分
精妙+1!我的那个的各类树表示方法还是弗如远甚……ListConvolve也是用的好!
2016年02月02日 23点02分
看了好久终于看明白了,漂亮!学习了!
2016年02月05日 13点02分
[勉强]哪位看懂的大哥能否来个分解加注释的版本呀?[疑问]
2016年02月06日 00点02分
level 12
顶起来
2016年02月06日 01点02分 9
level 12
eyofdu 楼主
在本帖子继续好了,希望有人能看见。
那么既然大家都要过年,系列题目就先不来了,来一道简单有趣的题目。但是有没有看起来如@新泽西管风琴 上一题那样好的解法呢?我不知道。
这次题目是self avoiding walks,一个人从某一个格点出发,在格点上随机行走。但是我们要求它的行走路线是不自相交的。研究这个路线的特性(开放~)
容我偷个懒,春节过后再来系列题目。
假如大家想提前做的话可以自己先百度一下DLA和CCA
2016年02月06日 08点02分 10
level 12
eyofdu 楼主
既然是开放题目,那么我先放上来我的解答也无所谓。
没想到什么好的方法,先这样吧……
size = 200;
loops = 10000;
l = {};
lg = {};
pbc[x_] := x - Floor[x - 1, size];
SetAttributes[pbc, Listable];
f[{l_, pos_, c_}] :=
Module[{dir = pbc[pos + # & /@ {{1, 0}, {-1, 0}, {0, 1}, {0, -1}}],
n, r}, n = Pick[dir, l[[##]] & @@@ dir, 0];
If[Length@n != 0,
r = RandomChoice@n; {{ReplacePart[l, r -> c], r, c + 1},
True}, {{l, c}, False}]]
t0 = SessionTime[];
Dynamic[TemplateApply[
StringTemplate[
"Time Used: ``\nFinished Loops: ``\nComplete Percntages: ``%"], \
{SessionTime[] - t0, t, 100 t/N@loops}]]
Do[t = k;
res = NestWhile[
f[#[[1]]] &, {{SparseArray[{}, {size, size}], {1, 1}, 1},
True}, #[[2]] &][[1]];
AppendTo[l, res[[2]]];
AppendTo[lg, res[[1]]], {k, loops}]
(*为了开动态,这里没有用 Reap 和 Sow *)
Histogram[l, 100, AxesLabel -> {"Steps", "Count"}]
tg = Total@Sign@lg;
dat = ArrayFlatten[ConstantArray[tg, {2, 2}]];
MatrixPlot[dat]
ListPlot3D[dat,
PlotRange -> {{size/2, 3 size/2}, {size/2, 3 size/2}, Full}]
2016年02月06日 08点02分 11
level 12
eyofdu 楼主
emmm,既然有人要求写一个注释版本,那么我就写一个吧(我的想法和@新泽西管风琴自认为有些相似,所以秒懂,我猜理解的也更好一些吧~)。如果有误,请轻喷~
代码:
ParallelTable[{l,Total[NestWhile[
#-52UnitStep[#
-50]&[#+UnitStep[ListConvolve[.1Exp[-l^2Array[Norm[{#
#}-5.]^2&,{9,9}]],UnitStep@#
,{{5,5}}]-1~RandomReal~{256,256}](1-Unitize[
#+1])+UnitStep@#
]&,UnitStep[.001-1~RandomReal~{256,256}]-1,Max@#>=0&]+2,2]/256.^2},{l,1.2,1.7,.02}];
首先,我们要知道这里的数据是怎么存储的:这里有一个列表,里面存储着各个树的生死状态:-2代表这个位置是死树,-1代表这个位置是健康的树,正的数值代表这颗树还活着,数值体现了已经活了多少回合。
在了解了基本的数据存储方法之后我们分开来看:
首先,最内层,先记住这里的#就代指一个表征各个树状态的列表啦。:
ListConvolve[.1Exp[-l^2Array[Norm[{#
#}-5.]^2&,{9,9}]],UnitStep@#
,{{5,5}}]
UnitStep@# 的意思是把所有生病的树搞成1,别的搞成0。在此之后,我们用ListConvolve把每一个点变成旁边一个9*9的区域内每一个点对该点传染概率的贡献(所有变化之前是0的点要么是死掉了,要么没生病,都没有贡献,做卷积的时候贡献也正好是0,假如是生病的,那么就会是1,符合要求)。
获得了每个点的生病几率之后肯定要判断它是否会生病啊,于是下面做如下操作:
***前面那一堆***-1~RandomReal~{256,256}
也就是-RandomReal[1,{256,256}]
在此之后做一个UnitStep——正的变成1,否则是0,这就意味着假如RandomReal给出的数小于几率的数值,就会是1,也就是染病啦~
但是染病这个玩意只对健康的树有效啊,所以后面乘上了Unitize[#+1],意味着只有#上的值为-1,也就是代表健康的树才会为1.别的都是0,也就是乘完了是0啦。
然后再加上UnitStep@#和#,也就是原来生病的一定是多活了一回合了啦~,别的该怎么样怎么样啦。
后面再对这个结果进行操作:

#-52UnitStep[#
-50],假如没有大于等于50,UnitStep乖乖的是0,维持原来结果,一旦破了50,一下减52,一夜回到解放前。完成操作!
最后NestWhileList这个结果直到最大的值是-1为止(只有死掉的和没生病的)
取初始条件为以0.001为几率的列表。
再往后求一个和,获得有多少死掉的,外面套一个ParallelTable造个列表。画个图呗
2016年02月07日 05点02分 12
非常感谢!其实这个代码可读性并不好,解释起来还是比较困难的.
2016年02月07日 15点02分
这么详细,太感谢你啦
2016年02月07日 22点02分
[真棒][真棒][花心]
2016年02月07日 22点02分
@新泽西管风琴 是的,但是结构很好,设计很精妙,赞一个!能加个好友吗?
2016年02月08日 00点02分
level 6
在处理二维点的时候,用复数代替列表会让代码变短,很多时候也能提高效率。下一步位置的计算是在 NestWhile 的判断表达式里做的。
p 为当前位置,i 计数,d 为下一步可能的位置,max 为步数限制,循环的变量是路线。
Block[{p=0,i=0,d,max=1024},NestWhile[ReplacePart[#,++i->(p=RandomChoice@d)]&,0~Table~{max},Length[d=Complement[p+{1,-1,I,-I},#~Take~i]]>0&,1,max];i]~ParallelTable~{10000};
Histogram[%]
在我的笔记本上的运行时间是 3.3s。
2016年02月08日 03点02分 14
厉害厉害!
2016年02月08日 05点02分
但是在走的步数数量很多之后这个算法不会在搜索那四个位置的时候比较慢吗?相比之下在列表中取确定位置的值会更快吧[呵呵]
2016年02月08日 07点02分
或者association是不是也会快?
2016年02月08日 07点02分
@eyofdu 确实如此,用SparseArray存储已经走过的点其实是一种哈希表的方法。但是实际总的路径并不很长,哈希表的方法并不一定会有效率的提升。
2016年02月08日 11点02分
level 12
eyofdu 楼主
在系列题目开始之前再插入一段有趣的东西:
迷宫游戏和走迷宫
我们现在需要来制造最简单的迷宫:迷宫由一个2D或者3D列表来表征,0代表道路,1代表墙。
以2D情况为例,{{1,1,1,1,1},{1,0,1,0,1},{1.0,1,0,1},{1,0,0,0,1},{1,1,1,1,1}}再加上起点设定为{2,2},终点设定为{4,4},这就是一个最简单的迷宫。
首先,一个简单地问题:如何走迷宫?
两方面考虑:
1. 单纯要求快
2. 单纯要求代码简明
由于代码都不长,都写下好了。
第二个问题较为困难:如何制造迷宫?
我们为了保证迷宫不那么弱智,有如下几个要求:
1. 整个地图外边一圈墙。
2. 起点{2,2},终点{l-1,l-1},l为设定的维数。
3. 所有的0必须全部联结(不可以出现一个区域中的0完全不和别的区域相邻)
4. 从起点到终点路径唯一
5. 不可以出现只与一个0相连的1(例如:{{1,0,1},{1,0,1},{1,1,1},{1,1,1}}是不允许的,因为{3,2}位置的1只与一个0相连。这意味着这一条道路可以进一步的延伸而不产生其余影响)
BTW:不要最后生成的路径都是只能向右向下的偶~迷宫迷宫,必须要有难度啊!
最后拿上面的走迷宫程序走一下生成的迷宫,别待会儿自己的程序走不出来自己的迷宫偶~
2016年02月11日 08点02分 15
可能有点意思的迷宫还有一些要求,慢慢在评论里面加。
2016年02月11日 08点02分
2016年02月11日 08点02分
level 6
迷宫的本质是一个无向图。那么问题分为两部分:1.将以矩阵描述的迷宫转换成图;2.找到这个图中连接两个顶点的最短路径。为了简便,矩阵中的 1 表示空地,而 0 表示墙,可以在空地之间向四个方向走动。
函数中,m 表示迷宫,s 表示起点 (source),t 表示目的地 (target):
f[m_,s_,t_]:=FindShortestPath[Graph@Flatten@MapAt[Reverse,#<->#+{1,0}&/@Position[Rest@
# Most@#
,1]&/@{m,Transpose@m},{2,;;,;;}],s,t]
测试:
maze={{1,0,0,1,1,1,0,1},{1,0,1,1,0,0,0,1},{1,0,0,1,1,1,0,1},{1,1,1,1,0,1,1,1},{0,0,0,0,0,1,0,1},{1,1,1,1,1,1,0,1},{1,0,1,0,0,0,0,0},{1,1,1,1,1,1,1,1}};
f[maze,{1,1},{8,8}]
输出 {{1,1},{2,1},{3,1},{4,1},{4,2},{4,3},{4,4},{3,4},{3,5},{3,6},{4,6},{5,6},{6,6},{6,5},{6,4},{6,3},{7,3},{8,3},{8,4},{8,5},{8,6},{8,7},{8,8}}
2016年02月11日 12点02分 16
WatershedComponents与Thinning这两个函数都可以走迷宫,而且所有通路都可以找到,第一个超强。
2016年03月01日 14点03分
@jiaoeyushushu 嗯,确实
2016年03月02日 02点03分
回复 jiaoeyushushu :是的[呵呵]NeatExamples
2016年03月02日 13点03分
回复 jiaoeyushushu :Thinning和Pruning的话似乎只能做比较简单的迷宫?多个分支似乎就跪了,因为Pruning写的是去除最外分支
2016年03月02日 13点03分
level 12
eyofdu 楼主
思路和@新泽西管风琴几乎完全一样。但效率低了不少(1倍)(用了RotateRight而不是Most[
#]*Rest[#
])
放上来一个新的迷宫:
entrance = {2, 2};
exit = {8, 5};
grid = {{1, 1, 1, 1, 1, 1, 1, 1, 1}, {1, 0, 0, 0, 0, 0, 0, 0, 1}, {1,
0, 1, 0, 1, 0, 1, 0, 1}, {1, 0, 0, 1, 0, 0, 1, 0, 1}, {1, 1, 0, 1,
0, 1, 0, 0, 1}, {1, 0, 0, 0, 1, 0, 0, 0, 1}, {1, 0, 1, 0, 0, 1,
1, 1, 1}, {1, 0, 0, 1, 0, 0, 0, 0, 1}, {1, 1, 1, 1, 1, 1, 1, 1,
1}};
ArrayPlot@grid
FindPath[Join @@
Function[{t}, # <-> # - t & /@
Position[(1 - grid)*(1 - RotateRight[grid, t]), 1]] /@ {{0,
1}, {1, 0}}, entrance, exit][[1]]
2016年02月11日 14点02分 17
level 12
eyofdu 楼主
有人解决了制造迷宫的问题了吗?
2016年02月13日 15点02分 18
2016年02月14日 02点02分
level 12
eyofdu 楼主
又查了查,迷宫问题有很多种特定的算法……
http://www.wolfram.com/mathematica/new-in-10/enhanced-graphs-and-networks/maze-generator.html
这就是一个例子。但是在这个例子之中错误路径中也不会有环。
https://en.wikipedia.org/wiki/Maze_generation_algorithm
这是所有的算法。
2016年02月15日 07点02分 19
赞
2016年02月23日 05点02分
level 12
eyofdu 楼主
2016年02月15日 07点02分 20
1 2 尾页