import pathlib  # for filepath path tooling
import lzma  # to decompress the iCOM file

import numpy as np  # for array tooling
import matplotlib.pyplot as plt  # for plotting

"""
Explanation: Directly compare Monaco to Elekta Linac iCOM
This notebook uses the PyMedPhys library to compare the collected iCOM delivery data directly to the recorded plan within Monaco's tel files.

Description of method
There are a few key stages required to make this notebook work. Firstly the iCOM stream from the Linac needs to be recorded to disk, this iCOM stream then needs to be archived and grouped by patient.
Once these iCOM streams are running the plan needs to be delivered on the machine being listened to. From this point the iCOM stream is then able to be compared directly to the Monaco tel file.
This comparison is done by calculating an MU Density upon which the difference is reported on using the Gamma comparison tool.
iCOM stream service
A Windows service is created to listen to the iCOM stream. This service is made from a .bat file, an example of which can be seen at:

https://github.com/CCA-Physics/physics-server/blob/8a0954e5/RCCC/icom/harry_listening.bat

To create this service the nssm tool is used, with an example of its usage available at:

https://github.com/CCA-Physics/physics-server/blob/8a0954e5/RCCC/icom/services-install.bat

Warning: Take Note
A force closing of the service which is listening to the iCOM stream done in such a way that it is not able to properly close down the listening socket will cause the Linac being listened to raise an interlock which won't interupt the beam, but it will not let a new beam be delivered until the machine is logged out of and re-logged back in.
This can happen when the service is force killed, the host machine is force shutdown, or the network connection is abruptly disconnected. Normal shutdown of the service or machine should not have this effect.
Grouping iCOM stream by patient
For the iCOM stream to be easily indexed in the future deliveries are stored by patient id and name. This is done by setting up a .bat file to run on machine boot. An example .bat file that achieves this can be seen at:

https://github.com/CCA-Physics/physics-server/blob/8a0954e5/RCCC/icom/patient_archiving.bat

Reading the iCOM file, Monaco file, and comparing them
The resulting files are then loaded into a PyMedPhys Delivery object from which an MU Density can be calculated and used as a comparison and reporting tool.
These steps will be further expanded on below, prior to the lines of code that implement them.
Importing the required libraries
Third party libraries
End of explanation
"""


# Makes it so that any changes in pymedphys is automatically
# propagated into the notebook without needing a kernel reset.
from IPython.lib.deepreload import reload
%load_ext autoreload
%autoreload 2

import pymedphys

"""
Explanation: Importing PyMedPhys
PyMedPhys contains all the tooling to log the iCOM stream, read the iCOM stream, read Monaco tel files, create an MU Density, and compare the results using a Gamma implementation
End of explanation
"""


patient_id = '015112'

"""
Explanation: Patient ID Configuration
Provide here the Patient ID, this will be used to look up the iCOM file record as well as the Monaco tel file record.
End of explanation
"""


icom_directory = pathlib.Path(r'\\physics-server\iComLogFiles\patients')
monaco_directory = pathlib.Path(r'\\monacoda\FocalData\RCCC\1~Clinical')

"""
Explanation: File Path Configurations
Here is where the various root paths of the data has been stored.
Input directories
End of explanation
"""


output_directory = pathlib.Path(r'S:\Physics\Patient Specific Logfile Fluence')
pdf_directory = pathlib.Path(r'P:\Scanned Documents\RT\PhysChecks\Logfile PDFs')

"""
Explanation: Output directories
The output_directory below presents the location where raw results are saved for the permanent record.
The pdf_directory is designed to be a sub-directory of the Mosaiq eScan directory so that created PDFs can be easily imported into Mosaiq for reporting purposes.
End of explanation
"""


GRID = pymedphys.mudensity.grid()
COORDS = (GRID["jaw"], GRID["mlc"])

GAMMA_OPTIONS = {
    'dose_percent_threshold': 2,  # Not actually comparing dose though
    'distance_mm_threshold': 0.5,
    'local_gamma': True,
    'quiet': True,
    'max_gamma': 2,
}

