#!/usr/bin/env python

"""
Script report_protocol
==========================
This script is used for generating CSV files and their corresponding plots for
reporting a NPFC run on a given dataset.

The folder tree:
"""

# standard
import warnings
import logging
import argparse
from datetime import datetime
from collections import Counter
import re
from pathlib import Path
from collections import OrderedDict
# data handling
import json
import pandas as pd
from pandas import DataFrame
import networkx as nx
# data visualization
from matplotlib import pyplot as plt
import seaborn as sns
from matplotlib.ticker import FuncFormatter
from pylab import savefig
# chemoinformatics
import rdkit
from rdkit import Chem
# docs
from typing import List
from typing import Tuple
from rdkit import RDLogger
# custom libraries
import npfc
from npfc import utils
from npfc import load
from npfc import save
from npfc import fragment_combination


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ PRIVATE FUNCTIONS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #


def _print_title(title, logger, pad):
    pad = pad * 2
    logger.info("=" * pad)
    logger.info(title.upper().center(pad ))
    logger.info("=" * pad)


def _get_chunks(WD: str, pattern: str) -> List[str]:
    """This function returns the list of all chunks found in a given directory.
    It does not use the _000_ notation because it can also be applied to fragments.

    :param WD: the WD where all subfolders for each step are present (i.e. 'natural/dnp/data')
    :param pattern: the pattern that desired chunks contain (i.e. '1_input/data/*([0-9][0-9][0-9])?.sdf.gz')
    :return: the list of chunks
    """
    WD = Path(WD)
    pattern = re.compile(pattern)

    chunks = [str(x) for x in list(WD.glob("*"))]

    return sorted(list(filter(pattern.match, chunks)))


def _parse_std_chunks(chunks: List[str]) -> DataFrame:
    """Parse all output files of a category (passed, filtered or error) for the std step and return a corresponding a results summary.

    :param chunks: output files for a category of std results
    :return: summary DF with counts
    """
    # parse all files
    dfs = []
    for c in chunks:
        df = pd.read_csv(c, sep="|", compression="gzip").groupby("task").count()[['status']].rename({'status': 'Count'}, axis=1)
        if len(df.index) > 0:
            dfs.append(df)

    # if no case was found, return an empty dataframe
    if len(dfs) == 0:
        df = pd.DataFrame([], columns=["Count", "Category"])
        return df

    # concatenate all dataframes and compute the sum of all counts
    df = pd.concat(dfs)
    df["Category"] = df.index  # I don't know how to group by index!
    df = df.groupby("Category").sum()
    df["Category"] = df.index.map(lambda x: x.replace('filter_', ''))
    # df['Perc_status'] = df['Count'].map(lambda x: f"{x/tot_mols:.2%}")

    return df


def _step_is_valid(subdir: str, errors_data: OrderedDict, errors_log: OrderedDict) -> bool:
    """Check is the data and log output files in a subdirectory (=step) are valid (no missing chunks).

    :param subdir: either a string or an iterable of strings
    :param errors_data: a dictionary with key=subdir and values=missing chunks ids
    :param errors_log: same as errors_data but with log output instead of data
    :return: Return True if no missing chunk in step were identified, False otherwise
    """
    # make subdir into an iterable if not already one
    if isinstance(subdir, str):
        subdir = [subdir]
    # check for missing data/logs
    for step in subdir:
        if step in errors_data.keys():
            return False
        elif step in errors_data.keys():
            return False
    return True


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ FUNCTIONS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #


def save_barplot(df: DataFrame,
                 output_png: str,
                 x_name: str,
                 y_name: str,
                 title: str,
                 color: str,
                 x_label: str = None,
                 y_label: str = None,
                 rotate_x: int = 0,
                 perc_labels: str = None,
                 perc_label_size: int = 15,
                 fig_size: Tuple[int] = (24, 12),
                 force_order: bool = True,
                 ):
    """This function helps for computing automated barplots using seaborn.
    It sets up somewhat standardized figure output for a harmonized rendering.

    :param df: the DataFrame with data to plot
    :param output_png: the output png full file name
    :param x_name: DF column name to use for x-axis
    :param y_name: DF column name to use for y-axis
    :param x_label: the name to display on the plot for x-axis
    :param y_label: the name to display on the plot for y-axis
    :param rotate_x: rotate the x-axis ticks anticlock-wise
    :param perc_labels: DF column name to use display percentage labels above bars
    :param perc_label_size: annotation text size for perc_labels
    :param color: color to use for bars, theoritically could also be a list of colors
    :param fig_size: tuple of integers defining the plot dimensions (x, y)
    :param force_order: force the plot to display the bars in the provided order. I do not know why sometimes this is needed, sometimes not depending on the plot.
    :return: the figure in searborn format
    """
    # delete existing file for preventing stacking of plots
    p = Path(output_png)
    if p.exists():
        p.unlink()

    # general style for plotting
    # sns.set(rc={'figure.figsize': fig_size})
    plt.figure(figsize=fig_size)
    sns.set_style('whitegrid', {'axes.edgecolor': '0.2'})
    sns.set_context("paper", font_scale=2)

    # barplot
    if force_order:
        ax = sns.barplot(x=df[x_name], y=df[y_name], color=color, order=df[x_name])
    else:
        ax = sns.barplot(x=df[x_name], y=df[y_name], color=color)
    ax.set_title(title, fontsize=24, y=1.02)
    ax.tick_params(labelsize=20)
    ax.tick_params(axis='x', rotation=rotate_x)
    # ax.set_xticklabels(df[x_name])
    ax.set_xlabel(x_label,fontsize=25, labelpad=20)
    ax.set_ylabel(y_label,fontsize=25, labelpad=20)
    # if no entry, then y ticks get all confused
    if df[y_name].sum() == 0:
        ax.set_yticks((0, 1, 2, 3, 4, 5))
        ax.set_ylim((0, 5))

    # format y labels
    ylabels = ['{:,.0f}'.format(x) for x in ax.get_yticks()]
    ax.set_yticklabels(ylabels)
    # add percentage annotations
    if perc_labels is not None:
        for i in range(len(df.index)):
            row = df.iloc[i]
            ax.text(row.name, row[y_name], row[perc_labels], color='black', ha='center', fontdict={'fontsize': perc_label_size})
    # save
    figure = ax.get_figure()
    if rotate_x > 0:
        figure.subplots_adjust(bottom=0.2)
    figure.savefig(output_png, dpi=600)
    plt.clf()

    return figure


def save_kdeplot(df: DataFrame,
                 output_png: str,
                 x_name: str,
                 title: str,
                 color: str,
                 x_label: str = None,
                 y_label: str = None,
                 normalize_x: bool = True,
                 fig_size: Tuple[int] = (24, 12),
                 ):
    """This function helps for computing automated kdeplots using seaborn.
    It sets up somewhat standardized figure output for a harmonized rendering.

    :param df: the DataFrame with data to plot
    :param output_png: the output png full file name
    :param x_name: DF column name to use for x-axis
    :param x_label: the name to display on the plot for x-axis
    :param y_label: the name to display on the plot for y-axis
    :param color: color to use for bars, theoritically could also be a list of colors
    :param fig_size: tuple of integers defining the plot dimensions (x, y)
    :return: the figure in searborn format
    """
    # general style for plotting
    sns.set(rc={'figure.figsize': fig_size})
    sns.set_style('whitegrid', {'axes.edgecolor': '0.2'})
    sns.set_context("paper", font_scale=2)

    ax = sns.kdeplot(df[x_name], shade=True, label='', color=color)
    ax.set_title(title, fontsize=24, y=1.02)
    ax.tick_params(labelsize=20)
    ax.tick_params(axis='x', rotation=0)
    ax.set_xlim(0, 1)
    #ax.set_xticklabels(df[x_name])
    ax.set_xticklabels(['{:,.0%}'.format(x) for x in ax.get_xticks()])
    ax.set_xlabel(x_label,fontsize=25, labelpad=20)
    ax.set_ylabel(y_label,fontsize=25, labelpad=20)

    # save
    figure = ax.get_figure()
    figure.savefig(output_png, dpi=600)
    plt.clf()

    return figure


def _get_missing_chunkids(WD: Path, subfolders: list, chunks_ini_id: list, check_on: str) -> OrderedDict:
    """Iterate over defined subfolders in WD and compare the number of chunks found in sub-subfolder check_on (data or log) to check
    if files are missing. An OrderedDict is computed with for each key the subfolder and as values the list of missing chunk ids.

    :param WD: the WD
    :param subfolders: the list of subfolders to check (02_load, 03_deglyco, etc.). For any subfolder with '_std' in its name, then number of counted chunks if checking on data will be a third (because there are 3 outputs per chunk: passed, filtered and error).
    :param: chunks_ini_id: the list of chunk ids defined in the step chunk
    """
    # init
    errors = OrderedDict()
    steps = [s.split('_')[1] for s in subfolders]
    # no chunks
    if len(chunks_ini_id) == 0:  # so no reference, just look if there are the expected number of outputs per subfolder
        logger.info("NO CHUNK SUBFOLDER COULD BE FOUND, ASSUMING THIS IS A FRAGMENTS DATASET")
        for step, sf in zip(steps, subfolders):
            error = False
            sf_output_files = [s.name for s in WD.glob(f"{sf}/{check_on}/*")]
            if len(sf_output_files) < 1:
                error = True
            if check_on == 'data':
                if ('_std' in sf and len(sf_output_files) < 3) or (not '_std' in sf and len(sf_output_files) < 1):
                    error = True
                else:
                    error = False
            if error:
                errors[step] = ["Missing output files!"]
    else:
        # chunks
        for step, sf in zip(steps, subfolders):
            sf_chunks = [s.name for s in WD.glob(f"{sf}/{check_on}/*")]
            if sf == '02_load':
                sf_chunks_id = sorted([s.split("_")[-1] for s in sf_chunks])
            else:
                sf_chunks_id = sorted([s.split("_")[-2] for s in sf_chunks])
            # check number of chunks to see if there is a problem, 3 output chunks for std so speical case
            if check_on == 'data':
                if ('_std' in sf and (len(sf_chunks_id)/3) < len(chunks_ini_id)) or (not '_std' in sf and len(sf_chunks_id) < len(chunks_ini_id)):
                    error = True
                else:
                    error = False
            # no this problem with log files
            else:
                if len(sf_chunks_id) < len(chunks_ini_id):
                    error = True
                else:
                    error = False

            if error:
                errors[step] = [x for x in chunks_ini_id if x not in sf_chunks_id]

    return errors


def _get_subfolders(WD: Path) -> dict:
    d_subfolders = {}
    d_subfolders['excluded'] = ["00_raw", "01_chunk"]
    d_subfolders['included'] = sorted([d.name for d in list(WD.glob('[0-9][0-9]_*')) if d.is_dir() if d.name not in d_subfolders['excluded']])
    return d_subfolders


def find_missing_chunks(WD: Path) -> tuple:
    """This functions will look up the expected amount of chunks produced by the
    "1_input" step and verify if all chunks are present in subsequent tasks.
    If not, it will raise a ValueError.

    Subfolders 0_raw and 1_input are excluded from the checks.

    :param WD: Working Directory of the pipeline to check (parent folder of the smk file). Can also be a str object.
    :return: a tuple of OrderedDict with steps as keys and the list of missing chunk ids as values
    """
    # setup WD
    if not isinstance(WD, Path):
        WD = Path(WD)

    # initialize the references for comparing
    logger.info(f"COUNTING NUMBER OF EXPECTED CHUNKS BY USING CHUNK STEP AS REFERENCE")
    subdir = f"{str(WD)}/01_chunk/data"
    logger.info(f"CHUNK SUBFOLDER: {subdir}")
    pattern = f"*"
    chunks_ini = [s.name for s in list(Path(subdir).glob(pattern))]
    chunks_ini_id = [str(s).split("_")[-1].split(".")[0] for s in chunks_ini]
    num_chunks_ini = len(chunks_ini)
    if num_chunks_ini == 0:
        logger.error(f"NO CHUNK FILES COULD BE FOUND AT '{subdir}'!")

    logger.info(f"NUMBER OF EXPECTED CHUNKS FOR EACH SUBSTEP: {num_chunks_ini:,}")
    # define what subfolders are to be checked
    d_subfolders = _get_subfolders(WD)

    logger.info(f"EXCLUDING SUBFOLDER(S) {', '.join(d_subfolders['excluded'])}")
    logger.info(f"INCLUDING SUBFOLDER(S) {', '.join(d_subfolders['included'])}")

    # record errors in case of missing chunks
    logger.info(f"CHECKING FOR DATA FILES")
    errors_data = _get_missing_chunkids(WD, d_subfolders['included'], chunks_ini_id, 'data')
    logger.info(f"CHECKING FOR LOG FILES")
    errors_log = _get_missing_chunkids(WD, d_subfolders['included'], chunks_ini_id, 'log')

    return errors_data, errors_log

