优化一段代码 - mathematica吧

吧务

level 9

草红样楼主

{$, {data}} = Do[ With[{r = Sqrt[x^2 + y^2], th = ArcTan[x, y]},
If[(Cos[th - r] - Sin[th]) (r^4 - 2 r^2 Cos[2 th + 2.4] + 0.9) + (0.62 r)^1000 < 0, Sow@{x, y}]],
{x, -2, 2,0.005}, {y, -2, 2, 0.005}] // Reap;
(*ListPlot[data]*)
我的笔记本上运行大概17s，matlab版大概0.3s，貌似不好编译

2013年09月06日 13点09分 1

level 10

mm_酱

data = Compile[{}, Table[Evaluate@ With[{r = Sqrt[x^2 + y^2], th = ArcTan[x, y]}, If[(Cos[th - r] - Sin[th]) (r^4 - 2 r^2 Cos[2 th + 2.4] + 0.9) + (0.62 r)^1000 < 0, 1, 0]], {x, -2, 2, 0.0057}, {y, -2, 2, 0.0057}]];data[] // MatrixPlot

2013年09月06日 14点09分 3

xzcyr

把Evaluate去掉更快……

2013年09月07日 05点09分

乡绅地鼠

学习了.

2013年09月08日 00点09分

吧务

level 9

草红样楼主

编译出错原因是ArcTan的问题，试试：ArcTan @@@ Tuples[{0, 0.}, 2]
下面编译后运行时间约2.5s.
cf = Compile[{},
Do[With[{r = Sqrt[x^2 + y^2], th = ArcTan[x + 10^-9, y]},
If[(Cos[th - r] - Sin[th]) (r^4 - 2 r^2 Cos[2 th + 2.4] + 0.9) +
(0.62 r)^1000 < 0, Sow@{x, y}];
],
{x, -2, 2, 0.005}, {y, -2, 2, 0.005}]
];
{$, {data}} = cf[] // Reap;

2013年09月07日 12点09分 4

xzcyr

被你这么一说，好像@situxuming 之前讨论过这个？

2013年09月07日 13点09分

吧务

level 15

xzcyr

仔细一想，我们的方向是不是错了？：
With[{r = Sqrt[x^2 + y^2], th = ArcTan[x, y]}, RegionPlot[(Cos[th - r] - Sin[th]) (r^4 - 2 r^2 Cos[2 th + 24/10] + 9/10) + (62/100 r)^1000 < 0, {x, -2, 2}, {y, -2, 2}, PlotPoints -> 67]] // AbsoluteTiming
Matlab加上出图要几秒？

2013年09月08日 06点09分 5

草红样

这样更方便，但是要对比需要等价的代码，比的都是散点图，matlab里没有画不等式的函数

2013年09月08日 06点09分

草红样

RegionPlot[ Evaluate@With[{r = Sqrt[x^2 + y^2], th = ArcTan[x, y]}, (Cos[th - r] - Sin[th]) (r^4 - 2 r^2 Cos[2 th + 24/10] + 9/10) + (62/100 r)^1000 < 0], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 67] 这样稍快一点

2013年09月08日 06点09分

草红样

matlab纯计算时间0.3s，plot画图也稍快

2013年09月08日 06点09分

吧务

level 9

草红样楼主

这个速度比matlab更快0.25s：
r = Range[-2, 2, 0.005];
bin = Compile[{},
With[{y = r},
Table[With[{r = Sqrt[x^2 + y^2], th = ArcTan[x, y]},
UnitStep@ Minus[(Cos[th - r] - Sin[th]) (r^4 - 2 r^2 Cos[2 th + 2.4] +
0.9) + (0.62 r)^1000]], {x, y}] // Flatten]];
pts= Pick[Tuples[r, 2], bin[], 1]; // AbsoluteTiming
ListPlot[pts]

2013年09月08日 06点09分 6

xzcyr

……怎么想到的？而且貌似这里如果不利用几个函数的Listable属性的话，甚至没法编译？

2013年09月09日 02点09分

草红样

回复 xzcyr :matlab快就是因为用了向量化的操作（就是Listable），mma中要速度的话尽可能地也要用，还有Select和Position都是有点低效的函数，Pick，UnitStep，Unitize，Sign，Clip等需要多用

2013年09月09日 03点09分

flaycatc

为啥mma12运行都要0.3秒+，难道mma越更新越慢了么

2020年02月22日 08点02分

草红样

@flaycatc 应该是遇到奇点了，ArcTan那里改成ArcTan[x, y + $MachineEpsilon]

2020年02月26日 06点02分

吧务

level 9

草红样楼主

matlab代码：
tic;
[x,y]=meshgrid(-2:0.005:2);
r=(x.^2+y.^2).^0.5;
th=atan2(y,x);
i=(cos(th-r)-sin(th)).*(r.^4-2.*r.^2.*cos(2.*th+2.4)+0.9)+(0.62.*r).^1000<0;
toc;
plot(x(i), y(i), '.k')
axis('equal')

2013年09月08日 06点09分 7

level 6

zeroclyy

学习！

2013年09月08日 13点09分 8

吧务

level 15

xzcyr

基于这个帖子里的经验：http://mathematica.stackexchange.com/a/63372/1871
cf = Compile[{},
Table[With[{r = Sqrt[x^2 + y^2], th = ArcTan[x, y]},
If[(Cos[th - r] - Sin[th]) (r^4 - 2 r^2 Cos[2 th + 24/10] +
9/10) + (62/100 r)^1000 < 0, 1, 0]], {x, -2, 2,
0.005}, {y, -2, 2, 0.005}], CompilationTarget -> "C",
RuntimeOptions -> "Speed"];
data = cf[]; // AbsoluteTiming
ArrayPlot@pts

这个解法的速度仅略逊于6楼的解法，但是更加简明。一个值得注意的有趣之处是，C编译器把奇点给绕过去了，具体理由不明……

2015年01月03日 06点01分 10

吧务

level 9

草红样楼主

补充两个以前画的
Plot版：
Plot[{Sqrt[1-x^2],-Sqrt[1-x^2],Sqrt[x-x^2],-Sqrt[-x-x^2],Sqrt[0.005-(x-0.5)^2],-Sqrt[0.005-(x-0.5)^2],Sqrt[0.005-(x+0.5)^2],-Sqrt[0.005-(x+0.5)^2]},{x,-1,1},Axes->False,PlotStyle->Black,
Filling->{
2->{{4},Black},2->{Axis,Black},3->{Axis,Black},4->{Axis,White},5->{{6},White},7->{{8},Black}
},AspectRatio->Automatic]
Graphics版：
Graphics[{Disk[{0,0},1,{Pi/2,(3 Pi)/2}],Disk[{0,1/2},1/2],{White,Disk[{0,-(1/2)},1/2]},{White,Disk[{0,1/2},0.1]},{Disk[{0,-(1/2)},0.1]},Circle[]}]

2015年11月08日 04点11分 12