multithread optimizations

Author: D-McDougall

figures/ap_plot.py

@@ -13,8 +13,7 @@
 traces = {
     "matexp": {},
     "sparse": {},
-    "approx32": {},
-    "approx64": {},
+    "approx": {},
     "accuracy": {},
 }
 
@@ -42,7 +41,7 @@
     "Backwards Euler Method",
     "Approximate Matrix Exponential Method vs Time Step",
     "Approximate Matrix Exponential Method vs Accuracy",]
-methods = ["matexp", "sparse", "approx32", "accuracy"]
+methods = ["matexp", "sparse", "approx", "accuracy"]
 
 t_min, t_max = (3, 3.6)
 
@@ -67,7 +66,7 @@
     title = titles[index]
     method = methods[index]
     # 
-    axes.text(t_min+.1, 15, chr(ord("A") + index), ha='left', va='top',
+    axes.text(t_min+.05, 25, chr(ord("A") + index), ha='left', va='top',
               fontsize="large", weight="bold")
     num_traces = len(traces[method])
     for trace_index, (value, (t, v)) in enumerate(traces[method].items()):
@@ -83,7 +82,10 @@
     # axes.grid(which="major", axis='both', linestyle='solid', linewidth=1)
     axes.set_xlabel("time (ms)")
     axes.set_xlim(xmin=t_min, xmax=t_max)
-    axes.legend()
+    if index == 3:
+        axes.legend(loc='center right')
+    else:
+        axes.legend(loc='best')
 
 fig.savefig("ap_demo.png", dpi=600, bbox_inches='tight')
 if not os.environ.get('NOSHOW', ''): plt.show()

figures/ap_run.sh

@@ -32,7 +32,7 @@ max_errors=(
 )
 
 for err in "${max_errors[@]}"; do
-    python ap_sim.py approx .025 $err
+    python ap_sim.py approx .001 $err
 done
 
 python ap_plot.py

figures/err_run.sh

@@ -7,7 +7,8 @@ exec 2>&1
 
 CELLS=10000
 
-SEED=$RANDOM
+# SEED=$RANDOM
+SEED=12345
 
 # python -c "import numpy; print('\n'.join(str(x) for x in numpy.geomspace(.001, 1, 10)))"
 

matexp/__init__.py

@@ -11,12 +11,6 @@
 
 # Written by David McDougall, 2022-2026
 
-# Disable automatic multithreading
-import os
-os.environ['OPENBLAS_NUM_THREADS'] = '1'
-os.environ['MKL_NUM_THREADS'] = '1'
-os.environ['OMP_NUM_THREADS'] = '1'
-
 from .approx import Approx1D, Approx2D, MatrixSamples
 from .codegen import Codegen
 from .inputs import LinearInput, LogarithmicInput
@@ -25,14 +19,17 @@
 from pathlib import Path
 import multiprocessing
 import numpy as np
+import dill
 import time
+import os
 import sys
 
 __all__ = ('main', 'LinearInput', 'LogarithmicInput')
 
 _num_threads = len(os.sched_getaffinity(0))
 _thread_pool = None
-def _initialize_thread_pool(verbose):
+_derivative  = None
+def _initialize_thread_pool(model, verbose):
     global _thread_pool
     if verbose: print("Worker pool:", _num_threads, 'processes')
     # Manually delete any leftover shared memory files from a previous run.
@@ -42,15 +39,27 @@ def _initialize_thread_pool(verbose):
     else:
         pass # todo
     multiprocessing.set_start_method('spawn')
-    _thread_pool = multiprocessing.Pool(_num_threads)
+    _thread_pool = multiprocessing.Pool(
+            _num_threads,
+            _initialize_worker_process,
+            (dill.dumps(model.derivative),)) # Send the derivative function to every worker.
     return _thread_pool
 
+def _initialize_worker_process(derivative_pickle):
+    # Recv the derivative function.
+    global _derivative
+    _derivative = dill.loads(derivative_pickle)
+    # Disable automatic multithreading
+    os.environ['OPENBLAS_NUM_THREADS'] = '1'
+    os.environ['MKL_NUM_THREADS'] = '1'
+    os.environ['OMP_NUM_THREADS'] = '1'
+
 def main(nmodl_filename, inputs, time_step, temperature,
          error, target,
          outfile=None, verbose=False):
-    _initialize_thread_pool(verbose >= 2)
     # Read and process the NMODL file.
     model = LTI_Model(nmodl_filename, inputs, time_step, temperature)
+    _initialize_thread_pool(model, verbose >= 2)
     if   model.num_inputs == 1: OptimizerClass = Optimize1D
     elif model.num_inputs == 2: OptimizerClass = Optimize2D
     else: raise NotImplementedError('too many inputs.')

matexp/approx.py

@@ -4,7 +4,6 @@
 from itertools import pairwise, repeat
 import math
 import numpy as np
-import scipy.stats
 
 class MatrixSamples:
     def __init__(self, model, verbose=False):

matexp/lti_model.py

@@ -1,7 +1,6 @@
 from .nmodl_compiler import NMODL_Compiler
 from multiprocessing.shared_memory import SharedMemory
 from itertools import pairwise, repeat
-import dill
 import numpy as np
 import scipy.linalg
 
@@ -48,8 +47,10 @@ def make_deriv_matrix(self, inputs):
             boundaries = [num_samples * i // num_chunks for i in range(num_chunks + 1)]
             input_slices = [slice(*pair) for pair in pairwise(boundaries)]
             # 
-            derivative = dill.dumps(self.derivative)
-            args = (repeat(derivative), repeat(self.num_inputs), repeat(self.num_states), repeat(num_samples), input_slices)
+            args = (repeat(self.num_inputs),
+                    repeat(self.num_states),
+                    repeat(num_samples),
+                    input_slices)
             for _ in _thread_pool.map(self._compute_deriv, zip(*args), chunksize=1): pass
             # for _ in map(self._compute_deriv, zip(*args)): pass
             return np.ndarray(deriv_shape, dtype=np.float64, buffer=deriv_sm.buf).copy(), deriv_sm
@@ -60,8 +61,8 @@ def make_deriv_matrix(self, inputs):
 
     @staticmethod
     def _compute_deriv(args):
-        derivative, num_inputs, num_states, num_samples, input_slice = args
-        derivative = dill.loads(derivative)
+        num_inputs, num_states, num_samples, input_slice = args
+        from . import _derivative
         inputs_shape = (num_inputs, num_samples)
         deriv_shape = (num_samples, num_states, num_states)
         inputs_sm = SharedMemory('matexp_deriv_inputs', False)
@@ -74,7 +75,7 @@ def _compute_deriv(args):
         for col in range(num_states):
             state.fill(0.)
             state[col, :] = 1.
-            deriv[input_slice, :, col] = np.transpose(derivative(*chunk_inputs, *state))
+            deriv[input_slice, :, col] = np.transpose(_derivative(*chunk_inputs, *state))
         inputs_sm.close()
         deriv_sm.close()