def get_df_load(WD: str) -> DataFrame:
    """Get a DF summarizing the load step.

    :param WD: the directory of the load step
    :return: a DF summarizing the load step
    """
    logger.info("PREP -- COMPUTING LOAD RESULTS")
    # iterate over the log files to count status
    pattern = ".*([0-9]{3})?.log"
    chunks = _get_chunks(f"{WD}/log", pattern)
    logger.info(f"FOUND {len(chunks):,d} CHUNKS")
    # initiate counts
    num_passed = 0
    num_errors = 0
    for c in chunks:
        df = pd.read_csv(c, sep="@", header=None)  # char not found in the log file so we can extract all lines as one column
        records = df[df[0].str.contains("FAILURE")].iloc[0][0].split()
        total = int(records[6])
        errors = int(records[9])
        passed = int(df[df[0].str.contains("SAVED")].iloc[0][0].split()[6])
        logger.debug(f"{c} => PASSED: {passed}/{total} ({errors} ERRORS)")
        num_passed += passed
        num_errors += errors
    # create a dataframe with counts
    df = pd.DataFrame({'Category': ['loaded', 'cannot_load'], 'Count': [num_passed, num_errors]})
    logger.info(f"PREP -- RESULTS FOR LOADING MOLECULES:\n\n{df}\n")

    return df


def get_df_deglyco(WD: str) -> DataFrame:
    """Get a DF summarizing the deglycoslyation step.

    :param WD: the directory of the deglyco step
    :return: a DF summarizing the deglycoslyation step
    """
    if WD.endswith('/'):
        WD = WD[0:-1]
    logger.info("PREP -- COMPUTING DEGLYCO RESULTS")
    # iterate over the log files to count status
    pattern = ".*([0-9]{3})?_deglyco\.csv"
    chunks = _get_chunks(f"{WD}/log", pattern)
    logger.info(f"PREP -- FOUND {len(chunks):,d} CHUNKS")
    # initiate counts
    num_tot = 0
    num_unchanged = 0
    num_deglycosylated = 0
    num_failed = 0
    num_error = 0
    # count status
    for c in chunks:
        df = pd.read_csv(c, sep="|")
        num_tot += len(df.index)
        num_unchanged += len(df[df['status'] == 'unchanged'].index)
        num_deglycosylated += len(df[df['status'] == 'deglycosylated'].index)
        num_failed += len(df[df['status'] == 'failed'].index)
        num_error += len(df[df['status'] == 'error'].index)
    # create a dataframe with counts
    df = pd.DataFrame({'Category': ['deglycosylated', 'unchanged', 'failed', 'error'], 'Count': [num_deglycosylated, num_unchanged, num_failed, num_error]})
    logger.info(f"PREP -- RESULTS FOR DEGLYCO:\n\n{df}\n")

    return df


def get_df_std_passed(WD: str) -> DataFrame:
    """Get a DF summarizing the passed results of the sandardization step.

    :param WD: the directory of the std step
    :return: a DF summarizing passed results of the sandardization step
    """
    logger.info("PREP -- COMPUTING STD PASSED RESULTS")
    # iterate over the log files to count status
    pattern = ".*([0-9]{3})?_passed\.csv.gz"
    chunks = _get_chunks(f"{WD}/data", pattern)
    logger.info(f"PREP -- FOUND {len(chunks):,d} CHUNKS")
    # initiate counts
    num_passed = 0
    for c in chunks:
        df = pd.read_csv(c, sep="|")
        num_passed += len(df.index)
    df = pd.DataFrame({"Category": ['passed'], 'Count': [num_passed]})
    logger.info(f"PREP -- RESULTS FOR STD PASSED:\n\n{df}\n")

    return df


def get_df_std_filtered(WD: str) -> DataFrame:
    """Get a DF summarizing the filtered results of the sandardization step.

    :param WD: the directory of the std step
    :return: a DF summarizing filtered results of the sandardization step
    """
    logger.info("PREP -- COMPUTING STD FILTERED RESULTS")
    # iterate over the log files to count status
    pattern = ".*([0-9]{3})?_filtered\.csv.gz"
    chunks = _get_chunks(f"{WD}/data", pattern)
    logger.info(f"PREP -- FOUND {len(chunks):,d} CHUNKS")
    # initiate counts
    df = _parse_std_chunks(chunks)
    logger.info(f"RESULTS FOR STD FILTERED:\n\n{df}\n")

    return df


def get_df_std_error(WD: str) -> DataFrame:
    """Get a DF summarizing the error results of the sandardization step.

    :param WD: the directory of the std step
    :return: a DF summarizing error results of the sandardization step
    """
    logger.info("PREP -- COMPUTING STD ERROR RESULTS")
    # iterate over the log files to count status
    pattern = ".*([0-9]{3})?_error\.csv.gz"
    chunks = _get_chunks(f"{WD}/data", pattern)
    logger.info(f"PREP -- FOUND {len(chunks):,d} CHUNKS")
    # initiate counts
    df = _parse_std_chunks(chunks)
    logger.info(f"PREP -- RESULTS FOR STD ERROR:\n\n{df}\n")

    return df


def get_df_dedupl(WD: str) -> DataFrame:
    """Get a DF summarizing the results of the deduplication step.

    :param WD: the directory of the std step
    :return: a DF summarizing results of the deduplication step
    """
    logger.info("PREP -- COMPUTING DEDUPL RESULTS")
    # iterate over the log files to count status
    pattern = ".*([0-9]{3})?_dedupl\.log"
    chunks = _get_chunks(f"{WD}/log", pattern)
    logger.info(f"PREP -- FOUND {len(chunks):,d} CHUNKS")
    # initiate counts
    num_tot = 0
    num_passed = 0
    num_filtered = 0
    for c in chunks:
        df = pd.read_csv(c, sep="@", header=None)  # char not found in the log file so we can extract all lines as one column
        passed, total = [int(x) for x in df[df[0].str.contains("REMAINING MOLECULES")].iloc[0][0].split("MOLECULES:")[1].split("/")]
        num_passed += passed
        num_filtered += total - passed
        num_tot += total
    # create a dataframe with counts
    df = pd.DataFrame({'Category': ['unique', 'duplicate'], 'Count': [num_passed, num_filtered]})
    logger.info(f"PREP -- RESULTS FOR DEDUPL:\n\n{df}\n")

    return df


def get_dfs_prep(WD: str) -> Tuple[DataFrame]:
    """Get a list of DFs summarizing the whole preprocess superstep: load, deglyco, std and dedupl.

    - DF_deglyco is the detailed summary of deglycosylation appended with the number of mols that did not get processed because they could not be loaded (NA).
    - DF_prep_filtered is the detailed summary of std and dedupl
    - DF_prep_error is the detailed summary of std and load
    - DF_prep_all is the general summary with the final number of passed, filtered and error molecules.

    :param WD: the main directory of the dataset data (i.e. 'natural/dnp/data')
    :return: a list of DFs of interest: [DF_deglyco, DF_prep_filtered, DF_prep_error, DF_prep_all]
    """

    logger.info("PREP -- COMPUTE RESULTS FOR PREPROCESS")
    logger.info("PREP -- PROPRESS CONTAINS LOAD, DEGLYCO, STD AND DEDUPL STEPS")

    # define subfolders
    p = Path(WD)
    WD_LOAD = [str(x) for x in list(p.glob("*_load"))][0]
    WD_DEGLYCO = [str(x) for x in list(p.glob("*_deglyco"))][0]
    WD_STD = [str(x) for x in list(p.glob("*_std"))][0]
    WD_DEDUPL = [str(x) for x in list(p.glob("*_dedupl"))][0]

    # get dfs
    df_load = get_df_load(WD_LOAD)
    df_deglyco = get_df_deglyco(WD_DEGLYCO)
    df_std_passed = get_df_std_passed(WD_STD)
    df_std_filtered = get_df_std_filtered(WD_STD)
    df_std_error = get_df_std_error(WD_STD)
    df_dedupl = get_df_dedupl(WD_DEDUPL)

    # get total of molecules in input
    num_mols_tot = df_load['Count'].sum()
    # count not loaded molecules as well in deglyco
    num_mols_deglyco_na = df_load[df_load['Category'] == 'cannot_load']['Count']
    df_deglyco = pd.concat([df_deglyco, pd.DataFrame({'Category': ['NA'], 'Count': [num_mols_deglyco_na]})])
    df_deglyco.reset_index(inplace=True, drop=True)
    df_deglyco['Count'] = df_deglyco['Count'].astype(int)
    df_deglyco['Perc_Status'] = df_deglyco['Count'].map(lambda x: f"{x/num_mols_tot:.2%}")
    logger.info(f"PREP -- RESULTS FOR DEGLYCOSYLATION:\n\n{df_deglyco}\n")

    # gather all filtered molecules
    df_dedupl_dupl = df_dedupl[df_dedupl['Category'] == 'duplicate']
    num_dedupl_dupl = df_dedupl_dupl['Count'].sum()
    df_std_filtered = pd.concat([df_std_filtered, df_dedupl_dupl], sort=True)
    # count even unoccurred cases in df_std_filtered
    filters = ['empty', 'hac', 'molweight', 'nrings', 'medchem', 'timeout', 'duplicate']
    df_std_filtered.set_index('Category', inplace=True)
    df_std_filtered = df_std_filtered.reindex(filters)
    df_std_filtered.reset_index(inplace=True)
    df_std_filtered.fillna(0, inplace=True)
    df_std_filtered['Count'] = df_std_filtered['Count'].astype(int)
    df_std_filtered['Perc_Status'] = df_std_filtered['Count'].map(lambda x: f"{x/num_mols_tot:.2%}")
    logger.info(f"PREP -- RESULTS FOR STD_FILTERED:\n\n{df_std_filtered}\n")

    # gather all molecules that raised an error
    df_std_error = pd.concat([df_std_error, df_load[df_load['Category'] == 'cannot_load']], sort=True)
    # count even unoccurred cases in df_std_error
    errors = ['cannot_load', 'initiate_mol', 'disconnect_metal', 'sanitize', 'remove_isotopes', 'normalize', 'uncharge', 'canonicalize', 'remove_stereo']
    df_std_error.set_index('Category', inplace=True)
    df_std_error = df_std_error.reindex(errors)
    df_std_error.reset_index(inplace=True)
    df_std_error.fillna(0, inplace=True)
    df_std_error['Count'] = df_std_error['Count'].astype(int)
    df_std_error['Perc_Status'] = df_std_error['Count'].map(lambda x: f"{x/num_mols_tot:.2%}")
    logger.info(f"PREP -- RESULTS FOR STD_ERRORS:\n\n{df_std_error}\n")

    # general count for passed/filtered/errors
    num_tot_filtered = df_std_filtered['Count'].sum()
    num_tot_passed = df_std_passed['Count'].sum() - num_dedupl_dupl  # dedupl happens after std, so std contains passsed mols that get filtered
    num_tot_errors = df_std_error['Count'].sum()
    df_std_all = pd.DataFrame({'Category': ['passed', 'filtered', 'errors'], 'Count': [num_tot_passed, num_tot_filtered, num_tot_errors]})
    df_std_all['Perc_Status'] = df_std_all['Count'].map(lambda x: f"{x/num_mols_tot:.2%}")
    logger.info(f"PREP -- RESULTS FOR STD_ALL:\n\n{df_std_all}\n")

    return (df_std_all, df_std_filtered, df_std_error, df_deglyco)


def get_df_murcko(WD: Path) -> DataFrame:
    """Get a DF summarizing the results of the Murcko Scaffold extract step.
    It looks into 2 subfolders:
        -

    :param WD: the main directory of the dataset data (i.e. 'natural/dnp/data')
    :return: a DF summarizing results of the murcko scaffold extraction step
    """
    pass