"""
Explanation: MU Density and Gamma configuration
MU Density and Gamma have a range of options available to them. Here are is where those options are configured.
End of explanation
"""


all_tel_paths = list(monaco_directory.glob(f'*~{patient_id}/plan/*/tel.1'))
all_tel_paths

plan_names_to_choose_from = [
    path.parent.name for path in all_tel_paths
]
plan_names_to_choose_from

"""
Explanation: Choosing from the available Monaco Plans and iCOM Deliveries
For the Patient ID provided the following two cells list out all the tel files and iCOM delivery files that were found. Run both of these cells below and then choose the appropriate plan and delivery in the third cell.
Monaco Plans
End of explanation
"""


icom_deliveries = list(icom_directory.glob(f'{patient_id}_*/*.xz'))
icom_deliveries

icom_files_to_choose_from = [
    path.stem for path in icom_deliveries
]
icom_files_to_choose_from

"""
Explanation: iCOM Deliveries
End of explanation
"""


monaco_plan_name = 'LeftIlium1'  # plan directory name
icom_delivery = '20200213_133208'  # iCOM timestamp

"""
Explanation: Plan and iCOM choice
Provide the directory name for the monaco plan and the timestamp for the iCOM delivery which you wish to compare.
End of explanation
"""


tel_path = list(monaco_directory.glob(f'*~{patient_id}/plan/{monaco_plan_name}/tel.1'))[-1]
tel_path

icom_path = list(icom_directory.glob(f'{patient_id}_*/{icom_delivery}.xz'))[-1]
icom_path

"""
Explanation: Resulting paths found matching provided search query
End of explanation
"""


with lzma.open(icom_path, 'r') as f:
    icom_stream = f.read()

"""
Explanation: Loading the iCOM file
In order to save space on disk the patient iCOM streams are compressed using lzma. The following code opens these compressed files and gets the raw iCOM stream.
End of explanation
"""


# Print out available methods and attributes on the Delivery object
[command for command in dir(pymedphys.Delivery) if not command.startswith('_')]

"""
Explanation: Delivery Objects
Within PyMedPhys there is a Delivery object. This object can be created from a range of sources such as RT Plan DICOM files, Mosaiq SQL queries, iCOM data, trf log files, as well as Monaco tel files.
From this Delivery object many tasks can be undergone. The available methods and attributes on the Delivery object are given below:
End of explanation
"""


delivery_icom = pymedphys.Delivery.from_icom(icom_stream)
delivery_tel = pymedphys.Delivery.from_monaco(tel_path)

"""
Explanation: Creating the Delivery Objects
We can create two of these Delivery objects, one from the iCOM stream, and the other from the Monaco tel file.
End of explanation
"""


mudensity_icom = delivery_icom.mudensity()
mudensity_tel = delivery_tel.mudensity()

"""
Explanation: Using the Delivery Objects
Once we have two Delivery objects we can calculate the MU Density of these. Of note, this same method of using the Delivery object can also be employed to compare to an RT Plan file, Elekta Linac TRF log file, or Mosaiq SQL.
End of explanation
"""


def to_tuple(array):
    return tuple(map(tuple, array))

gamma = pymedphys.gamma(
    COORDS,
    to_tuple(mudensity_tel),
    COORDS,
    to_tuple(mudensity_icom),
    **GAMMA_OPTIONS
)

"""
Explanation: Calculating Gamma
PyMedPhys also has within it tooling to calculate Gamma. This is done below.
End of explanation
"""


def plot_gamma_hist(gamma, percent, dist):
    valid_gamma = gamma[~np.isnan(gamma)]

    plt.hist(valid_gamma, 50, density=True)
    pass_ratio = np.sum(valid_gamma <= 1) / len(valid_gamma)

    plt.title(
        "Local Gamma ({0}%/{1}mm) | Percent Pass: {2:.2f} % | Max Gamma: {3:.2f}".format(
            percent, dist, pass_ratio * 100, np.max(valid_gamma)
        )
    )

