Repository: UC-Davis-molecular-computing/nuad (language: python)

Issue:
implement `SearchParameters.report_only_violations`
When implementing #92, we lost the ability to report constraint evaluations that passed, since only violations are stored in the `ViolationSet` that is saved between the constraint evaluation and reporting code.

One way to fix this is to rename the classes `Violation` and `ViolationSet` to `Evaluation` and `EvaluationSet`, so that each object in it represents a single evaluation of a `Constraint` on a design `Part`, whether the `Constraint` was violated or not.

`Evaluation` could essentially be the same class as `Violation`, but with an extra Boolean field `violated`. However, a significant number of parts of the code need to change to implement this, since currently it is assumed in many parts of code that only violations need to be stored.
implement parallelism for SingularConstraints
When addressing issue #93, we removed the support for parallel processing of `SingularConstraint`'s.

Re-implement it, and then address issue #132 if it still applies (i.e., if the parallelism still seems slower than it should be).
closes #213: display constraints report graphically in jupyter notebook



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 690fc630..e866a003 100644
--- a/README.md
+++ b/README.md
@@ -53,6 +53,8 @@ To install nuad, you can either install it using pip (the slightly simpler optio
 
 - git
 
+  This method has more steps, but it might be preferable if you want to use a new feature that is not on the main branch: the package installed by the pip instructions above will install the version currently on the main branch.
+
   1. Download the git repo, by one of two methods:
       - Install [git](https://git-scm.com/downloads) if necessary, then type 
     
@@ -88,30 +90,30 @@ Install NUPACK (version 4) and ViennaRNA following their installation instructio
 
 After installing ViennaRNA, it may be necessary to add its executables directory (the directory containing executable programs such as RNAduplex) to your `PATH` environment variable. (Similarly to how the `PYTHONPATH` variable is adjusted above.) NUPACK 4 does not come with an executable, so this step is unnecessary; it is called directly from within Python.
 
-To test that NUPACK 4 is installed correctly, run `python3 -m pip show nupack`.
+<!-- To test that NUPACK 4 is installed correctly, run `python3 -m pip show nupack`. -->
 
-To test that ViennaRNA is installed correctly, type `RNAduplex` at the command line.
+<!-- To test that ViennaRNA is installed correctly, type `RNAduplex` at the command line. -->
 
-Test NUPACK and ViennaRNA are available from within nuad by typing `python` at the command line, then typing `import nuad`. It should import without errors:
+Type `python` at the command line, then type `import nuad`. It should import without errors:
 
-    ```python
-    $ python
-    Python 3.7.5 (default, Nov  7 2019, 10:50:52)
-    [GCC 8.3.0] on linux
-    Type "help", "copyright", "credits" or "license" for more information.
-    >>> import nuad
-    >>>
-    ```
+```python
+$ python
+Python 3.7.5 (default, Nov  7 2019, 10:50:52)
+[GCC 8.3.0] on linux
+Type "help", "copyright", "credits" or "license" for more information.
+>>> import nuad
+>>>
+```
 
-To test that NUPACK and ViennaRNA can each be called from within the Python library (note that if you do not install NUPACK and/or ViennaRNA, then only a subset of the following will succeed):
+To test that NUPACK and ViennaRNA can each be called from within the Python library (note that if you do not install NUPACK and/or ViennaRNA, then this will fail):
 
-    ```python
-    >>> import nuad.vienna_nupack as nv
-    >>> nv.pfunc('GCGCGCGCGC')  # test NUPACK 4
-    -1.9079766874655928
-    >>> nv.rna_duplex_multiple([('GCGCGCGCGC', 'GCGCGCGCGC')]) # test ViennaRNA
-    [-15.7]
-    ```
+```python
+>>> import nuad.vienna_nupack as nv
+>>> nv.pfunc('GCGCGCGCGC')  # test NUPACK 4
+-1.9079766874655928
+>>> nv.rna_duplex_multiple([('GCGCGCGCGC', 'GCGCGCGCGC')]) # test ViennaRNA
+[-15.7]
+```
 
 ## Data model
 nuad allows one to go from a design with abstract "domains", such as `a`, `a*`, `b`, `b*`, to concrete DNA sequences, for example, 
diff --git a/examples/many_strands_no_common_domains.py b/examples/many_strands_no_common_domains.py
index 2b980c0e..26095b7b 100644
--- a/examples/many_strands_no_common_domains.py
+++ b/examples/many_strands_no_common_domains.py
@@ -139,19 +139,15 @@ def main() -> None:
     strand_pairs_rna_duplex_constraint = nc.rna_duplex_strand_pairs_constraint(
         threshold=-1.0, temperature=52, short_description='StrandPairRNA', parallel=parallel)
 
-    strand_individual_ss_constraint = nc.nupack_strand_complex_free_energy_constraint(
+    strand_individual_ss_constraint = nc.nupack_strand_free_energy_constraint(
         threshold=-1.0, temperature=52, short_description='StrandSS', parallel=parallel)
 
-    strand_individual_ss_constraint2 = nc.nupack_strand_complex_free_energy_constraint(
-        threshold=-1.0, temperature=52, short_description='StrandSS2', parallel=parallel)
-
-    strand_pair_nupack_constraint = nc.nupack_strand_pairs_constraint(
+    strand_pair_nupack_constraint = nc.nupack_strand_pair_constraint(
         threshold=3.0, temperature=52, short_description='StrandPairNUPACK', parallel=parallel, weight=0.1)
 
     params = ns.SearchParameters(constraints=[
         # domain_nupack_ss_constraint,
-        # strand_individual_ss_constraint,
-        # strand_individual_ss_constraint2,
+        strand_individual_ss_constraint,
         strand_pairs_rna_duplex_constraint,
         # strand_pair_nupack_constraint,
         # domain_pair_nupack_constraint,
@@ -163,11 +159,12 @@ def main() -> None:
         restart=args.restart,
         random_seed=random_seed,
         max_iterations=None,
-        save_sequences_for_all_updates=True,
-        save_report_for_all_updates=True,
-        save_design_for_all_updates=True,
+        # save_sequences_for_all_updates=True,
+        # save_report_for_all_updates=True,
+        # save_design_for_all_updates=True,
         force_overwrite=True,
-        scrolling_output=False,
+        # log_time=True,
+        # scrolling_output=False,
         # report_only_violations=False,
     )
     ns.search_for_dna_sequences(design, params)
diff --git a/examples/sample_designer.py b/examples/sample_designer.py
index f1d432ae..dde60b22 100644
--- a/examples/sample_designer.py
+++ b/examples/sample_designer.py
@@ -121,7 +121,7 @@ def main() -> None:
         threshold=-7.0, temperature=52, short_description='StrandPairCompl',
         pairs=[(strand0, strand2), (strand0, strand3), (strand1, strand2), (strand1, strand3)]
     )
-    strand_individual_ss_constraint = nc.nupack_strand_complex_free_energy_constraint(
+    strand_individual_ss_constraint = nc.nupack_strand_free_energy_constraint(
         threshold=-0.0, temperature=52, short_description='StrandSS')
 
     params = ns.SearchParameters(
diff --git a/examples/sst_canvas.py b/examples/sst_canvas.py
index 5a579c84..44538999 100644
--- a/examples/sst_canvas.py
+++ b/examples/sst_canvas.py
@@ -206,7 +206,7 @@ class Thresholds:
 
 
 def create_constraints(design: nc.Design, thresholds: Thresholds) -> List[nc.Constraint]:
-    strand_individual_ss_constraint = nc.nupack_strand_complex_free_energy_constraint(
+    strand_individual_ss_constraint = nc.nupack_strand_free_energy_constraint(
         threshold=thresholds.tile_ss, temperature=thresholds.temperature, short_description='StrandSS')
 
     # This reduces the number of times we have to create these sets from quadratic to linear
diff --git a/notebooks/Untitled.ipynb b/notebooks/Untitled.ipynb
new file mode 100644
index 00000000..0f6b8fe2
--- /dev/null
+++ b/notebooks/Untitled.ipynb
@@ -0,0 +1,90 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "cf3567b8-b41b-4ce0-aa83-f8cbd4dd45b3",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(0.796078431372549, 0.3764705882352941, 0.08235294117647059)"
+      ]
+     },
+     "execution_count": 1,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(203/255, 96/255, 21/255)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "1beec0c0-7e6f-4da6-a325-a4db11554a83",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(0.7490196078431373, 0.3411764705882353, 0.0)"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(191/255, 87/255, 0/255)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "4ff43f83-c941-4a79-8b74-c97c70814645",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 1,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "t = (2,3)\n",
+    "s = tuple(t)\n",
+    "s is t"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/notebooks/nuad_parallel_time_trials.ipynb b/notebooks/nuad_parallel_time_trials.ipynb
new file mode 100644
index 00000000..d6107c17
--- /dev/null
+++ b/notebooks/nuad_parallel_time_trials.ipynb
@@ -0,0 +1,650 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "6fc8df87-3dea-4b7a-856d-144ab5d7d917",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "10\n"
+     ]
+    }
+   ],
+   "source": [
+    "from multiprocessing.pool import ThreadPool\n",
+    "from pathos.pools import ProcessPool\n",
+    "import psutil\n",
+    "\n",
+    "num_cores = psutil.cpu_count(logical=False)\n",
+    "# num_cores = psutil.cpu_count(logical=True)\n",
+    "print(num_cores)\n",
+    "thread_pool = ThreadPool(processes=num_cores)\n",
+    "# process_pool = ProcessPool(processes=num_cores)\n",
+    "process_pool = ProcessPool(nodes=num_cores)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "69e37ca8-911f-4186-ab76-e5299763ab95",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "using 10 processes in ProcessPool\n",
+      "averaging each run over 5 trials\n",
+      "################################################################################\n",
+      "#  seq length = 10\n",
+      "#  ************************************************************\n",
+      "#  * num seqs = 100\n",
+      "#  * ave seq time = 0.08 sec\n",
+      "#  * ave par time = 0.28 sec\n",
+      "#  * seq/par      = 0.30\n",
+      "#\n",
+      "\n",
+      "################################################################################\n",
+      "#  seq length = 100\n",
+      "#  ************************************************************\n",
+      "#  * num seqs = 100\n",
+      "#  * ave seq time = 0.25 sec\n",
+      "#  * ave par time = 0.08 sec\n",
+      "#  * seq/par      = 2.95\n",
+      "#\n",
+      "\n",
+      "################################################################################\n",
+      "#  seq length = 125\n",
+      "#  ************************************************************\n",
+      "#  * num seqs = 100\n",
+      "#  * ave seq time = 0.37 sec\n",
+      "#  * ave par time = 0.09 sec\n",
+      "#  * seq/par      = 3.91\n",
+      "#\n",
+      "\n",
+      "################################################################################\n",
+      "#  seq length = 150\n",
+      "#  ************************************************************\n",
+      "#  * num seqs = 100\n",
+      "#  * ave seq time = 0.64 sec\n",
+      "#  * ave par time = 0.13 sec\n",
+      "#  * seq/par      = 4.85\n",
+      "#\n",
+      "\n",
+      "################################################################################\n",
+      "#  seq length = 175\n",
+      "#  ************************************************************\n",
+      "#  * num seqs = 100\n",
+      "#  * ave seq time = 0.87 sec\n",
+      "#  * ave par time = 0.19 sec\n",
+      "#  * seq/par      = 4.66\n",
+      "#\n",
+      "\n",
+      "################################################################################\n",
+      "#  seq length = 200\n",
+      "#  ************************************************************\n",
+      "#  * num seqs = 100\n",
+      "#  * ave seq time = 1.21 sec\n",
+      "#  * ave par time = 0.27 sec\n",
+      "#  * seq/par      = 4.57\n",
+      "#\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# pfunc\n",
+    "import nuad.vienna_nupack as nv\n",
+    "from pprint import pprint\n",
+    "from timeit import timeit\n",
+    "from importlib import reload\n",
+    "reload(nv)\n",
+    "\n",
+    "nv.pfunc('ACGTACGTAGCTA') # call once so overhead of loading Model isn't being measured in any call below\n",
+    "\n",
+    "num_trials = 5\n",
+    "print(f'using {process_pool.ncpus} processes in ProcessPool')\n",
+    "print(f'averaging each run over {num_trials} trials')\n",
+    "# for seq_length in [10, 15, 20, 30, 40, 50, 75, 100, 125, 150]:\n",
+    "for seq_length in [10, 100, 125, 150, 175, 200]:\n",
+    "# for seq_length in [200]:\n",
+    "    print('#'*80)\n",
+    "    print(f'#  seq length = {seq_length}')\n",
+    "    #for num_seqs in [1, 2, 4, 8, 10, 15, 20, 30, 40, 50, 75, 100, 200, 300, 400, 500]: # for lower seq_length\n",
+    "    # for num_seqs in [1, 2, 4, 8, 10, 15, 20]:\n",
+    "    for num_seqs in [100]:\n",
+    "        parallel_seqs = [nv.random_dna_seq(seq_length) for _ in range(num_seqs)]\n",
+    "        sequential_seqs = [nv.random_dna_seq(seq_length) for _ in range(num_seqs)]\n",
+    "        print(f'#  {\"*\"*60}')\n",
+    "        print(f'#  * num seqs = {num_seqs}')\n",
+    "\n",
+    "        def pfunc_sequential():\n",
+    "            energies = tuple(nv.pfunc(seq) for seq in sequential_seqs)\n",
+    "            return energies\n",
+    "        \n",
+    "        time_s = timeit(pfunc_sequential, number=num_trials)\n",
+    "        time_s_ave = time_s / num_trials\n",
+    "        print(f'#  * ave seq time = {time_s_ave:.2f} sec')\n",
+    "\n",
+    "        time_p = timeit(lambda: nv.pfunc_parallel(process_pool, parallel_seqs), number=num_trials)\n",
+    "        time_p_ave = time_p / num_trials\n",
+    "        print(f'#  * ave par time = {time_p_ave:.2f} sec')\n",
+    "\n",
+    "        print(f'#  * seq/par      = {time_s_ave/time_p_ave:.2f}')\n",
+    "        print('#')\n",
+    "    print()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "id": "3c022c07-8f3f-4b51-b833-0a642e5042dd",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "using 10 processes in ProcessPool\n",
+      "averaging each run over 3 trials\n",
+      "################################################################################\n",
+      "#  seq length = 10\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 100\n",
+      "#  * ave seq time = 0.48 sec\n",
+      "#  * ave par time = 0.42 sec\n",
+      "#  * seq/par      = 1.15\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 200\n",
+      "#  * ave seq time = 0.40 sec\n",
+      "#  * ave par time = 0.32 sec\n",
+      "#  * seq/par      = 1.25\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 300\n",
+      "#  * ave seq time = 0.29 sec\n",
+      "#  * ave par time = 0.35 sec\n",
+      "#  * seq/par      = 0.84\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 400\n",
+      "#  * ave seq time = 0.38 sec\n",
+      "#  * ave par time = 0.36 sec\n",
+      "#  * seq/par      = 1.04\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 500\n",
+      "#  * ave seq time = 0.39 sec\n",
+      "#  * ave par time = 0.33 sec\n",
+      "#  * seq/par      = 1.19\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 1000\n",
+      "#  * ave seq time = 0.39 sec\n",
+      "#  * ave par time = 0.38 sec\n",
+      "#  * seq/par      = 1.02\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 2000\n",
+      "#  * ave seq time = 0.39 sec\n",
+      "#  * ave par time = 0.36 sec\n",
+      "#  * seq/par      = 1.09\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 5000\n",
+      "#  * ave seq time = 0.45 sec\n",
+      "#  * ave par time = 0.48 sec\n",
+      "#  * seq/par      = 0.93\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 10000\n",
+      "#  * ave seq time = 0.57 sec\n",
+      "#  * ave par time = 1.16 sec\n",
+      "#  * seq/par      = 0.49\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 20000\n",
+      "#  * ave seq time = 0.71 sec\n",
+      "#  * ave par time = 1.02 sec\n",
+      "#  * seq/par      = 0.69\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 50000\n",
+      "#  * ave seq time = 1.17 sec\n",
+      "#  * ave par time = 1.04 sec\n",
+      "#  * seq/par      = 1.13\n",
+      "#\n",
+      "\n",
+      "################################################################################\n",
+      "#  seq length = 15\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 100\n",
+      "#  * ave seq time = 0.38 sec\n",
+      "#  * ave par time = 0.42 sec\n",
+      "#  * seq/par      = 0.91\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 200\n",
+      "#  * ave seq time = 0.40 sec\n",
+      "#  * ave par time = 0.38 sec\n",
+      "#  * seq/par      = 1.04\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 300\n",
+      "#  * ave seq time = 0.34 sec\n",
+      "#  * ave par time = 0.40 sec\n",
+      "#  * seq/par      = 0.85\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 400\n",
+      "#  * ave seq time = 0.40 sec\n",
+      "#  * ave par time = 0.38 sec\n",
+      "#  * seq/par      = 1.05\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 500\n",
+      "#  * ave seq time = 0.39 sec\n",
+      "#  * ave par time = 0.40 sec\n",
+      "#  * seq/par      = 0.98\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 1000\n",
+      "#  * ave seq time = 0.47 sec\n",
+      "#  * ave par time = 0.32 sec\n",
+      "#  * seq/par      = 1.47\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 2000\n",
+      "#  * ave seq time = 0.49 sec\n",
+      "#  * ave par time = 0.43 sec\n",
+      "#  * seq/par      = 1.15\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 5000\n",
+      "#  * ave seq time = 0.57 sec\n",
+      "#  * ave par time = 0.69 sec\n",
+      "#  * seq/par      = 0.82\n",
+      "#\n",
+      "#  ************************************************************\n",
+      "#  * num_pairs = 10000\n",
+      "#  * ave seq time = 0.62 sec\n"
+     ]
+    },
+    {
+     "ename": "KeyboardInterrupt",
+     "evalue": "",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
+      "Input \u001b[0;32mIn [19]\u001b[0m, in \u001b[0;36m<cell line: 11>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     20\u001b[0m time_s_ave \u001b[38;5;241m=\u001b[39m time_s \u001b[38;5;241m/\u001b[39m num_trials\n\u001b[1;32m     21\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m#  * ave seq time = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mtime_s_ave\u001b[38;5;132;01m:\u001b[39;00m\u001b[38;5;124m.2f\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m sec\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m---> 23\u001b[0m time_p \u001b[38;5;241m=\u001b[39m \u001b[43mtimeit\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43;01mlambda\u001b[39;49;00m\u001b[43m:\u001b[49m\u001b[43m \u001b[49m\u001b[43mnv\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrna_duplex_multiple_parallel\u001b[49m\u001b[43m(\u001b[49m\u001b[43mthread_pool\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpairs\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnumber\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnum_trials\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     24\u001b[0m time_p_ave \u001b[38;5;241m=\u001b[39m time_p \u001b[38;5;241m/\u001b[39m num_trials\n\u001b[1;32m     25\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m#  * ave par time = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mtime_p_ave\u001b[38;5;132;01m:\u001b[39;00m\u001b[38;5;124m.2f\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m sec\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
+      "File \u001b[0;32m/usr/lib/python3.8/timeit.py:233\u001b[0m, in \u001b[0;36mtimeit\u001b[0;34m(stmt, setup, timer, number, globals)\u001b[0m\n\u001b[1;32m    230\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mtimeit\u001b[39m(stmt\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpass\u001b[39m\u001b[38;5;124m\"\u001b[39m, setup\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpass\u001b[39m\u001b[38;5;124m\"\u001b[39m, timer\u001b[38;5;241m=\u001b[39mdefault_timer,\n\u001b[1;32m    231\u001b[0m            number\u001b[38;5;241m=\u001b[39mdefault_number, \u001b[38;5;28mglobals\u001b[39m\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m    232\u001b[0m     \u001b[38;5;124;03m\"\"\"Convenience function to create Timer object and call timeit method.\"\"\"\u001b[39;00m\n\u001b[0;32m--> 233\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mTimer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstmt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msetup\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtimer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mglobals\u001b[39;49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtimeit\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnumber\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m/usr/lib/python3.8/timeit.py:177\u001b[0m, in \u001b[0;36mTimer.timeit\u001b[0;34m(self, number)\u001b[0m\n\u001b[1;32m    175\u001b[0m gc\u001b[38;5;241m.\u001b[39mdisable()\n\u001b[1;32m    176\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 177\u001b[0m     timing \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43minner\u001b[49m\u001b[43m(\u001b[49m\u001b[43mit\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtimer\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    178\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n\u001b[1;32m    179\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m gcold:\n",
+      "File \u001b[0;32m<timeit-src>:6\u001b[0m, in \u001b[0;36minner\u001b[0;34m(_it, _timer, _stmt)\u001b[0m\n",
+      "Input \u001b[0;32mIn [19]\u001b[0m, in \u001b[0;36m<lambda>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     20\u001b[0m time_s_ave \u001b[38;5;241m=\u001b[39m time_s \u001b[38;5;241m/\u001b[39m num_trials\n\u001b[1;32m     21\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m#  * ave seq time = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mtime_s_ave\u001b[38;5;132;01m:\u001b[39;00m\u001b[38;5;124m.2f\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m sec\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m---> 23\u001b[0m time_p \u001b[38;5;241m=\u001b[39m timeit(\u001b[38;5;28;01mlambda\u001b[39;00m: \u001b[43mnv\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrna_duplex_multiple_parallel\u001b[49m\u001b[43m(\u001b[49m\u001b[43mthread_pool\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpairs\u001b[49m\u001b[43m)\u001b[49m, number\u001b[38;5;241m=\u001b[39mnum_trials)\n\u001b[1;32m     24\u001b[0m time_p_ave \u001b[38;5;241m=\u001b[39m time_p \u001b[38;5;241m/\u001b[39m num_trials\n\u001b[1;32m     25\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m#  * ave par time = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mtime_p_ave\u001b[38;5;132;01m:\u001b[39;00m\u001b[38;5;124m.2f\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m sec\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
+      "File \u001b[0;32m/mnt/c/Dropbox/git/nuad/nuad/vienna_nupack.py:373\u001b[0m, in \u001b[0;36mrna_duplex_multiple_parallel\u001b[0;34m(thread_pool, pairs, logger, temperature, parameters_filename)\u001b[0m\n\u001b[1;32m    370\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m calculate_energies_sequential(pairs)\n\u001b[1;32m    372\u001b[0m lists_of_sequence_pairs \u001b[38;5;241m=\u001b[39m nc\u001b[38;5;241m.\u001b[39mchunker(pairs, num_chunks\u001b[38;5;241m=\u001b[39mnum_cores)\n\u001b[0;32m--> 373\u001b[0m lists_of_energies \u001b[38;5;241m=\u001b[39m \u001b[43mthread_pool\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmap\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcalculate_energies_sequential\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlists_of_sequence_pairs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    374\u001b[0m energies \u001b[38;5;241m=\u001b[39m nc\u001b[38;5;241m.\u001b[39mflatten(lists_of_energies)\n\u001b[1;32m    375\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mtuple\u001b[39m(energies)\n",
+      "File \u001b[0;32m/usr/lib/python3.8/multiprocessing/pool.py:364\u001b[0m, in \u001b[0;36mPool.map\u001b[0;34m(self, func, iterable, chunksize)\u001b[0m\n\u001b[1;32m    359\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmap\u001b[39m(\u001b[38;5;28mself\u001b[39m, func, iterable, chunksize\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m    360\u001b[0m     \u001b[38;5;124;03m'''\u001b[39;00m\n\u001b[1;32m    361\u001b[0m \u001b[38;5;124;03m    Apply `func` to each element in `iterable`, collecting the results\u001b[39;00m\n\u001b[1;32m    362\u001b[0m \u001b[38;5;124;03m    in a list that is returned.\u001b[39;00m\n\u001b[1;32m    363\u001b[0m \u001b[38;5;124;03m    '''\u001b[39;00m\n\u001b[0;32m--> 364\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_map_async\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43miterable\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmapstar\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mchunksize\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m/usr/lib/python3.8/multiprocessing/pool.py:765\u001b[0m, in \u001b[0;36mApplyResult.get\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    764\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget\u001b[39m(\u001b[38;5;28mself\u001b[39m, timeout\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m--> 765\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    766\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mready():\n\u001b[1;32m    767\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTimeoutError\u001b[39;00m\n",
+      "File \u001b[0;32m/usr/lib/python3.8/multiprocessing/pool.py:762\u001b[0m, in \u001b[0;36mApplyResult.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    761\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mwait\u001b[39m(\u001b[38;5;28mself\u001b[39m, timeout\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m--> 762\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_event\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m/usr/lib/python3.8/threading.py:558\u001b[0m, in \u001b[0;36mEvent.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    556\u001b[0m signaled \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_flag\n\u001b[1;32m    557\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m signaled:\n\u001b[0;32m--> 558\u001b[0m     signaled \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_cond\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    559\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m signaled\n",
+      "File \u001b[0;32m/usr/lib/python3.8/threading.py:302\u001b[0m, in \u001b[0;36mCondition.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    300\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:    \u001b[38;5;66;03m# restore state no matter what (e.g., KeyboardInterrupt)\u001b[39;00m\n\u001b[1;32m    301\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m timeout \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 302\u001b[0m         \u001b[43mwaiter\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43macquire\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    303\u001b[0m         gotit \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[1;32m    304\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n",
+      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Process ForkPoolWorker-33:\n",
+      "Process ForkPoolWorker-31:\n",
+      "Process ForkPoolWorker-32:\n",
+      "Process ForkPoolWorker-46:\n",
+      "Process ForkPoolWorker-43:\n",
+      "Process ForkPoolWorker-41:\n",
+      "Process ForkPoolWorker-42:\n",
+      "Process ForkPoolWorker-36:\n",
+      "Process ForkPoolWorker-38:\n",
+      "Process ForkPoolWorker-45:\n",
+      "Process ForkPoolWorker-35:\n",
+      "Process ForkPoolWorker-37:\n",
+      "Process ForkPoolWorker-39:\n",
+      "Process ForkPoolWorker-34:\n",
+      "Process ForkPoolWorker-44:\n",
+      "Traceback (most recent call last):\n",
+      "Process ForkPoolWorker-40:\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "Traceback (most recent call last):\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "Traceback (most recent call last):\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "Traceback (most recent call last):\n",
+      "Traceback (most recent call last):\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 315, in _bootstrap\n",
+      "    self.run()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/process.py\", line 108, in run\n",
+      "    self._target(*self._args, **self._kwargs)\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "KeyboardInterrupt\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/pool.py\", line 114, in worker\n",
+      "    task = get()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "KeyboardInterrupt\n",
+      "KeyboardInterrupt\n",
+      "KeyboardInterrupt\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "KeyboardInterrupt\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "KeyboardInterrupt\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "KeyboardInterrupt\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/queues.py\", line 358, in get\n",
+      "    with self._rlock:\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "KeyboardInterrupt\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "KeyboardInterrupt\n",
+      "KeyboardInterrupt\n",
+      "KeyboardInterrupt\n",
+      "  File \"/home/ddoty/.local/lib/python3.8/site-packages/multiprocess/synchronize.py\", line 101, in __enter__\n",
+      "    return self._semlock.__enter__()\n",
+      "KeyboardInterrupt\n",
+      "KeyboardInterrupt\n",
+      "KeyboardInterrupt\n",
+      "KeyboardInterrupt\n",
+      "KeyboardInterrupt\n"
+     ]
+    }
+   ],
+   "source": [
+    "# rna_duplex_multiple\n",
+    "import nuad.vienna_nupack as nv\n",
+    "from pprint import pprint\n",
+    "from timeit import timeit\n",
+    "from importlib import reload\n",
+    "reload(nv)\n",
+    "\n",
+    "num_trials = 3\n",
+    "print(f'using {thread_pool.??} threads in ThreadPool')\n",
+    "print(f'averaging each run over {num_trials} trials')\n",
+    "for seq_length in [10, 15, 20, 30, 40, 50, 75, 100]:\n",
+    "    print('#'*80)\n",
+    "    print(f'#  seq length = {seq_length}')\n",
+    "    for num_pairs in [100, 200, 300, 400, 500, 1000, 2000, 5000, 10000, 20000, 50000]:\n",
+    "        pairs = [(nv.random_dna_seq(seq_length), nv.random_dna_seq(seq_length)) for _ in range(num_pairs)]\n",
+    "        print(f'#  {\"*\"*60}')\n",
+    "        print(f'#  * num_pairs = {num_pairs}')\n",
+    "\n",
+    "        time_s = timeit(lambda: nv.rna_duplex_multiple(pairs), number=num_trials)\n",
+    "        time_s_ave = time_s / num_trials\n",
+    "        print(f'#  * ave seq time = {time_s_ave:.2f} sec')\n",
+    "\n",
+    "        time_p = timeit(lambda: nv.rna_duplex_multiple_parallel(thread_pool, pairs), number=num_trials)\n",
+    "        time_p_ave = time_p / num_trials\n",
+    "        print(f'#  * ave par time = {time_p_ave:.2f} sec')\n",
+    "\n",
+    "        print(f'#  * seq/par      = {time_s_ave/time_p_ave:.2f}')\n",
+    "        print('#')\n",
+    "    print()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 37,
+   "id": "cbb8a6fd-1890-4988-85bb-396fd094174a",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(-4.4, -4.3, -4.6, -2.0, -5.4, -4.1, -3.2, -6.1, -3.6, -4.9, -7.3, -5.9, -2.8, -4.5, -6.6, -9.4, -8.6, -12.4, -3.0, -3.6, -5.2, -9.8, -3.7, -4.6, -5.3, -6.5, -6.8, -4.7, -7.9, -6.9, -2.8, -6.1, -5.4, -3.3, -3.3, -5.2, -5.3, -4.7, -7.5, -4.3, -8.2, -4.1, -4.4, -5.2, -5.0, -5.6, -2.9, -3.1, -3.5, -7.5, -5.1, -9.8, -6.8, -7.9, -3.7, -10.0, -7.6, -3.2, -4.2, -4.0, -4.5, -8.1, -6.3, -6.5, -5.3, -3.3, -10.7, -5.3, -7.1, -4.2, -4.8, -2.0, -2.0, -3.3, -7.1, -3.3, -3.5, -5.4, -2.1, -5.4, -4.9, -5.6, -9.0, -7.8, -3.7, -5.7, -3.8, -6.0, -1.9, -5.5, -4.7, -4.7, -10.3, -4.6, -5.0, -6.2, -3.8, -7.3, -9.5, -3.2, -5.0, -3.6, -3.0, -2.6, -5.0, -7.4, -4.1, -7.9, -3.3, -5.4, -8.4, -6.2, -6.3, -3.8, -4.5, -3.8, -5.4, -5.9, -5.7, -3.5, -6.6, -6.4, -4.8, -2.1, -3.7, -3.8, -6.5, -3.2, -3.0, -5.0, -4.4, -5.2, -10.7, -4.0, -5.6, -6.6, -5.8, -4.6, -4.7, -5.6, -5.7, -3.1, -5.8, -5.5, -6.8, -6.5, -7.0, -6.2, -3.9, -4.6, -3.5, -2.9, -2.7, -4.9, -6.6, -8.6, -8.4, -4.7, -6.3, -11.6, -3.0, -2.7, -4.3, -4.2, -4.3, -5.6, -3.1, -6.6, -3.8, -3.9, -5.4, -3.0, -5.6, -4.5, -3.3, -3.9, -2.8, -3.9, -7.6, -8.4, -5.3, -7.4, -3.3, -7.7, -4.0, -4.8, -5.5, -5.8, -4.6, -6.1, -3.7, -7.4, -4.1, -2.4, -3.3, -3.2, -7.9, -4.2, -5.9, -7.0, -6.5, -3.6, -4.8, -7.7, -2.3, -5.9, -9.0, -6.9, -5.0, -1.7, -6.4, -2.3, -7.6, -8.3, -7.0, -4.9, -8.6, -6.8, -4.7, -2.9, -10.6, -9.8, -2.2, -8.3, -3.1, -5.2, -7.6, -3.4, -5.7, -7.1, -5.4, -6.5, -4.4, -3.2, -3.7, -3.7, -5.9, -5.6, -4.3, -4.3, -3.0, -8.3, -3.2, -5.7, -7.5, -3.7, -8.3, -8.6, -7.4, -7.2, -5.7, -8.0, -8.9, -4.2, -4.7, -2.1, -6.5, -6.5, -4.5, -2.3, -5.4, -5.4, -6.1, -2.2, -5.8, -4.2, -6.0, -7.1, -3.3, -4.8, -3.9, -5.6, -8.2, -4.1, -7.0, -2.4, -4.1, -3.6, -2.9, -3.5, -6.7, -4.5, -8.4, -7.1, -6.3, -6.0, -4.8, -3.1, -4.6, -2.5, -3.2, -1.9, -8.1, -7.7, -7.8, -2.5, -4.7, -5.9, -3.9, -2.7, -3.0, -2.8, -3.4, -6.8, -5.1, -5.8, -5.2, -7.1, -9.3, -5.6, -5.7, -3.1, -11.9, -4.4, -6.5, -4.6, -8.8, -2.1, -5.1, -3.4, -7.0, -4.4, -5.4, -4.9, -6.5, -2.9, -4.7, -6.2, -4.0, -4.7, -3.3, -5.9, -3.6, -6.6, -4.7, -6.4, -3.6, -1.7, -5.2, -2.6, -4.4, -6.5, -5.7, -6.1, -3.0, -3.9, -5.1, -3.1, -5.2, -6.3, -2.6, -3.8, -3.2, -2.9, -5.2, -7.9, -8.8, -3.1, -5.6, -7.4)\n"
+     ]
+    }
+   ],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "ef46573a-55b2-45d8-ba63-33321d42451f",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "ph.cpu_count()=8\n",
+      "ps.cpu_count()=8\n",
+      "psutil.cpu_count(logical=True)=8\n",
+      "psutil.cpu_count(logical=False)=4\n"
+     ]
+    },
+    {
+     "ename": "AttributeError",
+     "evalue": "'function' object has no attribute 'cache_clear'",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
+      "Input \u001b[0;32mIn [6]\u001b[0m, in \u001b[0;36m<cell line: 9>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mpsutil\u001b[38;5;241m.\u001b[39mcpu_count(logical\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\u001b[38;5;132;01m=}\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m      8\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mpsutil\u001b[38;5;241m.\u001b[39mcpu_count(logical\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\u001b[38;5;132;01m=}\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m----> 9\u001b[0m \u001b[43mpsutil\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcpu_count\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcache_clear\u001b[49m()\n",
+      "\u001b[0;31mAttributeError\u001b[0m: 'function' object has no attribute 'cache_clear'"
+     ]
+    }
+   ],
+   "source": [
+    "import pathos.helpers as ph\n",
+    "import psutil as psutil\n",
+    "from inspect import getsource\n",
+    "\n",
+    "print(f'{ph.cpu_count()=}')\n",
+    "print(f'{ps.cpu_count()=}')\n",
+    "print(f'{psutil.cpu_count(logical=True)=}')\n",
+    "print(f'{psutil.cpu_count(logical=False)=}')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "67c2455f-645d-4a75-bac1-654816e28f8d",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[[1, 2], [3, 4], [5, 6], [7, 8], [9]]"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "import nuad.constraints as nc\n",
+    "from importlib import reload\n",
+    "reload(nc)\n",
+    "\n",
+    "nc.chunker(('A', 'B', 'C', 'D', 'E', 'F', 'G'), num_chunks=4)\n",
+    "num = 9\n",
+    "nc.chunker(list(range(1,num+1)), num_chunks=8)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/nuad/constraints.py b/nuad/constraints.py
index c8c823a2..6cf24939 100644
--- a/nuad/constraints.py
+++ b/nuad/constraints.py
@@ -22,6 +22,7 @@
 import os
 import math
 import json
+from decimal import Decimal
 from typing import List, Set, Optional, Dict, Callable, Iterable, Tuple, Union, Collection, TypeVar, Any, \
     cast, Generic, DefaultDict, FrozenSet, Iterator, Sequence, Type
 from dataclasses import dataclass, field, InitVar
@@ -32,7 +33,9 @@
 from multiprocessing.pool import ThreadPool
 from numbers import Number
 from enum import Enum, auto, unique
+import functools
 
+import pint
 import numpy as np  # noqa
 from ordered_set import OrderedSet
 
@@ -43,6 +46,10 @@
 import nuad.modifications as nm
 from nuad.json_noindent_serializer import JSONSerializable, json_encode, NoIndent
 
+from pint import UnitRegistry
+
+ureg = UnitRegistry()
+
 # need typing_extensions package prior to Python 3.8 to get Protocol object
 try:
     from typing import Protocol
@@ -2113,7 +2120,7 @@ def domains_not_substrings_of_each_other_constraint(
 
     # def evaluate(s1: str, s2: str, domain1: Optional[Domain], domain2: Optional[Domain]) -> float:
     def evaluate(seqs: Tuple[str, ...],
-                 domains: Optional[Tuple[Domain, Domain]]) -> Tuple[float, str]:  # noqa
+                 domains: Optional[Tuple[Domain, Domain]]) -> Result:  # noqa
         s1, s2 = seqs
         if len(s1) > len(s2):
             s1, s2 = s2, s1
@@ -2135,7 +2142,8 @@ def evaluate(seqs: Tuple[str, ...],
                 else:
                     summary = msg
                 score += 1.0
-        return score, summary
+
+        return Result(excess=score, summary=summary)
 
     return DomainPairConstraint(description='domains not substrings of each other',
                                 short_description=short_description,
@@ -2252,14 +2260,16 @@ class Strand(Part, JSONSerializable, Generic[StrandLabel, DomainLabel]):
     """Hash value of _domain_names_concatenated; cached for efficiency."""
 
     idt: Optional[IDTFields] = None
-    """Fields used when ordering strands from the synthesis company IDT 
+    """
+    Fields used when ordering strands from the synthesis company IDT 
     (Integrated DNA Technologies, Coralville, IA). If present (i.e., not equal to :const:`None`)
     then the method :py:meth:`Design.write_idt_bulk_input_file` can be called to automatically
     generate an text file for ordering strands in test tubes: 
     https://www.idtdna.com/site/order/oligoentry,
     as can the method :py:meth:`Design.write_idt_plate_excel_file` for writing a Microsoft Excel 
     file that can be uploaded to IDT's website for describing DNA sequences to be ordered in 96-well
-    or 384-well plates."""
+    or 384-well plates.
+    """
 
     _name: Optional[str] = None
     """Optional name of strand."""
@@ -3147,6 +3157,34 @@ def to_json_serializable(self, suppress_indent: bool = True) -> Dict[str, Any]:
 
         return dct
 
+    def write_design_file(self, directory: str = '.', filename: Optional[str] = None,
+                          extension: str = 'json') -> None:
+        """
+        Write text file representing this :any:`Design`,
+        suitable for reading by the scadnano web interface,
+        with the output file having the same name as the running script but with ``.py`` changed to
+        :attr:`default_scadnano_file_extension`,
+        unless `filename` is explicitly specified.
+        For instance, if the script is named ``my_design.py``,
+        then the design will be written to ``my_design.json``.
+        If `extension` is specified (but `filename` is not), then the design will be written to
+        ``my_design.<extension>``
+
+        The string written is that returned by :meth:`Design.to_json`.
+
+        :param directory:
+            directory in which to put file (default: current working directory)
+        :param filename:
+            filename (default: name of script with ``.py`` replaced by
+            ``.sc``).
+            Mutually exclusive with `extension`
+        :param extension:
+            extension for filename (default: ``.sc``)
+            Mutually exclusive with `filename`
+        """
+        content = self.to_json()
+        sc.write_file_same_name_as_running_python_script(content, extension, directory, filename)
+
     @staticmethod
     def from_design_file(filename: str,
                          strand_label_decoder: Callable[[Any], StrandLabel] = lambda label: label,
@@ -3738,7 +3776,8 @@ def get_group_name_from_strand_label(sc_strand: Strand) -> Any:
             raise AssertionError(f'label does not have either an attribute or a dict key "{group_key}"')
 
     def assign_sequences_to_scadnano_design(self, sc_design: sc.Design[StrandLabel, DomainLabel],
-                                            ignored_strands: Iterable[Strand] = ()) -> None:
+                                            ignored_strands: Iterable[Strand] = (),
+                                            overwrite: bool = False) -> None:
         """
         Assigns sequences from this :any:`Design` into `sc_design`.
 
@@ -3752,6 +3791,9 @@ def assign_sequences_to_scadnano_design(self, sc_design: sc.Design[StrandLabel,
             a scadnano design.
         :param ignored_strands:
             strands in the scadnano design that are to be ignored by the sequence designer.
+        :param overwrite:
+            if True, overwrites existing sequences; otherwise gives an error if an existing sequence
+            disagrees with the newly assigned sequence
         """
 
         # filter out ignored strands
@@ -3770,9 +3812,9 @@ def assign_sequences_to_scadnano_design(self, sc_design: sc.Design[StrandLabel,
 
         for sc_strand in sc_strands_to_include:
             domain_names = [domain.name for domain in sc_strand.domains]
-            if sc_strand.dna_sequence is None:
+            if sc_strand.dna_sequence is None or overwrite:
                 assert None not in domain_names
-                self._assign_to_strand_with_no_sequence(sc_strand, sc_design)
+                self._assign_to_strand_with_without_checking_existing_sequence(sc_strand, sc_design)
             elif None not in domain_names:
                 self._assign_to_strand_with_partial_sequence(sc_strand, sc_design, sc_domain_name_tuples)
             else:
@@ -3857,9 +3899,11 @@ def assign_modifications_to_scadnano_design(self, sc_design: sc.Design[StrandLab
                                      f'{sc_strand.modifications_int[offset]}')
                 sc_strand.modifications_int[offset] = mod_int.to_scadnano_modification()
 
-    def _assign_to_strand_with_no_sequence(self,
-                                           sc_strand: sc.Strand[StrandLabel, DomainLabel],
-                                           sc_design: sc.Design[StrandLabel, DomainLabel]) -> None:
+    def _assign_to_strand_with_without_checking_existing_sequence(
+            self,
+            sc_strand: sc.Strand[StrandLabel, DomainLabel],
+            sc_design: sc.Design[StrandLabel, DomainLabel]
+    ) -> None:
         # check types
         if not isinstance(sc_design, sc.Design):
             raise TypeError(f'sc_design must be an instance of scadnano.Design, but it is {type(sc_design)}')
@@ -3881,8 +3925,7 @@ def _assign_to_strand_with_no_sequence(self,
             domain_sequence = dsd_domain.concrete_sequence(starred)
             sequence_list.append(domain_sequence)
         strand_sequence = ''.join(sequence_list)
-        sc_design.assign_dna(strand=sc_strand, sequence=strand_sequence, assign_complement=False,
-                             check_length=True)
+        sc_strand.set_dna_sequence(strand_sequence)
 
     @staticmethod
     def _assign_to_strand_with_partial_sequence(sc_strand: sc.Strand[StrandLabel, DomainLabel],
@@ -4115,10 +4158,143 @@ def _raise_unreachable():
     raise AssertionError('This should be unreachable')
 
 
+@dataclass
+class Result(Generic[DesignPart]):
+    """
+    A :any:`Result` is returned from the function :data:`SingularConstraint.evaluate`, and a list of
+    :any:`Result`'s is returned from the function :data:`BulkConstraint.evaluate_bulk`, describing the
+    result of evaluating the constraint on the design "part".
+
+    A :any:`Result` must have an "excess" and "summary" specified.
+
+    Optionally one may also specify a "value", which helps in graphically displaying the results of
+    evaluating constraints using the function :meth:`display_report`.
+
+    For example, if the constraint checks that the NUPACK complex free energy of a strand is at least
+    -2.5 kcal/mol, and a strand has energy -3.4 kcal/mol, then the following are sensible values for
+    these fields:
+
+    - ``value`` = ``-3.4``  or  ``"-3.4 kcal/mol"``  or  ``pint.Quantity(Decimal(-3.4), "kcal/mol")``
+    - ``excess`` = ``-0.9``
+    - ``summary`` = ``"-3.4 kcal/mol"``
+    """
+
+    excess: float
+    """
+    The excess is a nonnegative value that is turned into a score, and the search minimizes the total score 
+    of all constraint evaluations. Setting this to 0 (or a negative value) means the constraint 
+    is satisfied, and setting it to a positive value means the constraint is violated. The interpretation
+    is that the larger `excess` is, the more the constraint is violated.
+    
+    For example, a common value for excess is the amount by which the NUPACK complex free energy exceeds
+    a threshold.
+    """
+
+    summary: str = ''
+    """
+    This string is displayed in the text report on constraints, after the name of the "part" (e.g.,
+    strand, pair of domains, pair of strands).
+    """
+
+    value: Optional[pint.Quantity[Decimal]] = None
+    """
+    If this is a "numeric" constraint, i.e., checking some number such as the complex free energy of a 
+    strand and comparing it to a threshold, this is the "raw" value. It is optional, but if specified,
+    then the raw values can be plotted in a Jupyter notebook by the function :meth:`display_report`.
+    
+    If a ``float``, then no units are assumed. If it is a ``str``, then it is assumed that it can be 
+    passed to the constructor pint.Quantity and interpreted as a value with units, e.g., the string
+    "-3.4 kcal/mol".
+    """
+
+    score: float = field(init=False)
+    """
+    Set by the search algorithm based on :data:`Result.excess` as well as other data such as the 
+    constraint's weight and the :data:`SearchParameters.score_transfer_function`.
+    """
+
+    part: DesignPart = field(init=False)
+    """
+    Set by the search algorithm based on the part that was evaluated.
+    """
+
+    def __init__(self,
+                 excess: float,
+                 summary: Optional[str] = None,
+                 value: Optional[Union[float, str, pint.Quantity[Decimal]]] = None) -> None:
+        self.excess = excess
+        if summary is None:
+            if value is None:
+                raise ValueError('at least one of value or summary must be specified')
+            self.summary = str(value)
+        else:
+            self.summary = summary
+        if value is not None:
+            self.value = parse_and_normalize_quantity(value)
+
+        self.score = 0.0
+        self.part = None  # type:ignore
+
+
+def parse_and_normalize_quantity(quantity: Union[float, int, str, pint.Quantity]) \
+        -> pint.Quantity[Decimal]:
+    if isinstance(quantity, (str, float, int)):
+        quantity = ureg(quantity)
+    quantity = normalize_quantity(quantity)
+    return quantity
+
+
+def Q_(qty: int | str | Decimal | float, unit: str | pint.Unit) -> pint.Quantity[Decimal]:
+    # Convenient constructor for units, eg, :code:`Q_(5.0, 'nM')`.
+    # Ensures that the quantity is a Decimal.
+    if isinstance(qty, Decimal):
+        return ureg.Quantity(qty, unit)
+    else:
+        # we convert to string to avoid floating-point weirdness. For example
+        # ureg.Quantity(Decimal(-2.1), 'kcal/mol') gives
+        #   -2.100000000000000088817841970012523233890533447265625 kilocalorie / mole,
+        # but ureg.Quantity(Decimal(str(-2.1)), 'kcal/mol') gives
+        #   -2.1 kilocalorie / mole,
+        qty_str = str(qty)
+        return ureg.Quantity(Decimal(qty_str), unit)
+
+
+def normalize_quantity(quantity: pint.Quantity, compact: bool = False) -> pint.Quantity[Decimal]:
+    """
+    Normalize `quantity` so that it has a Decimal madnitude,
+    is "compact" if specified (uses units within the correct "3 orders of magnitude":
+    https://pint.readthedocs.io/en/0.18/tutorial.html#simplifying-units)
+    and eliminate trailing zeros.
+
+    :param quantity:
+        a pint Quantity[Decimal]
+    :param compact:
+        whether to change units to make compact (within correct 3 orders of magnitude, e.g.,
+        30 kg instead of 30,000 g)
+    :return:
+        `quantity` normalized to be compact and without trailing zeros.
+    """
+    if not isinstance(quantity.magnitude, Decimal):
+        quantity = Q_(quantity.magnitude, quantity.units)
+    if compact:
+        quantity = quantity.to_compact()
+    mag_int = quantity.magnitude.to_integral()
+    if mag_int == quantity.magnitude:
+        # can be represented exactly as integer, so return that;
+        # quantity.magnitude.normalize() would use scientific notation in this case, which we don't want
+        quantity = Q_(mag_int, quantity.units)
+    else:
+        # is not exact integer, so normalize will return normal float literal such as 10.2
+        # and not scientific notation like it would for an integer
+        mag_norm = quantity.magnitude.normalize()
+        quantity = Q_(mag_norm, quantity.units)
+    return quantity
+
+
 @dataclass(eq=False)
 class SingularConstraint(Constraint[DesignPart], Generic[DesignPart], ABC):
     evaluate: Callable[[Tuple[str, ...], Optional[DesignPart]],
-                       Tuple[float, str]] = lambda _: _raise_unreachable()
+                       Result[DesignPart]] = lambda _: _raise_unreachable()
     """
     Essentially a wrapper for a function that evaluates the :any:`Constraint`. 
     It takes as input a tuple of DNA sequences 
@@ -4126,11 +4302,14 @@ class SingularConstraint(Constraint[DesignPart], Generic[DesignPart], ABC):
     :any:`Domain`, :any:`Strand`, :any:`DomainPair`, :any:`StrandPair`, or :any:`Complex`
     (the latter being an alias for arbitrary-length tuple of :any:`Strand`'s).
 
-    The second argument will be None if :data:`Constraint.parallel` is True (since it's more expensive
-    to serialize the :any:`Domain` and :any:`Strand` objects than strings for passing data to processes
-    executing in parallel). Thus, if the :any:`Constraint` needs to use more data about the :any:`Part`
-    than just its DNA sequence, by accessing the second argument, :data:`Constraint.parallel` should be set
-    to False.
+    The second argument will be None if :data:`SingularConstraint.parallel` is True 
+    (since it's more expensive to serialize the :any:`Domain` and :any:`Strand` objects than strings for 
+    passing data to processes executing in parallel).
+     
+    Thus, if the :any:`Constraint` needs to use more data about the :any:`Part` than just its DNA sequence, 
+    by accessing the second argument, :data:`Constraint.parallel` should be set to False.
+    
+    It should return a :any:`Result` object.
     """
 
     parallel: bool = False
@@ -4151,7 +4330,7 @@ def __post_init__(self) -> None:
         if self.weight <= 0:
             raise ValueError(f'weight must be positive but it is {self.weight}')
 
-    def call_evaluate(self, seqs: Tuple[str, ...], part: Optional[DesignPart]) -> Tuple[float, str]:
+    def call_evaluate(self, seqs: Tuple[str, ...], part: Optional[DesignPart]) -> Result[DesignPart]:
         """
         Evaluates this :any:`Constraint` using function :data:`SingularConstraint.evaluate`
         supplied in constructor.
@@ -4171,27 +4350,28 @@ def call_evaluate(self, seqs: Tuple[str, ...], part: Optional[DesignPart]) -> Tu
             `excess` might be the difference 1.5 between the energy and the threshold,
             and `summary` might be the string "-2.5 kcal/mol".
         """
-        excess, summary = (self.evaluate)(seqs, part)  # noqa
-        if excess < 0.0:
-            excess = 0.0
-        score = self.weight * self.score_transfer_function(excess)
-        return score, summary
+        result = (self.evaluate)(seqs, part)  # noqa
+        if result.excess < 0.0:
+            result.excess = 0.0
+        result.score = self.weight * self.score_transfer_function(result.excess)
+        result.part = part
+        return result
 
 
 @dataclass(eq=False)
 class BulkConstraint(Constraint[DesignPart], Generic[DesignPart], ABC):
     evaluate_bulk: Callable[[Sequence[DesignPart]],
-                            List[Tuple[DesignPart, float, str]]] = lambda _: _raise_unreachable()
+                            List[Result]] = lambda _: _raise_unreachable()
 
-    def call_evaluate_bulk(self, parts: Sequence[DesignPart]) -> List[Tuple[DesignPart, float, str]]:
-        results: List[Tuple[DesignPart, float, str]] = (self.evaluate_bulk)(parts)  # noqa
-        results_with_scores_transferred_and_weighted: List[Tuple[DesignPart, float, str]] = []
-        for part, excess, summary in results:
-            if excess < 0.0:
-                excess = 0.0
-            score = self.weight * self.score_transfer_function(excess)
-            results_with_scores_transferred_and_weighted.append((part, score, summary))
-        return results_with_scores_transferred_and_weighted
+    def call_evaluate_bulk(self, parts: Sequence[DesignPart]) -> List[Result]:
+        results: List[Result[DesignPart]] = (self.evaluate_bulk)(parts)  # noqa
+        # apply weight and transfer scores
+        for result, part in zip(results, parts):
+            if result.excess < 0.0:
+                result.excess = 0.0
+            result.score = self.weight * self.score_transfer_function(result.excess)
+            result.part = part
+        return results
 
 
 _no_summary_string = "No summary for this constraint. " \
@@ -4424,15 +4604,15 @@ def __post_init__(self) -> None:
             raise ValueError('_evaluate_design should be specified in a DesignConstraint')
 
     def call_evaluate_design(self, design: Design, domains_changed: Iterable[Domain]) \
-            -> List[Tuple[DesignPart, float, str]]:
+            -> List[Result]:
         results = (self._evaluate_bulk)(design, domains_changed)  # noqa
-        results_with_scores_transferred_and_weighted = []
-        for part, excess, summary in results:
-            if excess < 0.0:
-                excess = 0.0
-            weighted_score = self.weight * self.score_transfer_function(excess)
-            results_with_scores_transferred_and_weighted.append((part, weighted_score, summary))
-        return results_with_scores_transferred_and_weighted
+        # apply weight and transfer scores
+        for result in zip(results):
+            if result.excess < 0.0:
+                result.excess = 0.0
+            result.score = self.weight * self.score_transfer_function(result.excess)
+            result.part = design
+        return results
 
     @staticmethod
     def part_name() -> str:
@@ -4562,11 +4742,12 @@ def nupack_domain_complex_free_energy_constraint(
     """
     _check_nupack_installed()
 
-    def evaluate(seqs: Tuple[str], _: Optional[Domain]) -> Tuple[float, str]:
+    def evaluate(seqs: Tuple[str], _: Optional[Domain]) -> Result:
         sequence = seqs[0]
-        energy = nv.complex_free_energy_single_strand(sequence, temperature, sodium, magnesium)
-        excess = threshold - energy
-        return max(0.0, excess), f'{energy:6.2f} kcal/mol'
+        energy = nv.free_energy_single_strand(sequence, temperature, sodium, magnesium)
+        excess = max(0.0, threshold - energy)
+        value = f'{energy:6.2f} kcal/mol'
+        return Result(excess=excess, value=value)
 
     if description is None:
         description = f'NUPACK secondary structure of domain exceeds {threshold} kcal/mol'
@@ -4583,7 +4764,7 @@ def evaluate(seqs: Tuple[str], _: Optional[Domain]) -> Tuple[float, str]:
                             domains=domains)
 
 
-def nupack_strand_complex_free_energy_constraint(
+def nupack_strand_free_energy_constraint(
         threshold: float,
         temperature: float = nv.default_temperature,
         sodium: float = nv.default_sodium,
@@ -4596,7 +4777,8 @@ def nupack_strand_complex_free_energy_constraint(
         strands: Optional[Iterable[Strand]] = None) -> StrandConstraint:
     """
     Returns constraint that checks individual :any:`Strand`'s for excessive interaction using
-    NUPACK's pfunc.
+    NUPACK's pfunc. This is the so-called "complex free energy":
+    https://docs.nupack.org/definitions/#complex-free-energy
 
     NUPACK 4 must be installed. Installation instructions can be found at
     https://piercelab-caltech.github.io/nupack-docs/start/.
@@ -4629,14 +4811,15 @@ def nupack_strand_complex_free_energy_constraint(
     """
     _check_nupack_installed()
 
-    def evaluate(seqs: Tuple[str], _: Optional[Strand]) -> Tuple[float, str]:
+    def evaluate(seqs: Tuple[str], _: Optional[Strand]) -> Result:
         sequence = seqs[0]
-        energy = nv.complex_free_energy_single_strand(sequence, temperature, sodium, magnesium)
-        excess = threshold - energy
-        return max(0.0, excess), f'{energy:6.2f} kcal/mol'
+        energy = nv.free_energy_single_strand(sequence, temperature, sodium, magnesium)
+        excess = max(0.0, threshold - energy)
+        value = f'{energy:6.2f} kcal/mol'
+        return Result(excess=excess, value=value)
 
     if description is None:
-        description = f'NUPACK secondary structure of strand exceeds {threshold} kcal/mol'
+        description = f'strand NUPACK energy >= {threshold} kcal/mol at {temperature}C'
 
     if strands is not None:
         strands = tuple(strands)
@@ -4713,7 +4896,7 @@ def binding_closure(seq_pair: Tuple[str, str]) -> float:
                           sodium=sodium, magnesium=magnesium)
 
     # def evaluate(seq1: str, seq2: str, domain1: Optional[Domain], domain2: Optional[Domain]) -> float:
-    def evaluate(seqs: Tuple[str, ...], domain_pair: Optional[DomainPair]) -> Tuple[float, str]:
+    def evaluate(seqs: Tuple[str, ...], domain_pair: Optional[DomainPair]) -> Result:
         seq1, seq2 = seqs
         name_pairs = [(None, None)] * 4
         if domain_pair is not None:
@@ -4749,16 +4932,17 @@ def evaluate(seqs: Tuple[str, ...], domain_pair: Optional[DomainPair]) -> Tuple[
 
         max_excess = max(excesses)
 
-        max_name_length = max(len(name) for name in _flatten(name_pairs))
+        max_name_length = max(len(name) for name in flatten(name_pairs))
         lines_and_energies = [(f'{name1:{max_name_length}}, '
                                f'{name2:{max_name_length}}: '
                                f' {energy:6.2f} kcal/mol', energy)
                               for (name1, name2), energy in zip(name_pairs, energies)]
         lines_and_energies.sort(key=lambda line_and_energy: line_and_energy[1])
         lines = [line for line, _ in lines_and_energies]
-        msg = '\n  ' + '\n  '.join(lines)
+        summary = '\n  ' + '\n  '.join(lines)
 
-        return max(0.0, max_excess), msg
+        max_excess = max(0.0, max_excess)
+        return Result(excess=max_excess, summary=summary, value=max_excess)
 
     if pairs is not None:
         pairs = tuple(pairs)
@@ -4772,7 +4956,93 @@ def evaluate(seqs: Tuple[str, ...], domain_pair: Optional[DomainPair]) -> Tuple[
                                 pairs=pairs)
 
 
-def nupack_strand_pairs_constraint(
+def nupack_strand_pair_constraint_by_number_matching_domains(
+        thresholds: Dict[int, float],
+        temperature: float = nv.default_temperature,
+        sodium: float = nv.default_sodium,
+        magnesium: float = nv.default_magnesium,
+        weight: float = 1.0,
+        score_transfer_function: Callable[[float], float] = default_score_transfer_function,
+        descriptions: Optional[Dict[int, str]] = None,
+        short_descriptions: Optional[Dict[int, str]] = None,
+        parallel: bool = False,
+        strands: Optional[Iterable[Strand]] = None,
+        pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
+) -> List[StrandPairConstraint]:
+    """
+    Convenience function for creating many constraints as returned by
+    :meth:`nupack_strand_pair_constraint`, one for each threshold specified in parameter `thresholds`,
+    based on number of matching (complementary) domains between pairs of strands.
+
+    Optional parameters `description` and `short_description` are also dicts keyed by the same keys.
+
+    Exactly one of `strands` or `pairs` must be specified. If `strands`, then all pairs of strands
+    (including a strand with itself) will be checked; otherwise only those pairs in `pairs` will be checked.
+
+    It is also common to set different thresholds according to the lengths of the strands.
+    This can be done by calling :meth:`strand_pairs_by_lengths` to separate first by lengths
+    in a dict mapping length pairs to strand pairs,
+    then calling this function once for each (key, value) in that dict, giving the value
+    (which is a list of pairs of strands) as the `pairs` parameter to this function.
+
+    Args:
+        thresholds: Energy thresholds in kcal/mol. If `k` domains are complementary between the strands,
+                    then use threshold `thresholds[k]`.
+        temperature: Temperature in Celsius.
+        sodium: concentration of Na+ in molar
+        magnesium: concentration of Mg++ in molar
+        weight: How much to weigh this :any:`Constraint`.
+        score_transfer_function: See :py:data:`Constraint.score_transfer_function`.
+        descriptions: Long descriptions of constraint suitable for putting into constraint report.
+        short_descriptions: Short descriptions of constraint suitable for logging to stdout.
+        parallel: Whether to test the each pair of :any:`Strand`'s in parallel.
+        strands: Pairs of :any:`Strand`'s to compare; if not specified, checks all pairs in `pairs`.
+                 Mutually exclusive with `pairs`.
+        pairs: Pairs of :any:`Strand`'s to compare; if not specified, checks all pairs in `strands`,
+               including each strand with itself.
+               Mutually exclusive with `strands`.
+
+    Returns:
+        list of constraints, one per threshold in `thresholds`
+    """
+    # ignoring the type error due to this issue: https://github.com/python/mypy/issues/1484
+    # Seems functools.partial with keyword arguments isn't supported well in mypy
+    nupack_strand_pair_constraint_partial: _StrandPairsConstraintCreator = \
+        functools.partial(nupack_strand_pair_constraint, sodium=sodium, magnesium=magnesium)  # type: ignore
+
+    if descriptions is None:
+        descriptions = {
+            num_matching: _pair_default_description('strand', 'NUPACK', threshold, temperature) +
+                          f'\nfor strands with {num_matching} complementary '
+                          f'{"domain" if num_matching==1 else "domains"}'
+            for num_matching, threshold in thresholds.items()
+        }
+
+    if short_descriptions is None:
+        short_descriptions = {
+            num_matching: f'NUPACKpair{num_matching}comp'
+            for num_matching, threshold in thresholds.items()
+        }
+
+    return _strand_pairs_constraints_by_number_matching_domains(
+        constraint_creator=nupack_strand_pair_constraint_partial,
+        thresholds=thresholds,
+        temperature=temperature,
+        weight=weight,
+        score_transfer_function=score_transfer_function,
+        descriptions=descriptions,
+        short_descriptions=short_descriptions,
+        parallel=parallel,
+        strands=strands,
+        pairs=pairs,
+    )
+
+
+def _pair_default_description(part_name: str, func_name: str, threshold: float, temperature: float) -> str:
+    return f'{part_name} pair {func_name} energy >= {threshold} kcal/mol at {temperature}C'
+
+
+def nupack_strand_pair_constraint(
         threshold: float,
         temperature: float = nv.default_temperature,
         sodium: float = nv.default_sodium,
@@ -4786,7 +5056,7 @@ def nupack_strand_pairs_constraint(
 ) -> StrandPairConstraint:
     """
     Returns constraint that checks given pairs of :any:`Strand`'s for excessive interaction using
-    NUPACK's pfunc executable.
+    NUPACK's pfunc function.
 
     NUPACK 4 must be installed. Installation instructions can be found at
     https://piercelab-caltech.github.io/nupack-docs/start/.
@@ -4819,13 +5089,14 @@ def nupack_strand_pairs_constraint(
     _check_nupack_installed()
 
     if description is None:
-        description = f'NUPACK binding energy of strand pair exceeds {threshold} kcal/mol'
+        description = _pair_default_description('strand', 'NUPACK', threshold, temperature)
 
-    def evaluate(seqs: Tuple[str, ...], _: Optional[StrandPair]) -> Tuple[float, str]:
+    def evaluate(seqs: Tuple[str, ...], _: Optional[StrandPair]) -> Result:
         seq1, seq2 = seqs
         energy = nv.binding(seq1, seq2, temperature=temperature, sodium=sodium, magnesium=magnesium)
-        excess = threshold - energy
-        return max(0.0, excess), f'{energy:6.2f} kcal/mol'
+        excess = max(0.0, threshold - energy)
+        value = f'{energy:6.2f} kcal/mol'
+        return Result(excess=excess, value=value)
 
     if pairs is not None:
         pairs = tuple(pairs)
@@ -4837,12 +5108,12 @@ def evaluate(seqs: Tuple[str, ...], _: Optional[StrandPair]) -> Tuple[float, str
                                 parallel=parallel,
                                 pairs=pairs,
                                 evaluate=evaluate,
-                                # summary=summary
                                 )
 
 
-def chunker(sequence: Sequence[T], chunk_length: Optional[int] = None, num_chunks: Optional[int] = None) -> \
-        List[List[T]]:
+def chunker(sequence: Sequence[T],
+            chunk_length: Optional[int] = None,
+            num_chunks: Optional[int] = None) -> List[List[T]]:
     """
     Collect data into fixed-length chunks or blocks, e.g., chunker('ABCDEFG', 3) --> ABC DEF G
 
@@ -4956,19 +5227,18 @@ def rna_duplex_domain_pairs_constraint(
     _check_vienna_rna_installed()
 
     if description is None:
-        description = f'RNAduplex energy for some domain pairs exceeds {threshold} kcal/mol'
+        description = _pair_default_description('domain', 'RNAduplex', threshold, temperature)
 
-    def evaluate_bulk(domain_pairs: Iterable[DomainPair]) -> List[Tuple[DomainPair, float, str]]:
+    def evaluate_bulk(domain_pairs: Iterable[DomainPair]) -> List[Result]:
         sequence_pairs, _, _ = _all_pairs_domain_sequences_complements_names_from_domains(domain_pairs)
-        pairs_scores_summaries: List[Tuple[DomainPair, float, str]] = []
         energies = nv.rna_duplex_multiple(sequence_pairs, logger, temperature, parameters_filename)
 
+        results = []
         for pair, energy in zip(domain_pairs, energies):
             excess = threshold - energy
-            summary = f'{energy:6.2f} kcal/mol'
-            pair_score_summary = (pair, excess, summary)
-            pairs_scores_summaries.append(pair_score_summary)
-        return pairs_scores_summaries
+            value = f'{energy:6.2f} kcal/mol'
+            results.append(Result(excess=excess, value=value))
+        return results
 
     pairs_tuple = None
     if pairs is not None:
@@ -4982,19 +5252,6 @@ def evaluate_bulk(domain_pairs: Iterable[DomainPair]) -> List[Tuple[DomainPair,
                                  pairs=pairs_tuple)
 
 
-# def _normalize_threshold_as_dict(threshold: Union[float, Dict[int, float]],
-#                                  keys: Iterable[int]) -> Dict[int, float]:
-#     # normalize threshold to be a dict with keys `keys` if it is not already, mapping all of the keys
-#     # to its float value if it is a float
-#     if isinstance(threshold, float):
-#         thresholds: Dict[int, float] = {}
-#         for num_comp_domains in keys:
-#             thresholds[num_comp_domains] = threshold
-#         return thresholds
-#     else:
-#         return threshold
-
-
 def _populate_strand_list_and_pairs(strands: Optional[Iterable[Strand]],
                                     pairs: Optional[Iterable[Tuple[Strand, Strand]]]) \
         -> Tuple[List[Strand], List[Tuple[Strand, Strand]]]:
@@ -5085,29 +5342,36 @@ def strand_pairs_by_number_matching_domains(*, strands: Optional[Iterable[Strand
     return strand_pairs
 
 
-class _StrandPairsConstraintCreator(Protocol):
+SPC = TypeVar('SPC',
+              StrandPairConstraint,
+              StrandPairsConstraint)
+
+
+class _StrandPairsConstraintCreator(Protocol[SPC]):
     # Used to specify type of function that
+    #   nupack_strand_pair_constraint
+    # or
     #   rna_duplex_strand_pairs_constraints_by_number_matching_domains
-    # and
+    #   and
     #   rna_cofold_strand_pairs_constraints_by_number_matching_domains
-    # both are. See https://mypy.readthedocs.io/en/stable/protocols.html#callback-protocols
+    # are. See https://mypy.readthedocs.io/en/stable/protocols.html#callback-protocols
     # The Protocol class seems to be available in the typing module, even though the above
     # documentation seems to indicate it is only in typing_extensions?
-    def __call__(self,
+    def __call__(self, *,
                  threshold: float,
                  temperature: float = nv.default_temperature,
                  weight: float = 1.0,
                  score_transfer_function: Callable[[float], float] = default_score_transfer_function,
                  description: Optional[str] = None,
-                 short_description: str = 'rna_dup_strand_pairs',
+                 short_description: str,
                  parallel: bool = False,
                  pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
-                 parameters_filename: str = nv.default_vienna_rna_parameter_filename,
-                 ) -> StrandPairsConstraint: ...
+                 ) -> SPC: ...
 
 
 def _strand_pairs_constraints_by_number_matching_domains(
-        constraint_creator: _StrandPairsConstraintCreator,
+        *,
+        constraint_creator: _StrandPairsConstraintCreator[SPC],
         thresholds: Dict[int, float],
         temperature: float = nv.default_temperature,
         weight: float = 1.0,
@@ -5117,12 +5381,12 @@ def _strand_pairs_constraints_by_number_matching_domains(
         parallel: bool = False,
         strands: Optional[Iterable[Strand]] = None,
         pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
-        parameters_filename: str = nv.default_vienna_rna_parameter_filename,
-):
+) -> List[SPC]:
     # function to share common code between
     #   rna_duplex_strand_pairs_constraints_by_number_matching_domains
     # and
     #   rna_cofold_strand_pairs_constraints_by_number_matching_domains
+
     if strands is None and pairs is None:
         raise ValueError('exactly one of strands or pairs must be specified, but neither is')
     elif strands is not None and pairs is not None:
@@ -5134,29 +5398,22 @@ def _strand_pairs_constraints_by_number_matching_domains(
 
     pairs_by_matching_domains = strand_pairs_by_number_matching_domains(pairs=pairs)
     keys = set(pairs_by_matching_domains.keys())
-    if thresholds is not None:
-        thres_keys = set(thresholds.keys())
-        if keys != thres_keys:
-            raise ValueError(f'The keys of parameter thresholds must be exactly {sorted(list(keys))}, '
-                             'which is the set of integers representing the number of matching domains '
-                             'across all pairs of Strands in the parameter pairs, '
-                             f'but instead the thresholds.keys() is {sorted(list(thres_keys))}')
+    thres_keys = set(thresholds.keys())
+    if keys != thres_keys:
+        raise ValueError(f'The keys of parameter thresholds must be exactly {sorted(list(keys))}, '
+                         'which is the set of integers representing the number of matching domains '
+                         'across all pairs of Strands in the parameter pairs, '
+                         f'but instead the thresholds.keys() is {sorted(list(thres_keys))}')
 
-    constraints: List[StrandPairsConstraint] = []
+    constraints: List[SPC] = []
 
     for num_matching_domains, threshold in thresholds.items():
         pairs_with_matching_domains = pairs_by_matching_domains[num_matching_domains]
 
-        if descriptions is None or num_matching_domains not in descriptions:
-            description = f'RNAduplex energy for pairs strand with {num_matching_domains} ' \
-                          f'matching domains exceeds {threshold} kcal/mol at {temperature} C'
-        else:
-            description = descriptions[num_matching_domains]
-
-        if short_descriptions is None or num_matching_domains not in short_descriptions:
-            short_description = f'RNADup{num_matching_domains}Comp'
-        else:
-            short_description = short_descriptions[num_matching_domains]
+        description = None if descriptions is None \
+            else descriptions.get(num_matching_domains)
+        short_description = None if short_descriptions is None \
+            else short_descriptions.get(num_matching_domains)
 
         constraint = constraint_creator(
             threshold=threshold,
@@ -5167,7 +5424,6 @@ def _strand_pairs_constraints_by_number_matching_domains(
             short_description=short_description,
             parallel=parallel,
             pairs=pairs_with_matching_domains,
-            parameters_filename=parameters_filename,
         )
         constraints.append(constraint)
 
@@ -5175,6 +5431,7 @@ def _strand_pairs_constraints_by_number_matching_domains(
 
 
 def rna_cofold_strand_pairs_constraints_by_number_matching_domains(
+        *,
         thresholds: Dict[int, float],
         temperature: float = nv.default_temperature,
         weight: float = 1.0,
@@ -5184,14 +5441,24 @@ def rna_cofold_strand_pairs_constraints_by_number_matching_domains(
         parallel: bool = False,
         strands: Optional[Iterable[Strand]] = None,
         pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
-        parameters_filename: str = nv.default_vienna_rna_parameter_filename) \
-        -> List[StrandPairsConstraint]:
+        parameters_filename: str = nv.default_vienna_rna_parameter_filename
+) -> List[StrandPairsConstraint]:
     """
     Similar to :meth:`rna_duplex_strand_pairs_constraints_by_number_matching_domains`
     but creates constraints as returned by :meth:`rna_cofold_strand_pairs_constraint`.
     """
+    rna_cofold_with_parameters_filename: _StrandPairsConstraintCreator = \
+        functools.partial(rna_cofold_strand_pairs_constraint,  # type:ignore
+                          parameters_filename=parameters_filename)
+    if descriptions is None:
+        descriptions = {
+            num_matching: _pair_default_description('strand', 'RNAcofold', threshold, temperature) +
+                          f'\nfor strands with {num_matching} complementary '
+                          f'{"domain" if num_matching==1 else "domains"}'
+            for num_matching, threshold in thresholds.items()
+        }
     return _strand_pairs_constraints_by_number_matching_domains(
-        rna_cofold_strand_pairs_constraint,
+        constraint_creator=rna_cofold_with_parameters_filename,
         thresholds=thresholds,
         temperature=temperature,
         weight=weight,
@@ -5201,11 +5468,11 @@ def rna_cofold_strand_pairs_constraints_by_number_matching_domains(
         parallel=parallel,
         strands=strands,
         pairs=pairs,
-        parameters_filename=parameters_filename,
     )
 
 
 def rna_duplex_strand_pairs_constraints_by_number_matching_domains(
+        *,
         thresholds: Dict[int, float],
         temperature: float = nv.default_temperature,
         weight: float = 1.0,
@@ -5215,8 +5482,8 @@ def rna_duplex_strand_pairs_constraints_by_number_matching_domains(
         parallel: bool = False,
         strands: Optional[Iterable[Strand]] = None,
         pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
-        parameters_filename: str = nv.default_vienna_rna_parameter_filename) \
-        -> List[StrandPairsConstraint]:
+        parameters_filename: str = nv.default_vienna_rna_parameter_filename
+) -> List[StrandPairsConstraint]:
     """
     Convenience function for creating many constraints as returned by
     :meth:`rna_duplex_strand_pairs_constraint`, one for each threshold specified in parameter `thresholds`,
@@ -5253,8 +5520,26 @@ def rna_duplex_strand_pairs_constraints_by_number_matching_domains(
     Returns:
         list of constraints, one per threshold in `thresholds`
     """
+    rna_duplex_with_parameters_filename: _StrandPairsConstraintCreator = \
+        functools.partial(rna_duplex_strand_pairs_constraint,  # type:ignore
+                          parameters_filename=parameters_filename)
+
+    if descriptions is None:
+        descriptions = {
+            num_matching: _pair_default_description('strand', 'RNAduplex', threshold, temperature) +
+                          f'\nfor strands with {num_matching} complementary '
+                          f'{"domain" if num_matching==1 else "domains"}'
+            for num_matching, threshold in thresholds.items()
+        }
+
+    if short_descriptions is None:
+        short_descriptions = {
+            num_matching: f'RNAdup{num_matching}comp'
+            for num_matching, threshold in thresholds.items()
+        }
+
     return _strand_pairs_constraints_by_number_matching_domains(
-        rna_duplex_strand_pairs_constraint,
+        constraint_creator=rna_duplex_with_parameters_filename,
         thresholds=thresholds,
         temperature=temperature,
         weight=weight,
@@ -5264,11 +5549,11 @@ def rna_duplex_strand_pairs_constraints_by_number_matching_domains(
         parallel=parallel,
         strands=strands,
         pairs=pairs,
-        parameters_filename=parameters_filename,
     )
 
 
 def rna_duplex_strand_pairs_constraint(
+        *,
         threshold: float,
         temperature: float = nv.default_temperature,
         weight: float = 1.0,
@@ -5277,8 +5562,8 @@ def rna_duplex_strand_pairs_constraint(
         short_description: str = 'rna_dup_strand_pairs',
         parallel: bool = False,
         pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
-        parameters_filename: str = nv.default_vienna_rna_parameter_filename) \
-        -> StrandPairsConstraint:
+        parameters_filename: str = nv.default_vienna_rna_parameter_filename
+) -> StrandPairsConstraint:
     """
     Returns constraint that checks given pairs of :any:`Strand`'s for excessive interaction using
     Vienna RNA's RNAduplex executable.
@@ -5314,38 +5599,37 @@ def rna_duplex_strand_pairs_constraint(
     _check_vienna_rna_installed()
 
     if description is None:
-        description = f'RNAduplex energy for some strand pairs exceeds ' \
-                      f'{threshold} kcal/mol at {temperature} C'
+        description = _pair_default_description('strand', 'RNAduplex', threshold, temperature)
 
-    num_threads = max(cpu_count() - 1, 1)  # this seems to be slightly faster than using all cores
+    num_cores = max(cpu_count(), 1)
 
     # we use ThreadPool instead of pathos because we're farming this out to processes through
     # subprocess module anyway, no need for pathos to boot up separate processes or serialize through dill
-    thread_pool = ThreadPool(processes=num_threads)
-
-    def calculate_energies_unparallel(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]:
-        return nv.rna_duplex_multiple(sequence_pairs, logger, temperature, parameters_filename)
+    thread_pool = ThreadPool(processes=num_cores)
 
-    def calculate_energies(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]:
-        if parallel and len(sequence_pairs) > 1:
-            lists_of_sequence_pairs = chunker(sequence_pairs, num_chunks=num_threads)
-            lists_of_energies = thread_pool.map(calculate_energies_unparallel, lists_of_sequence_pairs)
-            energies = _flatten(lists_of_energies)
+    def calculate_energies(seq_pairs: Sequence[Tuple[str, str]]) -> Tuple[float]:
+        if parallel:
+            energies = nv.rna_duplex_multiple_parallel(thread_pool, seq_pairs, logger, temperature,
+                                                       parameters_filename)
         else:
-            energies = calculate_energies_unparallel(sequence_pairs)
+            energies = nv.rna_duplex_multiple(seq_pairs, logger, temperature, parameters_filename)
         return energies
 
-    def evaluate(strand_pairs: Iterable[StrandPair]) -> List[Tuple[StrandPair, float, str]]:
+    def evaluate_bulk(strand_pairs: Iterable[StrandPair]) -> List[Result]:
         sequence_pairs = [(pair.strand1.sequence(), pair.strand2.sequence()) for pair in strand_pairs]
         energies = calculate_energies(sequence_pairs)
-        pairs_scores_summaries: List[Tuple[StrandPair, float, str]] = []
 
+        for pair, energy in zip(sequence_pairs, energies):
+            if energy > 100:
+                print(f'energy = {energy};  pair = {pair}')
+
+        results = []
         for pair, energy in zip(strand_pairs, energies):
             excess = threshold - energy
-            summary = f'{energy:6.2f} kcal/mol'
-            pair_score_summary = (pair, excess, summary)
-            pairs_scores_summaries.append(pair_score_summary)
-        return pairs_scores_summaries
+            value = f'{energy:6.2f} kcal/mol'
+            result = Result(excess=excess, value=value)
+            results.append(result)
+        return results
 
     pairs_tuple = None
     if pairs is not None:
@@ -5355,7 +5639,7 @@ def evaluate(strand_pairs: Iterable[StrandPair]) -> List[Tuple[StrandPair, float
                                  short_description=short_description,
                                  weight=weight,
                                  score_transfer_function=score_transfer_function,
-                                 evaluate_bulk=evaluate,
+                                 evaluate_bulk=evaluate_bulk,
                                  pairs=pairs_tuple)
 
 
@@ -5377,6 +5661,7 @@ def energy_excess_domains(energy: float,
 
 
 def rna_cofold_strand_pairs_constraint(
+        *,
         threshold: float,
         temperature: float = nv.default_temperature,
         weight: float = 1.0,
@@ -5385,8 +5670,8 @@ def rna_cofold_strand_pairs_constraint(
         short_description: str = 'rna_dup_strand_pairs',
         parallel: bool = False,
         pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
-        parameters_filename: str = nv.default_vienna_rna_parameter_filename) \
-        -> StrandPairsConstraint:
+        parameters_filename: str = nv.default_vienna_rna_parameter_filename
+) -> StrandPairsConstraint:
     """
     Returns constraint that checks given pairs of :any:`Strand`'s for excessive interaction using
     Vienna RNA's RNAduplex executable.
@@ -5431,22 +5716,21 @@ def calculate_energies(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]
         if parallel and len(sequence_pairs) > 1:
             lists_of_sequence_pairs = chunker(sequence_pairs, num_chunks=num_threads)
             lists_of_energies = thread_pool.map(calculate_energies_unparallel, lists_of_sequence_pairs)
-            energies = _flatten(lists_of_energies)
+            energies = flatten(lists_of_energies)
         else:
             energies = calculate_energies_unparallel(sequence_pairs)
         return energies
 
-    def evaluate(strand_pairs: Iterable[StrandPair]) -> List[Tuple[StrandPair, float, str]]:
+    def evaluate_bulk(strand_pairs: Iterable[StrandPair]) -> List[Result]:
         sequence_pairs = [(pair.strand1.sequence(), pair.strand2.sequence()) for pair in strand_pairs]
         energies = calculate_energies(sequence_pairs)
-        pairs_scores_summaries: List[Tuple[StrandPair, float, str]] = []
 
+        results = []
         for pair, energy in zip(strand_pairs, energies):
             excess = threshold - energy
-            summary = f'{energy:6.2f} kcal/mol'
-            pair_score_summary = (pair, excess, summary)
-            pairs_scores_summaries.append(pair_score_summary)
-        return pairs_scores_summaries
+            value = f'{energy:6.2f} kcal/mol'
+            results.append(Result(excess=excess, value=value))
+        return results
 
     pairs_tuple = None
     if pairs is not None:
@@ -5456,7 +5740,7 @@ def evaluate(strand_pairs: Iterable[StrandPair]) -> List[Tuple[StrandPair, float
                                  short_description=short_description,
                                  weight=weight,
                                  score_transfer_function=score_transfer_function,
-                                 evaluate_bulk=evaluate,
+                                 evaluate_bulk=evaluate_bulk,
                                  pairs=pairs_tuple)
 
 
@@ -5494,9 +5778,9 @@ def _all_pairs_domain_sequences_complements_names_from_domains(
     return sequence_pairs, names, domains
 
 
-def _flatten(list_of_lists: Iterable[Iterable[T]]) -> List[T]:
+def flatten(list_of_lists: Iterable[Iterable[T]]) -> Tuple[T]:
     #  Flatten one level of nesting
-    return list(itertools.chain.from_iterable(list_of_lists))
+    return tuple(itertools.chain.from_iterable(list_of_lists))
 
 
 #########################################################################################
@@ -7351,7 +7635,7 @@ def nupack_complex_base_pair_probability_constraint(
     if description is None:
         description = 'Base pair probability of complex'
 
-    def evaluate(seqs: Tuple[str, ...], strand_complex_: Complex) -> Tuple[float, str]:
+    def evaluate(seqs: Tuple[str, ...], strand_complex_: Complex) -> Result:
         assert len(seqs) == len(strand_complex)
         bps = _violation_base_pairs(strand_complex_)
         err_sq = 0.0
@@ -7374,7 +7658,8 @@ def evaluate(seqs: Tuple[str, ...], strand_complex_: Complex) -> Tuple[float, st
                     f'\t{i},{j}: {math.floor(100 * p)}% '
                     f'(<{round(100 * base_type_probability_threshold[t])}% [{t}])')
             summary = '\n'.join(summary_list)
-        return err_sq, summary
+
+        return Result(excess=err_sq, summary=summary, value=err_sq)
 
     def _violation_base_pairs(strand_complex_: Complex) -> List[_BasePair]:
         nupack_complex_result = nv.nupack_complex_base_pair_probabilities(strand_complex_,
diff --git a/nuad/search.py b/nuad/search.py
index e508ec77..153066ff 100644
--- a/nuad/search.py
+++ b/nuad/search.py
@@ -28,6 +28,16 @@
 import datetime
 from functools import lru_cache
 
+try:
+    from typing import Literal
+except ImportError:
+    try:
+        from typing_extensions import Literal
+    except ImportError as err:
+        print('If you are using Python prior to version 3.8, you need to install the typing_extensions\n'
+              'package in order to use nuad. Exiting.')
+        sys.exit(-1)
+
 from tabulate import tabulate
 import numpy.random
 from ordered_set import OrderedSet
@@ -61,11 +71,11 @@
 from nuad.stopwatch import Stopwatch
 
 
-def new_process_pool(cpu_count: int) -> pathos.multiprocessing.Pool:
-    return pathos.multiprocessing.Pool(processes=cpu_count)
+def new_process_pool(cpu_count: int) -> pathos.pools.ProcessPool:
+    return pathos.pools.ProcessPool(processes=cpu_count)
 
 
-_process_pool: pathos.multiprocessing.Pool = new_process_pool(nc.cpu_count())
+_process_pool: Optional[pathos.pools.ProcessPool] = None
 
 log_names_of_domains_and_strands_checked = False
 pprint_indent = 4
@@ -78,8 +88,8 @@ def default_output_directory() -> str:
 # This function takes a lot of time if we don't cache results; but there's not too many different
 # combinations of inputs so it's worth it to maintain a cache.
 @lru_cache()
-def find_parts_to_check(constraint: nc.Constraint, design: nc.Design,
-                        domains_changed: Optional[Tuple[Domain]]) -> Sequence[nc.DesignPart]:
+def find_parts_to_check(constraint: nc.Constraint[DesignPart], design: nc.Design,
+                        domains_changed: Optional[Tuple[Domain]]) -> Tuple[DesignPart, ...]:
     if domains_changed is not None:
         domains_changed_full: OrderedSet[Domain] = OrderedSet(domains_changed)
         for domain in domains_changed:
@@ -92,7 +102,7 @@ def find_parts_to_check(constraint: nc.Constraint, design: nc.Design,
         if len(domains_changed_full) > len(domains_changed):
             domains_changed = tuple(domains_changed_full)
 
-    parts_to_check: Sequence[nc.DesignPart]
+    parts_to_check: Tuple[DesignPart, ...]
     if isinstance(constraint, ConstraintWithDomains):
         parts_to_check = _determine_domains_to_check(design, domains_changed, constraint)
     elif isinstance(constraint, ConstraintWithStrands):
@@ -104,9 +114,11 @@ def find_parts_to_check(constraint: nc.Constraint, design: nc.Design,
     elif isinstance(constraint, ConstraintWithComplexes):
         parts_to_check = _determine_complexes_to_check(domains_changed, constraint)
     elif isinstance(constraint, nc.DesignConstraint):
-        parts_to_check = []  # not used when checking DesignConstraint
+        parts_to_check = tuple()  # not used when checking DesignConstraint
     else:
-        raise NotImplementedError()
+        raise AssertionError('should be unreachable; type of constraint not recognized:\n'
+                             f'type(constraint) = {type(constraint)}\n'
+                             f'constraint = {constraint}')
 
     return parts_to_check
 
@@ -136,7 +148,7 @@ def _at_least_one_domain_unfixed(pair: Tuple[Domain, Domain]) -> bool:
 
 def _determine_domains_to_check(design: Design,
                                 domains_changed: Optional[Tuple[Domain]],
-                                constraint: ConstraintWithDomains) -> Sequence[Domain]:
+                                constraint: ConstraintWithDomains) -> Tuple[Domain, ...]:
     """
     Determines domains to check in `all_domains`.
     If `domains_new` is None, then this is all that are not fixed if constraint.domains
@@ -145,20 +157,20 @@ def _determine_domains_to_check(design: Design,
     it is just those in `domains_new` that appear in `all_domains`.
     """
     # either all pairs, or just constraint.pairs if specified
-    domains_to_check_if_domain_changed_none = design.domains \
+    domains_to_check_if_domain_changed_none = tuple(design.domains) \
         if constraint.domains is None else constraint.domains
 
     # filter out those not containing domain_change if specified
-    domains_to_check = list(domains_to_check_if_domain_changed_none) if domains_changed is None \
-        else [domain for domain in domains_to_check_if_domain_changed_none
-              if domain in domains_changed]
+    domains_to_check = tuple(domains_to_check_if_domain_changed_none) if domains_changed is None \
+        else tuple(domain for domain in domains_to_check_if_domain_changed_none
+                   if domain in domains_changed)
 
     return domains_to_check
 
 
 def _determine_strands_to_check(design: Design,
                                 domains_changed: Optional[Tuple[Domain]],
-                                constraint: ConstraintWithStrands) -> Sequence[Strand]:
+                                constraint: ConstraintWithStrands) -> Tuple[Strand, ...]:
     """
     Similar to _determine_domains_to_check but for strands.
     """
@@ -177,7 +189,7 @@ def _determine_strands_to_check(design: Design,
                     strands_to_check.append(strand)
                     break
 
-    return strands_to_check
+    return tuple(strands_to_check)
 
 
 def _determine_domain_pairs_to_check(design: Design,
@@ -308,7 +320,7 @@ def _strands_containing_domains(domains: Optional[Iterable[Domain]], strands: Li
 _empty_frozen_set: FrozenSet = frozenset()
 
 
-def _independent_domains_in_part(part: nc.DesignPart, exclude_fixed: bool) -> Tuple[Domain]:
+def _independent_domains_in_part(part: DesignPart, exclude_fixed: bool) -> Tuple[Domain]:
     """
     :param part:
         DesignPart (e.g., :any:`Strand`, :any:`Domani`, Tuple[:any:`Strand`, :any:`Strand`])
@@ -414,11 +426,8 @@ def _write_sequences(design: Design, params: SearchParameters, directories: _Dir
 
 def _write_report(params: SearchParameters, directories: _Directories,
                   num_new_optimal_padded: str, eval_set: EvaluationSet) -> None:
-    content = f'''\
-Report on constraints
-=====================
-''' + summary_of_constraints(params.constraints, params.report_only_violations,
-                             eval_set=eval_set)
+    content = summary_of_constraints(params.constraints, params.report_only_violations,
+                                     eval_set=eval_set)
 
     best_filename = directories.best_report_full_filename_noext()
     idx_filename = directories.indexed_report_full_filename_noext(num_new_optimal_padded) \
@@ -669,9 +678,6 @@ class SearchParameters:
 
     report_only_violations: bool = True
     """
-    NOTE: due to a recent change in how the search is conducted, setting this parameter to False is not
-    currently supported. 
-    
     If True, does not give report on each constraint that was satisfied; only reports violations
     and summary of all constraint checks of a certain type (e.g., how many constraint checks there were).
     """
@@ -976,12 +982,15 @@ def search_for_dna_sequences(design: nc.Design, params: SearchParameters) -> Non
 
         _log_constraint_summary(params=params, eval_set=eval_set,
                                 iteration=iteration, num_new_optimal=num_new_optimal)
+        print()
 
     finally:
         # if sys.platform != 'win32':
         #     _pfunc_killall()
-        _process_pool.close()  # noqa
-        _process_pool.terminate()
+        global _process_pool
+        if _process_pool is not None:
+            _process_pool.close()  # noqa
+            _process_pool.terminate()
 
     if directories.debug_file_handler is not None:
         nc.logger.removeHandler(directories.debug_file_handler)  # noqa
@@ -1008,40 +1017,6 @@ def _done(iteration: int, params: SearchParameters, eval_set: EvaluationSet) ->
     return True
 
 
-def create_report(design: nc.Design, constraints: Iterable[Constraint],
-                  report_only_violations: bool = False) -> str:
-    """
-    Returns string containing report of how well `design` does according to `constraints`, assuming
-    `design` has sequences assigned to it, for example, if it was read using
-    :meth:`constraints.Design.from_design_file`
-    from a design.json file writte as part of a call to :meth:`search_for_dna_sequences`.
-
-    The report is the same format as written to the reports generated when calling
-    :meth:`search_for_dna_sequences`.
-    Unfortunately this means it suffers the limitation that it currently only prints a summary
-    for *violations* of constraints;
-    see https://github.com/UC-Davis-molecular-computing/dsd/issues/134
-
-    :param design:
-        the :any:`constraints.Design`, with sequences assigned to all :any:`Domain`'s
-    :param constraints:
-        the list of :any:`constraints.Constraint`'s to evaluate in the report
-    :param report_only_violations:
-        if True, lists only violations of constraints
-    :return:
-        string describing a report of how well `design` does according to `constraints`
-    """
-    evaluation_set = EvaluationSet(constraints, False)
-    evaluation_set.evaluate_all(design)
-
-    content = f'''\
-Report on constraints
-=====================
-''' + summary_of_constraints(constraints, report_only_violations, eval_set=evaluation_set)
-
-    return content
-
-
 def _check_cpu_count(cpu_count: int) -> None:
     # alters number of threads in ThreadPool if cpu count changed. (Lets us hot-swap CPUs, e.g.,
     # in Amazon web services, without stopping the program.)
@@ -1050,8 +1025,9 @@ def _check_cpu_count(cpu_count: int) -> None:
                     f'\nallocating new ThreadPool')
         cpu_count = nc.cpu_count()
         global _process_pool
-        _process_pool.close()
-        _process_pool.terminate()
+        if _process_pool is not None:
+            _process_pool.close()
+            _process_pool.terminate()
         _process_pool = new_process_pool(cpu_count)
 
 
@@ -1278,25 +1254,26 @@ def _find_highest_index_in_directory(directory: str, filename_start: str, ext: s
     return max_index
 
 
-def _pfunc_killall() -> None:
-    import subprocess
-    delim = '#'
-    logger.warning('\n' + delim * 79)
-    logger.warning('# attempting to kill all pfunc processes with `killall pfunc Pfunc`')
-    command = ['killall', 'pfunc', 'Pfunc']
-    process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
-    output, error = process.communicate()
-    if output is None:
-        output = 'None'
-    if error is None:
-        error = 'None'
-    output_str = output.decode() if isinstance(output, bytes) else output
-    error_str = error.decode() if isinstance(error, bytes) else error
-    output_str = textwrap.indent(output_str, delim + '   ')
-    error_str = textwrap.indent(error_str, delim + '   ')
-    logger.warning(f'{delim} output from killall command:\n{output_str}')
-    logger.warning(f'{delim} error from killall command:\n{error_str}')
-    logger.warning('#' * 79)
+# keeping this code around in case we want it for killing parallel processes if necessary
+# def _pfunc_killall() -> None:
+#     import subprocess
+#     delim = '#'
+#     logger.warning('\n' + delim * 79)
+#     logger.warning('# attempting to kill all pfunc processes with `killall pfunc Pfunc`')
+#     command = ['killall', 'pfunc', 'Pfunc']
+#     process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+#     output, error = process.communicate()
+#     if output is None:
+#         output = 'None'
+#     if error is None:
+#         error = 'None'
+#     output_str = output.decode() if isinstance(output, bytes) else output
+#     error_str = error.decode() if isinstance(error, bytes) else error
+#     output_str = textwrap.indent(output_str, delim + '   ')
+#     error_str = textwrap.indent(error_str, delim + '   ')
+#     logger.warning(f'{delim} output from killall command:\n{output_str}')
+#     logger.warning(f'{delim} error from killall command:\n{error_str}')
+#     logger.warning('#' * 79)
 
 
 n_in_last_n_calls = 50
@@ -1315,11 +1292,11 @@ def _log_time(stopwatch: Stopwatch, include_median: bool = False) -> None:
             med_time = statistics.median(time_last_n_calls)
             content += f' median: {med_time:.1f} ms |'
         content_width = len(content)
-        logger.info('-' * content_width + '\n' + content)
+        logger.info('\n' + ('-' * content_width) + '\n' + content)
     else:
         # skip appending first time, since it is much larger and skews the average
         content = f'| time for first call: {stopwatch.milliseconds_str()} ms |'
-        logger.info('-' * len(content) + '\n' + content)
+        logger.info('\n' + ('-' * len(content)) + '\n' + content)
         time_last_n_calls_available = True
 
 
@@ -1341,8 +1318,11 @@ def _log_constraint_summary(*, params: SearchParameters,
         row1 = ['iteration', 'update', 'opt score', 'new score'] + [f'{constraint.short_description}'
                                                                     for constraint in params.constraints]
         header = tabulate([row1], tablefmt='github')
+        if params.scrolling_output and iteration > 0:
+            print()
         print(header)
 
+
     def _dec(score_: float) -> int:
         # how many decimals after decimal point to use given the score
         return max(1, math.ceil(math.log(1 / score_, 10)) + 2) if score_ > 0 else 1
@@ -1383,8 +1363,8 @@ def _dec(score_: float) -> int:
     table_str = tabulate(table, tablefmt='github', numalign='right', stralign='right')
     table_str = _remove_first_lines_from_string(table_str, 2)
     # logger.info(table_str)
-    line_end = '\n' if params.scrolling_output else '\r'
-    print(table_str, end=line_end)
+    first_newline = '' if params.scrolling_output else '\r'
+    print(first_newline + table_str, end='')
 
 
 def assign_sequences_to_domains_randomly_from_pools(design: Design,
@@ -1682,6 +1662,37 @@ def calculate_initial_score_gap(self, design: Design, domains_new: Tuple[Domain]
                 score_gap += ev.score
         return score_gap
 
+    @staticmethod
+    def evaluate_singular_constraint_parallel(constraint: SingularConstraint[DesignPart],
+                                              parts: Tuple[nc.DesignPart, ...],
+                                              score_gap: float) \
+            -> Tuple[List[Tuple[nc.DesignPart, float, str]], float]:
+        if len(parts) == 0:
+            return tuple()
+
+        num_cpus = nc.cpu_count()
+
+        parts_chunks = nc.chunker(parts, num_chunks=num_cpus)
+
+        def call_evaluate_sequential(partz: Tuple[nc.DesignPart]) -> List[Tuple[nc.DesignPart, float, str]]:
+            partz_scores_summaries: List[Tuple[nc.DesignPart, float, str]] = []
+            for part in partz:
+                seqs = tuple(indv_part.sequence() for indv_part in part.individual_parts())
+                score, summary = constraint.call_evaluate(seqs, part)
+                partz_scores_summaries.append((part, score, summary))
+            return partz_scores_summaries
+
+        global _process_pool
+        if _process_pool is None:
+            _process_pool = new_process_pool(num_cpus)
+
+        lists_of_violating_parts_scores_summaries = _process_pool.map(call_evaluate_sequential, parts_chunks)
+        parts_scores_summaries = [elt
+                                  for elts in lists_of_violating_parts_scores_summaries
+                                  for elt in elts]
+
+        return parts_scores_summaries, score_gap
+
     def evaluate_constraint(self,
                             constraint: Constraint[DesignPart],
                             design: Design,  # only used with DesignConstraint
@@ -1690,37 +1701,38 @@ def evaluate_constraint(self,
                             ) -> float:
         # returns score gap = score(old evals) - score(new evals);
         # if gap > 0, then new evals haven't added up to
-        assert ((score_gap is None and domains_new is None) or
-                (score_gap is not None and domains_new is not None))
+        if score_gap is not None:
+            assert domains_new is not None  # ensure we are only using score gap when doing new evaluation
 
         parts = find_parts_to_check(constraint, design, domains_new)
 
         # measure violations of constraints and collect in list of triples (part, score, summary)
-        violating_parts_scores_summaries: List[Tuple[DesignPart, float, str]] = []
+        results: List[nc.Result] = []
         if isinstance(constraint, SingularConstraint):
             if not constraint.parallel or len(parts) == 1 or nc.cpu_count() == 1:
                 for part in parts:
                     seqs = tuple(indv_part.sequence() for indv_part in part.individual_parts())
-                    score, summary = constraint.call_evaluate(seqs, part)
-                    violating_parts_scores_summaries.append((part, score, summary))
-                    if score > 0.0:
+                    result = constraint.call_evaluate(seqs, part)
+                    results.append(result)
+                    if result.score > 0.0:
                         if score_gap is not None:
-                            score_gap -= score
+                            score_gap -= result.score
                             if _is_significantly_greater(0.0, score_gap):
                                 break
             else:
-                raise NotImplementedError('TODO: implement parallelization')
+                violating_parts_scores_summaries, score_gap = \
+                    EvaluationSet.evaluate_singular_constraint_parallel(constraint, parts, score_gap)
 
         elif isinstance(constraint, (BulkConstraint, DesignConstraint)):
             if isinstance(constraint, DesignConstraint):
-                violating_parts_scores_summaries = constraint.call_evaluate_design(design, domains_new)
+                results = constraint.call_evaluate_design(design, domains_new)
             else:
                 # XXX: I don't understand the mypy error on the next line
-                violating_parts_scores_summaries = constraint.call_evaluate_bulk(parts)  # type: ignore
+                results = constraint.call_evaluate_bulk(parts)  # type: ignore
 
             # we can't quit this function early,
             # but we can let the caller know to stop evaluating constraints
-            total_score = sum(score for _, score, _ in violating_parts_scores_summaries)
+            total_score = sum(result.score for result in results)
             if score_gap is not None:
                 score_gap -= total_score
         else:
@@ -1738,17 +1750,18 @@ def evaluate_constraint(self,
             domain_to_evals = self.domain_to_evaluations_new
             domain_to_viols = self.domain_to_violations_new
 
-        for part, score, summary in violating_parts_scores_summaries:
-            domains = _independent_domains_in_part(part, exclude_fixed=False)
-            evaluation = Evaluation(constraint=constraint, part=part, domains=domains,
-                                    score=score, summary=summary, violated=score > 0)
+        for result in results:
+            domains = _independent_domains_in_part(result.part, exclude_fixed=False)
+            evaluation = Evaluation(constraint=constraint, part=result.part, domains=domains,
+                                    score=result.score, summary=result.summary, violated=result.score > 0,
+                                    result=result)
 
-            evals_of_constraint[part] = evaluation
+            evals_of_constraint[result.part] = evaluation
             for domain in domains:
                 domain_to_evals[domain].append(evaluation)
 
             if evaluation.violated:
-                viols_of_constraint[part] = evaluation
+                viols_of_constraint[result.part] = evaluation
                 for domain in domains:
                     domain_to_viols[domain].append(evaluation)
 
@@ -1909,7 +1922,7 @@ def _assert_violations_are_accurate(evaluations: Dict[Constraint, Dict[nc.Part,
             assert not ev.violated
 
 
-@dataclass(frozen=True)
+@dataclass
 class Evaluation(Generic[DesignPart]):
     # Represents a violation of a single :any:`Constraint` in a :any:`Design`.
     # The "part" of the :any:`Design` that was evaluated for the constraint is generic type `DesignPart`
@@ -1931,8 +1944,10 @@ class Evaluation(Generic[DesignPart]):
 
     score: float
 
+    result: nc.Result
+
     def __init__(self, constraint: Constraint, violated: bool, part: DesignPart, domains: Iterable[Domain],
-                 score: float, summary: str) -> None:
+                 score: float, summary: str, result: nc.Result) -> None:
         # :param constraint:
         #     :any:`Constraint` that was violated to result in this
         # :param domains:
@@ -1940,13 +1955,22 @@ def __init__(self, constraint: Constraint, violated: bool, part: DesignPart, dom
         # :param score:
         #     total "score" of this violation, typically something like an excess energy over a
         #     threshold, squared, multiplied by the :data:`Constraint.weight`
-        object.__setattr__(self, 'constraint', constraint)
-        object.__setattr__(self, 'violated', violated)
-        object.__setattr__(self, 'part', part)
-        domains_frozen = frozenset(domains)
-        object.__setattr__(self, 'domains', domains_frozen)
-        object.__setattr__(self, 'score', score)
-        object.__setattr__(self, 'summary', summary)
+        self.constraint = constraint
+        self.violated = violated
+        self.part = part
+        self.domains = frozenset(domains)
+        self.score = score
+        self.summary = summary
+        self.result = result
+
+        # object.__setattr__(self, 'constraint', constraint)
+        # object.__setattr__(self, 'violated', violated)
+        # object.__setattr__(self, 'part', part)
+        # domains_frozen = frozenset(domains)
+        # object.__setattr__(self, 'domains', domains_frozen)
+        # object.__setattr__(self, 'score', score)
+        # object.__setattr__(self, 'summary', summary)
+        # object.__setattr__(self, 'result', result)
 
     def __repr__(self) -> str:
         return f'Evaluation({self.constraint.short_description}, score={self.score:.2f}, ' \
@@ -1967,13 +1991,106 @@ def __eq__(self, other):
 ####################################################################################
 # report generating functions
 
+
+def create_constraints_report(design: nc.Design, constraints: Iterable[Constraint],
+                              report_only_violations: bool = False,
+                              include_only_with_values: bool = False) -> ConstraintsReport:
+    """
+    Returns :any:`ConstraintsReport`, where its :data:`ConstraintsReport.reports` field has one
+    :any:`ConstraintReport` for each :any:`Constraint` in `constraints`,
+    indicating how well `design` does according to `constraints`, assuming
+    `design` has sequences assigned to it, for example, if it was read using
+    :meth:`constraints.Design.from_design_file`
+    from a design.json file written as part of a call to :meth:`search_for_dna_sequences`.
+
+    The report contains the same information as written in the return value of
+    :meth:`summary_of_constraints`, but in a more structred format using :any:`ConstraintReport` objects
+    rather than just text. (:meth:`summary_of_constraints` calls the method :meth:`ConstraintReport.content`
+    on each :any:`ConstraintReport` object).
+
+    :param design:
+        the :any:`constraints.Design`, with sequences assigned to all :any:`Domain`'s
+    :param constraints:
+        the list of :any:`constraints.Constraint`'s to evaluate in the report
+    :param report_only_violations:
+        if True, lists only violations of constraints
+    :param include_only_with_values:
+        if True, lists only violations of constraints
+    :return:
+        :any:`ConstraintsReport` describing a report of how well `design` does
+        according to `constraints`
+    """
+    eval_set = EvaluationSet(constraints, False)
+    eval_set.evaluate_all(design)
+
+    reports = [ConstraintReport(constraint, eval_set, report_only_violations) for constraint in constraints]
+
+    constraints_report = ConstraintsReport(reports=reports,
+                                           total_score=eval_set.total_score,
+                                           total_score_fixed=eval_set.total_score_fixed,
+                                           total_score_nonfixed=eval_set.total_score_nonfixed,
+                                           num_evaluations=eval_set.num_evaluations,
+                                           num_violations=eval_set.num_violations,
+                                           )
+    if include_only_with_values:
+        constraints_report.filter_no_values()
+
+    return constraints_report
+
+
+def create_text_report(design: nc.Design, constraints: Iterable[Constraint],
+                       report_only_violations: bool = False) -> str:
+    """
+    Returns text string containing report of how well `design` does according to `constraints`, assuming
+    `design` has sequences assigned to it, for example, if it was read using
+    :meth:`constraints.Design.from_design_file`
+    from a design.json file written as part of a call to :meth:`search_for_dna_sequences`.
+
+    The report is the same format as written to the reports generated when calling
+    :meth:`search_for_dna_sequences`.
+
+    :param design:
+        the :any:`constraints.Design`, with sequences assigned to all :any:`Domain`'s
+    :param constraints:
+        the list of :any:`constraints.Constraint`'s to evaluate in the report
+    :param report_only_violations:
+        if True, lists only violations of constraints
+    :return:
+        string describing a report of how well `design` does according to `constraints`
+    """
+    constraints_report = create_constraints_report(design, constraints, report_only_violations)
+
+    summaries = [report.content() for report in constraints_report.reports]
+
+    score = constraints_report.total_score
+    score_unfixed = constraints_report.total_score_nonfixed
+    score_total_summary = f'total score of constraint violations: {score:.2f}'
+    score_unfixed_summary = f'total score of unfixed constraint violations: {score_unfixed:.2f}'
+
+    summary = (f'total evaluations: {constraints_report.num_evaluations}\n'
+               f'total violations: {constraints_report.num_violations}\n'
+               + score_total_summary + '\n'
+               + (score_unfixed_summary + '\n\n' if score_unfixed != score else '\n')
+               + '\n\n'.join(summaries))
+
+    return f'''\
+Report on constraints
+=====================
+''' + summary
+
+
+# does essentially the same thing as create_text_report, but is used internally in the search
+# assuming we already have an evaluation set with constraints evaluated; maybe there's an elegant
+# way to share code, but currently ConstraintsReport does not contain the EvaluationSet that is
+# used to populate it.
 def summary_of_constraints(constraints: Iterable[Constraint], report_only_violations: bool,
                            eval_set: EvaluationSet) -> str:
     summaries: List[str] = []
 
     # other constraints
     for constraint in constraints:
-        summary = summary_of_constraint(constraint, report_only_violations, eval_set)
+        report = ConstraintReport(constraint, eval_set, report_only_violations)
+        summary = report.content()
         summaries.append(summary)
 
     score = eval_set.total_score
@@ -1982,117 +2099,225 @@ def summary_of_constraints(constraints: Iterable[Constraint], report_only_violat
     score_unfixed_summary = f'total score of unfixed constraint violations: {score_unfixed:.2f}'
 
     summary = (f'total evaluations: {eval_set.num_evaluations}\n'
-               f'total violations: {eval_set.num_violations}'
+               f'total violations: {eval_set.num_violations}\n'
                + score_total_summary + '\n'
                + (score_unfixed_summary + '\n\n' if score_unfixed != score else '\n')
-               + '\n'.join(summaries))
+               + '\n\n'.join(summaries))
+
+    return f'''\
+Report on constraints
+=====================
+''' + summary
 
-    return summary
 
+@dataclass
+class DisplayOptions:
+    """
+    Options for displaying graphical report in form of histogram of :data:`Report.value` for each
+    constraint evaluation.
+    """
 
-def summary_of_constraint(constraint: Constraint, report_only_violations: bool,
-                          evaluation_set: EvaluationSet) -> str:
-    if isinstance(constraint, (DomainConstraint, StrandConstraint,
-                               DomainPairConstraint, StrandPairConstraint, ComplexConstraint,
-                               DomainsConstraint, StrandsConstraint,
-                               DomainPairsConstraint, StrandPairsConstraint, ComplexesConstraint)):
-        summaries = []
-        num_evals = evaluation_set.num_evaluations_of(constraint)
-        num_viols = evaluation_set.num_violations_of(constraint)
-        num_violations = 0
-        part_type_name = constraint.part_name()
+    bins: int = 10
+    """Number of bins for histogram in report."""
 
-        evals_nonfixed = evaluation_set.evaluations_nonfixed_of_constraint(constraint, report_only_violations)
-        evals_fixed = evaluation_set.evaluations_fixed_of_constraint(constraint, report_only_violations)
 
-        some_fixed_evals = len(evals_fixed) > 0
+def _value_from_constraint_dict(dct: Union[V, Dict[Union[str, Constraint], V]],
+                                constraint: Constraint, default_value: V, klass: type) -> V:
+    # if dct is of type the value we want to return, just return it
+    if isinstance(dct, klass):
+        return dct
 
-        for evals, header_name in [(evals_nonfixed, f"unfixed {part_type_name}s"),
-                                   (evals_fixed, f"fixed {part_type_name}s")]:
-            if len(evals) == 0:
-                continue
+    # otherwise look up value in dictionary
+    assert isinstance(dct, dict)
+    value = default_value
+    if constraint in dct:
+        value = dct[constraint]
+    elif constraint.short_description in dct:
+        value = dct[constraint.short_description]
+    elif constraint.description in dct:
+        value = dct[constraint.description]
+    return value
 
-            max_part_name_length = max(len(violation.part.name) for violation in evals)
-            num_violations += len(evals)
 
-            lines_and_scores: List[Tuple[str, float]] = []
-            for ev in evals:
-                line = f'{part_type_name} {ev.part.name:{max_part_name_length}}: ' \
-                       f'{ev.summary};  score: {ev.score:.2f}'
-                lines_and_scores.append((line, ev.score))
+_default_num_bins = 10
 
-            lines_and_scores.sort(key=lambda line_and_score: line_and_score[1], reverse=True)
 
-            lines = (line for line, _ in lines_and_scores)
-            content = '\n'.join(lines)
+def display_report(design: nc.Design, constraints: Iterable[Constraint],
+                   report_only_violations: bool = False,
+                   layout: Literal['horz', 'vert'] = 'vert',
+                   xlims: Union[Optional[Tuple[float, float]],
+                                Dict[Union[str, Constraint], Optional[Tuple[float, float]]]] = None,
+                   yscale: Literal['log', 'linear', 'symlog'] = 'linear',
+                   bins: Union[int, Dict[Union[str, Constraint], int]] = _default_num_bins) -> None:
+    """
+    When run in a Jupyter notebook cell, creates a :any:`ConstraintsReport` (the one returned from
+    :meth:`create_report`) and displays its data graphically in the notebook using matplotlib.
 
-            # only put header to distinguish fixed from unfixed violations if there are some fixed
-            full_header = _small_header(header_name, "=") if some_fixed_evals else ''
-            summary = full_header + f'\n{content}\n'
-            summaries.append(summary)
+    :param design:
+        the :any:`constraints.Design`, with sequences assigned to all :any:`Domain`'s
+    :param constraints:
+        the list of :any:`constraints.Constraint`'s to evaluate in the report
+    :param report_only_violations:
+        if True, lists only violations of constraints
+    :param layout:
+        layout of plots. If 'horz', they will be laid out horizontally, which is smaller if you have several
+        constraints, but might make it more useful for compare results of different choices of sequences
+        if you call :meth:`display_report` repeatedly for different sequences
+    :param xlims:
+        If specified, is either a single pair of floats to use as the argument to matplotlib.pyplot.xlim:
+        https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.xlim.html
+        or to specify different xlim values for different constraints, can be a
+        dict mapping :any:`Constraint` (or for conveience, string in the fields
+        :data:`Constraint.short_description` or :data:`Constraint.description`) to an xlim pair of floats.
+    :param yscale:
+        same as argument `value` to matplotlib.pyplot.yscale:
+        https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.yscale.html
+    :param bins:
+        same as argument `bins` to matplotlib.pyplot.hist:
+        https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.yscale.html
+        Can either be a single int to apply to all charts, or a dict mapping a constraint to an int,
+        with the same rules as `xlims`
+    """
+    import matplotlib.pyplot as plt
+
+    if xlims is None:
+        xlims = {}
+
+    assert layout in ['horz', 'vert']
+    constraints_report = create_constraints_report(design, constraints, report_only_violations,
+                                                   include_only_with_values=True)
+
+    num_figs = len(constraints_report.reports)
+
+    for i, report in enumerate(constraints_report.reports):
+        quantities = [ev.result.value for ev in report.evaluations if ev.result.value is not None]
+        if len(quantities) == 0:
+            # print(f'skipping constraint "{report.constraint.description}" since its evaluations have no values')
+            continue
+
+        # convert pint.Quantity to unitless magnitude to avoid UnitStrippedWarning when calling py.hist
+        values = [q.magnitude for q in quantities]
+
+        if layout == 'horz':
+            plt.subplot(1, num_figs, i + 1)
+        weights = np.ones(len(values)) / len(values) if yscale == 'linear' else None
+        # the histogram of the data
+        # weights is to display as percentage: https://stackoverflow.com/a/51477080
+
+        num_bins = _value_from_constraint_dict(bins, report.constraint, _default_num_bins, int)
+        _, __, ___ = plt.hist(values, num_bins, density=True, facecolor='g', alpha=0.75, weights=weights, )
+
+        plt.yscale(yscale)
+        if yscale == 'linear':
+            from matplotlib.ticker import PercentFormatter
+            plt.gca().yaxis.set_major_formatter(PercentFormatter(1))
+
+        # see if user set custom x limits for this constraint
+        # not sure why getting mypy error on next line
+        xlim = _value_from_constraint_dict(xlims, report.constraint, None, tuple)  # type:ignore
+        if xlim != (-1, -1):
+            plt.xlim(xlim)
+
+        # label x-axis with units (e.g., kilocalorie / mole)
+        unit = str(quantities[0].units)
+        plt.xlabel(unit)
+
+        plt.title(report.constraint.description)
+
+        if layout == 'vert':
+            plt.show()
+
+    if layout == 'horz':
+        plt.tight_layout(rect=(0, 0, max(1, num_figs), 1))
+        plt.show()
 
-        if report_only_violations and num_violations != num_viols:
-            assert num_violations == num_viols
 
-        content = ''.join(summaries)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_viols, num_evaluations=num_evals)
+@dataclass
+class ConstraintsReport:
+    """Represents a report on how well a design did on all constraints."""
 
-    elif isinstance(constraint, DesignConstraint):
-        raise NotImplementedError()
-    else:
-        content = f'skipping summary of constraint {constraint.description}; ' \
-                  f'unrecognized type {type(constraint)}'
-        report = ConstraintReport(constraint=constraint, content=content, num_violations=0, num_evaluations=0)
+    reports: List[ConstraintReport]
+    """Has one :any:`ConstraintReport` per :any:`Constraint` evaluated."""
 
-    summary = add_header_to_content_of_summary(report, evaluation_set)
-    return summary
+    num_evaluations: int
+    """Total number of evaluations of all :any:`Constraint`'s."""
 
+    num_violations: int
+    """Total number of evaluations of all :any:`Constraint`'s that were violated."""
 
-def add_header_to_content_of_summary(report: ConstraintReport, eval_set: EvaluationSet) -> str:
-    score = eval_set.score_of_constraint(report.constraint, True)
-    score_unfixed = eval_set.score_of_constraint_nonfixed(report.constraint, True)
+    total_score: float
+    """
+    Total "score" of evaluations of all :any:`Constraint`'s (the score is what the search done by 
+    :meth:`search_for_dna_sequences` is trying to minimize).
+    """
 
-    if score != score_unfixed:
-        summary_score_unfixed = f'\n* unfixed score of violations: {score_unfixed:.2f}'
-    else:
-        summary_score_unfixed = None
+    total_score_nonfixed: float
+    """
+    Total "score" of evaluations of all :any:`Constraint`'s that are not "fixed".
+    
+    It should obey
+    :data:`ConstraintsReport.total_score` == :data:`ConstraintsReport.total_score_nonfixed` + 
+    :data:`ConstraintsReport.total_score_fixed`
+    """
 
-    indented_content = textwrap.indent(report.content, '  ')
+    total_score_fixed: float
+    """
+    Total "score" of evaluations of all :any:`Constraint`'s that are "fixed", meaning that 
+    all domains involved in the evaluation of the constraint have :data:`Domain.fixed` == True,
+    so the search will not change the sequences of those domains (so it is not possible to satisfy 
+    the constraint for that evaluation).
+    
+    It should obey
+    :data:`ConstraintsReport.total_score` == :data:`ConstraintsReport.total_score_nonfixed` + 
+    :data:`ConstraintsReport.total_score_fixed`
+    """
 
-    summary = f'''
-**{"*" * len(report.constraint.description)}
-* {report.constraint.description}
-* evaluations: {report.num_evaluations}
-* violations:  {report.num_violations}
-* score of violations: {score:.2f}{"" if summary_score_unfixed is None else summary_score_unfixed}
-{indented_content}'''
+    def filter_no_values(self) -> None:
+        """
+        Removes any reports on constraints that have no values. This is useful for plotting numeric
+        constraints without having non-numeric constraints (such "domains_not_substrings") in the way.
+        """
+        reports_with_values = []
+        score_removed = 0
+        score_nonfixed_removed = 0
+        score_fixed_removed = 0
+        num_evals_removed = 0
+        num_viols_removed = 0
+        for report in self.reports:
+            if has_values(report):
+                reports_with_values.append(report)
+            else:
+                score_removed += report.score
+                score_nonfixed_removed += report.score_nonfixed
+                score_fixed_removed += report.score_fixed
+                num_evals_removed += report.num_evaluations
+                num_viols_removed += report.num_violations
 
-    return summary
+        self.reports = reports_with_values
+        self.total_score -= score_removed
+        self.total_score_nonfixed -= score_nonfixed_removed
+        self.total_score_fixed -= score_fixed_removed
 
 
-def _small_header(header: str, delim: str) -> str:
-    width = len(header)
-    return f'\n{header}\n{delim * width}'
+def has_values(report: ConstraintReport) -> bool:
+    for ev in report.evaluations:
+        if ev.result.value is not None:
+            return True
+    return False
 
 
 @dataclass
-class ConstraintReport:
+class ConstraintReport(Generic[DesignPart]):
     """
     Represents a report on how well a design did on a constraint.
     """
 
-    constraint: Optional['Constraint']
+    constraint: nc.Constraint[DesignPart]
     """
     The :any:`Constraint` to report on. This can be None if the :any:`Constraint` object is not available
     at the time the :meth:`Constraint.generate_summary` function is defined. If so it will be
     automatically inserted by the report generating code."""
 
-    content: str
-    """
-    Summary of constraint information on the :any:`Design`.
-    """
-
     num_violations: int
     """
     Total number of "parts" of the :any:`Design` (e.g., :any:`Strand`'s, pairs of :any:`Domain`'s) that
@@ -2104,3 +2329,102 @@ class ConstraintReport:
     Total number of "parts" of the :any:`Design` (e.g., :any:`Strand`'s, pairs of :any:`Domain`'s) that
     were checked against the constraint.
     """
+
+    score: float
+    score_nonfixed: float
+    score_fixed: float
+
+    evaluations: List[Evaluation]
+    evaluations_nonfixed: List[Evaluation]
+    evaluations_fixed: List[Evaluation]
+
+    report_only_violations: bool
+
+    def __init__(self, constraint: nc.Constraint[DesignPart],
+                 evaluation_set: EvaluationSet, report_only_violations: bool) -> None:
+
+        if not isinstance(constraint, (DomainConstraint, StrandConstraint,
+                                       DomainPairConstraint, StrandPairConstraint, ComplexConstraint,
+                                       DomainsConstraint, StrandsConstraint,
+                                       DomainPairsConstraint, StrandPairsConstraint, ComplexesConstraint,
+                                       DesignConstraint)):
+            raise NotImplementedError(f'unrecognized type {type(constraint)}')
+
+        self.constraint = constraint
+        self.report_only_violations = report_only_violations
+        self.evaluations = evaluation_set.evaluations_of_constraint(constraint,
+                                                                    violations=report_only_violations)
+        self.num_evaluations = evaluation_set.num_evaluations_of(constraint)
+        self.num_violations = evaluation_set.num_violations_of(constraint)
+        self.score = evaluation_set.score_of_constraint(constraint, True)
+        self.score_nonfixed = evaluation_set.score_of_constraint_nonfixed(constraint, True)
+        self.score_fixed = evaluation_set.score_of_constraint_fixed(constraint, True)
+
+        self.evaluations_nonfixed = evaluation_set.evaluations_nonfixed_of_constraint(constraint,
+                                                                                      report_only_violations)
+        self.evaluations_fixed = evaluation_set.evaluations_fixed_of_constraint(constraint,
+                                                                                report_only_violations)
+
+    def header(self) -> str:
+        if self.score != self.score_nonfixed:
+            summary_score_unfixed = f'\n* unfixed score of violations: {self.score_nonfixed:.2f}'
+        else:
+            summary_score_unfixed = None
+
+        summary_unfixed_content = "" if summary_score_unfixed is None else summary_score_unfixed
+        summary = f'''\
+**{"*" * len(self.constraint.description)}
+* {self.constraint.description}
+* evaluations: {self.num_evaluations}
+* violations:  {self.num_violations}
+* score of violations: {self.score:.2f}{summary_unfixed_content}'''
+
+        return summary
+
+    def content_no_header(self) -> str:
+        part_type_name = self.constraint.part_name()
+        some_fixed_evals = len(self.evaluations_fixed) > 0
+
+        summaries = []
+        num_violations_counted = 0
+        for evals, header_name in [(self.evaluations_nonfixed, f"unfixed {part_type_name}s"),
+                                   (self.evaluations_fixed, f"fixed {part_type_name}s")]:
+            if len(evals) == 0:
+                continue
+
+            max_part_name_length = max(len(violation.part.name) for violation in evals)
+            num_violations_counted += len(evals)
+
+            lines_and_scores: List[Tuple[str, float]] = []
+            for ev in evals:
+                line = f'{part_type_name} {ev.part.name:{max_part_name_length}}: ' \
+                       f'{ev.summary};  score: {ev.score:.2f}'
+                lines_and_scores.append((line, ev.score))
+
+            lines_and_scores.sort(key=lambda line_and_score: line_and_score[1], reverse=True)
+
+            lines = (line for line, _ in lines_and_scores)
+            content = '\n'.join(lines)
+
+            # only put header to distinguish fixed from unfixed violations if there are some fixed
+            full_header = _small_header(header_name, "=") if some_fixed_evals else ''
+            summary = full_header + f'\n{content}\n'
+            summaries.append(summary)
+
+        if self.report_only_violations:
+            assert num_violations_counted == self.num_violations
+        else:
+            assert num_violations_counted == self.num_evaluations
+
+        return '\n'.join(summaries)
+
+    def content(self) -> str:
+        header = self.header()
+        content_no_header = self.content_no_header()
+        indented_content = textwrap.indent(content_no_header, '  ')
+        return header + '\n' + indented_content
+
+
+def _small_header(header: str, delim: str) -> str:
+    width = len(header)
+    return f'\n{header}\n{delim * width}'
diff --git a/nuad/vienna_nupack.py b/nuad/vienna_nupack.py
index cf80e448..fce84843 100644
--- a/nuad/vienna_nupack.py
+++ b/nuad/vienna_nupack.py
@@ -9,7 +9,6 @@
 that is much faster than calling :meth:`pfunc` on the same pair of strands).
 """  # noqa
 
-import collections
 import math
 import itertools
 import os
@@ -17,13 +16,13 @@
 import random
 import subprocess as sub
 import sys
-from functools import lru_cache
 from multiprocessing.pool import ThreadPool
-from typing import Sequence, Union, Tuple, List, Dict, Optional, cast, Deque
+from pathos.pools import ProcessPool
+from typing import Sequence, Union, Tuple, List, Optional, Iterable
 
 import numpy as np
 
-import nuad.constraints as dc
+import nuad.constraints as nc
 
 os_is_windows = sys.platform == 'win32'
 
@@ -65,7 +64,6 @@ def calculate_strand_association_penalty(temperature: float, num_seqs: int) -> f
     return adjust * (num_seqs - 1)
 
 
-@lru_cache(maxsize=10_000)
 def pfunc(seqs: Union[str, Tuple[str, ...]],
           temperature: float = default_temperature,
           sodium: float = default_sodium,
@@ -108,8 +106,7 @@ def pfunc(seqs: Union[str, Tuple[str, ...]],
     :return:
         complex free energy ("delta G") of ordered complex with strands in given cyclic permutation
     """
-    if isinstance(seqs, str):
-        seqs = (seqs,)
+    seqs = tupleize(seqs)
 
     try:
         from nupack import pfunc as nupack_pfunc  # type: ignore
@@ -134,7 +131,55 @@ def pfunc(seqs: Union[str, Tuple[str, ...]],
     return dg
 
 
-def nupack_complex_base_pair_probabilities(strand_complex: 'dc.Complex',  # circular import causes problems
+def tupleize(seqs: Union[str, Iterable[str]]) -> Tuple[str, ...]:
+    return (seqs,) if isinstance(seqs, str) else tuple(seqs)
+
+
+def pfunc_parallel(
+        pool: ProcessPool,
+        all_seqs: Sequence[Union[str, Tuple[str, ...]]],
+        temperature: float = default_temperature,
+        sodium: float = default_sodium,
+        magnesium: float = default_magnesium,
+        strand_association_penalty: bool = True,
+) -> Tuple[float]:
+    num_seqs = len(all_seqs)
+    if num_seqs == 0:
+        return tuple()
+
+    all_seqs = tuple(tupleize(seqs) for seqs in all_seqs)
+
+    first_seqs = all_seqs[0]
+
+    bases = sum(len(seq) for seq in first_seqs)
+    num_cores = nc.cpu_count(logical=True)
+
+    # these thresholds were measured empirically; see notebook nuad_parallel_time_trials.ipynb
+    call_sequential = (len(all_seqs) == 1
+                       or (bases <= 30 and num_seqs <= 50)
+                       or (bases <= 40 and num_seqs <= 40)
+                       or (bases <= 50 and num_seqs <= 30)
+                       or (bases <= 75 and num_seqs <= 20)
+                       or (bases <= 100 and num_seqs <= 10)
+                       or (bases <= 125 and num_seqs <= 4)
+                       or (bases <= 150 and num_seqs <= 3)
+                       or (num_seqs <= 1)
+                       )
+
+    def calculate_energies_sequential(all_tuples: Sequence[Tuple[str, ...]]) -> Tuple[float]:
+        return tuple(pfunc(seqs, temperature, sodium, magnesium, strand_association_penalty)
+                     for seqs in all_tuples)
+
+    if call_sequential:
+        return calculate_energies_sequential(all_seqs)
+
+    lists_of_sequence_pairs = nc.chunker(all_seqs, num_chunks=num_cores)
+    lists_of_energies = pool.map(calculate_energies_sequential, lists_of_sequence_pairs)
+    energies = nc.flatten(lists_of_energies)
+    return tuple(energies)
+
+
+def nupack_complex_base_pair_probabilities(strand_complex: 'nc.Complex',  # circular import causes problems
                                            temperature: float = default_temperature,
                                            sodium: float = default_sodium,
                                            magnesium: float = default_magnesium) -> np.ndarray:
@@ -226,33 +271,12 @@ def call_subprocess(command_strs: List[str], user_input: str) -> Tuple[str, str]
     return output_decoded, stderr_decoded
 
 
-# https://github.com/python/cpython/blob/42336def77f53861284336b3335098a1b9b8cab2/Lib/functools.py#L485
-_sentinel = object()
-_max_size = 1_000_000
-_rna_duplex_cache: Dict[Tuple[Tuple[str, str], float, str], float] = {}
-_rna_duplex_queue: Deque[Tuple[Tuple[str, str], float, str]] = collections.deque(maxlen=_max_size)
-
-
-def _rna_duplex_single_cacheable(seq_pair: Tuple[str, str], temperature: float,
-                                 parameters_filename: str) -> Optional[float]:
-    # return None if arguments are not in _rna_duplex_cache,
-    # otherwise return cached energy in _rna_duplex_cache
-    key = (seq_pair, temperature, parameters_filename)
-    result: Union[object, float] = _rna_duplex_cache.get(key, _sentinel)
-    if result is _sentinel:
-        return None
-    else:
-        result_float: float = cast(float, result)
-        return result_float
-
-
-def rna_duplex_multiple(seq_pairs: Sequence[Tuple[str, str]],
+def rna_duplex_multiple(pairs: Sequence[Tuple[str, str]],
                         logger: logging.Logger = logging.root,
                         temperature: float = default_temperature,
                         parameters_filename: str = default_vienna_rna_parameter_filename,
-                        # cache: bool = True,
-                        cache: bool = False,  # off until I implement LRU queue to bound cache size
-                        ) -> List[float]:
+                        max_energy: float = 0.0,
+                        ) -> Tuple[float]:
     """
     Calls RNAduplex (from ViennaRNA package: https://www.tbi.univie.ac.at/RNA/)
     on a list of pairs, specifically:
@@ -260,7 +284,7 @@ def rna_duplex_multiple(seq_pairs: Sequence[Tuple[str, str]],
     where seqi is a string over {A,C,T,G}. Temperature is in Celsius.
     Returns a list (in the same order as seqpairs) of free energies.
 
-    :param seq_pairs:
+    :param pairs:
         sequence (list or tuple) of pairs of DNA sequences
     :param logger:
         logger to use for printing error messages
@@ -268,13 +292,16 @@ def rna_duplex_multiple(seq_pairs: Sequence[Tuple[str, str]],
         temperature in Celsius
     :param parameters_filename:
         name of parameters file for NUPACK
-    :param cache:
-        Whether to cache results to save on number of sequences to give to RNAduplex.
+    :param max_energy:
+        This is the maximum energy possible to assign. If RNAduplex reports any energies larger than this,
+        they will be changed to `max_energy`. This is useful in case two sequences have no possible
+        base pairs between them (e.g., CCCC and TTTT), in which case RNAduplex assigns a free energy
+        of 100000 (perhaps its approximation of infinity). But for meaningful comparison and particularly
+        for graphing energies, it's nice if there's not some value several orders of magnitude larger
+        than all the rest.
     :return:
-        list of free energies, in the same order as `seq_pairs`
+        list of free energies, in the same order as `pairs`
     """
-    # print(f'rna_duplex_multiple.lru_cache = {rna_duplex_multiple.cache_info()}')
-
     # NB: the string NA_parameter_set needs to be exactly the intended filename;
     # e.g. any extra whitespace characters cause RNAduplex to default to RNA parameter set
     # without warning the user!
@@ -285,15 +312,6 @@ def rna_duplex_multiple(seq_pairs: Sequence[Tuple[str, str]],
     # https://stackoverflow.com/questions/10174211/how-to-make-an-always-relative-to-current-module-file-path
 
     # fill in cached energies and determine which indices still need to have their energies calculated
-    if cache:
-        energies = [_rna_duplex_single_cacheable(seq_pair, temperature, parameters_filename)
-                    for seq_pair in seq_pairs]
-    else:
-        energies = [None] * len(seq_pairs)
-    idxs_to_calculate = [i for i, energy in enumerate(energies) if energy is None]
-    idxs_to_calculate_set = set(idxs_to_calculate)
-    seq_pairs_to_calculate = [seq_pair for i, seq_pair in enumerate(seq_pairs) if i in idxs_to_calculate_set]
-
     full_parameters_filename = os.path.join(os.path.dirname(__file__),
                                             parameter_set_directory, parameters_filename)
 
@@ -304,7 +322,7 @@ def rna_duplex_multiple(seq_pairs: Sequence[Tuple[str, str]],
                                '--noGU', '−−noconv']
 
     # DNA sequences to type after RNAduplex starts up
-    user_input = '\n'.join(f'{seq_pair[0]}\n{seq_pair[1]}' for seq_pair in seq_pairs_to_calculate) + '\n@\n'
+    user_input = '\n'.join(f'{seq1}\n{seq2}' for seq1, seq2 in pairs) + '\n@\n'
 
     output, error = call_subprocess(command_strs, user_input)
 
@@ -315,31 +333,61 @@ def rna_duplex_multiple(seq_pairs: Sequence[Tuple[str, str]],
                              'is a different error that I don\'t know how to handle. Exiting...'
                              f'\nerror:\n{error}')
 
-    energies_to_calculate: List[float] = []
     lines = output.split('\n')
-    for line in lines[:-1]:
-        energies_to_calculate.append(float(line.split(':')[1].split('(')[1].split(')')[0]))
-    if len(lines) - 1 != len(seq_pairs_to_calculate):
-        raise ValueError(f'lengths do not match: #lines:{len(lines) - 1} #seqpairs:{len(seq_pairs)}')
+    if len(lines) - 1 != len(pairs):
+        raise ValueError(f'lengths do not match: #lines:{len(lines) - 1} #seqpairs:{len(pairs)}')
 
-    assert len(energies_to_calculate) == len(seq_pairs_to_calculate)
+    energies = []
+    for line in lines[:-1]:
+        energy = float(line.split(':')[1].split('(')[1].split(')')[0])
+        energy = min(energy, max_energy)
+        energies.append(energy)
 
-    # put calculated energies into list to return alongside cached energies
-    assert len(idxs_to_calculate) == len(energies_to_calculate)
-    for i, energy, seq_pair in zip(idxs_to_calculate, energies_to_calculate, seq_pairs_to_calculate):
-        energies[i] = energy
-        if cache:
-            # clear out oldest cache key if _rna_duplex_queue is full
-            if len(_rna_duplex_queue) == _rna_duplex_queue.maxlen:
-                lru_item = _rna_duplex_queue[0]
-                if lru_item in _rna_duplex_cache:
-                    del _rna_duplex_cache[lru_item]
+    return tuple(energies)
 
-            key = (seq_pair, temperature, parameters_filename)
-            _rna_duplex_cache[key] = energy
-            _rna_duplex_queue.append(key)
 
-    return energies
+def rna_duplex_multiple_parallel(
+        thread_pool: ThreadPool,
+        pairs: Sequence[Tuple[str, str]],
+        logger: logging.Logger = logging.root,
+        temperature: float = default_temperature,
+        parameters_filename: str = default_vienna_rna_parameter_filename,
+        max_energy: float = 0.0,
+) -> Tuple[float]:
+    """
+    Parallel version of :meth:`rna_duplex_multiple`. TODO document this
+    """
+    num_pairs = len(pairs)
+    if num_pairs == 0:
+        return tuple()
+
+    bases = len(pairs[0][0] + pairs[0][1])
+    num_cores = nc.cpu_count(logical=True)
+
+    # these thresholds were measured empirically; see notebook nuad_parallel_time_trials.ipynb
+    call_sequential = (len(pairs) == 1
+                       or (bases <= 10 and num_pairs <= 20000)
+                       or (bases <= 15 and num_pairs <= 10000)
+                       or (bases <= 20 and num_pairs <= 5000)
+                       or (bases <= 30 and num_pairs <= 2000)
+                       or (bases <= 40 and num_pairs <= 1000)
+                       or (bases <= 50 and num_pairs <= 800)
+                       or (bases <= 75 and num_pairs <= 200)
+                       or (bases <= 100 and num_pairs <= 150)
+                       or (num_pairs < num_cores)
+                       )
+
+    def calculate_energies_sequential(seq_pairs: Sequence[Tuple[str, str]]) -> Tuple[float]:
+        return rna_duplex_multiple(pairs=seq_pairs, logger=logger, temperature=temperature,
+                                   parameters_filename=parameters_filename, max_energy=max_energy)
+
+    if call_sequential:
+        return calculate_energies_sequential(pairs)
+
+    lists_of_sequence_pairs = nc.chunker(pairs, num_chunks=num_cores)
+    lists_of_energies = thread_pool.map(calculate_energies_sequential, lists_of_sequence_pairs)
+    energies = nc.flatten(lists_of_energies)
+    return tuple(energies)
 
 
 def _fix_filename_windows(parameters_filename: str) -> str:
@@ -423,7 +471,7 @@ def wc(seq: str) -> str:
     return seq.translate(_wctable)[::-1]
 
 
-def complex_free_energy_single_strand(
+def free_energy_single_strand(
         seq: str, temperature: float = default_temperature, sodium: float = default_sodium,
         magnesium: float = default_magnesium) -> float:
     """Computes the complex free energy of a single strand.
@@ -469,7 +517,7 @@ def binding(seq1: str, seq2: str, *, temperature: float = default_temperature,
     if seq1 > seq2:
         seq1, seq2 = seq2, seq1
     return pfunc((seq1, seq2), temperature, sodium, magnesium) - (
-        pfunc(seq1, temperature, sodium, magnesium) + pfunc(seq2, temperature, sodium, magnesium))
+            pfunc(seq1, temperature, sodium, magnesium) + pfunc(seq2, temperature, sodium, magnesium))
 
 
 def random_dna_seq(length: int, bases: Sequence = 'ACTG') -> str:
@@ -572,7 +620,7 @@ def domain_pairwise_concatenated_no_sec_struct(seq: str, seqs: Sequence[str], te
         wc_seq = wc(seq)
         wc_altseq = wc(altseq)
         if parallel:
-            results = [global_thread_pool.apply_async(complex_free_energy_single_strand,
+            results = [global_thread_pool.apply_async(free_energy_single_strand,
                                                       args=(seq1 + seq2, temperature, sodium, magnesium)) for
                        (seq1, seq2) in
                        [(seq, altseq),
@@ -590,30 +638,30 @@ def domain_pairwise_concatenated_no_sec_struct(seq: str, seqs: Sequence[str], te
                 return False
             energy_sum += sum(energies)
         else:
-            seq_alt = complex_free_energy_single_strand(seq + altseq, temperature, sodium, magnesium)
+            seq_alt = free_energy_single_strand(seq + altseq, temperature, sodium, magnesium)
             if seq_alt > concat:
                 return False
-            seq_wcalt = complex_free_energy_single_strand(seq + wc_altseq, temperature, sodium, magnesium)
+            seq_wcalt = free_energy_single_strand(seq + wc_altseq, temperature, sodium, magnesium)
             if seq_wcalt > concat:
                 return False
-            wcseq_alt = complex_free_energy_single_strand(wc_seq + altseq, temperature, sodium, magnesium)
+            wcseq_alt = free_energy_single_strand(wc_seq + altseq, temperature, sodium, magnesium)
             if wcseq_alt > concat:
                 return False
-            wcseq_wcalt = complex_free_energy_single_strand(wc_seq + wc_altseq, temperature, sodium,
-                                                            magnesium)
+            wcseq_wcalt = free_energy_single_strand(wc_seq + wc_altseq, temperature, sodium,
+                                                    magnesium)
             if wcseq_wcalt > concat:
                 return False
-            alt_seq = complex_free_energy_single_strand(altseq + seq, temperature, sodium, magnesium)
+            alt_seq = free_energy_single_strand(altseq + seq, temperature, sodium, magnesium)
             if alt_seq > concat:
                 return False
-            alt_wcseq = complex_free_energy_single_strand(altseq + wc_seq, temperature, sodium, magnesium)
+            alt_wcseq = free_energy_single_strand(altseq + wc_seq, temperature, sodium, magnesium)
             if alt_wcseq > concat:
                 return False
-            wcalt_seq = complex_free_energy_single_strand(wc_altseq + seq, temperature, sodium, magnesium)
+            wcalt_seq = free_energy_single_strand(wc_altseq + seq, temperature, sodium, magnesium)
             if wcalt_seq > concat:
                 return False
-            wcalt_wcseq = complex_free_energy_single_strand(wc_altseq + wc_seq, temperature, sodium,
-                                                            magnesium)
+            wcalt_wcseq = free_energy_single_strand(wc_altseq + wc_seq, temperature, sodium,
+                                                    magnesium)
             if wcalt_wcseq > concat:
                 return False
             energy_sum += (seq_alt + seq_wcalt + wcseq_alt + wcseq_wcalt +
diff --git a/requirements.txt b/requirements.txt
index df8127d3..2b6799f9 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -6,5 +6,7 @@ xlwt
 xlrd
 nupack
 tabulate
+pint
+matplotlib
 scadnano
 typing_extensions; python_version < "3.8"


Test patch:
diff --git a/examples/hamming_dist_test.py b/examples/hamming_dist_test.py
index 043281b2..2435b920 100644
--- a/examples/hamming_dist_test.py
+++ b/examples/hamming_dist_test.py
@@ -19,7 +19,7 @@ def main():
         for domain in strand.domains:
             domain.pool = domain_pool
     params = ns.SearchParameters(
-        constraints=[nc.nupack_strand_complex_free_energy_constraint(
+        constraints=[nc.nupack_strand_free_energy_constraint(
             threshold=-1.5, temperature=52, short_description='StrandSS', parallel=False)],
         out_directory='output/hamming_dist_test',
         # weigh_violations_equally=True