diff --git a/CHANGELOG.md b/CHANGELOG.md index 55c98b2d..f1356e5c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,10 @@ - Make the code compatible with Python 3.12 [#332](https://github.com/litebird/litebird_sim/pull/332) +- plot_fp.py which visualizes focal plane and `DetectorInfo` is implemented. +Also it can generate a dector list file by clicking visualized detectors. +The function is executable by: `python -m litebird_sim.plot_fp` [#345](https://github.com/litebird/litebird_sim/pull/345) + # Version 0.13.0 - **Breaking change**: new API for pointing computation [#319](https://github.com/litebird/litebird_sim/pull/319). Here is a in-depth list of all the breaking changes in this PR: diff --git a/README.md b/README.md index 33ce83c4..c5306812 100644 --- a/README.md +++ b/README.md @@ -112,6 +112,20 @@ command: make.bat html ``` +### Focal plane visualizer + +We can visualize detectors in the focal plane by: +``` +python -m litebird_sim.plot_fp +``` +This software loads the IMo which is installed in the machine you are using. +As the conversation unfolds, an interactive Matplotlib window will appear. +Detectors corresponding to the specified channels are represented as blue dots. +Clicking on a dot reveals the `DetectorInfo` for that detector in real time, highlighted with a red star. +Additionally, if you agree during the conversation to generate a detector file, +a list of starred detectors will be saved into a text file at the designated location after you closed the plot. + +![plot_fp_usage](https://private-user-images.githubusercontent.com/83496454/388083441-93876c66-8b61-40c0-b23f-79a6998e5e33.gif?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MzIxMDQ2MDcsIm5iZiI6MTczMjEwNDMwNywicGF0aCI6Ii84MzQ5NjQ1NC8zODgwODM0NDEtOTM4NzZjNjYtOGI2MS00MGMwLWIyM2YtNzlhNjk5OGU1ZTMzLmdpZj9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDExMjAlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQxMTIwVDEyMDUwN1omWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTZlOTk2M2M3NDgyZjUyY2M4NjFiOGE3ZjliM2RhMWU1YWU2NjFkNmM1ZGQxYzAzM2Y3ZDVmOGI3ZTdhMTAwNjEmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0In0.MoAArSYcXhQgVYaer_O72XHIamO_Ex6B5ITuyVgB8_g) ## Roadmap diff --git a/docs/source/part2.rst b/docs/source/part2.rst index de346164..cb897652 100644 --- a/docs/source/part2.rst +++ b/docs/source/part2.rst @@ -7,6 +7,7 @@ Structure of the framework simulations.rst imo.rst detectors.rst + plot_fp.rst observations.rst data_layout.rst profiling.rst diff --git a/docs/source/plot_fp.rst b/docs/source/plot_fp.rst new file mode 100644 index 00000000..54a49212 --- /dev/null +++ b/docs/source/plot_fp.rst @@ -0,0 +1,44 @@ +.. _plot_fp: + +Focal plane visualization +=========== + +We can visualize detectors in the focal plane by: + +.. code-block:: text + + python -m litebird_sim.plot_fp + +This software loads the IMo, which is installed on the machine you are using. +As the conversation unfolds, an interactive Matplotlib window will appear. + +Detectors corresponding to the specified channels are represented as blue dots. +Clicking on a dot reveals the `DetectorInfo` for that detector in real time, highlighted with a red star. + +Additionally, if you agree during the conversation to generate a detector file, +a list of starred detectors will be saved into a text file at the designated location after you close the plot. + +The format of the detector file is as follows: + ++------------+---------+---------+------------+------------+-----------------------+ +| Telescope | Channel | IMO_NET | Number_det | Scaled_NET | Detector_name | ++------------+---------+---------+------------+------------+-----------------------+ +| LFT | L1-040 | 114.63 | 2/48 | 13.51 | `000_003_003_UB_040_T`| ++------------+---------+---------+------------+------------+-----------------------+ +| LFT | L1-040 | 114.63 | 2/48 | 13.51 | `000_003_003_UB_040_B`| ++------------+---------+---------+------------+------------+-----------------------+ + +The description of each column is as follows: + +- `Telescope`: The telescope name. +- `Channel`: The channel name. +- `IMO_NET`: The NET of the detector in IMo. +- `Number_det`: :math:`N_{\text{selected}}/N_{\text{total}}` where :math:`N_{\text{selected}}` is the number of selected detectors by clicking and :math:`N_{\text{total}}` is the total number of detectors in the channel. +- `Scaled_NET`: The scaled NET of the detectors is given by the following equation: + + .. math:: + + \text{Scaled NET} = \text{NET}_{\text{IMO}} \sqrt{\frac{t_{\text{obs}}}{3} \frac{N_{\text{selected}}}{N_{\text{total}}}} + + where :math:`t_{\text{obs}}` is the observation time in years that you can specify in the conversation. The factor of 3 is the nominal observation time in years. +- `Detector_name`: The detector name. diff --git a/litebird_sim/plot_fp.py b/litebird_sim/plot_fp.py new file mode 100644 index 00000000..cf265560 --- /dev/null +++ b/litebird_sim/plot_fp.py @@ -0,0 +1,353 @@ +# -*- encoding: utf-8 -*- + +import numpy as np +import healpy as hp +import matplotlib.pyplot as plt +import os +from pathlib import Path +import toml +from rich import print +from .quaternions import ( + quat_left_multiply, + quat_right_multiply, + quat_rotation_z, +) +from .simulations import Simulation +from .imo import Imo +from .detectors import DetectorInfo, FreqChannelInfo + +CONFIG_FILE_PATH = Path.home() / ".config" / "litebird_imo" / "imo.toml" + + +class DetectorInfoViewer: + def __init__(self): + self.selected_detector_list = [] + self.scatter = [] + self.fig = None + self.ax1 = None + self.ax2 = None + self.info_box = None + self.x_tot = None + self.y_tot = None + self.x_ch = None + self.y_ch = None + self.ndets_in_channel = None + self.channel_dets_list = [] + self.total_dets_list = [] + self.base_path = None + self.telescope = None + self.channel = None + self.imo_version = None + + def get_det_xy(self, det_info, telescope): + """Get the x, y coordinates of the detector on the focal plane. + + Args: + det_info (DetectorInfo): Detector information. + + Returns: + x (float): The x coordinate of the detector on the focal plane. + y (float): The y coordinate of the detector on the focal plane. + """ + q = det_info.quat.quats[0] + r = np.array([0.0, 0.0, 1.0, 0.0]) + if telescope == "LFT": + # Rotate the FPU 180 degrees for projection if LFT is selected + q_z = quat_rotation_z(np.deg2rad(180)) + for q_val in [q, q_z]: + q_conj = np.array([-q_val[0], -q_val[1], -q_val[2], q_val[3]]) + quat_left_multiply(r, *q_val) + quat_right_multiply(r, *q_conj) + else: + q_conj = np.array([-q[0], -q[1], -q[2], q[3]]) + quat_left_multiply(r, *q) + quat_right_multiply(r, *q_conj) + theta, phi = hp.vec2ang(r[:3]) + x = np.rad2deg(theta) * np.cos(phi) + y = np.rad2deg(theta) * np.sin(phi) + return x, y + + def gen_info_text(self, detector): + """Generate the information text of the detector. + + Args: + detector (DetectorInfo): Detector information. + """ + info_text = rf""" + Detector info. + $\cdot~$name: {detector.name} + $\cdot~$wafer: {detector.wafer} + $\cdot~$pixel: {detector.pixel} + $\cdot~$pixtype: {detector.pixtype} + $\cdot~$channel: {detector.channel} + $\cdot~$sampling_rate_hz: {detector.sampling_rate_hz} + $\cdot~$fwhm_arcmin: {detector.fwhm_arcmin} + $\cdot~$ellipticity: {detector.ellipticity} + $\cdot~$bandcenter_ghz: {detector.bandcenter_ghz} + $\cdot~$net_ukrts: {detector.net_ukrts} + $\cdot~$pol_sensitivity_ukarcmin: {detector.pol_sensitivity_ukarcmin} + $\cdot~$fknee_mhz: {detector.fknee_mhz} + $\cdot~$fmin_hz: {detector.fmin_hz} + $\cdot~$alpha: {detector.alpha} + $\cdot~$pol: {detector.pol} + $\cdot~$orient: {detector.orient} + $\cdot~$quat: {detector.quat.quats[0]} + """ + return info_text + + def gen_detsinfo_text(self, det1, det2): + """Generate the information text of the detector. + + Args: + detector (DetectorInfo): Detector information. + """ + info_text = rf""" +Detector info. + name: {det1.name}, {det2.name} + pol: {det1.pol}, {det2.pol} + quat: {det1.quat.quats[0]} + : {det2.quat.quats[0]} + orient: {det1.orient} + wafer: {det1.wafer} + pixel: {det1.pixel} + pixtype: {det1.pixtype} + channel: {det1.channel} + sampling_rate_hz: {det1.sampling_rate_hz} + fwhm_arcmin: {det1.fwhm_arcmin} + ellipticity: {det1.ellipticity} + bandcenter_ghz: {det1.bandcenter_ghz} + net_ukrts: {det1.net_ukrts} + pol_sensitivity_ukarcmin: {det1.pol_sensitivity_ukarcmin} + fknee_mhz: {det1.fknee_mhz} + fmin_hz: {det1.fmin_hz} + alpha: {det1.alpha} + """ + return info_text + + def generate_dets_list_file( + self, filename, selected_detector_list, duration_yr=1.0 + ): + """Generate a text file with the selected detector list. + The NET on detector is scaled by the`scaling_factor`: + :math:`\sqrt{\frac{duration_yr}{3} \times \frac{N_{\rm dets}^{\rm e2e}}{N_{\rm dets}^{\rm ch}}}` + + Args: + filename (str): The name of the file to be created. + selected_detector_list (list): A list of selected detectors. + duration_yr (float): The duration of the end-to-end simulation in years. + """ + header = "# Telescope Channel IMO_NET Number_det Scaled_NET Detector_name\n" + selected_number_of_dets = len(selected_detector_list) + scaling_factor = np.sqrt( + duration_yr / 3.0 * selected_number_of_dets / self.ndets_in_channel + ) + with open(os.path.join(self.base_path, filename), "w") as file: + file.write(header) + for detector in selected_detector_list: + scaled_net = np.round(detector.net_ukrts * scaling_factor, 2) + line = f"{self.telescope}\t\t{self.channel}\t\t{detector.net_ukrts}\t\t{selected_number_of_dets}/{self.ndets_in_channel}\t\t{scaled_net}\t\t{detector.name}\n" + file.write(line) + print(f"[green]The {filename} is generated.[/green]") + print(f"[green]Location:[/green] [cyan]{self.base_path}/{filename}[/cyan]") + + def on_plot_click(self, event): + """Select the detector by clicking on the plot. + + Args: + event (matplotlib.backend_bases.MouseEvent): The mouse event. + """ + if not event.inaxes: + return + else: + blue = "#1f77b4" + red = "#b41f44" + # distance between the clicked point and the detector + distance = ( + (self.x_ch - event.xdata) ** 2 + (self.y_ch - event.ydata) ** 2 + ) ** 0.5 + if np.min(distance) < 0.3: + sorted_indices = np.argsort(distance) + indices = [sorted_indices[0], sorted_indices[1]] + dets_list = [] + for idx in indices: + detector = self.channel_dets_list[idx] + dets_list.append(detector) + if self.scatter[idx].get_markerfacecolor() == blue: + self.scatter[idx].set_markerfacecolor(red) + self.scatter[idx].set_markeredgecolor(red) + self.scatter[idx].set_marker("*") + self.scatter[idx].set_markersize(12) + self.selected_detector_list.append(detector) + elif self.scatter[idx].get_markerfacecolor() == red: + self.scatter[idx].set_markerfacecolor(blue) + self.scatter[idx].set_markeredgecolor(blue) + self.scatter[idx].set_marker("o") + self.scatter[idx].set_markersize(8) + if detector in self.selected_detector_list: + self.selected_detector_list.remove(detector) + info_text = self.gen_detsinfo_text(dets_list[0], dets_list[1]) + self.info_box.set_text(info_text) + self.fig.canvas.draw() + + def ask_yes_or_no(self): + while True: + ans = input().lower() + if ans == "y": + print("[green]Create a detector list file.[/green]") + break + elif ans == "n": + print("[green]No detector list file will be created.[/green]") + break + else: + print("[ref]Invalid input. Please enter 'y' or 'n'.[/red]") + print("[green]Do you want to make a detector list file? [y/n][/green]") + return ans + + def extract_location_from_toml(self, file_path): + """Extract the location of IMo from the toml file. + + Args: + file_path (str): The path to the toml file. + """ + with open(file_path, "r") as file: + data = toml.load(file) + loc = data["repositories"][0]["location"] + return loc + + def main(self): + if not CONFIG_FILE_PATH.exists(): + imo = Imo() + self.imo_version = "vPTEP" + else: + IMO_ROOT_PATH = self.extract_location_from_toml(CONFIG_FILE_PATH) + imo = Imo(flatfile_location=os.path.join(IMO_ROOT_PATH, "schema.json")) + versions = list(imo.imoobject.releases.keys()) + versions_with_idx = [f"({i+1}). {ver}" for i, ver in enumerate(versions)] + print( + f"[green]Available IMO versions:[/green] [cyan]{versions_with_idx}[/cyan]" + ) + print("[green]Input IMO version's number: [/green]") + version_idx = input() + self.imo_version = versions[int(version_idx) - 1] + + sim = Simulation(random_seed=None, imo=imo) + + print( + "[green]Input telescope's number:[/green] [cyan]['(1). LFT', '(2). MFT', '(3). HFT'][/cyan]" + ) + telescope_id = input() + + if telescope_id == "1": + self.telescope = "LFT" + elif telescope_id == "2": + self.telescope = "MFT" + elif telescope_id == "3": + self.telescope = "HFT" + else: + raise ValueError("Invalid telescope number. Please input 1, 2 or 3.") + + inst_info = sim.imo.query( + f"/releases/{self.imo_version}/satellite/{self.telescope}/instrument_info" + ) + channel_list = inst_info.metadata["channel_names"] + # add index to the channel list + channel_list_with_idx = [ + f"({i+1}). {channel}" for i, channel in enumerate(channel_list) + ] + + print( + f"[green]The availavle channels are:[/green] [cyan]{channel_list_with_idx}[/cyan]" + ) + print("[green]Input the channel's number:[/green]") + channel_idx = input() + self.channel = channel_list[int(channel_idx) - 1] + + print("[green]Do you want to make a detector list file? (y/n) [/green]") + ans = self.ask_yes_or_no() + if ans == "y": + print( + "[green]Input mission duration to define a scaling factor for NET (unit: yr):[/green]" + ) + duration_yr = float(input()) + print("[green]Specify the directory to save:[/green]") + self.base_path = input() + if self.base_path == "": + self.base_path = "./" + print("[green]The file will be saved in the current directory.[/green]") + if self.base_path.endswith("/"): + self.base_path = self.base_path[:-1] + + for ch in channel_list: + channel_info = FreqChannelInfo.from_imo( + imo=imo, + url=f"/releases/{self.imo_version}/satellite/{self.telescope}/{ch}/channel_info", + ) + for detector_name in channel_info.detector_names: + det = DetectorInfo.from_imo( + imo=imo, + url=f"/releases/{self.imo_version}/satellite/{self.telescope}/{ch}/{detector_name}/detector_info", + ) + if self.channel == ch: + self.channel_dets_list.append(det) + else: + self.total_dets_list.append(det) + + self.x_tot = np.zeros(len(self.total_dets_list)) + self.y_tot = np.zeros(len(self.total_dets_list)) + self.x_ch = np.zeros(len(self.channel_dets_list)) + self.y_ch = np.zeros(len(self.channel_dets_list)) + self.ndets_in_channel = len(self.channel_dets_list) + + # Make a figure + self.fig = plt.figure(figsize=(10, 8)) + self.ax1 = self.fig.add_subplot(1, 2, 1, aspect="equal") + self.ax1.set_title(self.channel) + self.ax1.plot() + self.scatter = [] + + for i, det in enumerate(self.total_dets_list): + self.x_tot[i], self.y_tot[i] = self.get_det_xy(det, self.telescope) + self.ax1.scatter(self.x_tot, self.y_tot, marker="x", s=25, color="black") + + for i, det in enumerate(self.channel_dets_list): + self.x_ch[i], self.y_ch[i] = self.get_det_xy(det, self.telescope) + self.scatter.append( + self.ax1.plot( + self.x_ch[i], self.y_ch[i], "o", markersize=8, color="#1f77b4" + )[0] + ) + + self.ax1.set_xlabel(r"$\theta\cos(\phi)$ [degrees]") + self.ax1.set_ylabel(r"$\theta\sin(\phi)$ [degrees]") + + self.ax2 = self.fig.add_subplot(1, 2, 2, aspect="equal") + self.info_box = self.ax2.text( + 0.02, + 0.98, + "", + transform=self.ax2.transAxes, + va="top", + ha="left", + bbox=dict(boxstyle="round", facecolor="wheat", alpha=0.5), + ) + self.ax2.set_axis_off() + + print("[blue]Click is available...[/blue]") + self.fig.canvas.mpl_connect("button_press_event", self.on_plot_click) + plt.tight_layout() + plt.show() + + # Save the detector list file. + if ans == "y": + filename = "detectors_" + self.telescope + "_" + self.channel + "_T+B.txt" + self.selected_detector_list = sorted( + self.selected_detector_list, key=lambda detector: detector.name + ) + self.generate_dets_list_file( + filename, self.selected_detector_list, duration_yr + ) + + +if __name__ == "__main__": + viewer = DetectorInfoViewer() + viewer.main()