|
1 | 1 | """ |
2 | | -Histogram aggregation with HistogramsAggregator |
| 2 | +Histogram aggregation with HistogramAggregator |
3 | 3 | =============================================== |
4 | 4 |
|
5 | 5 | This tutorial shows how to: |
6 | 6 |
|
7 | 7 | 1. Build an event table with camera-like data (images and peak times) and some invalid values. |
8 | | -2. Configure and run HistogramsAggregator in chunks. |
9 | | -3. Access counts, bin edges, and valid-event counts (n_events). |
| 8 | +2. Configure and run HistogramAggregator in chunks. |
| 9 | +3. Access histogram counts, bin edges, summary statistics, and valid-event counts (n_events). |
10 | 10 | 4. Plot one pixel histogram from the selected chunks and both gain channels for both image and peak_time columns. |
| 11 | +5. Overlay mean, median, and std on top of the histogram curves. |
11 | 12 | """ |
12 | 13 |
|
13 | 14 | import matplotlib.pyplot as plt |
14 | 15 | import numpy as np |
15 | | -import hist |
| 16 | +from matplotlib.lines import Line2D |
| 17 | +from matplotlib.patches import Patch |
16 | 18 | from astropy.table import Table |
17 | 19 | from astropy.time import Time |
18 | 20 | from traitlets.config import Config |
19 | 21 |
|
20 | | -from ctapipe.monitoring.aggregator import HistogramsAggregator |
| 22 | +from ctapipe.monitoring.aggregator import HistogramAggregator |
21 | 23 |
|
22 | 24 |
|
23 | 25 | # ------------------------------------------------------------------- |
|
36 | 38 | ) |
37 | 39 | event_ids = np.arange(n_events) |
38 | 40 | images = rng.normal(loc=77.0, scale=10.0, size=(n_events, n_channels, n_pixels)) |
39 | | -peak_time = rng.normal(loc=20.0, scale=5.0, size=(n_events, n_channels, n_pixels)) |
| 41 | +peak_time = rng.normal(loc=20.0, scale=2.0, size=(n_events, n_channels, n_pixels)) |
40 | 42 |
|
41 | 43 | # Add a few invalid values to demonstrate n_events behavior. |
42 | 44 | images[3, 0, 10] = np.nan |
|
60 | 62 | # ------------------------------------------------------------------- |
61 | 63 | config_image = Config( |
62 | 64 | { |
63 | | - "HistogramsAggregator": { |
| 65 | + "HistogramAggregator": { |
64 | 66 | "chunking_type": "SizeChunking", |
| 67 | + "hist_axis_dict": { |
| 68 | + "axis_class_name": "Regular", |
| 69 | + "kwargs": { |
| 70 | + "bins": 50, |
| 71 | + "start": 40.0, |
| 72 | + "stop": 110.0, |
| 73 | + "name": "value", |
| 74 | + }, |
| 75 | + }, |
65 | 76 | }, |
66 | 77 | "SizeChunking": {"chunk_size": 1000}, |
67 | 78 | } |
68 | 79 | ) |
69 | 80 |
|
70 | | -aggregator_image = HistogramsAggregator( |
71 | | - hist.axis.Regular(50, 40.0, 110.0, name="value"), |
72 | | - config=config_image, |
73 | | -) |
| 81 | +aggregator_image = HistogramAggregator(config=config_image) |
74 | 82 | result = aggregator_image( |
75 | 83 | table=table, |
76 | 84 | col_name="image", |
|
79 | 87 |
|
80 | 88 | config_peak_time = Config( |
81 | 89 | { |
82 | | - "HistogramsAggregator": { |
| 90 | + "HistogramAggregator": { |
83 | 91 | "chunking_type": "SizeChunking", |
| 92 | + "hist_axis_dict": { |
| 93 | + "axis_class_name": "Regular", |
| 94 | + "kwargs": { |
| 95 | + "bins": 50, |
| 96 | + "start": 2.0, |
| 97 | + "stop": 38.0, |
| 98 | + "name": "value", |
| 99 | + }, |
| 100 | + }, |
84 | 101 | }, |
85 | 102 | "SizeChunking": {"chunk_size": 1000}, |
86 | 103 | } |
87 | 104 | ) |
88 | 105 |
|
89 | | -aggregator_peak_time = HistogramsAggregator( |
90 | | - hist.axis.Regular(50, 2.0, 38.0, name="value"), |
91 | | - config=config_peak_time, |
92 | | -) |
| 106 | +aggregator_peak_time = HistogramAggregator(config=config_peak_time) |
93 | 107 | result_peak_time = aggregator_peak_time( |
94 | 108 | table=table, |
95 | 109 | col_name="peak_time", |
96 | 110 | masked_elements_of_sample=masked_elements_of_sample, |
97 | 111 | ) |
98 | 112 |
|
99 | 113 | print(f"Number of chunks: {len(result)}") |
100 | | -print(f"counts shape per chunk: {result[0]['counts'].shape}") |
101 | | -print(f"edges shape per chunk: {result[0]['edges'].shape}") |
| 114 | +print(f"histogram shape per chunk: {result[0]['histogram'].shape}") |
| 115 | +print(f"edges shape per chunk: {result[0]['meta']['bin_edges'].shape}") |
102 | 116 | print(f"n_events shape per chunk: {result[0]['n_events'].shape}") |
103 | 117 |
|
104 | 118 |
|
|
110 | 124 |
|
111 | 125 | fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True) |
112 | 126 | for chunk_index, ax in enumerate(axes): |
113 | | - edges = result[chunk_index]["edges"] |
| 127 | + edges = result[chunk_index]["meta"]["bin_edges"] |
| 128 | + channel_handles = [] |
114 | 129 |
|
115 | 130 | for channel_index in range(n_channels): |
116 | | - counts = result[chunk_index]["counts"][:, channel_index, pixel_index] |
| 131 | + counts = result[chunk_index]["histogram"][:, channel_index, pixel_index] |
117 | 132 | valid_events = result[chunk_index]["n_events"][channel_index, pixel_index] |
118 | | - ax.step( |
| 133 | + mean_val = result[chunk_index]["mean"][channel_index, pixel_index] |
| 134 | + median_val = result[chunk_index]["median"][channel_index, pixel_index] |
| 135 | + std_val = result[chunk_index]["std"][channel_index, pixel_index] |
| 136 | + |
| 137 | + line = ax.step( |
119 | 138 | edges[:-1], |
120 | 139 | counts, |
121 | 140 | where="post", |
122 | 141 | label=f"{gain_label[channel_index]} (n_events={valid_events})", |
| 142 | + )[0] |
| 143 | + channel_handles.append(line) |
| 144 | + color = line.get_color() |
| 145 | + |
| 146 | + ax.axvline(mean_val, color=color, linestyle="--", linewidth=1.2) |
| 147 | + ax.axvline(median_val, color=color, linestyle=":", linewidth=1.2) |
| 148 | + ax.axvspan( |
| 149 | + mean_val - std_val, |
| 150 | + mean_val + std_val, |
| 151 | + color=color, |
| 152 | + alpha=0.12, |
123 | 153 | ) |
124 | 154 |
|
125 | 155 | ax.set_title(f"Chunk {chunk_index}, pixel {pixel_index}") |
126 | 156 | ax.set_xlabel("image value") |
127 | 157 | ax.set_ylabel("Counts") |
128 | | - ax.legend(loc="upper right", fontsize=8) |
| 158 | + stat_handles = [ |
| 159 | + Line2D([0], [0], color="black", linestyle="--", linewidth=1.2, label="Mean"), |
| 160 | + Line2D([0], [0], color="black", linestyle=":", linewidth=1.2, label="Median"), |
| 161 | + Patch(facecolor="gray", alpha=0.12, label="Mean ± Std"), |
| 162 | + ] |
| 163 | + ax.legend(handles=channel_handles + stat_handles, loc="upper left", fontsize=8) |
129 | 164 |
|
130 | 165 | plt.show() |
131 | 166 |
|
|
135 | 170 | # ------------------------------------------------------------------- |
136 | 171 | fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True) |
137 | 172 | for chunk_index, ax in enumerate(axes): |
138 | | - edges = result_peak_time[chunk_index]["edges"] |
| 173 | + edges = result_peak_time[chunk_index]["meta"]["bin_edges"] |
| 174 | + channel_handles = [] |
139 | 175 |
|
140 | 176 | for channel_index in range(n_channels): |
141 | | - counts = result_peak_time[chunk_index]["counts"][:, channel_index, pixel_index] |
| 177 | + counts = result_peak_time[chunk_index]["histogram"][ |
| 178 | + :, channel_index, pixel_index |
| 179 | + ] |
142 | 180 | valid_events = result_peak_time[chunk_index]["n_events"][ |
143 | 181 | channel_index, pixel_index |
144 | 182 | ] |
145 | | - ax.step( |
| 183 | + mean_val = result_peak_time[chunk_index]["mean"][channel_index, pixel_index] |
| 184 | + median_val = result_peak_time[chunk_index]["median"][channel_index, pixel_index] |
| 185 | + std_val = result_peak_time[chunk_index]["std"][channel_index, pixel_index] |
| 186 | + |
| 187 | + line = ax.step( |
146 | 188 | edges[:-1], |
147 | 189 | counts, |
148 | 190 | where="post", |
149 | 191 | label=f"{gain_label[channel_index]} (n_events={valid_events})", |
| 192 | + )[0] |
| 193 | + channel_handles.append(line) |
| 194 | + color = line.get_color() |
| 195 | + |
| 196 | + ax.axvline(mean_val, color=color, linestyle="--", linewidth=1.2) |
| 197 | + ax.axvline(median_val, color=color, linestyle=":", linewidth=1.2) |
| 198 | + ax.axvspan( |
| 199 | + mean_val - std_val, |
| 200 | + mean_val + std_val, |
| 201 | + color=color, |
| 202 | + alpha=0.12, |
150 | 203 | ) |
151 | 204 |
|
152 | 205 | ax.set_title(f"Peak Time - Chunk {chunk_index}, pixel {pixel_index}") |
153 | 206 | ax.set_xlabel("peak_time value") |
154 | 207 | ax.set_ylabel("Counts") |
155 | | - ax.legend(loc="upper right", fontsize=8) |
| 208 | + stat_handles = [ |
| 209 | + Line2D([0], [0], color="black", linestyle="--", linewidth=1.2, label="Mean"), |
| 210 | + Line2D([0], [0], color="black", linestyle=":", linewidth=1.2, label="Median"), |
| 211 | + Patch(facecolor="gray", alpha=0.12, label="Mean ± Std"), |
| 212 | + ] |
| 213 | + ax.legend(handles=channel_handles + stat_handles, loc="upper left", fontsize=8) |
156 | 214 |
|
157 | 215 | plt.show() |
0 commit comments