diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index 5e32f45444c981b9bb671caac174815a30eaf2dd..a8d2f90ccffd827ff38366e61337bb91c969540c 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -1044,7 +1044,7 @@ def merge_multiplicities(result_array, max_multiplicity, verbosity=0):
     verbose_print(
         "\tM&S: merge_multiplicities: Merging multiplicities with initial array {} and max multiplicity {}".format(result_array, max_multiplicity),
         verbosity,
-        0,
+        _MS_VERBOSITY_LEVEL,
     )
     for i in range(max_multiplicity, len(result_array))[::-1]:
         result_array[i-1] += result_array[i]
@@ -1052,7 +1052,7 @@ def merge_multiplicities(result_array, max_multiplicity, verbosity=0):
     verbose_print(
         "\tM&S: merge_multiplicities: Merging multiplicities to new array {}".format(result_array),
         verbosity,
-        0,
+        _MS_VERBOSITY_LEVEL,
     )
 
     return result_array
@@ -1085,9 +1085,11 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
     # Use the global Moecache
     global Moecache
 
+    multiplicity_modulator_array = np.array(options["multiplicity_modulator"]) # Modulator array
+
     # Check for length
-    if not len(options["multiplicity_modulator"])==4:
-        msg = "Multiplicity modulator has to have 4 elements. Now it is {}, len: {}".format(options["multiplicity_modulator"], len(options["multiplicity_modulator"]))
+    if not len(multiplicity_modulator_array)==4:
+        msg = "Multiplicity modulator has to have 4 elements. Now it is {}, len: {}".format(multiplicity_modulator_array, len(multiplicity_modulator_array))
         verbose_print(
             msg,
             verbosity,
@@ -1098,7 +1100,9 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
     # Set up some arrays
     full_fractions_array = np.zeros(4) # Meant to contain the real fractions
     weighted_fractions_array = np.zeros(4) # Meant to contain the fractions multiplied by the multiplicity modulator
-    multiplicity_modulator_array = np.array(options["multiplicity_modulator"]) # Modulator array
+
+    # Get max multiplicity
+    max_multiplicity = get_max_multiplicity(multiplicity_modulator_array)
 
     # ... it's better to interpolate the multiplicity and then
     # use a Poisson distribution to calculate the fractions
@@ -1120,7 +1124,6 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
 
         # Fill the multiplicity array
         for n in range(4):
-            new_result[n] = poisson(multiplicity, n, 3, verbosity)
             full_fractions_array[n] = poisson(multiplicity, n, 3, verbosity)
 
         # Normalize it so it fills to one when taking all the multiplicities:
@@ -1161,8 +1164,7 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
             verbosity,
             _MS_VERBOSITY_LEVEL,
         )
-
-        pass
+        result = full_fractions_array
 
     elif options["normalize_multiplicities"] == "norm":
         # Multiply the full_multiplicity_fraction array by the multiplicity_multiplier_array, creating a weighted fractions array
@@ -1178,9 +1180,6 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
         )
 
     elif options["normalize_multiplicities"] == "merge":
-        # Get max multiplicity
-        max_multiplicity = get_max_multiplicity(options["multiplicity_modulator"])
-
         # We first take the full multiplicity array
         # (i.e. not multiplied by multiplier) and do the merging
         result = merge_multiplicities(full_fractions_array, max_multiplicity, verbosity=verbosity)
diff --git a/binarycpython/utils/functions.py b/binarycpython/utils/functions.py
index 0df3abe35e29c35b7ab96d72932c55bf22ffa44d..a8b54fa5865e6d6e40fa5ab8d56d7dcfc1782a01 100644
--- a/binarycpython/utils/functions.py
+++ b/binarycpython/utils/functions.py
@@ -31,7 +31,33 @@ import py_rinterpolate
 ########################################################
 # Unsorted
 ########################################################
