• Facebook
  • Twitter
  • Reddit
  • StumbleUpon
  • Digg
  • email

#!/usr/bin/env python
 
# plot_flame_dist.py
 
# By Nathan C. Hearn
#    September 24, 2008
#
# Plotting utility for flame_dist_analysis output
#
# Requires QuickFlash/utils to be in PYTHONPATH environment variable
 
 
try :
    import psyco
except ImportError :
    pass
 
 
data_start_keyword = "DistribData"
default_filename_extension = "png"
 
png_filename_extension = "png"
png_matplotlib_backend = "Agg"
 
svg_filename_extension = "svg"
svg_matplotlib_backend = "SVG"
 
ps_filename_extension = "ps"
ps_matplotlib_backend = "PS"
 
pdf_filename_extension = "pdf"
pdf_matplotlib_backend = "PDF"
 
mean_value_suffix = "_Mean_Value"
dispersion_suffix = "_Dispersion"
rms_suffix = "_RMS_Difference"
 
density_suffix = "_Density"  ## Initially, a double-underscore was used
area_suffix = "_Area"
 
 
def read_file(filename) :
 
    from plot_gen import read_keywords_table
 
    ## Read file data and attributes
 
    keyword_values, table_dict \
        = read_keywords_table(filename, table_start_line=data_start_keyword,
                              ignore_leading_comment_mark=True)
 
    return keyword_values, table_dict
 
 
def build_plot(output_file, matplotlib_backend, plot_title, x_axis_title,
               min_x, max_x, y_axis_title, pos_array_list, data_array_list,
               legend_list) :
 
    from plot_gen import Plot
 
    num_lists = len(pos_array_list)
 
    if (len(data_array_list) != num_lists) :
        raise RuntimeError("List count mismatch")
 
    if (len(legend_list) != num_lists) :
        raise RuntimeError("List count mismatch")
 
    p = Plot(plot_title, x_axis_title, y_axis_title)
 
    data_bounds_set = False
 
    min_y = None
    max_y = None
 
    for list_index in range(num_lists) :
 
        data_array = data_array_list[list_index]
 
        local_min_y = min(data_array)
        local_max_y = max(data_array)
 
        if (not data_bounds_set) :
 
            min_y = local_min_y
            max_y = local_max_y
 
            data_bounds_set = True
 
        else :
 
            if (local_min_y < min_y) :
                min_y = local_min_y
 
            if (local_max_y > max_y) :
                max_y = local_max_y
 
        p.add_plot(pos_array_list[list_index], data_array,
                   legend_text=legend_list[list_index])
 
    p.set_limits(min_x=min_x, max_x=max_x, min_y=min_y, max_y=max_y)
 
    p.show_plot(output_file, matplotlib_backend=matplotlib_backend)
 
 
def scale_data(data_array, scale_factor) :
 
    new_array = list()
 
    for elem in data_array :
        new_array.append(elem * scale_factor)
 
    return new_array
 
 
def offset_data(data_array, offset_value) :
 
    new_array = list()
 
    for elem in data_array :
        new_array.append(elem + offset_value)
 
    return new_array
 
 
