diff options
| -rw-r--r-- | sem4/hpp/miniproject/multiproc.py | 137 | ||||
| -rw-r--r-- | sem4/hpp/miniproject/naive.py | 30 | ||||
| -rw-r--r-- | sem4/hpp/miniproject/optimised.py | 66 | ||||
| -rw-r--r-- | sem4/hpp/miniproject/runner.py | 82 |
4 files changed, 188 insertions, 127 deletions
diff --git a/sem4/hpp/miniproject/multiproc.py b/sem4/hpp/miniproject/multiproc.py index be39326..96fe788 100644 --- a/sem4/hpp/miniproject/multiproc.py +++ b/sem4/hpp/miniproject/multiproc.py @@ -1,107 +1,88 @@ -#!/usr/bin/env python3 - import numpy as np -import matplotlib.pyplot as plt -import time import multiprocessing as mp -from multiprocessing import RawArray + +# Local import +from optimised import numpylota # c-mesh limits limitre = ( -2, 1 ) limitim = ( -1.5, 1.5 ) +def createmangelctx(pre, pim, T, l, savez): + """ + Calculates the context used when calculating mangelbrot -def worker(gridchunk, s, step, T, l): - global rs - # Preallocate z array - rschunk = rs[s : s + step] + :param pre: Number of real numbers in cmatrix + :param pim: Number of imaginary numbers in cmatrix + :param T: Mangelbrot threshold + :param l: Iterations + """ + sre = ( limitre[1] - limitre[0] ) / (pre-1) - z = np.zeros(rschunk.shape) - print(f"rschunk.shape: {rschunk.shape}, gridchunk.shape: {gridchunk.shape}, T: {T}, l: {l}") + ctx = (sre, pim, T, l, savez) + + return ctx - # Calculate ι for all complex numbers - for i in range(l): - # This will generate warnings for some of the values rising above T. - # Because these values are above T they are not used, thus the warnings - # can be ignored - z = z*z + gridchunk +def mangelstep(ctx, re): + """ + Calculates a single mangelbrot row. + + :param ctx: Context containing information about mangelbrot + :param re: Row to calculate + """ - # This will generate 1 in all the places - # where z < T and zeros elsewhere - below = (np.abs(z) < T) + # Unpack context + (sre, pim, T, l, savez) = ctx - # Add this to the result - # Because the ones that pass T are 0 - # they will stop counting. - # - # If a specific z never reaches >= T its value in rs will - # be l - rschunk += below + # Create c-mesh row + im = np.linspace(limitim[0], limitim[1], pim, dtype=complex) + np.add(1j * im, limitre[0] + sre * re, im) - np.divide(rschunk, l, out=rschunk) + # Calculate ι + rs = numpyiota(im, T, l, savez) + if savez: + return rs + else: + # It takes some time to unpack a list of tuples + # So if no savez do not return a tuple + return rs[0] -def mangel(pre, pim, T, l, workers): +def mangel(pre, pim, T, l, savez): """ - Calculate the mangelbrot image with multiple processes + Calculate the mangelbrot image (pre, pim) discribes the image size. Use T and l to tune the mangelbrot - - This will split the image in horizontal parts and distribute it - between the workers. - Because the result array is row major, data will be nicely together if - the workers work with rows not columns. - - Pre must be devisible by workers. - - The result is saved in rs. Sorry couln't get numpy references through to the process as arguments + This function uses the global variables limitre and limitim to determine + the c-mesh range. :param pre: Number of real numbers used :param pim: Number of imaginary numbers :param T: Mangelbrot threshold :param l: Iterations - :param workers: Number of workers. + :param savez: Return z as the second element of returned tuple """ - # Used to calculate c-mesh - re = np.linspace(limitre[0], limitre[1], pre) - im = np.linspace(limitim[0], limitim[1], pim) - - # Calculate c-mesh - grid = np.add.outer(re, 1j * im) - - # Calculate the partition variables - step = int(pre / workers) - - - # Loop chunks and start the workers - wl = [] - for s in range(0, pre, step): - gridchunk = grid[s : s + step] - p = mp.Process(target=worker, args=(gridchunk, s, step, T, l)) - wl.append(p) - p.start() - - # Wait for them to be done - for p in wl: - p.join() + ctx = createmangelctx(pre, pim, T, l, savez) - return rs + # Number of processes + procs = mp.cpu_count() -rs = np.full((500, 500), 0.5) + pool = mp.Pool(processes=procs) -start = time.time() -arr = mangel(500, 500, 2, 100, 1) -end = time.time() + args = ((ctx, i) for i in range(pre)) + result = pool.starmap_async(mangelstep, args) + pool.close() + pool.join() -plt.imshow(arr, cmap=plt.cm.hot, vmin=0, vmax=1) -plt.savefig("mult.png") + result = result.get() + if savez: + z = [None] * pre + rs = [None] * pre -print(f"Took {end - start} seconds") -""" - p = mp.Process(target=worker, args=(rs[s:s + step], grid[s:s + step], T, l)) - wl.append(p) - p.start() + # Unzip rs + for i, r in enumerate(result): + rs[i] = r[0] + z[i] = r[1] - # Wait for them to be done - for p in wl: - p.join() -""" + return (np.vstack(rs), np.vstack(z)) + else: + return (np.vstack(result), None) diff --git a/sem4/hpp/miniproject/naive.py b/sem4/hpp/miniproject/naive.py index 765878a..24e366d 100644 --- a/sem4/hpp/miniproject/naive.py +++ b/sem4/hpp/miniproject/naive.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python3 import numpy as np import matplotlib.pyplot as plt import time @@ -7,9 +6,10 @@ import time limitre = ( -2, 1 ) limitim = ( -1.5, 1.5 ) -def lota(c, T, l): +def iota(c, T, l): """ Implement the ι function used in mangelbrot + Also devides by l :param c: Complex number from the c-mesh :param T: Mangelbrot threshold @@ -22,12 +22,12 @@ def lota(c, T, l): # Check if we found or z if np.abs(z) > T: - return i + return (i / l, z) # If we did not find z, use l - return l + return (l / l, z) -def mangel(pre, pim, T, l): +def mangel(pre, pim, T, l, savez): """ Calculate the mangelbrot image (pre, pim) discribes the image size. Use T and l to tune the mangelbrot @@ -38,10 +38,12 @@ def mangel(pre, pim, T, l): :param pim: Number of imaginary numbers :param T: Mangelbrot threshold :param l: Iterations + :param savez: Return z as the second element of returned tuple """ # Preallocate result array rs = np.zeros((pre, pim)) + z = np.empty((pre, pim), dtype=complex) # Calculate scaling variables sre = ( limitre[1] - limitre[0] ) / (pre-1) @@ -54,17 +56,9 @@ def mangel(pre, pim, T, l): c = limitre[0] + limitim[0] * 1j + sre * re + 1j * sim * im # Calculate the ι - rs[re,im] = lota(c, T, l) / l + (rs[re,im], z[re, im]) = iota(c, T, l) - return rs - - -start = time.time() -arr = mangel(500, 500, 2, 100) -end = time.time() - -plt.imshow(arr, cmap=plt.cm.hot, vmin=0, vmax=1) -plt.savefig("test.png") -plt.savefig("test.pdf") - -print(f"Took {end - start} seconds") + if savez: + return (rs, z) + else: + return (rs, None) diff --git a/sem4/hpp/miniproject/optimised.py b/sem4/hpp/miniproject/optimised.py index 2d128c9..299d25e 100644 --- a/sem4/hpp/miniproject/optimised.py +++ b/sem4/hpp/miniproject/optimised.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python3 - import numpy as np import matplotlib.pyplot as plt import time @@ -8,33 +6,23 @@ import time limitre = ( -2, 1 ) limitim = ( -1.5, 1.5 ) -def mangel(pre, pim, T, l): +def numpyiota(grid, T, l, savez): """ - Calculate the mangelbrot image - (pre, pim) discribes the image size. Use T and l to tune the mangelbrot - This function uses the global variables limitre and limitim to determine - the c-mesh range. - - :param pre: Number of real numbers used - :param pim: Number of imaginary numbers - :param T: Mangelbrot threshold - :param l: Iterations + Calculates the ι using numpy arrays. + This works with both numpy.arrays and numpy.ndarrays + Also devides by l + + :param grid: c-mesh. Shape will be used to make result + :param savez: Return z as the second element of returned tuple """ # Preallocate result array and z array - rs = np.zeros((pre, pim)) - z = np.zeros((pre, pim)) - - # Used to calculate c-mesh - re = np.linspace(limitre[0], limitre[1], pre) - im = np.linspace(limitim[0], limitim[1], pim) - - # Calculate c-mesh - grid = np.add.outer(re, 1j * im) + rs = np.zeros(grid.shape) + z = np.zeros(grid.shape) # Calculate ι for all complex numbers for i in range(l): # This will generate warnings for some of the values rising above T. - # Because these values are above T they are not used, thus the warnings + # Because these values become NAN they are not used, thus the warnings # can be ignored z = z*z + grid @@ -51,15 +39,31 @@ def mangel(pre, pim, T, l): rs += below rs /= l - - return rs -if __name__ == "__main__": - start = time.time() - arr = mangel(500, 500, 2, 100) - end = time.time() + if not savez: + z = None + return (rs, z) - plt.imshow(arr, cmap=plt.cm.hot, vmin=0, vmax=1) - plt.savefig("opt.png") +def mangel(pre, pim, T, l, savez): + """ + Calculate the mangelbrot image + (pre, pim) discribes the image size. Use T and l to tune the mangelbrot + This function uses the global variables limitre and limitim to determine + the c-mesh range. + + :param pre: Number of real numbers used + :param pim: Number of imaginary numbers + :param T: Mangelbrot threshold + :param l: Iterations + :param savez: Return z as the second element of returned tuple + """ + + # Used to calculate c-mesh + re = np.linspace(limitre[0], limitre[1], pre) + im = np.linspace(limitim[0], limitim[1], pim) + + # Calculate c-mesh + grid = np.add.outer(re, 1j * im) - print(f"Took {end - start} seconds") + # Calculate ι + return numpyiota(grid, T, l, savez) diff --git a/sem4/hpp/miniproject/runner.py b/sem4/hpp/miniproject/runner.py new file mode 100644 index 0000000..1565bd5 --- /dev/null +++ b/sem4/hpp/miniproject/runner.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 + +import matplotlib.pyplot as plt +import argparse +import os.path as path +import numpy as np +import time +from pathlib import Path +from PIL import Image +import warnings +import h5py + +# Many of the implementations overflow some of the values +# This creates warnings which can be ignored +warnings.simplefilter("ignore") + +# Implementations +impls = ["naive", "optimised", "multiproc"] + +parser = argparse.ArgumentParser(description="runs the different mangelbrot implementations") +parser.add_argument("impl", help="which implementation", default='all', choices=impls + ["all"], nargs='*') +parser.add_argument("--numreal", "--nr", type=int, help="number of real numbers", default=500) +parser.add_argument("--numimag", "--ni", type=int, help="number of imaginary numbers", default=500) +parser.add_argument("--iterations", "-i", type=int, help="mangelbrot iterations", default=100) +parser.add_argument("--threshold", "-t", type=int, help="mangelbrot threshold", default=2) +parser.add_argument("--outputfolder", "-o", help="output folder", default='.') +parser.add_argument("--createpng", help="create png images", action="store_true") +parser.add_argument("--skipz", "-z", help="do not export z", action="store_true") + +args = parser.parse_args() + +# If no implementations are specified run all +if args.impl == "all" or "all" in args.impl: + args.impl = impls + +def runtest(mangel, name): + """ + Run test on specified mangel function. + `mangel` function must be of form + + ``` + mangel(pre, pim, T, l) + ``` + + Uses the global args for options. + """ + + # Run mangel and measure runtime + start = time.time() + (rs, z) = mangel(args.numreal, args.numimag, args.threshold, args.iterations, not args.skipz) + end = time.time() + + plt.imshow(rs, cmap=plt.cm.hot, vmin=0, vmax=1) + plt.savefig(path.join(args.outputfolder, name + ".pdf")) + if args.createpng: + # Convert o simple b/w image + formatted = (rs * 255).astype('uint8') + + img = Image.fromarray(formatted, mode='L') + img.save(path.join(args.outputfolder, name + ".png")) + + # Save data as h5py + h5file = h5py.File(path.join(args.outputfolder, name + ".h5"), "w") + h5file.create_dataset("result", data=rs) + if not args.skipz: + h5file.create_dataset("z", data=z) + h5file.close() + + return end - start + +# Make sure output folder exists +Path(args.outputfolder).mkdir(parents=True, exist_ok=True) + +print(f"Running implementations {args.impl}") + +for impl in args.impl: + # Import module + m = __import__(impl) + + print(f"Running {impl}") + dur = runtest(m.mangel, impl) + print(f"{impl} took {dur} seconds") |