def get_df_subset(WD: Path) -> DataFrame:
    """Get a DF summarizing the results of the subset step.

    At the moment only one subset is recognized (i.e. 'subset' subfolder in WD).


    :param WD: the main directory of the dataset data (i.e. 'natural/dnp/data')
    :return: a DF summarizing results of the murcko subset step
    """
    logger.info("SUB -- COMPUTING RESULTS FOR SUBSET")
    if not isinstance(WD, Path):
        WD = Path(WD)
    # parse results before fragment search
    WD_SUBSET = [str(x) for x in list(WD.glob("*_subset"))][0]    # get latest step
    pattern = ".*([0-9]{3})?.log"
    chunks = _get_chunks(f"{WD_SUBSET}/log", pattern)
    num_passed = 0
    num_filtered = 0
    num_total = 0
    for c in chunks:
        df = pd.read_csv(c, sep="@", header=None)  # char not found in the log file so we can extract all lines as one column
        records = df[df[0].str.contains("NUMBER OF REMAINING RECORDS IN SUBSET")].iloc[0][0].split()[-1].split('/')
        passed = int(records[0])
        total = int(records[1])
        num_filtered += total - passed
        num_passed += passed
        num_total += total
    # create a dataframe with counts
    df = pd.DataFrame({'Category': ['passed', 'filtered'], 'Count': [num_passed, num_filtered]})
    df['Perc_Status'] = df['Count'].map(lambda x: f"{x/num_total:.2%}")
    logger.info(f"SUB -- RESULTS FOR SUBSETTING {num_total:,} MOLECULES:\n\n{df}\n")

    return df


def get_dfs_fs(WD: Path) -> Tuple[DataFrame]:
    """Get a list of DFs summarizing the Fragment Search step.

    1. DF_NHITS is the summary of the number of fragments found per molecule
    2. DF_NHITS_U is the summary of the number of unique fragments found per molecule
    3. DF_FRAG_RATIO lists for each molecule the proportion of it that is matched with natural fragments
    4. DF_FRAG_RATIO_U lists for each molecule the proportion of it that is matched with natural fragments, counting each type of fragment only once
    5. DF_TOP_FRAGS is the list of all fragments found, ranked by decreasing occurrence
    6. DF_TOP_FRAGS_U is the list of all fragments found, ranked by decreasing occurrence,  but each type of fragment is counted only once per molecule,

    :param WD: the main directory of the dataset data (i.e. 'natural/dnp/data')
    :return: a list of DFs of interest: [DF_NHITS, DF_NHITS_U, DF_FRAG_RATIO, DF_FRAG_RATIO_U, DF_TOP_FRAGS, DF_TOP_FRAGS_U]
    """
    logger.info("FS -- COMPUTING RESULTS FOR FRAGMENT SEARCH")
    if not isinstance(WD, Path):
        WD = Path(WD)
    # parse results before fragment search
    logger.info("FS -- COMPUTING TOTAL NUMBER OF MOLECULES BEFORE FRAGMENT SEARCH")
    WD_FSEARCH = [str(x) for x in list(WD.glob("*_fsearch"))][0]    # get latest step
    fsearch_dirname = Path(WD_FSEARCH).name
    pstep_direname = list(WD.glob(f"{str(int(fsearch_dirname.split('_')[0]) - 1).zfill(2)}_*"))[0].name
    logger.info(f"FS -- COUNTING MOLECULES IN PREVIOUS STEP='{pstep_direname}'")

    # previous step to assess how many (1 row = 1 mol!)
    chunks_pstep = sorted([str(x) for x in list(WD.glob(f"{pstep_direname}/data/*_[0-9][0-9][0-9]_*.csv.gz"))])
    n_tot = sum([len(pd.read_csv(x, sep='|', compression='gzip')) for x in chunks_pstep])
    logger.info(f"FS -- TOTAL NUMBER OF MOLECULES IN PREVIOUS STEP: {n_tot:,}")
    logger.info(f"FS -- NOW INVESTIGATING RESULTS IN WD='{fsearch_dirname}'")

    # parse chunks
    chunks_fsearch = sorted([str(x) for x in list(WD.glob(f"{fsearch_dirname}/data/*_[0-9][0-9][0-9]_*.csv.gz"))])

    logger.info(f"FS -- STARTING CHUNK ITERATION...")
    n_fs_nhits_tot = 0
    n_fs_mols_tot = 0
    n_fs_nhits_tot_u = 0
    dfs_fs_nhits = []
    dfs_top_frags = []
    dfs_frag_ratio = []
    dfs_fs_nhits_u = []
    dfs_top_frags_u = []
    for x in chunks_fsearch:  # chunk-wise iteration to save up memory
        df_fsearch = load.file(x)
        df_fsearch['_aidxf'] = df_fsearch['_aidxf'].map(list)
        df_smiles = df_fsearch[['idf', 'mol_frag']].drop_duplicates(subset=['idf'])
        df_smiles['smiles_frag'] = df_smiles['mol_frag'].map(Chem.MolToSmiles)
        df_smiles = df_smiles[['idf', 'smiles_frag', 'mol_frag']]
        df_fsearch['hac_mol'] = df_fsearch['mol'].map(lambda x: x.GetNumAtoms())

        # fragment hits per mol
        groups = df_fsearch.groupby('idm')
        n_fs_nhits_tot += len(df_fsearch.index)
        n_fs_mols_tot += len(groups)  # there are no duplicate ids at this point, each chunk contain unique and non-redundant mols
        df_fs_nhits = groups.count()[['idf']].reset_index().rename({'idf': 'NumFrags'}, axis=1).groupby('NumFrags').count().rename({'idm': 'Count'}, axis=1).reset_index()  # this counts >0 hits
        dfs_fs_nhits.append(df_fs_nhits)

        # top fragments
        df_top_frags = df_fsearch[['idf', 'idm']].groupby('idf').count().reset_index().rename({'idm': 'Count'}, axis=1).sort_values('Count', ascending=False).reset_index(drop=True)
        dfs_top_frags.append(df_top_frags)

        # fragment ratio per molecule (to do last)
        df_fsearch['hac_mol'] = df_fsearch['mol'].map(lambda x: x.GetNumAtoms())
        df_frag_ratio = df_fsearch.copy()
        df_frag_ratio = groups.agg({'_aidxf': 'sum', 'hac_mol': 'first'}).reset_index()
        df_frag_ratio['_aidxf'] = df_frag_ratio['_aidxf'].map(lambda x: len(set(x)))
        df_frag_ratio.rename({'_aidxf': 'hac_frags'}, axis=1, inplace=True)
        df_frag_ratio['frag_ratio'] = df_frag_ratio['hac_frags'] / df_frag_ratio['hac_mol']
        dfs_frag_ratio.append(df_frag_ratio)

        # unique fragments per molecule only
        df_fsearch_u = df_fsearch.groupby(['idm', 'idf']).first().reset_index()
        n_fs_nhits_tot_u += len(df_fsearch_u.index)

        # unique fragment hits per mol
        df_fs_nhits_u = df_fsearch_u.groupby('idm').count()[['idf']].reset_index().rename({'idf': 'NumFrags'}, axis=1).groupby('NumFrags').count().rename({'idm': 'Count'}, axis=1).reset_index()  # this counts >0 hits
        dfs_fs_nhits_u.append(df_fs_nhits_u)

        # top unique fragments
        df_top_frags_u = df_fsearch_u[['idf', 'idm']].groupby('idf').count().reset_index().rename({'idm': 'Count'}, axis=1).sort_values('Count', ascending=False).reset_index(drop=True)
        dfs_top_frags_u.append(df_top_frags_u)

    logger.info(f"FS -- COMPLETED CHUNK ITERATION")
    logger.info(f"FS -- TOTAL NUMBER OF FRAGMENTS HITS={n_fs_nhits_tot:,}")
    logger.info(f"FS -- TOTAL NUMBER OF MOLECULES={n_fs_mols_tot:,}")
    logger.info(f"FS -- AVERAGE NUMBER OF FRAGMENTS PER MOL={n_fs_nhits_tot/n_fs_mols_tot:,.2f}")

    # count used later for fragment hit counts
    n_fs_nhits_0 = n_tot - n_fs_mols_tot  # this counts nhits = 0

    # fragment hits per mol
    logger.info(f"FS -- INVESTIGATING THE NUMBER OF FRAGMENT HITS PER MOLECULE")
    df_fs_nhits = pd.concat(dfs_fs_nhits).groupby('NumFrags').sum().reset_index().sort_values('NumFrags')  # concatenate counts because they do not take any memory anymore
    df_fs_nhits = pd.concat([pd.DataFrame({'NumFrags': [0], 'Count': [n_fs_nhits_0]}), df_fs_nhits]).sort_values('NumFrags').reset_index(drop=True)  # now we have all nhits
    df_fs_nhits['Perc_Mols'] = df_fs_nhits['Count'].map(lambda x: f"{x / n_fs_mols_tot:.2%}")
    logger.info(f"FS -- RESULTS FOR THE NUMBER OF FRAGMENT HITS PER MOLECULE\n\n{df_fs_nhits}\n")

    # top fragments
    logger.info(f"FS -- INVESTIGATING THE TOP FRAGMENTS")
    df_top_frags = pd.concat(dfs_top_frags)
    df_top_frags = df_top_frags.groupby('idf').sum().reset_index().sort_values('Count', ascending=False).reset_index(drop=True)
    df_top_frags['Rank'] = df_top_frags.index + 1
    df_top_frags['Perc_FHits'] = df_top_frags['Count'].map(lambda x: f"{x / n_fs_nhits_tot:.2%}")
    df_top_frags = df_top_frags.merge(df_smiles, how='inner', on='idf')
    logger.info(f"FS -- TOTAL NUMBER OF FRAGMENTS={len(df_top_frags):,}")
    logger.info(f"FS -- RESULTS FOR THE TOP FRAGMENTS\n\n{df_top_frags.drop('mol_frag', axis=1)}\n")

    # fragment ratio per molecule
    logger.info(f"FS -- INVESTIGATING THE RATIO OF FRAGMENT PER MOLECULE")
    df_frag_ratio = pd.concat(dfs_frag_ratio)
    logger.info(f"FS -- RESULTS FOR THE RATIO OF FRAGMENT PER MOLECULE\n\n{df_frag_ratio}\n")

    # unique fragment hits per mol
    logger.info(f"FS -- NOW COMPUTING RESULTS WHEN COUNTING ONLY UNIQUE FRAGMENTS PER MOLECULE")
    logger.info(f"FS -- TOTAL NUMBER OF UNIQUE FRAGMENTS HITS={n_fs_nhits_tot_u:,}")
    logger.info(f"FS -- INVESTIGATING THE NUMBER OF UNIQUE FRAGMENT HITS PER MOLECULE")
    df_fs_nhits_u = pd.concat(dfs_fs_nhits_u).groupby('NumFrags').sum().reset_index().sort_values('NumFrags')  # concatenate counts because they do not take any memory anymore
    df_fs_nhits_u = pd.concat([pd.DataFrame({'NumFrags': [0], 'Count': [n_fs_nhits_0]}), df_fs_nhits_u]).sort_values('NumFrags').reset_index(drop=True)  # now we have all nhits
    df_fs_nhits_u['Perc_Mols'] = df_fs_nhits_u['Count'].map(lambda x: f"{x / n_fs_nhits_tot_u:.2%}")
    logger.info(f"FS -- TOTAL NUMBER OF UNIQUE FRAGMENT HITS={n_fs_nhits_tot_u:,}")
    logger.info(f"FS -- AVERAGE NUMBER OF UNIQUE FRAGMENTS PER MOL={n_fs_nhits_tot_u/n_fs_mols_tot:,.2f}")
    logger.info(f"FS -- RESULTS FOR THE NUMBER OF UNIQUE FRAGMENT HITS PER MOLECULE\n\n{df_fs_nhits_u}\n")

    # top unique fragments
    logger.info(f"FS -- INVESTIGATING THE TOP UNIQUE FRAGMENTS")
    df_top_frags_u = pd.concat(dfs_top_frags_u)
    df_top_frags_u = df_top_frags_u.groupby('idf').sum().reset_index().sort_values('Count', ascending=False).reset_index(drop=True)
    df_top_frags_u['Rank'] = df_top_frags_u.index + 1
    df_top_frags_u['Perc_FHits'] = df_top_frags_u['Count'].map(lambda x: f"{x / n_fs_nhits_tot_u:.2%}")
    df_top_frags_u = df_top_frags_u.merge(df_smiles, how='inner', on='idf')
    logger.info(f"FS -- RESULTS FOR THE TOP UNIQUE FRAGMENTS\n\n{df_top_frags_u.drop('mol_frag', axis=1)}\n")


    # no unique fragment ratio because I cannot think of a good way to decide what atom index to retain

    return (df_fs_nhits, df_fs_nhits_u, df_frag_ratio, df_top_frags, df_top_frags_u)



