Finished implementations and runner

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")