Julia Set#
In the notebook I present different way to accelerate python code. This is a modified version from Loic Gouarin
The test case is the computation of the Julia set wikipedia
import os, sys
if sys.platform == 'darwin':
os.environ['CC'] = 'gcc-10'
os.environ['CXX'] = 'g++-10'
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.rcParams['figure.figsize'] = (6,6)
Initialize Mesh#
nx, ny = 512, 512 # mesh
lim, maxit = 400, 2000 # limits
vmin, vmax = 0, 200
x = np.linspace(-1.6, 1.6, nx)
y = np.linspace(-1.6, 1.6, ny)
c = -0.772691322542185 + 0.124281466072787j
Pure Python#
def juliaset_python(x, y, c, lim, maxit):
"""
returns Julia set
"""
julia = np.zeros((x.size, y.size))
for i in range(x.size):
for j in range(y.size):
z = x[i] + 1j*y[j]
ite = 0
while abs(z) < lim and ite < maxit:
z = z**2 + c
ite += 1
julia[j, i] = ite
return julia
def plot_julia_set(julia):
plt.figure(figsize=(6,6))
plt.imshow(julia, cmap = cm.Greys, vmin=vmin, vmax=vmax)
plot_julia_set(juliaset_python(x, y, c, lim, maxit))
numba#
Numba will accelerate the pure python function just with
the decorator @jit. Numba does everything for you.
from numba import jit
@jit(nopython=True, parallel=True)
def juliaset_numba(x, y, c, lim, maxit):
julia = np.zeros((x.size, y.size))
lim2 = lim*lim
c = complex(c) # needed for numba
for j in range(y.size):
for i in range(x.size):
z = complex(x[i], y[j])
ite = 0
while (z.real*z.real + z.imag*z.imag) < lim2 and ite < maxit:
z = z*z + c
ite += 1
julia[j, i] = ite
return julia
plot_julia_set(juliaset_numba(x, y, c, lim, maxit))
PyJulia#
PyJulia is a python module to import Julia function in your Python session. You can also run Julia code in a middle of a Jupyter notebook with a Python kernel.
To use pyjulia you need to install Julia and install PyCall.jl and REPL
julia> using Pkg
julia> ENV["PYTHON"] = "/usr/local/bin/python3.7"
julia> Pkg.add("PyCall")
julia> Pkg.add("REPL")
julia> Pkg.build("PyCall")
print the value of sys.executable to know the python path. But the cell above could do the job.
import julia
julia.install()
from julia.api import Julia
jl = Julia(compiled_modules=False)
[ Info: Julia version info
Julia Version 1.10.4
Commit 48d4fd4843 (2024-06-04 10:41 UTC)
Build Info:
https://github.com/conda-forge/julia-feedstock
Note: This is an unofficial build, please report bugs to the project
responsible for this build and not to the Julia project unless you can
reproduce the issue using official builds available at https://julialang.org/downloads
Platform Info:
OS: Linux (x86_64-conda-linux-gnu)
Ubuntu 22.04.5 LTS
uname: Linux 6.5.0-1025-azure #26~22.04.1-Ubuntu SMP Thu Jul 11 22:33:04 UTC 2024 x86_64 x86_64
CPU: AMD EPYC 7763 64-Core Processor:
speed user nice sys idle irq
#1 2684 MHz 477 s 0 s 130 s 3470 s 0 s
#2 2670 MHz 670 s 0 s 151 s 3302 s 0 s
#3 3199 MHz 322 s 0 s 144 s 3476 s 0 s
#4 3145 MHz 292 s 0 s 139 s 3707 s 0 s
Memory: 15.606491088867188 GB (14381.24609375 MB free)
Uptime: 421.27 sec
Load Avg: 1.42 0.73 0.3
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_CONDAPKG_BACKEND_BACKUP =
JULIA_CONDAPKG_BACKEND = System
JULIA_LOAD_PATH = @:@python-fortran:@stdlib
JULIA_LOAD_PATH_BACKUP =
JULIA_DEPOT_PATH_BACKUP =
JULIA_PROJECT_BACKUP =
JULIA_CONDAPKG_EXE = /usr/share/miniconda/bin/conda
JULIA_SSL_CA_ROOTS_PATH = /usr/share/miniconda/envs/python-fortran/ssl/cacert.pem
JULIA_DEPOT_PATH = /usr/share/miniconda/envs/python-fortran/share/julia:
LD_GOLD = /usr/share/miniconda/envs/python-fortran/bin/x86_64-conda-linux-gnu-ld.gold
LD_LIBRARY_PATH = /opt/hostedtoolcache/Python/3.9.20/x64/lib
JULIA_PROJECT = @python-fortran
JULIA_SSL_CA_ROOTS_PATH_BACKUP =
JULIA_CONDAPKG_EXE_BACKUP =
SELENIUM_JAR_PATH = /usr/share/java/selenium-server.jar
JAVA_HOME_11_X64 = /usr/lib/jvm/temurin-11-jdk-amd64
PKG_CONFIG_PATH = /opt/hostedtoolcache/Python/3.9.20/x64/lib/pkgconfig
GITHUB_PATH = /home/runner/work/_temp/_runner_file_commands/add_path_2bd343cc-b8d3-4016-802a-2940f75c6760
JAVA_HOME = /usr/lib/jvm/temurin-11-jdk-amd64
GRADLE_HOME = /usr/share/gradle-8.10.2
XDG_CONFIG_HOME = /home/runner/.config
ANT_HOME = /usr/share/ant
JAVA_HOME_8_X64 = /usr/lib/jvm/temurin-8-jdk-amd64
HOMEBREW_CLEANUP_PERIODIC_FULL_DAYS = 3650
DEPLOYMENT_BASEPATH = /opt/runner
ANDROID_NDK_LATEST_HOME = /usr/local/lib/android/sdk/ndk/27.1.12297006
CXXFLAGS = -fvisibility-inlines-hidden -fmessage-length=0 -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/python-fortran/include
DEBUG_CXXFLAGS = -fvisibility-inlines-hidden -fmessage-length=0 -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-all -fno-plt -Og -g -Wall -Wextra -fvar-tracking-assignments -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/python-fortran/include
LDFLAGS = -Wl,-O2 -Wl,--sort-common -Wl,--as-needed -Wl,-z,relro -Wl,-z,now -Wl,--disable-new-dtags -Wl,--gc-sections -Wl,--allow-shlib-undefined -Wl,-rpath,/usr/share/miniconda/envs/python-fortran/lib -Wl,-rpath-link,/usr/share/miniconda/envs/python-fortran/lib -L/usr/share/miniconda/envs/python-fortran/lib
HOME = /home/runner
DEBUG_CFLAGS = -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-all -fno-plt -Og -g -Wall -Wextra -fvar-tracking-assignments -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/python-fortran/include
CMAKE_PREFIX_PATH = /usr/share/miniconda/envs/python-fortran:/usr/share/miniconda/envs/python-fortran/x86_64-conda-linux-gnu/sysroot/usr
CONDA_JL_HOME = /usr/share/miniconda/envs/python-fortran
CPPFLAGS = -DNDEBUG -D_FORTIFY_SOURCE=2 -O2 -isystem /usr/share/miniconda/envs/python-fortran/include
JAVA_HOME_21_X64 = /usr/lib/jvm/temurin-21-jdk-amd64
GITHUB_EVENT_PATH = /home/runner/work/_temp/_github_workflow/event.json
JAVA_HOME_17_X64 = /usr/lib/jvm/temurin-17-jdk-amd64
ANDROID_NDK_HOME = /usr/local/lib/android/sdk/ndk/27.1.12297006
CONDA_JL_HOME_BACKUP =
HOMEBREW_NO_AUTO_UPDATE = 1
ANDROID_HOME = /usr/local/lib/android/sdk
PIPX_HOME = /opt/pipx
DEBUG_CPPFLAGS = -D_DEBUG -D_FORTIFY_SOURCE=2 -Og -isystem /usr/share/miniconda/envs/python-fortran/include
LEIN_HOME = /usr/local/lib/lein
PATH = /usr/share/miniconda/envs/python-fortran/bin:/usr/share/miniconda/condabin:/usr/share/miniconda/condabin:/opt/hostedtoolcache/Python/3.9.20/x64/bin:/opt/hostedtoolcache/Python/3.9.20/x64:/snap/bin:/home/runner/.local/bin:/opt/pipx_bin:/home/runner/.cargo/bin:/home/runner/.config/composer/vendor/bin:/usr/local/.ghcup/bin:/home/runner/.dotnet/tools:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin:/home/runner/.dotnet/tools
CFLAGS = -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /usr/share/miniconda/envs/python-fortran/include
SWIFT_PATH = /usr/share/swift/usr/bin
TERM = xterm-color
[ Info: Julia executable: /usr/share/miniconda/envs/python-fortran/bin/julia
[ Info: Trying to import PyCall...
┌ Error: `import PyCall` failed
│ exception =
│ ArgumentError: Package PyCall [438e738f-606a-5dbb-bf0a-cddfbfd45ab0] is required but does not seem to be installed:
│ - Run `Pkg.instantiate()` to install all recorded dependencies.
│
│ Stacktrace:
│ [1] _require(pkg::Base.PkgId, env::Nothing)
│ @ Base ./loading.jl:1926
│ [2] __require_prelocked(uuidkey::Base.PkgId, env::Nothing)
│ @ Base ./loading.jl:1812
│ [3] #invoke_in_world#3
│ @ ./essentials.jl:926 [inlined]
│ [4] invoke_in_world
│ @ ./essentials.jl:923 [inlined]
│ [5] _require_prelocked
│ @ ./loading.jl:1803 [inlined]
│ [6] _require_prelocked
│ @ ./loading.jl:1802 [inlined]
│ [7] macro expansion
│ @ ./lock.jl:267 [inlined]
│ [8] require(uuidkey::Base.PkgId)
│ @ Base ./loading.jl:1797
│ [9] top-level scope
│ @ /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/julia/install.jl:36
│ [10] include(mod::Module, _path::String)
│ @ Base ./Base.jl:495
│ [11] exec_options(opts::Base.JLOptions)
│ @ Base ./client.jl:318
│ [12] _start()
│ @ Base ./client.jl:552
└ @ Main /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/julia/install.jl:38
[ Info: Installing PyCall...
Resolving package versions...
Updating `/usr/share/miniconda/envs/python-fortran/share/julia/environments/python-fortran/Project.toml`
[438e738f] + PyCall v1.96.4
Updating `/usr/share/miniconda/envs/python-fortran/share/julia/environments/python-fortran/Manifest.toml`
[8f4d0f93] + Conda v1.10.2
[682c06a0] + JSON v0.21.4
[1914dd2f] + MacroTools v0.5.13
[69de0a69] + Parsers v2.8.1
[aea7be01] + PrecompileTools v1.2.1
[21216c6a] + Preferences v1.4.3
[438e738f] + PyCall v1.96.4
[81def892] + VersionParsing v1.3.0
[0dad84c5] + ArgTools v1.1.1
[56f22d72] + Artifacts
[2a0f44e3] + Base64
[ade2ca70] + Dates
[f43a241f] + Downloads v1.6.0
[7b1f6079] + FileWatching
[b27032c2] + LibCURL v0.6.4
[8f399da3] + Libdl
[37e2e46d] + LinearAlgebra
[d6f4376e] + Markdown
[a63ad114] + Mmap
[ca575930] + NetworkOptions v1.2.0
[de0858da] + Printf
[9a3f8284] + Random
[ea8e919c] + SHA v0.7.0
[9e88b42a] + Serialization
[fa267f1f] + TOML v1.0.3
[cf7118a7] + UUIDs
[4ec0a83e] + Unicode
[e66e0078] + CompilerSupportLibraries_jll v1.1.1+0
[deac9b47] + LibCURL_jll v8.4.0+0
[29816b5a] + LibSSH2_jll v1.11.0+1
[c8ffd9c3] + MbedTLS_jll v2.28.2+1
[14a3606d] + MozillaCACerts_jll v2023.1.10
[4536629a] + OpenBLAS_jll v0.3.23+4
[83775a58] + Zlib_jll v1.2.13+1
[8e850b90] + libblastrampoline_jll v5.8.0+1
[8e850ede] + nghttp2_jll v1.52.0+1
Precompiling project...
✓ CompilerSupportLibraries_jll
✓ VersionParsing
✓ Preferences
✓ PrecompileTools
✓ MacroTools
✓ Parsers
✓ JSON
✓ Conda
✓ PyCall
9 dependencies successfully precompiled in 28 seconds. 1 already precompiled.
Precompiling PyCall...
Precompiling PyCall... DONE
PyCall is installed and built successfully.
%%file julia_set.jl
function escapetime(z, c, lim, maxit)
for n = 1:maxit
if abs(z) > lim
return n-1
end
z = z*z + c
end
return maxit
end
function juliaset_julia(x :: Vector{Float64}, y :: Vector{Float64},
c :: Complex, lim , maxit )
nx = length(x)
ny = length(y)
julia = zeros(Float64, (nx, ny))
Threads.@sync for i in eachindex(x)
Threads.@spawn for j in eachindex(y)
@inbounds z = x[i] + 1im * y[j]
@inbounds julia[j, i] = escapetime(z, c, lim, maxit)
end
end
return julia
end
Writing julia_set.jl
from julia import Main
juliaset_julia = Main.include("julia_set.jl")
plot_julia_set(juliaset_julia(x, y, c, lim, maxit))
Pythran#
Pythran is a Python-to-C++ translator
Add a comment line before your python function and it runs much faster.
Configuration#
~/.pythranrc file on macos (gcc is installed with hombrew and pythran with pip)
[compiler]
include_dirs=/usr/local/opt/openblas/include
library_dirs=/usr/local/opt/openblas/lib
blas=openblas
CXX=g++-10
CC=gcc-10
%load_ext pythran.magic
%%pythran -fopenmp
import numpy as np
#pythran export juliaset_pythran(float64[], float64[],complex, int, int)
def juliaset_pythran(x, y, c, lim, maxit):
"""
returns Julia set
"""
juliap = np.zeros((x.size, y.size), dtype=np.int32)
#omp parallel for private(z, ite)
for j in range(y.size):
for i in range(x.size):
z = x[i] + 1j*y[j]
ite = 0
while abs(z) < lim and ite < maxit:
z = z**2 + c
ite += 1
juliap[j, i] = ite
return juliap
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[15], line 1
----> 1 get_ipython().run_cell_magic('pythran', '-fopenmp', '\nimport numpy as np\n\n#pythran export juliaset_pythran(float64[], float64[],complex, int, int)\ndef juliaset_pythran(x, y, c, lim, maxit):\n """ \n returns Julia set\n """\n juliap = np.zeros((x.size, y.size), dtype=np.int32)\n\n #omp parallel for private(z, ite)\n for j in range(y.size):\n for i in range(x.size):\n z = x[i] + 1j*y[j]\n ite = 0\n while abs(z) < lim and ite < maxit:\n z = z**2 + c\n ite += 1\n juliap[j, i] = ite\n\n return juliap\n')
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/IPython/core/interactiveshell.py:2517, in InteractiveShell.run_cell_magic(self, magic_name, line, cell)
2515 with self.builtin_trap:
2516 args = (magic_arg_s, cell)
-> 2517 result = fn(*args, **kwargs)
2519 # The code below prevents the output from being displayed
2520 # when using magics with decorator @output_can_be_silenced
2521 # when the last Python token in the expression is a ';'.
2522 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/pythran/magic.py:67, in PythranMagics.pythran(self, line, cell)
65 m.update(cell.encode('utf-8'))
66 module_name = "pythranized_" + m.hexdigest()
---> 67 module_path = pythran.compile_pythrancode(module_name, cell, **kwargs)
68 loader = importlib.machinery.ExtensionFileLoader(module_name, module_path)
69 spec = importlib.machinery.ModuleSpec(name=module_name, loader=loader,
70 origin=module_path)
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/pythran/toolchain.py:465, in compile_pythrancode(module_name, pythrancode, specs, opts, cpponly, pyonly, output_file, module_dir, report_times, **kwargs)
463 # Compile to binary
464 try:
--> 465 output_file = compile_cxxcode(module_name,
466 str(module),
467 output_binary=output_file,
468 **kwargs)
469 except PythranCompileError:
470 logger.warning("Compilation error, "
471 "trying hard to find its origin...")
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/pythran/toolchain.py:402, in compile_cxxcode(module_name, cxxcode, output_binary, keep_temp, **kwargs)
400 # Get a temporary C++ file to compile
401 fdpath = _write_temp(cxxcode, '.cpp')
--> 402 output_binary = compile_cxxfile(module_name, fdpath,
403 output_binary, **kwargs)
404 if not keep_temp:
405 # remove tempfile
406 os.remove(fdpath)
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/pythran/toolchain.py:347, in compile_cxxfile(module_name, cxxfile, output_binary, **kwargs)
342 extension = PythranExtension(module_name,
343 [cxxfile],
344 **kwargs)
346 try:
--> 347 setup(name=module_name,
348 ext_modules=[extension],
349 cmdclass={"build_ext": PythranBuildExt},
350 # fake CLI call
351 script_name='setup.py',
352 script_args=['--verbose'
353 if logger.isEnabledFor(logging.INFO)
354 else '--quiet',
355 'build_ext',
356 '--build-lib', builddir,
357 '--build-temp', buildtmp]
358 )
359 except SystemExit as e:
360 raise PythranCompileError(str(e))
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/__init__.py:117, in setup(**attrs)
115 # Make sure we have any requirements needed to interpret 'attrs'.
116 _install_setup_requires(attrs)
--> 117 return distutils.core.setup(**attrs)
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/core.py:183, in setup(**attrs)
181 # And finally, run all the commands found on the command line.
182 if ok:
--> 183 return run_commands(dist)
185 return dist
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/core.py:199, in run_commands(dist)
192 """Given a Distribution object run all the commands,
193 raising ``SystemExit`` errors in the case of failure.
194
195 This function assumes that either ``sys.argv`` or ``dist.script_args``
196 is already set accordingly.
197 """
198 try:
--> 199 dist.run_commands()
200 except KeyboardInterrupt:
201 raise SystemExit("interrupted")
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/dist.py:954, in Distribution.run_commands(self)
949 """Run each command that was seen on the setup script command line.
950 Uses the list of commands found and cache of command objects
951 created by 'get_command_obj()'.
952 """
953 for cmd in self.commands:
--> 954 self.run_command(cmd)
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/dist.py:950, in Distribution.run_command(self, command)
946 self.set_defaults()
947 # Postpone defaults until all explicit configuration is considered
948 # (setup() args, config files, command line and plugins)
--> 950 super().run_command(command)
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/dist.py:973, in Distribution.run_command(self, command)
971 cmd_obj = self.get_command_obj(command)
972 cmd_obj.ensure_finalized()
--> 973 cmd_obj.run()
974 self.have_run[command] = True
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/command/build_ext.py:359, in build_ext.run(self)
356 self.compiler.set_link_objects(self.link_objects)
358 # Now actually compile and link everything.
--> 359 self.build_extensions()
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/command/build_ext.py:476, in build_ext.build_extensions(self)
474 self._build_extensions_parallel()
475 else:
--> 476 self._build_extensions_serial()
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/command/build_ext.py:502, in build_ext._build_extensions_serial(self)
500 for ext in self.extensions:
501 with self._filter_build_errors(ext):
--> 502 self.build_extension(ext)
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/pythran/dist.py:121, in PythranBuildExtMixIn.build_extension(self, ext)
118 self.compiler.compiler_so[i] = 'x86_64'
120 try:
--> 121 return super(PythranBuildExtMixIn, self).build_extension(ext)
122 finally:
123 # Revert compiler settings
124 for key in prev.keys():
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/setuptools/_distutils/command/build_ext.py:557, in build_ext.build_extension(self, ext)
554 for undef in ext.undef_macros:
555 macros.append((undef,))
--> 557 objects = self.compiler.compile(
558 sources,
559 output_dir=self.build_temp,
560 macros=macros,
561 include_dirs=ext.include_dirs,
562 debug=self.debug,
563 extra_postargs=extra_args,
564 depends=ext.depends,
565 )
567 # XXX outdated variable, kept here in case third-part code
568 # needs it.
569 self._built_objects = objects[:]
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/numpy/distutils/ccompiler.py:91, in replace_method.<locals>.<lambda>(self, *args, **kw)
89 def replace_method(klass, method_name, func):
90 # Py3k does not have unbound method anymore, MethodType does not work
---> 91 m = lambda self, *args, **kw: func(self, *args, **kw)
92 setattr(klass, method_name, m)
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/numpy/distutils/ccompiler.py:274, in CCompiler_compile(self, sources, output_dir, macros, include_dirs, debug, extra_preargs, extra_postargs, depends)
237 """
238 Compile one or more source files.
239
(...)
270
271 """
272 global _job_semaphore
--> 274 jobs = get_num_build_jobs()
276 # setup semaphore to not exceed number of compile jobs when parallelized at
277 # extension level (python >= 3.5)
278 with _global_lock:
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/numpy/distutils/misc_util.py:91, in get_num_build_jobs()
76 def get_num_build_jobs():
77 """
78 Get number of parallel build jobs set by the --parallel command line
79 argument of setup.py
(...)
89
90 """
---> 91 from numpy.distutils.core import get_distribution
92 try:
93 cpu_count = len(os.sched_getaffinity(0))
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/numpy/distutils/core.py:24
22 from numpy.distutils.extension import Extension # noqa: F401
23 from numpy.distutils.numpy_distribution import NumpyDistribution
---> 24 from numpy.distutils.command import config, config_compiler, \
25 build, build_py, build_ext, build_clib, build_src, build_scripts, \
26 sdist, install_data, install_headers, install, bdist_rpm, \
27 install_clib
28 from numpy.distutils.misc_util import is_sequence, is_string
30 numpy_cmdclass = {'build': build.build,
31 'build_src': build_src.build_src,
32 'build_scripts': build_scripts.build_scripts,
(...)
44 'bdist_rpm': bdist_rpm.bdist_rpm,
45 }
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/numpy/distutils/command/config.py:19
17 import distutils
18 from numpy.distutils.exec_command import filepath_from_subprocess_output
---> 19 from numpy.distutils.mingw32ccompiler import generate_manifest
20 from numpy.distutils.command.autodist import (check_gcc_function_attribute,
21 check_gcc_function_attribute_with_intrinsics,
22 check_gcc_variable_attribute,
(...)
25 check_restrict,
26 check_compiler_gcc)
28 LANG_EXT['f77'] = '.f'
File /usr/share/miniconda/envs/python-fortran/lib/python3.9/site-packages/numpy/distutils/mingw32ccompiler.py:27
25 import distutils.cygwinccompiler
26 from distutils.unixccompiler import UnixCCompiler
---> 27 from distutils.msvccompiler import get_build_version as get_build_msvc_version
28 from distutils.errors import UnknownFileError
29 from numpy.distutils.misc_util import (msvc_runtime_library,
30 msvc_runtime_version,
31 msvc_runtime_major,
32 get_build_architecture)
ModuleNotFoundError: No module named 'distutils.msvccompiler'
plot_julia_set(juliaset_pythran(x, y, c, lim, maxit))
Fortran#
You need numpy and fortran-magic
%load_ext fortranmagic
On my mac i need to change compilers install with homebrew
Complex computation inside the loop are avoided on purpose. It takes time even with fortran.
%%fortran
subroutine juliaset_fortran(x, y, c, lim, maxit, julia)
real(8), intent(in) :: x(:)
real(8), intent(in) :: y(:)
complex(8), intent(in) :: c
real(8), intent(in) :: lim
integer, intent(in) :: maxit
integer, intent(out) :: julia(size(x),size(y))
real(8) :: zr, zi, limsq, cr, ci, tmp
integer :: ite, nx, ny
nx = size(x)
ny = size(y)
limsq = lim * lim
cr = real(c)
ci = imag(c)
do i = 1, nx
do j = 1, ny
zr = x(i)
zi = y(j)
ite = 0
do while (zr*zr+zi*zi < limsq .and. ite < maxit)
tmp = zr*zr - zi*zi
zi = 2*zr*zi + ci
zr = tmp + cr
ite = ite + 1
end do
julia(j, i) = ite
end do
end do
end subroutine juliaset_fortran
plot_julia_set(juliaset_fortran(x, y, c, lim, maxit))
Fortran with OpenMP#
On Linux you don’t need to add -L/usr/local/lib
This is the same function from above with some openmp directives
%%fortran --f90flags "-fopenmp" --opt "-O3" --extra "-lgomp"
subroutine juliaset_fortran_omp(x, y, c, lim, maxit, julia)
real(8), intent(in) :: x(:)
real(8), intent(in) :: y(:)
complex(8), intent(in) :: c
real(8), intent(in) :: lim
integer, intent(in) :: maxit
integer, intent(out) :: julia(size(x),size(y))
real(8) :: zr, zi, limsq, cr, ci, tmp
integer :: ite, nx, ny
nx = size(x)
ny = size(y)
limsq = lim * lim
cr = real(c)
ci = imag(c)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP FIRSTPRIVATE(nx,ny,x,y,c,limsq,maxit,cr,ci) &
!$OMP PRIVATE(i,j,ite,zr,zi, tmp) &
!$OMP SHARED(julia)
!$OMP DO SCHEDULE(DYNAMIC)
do i = 1, nx
do j = 1, ny
zr = x(i)
zi = y(j)
ite = 0
do while (zr*zr+zi*zi < limsq .and. ite < maxit)
tmp = zr*zr - zi*zi
zi = 2*zr*zi + ci
zr = tmp + cr
ite = ite + 1
end do
julia(j, i) = ite
end do
end do
!$OMP END PARALLEL
end subroutine juliaset_fortran_omp
plot_julia_set(juliaset_fortran_omp(x, y, c, lim, maxit))
Numpy#
Vectorize version with numpy. You could get some overflow warnings
import itertools
def juliaset_numpy(x, y, c, lim, maxit):
julia = np.zeros((x.size, y.size), dtype=np.int32)
zx = x[np.newaxis, :]
zy = y[:, np.newaxis]
z = zx + zy*1j
for ite in itertools.count():
z = z**2 + c
mask = np.logical_not(julia) & (np.abs(z) >= lim)
julia[mask] = ite
if np.all(julia) or ite > maxit:
return julia
plot_julia_set(juliaset_numpy(x, y, c, lim, maxit))
Cython#
Install Cython
Use %%cython -a to highlight the Python-C and C-Python conversions.
Cython is not Python and not C, it is another language :-)
%load_ext cython
%%cython
import numpy as np
import cython
@cython.boundscheck(False)
@cython.wraparound(False)
def juliaset_cython(double [:] x, double [:] y, double complex c, double lim, int maxit):
cdef:
int [:, ::1] julia = np.zeros((x.size, y.size), dtype = np.int32)
double tmp, zr, zi, lim2 = lim*lim
double cr = c.real, ci = c.imag
int ite, i, j, nx=x.size, ny=y.size
for i in range(nx):
for j in range(ny):
zr = x[i]
zi = y[j]
ite = 0
while (zr*zr + zi*zi) < lim2 and ite < maxit:
zr, zi = zr*zr - zi*zi + cr, 2*zr*zi + ci
ite += 1
julia[j, i] = ite
return julia
plot_julia_set(juliaset_cython(x, y, c, lim, maxit))
As f2py we can use openmp with the Cython prange function
%%cython -f -c-fopenmp --link-args=-fopenmp
import numpy as np
import cython
from cython.parallel import prange
from libc.stdlib cimport malloc, free
@cython.boundscheck(False)
@cython.wraparound(False)
def juliaset_cython_omp(double [:] x, double [:] y, double complex c, double lim, int maxit):
cdef:
int [:, ::1] julia = np.zeros((x.size, y.size), dtype = np.int32)
double tmp, zr, zi, lim2 = lim*lim
double cr = c.real, ci = c.imag
int i, j, nx=x.size, ny=y.size
int *ite
for j in prange(ny, nogil=True, schedule='dynamic'):
ite = <int *> malloc(sizeof(int))
for i in range(nx):
zr = x[i]
zi = y[j]
ite[0] = 0
while (zr*zr + zi*zi) < lim2 and ite[0] < maxit:
zr, zi = zr*zr - zi*zi + cr, 2*zr*zi + ci
ite[0] += 1
julia[j, i] = ite[0]
free(ite)
return julia
plot_julia_set(juliaset_cython_omp(x, y, c, lim, maxit))
Set number of threads used for parallel functions#
%env OMP_NUM_THREADS=4
%env JULIA_NUM_THREADS=4
env: OMP_NUM_THREADS=4
env: JULIA_NUM_THREADS=4
import pandas as pd
from collections import defaultdict
from tqdm import tqdm_notebook as tqdm
results = defaultdict(list)
nx, ny = 1024, 1024 # increase mesh size
x = np.linspace(-1.6, 1.6, nx)
y = np.linspace(-1.6, 1.6, ny)
functions = [juliaset_numpy,
juliaset_fortran,
juliaset_fortran_omp,
juliaset_cython,
juliaset_cython_omp,
juliaset_numba,
juliaset_pythran,
juliaset_julia]
for f in tqdm(functions):
_ = %timeit -q -n 1 -o f(x, y, c, lim, maxit)
results['etime'] += [_.best]
results = pd.DataFrame(results, index=list(map(lambda f:f.__name__[9:],functions)))
results["speed_up"] = [results.etime["numpy"]/t for t in results.etime]
results.sort_values(by="speed_up",axis=0)
| etime | speed_up | |
|---|---|---|
| numpy | 3.744096 | 1.000000 |
| pythran | 0.059130 | 63.320112 |
| numba | 0.056067 | 66.778523 |
| fortran | 0.052708 | 71.035016 |
| cython | 0.052008 | 71.990972 |
| julia | 0.025038 | 149.536701 |
| cython_omp | 0.009314 | 401.970603 |
| fortran_omp | 0.007442 | 503.092445 |