new interpolator

bin/posydon-setup-grid

@@ -83,6 +83,8 @@ def find_inlist_from_scenario(source, gitcommit, system_type):
         else:
             Pwarn("git repository is already there, using that",
                   "OverwriteWarning")
+            Pwarn("git repository is already there, using that",
+                  "OverwriteWarning")
 
         inlists_dir = '{0}/.posydon_mesa_inlists'.format(os.environ['HOME'])
         branch = gitcommit.split('-')[0]
@@ -112,6 +114,8 @@ def find_inlist_from_scenario(source, gitcommit, system_type):
         else:
             Pwarn("git repository is already there, using that",
                   "OverwriteWarning")
+            Pwarn("git repository is already there, using that",
+                  "OverwriteWarning")
 
         inlists_dir = '{0}/.user_mesa_inlists'.format(os.environ['HOME'])
         branch = gitcommit.split('-')[0]
@@ -146,6 +150,7 @@ def find_inlist_from_scenario(source, gitcommit, system_type):
                                  stderr = subprocess.PIPE
                              )
     proc.wait()
+    proc.wait()
 
     # if this is looking at posydon defaults, all posydon defaults build from default common inlists
     if source == 'posydon':

posydon/active_learning/psy_cris/run_params/run_psycris_sequence.py

@@ -144,3 +144,6 @@ def main():
 
 if __name__ == "__main__":
     main()
+
+if __name__ == "__main__":
+    main()

posydon/binary_evol/MESA/step_mesa.py

@@ -280,10 +280,10 @@ def load_Interp(self, filename):
 
         # Check if interpolation files exist
         if not os.path.exists(filename):
-            data_download()
+            data_download() #TODO: specify dataset
 
         # Load interpolator
-        self._Interp = IFInterpolator()
+        self._Interp = IFInterpolator(load = True)
         self._Interp.load(filename=filename)
 
     def close(self):

posydon/binary_evol/SN/profile_collapse.py

