|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +""" |
| 4 | +Compute a box size correction file that can be used as the 'box_size_correction' |
| 5 | +argument for an autoplotter plot. |
| 6 | +
|
| 7 | +Usage: |
| 8 | + velociraptor-compute-box-size-correction \ |
| 9 | + smallbox largebox plotname plottype output |
| 10 | +
|
| 11 | +with: |
| 12 | + - smallbox/largebox: data*.yml output file from a pipeline run |
| 13 | + - plotname: Name of a particular plot in the data*.yml files |
| 14 | + - plottype: Type of plot (currently supported: mass_function) |
| 15 | + - output: Name of an output .yml file. If the .yml extension is missing, it is |
| 16 | + added. |
| 17 | +""" |
| 18 | + |
| 19 | +import argparse |
| 20 | +import os |
| 21 | +import yaml |
| 22 | +import numpy as np |
| 23 | +import scipy.interpolate as interpol |
| 24 | + |
| 25 | + |
| 26 | +def velociraptor_compute_box_size_correction(): |
| 27 | + argparser = argparse.ArgumentParser("Compute the box size correction for a plot.") |
| 28 | + argparser.add_argument( |
| 29 | + "smallbox", help="Pipeline output for the small box that needs to be corrected." |
| 30 | + ) |
| 31 | + argparser.add_argument( |
| 32 | + "largebox", help="Pipeline output for the large box that we want to correct to." |
| 33 | + ) |
| 34 | + argparser.add_argument("plotname", help="Name of the plot that we want to correct.") |
| 35 | + argparser.add_argument("plottype", help="Type of the plot we want to correct.") |
| 36 | + argparser.add_argument( |
| 37 | + "output", help="Name of the output file that will store the correction." |
| 38 | + ) |
| 39 | + args = argparser.parse_args() |
| 40 | + |
| 41 | + if not args.plottype in ["mass_function"]: |
| 42 | + raise AttributeError( |
| 43 | + f"Cannot compute box size correction for plot type {args.plottype}!" |
| 44 | + ) |
| 45 | + |
| 46 | + log_x = False |
| 47 | + log_y = False |
| 48 | + if args.plottype in ["mass_function"]: |
| 49 | + log_x = True |
| 50 | + log_y = True |
| 51 | + |
| 52 | + small_box = args.smallbox |
| 53 | + large_box = args.largebox |
| 54 | + for file in [args.smallbox, args.largebox]: |
| 55 | + if not os.path.exists(file): |
| 56 | + raise AttributeError(f"File {file} could not be found!") |
| 57 | + |
| 58 | + output_file = args.output |
| 59 | + if not output_file.endswith(".yml"): |
| 60 | + output_file += ".yml" |
| 61 | + try: |
| 62 | + open(output_file, "w").close() |
| 63 | + except: |
| 64 | + raise AttributeError(f"Can not write to {output_file}!") |
| 65 | + |
| 66 | + with open(args.smallbox, "r") as handle: |
| 67 | + small_box = yaml.safe_load(handle) |
| 68 | + with open(args.largebox, "r") as handle: |
| 69 | + large_box = yaml.safe_load(handle) |
| 70 | + |
| 71 | + try: |
| 72 | + small_box_data = small_box[args.plotname]["lines"] |
| 73 | + except: |
| 74 | + raise AttributeError(f"Could not find {args.plotname} in {args.smallbox}!") |
| 75 | + try: |
| 76 | + large_box_data = large_box[args.plotname]["lines"] |
| 77 | + except: |
| 78 | + raise AttributeError(f"Could not find {args.plotname} in {args.largebox}!") |
| 79 | + |
| 80 | + try: |
| 81 | + small_box_plot_data = small_box_data[args.plottype] |
| 82 | + except: |
| 83 | + raise AttributeError( |
| 84 | + f"{args.plottype} not found in plot {args.plotname} in {args.smallbox}!" |
| 85 | + ) |
| 86 | + try: |
| 87 | + large_box_plot_data = large_box_data[args.plottype] |
| 88 | + except: |
| 89 | + raise AttributeError( |
| 90 | + f"{args.plottype} not found in plot {args.plotname} in {args.largebox}!" |
| 91 | + ) |
| 92 | + |
| 93 | + small_box_x = small_box_plot_data["centers"] |
| 94 | + small_box_y = small_box_plot_data["values"] |
| 95 | + large_box_x = large_box_plot_data["centers"] |
| 96 | + large_box_y = large_box_plot_data["values"] |
| 97 | + |
| 98 | + if log_x: |
| 99 | + small_box_x = np.log10(small_box_x) |
| 100 | + large_box_x = np.log10(large_box_x) |
| 101 | + |
| 102 | + if log_y: |
| 103 | + small_box_y = np.log10(small_box_y) |
| 104 | + large_box_y = np.log10(large_box_y) |
| 105 | + |
| 106 | + small_spline = interpol.InterpolatedUnivariateSpline(small_box_x, small_box_y) |
| 107 | + large_spline = interpol.InterpolatedUnivariateSpline(large_box_x, large_box_y) |
| 108 | + |
| 109 | + xmin = max(small_box_x.min(), large_box_x.min()) |
| 110 | + xmax = min(small_box_x.max(), large_box_x.max()) |
| 111 | + x_range = np.linspace(xmin, xmax, 100) |
| 112 | + small_y_range = small_spline(x_range) |
| 113 | + large_y_range = large_spline(x_range) |
| 114 | + |
| 115 | + if log_y: |
| 116 | + small_y_range = 10.0 ** small_y_range |
| 117 | + large_y_range = 10.0 ** large_y_range |
| 118 | + |
| 119 | + correction = large_y_range / small_y_range |
| 120 | + |
| 121 | + correction_data = {} |
| 122 | + correction_data["plot_name"] = args.plotname |
| 123 | + correction_data["plot_type"] = args.plottype |
| 124 | + correction_data["is_log_x"] = True |
| 125 | + correction_data["x_units"] = small_box_plot_data["centers_units"] |
| 126 | + correction_data["x_limits"] = np.array([xmin, xmax]).tolist() |
| 127 | + correction_data["x"] = x_range.tolist() |
| 128 | + correction_data["y"] = correction.tolist() |
| 129 | + with open(output_file, "w") as handle: |
| 130 | + yaml.safe_dump(correction_data, handle) |
| 131 | + |
| 132 | + |
| 133 | +if __name__ == "__main__": |
| 134 | + velociraptor_compute_box_size_correction() |
0 commit comments