/
mandelbrot.html
270 lines (250 loc) · 17.5 KB
/
mandelbrot.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
<!DOCTYPE html>
<!-- Generated by pkgdown: do not edit by hand --><html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Improve R’s Performance using JuliaCall with Mandel Set Example • JuliaCall</title>
<!-- jquery --><script src="https://code.jquery.com/jquery-3.1.0.min.js" integrity="sha384-nrOSfDHtoPMzJHjVTdCopGqIqeYETSXhZDFyniQ8ZHcVy08QesyHcnOUpMpqnmWq" crossorigin="anonymous"></script><!-- Bootstrap --><link href="https://maxcdn.bootstrapcdn.com/bootswatch/3.3.7/cerulean/bootstrap.min.css" rel="stylesheet" crossorigin="anonymous">
<script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js" integrity="sha384-Tc5IQib027qvyjSMfHjOMaLkfuWVxZxUPnCJA7l2mCWNIpG9mGCD8wGNIcPD7Txa" crossorigin="anonymous"></script><!-- Font Awesome icons --><link href="https://maxcdn.bootstrapcdn.com/font-awesome/4.6.3/css/font-awesome.min.css" rel="stylesheet" integrity="sha384-T8Gy5hrqNKT+hzMclPo118YTQO6cYprQmhrYwIiQ/3axmI1hQomh7Ud2hPOy8SP1" crossorigin="anonymous">
<!-- pkgdown --><link href="../pkgdown.css" rel="stylesheet">
<script src="../jquery.sticky-kit.min.js"></script><script src="../pkgdown.js"></script><!-- mathjax --><script src="https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script><!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]-->
</head>
<body>
<div class="container template-vignette">
<header><div class="navbar navbar-default navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="../index.html">JuliaCall</a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
<li>
<a href="..//index.html">
<span class="fa fa-home fa-lg"></span>
</a>
</li>
<li>
<a href="../reference/index.html">Reference</a>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
Articles
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="../articles/Julia_in_RMarkdown.html">Julia in RMarkdown</a>
</li>
<li>
<a href="../articles/JuliaCall_in_Jupyter_R_Notebook.html">JuliaCall in Jupyter R Notebook</a>
</li>
<li>
<a href="../articles/mandelbrot.html">Improve R's Performance using JuliaCall with Mandel Set Example</a>
</li>
</ul>
</li>
<li>
<a href="../news/index.html">News</a>
</li>
</ul>
<ul class="nav navbar-nav navbar-right">
<li>
<a href="https://github.com/Non-Contradiction/JuliaCall">
<span class="fa fa-github fa-lg"></span>
</a>
</li>
</ul>
</div>
<!--/.nav-collapse -->
</div>
<!--/.container -->
</div>
<!--/.navbar -->
</header><div class="row">
<div class="col-md-9">
<div class="page-header toc-ignore">
<h1>Improve R’s Performance using JuliaCall with Mandel Set Example</h1>
</div>
<div class="contents">
<div id="what-is-mandelbrot-set" class="section level2">
<h2 class="hasAnchor">
<a href="#what-is-mandelbrot-set" class="anchor"></a>What is Mandelbrot Set?</h2>
<p>Mandelbrot set is the set of complex numbers <span class="math inline">\(c\)</span> for which the sequence <span class="math inline">\(\{0, f_{c}(0), f_{c}(f_{c}(0)), \ldots, f^{(n)}_{c}(0), \ldots\}\)</span> remains bounded in absolute value where function <span class="math inline">\(f_{c}(z)=z^{2}+c\)</span>. For more detailed introduction, you could go to <a href="https://en.wikipedia.org/wiki/Mandelbrot_set">Wikipedia page for Mandelbrot set</a>.</p>
<p>The calculation for Mandelbrot set involves heavy computation, and here we will use JuliaCall to speed up the computation by rewriting R function in Julia, and we will also see some easy and useful tips in writing higher performance Julia code. In this example, we will see that JuliaCall actually brings <strong>more than 100 times speedup</strong> of the calculation.</p>
</div>
<div id="pure-r-implementation" class="section level2">
<h2 class="hasAnchor">
<a href="#pure-r-implementation" class="anchor"></a>Pure R Implementation</h2>
<p>The R implementation for the calculation of Mandelbrot set is quite straightforward. Note that it could be proved that if absolute value of one item in the sequence is greater than 2, the sequence <span class="math inline">\(\{0, f_{c}(0), f_{c}(f_{c}(0)), \ldots, f^{(n)}_{c}(0), \ldots\}\)</span> will be divergent.</p>
<p>In the <code>mandelbrot</code> function, we keep the iteration number that the absolute value of the item is greater than 2 until the maximum iteration times.</p>
<p>And in the <code>mandelbrotImage</code> function, for two given sequences of <span class="math inline">\(x\)</span> and <span class="math inline">\(y\)</span>, we calculate the value of <code>mandelbrot</code> function divided by the maximum iteration times on every grid point.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">mandelbrot <-<span class="st"> </span><span class="cf">function</span>(c, <span class="dt">iterate_max =</span> <span class="dv">500</span>){
z <-<span class="st"> </span>0i
<span class="cf">for</span> (i <span class="cf">in</span> <span class="dv">1</span><span class="op">:</span>iterate_max) {
z <-<span class="st"> </span>z <span class="op">^</span><span class="st"> </span><span class="dv">2</span> <span class="op">+</span><span class="st"> </span>c
<span class="cf">if</span> (<span class="kw">abs</span>(z) <span class="op">></span><span class="st"> </span><span class="fl">2.0</span>) {
<span class="kw">return</span>(i)
}
}
iterate_max
}
mandelbrotImage <-<span class="st"> </span><span class="cf">function</span>(xs, ys, <span class="dt">iterate_max =</span> <span class="dv">500</span>){
<span class="kw">sapply</span>(ys, <span class="cf">function</span>(y) <span class="kw">sapply</span>(xs, <span class="cf">function</span>(x) <span class="kw">mandelbrot</span>(x <span class="op">+</span><span class="st"> </span>y <span class="op">*</span><span class="st"> </span>1i, <span class="dt">iterate_max =</span> iterate_max))) <span class="op">/</span><span class="st"> </span>iterate_max
}</code></pre></div>
</div>
<div id="julia-implementation-in-r-using-juliacall" class="section level2">
<h2 class="hasAnchor">
<a href="#julia-implementation-in-r-using-juliacall" class="anchor"></a>Julia Implementation in R using JuliaCall</h2>
<p>To use Julia from JuliaCall, we first need to do some necessary setup work:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(JuliaCall)
<span class="kw"><a href="../reference/julia_setup.html">julia_setup</a></span>()
<span class="co">#> Julia at location /Applications/Julia-0.6.app/Contents/Resources/julia/bin will be used.</span>
<span class="co">#> Julia version 0.6.1 found.</span>
<span class="co">#> Julia initiation...</span>
<span class="co">#> Finish Julia initiation.</span>
<span class="co">#> Loading setup script for JuliaCall...</span>
<span class="co">#> Finish loading setup script for JuliaCall.</span></code></pre></div>
<p>And then we could just define Julia functions using <code>julia_command</code>, note that the syntax is quite similar to R and easy to understand.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="../reference/julia_command.html">julia_command</a></span>(<span class="st">"</span>
<span class="st">function mandelbrot(c, iterate_max = 500)</span>
<span class="st"> z = 0.0im</span>
<span class="st"> for i in 1:iterate_max</span>
<span class="st"> z = z ^ 2 + c</span>
<span class="st"> if abs2(z) > 4.0</span>
<span class="st"> return(i)</span>
<span class="st"> end</span>
<span class="st"> end</span>
<span class="st"> iterate_max</span>
<span class="st">end"</span>)
<span class="co">#> mandelbrot (generic function with 2 methods)</span>
<span class="kw"><a href="../reference/julia_command.html">julia_command</a></span>(<span class="st">"</span>
<span class="st">function mandelbrotImage(xs, ys, iterate_max = 500)</span>
<span class="st"> z = zeros(Float64, length(xs), length(ys))</span>
<span class="st"> for i in 1:length(xs)</span>
<span class="st"> for j in 1:length(ys)</span>
<span class="st"> z[i, j] = mandelbrot(xs[i] + ys[j] * im, iterate_max) / iterate_max</span>
<span class="st"> end</span>
<span class="st"> end</span>
<span class="st"> z</span>
<span class="st">end"</span>)
<span class="co">#> mandelbrotImage (generic function with 2 methods)</span></code></pre></div>
</div>
<div id="performance-comparison-between-r-and-julia-implementation" class="section level2">
<h2 class="hasAnchor">
<a href="#performance-comparison-between-r-and-julia-implementation" class="anchor"></a>Performance Comparison Between R and Julia Implementation</h2>
<p>This is the setting we use to compare the performance between R and Julia implementations.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">iterate_max <-<span class="st"> </span>1000L
centerx <-<span class="st"> </span><span class="fl">0.37522</span> <span class="co">#0.3750001200618655</span>
centery <-<span class="st"> </span><span class="op">-</span><span class="fl">0.22</span> <span class="co">#-0.2166393884377127</span>
step <-<span class="st"> </span><span class="fl">0.000002</span>
size <-<span class="st"> </span><span class="dv">125</span>
xs <-<span class="st"> </span><span class="kw">seq</span>(<span class="op">-</span>step <span class="op">*</span><span class="st"> </span>size, step <span class="op">*</span><span class="st"> </span>size, step) <span class="op">+</span><span class="st"> </span>centerx
ys <-<span class="st"> </span><span class="kw">seq</span>(<span class="op">-</span>step <span class="op">*</span><span class="st"> </span>size, step <span class="op">*</span><span class="st"> </span>size, step) <span class="op">+</span><span class="st"> </span>centery</code></pre></div>
<div id="time-for-pure-r-implementation" class="section level3">
<h3 class="hasAnchor">
<a href="#time-for-pure-r-implementation" class="anchor"></a>Time for Pure R Implementation</h3>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">system.time</span>(zR <-<span class="st"> </span><span class="kw">mandelbrotImage</span>(xs, ys, iterate_max))
<span class="co">#> 用户 系统 流逝 </span>
<span class="co">#> 17.512 0.109 17.908</span></code></pre></div>
</div>
<div id="time-for-julia-implementation-using-juliacall" class="section level3">
<h3 class="hasAnchor">
<a href="#time-for-julia-implementation-using-juliacall" class="anchor"></a>Time for Julia Implementation using JuliaCall</h3>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># A little warm up for Julia implementation</span>
<span class="kw">invisible</span>(<span class="kw"><a href="../reference/call.html">julia_call</a></span>(<span class="st">"mandelbrotImage"</span>, xs, ys, 2L))
<span class="kw">system.time</span>(zJL <-<span class="st"> </span><span class="kw"><a href="../reference/call.html">julia_call</a></span>(<span class="st">"mandelbrotImage"</span>, xs, ys, iterate_max))
<span class="co">#> 用户 系统 流逝 </span>
<span class="co">#> 0.196 0.001 0.199</span></code></pre></div>
<p>We could see that JuliaCall brings <strong>a lot of times speedup</strong> of the calculation, actually, we could see more speedup with larger problem scale, like <strong>100 times speedup</strong> or even more. I won’t show the result here because I don’t want to wait minutes for this RMarkdown document to be knited.</p>
</div>
</div>
<div id="tips-for-julia-performance" class="section level2">
<h2 class="hasAnchor">
<a href="#tips-for-julia-performance" class="anchor"></a>Tips for Julia Performance</h2>
<div id="write-small-functions" class="section level3">
<h3 class="hasAnchor">
<a href="#write-small-functions" class="anchor"></a>Write Small Functions</h3>
<p>A general advice in writing Julia (as well as R) is that you should write small functions which target at doing one thing. For example, it is possible to write <code>mandelbrot</code> and <code>mandelbrotImage</code> function together, but it is not a good practice. And the function call is also very cheap in Julia.</p>
</div>
<div id="type-statibility" class="section level3">
<h3 class="hasAnchor">
<a href="#type-statibility" class="anchor"></a>Type Statibility</h3>
<p>If you want to write high performance Julia code, you should write type stability functions. It means the variable in the functions should be of only one type.</p>
<p>For example, if you change the first line of <code>mandelbrot</code> functions like this:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="../reference/julia_command.html">julia_command</a></span>(<span class="st">"</span>
<span class="st">function mandelbrot1(c, iterate_max = 500)</span>
<span class="st"> z = 0 ## instead of z = 0.0im in the original example</span>
<span class="st"> for i in 1:iterate_max</span>
<span class="st"> z = z ^ 2 + c</span>
<span class="st"> if abs2(z) > 4.0</span>
<span class="st"> return(i)</span>
<span class="st"> end</span>
<span class="st"> end</span>
<span class="st"> iterate_max</span>
<span class="st">end"</span>)
<span class="co">#> mandelbrot1 (generic function with 2 methods)</span>
<span class="kw"><a href="../reference/julia_command.html">julia_command</a></span>(<span class="st">"</span>
<span class="st">function mandelbrotImage1(xs, ys, iterate_max = 500)</span>
<span class="st"> z = zeros(Float64, length(xs), length(ys))</span>
<span class="st"> for i in 1:length(xs)</span>
<span class="st"> for j in 1:length(ys)</span>
<span class="st"> z[i, j] = mandelbrot1(xs[i] + ys[j] * im, iterate_max) / iterate_max</span>
<span class="st"> end</span>
<span class="st"> end</span>
<span class="st"> z</span>
<span class="st">end"</span>)
<span class="co">#> mandelbrotImage1 (generic function with 2 methods)</span></code></pre></div>
<p>And we do the timing again:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># A little warm up for Julia implementation</span>
<span class="kw">invisible</span>(<span class="kw"><a href="../reference/call.html">julia_call</a></span>(<span class="st">"mandelbrotImage1"</span>, xs, ys, 2L))
<span class="kw">system.time</span>(zJL <-<span class="st"> </span><span class="kw"><a href="../reference/call.html">julia_call</a></span>(<span class="st">"mandelbrotImage1"</span>, xs, ys, iterate_max))
<span class="co">#> 用户 系统 流逝 </span>
<span class="co">#> 3.197 0.029 3.271</span></code></pre></div>
<p>We could see the function becomes much slower, because in the <code>mandelbrot1</code> function, <code>z</code> is an integer at the beginning, but becomes a complex number in the iteration. We could use <code>@code_warntype</code> or <code>code_warntype</code> and other tools provided by Julia to check about this problem, see <a href="https://docs.julialang.org/en/stable/manual/performance-tips/" class="uri">https://docs.julialang.org/en/stable/manual/performance-tips/</a> for more information.</p>
</div>
</div>
<div id="conclusion" class="section level2">
<h2 class="hasAnchor">
<a href="#conclusion" class="anchor"></a>Conclusion</h2>
<p>For the end of this article, let us have a look at our Mandel set!</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">image</span>(xs, ys, zR, <span class="dt">col =</span> <span class="kw">topo.colors</span>(<span class="dv">12</span>))
<span class="kw">image</span>(xs, ys, zJL, <span class="dt">col =</span> <span class="kw">topo.colors</span>(<span class="dv">12</span>))</code></pre></div>
<p><img src="mandelbrot_files/figure-html/unnamed-chunk-10-1.png" width="672"></p>
</div>
</div>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<div id="tocnav">
<h2 class="hasAnchor">
<a href="#tocnav" class="anchor"></a>Contents</h2>
<ul class="nav nav-pills nav-stacked">
<li><a href="#what-is-mandelbrot-set">What is Mandelbrot Set?</a></li>
<li><a href="#pure-r-implementation">Pure R Implementation</a></li>
<li><a href="#julia-implementation-in-r-using-juliacall">Julia Implementation in R using JuliaCall</a></li>
<li><a href="#performance-comparison-between-r-and-julia-implementation">Performance Comparison Between R and Julia Implementation</a></li>
<li><a href="#tips-for-julia-performance">Tips for Julia Performance</a></li>
<li><a href="#conclusion">Conclusion</a></li>
</ul>
</div>
</div>
</div>
<footer><div class="copyright">
<p>Developed by Changcheng Li.</p>
</div>
<div class="pkgdown">
<p>Site built with <a href="http://hadley.github.io/pkgdown/">pkgdown</a>.</p>
</div>
</footer>
</div>
</body>
</html>