def _prep_df_fcc(df_fcc: DataFrame, frags_file: str, compute_df_ffo: bool = False) -> tuple:
    """Prepare the df for fcc analysis. Used once for the fcc step and once for the fmap step.

    :param df_fcc: the fcc DF to process
    :return: a tuple of DFs: (DF_FCC, DF_FC)
    """
    # define categories
    categories = fragment_combination.get_fragment_combination_categories()
    # process DF
    df_fcc.rename({'abbrev': 'Category'}, axis=1, inplace=True)
    num_mols_tot = len(df_fcc.groupby('idm'))
    logger.info(f"FCC -- NUMBER OF MOLECULES IN TOTAL: {num_mols_tot:,}")
    n_fcc_tot = len(df_fcc.index)
    logger.info(f"FCC -- COUNTING FALSE POSITIVES (OVERLAPS)")
    df_fcc_ffo = df_fcc[df_fcc['Category'] == 'ffo'][['idm', 'idf1', 'idf2']]
    df_fcc = df_fcc[df_fcc['Category'] != 'ffo'][['idf1', 'idf2', 'Category', 'idm']]
    num_mols_tot_tp = len(df_fcc.groupby('idm'))
    logger.info(f"FCC -- NUMBER OF MOLECULES IN TP: {num_mols_tot_tp:,}")
    logger.info(f"FCC -- NUMBER OF MOLECULES IN FP: {num_mols_tot - num_mols_tot_tp:,}")

    n_fcc_tp = len(df_fcc.index)
    n_fcc_fp = n_fcc_tot - n_fcc_tp
    logger.info(f"FCC -- NUMBER OF FRAGMENT COMBINATIONS: {n_fcc_tp:,} TP +  {n_fcc_fp:,} FP = {n_fcc_tot:,} IN TOTAL")


    # number of overlaps per molecule
    if compute_df_ffo:
        logger.info(f"FCC -- INVESTIGATING OVERLAPS")
        df_fcc_ffo = df_fcc_ffo.groupby('idm').count().rename({'idf1' : 'NumOverlaps'}, axis=1).groupby('NumOverlaps').count().reset_index().rename({'idf2': 'Count'}, axis=1)
        df_fcc_ffo = pd.concat([DataFrame({'NumOverlaps': [0], 'Count': [num_mols_tot_tp]}), df_fcc_ffo]).reset_index(drop=True)
        df_fcc_ffo['Perc_Mols'] = df_fcc_ffo['Count'].map(lambda x: f"{x / num_mols_tot:.2%}")
        logger.info(f"FCC -- RESULTS: NUMBER OF OVERLAPS PER MOLECULE:\n\n{df_fcc_ffo}\n")

    # fragment combination categories counts
    logger.info(f"FCC -- FOCUSING ONLY ON THE {n_fcc_tp:,} TP")
    df_fcc_fcc_default = pd.DataFrame({'Category': categories, 'Count': [0] * len(categories)})
    df_fcc_fcc = df_fcc[['Category', 'idm']].groupby('Category').count().rename({'idm': 'Count'}, axis=1).reset_index()
    df_fcc_fcc = pd.concat([df_fcc_fcc, df_fcc_fcc_default]).groupby('Category').sum().T
    df_fcc_fcc = df_fcc_fcc[categories]

    df_fcc_fcc = df_fcc_fcc.T.reset_index().rename({'index': 'Category'}, axis=1)
    df_fcc_fcc['Perc'] = df_fcc_fcc['Count'].map(lambda x: f"{x / n_fcc_tp:.2%}")
    logger.info(f"FCC -- RESULTS FOR FCC COUNTS\n\n{df_fcc_fcc}\n")

    # top 10 fragment combinations
    df_fcc_fc = df_fcc[['idf1', 'idf2', 'Category', 'idm']].groupby(['idf1', 'idf2', 'Category']).count().rename({'idm': 'Count'}, axis=1).sort_values('Count', ascending=False).reset_index()
    n_fc_u = len(df_fcc_fc.index)
    logger.info(f"FCC -- NUMBER OF UNIQUE FC: {n_fc_u:,}")
    logger.info(f"FCC -- COMPUTING THE TOP FRAGMENT COMBINATIONS (FC)")
    df_fcc_fc['Perc'] = df_fcc_fc['Count'].map(lambda x: f"{x / n_fcc_tp:.2%}")
    df_fcc_fc['fc'] = df_fcc_fc['idf1'].astype(str) + '[' + df_fcc_fc['Category'] + ']' + df_fcc_fc['idf2'].astype(str)
    df_fcc_fc_top10 = df_fcc_fc.head(10)  # top 10 only
    n_fcc_tot_top10 = df_fcc_fc_top10['Count'].sum()
    logger.info(f"FCC -- TOP 10 FC: {n_fcc_tot_top10:,} ENTRIES ({n_fcc_tot_top10/n_fcc_tp:.2%})")
    # add fragments to the mix
    df_frags = load.file(frags_file)
    df_frags['idm'] = df_frags['idm'].astype(str)
    df_fcc_fc['idf1'] = df_fcc_fc['idf1'].astype(str)
    df_fcc_fc['idf2'] = df_fcc_fc['idf2'].astype(str)
    df_fcc_fc = pd.merge(df_fcc_fc, df_frags[['idm', 'mol']], left_on='idf1', right_on='idm', how='inner').rename({'mol': 'mol_frag_1'}, axis=1).drop('idm', axis=1)
    df_fcc_fc = pd.merge(df_fcc_fc, df_frags[['idm', 'mol']], left_on='idf2', right_on='idm', how='inner').rename({'mol': 'mol_frag_2'}, axis=1).drop('idm', axis=1).sort_values('Count', ascending=False).reset_index(drop=True)
    df_fcc_fc['Rank'] = df_fcc_fc.index + 1

    # display smiles in the log
    df_fcc_fcc_smiles = df_fcc_fc.copy()
    df_fcc_fcc_smiles['mol_frag_1'] = df_fcc_fcc_smiles['mol_frag_1'].map(Chem.MolToSmiles)
    df_fcc_fcc_smiles['mol_frag_2'] = df_fcc_fcc_smiles['mol_frag_2'].map(Chem.MolToSmiles)
    logger.info(f"FCC -- RESULTS FOR TOP 10 FC\n\n{df_fcc_fcc_smiles.head(10)}\n")

    if compute_df_ffo:
        return (df_fcc_fcc, df_fcc_fc, df_fcc_ffo)
    else:
        return (df_fcc_fcc, df_fcc_fc)


def get_dfs_fcc(WD: Path, frags_file: str, compute_df_ffo=False) -> Tuple[DataFrame]:
    """Get a list of DFs summarizing the Fragment Combination Classification step.

    - DF_FCC is the summary of the number of FCC found in the dataset for each category
    - DF_TOP_FC is the list of all Fragment Combinations found, ranked by the decreasing occurrence

    :param WD: the main directory of the dataset data (i.e. 'natural/dnp/data'). Can also be a string.
    :param frags_file: the fragment file used for fragment search (i.e. 'fragments/crms/data/07_gen2D/data/cmrs_gen2D.csv.gz')

    :return: a list of DFs of interest: [DF_FCC, DF_TOP_FC]
    """
    if not isinstance(WD, Path):
        WD = Path(WD)
    # check if col with activity in dataframe, if so return top10 most active fcc
    logger.info("FCC -- COMPUTING RESULTS FOR FRAGMENT COMBINATION CLASSIFICATION (FCC)")

    # parse results before fragment search
    logger.info("FCC -- COMPUTING TOTAL NUMBER OF MOLECULES BEFORE FCC")
    # get latest step
    fcc_dirname = list(WD.glob("*_fcc"))[0].name
    pstep_direname = list(WD.glob(f"{str(int(fcc_dirname.split('_')[0]) - 1).zfill(2)}_*"))[0].name
    logger.info(f"FCC -- COUNTING MOLECULES IN PREVIOUS STEP='{pstep_direname}'")

    # previous step to assess how many mols
    chunks_pstep = sorted([str(x) for x in list(WD.glob(f"{pstep_direname}/data/*_[0-9][0-9][0-9]_*.csv.gz"))])
    n_tot = sum([len(pd.read_csv(x, sep='|', compression='gzip').groupby('idm')) for x in chunks_pstep])
    logger.info(f"FCC -- TOTAL NUMBER OF MOLECULES IN PREVIOUS STEP: {n_tot:,}")
    logger.info(f"FCC -- INVESTIGATING RESULTS IN WD='{fcc_dirname}'")

    # parse chunks
    chunks_fcc = sorted([str(x) for x in list(WD.glob(f"{fcc_dirname}/data/*_[0-9][0-9][0-9]_*.csv.gz"))])
    df_fcc = pd.concat([load.file(x) for x in chunks_fcc])

    # data is processed in function

    return _prep_df_fcc(df_fcc, frags_file, compute_df_ffo=compute_df_ffo)