def main() :
 
    from sys import exit, stderr
    from os import path
    from optparse import OptionParser
 
 
    ## Set up the command line parser
 
    usage = str("%prog [ options ] output_prefix input_file1[,legend_1] ...")
    version = str("%prog 0.1")
 
    argparser = OptionParser(usage, version=version)
 
    output_extension_str = str("output-extension")
    argparser.add_option("--" + output_extension_str, type="string",
                         dest="output_extension",
                         help="Filename extension for output files")
    argparser.set_defaults(output_format=None)
 
    plot_title_str = str("plot-title")
    argparser.add_option("--" + plot_title_str, type="string",
                         dest="plot_title",
                         help="Text for plot titles")
    argparser.set_defaults(plot_title=None)
 
    time_legend_str = str("time-legend")
    argparser.add_option("--" + time_legend_str, action="store_true",
                         dest="time_legend",
                         help="Add time to plot legend")
    argparser.set_defaults(time_legend=False)
 
    center_flame_str = str("center-flame")
    argparser.add_option("--" + center_flame_str, type="string",
                         dest="center_flame",
                         help="Center flame (choose mass or mass-fraction)")
    argparser.set_defaults(center_flame=None)
 
    rescale_pos_str = str("position-rescale")
    argparser.add_option("--" + rescale_pos_str, type="string",
                         dest="position_rescale",
                         help="Comma-separated scale factor and unit name")
    argparser.set_defaults(position_rescale=None)
 
    opts, args = argparser.parse_args()
 
    argc = len(args)
 
    min_args = int(2)
 
    if (argc < min_args) :
        stderr.write("\n")
        argparser.print_help(file=stderr)
        stderr.write("\n")
        exit(1)
 
 
    ## Get the required arguments
 
    argptr = int(0)
 
    user_output_prefix = str(args[argptr])
    argptr += 1
 
    file_list = list()
    legend_list = list()
 
    while (argptr < argc) :
 
        file_legend_entry = str(args[argptr])
        argptr += 1
 
        file_legend_list = file_legend_entry.split(",")
 
        num_words = len(file_legend_list)
 
        if (num_words < 1) :
            stderr.write("\nError: Unable to parse data file list\n\n")
            exit(-1)
 
        filename = file_legend_list[0]
 
        legend_text = None
 
        if (num_words < 2) :
            legend_text_end = filename.rfind("_hdf5")
 
            if (legend_text_end < 1) :
                legend_text_end = len(filename)
 
            legend_text = filename[0:legend_text_end]
 
        elif (num_words == 2) :
 
            legend_text = file_legend_list[1]
 
        else :
 
            stderr.write("\nError: Improper format for filename-legend pair [ "
                         + file_legend_entry + " ]\n\n")
            exit(-1)
 
        file_list.append(filename)
        legend_list.append(legend_text)
 
    num_files = len(file_list)
 
 
    ## Read the optional arguments
 
    output_extension = opts.output_extension
 
    if (output_extension is None) :
 
        output_extension = default_filename_extension
 
    else :
 
        if (output_extension == png_filename_extension) :
 
            matplotlib_backend = png_matplotlib_backend
 
        elif (output_extension == svg_filename_extension) :
 
            matplotlib_backend = svg_matplotlib_backend
 
        elif (output_extension == ps_filename_extension) :
 
            matplotlib_backend = ps_matplotlib_backend
 
        elif (output_extension == pdf_filename_extension) :
 
            matplotlib_backend = pdf_matplotlib_backend
 
        else :
 
            stderr.write("\nError: Output extension [ " + output_extension
                         + " ] not recognized\n\n")
            exit(-1)
 
    user_plot_title = opts.plot_title
 
    plot_title_prefix = str()
 
    if (user_plot_title is not None) :
        plot_title_prefix = user_plot_title + " - "
 
    add_time_legend = opts.time_legend
 
    center_flame_method = opts.center_flame
 
    recenter_flame = False
 
    if (center_flame_method is not None) :
 
        if (not ((center_flame_method == "mass")
                 or (center_flame_method == "mass-fraction"))) :
            stderr.write("\nError: Flame centering method [ "
                         + center_flame_method + " ] not recognized\n\n")
            exit(-1)
 
        recenter_flame = True
 
    position_scale = 1.0
    title_unit_suffix = " (cm)"
 
    if (opts.position_rescale is not None) :
 
        words = opts.position_rescale.split(",")
 
        num_words = len(words)
 
        if ((num_words < 1) or (num_words > 2)) :
 
            stderr.write("\nError: Position scale requires comma-separated "
                         + "scale factor and (optional) unit name\n\n")
            exit(-1)
 
        position_scale = float(words[0])
 
        if (num_words > 1) :
            title_unit_suffix = " (" + words[1] + ")"
        else :
            title_unit_suffix = str()
 
    abs_pos_x_title = "Position" + title_unit_suffix
 
    ## Get other quantities
 
    output_prefix = user_output_prefix
 
    if (len(output_prefix) < 1) :
        stderr.write("\nError: Improper output_prefix\n\n")
        exit(-1)
 
    if (path.split(output_prefix)[1] != "") :
        if (output_prefix[-1] != "_") :
            output_prefix += "_"
 
    output_suffix = "." + output_extension
 
    ## Read the files
 
    table_list = list()
 
    sim_times_list = list()
 
    flame_center_mass_list = list()
    flame_center_mass_fraction_list = list()
 
    area_fractions_present_list = list()  ## booleans
    area_min_flame_value_list = list()
    area_min_flame_value_str_list = list()
    area_max_flame_values_list = list()
    area_max_flame_values_str_list = list()
 
    flame_min_pos = None
    flame_max_pos = None
 
    bin_array_min_pos = None
    bin_array_max_pos = None
 
    user_dataset_list = None
 
    rms_ref_value = None
 
    time_data_filename = output_prefix + "time_data.csv"
 
    time_data_file = open(time_data_filename, "w")
 
    time_data_file.write("Filename" + "\t" + "Sim_Time" + "\t"
                         + "Flame_Center_Mass_Fraction" + "\t"
                         + "Flame_Center_Mass" + "\t" + "Min_Flame_Pos"
                         + "\t" + "Max_Flame_Pos" + "\n")
 
    for filename in file_list :
 
        keyword_value, table = read_file(filename)
 
        sim_time = float(keyword_value.get_first_value("SimTime"))
 
        raw_flame_center_mass \
            = float(keyword_value.get_first_value("FlameCenterMass"))
 
        flame_center_mass \
            = raw_flame_center_mass * position_scale
 
        raw_flame_center_mass_fraction \
            = float(keyword_value.get_first_value("FlameCenterMassFract"))
 
        flame_center_mass_fraction \
            = raw_flame_center_mass_fraction * position_scale
 
        raw_local_flame_min_pos \
            = float(keyword_value.get_first_value("MinFlamePos"))
 
        local_flame_min_pos \
            = raw_local_flame_min_pos * position_scale
 
        raw_local_flame_max_pos \
            = float(keyword_value.get_first_value("MaxFlamePos"))
 
        local_flame_max_pos \
            = raw_local_flame_max_pos * position_scale
 
        time_data_file.write(filename + "\t" + str(sim_time) + "\t"
                             + str(raw_flame_center_mass_fraction) + "\t"
                             + str(raw_flame_center_mass) + "\t"
                             + str(raw_local_flame_min_pos) + "\t"
                             + str(raw_local_flame_max_pos) + "\n")
 
        if (flame_min_pos is None) :
 
            flame_min_pos = local_flame_min_pos
            flame_max_pos = local_flame_max_pos
 
        else :
 
            if (local_flame_min_pos < flame_min_pos) :
                flame_min_pos = local_flame_min_pos
 
            if (local_flame_max_pos > flame_max_pos) :
                flame_max_pos = local_flame_max_pos
 
        local_bin_array_min_pos \
            = float(keyword_value.get_first_value("BinArrayMinPos")) \
              * position_scale
 
        local_bin_array_max_pos \
            = float(keyword_value.get_first_value("BinArrayMaxPos")) \
              * position_scale
 
        if (bin_array_min_pos is None) :
 
            bin_array_min_pos = local_bin_array_min_pos
            bin_array_max_pos = local_bin_array_max_pos
 
        else :
 
            if (local_bin_array_min_pos < bin_array_min_pos) :
                bin_array_min_pos = local_bin_array_min_pos
 
            if (local_bin_array_max_pos > bin_array_max_pos) :
                bin_array_max_pos = local_bin_array_max_pos
 
        distrib_datasets_index \
            = keyword_value.find_first_pair("DistribDatasets")
 
        if (distrib_datasets_index is not None) :
 
            local_dataset_str \
                = keyword_value[distrib_datasets_index].get_value()
 
            local_dataset_list = local_dataset_str.split(",")
 
            local_dataset_list.sort()
 
            if (user_dataset_list is None) :
                user_dataset_list = local_dataset_list
            elif (local_dataset_list != user_dataset_list) :
                stderr.write("\nError: Incompatible user dataset lists\n\n")
                exit(-2)
 
        elif (user_dataset_list is None) :
 
            user_dataset_list = list()
 
        elif (len(user_dataset_list) > 0) :
 
            stderr.write("\nError: Incompatible user dataset lists\n\n")
 
        rms_diff_index \
            = keyword_value.find_first_pair("RootMeanSqDiffRefValue")
 
        if (rms_diff_index is not None) :
 
            local_ref_value = float(keyword_value[rms_diff_index].get_value())
 
            if (rms_ref_value is None) :
                rms_ref_value = local_ref_value
            elif (local_ref_value != rms_ref_value) :
                stderr.write("\nError: Incompatible RMS reference value\n\n")
 
        table_list.append(table)
 
        sim_times_list.append(sim_time)
        flame_center_mass_list.append(flame_center_mass)
        flame_center_mass_fraction_list.append(flame_center_mass_fraction)
 
        area_flame_min_index \
            = keyword_value.find_first_pair("AreaFlameMinvalue")
 
        if (area_flame_min_index is not None) :
 
            area_fractions_present_list.append(True)
 
            min_flame_value_str = keyword_value[area_flame_min_index]
 
            max_flame_values_str \
                = keyword_value.get_first_value("AreaFlameMaxValues")
 
            max_flame_values_str_list = max_flame_values_str.split(",")
 
            max_flame_values_str_list.sort()
 
            max_flame_values = list()
 
            for value_string in max_flame_values_str_list :
                max_flame_values.append(float(value_string))
 
            area_min_flame_value_list.append(float(min_flame_value_str))
            area_min_flame_value_str_list.append(min_flame_value_str)
 
            area_max_flame_values_list.append(max_flame_values)
            area_max_flame_values_str_list.append(max_flame_values_str)
 
        else :
 
            area_fractions_present_list.append(False)
 
            area_min_flame_value_list.append(None)
            area_min_flame_value_str_list.append(None)
 
            area_max_flame_values_list.append(None)
            area_max_flame_values_str_list.append(None)
 
 
    if (add_time_legend) :
 
        for file_index in range(num_files) :
 
            sim_time = sim_times_list[file_index]
 
            new_legend = legend_list[file_index] + " " \
                         + ("%(time)3.2f" % { "time" : sim_time })
 
            legend_list[file_index] = new_legend
 
    ## Get the position lists
 
    pos_array_list = list()
    pos_array_center_mass_list = list()
    pos_array_center_massfract_list = list()
 
    for file_index in range(num_files) :
 
        orig_pos_array = scale_data(table_list[file_index]["Bin_Center"],
                                    position_scale)
 
        center_mass = flame_center_mass_list[file_index]
        center_massfract = flame_center_mass_fraction_list[file_index]
 
        center_mass_offset = -1.0 * center_mass
        center_massfract_offset = -1.0 * center_massfract
 
        pos_array_list.append(orig_pos_array)
 
        pos_array_mass = offset_data(orig_pos_array, center_mass_offset)
 
        pos_array_center_mass_list.append(pos_array_mass)
 
        pos_array_massfract = offset_data(orig_pos_array,
                                          center_massfract_offset)
 
        pos_array_center_massfract_list.append(pos_array_massfract)
 
 
    ## Density distribution
 
    dens_output_file = output_prefix + "density" + output_suffix
 
    dens_data_list = list()
 
    for data_table in table_list :
        dens_data_list.append(data_table["Mean_Density"])
 
    dens_plot_title = plot_title_prefix + "Mean Density"
 
    dens_plot_y_title = "Density (g / cc)"
 
    build_plot(dens_output_file, matplotlib_backend, dens_plot_title,
               abs_pos_x_title, flame_min_pos, flame_max_pos,
               dens_plot_y_title, pos_array_list, dens_data_list, legend_list)
 
    ## Burned mass fraction
 
    burned_massfract_output_file = output_prefix + "burned_mass_fraction" \
                                   + output_suffix
 
    burned_massfract_list = list()
 
    for data_table in table_list :
        burned_massfract_list.append(data_table["Burned_Mass_Fraction"])
 
    burned_massfract_title = plot_title_prefix + "Burned Mass Fraction"
 
    burned_massfract_y_title = "Mass Fraction"
 
    build_plot(burned_massfract_output_file, matplotlib_backend,
               burned_massfract_title, abs_pos_x_title, flame_min_pos,
               flame_max_pos, burned_massfract_y_title, pos_array_list,
               burned_massfract_list, legend_list)
 
    ## User dataset distribution
 
    if (user_dataset_list is not None) :
        for user_dataset in user_dataset_list :
 
            ## Mean value
 
            mean_output_file = output_prefix + user_dataset + "_mean" \
                               + output_suffix
 
            mean_name = user_dataset + mean_value_suffix
 
            mean_value_list = list()
 
            for data_table in table_list :
                mean_value_list.append(data_table[mean_name])
 
            mean_title = plot_title_prefix + "Mean Value"
 
            mean_y_title = user_dataset + " Mean Value"
 
            build_plot(mean_output_file, matplotlib_backend, mean_title,
                       abs_pos_x_title, bin_array_min_pos, bin_array_max_pos,
                       mean_y_title, pos_array_list, mean_value_list,
                       legend_list)
 
            ## Dispersion
 
            dispersion_output_file = output_prefix + user_dataset \
                                     + "_dispersion" + output_suffix
 
            dispersion_name = user_dataset + dispersion_suffix
 
            dispersion_list = list()
 
            for data_table in table_list :
                dispersion_list.append(data_table[dispersion_name])
 
            dispersion_title = plot_title_prefix + "Dispersion"
 
            dispersion_y_title = user_dataset + " Dispersion"
 
            build_plot(dispersion_output_file, matplotlib_backend,
                       dispersion_title, abs_pos_x_title, bin_array_min_pos,
                       bin_array_max_pos, dispersion_y_title, pos_array_list,
                       dispersion_list, legend_list)
 
            ## RMS Value
 
            rms_output_file = output_prefix + user_dataset + "_rms_diff" \
                              + output_suffix
 
            rms_name = user_dataset + rms_suffix
 
            rms_list = list()
 
            for data_table in table_list :
                rms_list.append(data_table[rms_name])
 
            rms_title = plot_title_prefix
 
            rms_y_title = user_dataset
 
            if (rms_ref_value is None) :
                rms_title += "Population Standard Deviation"
                rms_y_title += " Population Standard Deviation"
            else :
                rms_title += "RMS Difference from " + str(rms_ref_value)
                rms_y_title += " RMS Difference from " + str(rms_ref_value)
 
            build_plot(rms_output_file, matplotlib_backend,
                       rms_title, abs_pos_x_title, bin_array_min_pos,
                       bin_array_max_pos, rms_y_title, pos_array_list,
                       rms_list, legend_list)
 
    ## Flame area distribution
 
    for file_index in range(num_files) :
        if (area_fractions_present_list[file_index]) :
 
            pass
 
 
if (__name__ == "__main__") :
    main()