-
+def convert_bytes(size):
+    for x in ['bytes', 'KB', 'MB', 'GB', 'TB']:
+        if size < 1024.0:
+            return "%3.1f %s" % (size, x)
+        size /= 1024.0
+
+    return size
+
+def get_size(obj, seen=None):
+    """Recursively finds size of objects"""
+    size = sys.getsizeof(obj)
+    if seen is None:
+        seen = set()
+    obj_id = id(obj)
+    if obj_id in seen:
+        return 0
+    # Important mark as seen *before* entering recursion to gracefully handle
+    # self-referential objects
+    seen.add(obj_id)
+    if isinstance(obj, dict):
+        size += sum([get_size(v, seen) for v in obj.values()])
+        size += sum([get_size(k, seen) for k in obj.keys()])
+    elif hasattr(obj, '__dict__'):
+        size += get_size(obj.__dict__, seen)
+    elif hasattr(obj, '__iter__') and not isinstance(obj, (str, bytes, bytearray)):
+        size += sum([get_size(i, seen) for i in obj])
+    return size
 
 def subtract_dicts(dict_1: dict, dict_2: dict) -> dict:
     """
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 614c652a60fe629c736dfb3cb405611c475cf14c..5c21d811c4d4b8edf72dafddcc31ff2ef4da486c 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -63,6 +63,8 @@ from binarycpython.utils.functions import (
     update_dicts,
     extract_ensemble_json_from_string,
     get_moe_distefano_dataset,
+    get_size,
+    convert_bytes,
 )
 
 import py_rinterpolate
@@ -81,6 +83,7 @@ from binarycpython.utils.distribution_functions import (
     fill_data,
     Moe_de_Stefano_2017_pdf,
     get_max_multiplicity,
+    distribution_constants,
 )
 
 from binarycpython.utils.spacing_functions import (
@@ -890,6 +893,17 @@ class Population:
             # Print current size
             # print("Current size: {}".format(save_que.qsize()))
 
+            # if system_number%10==0:
+            #     print("system_queue_filler: system_number: {}".format(system_number))
+            #     bytes_size_Moecache = get_size(Moecache)
+            #     print("\tsystem_queue_filler: Size moecache: {}".format(convert_bytes(bytes_size_Moecache)))
+
+            #     bytes_size_distribution_constants = get_size(distribution_constants)
+            #     print("\tsystem_queue_filler: Size distribution_constants: {}".format(convert_bytes(bytes_size_distribution_constants)))
+
+            #     bytes_size_self = get_size(dir(self))
+            #     print("\tsystem_queue_filler: Size dir(self): {}".format(convert_bytes(bytes_size_self)))
+
         # Send closing signal to workers. When they receive this they will terminate
         stream_logger.debug(f"Signaling stop to processes")  # DEBUG
         for _ in range(amt_cores):
@@ -1118,6 +1132,18 @@ class Population:
                 ) as f:
                     f.write("RUNNING")
 
+            # if system_number%10==0:
+            #     print("_process_run_population_grid: system_number: {}".format(system_number))
+            #     bytes_size_Moecache = get_size(Moecache)
+            #     print("\t_process_run_population_grid: Size moecache: {}".format(convert_bytes(bytes_size_Moecache)))
+
+            #     bytes_size_distribution_constants = get_size(distribution_constants)
+            #     print("\t_process_run_population_grid: Size distribution_constants: {}".format(convert_bytes(bytes_size_distribution_constants)))
+
+            #     bytes_size_self = get_size(dir(self))
+            #     print("\t_process_run_population_grid: Size dir(self): {}".format(convert_bytes(bytes_size_self)))
+
+
             # Combine that with the other settings
             full_system_dict = self.bse_options.copy()
             full_system_dict.update(system_dict)
@@ -3750,187 +3776,6 @@ eccentricity3=0
             Moecache[key].destroy()
             del Moecache[key]
 
-    def plot_MS_grid_period_distribution(self, options, mass):
-        """
-        Function to plot the period distribution along with its probability for a given mass.
-        """
-
-        import matplotlib.pyplot as plt
-
-        default_options = copy.deepcopy(moe_distefano_default_options)
-
-        # Take the option dictionary that was given and override.
-        options = update_dicts(default_options, options)
-        options['multiplicity'] = 2
-        sampling_range = const(
-            options["ranges"]["logP"][0],
-            options["ranges"]["logP"][1],
-            options["resolutions"]["logP"][0])
-
-        total_prob = 0
-        probs = []
-        step = sampling_range[1]-sampling_range[0]
-        for value in sampling_range:
-            options['M1'] = mass
-            options['P'] = 10**value
-            prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity'])['P']
-            probs.append(prob)
-            total_prob += step * prob
-
-        self._clean_interpolators()
-
-        plt.title("Period distribution for M={}. Total prob: {}".format(mass, total_prob))
-        plt.plot(sampling_range, probs)
-        plt.ylabel("Probability")
-        plt.xlabel("log10(period)")
-        plt.show()
-
-    def plot_MS_grid_q_distribution(self, options, mass, period):
-        """
-        Function to plot the q distribution along with its probability for a given mass.
-        """
-
-        import matplotlib.pyplot as plt
-
-        default_options = copy.deepcopy(moe_distefano_default_options)
-
-        # Take the option dictionary that was given and override.
-        options = update_dicts(default_options, options)
-        options['multiplicity'] = 2
-        sampling_range = const(
-            options["ranges"]["q"][0],
-            options["ranges"]["q"][1],
-            options["resolutions"]["M"][1])
-
-        total_prob = 0
-        probs = []
-        step = sampling_range[1]-sampling_range[0]
-        print(sampling_range)
-        options['M1'] = mass
-        options['P'] = period
-        for value in sampling_range:
-            options['M2'] = value * options['M1']
-            prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity']).get('q', 0)
-            probs.append(prob)
-            total_prob += step * prob
-
-        self._clean_interpolators()
 
-        plt.title("q distribution for M={} P={}. Total prob: {}".format(mass, period, total_prob))
-        plt.plot(sampling_range, probs)
-        plt.ylabel("Probability")
-        plt.xlabel("mass ratio q")
-        plt.show()
-
-    def plot_MS_grid_ecc_distribution(self, options, mass, period):
-        """
-        Function to plot the ecc distribution along with its probability for a given mass.
-        """
-
-        import matplotlib.pyplot as plt
-
-        default_options = copy.deepcopy(moe_distefano_default_options)
-
-        # Take the option dictionary that was given and override.
-        options = update_dicts(default_options, options)
-        options['multiplicity'] = 2
-        sampling_range = const(
-            options["ranges"]["ecc"][0],
-            options["ranges"]["ecc"][1],
-            options["resolutions"]["ecc"][0])
-
-        total_prob = 0
-        probs = []
-        step = sampling_range[1]-sampling_range[0]
-        print(sampling_range)
-        options['M1'] = mass
-        options['P'] = period
-        for value in sampling_range:
-            options['ecc'] = value
-            prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity']).get('ecc', 0)
-            probs.append(prob)
-            total_prob += step * prob
 
-        self._clean_interpolators()
-
-        plt.title("ecc distribution for M={} P={}. Total prob: {}".format(mass, period, total_prob))
-        plt.plot(sampling_range, probs)
-        plt.ylabel("Probability")
-        plt.xlabel("Eccentricity")
-        plt.show()
-
-
-    def plot_MS_grid_ecc_q_distribution(self, options, mass, period):
-        """
-        Function to plot the ecc distribution along with its probability for a given mass.
-        """
-        import numpy as np
-        import matplotlib.pyplot as plt
-
-        default_options = copy.deepcopy(moe_distefano_default_options)
-
-        # Take the option dictionary that was given and override.
-        options = update_dicts(default_options, options)
-        options['multiplicity'] = 2
-        sampling_range_ecc = const(
-            options["ranges"]["ecc"][0],
-            options["ranges"]["ecc"][1],
-            options["resolutions"]["ecc"][0])
-
-        sampling_range_q = const(
-            options["ranges"]["q"][0],
-            options["ranges"]["q"][1],
-            options["resolutions"]["M"][1])
-
-        total_prob = 0
-        probs_all = np.zeros(shape=(options["resolutions"]["ecc"][0], options["resolutions"]["M"][1]))
-        probs_ecc = np.zeros(shape=(options["resolutions"]["ecc"][0], options["resolutions"]["M"][1]))
-        probs_q = np.zeros(shape=(options["resolutions"]["ecc"][0], options["resolutions"]["M"][1]))
-
-        step_ecc = sampling_range_ecc[1]-sampling_range_ecc[0]
-        step_q = sampling_range_q[1]-sampling_range_q[0]
-
-        options['M1'] = mass
-        options['P'] = period
-
-        for i, ecc_value in enumerate(sampling_range_ecc):
-            for j, q_value in enumerate(sampling_range_q):
-                options['ecc'] = ecc_value
-                options['M2'] = q_value * options['M1']
-
-                prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity'])
-
-                probs_all[i][j] = prob.get('total_probdens', 0)
-                probs_ecc[i][j] = prob.get('ecc', 0)
-                probs_q[i][j] = prob.get('q', 0)
-
-        # Calculate total prob ecc
-        total_prob_ecc = 0
-        for ecc_value in sampling_range_ecc:
-            options['ecc'] = ecc_value
-            prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity'])
-            total_prob_ecc += prob.get('ecc', 0) * step_ecc
-
-        # Calculate total prob q:
-        total_prob_q = 0
-        for q_value in sampling_range_q:
-            options['M2'] = q_value * options['M1']
-            prob = Moe_de_Stefano_2017_pdf(options, verbosity=self.grid_options['verbosity'])
-            total_prob_q += prob.get('q', 0) * step_q   
-
-        # Clean the interpolators
-        self._clean_interpolators()
 
-        # Set up plot
-        fig, axes = plt.subplots(nrows=1, ncols=1)
-        extent = [options["ranges"]["ecc"][0] , options["ranges"]["ecc"][1], options["ranges"]["q"][0], options["ranges"]["q"][1]]
-        im = axes.imshow(probs_all, origin='lower', extent=extent)
-        plt.title("ecc distribution for M={} P={}. Total prob ecc: {} Total prob q: {}".format(mass, period, total_prob_ecc, total_prob_q))
-        plt.ylabel("Eccentricity")
-        plt.xlabel("Mass ratio q")
-
-        fig.subplots_adjust(right=0.85)
-        cbar_ax = fig.add_axes([0.88, 0.15, 0.04, 0.7])
-        cb = fig.colorbar(im, cax=cbar_ax)
-        cb.set_label("Probability")
-        plt.show()