Modifying the behaviour of the multiplicty fraction calculating function

binarycpython/utils/distribution_functions.py

@@ -1069,7 +1069,9 @@ def normalize_dict(result_dict, verbosity=0):
 
 def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
     """
-    Function that returns a list of multiplicity fractions for a given list
+    Function that creates a list of probability fractions and
+    normalizes and merges them according to the users choice.
+
 
     TODO: make an extrapolation functionality in this. log10(1.6e1)
     is low, we can probably go a bit further
@@ -1083,8 +1085,20 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
     # Use the global Moecache
     global Moecache
 
-    result = {}
-    new_result = {}
+    # 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"]))
+        verbose_print(
+            msg,
+            verbosity,
+            0,
+        )
+        raise ValueError(msg)
+
+    # 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
 
     # ... it's better to interpolate the multiplicity and then
     # use a Poisson distribution to calculate the fractions
@@ -1104,25 +1118,16 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
             [np.log10(options["M1"])]
         )[0]
 
+        # Fill the multiplicity array
         for n in range(4):
             new_result[n] = poisson(multiplicity, n, 3, verbosity)
-            result[n] = (
-                options["multiplicity_modulator"][n]
-                * poisson(multiplicity, n, 3, verbosity)
-                if options["multiplicity_modulator"][n] > 0
-                else 0
-            )
-        print("For mass M1={}".format(options['M1']))
-        print("Multiplicty modulator: {}".format(options["multiplicity_modulator"]))
-        print("Old results: {}".format(result))
+            full_fractions_array[n] = poisson(multiplicity, n, 3, verbosity)
 
-        # Normalize dict
-        new_result = normalize_dict(new_result, verbosity=verbosity)
-        sum_result = sum([new_result[key] for key in new_result])
-        print("New result: {} (sum = {})".format(new_result, sum_result))
+        # Normalize it so it fills to one when taking all the multiplicities:
+        full_fractions_array = full_fractions_array/np.sum(full_fractions_array)
 
         verbose_print(
-            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: using model {}".format("Poisson"),
+            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: using model {}: full_fractions_array: {}".format("Poisson", full_fractions_array),
             verbosity,
             _MS_VERBOSITY_LEVEL,
         )
@@ -1133,27 +1138,26 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
         # note that in this case, there are no quadruples: these
         # are combined with triples
 
+        # Fill with the raw values
         for n in range(3):
-            result[n] = (
+            full_fractions_array[n] = (
                 Moecache["rinterpolator_multiplicity"].interpolate(
                     [np.log10(options["M1"])]
-                )[n + 1]
-                * options["multiplicity_modulator"][n]
-            )
-        result[3] = 0.0  # no quadruples
+                )[n + 1])
+
+        # Set last value
+        full_fractions_array[3] = 0.0  # no quadruples
         verbose_print(
-            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: using model {}".format("data"),
+            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: using model {}: full_fractions_array: {}".format("data", full_fractions_array),
             verbosity,
             _MS_VERBOSITY_LEVEL,
         )
 
     # Normalisation:
     if options["normalize_multiplicities"] == "raw":
-        # TODO: Not sure what actually happens here. It eems the same as above
-        # do nothing : use raw Poisson predictions
-
+        # Don't multiply by the multiplicity_array, but do give a fractions array
         verbose_print(
-            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Not normalizing (using raw results)",
+            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Not normalizing (using raw results): results: {}".format(full_fractions_array),
             verbosity,
             _MS_VERBOSITY_LEVEL,
         )
@@ -1161,34 +1165,37 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
         pass
 
     elif options["normalize_multiplicities"] == "norm":
-        # simply normalize so all multiplicities add to 1
-        result = normalize_dict(result, verbosity=verbosity)
+        # Multiply the full_multiplicity_fraction array by the multiplicity_multiplier_array, creating a weighted fractions array
+        weighted_fractions_array = full_fractions_array * multiplicity_modulator_array
+
+        # Normalize this so it is intotal 1:
+        result = weighted_fractions_array/np.sum(weighted_fractions_array)
 
         verbose_print(
-            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Normalizing with {}".format("norm"),
+            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Normalizing with {}. result: {}".format("norm", result),
             verbosity,
             _MS_VERBOSITY_LEVEL,
         )
 
     elif options["normalize_multiplicities"] == "merge":
+        # Get max multiplicity
         max_multiplicity = get_max_multiplicity(options["multiplicity_modulator"])
 
-        # Merge the multiplicities
-        result = merge_multiplicities(result, max_multiplicity, verbosity=verbosity)
-
-        new_result = merge_multiplicities(new_result, max_multiplicity, verbosity=verbosity)
+        # 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)
 
-        # Normalising results again
-        result = normalize_dict(result, verbosity=verbosity)
+        # Then normalize to be sure
+        result = result/np.sum(result)
 
         verbose_print(
-            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Normalizing with {}, max_multiplicity={}".format("merge", max_multiplicity),
+            "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: Normalizing with {}, max_multiplicity={} result={}".format("merge", max_multiplicity, result),
             verbosity,
             _MS_VERBOSITY_LEVEL,
         )
 
     verbose_print(
-        "\tM&S: Multiplicity array: {}".format(str(result)),
+        "\tM&S: Moe_de_Stefano_2017_multiplicity_fractions: {}".format(str(result)),
         verbosity,
         _MS_VERBOSITY_LEVEL,
     )
@@ -1196,7 +1203,6 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
     # return array reference
     return result
 
-
 def build_q_table(options, m, p, verbosity=0):
     ############################################################
     #

binarycpython/utils/grid_options_defaults.py

@@ -603,7 +603,7 @@ combined (and there are NO quadruples).
 e.g. [1,0,0,0] for single stars only
      [0,1,0,0] for binary stars only
 
-defaults to [1,1,1,1] i.e. all types
+defaults to [1,1,0,0] i.e. singles and binaries
 """,
     "normalize_multiplicities": """
 given a mix of multiplities, you can either (noting that
@@ -628,13 +628,16 @@ here (S,B,T,Q) = appropriate modulator * model(S,B,T,Q) )
           Note: if multiplicity_modulator == [1,1,1,1] this
                 option does nothing (equivalent to 'raw').
 
-note: if you only set one multiplicity_modulator
-to 1, and all the others to 0, then normalizing
+Note: If you only set one multiplicity modulator to 1, and all others to zero,
+we will not multiply the fractions by the modulator. The grid will run fine (i.e. it will only run e.g. singles), but it will still take into account the real fraction
+This is unlike how this setting was previously implemented: normalizing
 will mean that you effectively have the same number
 of stars as single, binary, triple or quad (whichever
 is non-zero) i.e. the multiplicity fraction is ignored.
-This is probably not useful except for
-testing purposes or comparing to old grids.
+
+
+
+
 """,
     "q_low_extrapolation_method": """
 q extrapolation (below 0.15) method