@@ -393,7 +393,6 @@ def do_core_collapse_BH(star,
             The mass of the disk radiated away in M_sun.
         'BZ_jet_power_total' : float
             The total Blandford-Znajek jet power in erg/s.
-
         # Additional keys that are not used in the current implementation:
         # 'BZ_jet_power_array' : np.array(BZ_jet_power_array),
         #       Blandford-Znajek jet power at each shell collapse in erg/s
@@ -645,6 +644,26 @@ def do_core_collapse_BH(star,
         BZ_jet_power_array.append(BZ_power)
         BZ_jet_power_total += BZ_power
 
+        # calculate the potential BZ jet power at this moment of the collapse
+        # We assume full efficiency for the magnetic flux and a BH spin
+        # dependence of a^2 for the BH spin efficiency.
+        # just an energy total per collapse step.
+        BZ_power = BZ_jet_power(M_dot=dm_disk,
+                     eta_phi=1,
+                     eta_a=a_BH**2)
+        BZ_jet_power_array.append(BZ_power)
+        BZ_jet_power_total += BZ_power
+
+        # calculate the potential BZ jet power at this moment of the collapse
+        # We assume full efficiency for the magnetic flux and a BH spin
+        # dependence of a^2 for the BH spin efficiency.
+        # just an energy total per collapse step.
+        BZ_power = BZ_jet_power(M_dot=dm_disk,
+                     eta_phi=1,
+                     eta_a=a_BH**2)
+        BZ_jet_power_array.append(BZ_power)
+        BZ_jet_power_total += BZ_power
+
         # Append all quantities to the arrays
         J_accreted_array.append(J_BH)
         M_BH_array.append(M_BH)

posydon/interpolation/IF_interpolation.py

@@ -196,9 +196,6 @@ class relies on the BaseIFInterpolator class to perform the interpolation
     find_constraints_to_apply,
     sanitize_interpolated_quantities,
 )
-from posydon.interpolation.data_scaling import DataScaler
-from posydon.utils.posydonwarning import Pwarn
-from posydon.visualization.plot_defaults import DEFAULT_LABELS
 
 
 # INITIAL-FINAL INTERPOLATOR
@@ -273,13 +270,16 @@ def evaluate(self, binary, sanitization_verbose=False):
         """
         ynums = {}
         ycats = {}
-
+        # s = time.time()
         for interpolator in self.interpolators:
             ynum, ycat = interpolator.evaluate(binary, sanitization_verbose)
 
             ynums = {**ynums, **ynum}
             ycats = {**ycats, **ycat}
 
+        # e = time.time()
+        # print(f"Iterated over {len(self.interpolators)} interpolators in {e - s}")
+
         return ynums, ycats
 
 

posydon/interpolation/constraints.py

@@ -49,10 +49,14 @@
 
 import numpy as np
 
-from posydon.utils.common_functions import (
-    orbital_separation_from_period,
-    stefan_boltzmann_law,
-)
+CLASSIFICATION_KEYS = [
+    "S<*>_state",
+    "mt_hist",
+    "S<*>_MOD<n>_SN_type",
+    "S<*>_MOD<n>_CO_type"
+]
+
+N_MODELS = 11 # how many super nova models are there?
 
 # toggle this flag to enable/disable constraints (used for debugging)
 INTERPOLATION_CONSTRAINTS_ON = True
@@ -523,3 +527,63 @@ def sanitize_interpolated_quantities(fvalues, constraints, verbose=False):
                                         constraint["constraint"])
 
     return sanitized
+
+
+def mt_constraint(classes):
+
+    interpolation_class = classes["interpolation_class"]
+
+    if interpolation_class == "initial_MT":
+        classes["mt_hist"] == "ini_RLO"
+    elif interpolation_class == "no_MT":
+        classes["mt_hist"] = "no_RLO"
+    elif interpolation_class == "stable_MT":
+        pass
+    elif interpolation_class == "unstable_MT":
+        pass
+    elif interpolation_class == "stable_reverse_MT":
+        pass
+
+
+CLASS_CONSTRAINTS = {
+    "S<*>_state": None,
+    "mt_hist": mt_constraint,
+    "S<*>_MOD<n>_SN_type": None,
+    "S<*>_MOD<n>_CO_type": None
+}
+
+def apply_class_constraint(key_name, classes):
+
+    if key_name not in classes.keys():
+        return
+    else:
+        CLASS_CONSTRAINTS[key_name](classes)
+
+def sanitize_classes(classes, ):
+
+    assert(type(classes) == dict)
+
+    if "interpolation_class" not in classes.keys():
+        raise ValueError(
+            "Interpolation class must be present as a classified quantity to enforce classification constraints!"
+        )
+
+    for key in CLASSIFICATION_KEYS:
+        if "<*>" in key:
+
+            for star in range(2):
+                key_name = key.replace("<*>", f"{star}")
+
+                if "MOD<n>" in key_name:
+
+                    for model in range(N_MODELS):
+                        key_name = key_name.replace("<n>", f"{model}")
+
+                        apply_class_constraint(key_name, classes)
+
+                else:
+                    apply_class_constraint(key_name, classes)
+        else:
+
+            apply_class_constraint(key, classes)
+

posydon/interpolation/data_scaling.py

@@ -6,8 +6,27 @@
 ]
 
 
+import sys
+import warnings
+
 import numpy as np
 
+# Convert UserWarning to an error
+warnings.simplefilter("error", RuntimeWarning)
+
+eps = 1.0e-16
+
+SCALING_OPTIONS = [
+    "none",
+    "min_max",
+    "max_abs",
+    # "standardize",
+    "log_min_max", # has
+    # "neg_log_min_max", # has
+    "log_max_abs", # has
+    # "log_standardize", # has
+    # "neg_log_standardize" # has
+]
 
 class DataScaler:
     """Data Normalization class.
@@ -68,27 +87,28 @@ def fit(self, x, method='none', lower=-1.0, upper=1.0):
         if method == 'min_max':
             assert upper > lower, "upper must be greater than lower"
             self.lower, self.upper = lower, upper
-            self.params = [x.min(axis=0), x.max(axis=0)]
+            self.params = [np.nanmin(x, axis=0), np.nanmax(x, axis=0)]
         elif method == 'log_min_max':
             assert upper > lower, "upper must be greater than lower"
             self.lower, self.upper = lower, upper
-            self.params = [np.log10(x.min(axis=0)), np.log10(x.max(axis=0))]
+            self.params = [self.log(np.nanmin(x, axis=0)), self.log(np.nanmax(x, axis=0))]
+
         elif method == 'neg_log_min_max':
             assert upper > lower, "upper must be greater than lower"
             self.lower, self.upper = lower, upper
-            self.params = [np.log10((-x).min(axis=0)),
-                           np.log10((-x).max(axis=0))]
+            self.params = [self.log(np.nanmin(-x, axis=0)),
+                           self.log(np.nanmax(-x, axis=0))]
         elif method == 'max_abs':
-            self.params = [np.abs(x).max(axis=0)]
+            self.params = [np.nanmax(np.abs(x), axis=0)]
         elif method == 'log_max_abs':
-            self.params = [np.abs(np.log10(x)).max(axis=0)]
-        elif method == 'standarize':
-            self.params = [x.mean(axis=0), x.std(axis=0)]
-        elif method == 'log_standarize':
+            self.params = [np.nanmax(np.abs(self.log(x)), axis=0)]
+        elif method == 'standardize':
+            self.params = [np.nanmean(x, axis=0), np.nanstd(x, axis=0)]
+        elif method == 'log_standardize':
             # log will be computed in transform again
-            self.params = [np.log10(x).mean(axis=0), np.log10(x).std(axis=0)]
-        elif method == 'neg_log_standarize':    # log(-x)
-            self.params = [np.log10(-x).mean(axis=0), np.log10(-x).std(axis=0)]
+            self.params = [np.nanmean(self.log(x), axis=0), np.nanstd(self.log(x), axis=0)]
+        elif method == 'neg_log_standardize':    # log(-x)
+            self.params = [np.nanmean(self.log(-x), axis=0), np.nanstd(self.log(-x), axis=0)]
         elif method == 'log':
             self.params = []
         elif method == 'none':  # no transformation
@@ -124,26 +144,26 @@ def transform(self, x):
             x_t = ((x - self.params[0]) / (self.params[1] - self.params[0])
                    * (self.upper - self.lower) + self.lower)
         elif self.method == 'log_min_max':
-            x_t = ((np.log10(x) - self.params[0])
+            x_t = ((self.log(x) - self.params[0])
                    / (self.params[1] - self.params[0])
                    * (self.upper - self.lower) + self.lower)
         elif self.method == 'neg_log_min_max':
-            x_t = ((np.log10(-x) - self.params[0])
+            x_t = ((self.log(-x) - self.params[0])
                    / (self.params[1] - self.params[0])
                    * (self.upper - self.lower) + self.lower)
         elif self.method == 'max_abs':
             x_t = x / self.params[0]
         elif self.method == 'log_max_abs':
-            x_t = np.log10(x) / self.params[0]
-        elif self.method == 'standarize':
+            x_t = self.log(x) / self.params[0]
+        elif self.method == 'standardize':
             x_t = (x - self.params[0]) / self.params[1]
-        elif self.method == 'log_standarize':
+        elif self.method == 'log_standardize':
             # log will be computed in transform again
-            x_t = (np.log10(x) - self.params[0]) / self.params[1]
-        elif self.method == 'neg_log_standarize':
-            x_t = (np.log10(-x) - self.params[0]) / self.params[1]
+            x_t = (self.log(x) - self.params[0]) / self.params[1]
+        elif self.method == 'neg_log_standardize':
+            x_t = (self.log(-x) - self.params[0]) / self.params[1]
         elif self.method == 'log':
-            x_t = np.log10(x)
+            x_t = self.log(x)
         else:  # no transformation
             x_t = x
 
@@ -201,24 +221,38 @@ def inv_transform(self, x_t):
                  / (self.upper - self.lower)
                  * (self.params[1] - self.params[0]) + self.params[0])
         elif self.method == 'log_min_max':
-            x = 10 ** ((x_t - self.lower) / (self.upper - self.lower)
+            x = self.unlog((x_t - self.lower) / (self.upper - self.lower)
                        * (self.params[1] - self.params[0]) + self.params[0])
         elif self.method == 'neg_log_min_max':
-            x = -10 ** ((x_t - self.lower) / (self.upper - self.lower)
+            x = -self.unlog((x_t - self.lower) / (self.upper - self.lower)
                         * (self.params[1] - self.params[0]) + self.params[0])
         elif self.method == 'max_abs':
             x = x_t * self.params[0]
         elif self.method == 'log_max_abs':
-            x = 10 ** (x_t * self.params[0])
+            x = self.unlog(x_t * self.params[0])
         elif self.method == 'standarize':
             x = x_t * self.params[1] + self.params[0]
         elif self.method == 'log_standarize':
-            x = 10 ** (x_t * self.params[1] + self.params[0])
+            x = self.unlog(x_t * self.params[1] + self.params[0])
         elif self.method == 'neg_log_standarize':
-            x = -10 ** (x_t * self.params[1] + self.params[0])
+            x = -self.unlog(x_t * self.params[1] + self.params[0])
         elif self.method == 'log':
-            x = 10 ** x_t
+            x = self.unlog(x_t)
         else:  # no transformation
             x = x_t
 
         return x
+
+    def log(self, x):
+        logged = None
+        try:
+            logged = np.log10(x + eps)
+        except RuntimeWarning:
+            print(self.method)
+            print(x, np.isinf(x).any(), np.isnan(x).any(), (x < 0).any(), np.nanmin(x))
+            # sys.exit()
+
+        return logged
+
+    def unlog(self, x):
+        return (10 ** x) - eps