def plot_and_save_results(
    mudensity_tel,
    mudensity_icom,
    gamma,
    png_filepath,
    pdf_filepath,
    header_text="",
    footer_text="",
):
    diff = mudensity_icom - mudensity_tel
    largest_item = np.max(np.abs(diff))

    widths = [1, 1]
    heights = [0.3, 1, 1, 1, 0.1]
    gs_kw = dict(width_ratios=widths, height_ratios=heights)

    fig, axs = plt.subplots(5, 2, figsize=(10, 16), gridspec_kw=gs_kw)
    gs = axs[0, 0].get_gridspec()

    for ax in axs[0, 0:]:
        ax.remove()

    for ax in axs[1, 0:]:
        ax.remove()

    for ax in axs[4, 0:]:
        ax.remove()

    axheader = fig.add_subplot(gs[0, :])
    axhist = fig.add_subplot(gs[1, :])
    axfooter = fig.add_subplot(gs[4, :])

    axheader.axis("off")
    axfooter.axis("off")

    axheader.text(0, 0, header_text, ha="left", wrap=True, fontsize=30)
    axfooter.text(0, 1, footer_text, ha="left", va="top", wrap=True, fontsize=6)

    plt.sca(axs[2, 0])
    pymedphys.mudensity.display(GRID, mudensity_tel)
    axs[2, 0].set_title("Monaco Plan MU Density")

    plt.sca(axs[2, 1])
    pymedphys.mudensity.display(GRID, mudensity_icom)
    axs[2, 1].set_title("Recorded iCOM MU Density")

    plt.sca(axs[3, 0])
    pymedphys.mudensity.display(
        GRID, diff, cmap="seismic", vmin=-largest_item, vmax=largest_item
    )
    plt.title("iCOM - Monaco")

    plt.sca(axs[3, 1])
    pymedphys.mudensity.display(GRID, gamma, cmap="coolwarm", vmin=0, vmax=2)
    plt.title(
        "Local Gamma | "
        f"{GAMMA_OPTIONS['dose_percent_threshold']}%/"
        f"{GAMMA_OPTIONS['distance_mm_threshold']}mm")

    plt.sca(axhist)
    plot_gamma_hist(
        gamma, 
        GAMMA_OPTIONS['dose_percent_threshold'], 
        GAMMA_OPTIONS['distance_mm_threshold'])

    return fig

"""
Explanation: Create Plotting and Reporting Functions
So that we can view the result as well as create a PDF that can be stored within Mosaiq the following functions create these plots using the matplotlib library.
End of explanation
"""


results_dir = output_directory.joinpath(patient_id, tel_path.parent.name, icom_path.stem)
results_dir.mkdir(exist_ok=True, parents=True)

header_text = (
    f"Patient ID: {patient_id}\n"
    f"Plan Name: {tel_path.parent.name}\n"
)

footer_text = (
    f"tel.1 file path: {str(tel_path)}\n"
    f"icom file path: {str(icom_path)}\n"
    f"results path: {str(results_dir)}"
)

png_filepath = str(results_dir.joinpath("result.png").resolve())
pdf_filepath = str(pdf_directory.joinpath(f"{patient_id}.pdf").resolve())

fig = plot_and_save_results(
    mudensity_tel, mudensity_icom, 
    gamma, png_filepath, pdf_filepath, 
    header_text=header_text, footer_text=footer_text
)

fig.tight_layout()
plt.savefig(png_filepath, dpi=300)
plt.show()

"""
Explanation: Plotting and saving the report
Now that have our data, and have formatted our report as we wish let's create and save this report as a png.
End of explanation
"""


!magick convert "{png_filepath}" "{pdf_filepath}"

"""
Explanation: Converting PNG to PDF for importing into Mosaiq
To create a pdf, the just created png file can be converted to pdf. To do this the tool imagemagick needs to be installed on your system. If you install this now you will need to reset your Jupyter server in a new command prompt so that the magick command is available within your path.
End of explanation
"""


