|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +from argparse import ArgumentParser |
| 4 | +from distributed import Client, Future |
| 5 | +import numpy as np |
| 6 | +import os |
| 7 | +import sys |
| 8 | +import time |
| 9 | + |
| 10 | + |
| 11 | +def init_julia(re, im, n): |
| 12 | + '''Initialize the complex domain. |
| 13 | +
|
| 14 | + Positional arguments: |
| 15 | + re -- minimum and maximum real value as 2-tuple |
| 16 | + im -- minimum and maximum imaginary value as 2-tuple |
| 17 | + n -- number of real and imaginary points as 2-tuple |
| 18 | + ''' |
| 19 | + re_vals, im_vals = np.meshgrid( |
| 20 | + np.linspace(re[0], re[1], n[0]), |
| 21 | + np.linspace(im[0], im[1], n[1]) |
| 22 | + ) |
| 23 | + domain = re_vals + im_vals*1j |
| 24 | + return domain.flatten() |
| 25 | + |
| 26 | + |
| 27 | +def init_pyx(dask_worker): |
| 28 | + import pyximport |
| 29 | + pyximport.install() |
| 30 | + sys.path.insert(0, os.getcwd()) |
| 31 | +# sys.path.insert(0, '/scratch/leuven/301/vsc30140/julia_set/') |
| 32 | + from julia_cython import julia_set |
| 33 | + |
| 34 | + |
| 35 | +def init_omp_pyx(dask_worker): |
| 36 | + import pyximport |
| 37 | + pyximport.install() |
| 38 | + sys.path.insert(0, os.getcwd()) |
| 39 | +# sys.path.insert(0, '/scratch/leuven/301/vsc30140/julia_set/') |
| 40 | + from julia_cython_omp import julia_set |
| 41 | + |
| 42 | + |
| 43 | + |
| 44 | +if __name__ == '__main__': |
| 45 | + arg_parser = ArgumentParser(description='Compute julia set') |
| 46 | + arg_parser.add_argument('--re_min', type=float, default=-1.8, |
| 47 | + help='minimum real value') |
| 48 | + arg_parser.add_argument('--re_max', type=float, default=1.8, |
| 49 | + help='maximum real value') |
| 50 | + arg_parser.add_argument('--im_min', type=float, default=-1.8, |
| 51 | + help='minimum imaginary value') |
| 52 | + arg_parser.add_argument('--im_max', type=float, default=1.8, |
| 53 | + help='maximum imaginary value') |
| 54 | + arg_parser.add_argument('--max_norm', type=float, default=2.0, |
| 55 | + help='maximum complex norm for z') |
| 56 | + arg_parser.add_argument('--n_re', type=int, default=100, |
| 57 | + help='number of points on the real axis') |
| 58 | + arg_parser.add_argument('--n_im', type=int, default=100, |
| 59 | + help='number of points on the imaginary axis') |
| 60 | + arg_parser.add_argument('--max_iters', type=int, default=300, |
| 61 | + help='maximum number of iterations') |
| 62 | + arg_parser.add_argument('--implementation', default='python', |
| 63 | + choices=['python', 'cython', 'cython_omp'], |
| 64 | + help='implementation to use') |
| 65 | + arg_parser.add_argument('--partitions', type=int, default=100, |
| 66 | + help='number of partitions for dask workers') |
| 67 | + arg_parser.add_argument('--host', required=True, |
| 68 | + help='hostname of the dask scheduler') |
| 69 | + arg_parser.add_argument('--port', type=int, required=True, |
| 70 | + help='port of the dask scheduler') |
| 71 | + options = arg_parser.parse_args() |
| 72 | + client = Client(f'{options.host}:{options.port:d}') |
| 73 | + if options.implementation == 'python': |
| 74 | + from julia_python import julia_set |
| 75 | + elif options.implementation == 'cython': |
| 76 | + from julia_cython import julia_set |
| 77 | + client.register_worker_callbacks(init_pyx) |
| 78 | + elif options.implementation == 'cython_omp': |
| 79 | + from julia_cython_omp import julia_set |
| 80 | + client.register_worker_callbacks(init_omp_pyx) |
| 81 | + else: |
| 82 | + msg = '{0} version not implemented\n' |
| 83 | + sys.stderr.write(msg.format(options.implementation)) |
| 84 | + sys.exit(1) |
| 85 | + |
| 86 | + domain = init_julia( |
| 87 | + (options.re_min, options.re_max), |
| 88 | + (options.im_min, options.im_max), |
| 89 | + (options.n_re, options.n_im) |
| 90 | + ) |
| 91 | + domains = np.array_split(domain, options.partitions) |
| 92 | + iterations = np.array_split(np.zeros(options.n_re*options.n_im, |
| 93 | + dtype=np.int32), options.partitions) |
| 94 | + start_time = time.time() |
| 95 | + futures = client.map(julia_set, domains, iterations) |
| 96 | + results = client.gather(futures) |
| 97 | + end_time = time.time() |
| 98 | + print('compute time = {0:.6f} s'.format(end_time - start_time)) |
| 99 | + np.savetxt('julia.txt', np.concatenate(results).reshape(options.n_re, |
| 100 | + options.n_im)) |
0 commit comments