diff --git a/README.md b/README.md index 1520bd9..78fed46 100644 --- a/README.md +++ b/README.md @@ -26,8 +26,8 @@ Several options can be passed to mkipp.Kipp_Args(): class Kipp_Args: def __init__(self, logs_dirs = ['LOGS'], - profile_names = [], - history_names = [], + profile_paths = [], + history_paths = [], clean_data = True, extra_history_cols = [], identifier = "eps_nuc", @@ -57,10 +57,10 @@ class Kipp_Args: All arguments are optional, if not provided defaults are assigned Args: - logs_dir (List[str]): List of paths to MESA LOGS directories. If profile_names and - history_names are not provided, they are automatically generated from logs_dir. - profile_names (List[str]): List of paths to MESA profile files. - history_names (List[str]): List of paths to MESA history files. + logs_dir (List[str]): List of paths to MESA LOGS directories. If profile_paths and + history_paths are not provided, they are automatically generated from logs_dir. + profile_paths (List[str]): List of paths to MESA profile files. + history_paths (List[str]): List of paths to MESA history files. clean_data (bool): Clean history files in case of redos. If data has no redos, a small performance gain can be had by setting this to False. extra_history_cols (List[str]): Additional column names to be extracted from history files. @@ -134,6 +134,7 @@ The data from mixing regions and profiles can be extracted using lower level fun ```python import mkipp import kipp_data +import mesa_data import matplotlib.pyplot as plt from matplotlib.patches import PathPatch import numpy as np @@ -141,10 +142,9 @@ import numpy as np kipp_args = mkipp.Kipp_Args() fig = plt.figure() axis = plt.gca() -profile_names = ["LOGS/profile"+str(int(i))+".data" \ - for i in np.loadtxt("LOGS/profiles.index", skiprows = 1, usecols = (2,))] +profile_paths = mesa_data.get_profile_paths(["LOGS"]) #if data is distributed among several history.data files, you can provide them -history_names = ["LOGS/history.data"] +history_paths = ["LOGS/history.data"] #read profile data #kipp_data.get_xyz_data returns an object containing # xyz_data.xlims : limits of data in x coordinate @@ -152,7 +152,7 @@ history_names = ["LOGS/history.data"] # xyz_data.Y : 2D array of xaxis values of profile data # xyz_data.Z : 2D array of xaxis values of profile data # the last three can be used as inputs for matplotlib contour or contourf -xyz_data = kipp_data.get_xyz_data(profile_names, kipp_args.xaxis_divide, kipp_args) +xyz_data = kipp_data.get_xyz_data(profile_paths, kipp_args.xaxis_divide, kipp_args) #read mixing regions #kipp_data.get_mixing_zones returns an object containing # mixing_zones.zones : matplotlib Path objects for each mixing zone. @@ -162,7 +162,7 @@ xyz_data = kipp_data.get_xyz_data(profile_names, kipp_args.xaxis_divide, kipp_ar # mixing_zones.y_coords : y coordinates for points at the surface # mixing_zones.histories : mesa_data history files to access additional data # the last three can be used as inputs for matplotlib contour or contourf -mixing_zones = kipp_data.get_mixing_zones(history_names, kipp_args.xaxis_divide, xyz_data.xlims, kipp_args) +mixing_zones = kipp_data.get_mixing_zones(history_paths, kipp_args.xaxis_divide, xyz_data.xlims, kipp_args) # just plot convection, overshooting and semiconvection for i,zone in enumerate(mixing_zones.zones): color = "" diff --git a/example.py b/example.py index 0dd34d8..6c84bcb 100755 --- a/example.py +++ b/example.py @@ -1,6 +1,7 @@ #!/usr/bin/env python import mkipp import kipp_data +import mesa_data import matplotlib.pyplot as plt from matplotlib.patches import PathPatch import numpy as np @@ -35,10 +36,9 @@ kipp_args = mkipp.Kipp_Args() fig = plt.figure() axis = plt.gca() -profile_names = ["LOGS/profile"+str(int(i))+".data" \ - for i in np.loadtxt("LOGS/profiles.index", skiprows = 1, usecols = (2,))] +profile_paths = mesa_data.get_profile_paths(["LOGS"]) #if data is distributed among several history.data files, you can provide them -history_names = ["LOGS/history.data"] +history_paths = ["LOGS/history.data"] #read profile data #kipp_data.get_xyz_data returns an object containing # xyz_data.xlims : limits of data in x coordinate @@ -46,7 +46,7 @@ # xyz_data.Y : 2D array of xaxis values of profile data # xyz_data.Z : 2D array of xaxis values of profile data # the last three can be used as inputs for matplotlib contour or contourf -xyz_data = kipp_data.get_xyz_data(profile_names, kipp_args.xaxis_divide, kipp_args) +xyz_data = kipp_data.get_xyz_data(profile_paths, kipp_args.xaxis_divide, kipp_args) #read mixing regions #kipp_data.get_mixing_zones returns an object containing # mixing_zones.zones : matplotlib Path objects for each mixing zone. @@ -56,7 +56,7 @@ # mixing_zones.y_coords : y coordinates for points at the surface # mixing_zones.histories : mesa_data history files to access additional data # the last three can be used as inputs for matplotlib contour or contourf -mixing_zones = kipp_data.get_mixing_zones(history_names, kipp_args.xaxis_divide, xyz_data.xlims, kipp_args) +mixing_zones = kipp_data.get_mixing_zones(history_paths, kipp_args.xaxis_divide, xyz_data.xlims, kipp_args) # just plot convection, overshooting and semiconvection for i,zone in enumerate(mixing_zones.zones): color = "" diff --git a/kipp_data.py b/kipp_data.py index b711f53..39d1d80 100644 --- a/kipp_data.py +++ b/kipp_data.py @@ -102,12 +102,12 @@ def merge_zone(self, zone2): self.new_blocks.extend(zone2.new_blocks) #returns a list of matplotlib path objects with the mixing regions -def get_mixing_zones(history_names, xaxis_divide, xlims, kipp_args): +def get_mixing_zones(history_paths, xaxis_divide, xlims, kipp_args): # Get mixing zones and mass borders from history.data files print "Reading history data" mix_data = [] histories = [] - for history_name in history_names: + for history_name in history_paths: history = Mesa_Data(history_name, read_data = False) columns = [] for key in history.columns.keys(): @@ -217,7 +217,7 @@ def get_mixing_zones(history_names, xaxis_divide, xlims, kipp_args): Mixing_Zones = namedtuple('Mixing_Zones', 'zones mix_types x_coords y_coords histories') return Mixing_Zones(zones, mix_types, x_coords, y_coords, histories) -def get_xyz_data(profile_names, xaxis_divide, kipp_args): +def get_xyz_data(profile_paths, xaxis_divide, kipp_args): # Extract data from profiles max_x_coord = -1e99 min_x_coord = 1e99 @@ -227,7 +227,7 @@ def get_xyz_data(profile_names, xaxis_divide, kipp_args): if kipp_args.yaxis_normalize: max_y = star_mass = star_radius = 1.0 else: - for i,profile_name in enumerate(profile_names): + for i,profile_name in enumerate(profile_paths): try: prof = Mesa_Data(profile_name, only_read_header = True) if kipp_args.yaxis == "mass": @@ -248,13 +248,13 @@ def get_xyz_data(profile_names, xaxis_divide, kipp_args): columns.extend(kipp_args.extractor(\ kipp_args.identifier, kipp_args.log10_on_data, prof, return_data_columns = True)) # Initialize interpolation grid - Z_data_array = np.zeros((kipp_args.yresolution,len(profile_names))) + Z_data_array = np.zeros((kipp_args.yresolution,len(profile_paths))) # XY coordinates for data - X_data_array = np.zeros((kipp_args.yresolution,len(profile_names))) - Y_data_array = np.zeros((kipp_args.yresolution,len(profile_names))) + X_data_array = np.zeros((kipp_args.yresolution,len(profile_paths))) + Y_data_array = np.zeros((kipp_args.yresolution,len(profile_paths))) for j in range(kipp_args.yresolution): Y_data_array[j,:] = max_y * j / (kipp_args.yresolution-1) - for i,profile_name in enumerate(profile_names): + for i,profile_name in enumerate(profile_paths): try: prof = Mesa_Data(profile_name, read_data = False) prof.read_data(columns) diff --git a/mesa_data.py b/mesa_data.py index 4e8ec28..a12a691 100644 --- a/mesa_data.py +++ b/mesa_data.py @@ -71,6 +71,26 @@ def read_data(self, column_names, clean_data = True): for column in column_names: self.data[column] = ma.masked_array(self.data[column], mask=mask).compressed() - def get(self,key): return self.data[key] + +#reads the profiles.index files in the folders specified by the logs_dirs array and returns +#an array containing paths to the individual profile files, after cleaning up redos and backups +def get_profile_paths(logs_dirs = ["LOGS"]): + profile_paths = [] + for log_dir in logs_dirs: + model_number, paths = np.loadtxt(log_dir+"/profiles.index", skiprows = 1, usecols = (0,2), unpack = True) + mask = np.zeros(len(paths)) + max_model_number = model_number[-1] + #last entry is valid, start from there and remove repeats + for i in range(len(model_number)-2,-1,-1): + if model_number[i] >= max_model_number: + mask[i] = 1 + else: + max_model_number = model_number[i] + + if sum(mask) > 0: + paths = ma.masked_array(paths, mask=mask).compressed() + profile_paths.extend([log_dir+"/profile"+str(int(i))+".data" for i in paths]) + return profile_paths + diff --git a/mkipp.py b/mkipp.py index 99915a0..a33501f 100755 --- a/mkipp.py +++ b/mkipp.py @@ -52,8 +52,8 @@ def default_extractor(identifier, log10_on_data, prof, return_data_columns = Fal class Kipp_Args: def __init__(self, logs_dirs = ['LOGS'], - profile_names = [], - history_names = [], + profile_paths = [], + history_paths = [], clean_data = True, extra_history_cols = [], identifier = "eps_nuc", @@ -83,10 +83,10 @@ def __init__(self, All arguments are optional, if not provided defaults are assigned Args: - logs_dir (List[str]): List of paths to MESA LOGS directories. If profile_names and - history_names are not provided, they are automatically generated from logs_dir. - profile_names (List[str]): List of paths to MESA profile files. - history_names (List[str]): List of paths to MESA history files. + logs_dir (List[str]): List of paths to MESA LOGS directories. If profile_paths and + history_paths are not provided, they are automatically generated from logs_dir. + profile_paths (List[str]): List of paths to MESA profile files. + history_paths (List[str]): List of paths to MESA history files. clean_data (bool): Clean history files in case of redos. If data has no redos, a small performance gain can be had by setting this to False. extra_history_cols (List[str]): Additional column names to be extracted from history files. @@ -124,8 +124,8 @@ def __init__(self, """ self.logs_dirs = logs_dirs - self.profile_names = profile_names - self.history_names = history_names + self.profile_paths = profile_paths + self.history_paths = history_paths self.clean_data = clean_data self.extra_history_cols = extra_history_cols self.identifier = identifier @@ -163,17 +163,14 @@ def kipp_plot(kipp_args, axis=None): axis = fig.gca() #Fill profile and history names if unspecified - profile_names = kipp_args.profile_names - if len(profile_names) == 0: - profile_names = [] + profile_paths = kipp_args.profile_paths + if len(profile_paths) == 0: + profile_paths = get_profile_paths(logs_dirs = kipp_args.logs_dirs) + history_paths = kipp_args.history_paths + if len(history_paths) == 0: + history_paths = [] for log_dir in kipp_args.logs_dirs: - profile_names.extend(\ - [log_dir+"/profile"+str(int(i))+".data" for i in np.loadtxt(log_dir+"/profiles.index", skiprows = 1, usecols = (2,))]) - history_names = kipp_args.history_names - if len(history_names) == 0: - history_names = [] - for log_dir in kipp_args.logs_dirs: - history_names.append(log_dir+"/history.data") + history_paths.append(log_dir+"/history.data") xaxis_divide = kipp_args.xaxis_divide if kipp_args.xaxis == "star_age": @@ -184,7 +181,7 @@ def kipp_plot(kipp_args, axis=None): elif kipp_args.time_units == "Gyr": xaxis_divide = 1e9 - xyz_data = get_xyz_data(profile_names, xaxis_divide, kipp_args) + xyz_data = get_xyz_data(profile_paths, xaxis_divide, kipp_args) #Get levels if undefined levels = kipp_args.levels @@ -200,8 +197,8 @@ def kipp_plot(kipp_args, axis=None): cmap=kipp_args.contour_colormap, levels=levels, antialiased = False) #zones, mix_types, x_coords, y_coords, histories = get_mixing_zones(\ - # history_names, xaxis_divide, xyz_data.xlims, kipp_args) - mixing_zones = get_mixing_zones(history_names, xaxis_divide, xyz_data.xlims, kipp_args) + # history_paths, xaxis_divide, xyz_data.xlims, kipp_args) + mixing_zones = get_mixing_zones(history_paths, xaxis_divide, xyz_data.xlims, kipp_args) for i,zone in enumerate(mixing_zones.zones):