def get_df_fm(WD: Path, frags_file) -> DataFrame:
    """test
    """
    # define data
    WD_FMAP = [str(x) for x in list(WD.glob("*_fmap"))][0]    # get latest step
    pattern = ".*([0-9]{3})?.csv.gz"
    chunks = _get_chunks(f"{WD_FMAP}/data", pattern)
    categories = fragment_combination.get_fragment_combination_categories()

    # parse fragments data
    logger.info(f"FM -- RETRIEVING FRAGMENT DATA")
    df_frags = load.file(frags_file).rename({'idm': 'idf', 'mol': 'mol_frag'}, axis=1)
    df_frags['idf'] = df_frags['idf'].astype(str)
    df_frags['frag_hac'] = df_frags['mol_frag'].map(lambda x: x.GetNumAtoms())
    df_frags['smiles_frag'] = df_frags['mol_frag'].map(Chem.MolToSmiles)
    df_frags = df_frags[['idf', 'frag_hac', 'smiles_frag', 'mol_frag']]
    logger.info(f"FM -- RETRIEVED {len(df_frags.index)} FRAGMENTS")


    # initialize chunk iteration
    logger.info(f"FS -- STARTING CHUNK ITERATION...")
    num_tot_fm_map = 0
    num_tot_mol_map = 0
    n_fm_nhits_tot = 0
    n_fm_fcc = 0
    dfs_fm_nfmpermol = []
    dfs_fm_nhits = []
    dfs_fm_top_frags = []
    dfs_fm_frag_ratio = []
    dfs_fm_nhits_u = []
    dfs_fm_top_frags_u = []
    dfs_fm_fcc = []
    dfs_fm_fc = []

    # chunk iteration
    logger.setLevel(logging.WARNING)  # function has loggings at info level that would flood the log file because of iteration
    for chunk in chunks:

        # retrieve data
        df_fm = load.file(chunk, decode=['_fmap']).sort_values(["idm", "nfrags"], ascending=True) # df_fm already sorted for the best examples per case
        num_tot_fm_map += len(df_fm.index)
        groups = df_fm[['idm', 'fmid', 'nfrags']].groupby('idm')
        num_tot_mol_map += len(groups)

        # number of fragments map per molecule
        df_fm_nfmpermol = groups.count().rename({'fmid' : 'NumFM'}, axis=1).groupby('NumFM').count().rename({'nfrags': 'Count'}, axis=1).reset_index()
        dfs_fm_nfmpermol.append(df_fm_nfmpermol)

        # extract the maps
        df_edges = pd.concat([nx.convert_matrix.to_pandas_edgelist(df_fm.iloc[i]["_fmap"]) for i in range(len(df_fm.index))]).rename({'source': 'idf1', 'target': 'idf2'}, axis=1)
        df_edges.rename({'source': 'idf1', 'target': 'idf2'}, axis=1)

        # fs analysis

        # initialization of a common df useful for fs analysis
        df_fm['frags'] = df_fm['frags'].map(lambda x: x.replace("'", '"'))
        df_fm['frags'] = df_fm['frags'].map(lambda x: json.loads(x))
        df_fm['aidxs'] = df_fm['_d_aidxs'].map(lambda x: [v for l in x.values() for v in l])  # extract all values from dict: list of tuples
        df_fm['aidxs'] = df_fm['aidxs'].map(lambda x: [v for l in x for v in l])   # flatten the list and identify only unique atom indices
        groups = df_fm.groupby('idm')
        df_fm_grouped = groups.agg({'frags': 'sum'}).reset_index(drop=True)  # concatenate lists in the same group
        df_fm_grouped['frags'] = df_fm_grouped['frags'].map(lambda x: list(set(x)))  # count each occurrence of a fragment
        df_fm_grouped['n_frags'] =  df_fm_grouped['frags'].map(lambda x: len(x))

        # fragment hits per mol
        df_fm_nhits = df_fm_grouped[['n_frags', 'frags']].groupby('n_frags').count().reset_index().rename({'frags': 'Count', 'n_frags': 'NumFrags'}, axis=1)
        df_fm_nhits['Perc_Mols'] = df_fm_nhits['Count'].map(lambda x: f"{x / num_tot_mol_map:.2%}")
        dfs_fm_nhits.append(df_fm_nhits)

        # top fragments

        # process fm data
        df_fm_top_frags = df_fm_grouped['frags'].apply(pd.Series).reset_index().melt(id_vars='index').dropna()[['index', 'value']]  # ungroup values by frag id in list
        df_fm_top_frags['value'] = df_fm_top_frags['value'].map(lambda x: x.split(':')[0])
        df_fm_top_frags = df_fm_top_frags.groupby('value').count().reset_index().rename({'value': 'idf', 'index': 'Count'}, axis=1).sort_values('Count', ascending=False).reset_index(drop=True)  # count and sort idfs
        dfs_fm_top_frags.append(df_fm_top_frags)
        n_fm_nhits_tot += df_fm_top_frags['Count'].sum()

        # fragment ratio per molecule
        df_fm_frag_ratio = groups.agg({'aidxs': 'sum', 'hac_mol': 'first'}).reset_index()  # concatenate all aidxs obtained previously
        df_fm_frag_ratio['hac_frags'] = df_fm_frag_ratio['aidxs'].map(lambda x: len(set(x)))  # the length of atom indices is the number of hac in fragments
        df_fm_frag_ratio['frag_ratio'] = df_fm_frag_ratio['hac_frags'] / df_fm_frag_ratio['hac_mol']
        df_fm_frag_ratio.drop('aidxs', axis=1, inplace=True)
        dfs_fm_frag_ratio.append(df_fm_frag_ratio)

        # unique fragment hits per mol
        df_fm_grouped['frags_u'] = df_fm_grouped['frags'].map(lambda x: list(set([v.split(':')[0] for v in x])))
        df_fm_grouped['n_frags_u'] =  df_fm_grouped['frags_u'].map(lambda x: len(x))
        df_fm_nhits_u = df_fm_grouped[['n_frags_u', 'frags_u']].groupby('n_frags_u').count().reset_index().rename({'frags_u': 'Count', 'n_frags_u': 'NumFrags'}, axis=1)
        dfs_fm_nhits_u.append(df_fm_nhits_u)

        # top unique fragments
        df_fm_top_frags_u = df_fm_grouped['frags_u'].apply(pd.Series).reset_index().melt(id_vars='index').dropna()[['index', 'value']]  # ungroup values by frag id in list
        df_fm_top_frags_u = df_fm_top_frags_u.groupby('value').count().reset_index().rename({'value': 'idf', 'index': 'Count'}, axis=1).sort_values('Count', ascending=False).reset_index(drop=True)  # count and sort idfs
        dfs_fm_top_frags_u.append(df_fm_top_frags_u)

        # fcc analysis
        n_fm_fcc += len(df_edges.index)
        df_fm_fcc, df_fm_fc = _prep_df_fcc(df_edges, frags_file, compute_df_ffo=False)  # no ffo in maps because they got splitted up during fmapping
        dfs_fm_fcc.append(df_fm_fcc)
        dfs_fm_fc.append(df_fm_fc)

    logger.setLevel(logging.INFO)  # ok now back to normal
    logger.info(f"FS -- COMPLETED CHUNK ITERATION...")

    # fm_nfmpermol
    logger.info(f"FM -- RESULTS FOR THE NUMBER OF FM PER MOLECULE")
    logger.info(f"FM -- TOTAL NUMBER OF FRAGMENT MAPS: {num_tot_fm_map:,d}")
    logger.info(f"FM -- TOTAL NUMBER OF MOLECULES: {num_tot_mol_map:,d}")
    df_fm_nfmpermol = pd.concat(dfs_fm_nfmpermol).groupby('NumFM').sum().reset_index()
    df_fm_nfmpermol['Perc_Mols'] = df_fm_nfmpermol['Count'].map(lambda x: f"{x / num_tot_mol_map:.2%}")
    logger.info(f"FM -- RESULTS FOR THE NUMBER OF FM PER MOLECULE:\n\n{df_fm_nfmpermol}\n")

    # fm_nhits
    logger.info(f"FM -- INVESTIGATING FOR THE NUMBER OF FRAGMENT HITS PER MOLECULE")
    df_fm_nhits = pd.concat(dfs_fm_nhits).groupby('NumFrags').sum().reset_index().sort_values('NumFrags').reset_index(drop=True)
    df_fm_nhits['Perc_Mols'] = df_fm_nhits['Count'].map(lambda x: f"{x / num_tot_mol_map:.2%}")
    logger.info(f"FM -- RESULTS FOR THE NUMBER OF FRAGMENT HITS PER MOLECULE\n\n{df_fm_nhits}\n")

    # fm_top_frags
    logger.info(f"FM -- INVESTIGATING FOR THE TOP FRAGMENTS")
    df_fm_top_frags = pd.concat(dfs_fm_top_frags).groupby('idf').sum().reset_index().sort_values('Count', ascending=False)
    df_fm_top_frags['Rank'] = df_fm_top_frags.index + 1
    logger.info(f"FM -- TOTAL NUMBER OF FRAGMENT HITS={n_fm_nhits_tot:,}")
    df_fm_top_frags['Perc_FHits'] = df_fm_top_frags['Count'].map(lambda x: f"{x / n_fm_nhits_tot:.2%}")
    df_fm_top_frags['idf'] = df_fm_top_frags['idf'].astype(str)
    df_fm_top_frags = pd.merge(df_fm_top_frags, df_frags, on='idf', how='inner').drop('frag_hac', axis=1)
    logger.info(f"FM -- RESULTS FOR THE TOP FRAGMENTS\n\n{df_fm_top_frags.drop('mol_frag', axis=1)}\n")

    # fm_frag_ratio
    logger.info(f"FM -- INVESTIGATING FOR THE RATIO OF FRAGMENT PER MOLECULE")
    df_fm_frag_ratio = pd.concat(dfs_fm_frag_ratio)
    logger.info(f"FM -- RESULTS FOR THE RATIO OF FRAGMENT PER MOLECULE\n\n{df_fm_frag_ratio}\n")

    # fm_nhits_u
    logger.info(f"FM -- INVESTIGATING THE NUMBER OF UNIQUE FRAGMENT HITS PER MOLECULE")
    df_fm_nhits_u = pd.concat(dfs_fm_nhits_u).groupby('NumFrags').sum().reset_index().sort_values('NumFrags').reset_index(drop=True)
    df_fm_nhits_u['Perc_Mols'] = df_fm_nhits_u['Count'].map(lambda x: f"{x / num_tot_mol_map:.2%}")
    logger.info(f"FM -- RESULTS FOR THE NUMBER OF UNIQUE FRAGMENT HITS PER MOLECULE\n\n{df_fm_nhits_u}\n")

    # fm_top_frags_u
    logger.info(f"FM -- INVESTIGATING THE TOP UNIQUE FRAGMENTS")
    df_fm_top_frags_u = pd.concat(dfs_fm_top_frags_u).groupby('idf').sum().reset_index().sort_values('Count', ascending=False)
    df_fm_top_frags_u['Rank'] = df_fm_top_frags_u.index + 1
    df_fm_top_frags_u['Perc_FHits'] = df_fm_top_frags_u['Count'].map(lambda x: f"{x / n_fm_nhits_tot:.2%}")
    df_fm_top_frags_u['idf'] = df_fm_top_frags_u['idf'].astype(str)
    df_fm_top_frags_u = pd.merge(df_fm_top_frags_u, df_frags, on='idf', how='inner').drop('frag_hac', axis=1)
    logger.info(f"FM -- RESULTS FOR THE TOP UNIQUE FRAGMENTS\n\n{df_fm_top_frags_u.drop('mol_frag', axis=1)}\n")

    # fm_fcc
    logger.info(f"FM -- INVESTIGATING THE FCC COUNTS")
    df_fm_fcc = pd.concat(dfs_fm_fcc).groupby('Category').sum().T  # use transposition for sorting cols in predefined order
    df_fm_fcc = df_fm_fcc[categories].T.reset_index()  # once it is all good, transpose again to get rows in expected order
    df_fm_fcc['Perc'] = df_fm_fcc['Count'].map(lambda x: f"{x / n_fm_fcc:.2%}")
    logger.info(f"FM -- RESULTS FOR THE FCC COUNTS\n\n{df_fm_fcc}\n")

    # fm_fc
    logger.info(f"FM -- INVESTIGATING THE TOP FRAGMENT COMBINATIONS")
    df_fm_fc = pd.concat(dfs_fm_fc)
    df_fm_fc['mol_frag_1'] = df_fm_fc['mol_frag_1'].map(Chem.MolToSmiles)
    df_fm_fc['mol_frag_2'] = df_fm_fc['mol_frag_2'].map(Chem.MolToSmiles)
    df_fm_fc = df_fm_fc.groupby(['idf1', 'idf2', 'Category', 'fc', 'mol_frag_1', 'mol_frag_2']).sum().reset_index().sort_values('Count', ascending=False).reset_index(drop=True)
    df_fm_fc['Perc'] = df_fm_fc['Count'].map(lambda x: f"{x / n_fm_fcc:.2%}")
    df_fm_fc['Rank'] = df_fm_fc.index + 1
    logger.info(f"FM -- RESULT FOR THE TOP FRAGMENT COMBINATIONS\n\n{df_fm_fc}\n")


    return (df_fm_nhits, df_fm_nhits_u, df_fm_frag_ratio, df_fm_top_frags, df_fm_top_frags_u, df_fm_fcc, df_fm_fc, df_fm_nfmpermol)


def get_df_pnp(WD: Path) -> DataFrame:
    """Get a DF summarizing the results of the pnp step.

    At the moment only one subset is recognized (i.e. 'subset' subfolder in WD).


    :param WD: the main directory of the dataset data (i.e. 'natural/dnp/data')
    :return: a DF summarizing results of the murcko subset step
    """
    logger.info("PNP -- COMPUTING RESULTS FOR PNP")
    if not isinstance(WD, Path):
        WD = Path(WD)
    # parse results before fragment search
    WD_PNP = [str(x) for x in list(WD.glob("*_pnp"))][0]    # get latest step
    pattern = ".*([0-9]{3})?.log"
    chunks = _get_chunks(f"{WD_PNP}/log", pattern)
    num_pnp = 0
    num_non_pnp = 0
    num_total = 0
    for c in chunks:
        df = pd.read_csv(c, sep="@", header=None)  # char not found in the log file so we can extract all lines as one column
        pnp = int(df[df[0].str.contains("LIST OF PNPs")].iloc[0][0].split()[-1].replace('(', '').replace(')', ''))
        non_pnp = int(df[df[0].str.contains("LIST OF NON-PNPs")].iloc[0][0].split()[-1].replace('(', '').replace(')', ''))
        num_pnp += pnp
        num_non_pnp += non_pnp
    num_total = num_pnp + num_non_pnp
    # create a dataframe with counts
    df = pd.DataFrame({'Category': ['PNP', 'Non-PNP'], 'Count': [num_pnp, num_non_pnp]})
    df['Perc_Mols'] = df['Count'].map(lambda x: f"{x/num_total:.2%}")
    logger.info(f"PNP -- RESULTS FOR LABELLING PNPs IN {num_total:,} MOLECULES:\n\n{df}\n")

    return df


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ BEGIN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #


if __name__ == '__main__':

    d0 = datetime.now()
    parser = argparse.ArgumentParser(description="Compute all required files for analyzing FCC results", formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('wd', type=str, default='./', help="Working directory from the job to analyse (root folder, where the smk file is)")
    parser.add_argument('-f', '--frags', type=str, default=None, help="Fragment file used for substructure search")
    parser.add_argument('-d', '--dataset', type=str, default=None, help="Dataset name for using in the csv/png outputs in the report folder.")
    parser.add_argument('-p', '--prefix', type=str, default=None, help="Prefix used for output files in the data/log folders.")
    parser.add_argument('-c', '--color', type=str, default='', help="Prefix used for output files in the data/log folders.")
    parser.add_argument('--csv', type=str, default=False, help="Generate only CSV output files")
    parser.add_argument('--log', type=str, default='INFO', help="Specify level of logging. Possible values are: CRITICAL, ERROR, WARNING, INFO, DEBUG.")
    args = parser.parse_args()


    # logging
    logger = utils._configure_logger(args.log)
    logger.info("PLEASE RUN THIS COMMAND FROM THE DATA FOLDER OF THE DATASET TO PARSE (i.e. natural/dnp/data)")
    logger.info("RUNNING PLOT_PIPELINE_RESULTS")
    warnings.filterwarnings('ignore', category=pd.io.pytables.PerformanceWarning)  # if None is returned instead of a molecule, do not complain about mixed types
    pd.options.mode.chained_assignment = None  # disable pd.io.pytables.SettingWithCopyWarning
    pad = 40
    lg = RDLogger.logger()
    lg.setLevel(RDLogger.INFO)

    logger.info("LIBRARY VERSIONS:")
    logger.info("rdkit".ljust(pad) + f"{rdkit.__version__}")
    logger.info("pandas".ljust(pad) + f"{pd.__version__}")
    logger.info("npfc".ljust(pad) + f"{npfc.__version__}")
    logger.info("seaborn".ljust(pad) + f"{sns.__version__}")

    logger.info("ARGUMENTS:")
    logger.info("WD".ljust(pad) + f"{args.wd}")
    logger.info("DATASET".ljust(pad) + f"{args.dataset}")
    logger.info("PREFIX".ljust(pad) + f"{args.prefix}")
    logger.info("FRAGMENT FILE".ljust(pad) + f"{args.frags}")
    logger.info("COLOR".ljust(pad) + f"{args.color}")
    logger.info("COMPUTE CSV ONLY".ljust(pad) + f"{args.csv}")

    # I like to have my very own palette of colors but have a hard time remembering hexadecimal codes
    if args.color == 'gray':
        color = '#808080'
    elif args.color == 'green':
        color = '#2CA02C'
    elif args.color == 'blue':
        color = '#378EBF'
    elif args.color == 'red':
        color = '#EB5763'
    else:
        color = args.color

    # output folder
    p = Path(args.wd)
    output_folder = Path(str(p) + "/report")
    logger.info("OUTPUT_FOLDER".ljust(pad) + f"{output_folder}")
    if not output_folder.exists():
        output_folder.mkdir(parents=True, exist_ok=True)

    if args.frags is None:
        logger.error("WARNING! FRAGS IS NOT DEFINED, SO NO ANALYSIS FURTHER THAN POSTPROCESSING OR MURCKO")
    else:
        utils.check_arg_input_file(args.frags)
        frags_subdir = '_'.join(Path(args.frags).stem.split('_')[:-1])  # remove last suffix from frags file, so works with both crms_gen2D and crms_nobz_gen2D

    steps = ['find_missing_chunks', 'preprocess']
    subfolders = _get_subfolders(p)['included']
    if any([sf.split('_')[1] == 'subset' for sf in subfolders]):
        steps += ['subset']
    if any([sf.split('_')[1] == 'fsearch' for sf in subfolders]):
        steps += ['fsearch']
    if any([sf.split('_')[1] == 'fcc' for sf in subfolders]):
        steps += ['fcc']
    if any([sf.split('_')[1] == 'fmap' for sf in subfolders]):
        steps += ['fmap']
    if any([sf.split('_')[1] == 'pnp' for sf in subfolders]):
        steps += ['pnp']


    logger.info(f"STEPS TO INVESTIGATE: {', '.join(steps)}")


    # begin iteration, better than sequential code because this makes it possible to easily go to next step if any missing chunk
    for step in steps:


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ MISSING CHUNKS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #


        if step == 'find_missing_chunks':
            _print_title("CHECKING FOR MISSING CHUNKS", logger, pad)
            errors_data, errors_log = find_missing_chunks(p)
            if not errors_data and not errors_log:
                logger.info("ALL EXPECTED CHUNKS WERE FOUND!")
            else:
                if errors_data:
                    logger.error("MISSING CHUNKS WERE IDENTIFIED IN DATA, ANY STEP INVOLVING BELOW STEPS WILL BE SKIPPED!")
                    logger.error("RESULTS: IDS OF IDENTIFIED MISSING CHUNKS\n" + '\n'.join(f"{k}: {', '.join(v)}" for k, v in errors_data.items()))

                if errors_log:
                    logger.error("MISSING CHUNKS WERE IDENTIFIED IN LOG, ANY STEP INVOLVING BELOW STEPS WILL BE SKIPPED!")
                    logger.error("RESULTS: IDS OF IDENTIFIED MISSING CHUNKS\n" + '\n'.join(f"{k}: {', '.join(v)}" for k, v in errors_log.items()))


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ PREP ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #

        elif step == 'preprocess':

            _print_title("CHECKING RESULTS FROM PREP", logger, pad)
            # check if data is complete
            subdir = ['load', 'deglyco', 'std', 'dedupl']
            if not _step_is_valid(subdir, errors_data, errors_log):
                logger.error(f"MISSING DATA OR LOG FOR PREP... SKIPPING THIS STEP!")
                continue

            # define csv outputs
            output_csv_prep_overview = f"{output_folder}/{args.dataset}_prep_overview.csv"
            output_csv_prep_filtered = f"{output_folder}/{args.dataset}_prep_filtered.csv"
            output_csv_prep_error = f"{output_folder}/{args.dataset}_prep_error.csv"
            output_csv_prep_deglyco = f"{output_folder}/{args.dataset}_prep_deglyco.csv"
            # define png outputs
            output_png_prep_overview = output_csv_prep_overview.replace('.csv', '.png')
            output_png_prep_filtered = output_csv_prep_filtered.replace('.csv', '.png')
            output_png_prep_error = output_csv_prep_error.replace('.csv', '.png')
            output_png_prep_deglyco = output_csv_prep_deglyco.replace('.csv', '.png')
            logger.info("PREP -- OUTPUT_CSV_PREP_OVERVIEW".ljust(pad) + f"{output_csv_prep_overview}")
            logger.info("PREP -- OUTPUT_CSV_PREP_FILTERED".ljust(pad) + f"{output_csv_prep_filtered}")
            logger.info("PREP -- OUTPUT_CSV_PREP_ERROR".ljust(pad) + f"{output_csv_prep_error}")
            logger.info("PREP -- OUTPUT_CSV_PREP_DEGLYCO".ljust(pad) + f"{output_csv_prep_deglyco}")
            logger.info("PREP -- OUTPUT PNG FILES HAVE THE SAME FILE NAMES AS OUTPUT CSV FILES")
            # retrieve data
            output_csv_files = [output_csv_prep_overview, output_csv_prep_filtered,
                                output_csv_prep_error, output_csv_prep_deglyco,
                                ]
            output_png_files = [output_png_prep_overview, output_png_prep_filtered,
                                output_png_prep_error, output_png_prep_deglyco,
                                ]
            if all([Path(x).exists() for x in output_png_files]):
                logger.info(f"PREP -- ALL OUTPUT PNG FILES ARE ALREADY AVAILABLE, NOTHING TO DO!")
            elif all([Path(x).exists() for x in output_csv_files]):
                logger.info(f"PREP -- PARSING OUTPUT CSV FILES INSTEAD OF COMPUTING THEM")
                df_prep_overview = load.file(output_csv_prep_overview)
                df_prep_filtered = load.file(output_csv_prep_filtered)
                df_prep_error = load.file(output_csv_prep_error)
                df_prep_deglyco = load.file(output_csv_prep_deglyco)
            else:
                logger.info(f"PREP -- COMPUTING OUTPUT CSV FILES")
                dfs_prep =  get_dfs_prep(p)
                for df_prep, output_csv_fs in zip(dfs_prep, output_csv_files):
                    save.file(df_prep, output_csv_fs)
                df_prep_overview = dfs_prep[0]
                df_prep_filtered =  dfs_prep[1]
                df_prep_error =  dfs_prep[2]
                df_prep_deglyco =  dfs_prep[3]

            # skip plots if computing only CSV output files
            if args.csv:
                continue
            # plot output_png_prep_overview
            if Path(output_png_prep_overview).exists():
                logger.info("PREP -- OUTPUT_PNG_PREP_OVERVIEW".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("PREP -- OUTPUT_PNG_PREP_OVERVIEW".ljust(pad) + f"COMPUTING...")
                save_barplot(df_prep_overview,
                             output_png_prep_overview,
                             'Category',
                             'Count',
                             f"Preprocesssing of Molecules in {args.dataset} - Overview",
                             x_label='Category',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Status',
                             fig_size=(12, 12),
                             )
            # plot output_png_prep_filtered
            if Path(output_png_prep_filtered).exists():
                logger.info("PREP -- OUTPUT_PNG_PREP_FILTERED".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("PREP -- OUTPUT_PNG_PREP_FILTERED".ljust(pad) + f"COMPUTING...")
                save_barplot(df_prep_filtered,
                             output_png_prep_filtered,
                             'Category',
                             'Count',
                             f"Preprocesssing of Molecules in {args.dataset} - Filtered",
                             x_label='Category',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Status',
                             )
            # plot output_png_prep_error
            if Path(output_png_prep_error).exists():
                logger.info("PREP -- OUTPUT_PNG_PREP_ERROR".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("PREP -- OUTPUT_PNG_PREP_ERROR".ljust(pad) + f"COMPUTING...")
                save_barplot(df_prep_error,
                             output_png_prep_error,
                             'Category',
                             'Count',
                             f"Preprocesssing of Molecules in {args.dataset} - Error",
                             x_label='Category',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Status',
                             )
            # plot output_png_prep_deglyco
            if Path(output_png_prep_deglyco).exists():
                logger.info("PREP -- OUTPUT_PNG_PREP_DEGLYCO".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("PREP -- OUTPUT_PNG_PREP_DEGLYCO".ljust(pad) + f"COMPUTING...")
                save_barplot(df_prep_deglyco,
                             output_png_prep_deglyco,
                             'Category',
                             'Count',
                             f"Preprocesssing of Molecules in {args.dataset} - Deglycosylation",
                             x_label='Category',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Status',
                             fig_size=(12, 12),
                             )


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBSET ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #


        elif step == 'subset':
            _print_title("CHECKING RESULTS FROM SUBSET", logger, pad)
            # check if data is complete
            subdir = 'subset'
            if not _step_is_valid(subdir, errors_data, errors_log):
                logger.error(f"MISSING DATA OR LOG FOR SUBSET... SKIPPING THIS STEP!")
                continue
            # define outputs
            output_csv_sub_subset = f"{output_folder}/{args.dataset}_sub_subset.csv"
            output_png_sub_subset = output_csv_sub_subset.replace('.csv', '.png')
            logger.info("SUB -- OUTPUT_CSV_SUB_SUBSET".ljust(pad) + f"{output_csv_sub_subset}")
            logger.info("SUB -- OUTPUT PNG FILES HAVE THE SAME FILE NAMES AS OUTPUT CSV FILES")
            # retrieve data
            output_csv_files = [output_csv_sub_subset]
            output_png_files = [output_png_sub_subset]
            if all([Path(x).exists() for x in output_png_files]):
                logger.info(f"SUB -- ALL OUTPUT PNG FILES ARE ALREADY AVAILABLE, NOTHING TO DO!")
            elif all([Path(x).exists() for x in output_csv_files]):
                logger.info(f"SUB -- PARSING OUTPUT CSV FILES INSTEAD OF COMPUTING THEM")
                df_sub_subset = load.file(output_csv_sub_subset)
            else:
                logger.info(f"SUB -- COMPUTING OUTPUT CSV FILES")
                df_sub_subset =  get_df_subset(p)
                save.file(df_sub_subset, output_csv_sub_subset)

            # skip plots if computing only CSV output files
            if args.csv:
                continue
            # plot output_png_sub_subset
            if Path(output_png_sub_subset).exists():
                logger.info("SUB -- OUTPUT_PNG_SUB_SUBSET".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("SUB -- OUTPUT_PNG_SUB_SUBSET".ljust(pad) + f"COMPUTING...")
                save_barplot(df_sub_subset,
                             output_png_sub_subset,
                             'Category',
                             'Count',
                             f"Creating a Synthetic Subset of {args.dataset}",
                             x_label='Category',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Status',
                             fig_size = (12, 12),
                             )


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ FS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #

        elif step == 'fsearch':

            _print_title("CHECKING RESULTS FROM FS", logger, pad)
            # check if data is complete
            subdir = 'sub'
            if not _step_is_valid(subdir, errors_data, errors_log):
                logger.error(f"MISSING DATA OR LOG FOR FSEARCH... SKIPPING THIS STEP!")
                continue
            # define outputs
            output_csv_fs_nfragpermol = f"{output_folder}/{args.dataset}_fs_nfragpermol.csv"
            output_csv_fs_nfragpermol_u = f"{output_folder}/{args.dataset}_fs_nfragpermol_u.csv"
            output_csv_fs_fragmolcov = f"{output_folder}/{args.dataset}_fs_fragmolcov.csv"
            output_csv_fs_top10frags = f"{output_folder}/{args.dataset}_fs_top10frags.csv"
            output_csv_fs_top10frags_u = f"{output_folder}/{args.dataset}_fs_top10frags_u.csv"
            output_png_fs_nfragpermol = output_csv_fs_nfragpermol.replace('.csv', '.png')
            output_png_fs_nfragpermol_u = output_csv_fs_nfragpermol_u.replace('.csv', '.png')
            output_png_fs_fragmolcov = output_csv_fs_fragmolcov.replace('.csv', '.png')
            output_png_fs_top10frags = output_csv_fs_top10frags.replace('.csv', '.png')
            output_png_fs_top10frags_u = output_csv_fs_top10frags_u.replace('.csv', '.png')
            logger.info("FS -- OUTPUT_CSV_FS_NFRAGPERMOL".ljust(pad) + f"{output_csv_fs_nfragpermol}")
            logger.info("FS -- OUTPUT_CSV_FS_NFRAGPERMOL_U".ljust(pad) + f"{output_csv_fs_nfragpermol_u}")
            logger.info("FS -- OUTPUT_CSV_FS_FRAGMOLCOV".ljust(pad) + f"{output_csv_fs_fragmolcov}")
            logger.info("FS -- OUTPUT_CSV_FS_TOP10FRAGS".ljust(pad) + f"{output_csv_fs_top10frags}")
            logger.info("FS -- OUTPUT_CSV_FS_TOP10FRAGS_U".ljust(pad) + f"{output_csv_fs_top10frags_u}")
            logger.info("FS -- OUTPUT PNG FILES HAVE THE SAME FILE NAMES AS OUTPUT CSV FILES")
            # retrieve data
            output_csv_files = [output_csv_fs_nfragpermol, output_csv_fs_nfragpermol_u,
                                output_csv_fs_fragmolcov,
                                output_csv_fs_top10frags, output_csv_fs_top10frags_u,
                                ]
            output_png_files = [output_png_fs_nfragpermol, output_png_fs_nfragpermol_u,
                                output_png_fs_fragmolcov,
                                output_png_fs_top10frags, output_png_fs_top10frags_u,
                                ]
            if all([Path(x).exists() for x in output_png_files]):
                logger.info(f"FS -- ALL OUTPUT PNG FILES ARE ALREADY AVAILABLE, NOTHING TO DO!")
            elif all([Path(x).exists() for x in output_csv_files]):
                logger.info(f"FS -- PARSING OUTPUT CSV FILES INSTEAD OF COMPUTING THEM")
                df_fs_nfragpermol = load.file(output_csv_fs_nfragpermol)
                df_fs_nfragpermol_u = load.file(output_csv_fs_nfragpermol_u)
                df_fs_fragmolcov = load.file(output_csv_fs_fragmolcov)
                df_fs_top10frags = load.file(output_csv_fs_top10frags).head(10)
                df_fs_top10frags_u = load.file(output_csv_fs_top10frags_u).head(10)
            else:
                logger.info(f"FS -- COMPUTING OUTPUT CSV FILES")
                dfs_fs =  get_dfs_fs(p)
                for df_fs, output_csv_fs in zip(dfs_fs, output_csv_files):
                    save.file(df_fs, output_csv_fs)
                df_fs_nfragpermol = dfs_fs[0]
                df_fs_nfragpermol_u =  dfs_fs[1]
                df_fs_fragmolcov =  dfs_fs[2]
                df_fs_top10frags =  dfs_fs[3].head(10)
                df_fs_top10frags_u =  dfs_fs[4].head(10)

            # skip plots if computing only CSV output files
            if args.csv:
                continue

            # plot output_png_fs_nfragpermol
            if Path(output_png_fs_nfragpermol).exists():
                logger.info("FS -- OUTPUT_PNG_FS_NFRAGPERMOL".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FS -- OUTPUT_PNG_FS_NFRAGPERMOL".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fs_nfragpermol,
                             output_png_fs_nfragpermol,
                             'NumFrags',
                             'Count',
                             f"Number of Fragment Hits Per Molecule in {args.dataset}",
                             x_label='Number of Fragment Hits Per Molecule',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Mols',
                             )

            # plot output_png_fs_fragmolcov
            if Path(output_png_fs_fragmolcov).exists():
                logger.info("FS -- OUTPUT_PNG_FS_FRAGMOLCOV".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FS -- OUTPUT_PNG_FS_FRAGMOLCOV".ljust(pad) + f"COMPUTING...")
                save_kdeplot(df=df_fs_fragmolcov,
                             output_png=output_png_fs_fragmolcov,
                             x_name='frag_ratio',
                             title=f"Distribution of Molecule Coverage by Fragments in {args.dataset}",
                             x_label='Molecule Coverage by Fragments',
                             y_label='Kernel Density Estimate of the Number of Molecules',
                             color=color,
                             )

            # plot output_png_fs_top10frags
            if Path(output_png_fs_top10frags).exists():
                logger.info("FS -- OUTPUT_PNG_FS_TOP10FRAGS".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FS -- OUTPUT_PNG_FS_TOP10FRAGS".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fs_top10frags,
                             output_png_fs_top10frags,
                             'idf',
                             'Count',
                             f"Top 10 Fragments in {args.dataset}",
                             x_label='Fragment ID',
                             y_label='Count',
                             color=color,
                             rotate_x=60,
                             perc_labels='Perc_FHits',
                             force_order=True,
                             fig_size = (12, 12),
                             )

            # plot output_png_fs_nfragpermol_u
            if Path(output_png_fs_nfragpermol_u).exists():
                logger.info("FS -- OUTPUT_PNG_FS_NFRAGPERMOL_U".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FS -- OUTPUT_PNG_FS_NFRAGPERMOL_U".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fs_nfragpermol_u,
                             output_png_fs_nfragpermol_u,
                             'NumFrags',
                             'Count',
                             f"Number of Unique Fragment Hits Per Molecule in {args.dataset}",
                             x_label='Number of Unique Fragment Hits Per Molecule',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Mols',
                             )

            # plot output_png_fs_top10frags_u
            if Path(output_png_fs_top10frags_u).exists():
                logger.info("FS -- OUTPUT_PNG_FS_TOP10FRAGS_U".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FS -- OUTPUT_PNG_FS_TOP10FRAGS_U".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fs_top10frags_u,
                             output_png_fs_top10frags_u,
                             'idf',
                             'Count',
                             f"Top 10 Unique Fragments in {args.dataset}",
                             x_label='Fragment ID',
                             y_label='Count',
                             color=color,
                             rotate_x=60,
                             perc_labels='Perc_FHits',
                             force_order=True,
                             fig_size = (12, 12),
                             )

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ FCC ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #

        elif step == 'fcc':

            _print_title("CHECKING RESULTS FROM FCC", logger, pad)
            # check if data is complete
            subdir = 'fcc'
            if not _step_is_valid(subdir, errors_data, errors_log):
                logger.error(f"MISSING DATA OR LOG FOR PREP... SKIPPING THIS STEP!")
                continue
            # define outputs
            output_csv_fcc_fcc = f"{output_folder}/{args.dataset}_fcc_fcc.csv"
            output_csv_fcc_fc = f"{output_folder}/{args.dataset}_fcc_fc.csv"
            output_csv_fcc_ffo = f"{output_folder}/{args.dataset}_fcc_ffo.csv"
            output_png_fcc_fcc = output_csv_fcc_fcc.replace('.csv', '.png')
            output_png_fcc_fc = output_csv_fcc_fc.replace('.csv', '.png')
            output_png_fcc_ffo = output_csv_fcc_ffo.replace('.csv', '.png')
            logger.info("FCC -- OUTPUT_CSV_FCC_FCC".ljust(pad) + f"{output_csv_fcc_fcc}")
            logger.info("FCC -- OUTPUT_CSV_FCC_FC".ljust(pad) + f"{output_csv_fcc_fc}")
            logger.info("FCC -- OUTPUT_CSV_FCC_FFO".ljust(pad) + f"{output_csv_fcc_ffo}")
            logger.info("FCC -- OUTPUT PNG FILES HAVE THE SAME FILE NAMES AS OUTPUT CSV FILES")
            # retrieve data
            output_csv_files = [output_csv_fcc_fcc, output_csv_fcc_fc, output_csv_fcc_ffo]
            output_png_files = [output_png_fcc_fcc, output_png_fcc_fc, output_png_fcc_ffo]
            if all([Path(x).exists() for x in output_png_files]):
                logger.info(f"FCC -- ALL OUTPUT PNG FILES ARE ALREADY AVAILABLE, NOTHING TO DO!")
            elif all([Path(x).exists() for x in output_csv_files]):
                logger.info(f"FCC -- PARSING OUTPUT CSV FILES INSTEAD OF COMPUTING THEM")
                df_fcc_fcc = load.file(output_csv_fcc_fcc)
                df_fcc_fc = load.file(output_csv_fcc_fc, decode=False).head(10)  # pot top 10 and no need for mol objects at the moment
                df_fcc_ffo = load.file(output_csv_fcc_ffo)
            else:
                logger.info(f"FCC -- COMPUTING OUTPUT CSV FILES")
                if args.frags is None or not Path(args.frags).exists():
                    logging.error("FCC -- ERROR! COULD NOT FOUND FRAGMENT FILE! SKIPPING STEP...")
                    continue
                df_fcc_fcc, df_fcc_fc, df_fcc_ffo = get_dfs_fcc(p, args.frags, compute_df_ffo=True)
                save.file(df_fcc_fcc, output_csv_fcc_fcc)
                save.file(df_fcc_fc, output_csv_fcc_fc)
                save.file(df_fcc_ffo, output_csv_fcc_ffo)
                df_fcc_fc = df_fcc_fc.head(10)  # use only top 10 for plotting

            # skip plots if computing only CSV output files
            if args.csv:
                continue

            # plot output_png_fcc_fcc
            if Path(output_png_fcc_fcc).exists():
                logger.info("FCC -- OUTPUT_PNG_FCC_COUNTS".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FCC -- OUTPUT_PNG_FCC_COUNTS".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fcc_fcc,
                             output_png_fcc_fcc,
                             'Category',
                             'Count',
                             f"Fragment Combination Classification in {args.dataset}",
                             x_label='Fragment Combination Categories',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc',
                             )

            # plot output_png_fcc_fc
            if Path(output_png_fcc_fc).exists():
                logger.info("FCC -- OUTPUT_PNG_FC_COUNTS".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FCC -- OUTPUT_PNG_FC_COUNTS".ljust(pad) + f"COMPUTING...")  # for that one percentage labelling go crazy...!
                save_barplot(df_fcc_fc,
                             output_png_fcc_fc,
                             'fc',
                             'Count',
                             f"Top 20 Fragment Combinations in {args.dataset}",
                             x_label='Fragment Combinations',
                             y_label='Count',
                             color=color,
                             rotate_x=60,
                             perc_labels='Perc',
                             )

            # plot output_png_fcc_ffo
            if Path(output_png_fcc_ffo).exists():
                logger.info("FCC -- OUTPUT_PNG_FCC_FFO".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FCC -- OUTPUT_PNG_FCC_FFO".ljust(pad) + f"COMPUTING...")  # for that one percentage labelling go crazy...!
                save_barplot(df_fcc_ffo,
                             output_png_fcc_ffo,
                             'NumOverlaps',
                             'Count',
                             f"Number of Overlap per Molecule in {args.dataset}",
                             x_label='Number of Overlaps',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Mols',
                             )


# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ACT ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #





# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ FMAP ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #


        elif step == 'fmap':
            subdir = 'fmap'
            _print_title(f"CHECKING RESULTS FROM {subdir.upper()}", logger, pad)
            # check if data is complete
            if not _step_is_valid(subdir, errors_data, errors_log):
                logger.error(f"MISSING DATA OR LOG FOR {subdir.upper()}... SKIPPING THIS STEP!")
                continue
            # define outputs
            output_csv_fm_nfragpermol = f"{output_folder}/{args.dataset}_fm_nfragpermol.csv"
            output_csv_fm_nfragpermol_u = f"{output_folder}/{args.dataset}_fm_nfragpermol_u.csv"
            output_csv_fm_fragmolcov = f"{output_folder}/{args.dataset}_fm_fragmolcov.csv"
            output_csv_fm_top10frags = f"{output_folder}/{args.dataset}_fm_top10frags.csv"
            output_csv_fm_top10frags_u = f"{output_folder}/{args.dataset}_fm_top10frags_u.csv"
            output_csv_fm_fcc = f"{output_folder}/{args.dataset}_fm_fcc.csv"
            output_csv_fm_top10fc = f"{output_folder}/{args.dataset}_fm_fc.csv"
            output_csv_fm_nfmpermol = f"{output_folder}/{args.dataset}_fm_nfmpermol.csv"
            output_png_fm_nfragpermol = output_csv_fm_nfragpermol.replace('.csv', '.png')
            output_png_fm_nfragpermol_u = output_csv_fm_nfragpermol_u.replace('.csv', '.png')
            output_png_fm_fragmolcov = output_csv_fm_fragmolcov.replace('.csv', '.png')
            output_png_fm_top10frags = output_csv_fm_top10frags.replace('.csv', '.png')
            output_png_fm_top10frags_u = output_csv_fm_top10frags_u.replace('.csv', '.png')
            output_png_fm_fcc = output_csv_fm_fcc.replace('.csv', '.png')
            output_png_fm_top10fc = output_csv_fm_top10fc.replace('.csv', '.png')
            output_png_fm_nfmpermol = output_csv_fm_nfmpermol.replace('.csv', '.png')
            logger.info("FM -- OUTPUT_CSV_FM_NFRAGPERMOL".ljust(pad) + f"{output_csv_fm_nfragpermol}")
            logger.info("FM -- OUTPUT_CSV_FM_NFRAGPERMOL_U".ljust(pad) + f"{output_csv_fm_nfragpermol_u}")
            logger.info("FM -- OUTPUT_CSV_FM_FRAGMOLCOV".ljust(pad) + f"{output_csv_fm_fragmolcov}")
            logger.info("FM -- OUTPUT_CSV_FM_TOP10FRAGS".ljust(pad) + f"{output_csv_fm_top10frags}")
            logger.info("FM -- OUTPUT_CSV_FM_TOP10FRAGS_U".ljust(pad) + f"{output_csv_fm_top10frags_u}")
            logger.info("FM -- OUTPUT_CSV_FM_FCC".ljust(pad) + f"{output_csv_fm_fcc}")
            logger.info("FM -- OUTPUT_CSV_FM_TOP10FC".ljust(pad) + f"{output_csv_fm_top10fc}")
            logger.info("FM -- OUTPUT_CSV_FM_NFRAGPERMOL".ljust(pad) + f"{output_csv_fm_nfmpermol}")
            logger.info("FM -- OUTPUT PNG FILES HAVE THE SAME FILE NAMES AS OUTPUT CSV FILES")
            # retrieve data
            output_csv_files = [output_csv_fm_nfragpermol, output_csv_fm_nfragpermol_u,
                                output_csv_fm_fragmolcov,
                                output_csv_fm_top10frags, output_csv_fm_top10frags_u,
                                output_csv_fm_fcc, output_csv_fm_top10fc,
                                output_csv_fm_nfmpermol,
                                ]
            output_png_files = [output_png_fm_nfragpermol, output_png_fm_nfragpermol_u,
                                output_png_fm_fragmolcov,
                                output_png_fm_top10frags, output_png_fm_top10frags_u,
                                output_png_fm_fcc, output_png_fm_top10fc,
                                output_png_fm_nfmpermol,
                                ]
            if all([Path(x).exists() for x in output_png_files]):
                logger.info(f"FM -- ALL OUTPUT PNG FILES ARE ALREADY AVAILABLE, NOTHING TO DO!")
            elif all([Path(x).exists() for x in output_csv_files]):
                logger.info(f"FM -- PARSING OUTPUT CSV FILES INSTEAD OF COMPUTING THEM")
                df_fm_nfragpermol = load.file(output_csv_fm_nfragpermol)
                df_fm_nfragpermol_u = load.file(output_csv_fm_nfragpermol_u)
                df_fm_fragmolcov = load.file(output_csv_fm_fragmolcov)
                df_fm_top10frags = load.file(output_csv_fm_top10frags).head(10)
                df_fm_top10frags_u = load.file(output_csv_fm_top10frags_u).head(10)
                df_fm_fcc = load.file(output_csv_fm_fcc)
                df_fm_top10fc = load.file(output_csv_fm_top10fc).head(10)
                df_fm_nfmpermol = load.file(output_csv_fm_nfmpermol)
            else:
                logger.info(f"FM -- COMPUTING OUTPUT CSV FILES")
                if args.frags is None or not Path(args.frags).exists():
                    logging.error("FM -- ERROR! COULD NOT FOUND FRAGMENT FILE! SKIPPING STEP...")
                    continue
                dfs_fm =  get_df_fm(p, args.frags)
                for df_fm, output_csv_fm in zip(dfs_fm, output_csv_files):
                    save.file(df_fm, output_csv_fm)
                df_fm_nfragpermol = dfs_fm[0]
                df_fm_nfragpermol_u =  dfs_fm[1]
                df_fm_fragmolcov =  dfs_fm[2]
                df_fm_top10frags =  dfs_fm[3].head(10)
                df_fm_top10frags_u =  dfs_fm[4].head(10)
                df_fm_fcc = dfs_fm[5]
                df_fm_top10fc = dfs_fm[6].head(10)
                df_fm_nfmpermol = dfs_fm[7]
            # skip plots if computing only CSV output files
            if args.csv:
                continue

            # plot output_png_fm_nfragpermol
            if Path(output_png_fm_nfragpermol).exists():
                logger.info("FM -- OUTPUT_PNG_FM_NFRAGPERMOL".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FM -- OUTPUT_PNG_FM_NFRAGPERMOL".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fm_nfragpermol,
                             output_png_fm_nfragpermol,
                             'NumFrags',
                             'Count',
                             f"Number of Fragment Hits Per Molecule in {args.dataset}",
                             x_label='Number of Fragment Hits Per Molecule',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Mols',
                             )

            # plot output_png_fm_fragmolcov
            if Path(output_png_fm_fragmolcov).exists():
                logger.info("FM -- OUTPUT_PNG_FM_FRAGMOLCOV".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FM -- OUTPUT_PNG_FM_FRAGMOLCOV".ljust(pad) + f"COMPUTING...")
                save_kdeplot(df=df_fm_fragmolcov,
                             output_png=output_png_fm_fragmolcov,
                             x_name='frag_ratio',
                             title=f"Distribution of Molecule Coverage by Fragments in {args.dataset}",
                             x_label='Molecule Coverage by Fragments',
                             y_label='Kernel Density Estimate of the Number of Molecules',
                             color=color,
                             )

            # plot output_png_fm_top10frags
            if Path(output_png_fm_top10frags).exists():
                logger.info("FM -- OUTPUT_PNG_FM_TOP10FRAGS".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FM -- OUTPUT_PNG_FM_TOP10FRAGS".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fm_top10frags,
                             output_png_fm_top10frags,
                             'idf',
                             'Count',
                             f"Top 10 Fragments in {args.dataset}",
                             x_label='Fragment ID',
                             y_label='Count',
                             color=color,
                             rotate_x=60,
                             perc_labels='Perc_FHits',
                             force_order=True,
                             fig_size = (12, 12),
                             )

            # plot output_png_fm_nfragpermol_u
            if Path(output_png_fm_nfragpermol_u).exists():
                logger.info("FM -- OUTPUT_PNG_FM_NFRAGPERMOL_U".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FM -- OUTPUT_PNG_FM_NFRAGPERMOL_U".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fm_nfragpermol_u,
                             output_png_fm_nfragpermol_u,
                             'NumFrags',
                             'Count',
                             f"Number of Unique Fragment Hits Per Molecule in {args.dataset}",
                             x_label='Number of Unique Fragment Hits Per Molecule',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Mols',
                             )

            # plot output_png_fm_top10frags_u
            if Path(output_png_fm_top10frags_u).exists():
                logger.info("FM -- OUTPUT_PNG_FM_TOP10FRAGS_U".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FM -- OUTPUT_PNG_FM_TOP10FRAGS_U".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fm_top10frags_u,
                             output_png_fm_top10frags_u,
                             'idf',
                             'Count',
                             f"Top 10 Unique Fragments in {args.dataset}",
                             x_label='Fragment ID',
                             y_label='Count',
                             color=color,
                             rotate_x=60,
                             perc_labels='Perc_FHits',
                             force_order=True,
                             fig_size = (12, 12),
                             )

            # plot output_png_fm_fcc
            if Path(output_png_fm_fcc).exists():
                logger.info("FM -- OUTPUT_PNG_FM_COUNTS".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FM -- OUTPUT_PNG_FM_COUNTS".ljust(pad) + f"COMPUTING...")
                save_barplot(df_fm_fcc,
                             output_png_fm_fcc,
                             'Category',
                             'Count',
                             f"Fragment Combination Classification in {args.dataset}",
                             x_label='Fragment Combination Categories',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc',
                             )

            # plot output_png_fm_top10fc
            if Path(output_png_fm_top10fc).exists():
                logger.info("FM -- OUTPUT_PNG_FC_COUNTS".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("FM -- OUTPUT_PNG_FC_COUNTS".ljust(pad) + f"COMPUTING...")  # for that one percentage labelling go crazy...!
                save_barplot(df_fm_top10fc,
                             output_png_fm_top10fc,
                             'fc',
                             'Count',
                             f"Top 20 Fragment Combinations in {args.dataset}",
                             x_label='Fragment Combinations',
                             y_label='Count',
                             color=color,
                             rotate_x=60,
                             perc_labels='Perc',
                             )

            # plot output_png_fm_nfmpermol
            if Path(output_png_fm_nfmpermol).exists():
                logger.info("FM -- OUTPUT_PNG_FM_NFMPERMOL".ljust(pad) + "ALREADY DONE")
            else:
                logger.info("FM -- OUTPUT_PNG_FM_NFMPERMOL".ljust(pad) + "COMPUTING...")
                save_barplot(df_fm_nfmpermol,
                             output_png_fm_nfmpermol,
                             'NumFM',
                             'Count',
                             f"Number of Fragment Maps Per Molecule in {args.dataset}",
                             x_label='Number of Fragment Maps Per Molecule',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Mols',
                             )

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ PNP ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #


        elif step == 'pnp':
            _print_title("CHECKING RESULTS FROM PNP", logger, pad)
            # check if data is complete
            subdir = 'pnp'
            if not _step_is_valid(subdir, errors_data, errors_log):
                logger.error(f"MISSING DATA OR LOG FOR PNP... SKIPPING THIS STEP!")
                continue
            # define outputs
            output_csv_pnp_pnp = f"{output_folder}/{args.dataset}_pnp_pnp.csv"
            output_png_pnp_pnp = output_csv_pnp_pnp.replace('.csv', '.png')
            logger.info("PNP -- OUTPUT_CSV_PNP_PNP".ljust(pad) + f"{output_csv_pnp_pnp}")
            logger.info("PNP -- OUTPUT PNG FILES HAVE THE SAME FILE NAMES AS OUTPUT CSV FILES")
            # retrieve data
            output_csv_files = [output_csv_pnp_pnp]
            output_png_files = [output_png_pnp_pnp]
            if all([Path(x).exists() for x in output_png_files]):
                logger.info(f"PNP -- ALL OUTPUT PNG FILES ARE ALREADY AVAILABLE, NOTHING TO DO!")
            elif all([Path(x).exists() for x in output_csv_files]):
                logger.info(f"PNP -- PARSING OUTPUT CSV FILES INSTEAD OF COMPUTING THEM")
                df_pnp_pnp = load.file(output_csv_pnp_pnp)
            else:
                logger.info(f"PNP -- COMPUTING OUTPUT CSV FILES")
                df_pnp_pnp =  get_df_pnp(p)
                save.file(df_pnp_pnp, output_csv_pnp_pnp)

            # skip plots if computing only CSV output files
            if args.csv:
                continue
            # plot output_png_pnp_pnp
            if Path(output_png_pnp_pnp).exists():
                logger.info("PNP -- OUTPUT_PNG_PNP_PNP".ljust(pad) + f"ALREADY DONE")
            else:
                logger.info("PNP -- OUTPUT_PNG_PNP_PNP".ljust(pad) + f"COMPUTING...")
                save_barplot(df_pnp_pnp,
                             output_png_pnp_pnp,
                             'Category',
                             'Count',
                             f"Identification of PNPs in {args.dataset}",
                             x_label='Category',
                             y_label='Count',
                             color=color,
                             perc_labels='Perc_Mols',
                             fig_size = (12, 12),
                             )
