"""Evaluate model structure against reference.

eg.

  ost compare-structures \\
      --model <MODEL> \\
      --reference <REF> \\
      --output output.json \\
      --lddt \\
      --structural-checks \\
      --consistency-checks \\
      --molck \\
      --remove oxt hyd \\
      --map-nonstandard-residues

Here we describe how the parameters can be set to mimick a CAMEO evaluation
(as of August 2018).

CAMEO calls the lddt binary as follows:

  lddt \\
      -p <PARAMETER FILE> \\
      -f \\
      -a 15 \\
      -b 15 \\
      -r 15 \\
      <MODEL> \\
      <REF>

Only model structures are "Molck-ed" in CAMEO. The call to molck is as follows:

  molck \\
      --complib=<COMPOUND LIB> \\
      --rm=hyd,oxt,unk,nonstd \\
      --fix-ele \\
      --map-nonstd \\
      --out=<OUTPUT> \\
      <FILEPATH>

To be as much compatible with with CAMEO as possible one should call
compare-structures as follows:

  ost compare-structures \\
      --model <MODEL> \\
      --reference <REF> \\
      --output output.json \\
      --molck \\
      --remove oxt hyd unk nonstd \\
      --clean-element-column \\
      --map-nonstandard-residues \\
      --structural-checks \\
      --bond-tolerance 15.0 \\
      --angle-tolerance 15.0 \\
      --residue-number-alignment \\
      --consistency-checks \\
      --qs-score \\
      --lddt \\
      --inclusion-radius 15.0
"""

import os
import sys
import json
import argparse

import ost
from ost.io import (LoadPDB, LoadMMCIF, SavePDB, MMCifInfoBioUnit, MMCifInfo,
                    MMCifInfoTransOp, ReadStereoChemicalPropsFile, profiles)
from ost import PushVerbosityLevel
from ost.mol.alg import (qsscoring, Molck, MolckSettings, lDDTSettings,
                         CheckStructure, ResidueNamesMatch)
from ost.conop import (CompoundLib, SetDefaultLib, GetDefaultLib,
                       RuleBasedProcessor)
from ost.seq.alg.renumber import Renumber


def _GetDefaultShareFilePath(filename):
    """Look for filename in working directory and OST shared data path.
    :return: Path to valid file or None if not found.
    """
    # Try current directory
    cwd = os.path.abspath(os.getcwd())
    file_path = os.path.join(cwd, filename)
    if not os.path.isfile(file_path):
        try:
            file_path = os.path.join(ost.GetSharedDataPath(), filename)
        except RuntimeError:
            # Ignore errors here (caught later together with non-existing file)
            pass
        if not os.path.isfile(file_path):
            file_path = None
    # Either file_path is valid file path or None
    return file_path

def _GetDefaultParameterFilePath():
    # Try to get in default locations
    parameter_file_path = _GetDefaultShareFilePath("stereo_chemical_props.txt")
    if parameter_file_path is None:
        msg = (
            "Could not set default stereochemical parameter file. In "
            "order to use the default one please set $OST_ROOT "
            "environmental variable, run the script with OST binary or"
            " provide a local copy of 'stereo_chemical_props.txt' in "
            "CWD. Alternatively provide the path to the local copy.")
    else:
        msg = ""
    return parameter_file_path, msg

def _GetDefaultCompoundLibraryPath():
    # Try to get in default locations
    compound_library_path = _GetDefaultShareFilePath("compounds.chemlib")
    if compound_library_path is None:
        msg = (
            "Could not set default compounds library path. In "
            "order to use the default one please set $OST_ROOT "
            "environmental variable, run the script with OST binary or"
            " provide a local copy of 'compounds.chemlib' in CWD"
            ". Alternatively provide the path to the local copy.")
    else:
        msg = ""
    return compound_library_path, msg

def _ParseArgs():
    """Parse command-line arguments."""

    parser = argparse.ArgumentParser(
        formatter_class=argparse.RawTextHelpFormatter,
        description=__doc__,
        prog="ost compare-structures")

    #
    # Required arguments
    #

    group_required = parser.add_argument_group('required arguments')

    group_required.add_argument(
        "-m",
        "--model",
        dest="model",
        required=True,
        help=("Path to the model file."))
    group_required.add_argument(
        "-r",
        "--reference",
        dest="reference",
        required=True,
        help=("Path to the reference file."))

    #
    # General arguments
    #

    group_general = parser.add_argument_group('general arguments')

    group_general.add_argument(
        '-v',
        '--verbosity',
        type=int,
        default=3,
        help="Set verbosity level. Defaults to 3.")
    group_general.add_argument(
        "-o",
        "--output",
        dest="output",
        help=("Output file name. The output will be saved as a JSON file."))
    group_general.add_argument(
        "-d",
        "--dump-structures",
        dest="dump_structures",
        default=False,
        action="store_true",
        help=("Dump cleaned structures used to calculate all the scores as\n"
              "PDB files using specified suffix. Files will be dumped to the\n"
              "same location as original files."))
    group_general.add_argument(
        "-ds",
        "--dump-suffix",
        dest="dump_suffix",
        default=".compare.structures.pdb",
        help=("Use this suffix to dump structures.\n"
              "Defaults to .compare.structures.pdb."))
    group_general.add_argument(
        "-rs",
        "--reference-selection",
        dest="reference_selection",
        default="",
        help=("Selection performed on reference structures."))
    group_general.add_argument(
        "-ms",
        "--model-selection",
        dest="model_selection",
        default="",
        help=("Selection performed on model structures."))
    group_general.add_argument(
        "-ca",
        "--c-alpha-only",
        dest="c_alpha_only",
        default=False,
        action="store_true",
        help=("Use C-alpha atoms only. Equivalent of calling the action with\n"
              "'--model-selection=\"aname=CA\" "
              "--reference-selection=\"aname=CA\"'\noptions."))
    group_general.add_argument(
        "-ft",
        "--fault-tolerant",
        dest="fault_tolerant",
        default=False,
        action="store_true",
        help=("Fault tolerant parsing."))
    group_general.add_argument(
        "-cl",
        "--compound-library",
        dest="compound_library",
        default=None,
        help=("Location of the compound library file (compounds.chemlib).\n"
              "If not provided, the following locations are searched in this\n"
              "order: 1. Working directory, 2. OpenStructure standard library"
              "\nlocation."))

    #
    # Molecular check arguments
    #

    group_molck = parser.add_argument_group('molecular check arguments')

    group_molck.add_argument(
        "-ml",
        "--molck",
        dest="molck",
        default=False,
        action="store_true",
        help=("Run molecular checker to clean up input."))
    group_molck.add_argument(
        "-rm",
        "--remove",
        dest="remove",
        nargs="+",  # *, +, ?, N
        required=False,
        default=["hyd"],
        help=("Remove atoms and residues matching some criteria:\n"
              " * zeroocc - Remove atoms with zero occupancy\n"
              " * hyd - remove hydrogen atoms\n"
              " * oxt - remove terminal oxygens\n"
              " * nonstd - remove all residues not one of the 20\n"
              "            standard amino acids\n"
              " * unk - Remove unknown and atoms not following the\n"
              "         nomenclature\n"
              "Defaults to hyd."))
    group_molck.add_argument(
        "-ce",
        "--clean-element-column",
        dest="clean_element_column",
        default=False,
        action="store_true",
        help=("Clean up element column"))
    group_molck.add_argument(
        "-mn",
        "--map-nonstandard-residues",
        dest="map_nonstandard_residues",
        default=False,
        action="store_true",
        help=("Map modified residues back to the parent amino acid, for\n"
              "example MSE -> MET, SEP -> SER."))
    
    #
    # Structural check arguments
    #

    group_sc = parser.add_argument_group('structural check arguments')

    group_sc.add_argument(
        "-sc",
        "--structural-checks",
        dest="structural_checks",
        default=False,
        action="store_true",
        help=("Perform structural checks and filter input data."))
    group_sc.add_argument(
        "-p",
        "--parameter-file",
        dest="parameter_file",
        default=None,
        help=("Location of the stereochemical parameter file\n"
              "(stereo_chemical_props.txt).\n"
              "If not provided, the following locations are searched in this\n"
              "order: 1. Working directory, 2. OpenStructure standard library"
              "\nlocation."))
    group_sc.add_argument(
        "-bt",
        "--bond-tolerance",
        dest="bond_tolerance",
        type=float,
        default=12.0,
        help=("Tolerance in STD for bonds. Defaults to 12."))
    group_sc.add_argument(
        "-at",
        "--angle-tolerance",
        dest="angle_tolerance",
        type=float,
        default=12.0,
        help=("Tolerance in STD for angles. Defaults to 12."))

    #
    # Chain mapping arguments
    #

    group_cm = parser.add_argument_group('chain mapping arguments')

    group_cm.add_argument(
        "-c",
        "--chain-mapping",
        nargs="+",
        type=lambda x: x.split(":"),
        dest="chain_mapping",
        help=("Mapping of chains between the reference and the model.\n"
              "Each separate mapping consist of key:value pairs where key\n"
              "is the chain name in reference and value is the chain name in\n"
              "model."))
    group_cm.add_argument(
        "--qs-max-mappings-extensive",
        dest="qs_max_mappings_extensive",
        type=int,
        default=1000000,
        help=("Maximal number of chain mappings to test for 'extensive'\n"
              "chain mapping scheme which is used as a last resort if\n"
              "other schemes failed. The extensive chain mapping search\n"
              "must in the worst case check O(N!) possible mappings for\n"
              "complexes with N chains. Two octamers without symmetry\n"
              "would require 322560 mappings to be checked. To limit\n"
              "computations, no scores are computed if we try more than\n"
              "the maximal number of chain mappings. Defaults to 1000000."))

    #
    # Sequence alignment arguments
    #

    group_aln = parser.add_argument_group('sequence alignment arguments')

    group_aln.add_argument(
        "-cc",
        "--consistency-checks",
        dest="consistency_checks",
        default=False,
        action="store_true",
        help=("Take consistency checks into account. By default residue name\n"
              "consistency between a model-reference pair would be checked\n"
              "but only a warning message will be displayed and the script\n"
              "will continue to calculate scores. If this flag is ON, checks\n"
              "will not be ignored and if the pair does not pass the test\n"
              "all the scores for that pair will be marked as a FAILURE."))
    group_aln.add_argument(
        "-rna",
        "--residue-number-alignment",
        dest="residue_number_alignment",
        default=False,
        action="store_true",
        help=("Make alignment based on residue number instead of using\n"
              "a global BLOSUM62-based alignment."))

    #
    # QS score arguments
    #

    group_qs = parser.add_argument_group('QS score arguments')

    group_qs.add_argument(
        "-qs",
        "--qs-score",
        dest="qs_score",
        default=False,
        action="store_true",
        help=("Calculate QS-score."))
    group_qs.add_argument(
        "--qs-rmsd",
        dest="qs_rmsd",
        default=False,
        action="store_true",
        help=("Calculate CA RMSD between shared CA atoms of mapped chains.\n"
              "This uses a superposition using all mapped chains which\n"
              "minimizes the CA RMSD."))

    #
    # lDDT score arguments
    #

    group_lddt = parser.add_argument_group('lDDT score arguments')

    group_lddt.add_argument(
        "-l",
        "--lddt",
        dest="lddt",
        default=False,
        action="store_true",
        help=("Calculate lDDT."))
    group_lddt.add_argument(
        "-ir",
        "--inclusion-radius",
        dest="inclusion_radius",
        type=float,
        default=15.0,
        help=("Distance inclusion radius for lDDT. Defaults to 15 A."))
    group_lddt.add_argument(
        "-ss",
        "--sequence-separation",
        dest="sequence_separation",
        type=int,
        default=0,
        help=("Sequence separation. Only distances between residues whose\n"
              "separation is higher than the provided parameter are\n"
              "considered when computing the score. Defaults to 0."))
    group_lddt.add_argument(
        "-spr",
        "--save-per-residue-scores",
        dest="save_per_residue_scores",
        default=False,
        action="store_true",
        help=(""))

    # Print full help is no arguments provided
    if len(sys.argv) == 1:
        parser.print_help(sys.stderr)
        sys.exit(1)

    opts = parser.parse_args()
    # Set chain mapping
    if opts.chain_mapping is not None:
        try:
            opts.chain_mapping = dict(opts.chain_mapping)
        except ValueError:
            parser.error(
                "Cannot parse chain mapping into dictionary. The "
                "correct format is: key:value [key2:value2 ...].")
    
    # Check parameter file if structural checks are on
    if opts.structural_checks:
        if opts.parameter_file is None:
            # try to get default if none provided
            opts.parameter_file, msg = _GetDefaultParameterFilePath()
            if msg:
                parser.error(msg)
        else:
            # if provided it must exist
            if not os.path.isfile(opts.parameter_file):
                parser.error("Parameter file %s does not exist." \
                             % opts.parameter_file)

    # Check compound library path (always required!)
    if opts.compound_library is None:
        # try to get default if none provided
        opts.compound_library, msg = _GetDefaultCompoundLibraryPath()
        if msg:
            parser.error(msg)
    else:
        # if provided it must exist
        if not os.path.isfile(opts.compound_library):
            parser.error("Compounds library file %s does not exist." \
                         % opts.compound_library)

    # Check model and reference paths
    if not os.path.isfile(opts.model):
        parser.error("Model file %s does not exist." % opts.model)
    if not os.path.isfile(opts.reference):
        parser.error("Reference file %s does not exist." % opts.reference)

    return opts


def _SetCompoundsChemlib(path_to_chemlib):
    """Set default compound library for OST."""
    # NOTE: This is adapted from ProMod3 code and should in the future be doable
    #       with some shared OST code!
    compound_lib = CompoundLib.Load(path_to_chemlib)
    SetDefaultLib(compound_lib)
    processor = RuleBasedProcessor(compound_lib)
    for profile_name in profiles:
        profiles[profile_name].processor = processor.Copy()


def _RevertChainNames(ent):
    """Revert chain names to original names.

    By default the first chain with given name will not have any number
    attached to it ie. if there are two chains mapping to chain A the resulting
    chain names will be: A and A2.
    """
    editor = ent.EditXCS()
    suffix = "_tmp"  # just a suffix for temporary chain name
    separator = ""  # dot causes selection error
    used_names = dict()
    reverted_chains = dict()
    for chain in ent.chains:
        try:
            original_name = chain.GetStringProp("original_name")
        except Exception as ex:
            ost.LogError("Cannot revert chain %s back to original: %s" % (
                chain.name,
                str(ex)))
            reverted_chains[chain.name] = chain.name
            editor.RenameChain(chain, chain.name + suffix)
            continue
        new_name = original_name
        if new_name not in used_names:
            used_names[original_name] = 2
            reverted_chains[chain.name] = new_name
            editor.RenameChain(chain, chain.name + suffix)
        else:
            new_name = "%s%s%i" % (original_name,
                                   separator,
                                   used_names[original_name])
            reverted_chains[chain.name] = new_name
            editor.RenameChain(chain, chain.name + suffix)
            used_names[original_name] += 1
    for chain in ent.chains:
        editor.RenameChain(chain, reverted_chains[chain.name[:-len(suffix)]])
    rev_out = ["%s -> %s" % (on, nn) for on, nn in list(reverted_chains.items())]
    ost.LogInfo("Reverted chains: %s" % ", ".join(rev_out))


def _CheckConsistency(alignments, log_error):
    is_cons = True
    for alignment in alignments:
        ref_chain = Renumber(alignment.GetSequence(0)).CreateFullView()
        mdl_chain = Renumber(alignment.GetSequence(1)).CreateFullView()
        new_is_cons = ResidueNamesMatch(mdl_chain, ref_chain, log_error)
        is_cons = is_cons and new_is_cons
    return is_cons


def _GetAlignmentsAsFasta(alignments):
    """Get the alignments as FASTA formated string.

    :param alignments: Alignments
    :type alignments: list of AlignmentHandle
    :returns: list of alignments in FASTA format
    :rtype: list of strings
    """
    strings = list()
    for alignment in alignments:
        aln_str = ">reference:%s\n%s\n>model:%s\n%s" % (
            alignment.GetSequence(0).name,
            alignment.GetSequence(0).GetString(),
            alignment.GetSequence(1).name,
            alignment.GetSequence(1).GetString())
        strings.append(aln_str)
    return strings


def _ReadStructureFile(path, c_alpha_only=False, fault_tolerant=False,
                       selection=""):
    """Safely read structure file into OST entities (split by biounit).

    The function can read both PDB and mmCIF files.
    
    :param path: Path to the file.
    :type path: :class:`str`
    :returns: list of entities
    :rtype: :class:`list` of :class:`~ost.mol.EntityHandle`
    """

    def _Select(entity):
        if selection:
            ost.LogInfo("Selecting %s" % selection)
            ent_view = entity.Select(selection)
            entity = mol.CreateEntityFromView(ent_view, False)
        return entity

    entities = list()
    if not os.path.isfile(path):
        raise IOError("%s is not a file" % path)

    # Determine file format from suffix.
    ext = path.split(".")
    if ext[-1] == "gz":
        ext = ext[:-1]
    if len(ext) <= 1:
        raise RuntimeError(f"Could not determine format of file {path}.")
    sformat = ext[-1].lower()

    if sformat in ["pdb"]:
        entity = LoadPDB(
            path,
            fault_tolerant=fault_tolerant,
            calpha_only=c_alpha_only)
        if not entity.IsValid() or len(entity.residues) == 0:
            raise IOError("Provided file does not contain valid entity.")
        entity.SetName(os.path.basename(path))
        entity = _Select(entity)
        entities.append(entity)
    elif sformat in ["cif", "mmcif"]:
        try:
            tmp_entity, cif_info = LoadMMCIF(
                path,
                info=True,
                fault_tolerant=fault_tolerant,
                calpha_only=c_alpha_only)
            if len(cif_info.biounits) == 0:
                tbu = MMCifInfoBioUnit()
                tbu.id = 'ASU'
                tbu.details = 'asymmetric unit'
                for chain in tmp_entity.chains:
                    tbu.AddChain(str(chain))
                tinfo = MMCifInfo()
                tops = MMCifInfoTransOp()
                tinfo.AddOperation(tops)
                tbu.AddOperations(tinfo.GetOperations())
                entity = tbu.PDBize(tmp_entity, min_polymer_size=0)
                entity.SetName(os.path.basename(path) + ".au")
                _RevertChainNames(entity)
                entity = _Select(entity)
                entities.append(entity)
            elif len(cif_info.biounits) > 1:
                for i, biounit in enumerate(cif_info.biounits, 1):
                    entity = biounit.PDBize(tmp_entity, min_polymer_size=0)
                    if not entity.IsValid():
                        raise IOError(
                            "Provided file does not contain valid entity.")
                    entity.SetName(os.path.basename(path) + "." + str(i))
                    _RevertChainNames(entity)
                    entity = _Select(entity)
                    entities.append(entity)
            else:
                biounit = cif_info.biounits[0]
                entity = biounit.PDBize(tmp_entity, min_polymer_size=0)
                if not entity.IsValid():
                    raise IOError(
                        "Provided file does not contain valid entity.")
                entity.SetName(os.path.basename(path))
                _RevertChainNames(entity)
                entity = _Select(entity)
                entities.append(entity)

        except Exception:
            raise
    else:
        raise RuntimeError(f"Unsupported file extension found for file {path}.")

    return entities


def _MolckEntity(entity, options):
    """Molck the entity."""
    lib = GetDefaultLib()
    to_remove = tuple(options.remove)

    ms = MolckSettings(rm_unk_atoms="unk" in to_remove,
                       rm_non_std="nonstd" in to_remove,
                       rm_hyd_atoms="hyd" in to_remove,
                       rm_oxt_atoms="oxt" in to_remove,
                       rm_zero_occ_atoms="zeroocc" in to_remove,
                       colored=False,
                       map_nonstd_res=options.map_nonstandard_residues,
                       assign_elem=options.clean_element_column)
    Molck(entity, lib, ms)


def _Main():
    """Do the magic."""
    #
    # Setup
    opts = _ParseArgs()
    PushVerbosityLevel(opts.verbosity)
    _SetCompoundsChemlib(opts.compound_library)
    #
    # Read the input files
    ost.LogInfo("#" * 80)
    ost.LogInfo("Reading input files (fault_tolerant=%s)" %
                str(opts.fault_tolerant))
    ost.LogInfo(" --> reading model from %s" % opts.model)
    models = _ReadStructureFile(
        opts.model,
        c_alpha_only=opts.c_alpha_only,
        fault_tolerant=opts.fault_tolerant,
        selection=opts.model_selection)
    ost.LogInfo(" --> reading reference from %s" % opts.reference)
    references = _ReadStructureFile(
        opts.reference,
        c_alpha_only=opts.c_alpha_only,
        fault_tolerant=opts.fault_tolerant,
        selection=opts.reference_selection)
    # molcking
    if opts.molck:
        ost.LogInfo("#" * 80)
        ost.LogInfo("Cleaning up input with Molck")
        for reference in references:
            _MolckEntity(reference, opts)
        for model in models:
            _MolckEntity(model, opts)
    # restrict to peptides (needed for CheckStructure anyways)
    for i in range(len(references)):
        references[i] = references[i].Select("peptide=true")
    for i in range(len(models)):
        models[i] = models[i].Select("peptide=true")
    # structure checking
    if opts.structural_checks:
        ost.LogInfo("#" * 80)
        ost.LogInfo("Performing structural checks")
        stereochemical_parameters = ReadStereoChemicalPropsFile(
            opts.parameter_file)
        ost.LogInfo(" --> for reference(s)")
        for reference in references:
            ost.LogInfo("Checking %s" % reference.GetName())
            CheckStructure(reference,
                           stereochemical_parameters.bond_table,
                           stereochemical_parameters.angle_table,
                           stereochemical_parameters.nonbonded_table,
                           opts.bond_tolerance,
                           opts.angle_tolerance)
        ost.LogInfo(" --> for model(s)")
        for model in models:
            ost.LogInfo("Checking %s" % model.GetName())
            CheckStructure(model,
                           stereochemical_parameters.bond_table,
                           stereochemical_parameters.angle_table,
                           stereochemical_parameters.nonbonded_table,
                           opts.bond_tolerance,
                           opts.angle_tolerance)
    if len(models) > 1 or len(references) > 1:
        ost.LogInfo("#" * 80)
        ost.LogInfo(
            "Multiple complexes mode ON. All combinations will be tried.")

    result = {
        "result": {},
        "options": vars(opts)}
    result["options"]["cwd"] = os.path.abspath(os.getcwd())
    #
    # Perform scoring
    skipped = list()
    for model in models:
        model_name = model.GetName()
        model_results = dict()
        for reference in references:
            reference_name = reference.GetName()
            reference_results = {
                "info": dict()}
            ost.LogInfo("#" * 80)
            ost.LogInfo("Comparing %s to %s" % (
                model_name,
                reference_name))
            qs_scorer = qsscoring.QSscorer(reference,
                                           model,
                                           opts.residue_number_alignment)
            qs_scorer.max_mappings_extensive = opts.qs_max_mappings_extensive
            if opts.chain_mapping is not None:
                ost.LogInfo(
                    "Using custom chain mapping: %s" % str(
                        opts.chain_mapping))
                qs_scorer.chain_mapping = opts.chain_mapping
            else:
                try:
                  qs_scorer.chain_mapping  # just to initialize it
                except qsscoring.QSscoreError as ex:
                  ost.LogError('Chain mapping failed:', str(ex))
                  ost.LogError('Skipping comparison')
                  continue
            ost.LogInfo("-" * 80)
            ost.LogInfo("Checking consistency between %s and %s" % (
                        model_name, reference_name))
            is_cons = _CheckConsistency(
                qs_scorer.alignments,
                opts.consistency_checks)
            reference_results["info"]["residue_names_consistent"] = is_cons
            reference_results["info"]["mapping"] = {
                "chain_mapping": qs_scorer.chain_mapping,
                "chain_mapping_scheme": qs_scorer.chain_mapping_scheme,
                "alignments": _GetAlignmentsAsFasta(qs_scorer.alignments)}
            skip_score = False
            if opts.consistency_checks:
                if not is_cons:
                    msg = (("Residue names in model %s and in reference "
                            "%s are inconsistent.") % (
                                model_name,
                                reference_name))
                    ost.LogError(msg)
                    skip_score = True
                    skipped.append(skip_score)
                else:
                    ost.LogInfo("Consistency check: OK")
                    skipped.append(False)
            else:
                skipped.append(False)
                if not is_cons:
                    msg = (("Residue names in model %s and in reference "
                            "%s are inconsistent.\nThis might lead to "
                            "corrupted results.") % (
                                model_name,
                                reference_name))
                    ost.LogWarning(msg)
                else:
                    ost.LogInfo("Consistency check: OK")
            if opts.qs_rmsd:
                ost.LogInfo("-" * 80)
                if skip_score:
                    ost.LogInfo(
                        "Skipping QS-RMSD because consistency check failed")
                    reference_results["qs_rmsd"] = {
                        "status": "FAILURE",
                        "error": "Consistency check failed."}
                else:
                    ost.LogInfo("Computing QS-RMSD")
                    try:
                        reference_results["qs_rmsd"] = {
                            "status": "SUCCESS",
                            "error": "",
                            "ca_rmsd": qs_scorer.superposition.rmsd}
                    except qsscoring.QSscoreError as ex:
                        ost.LogError('QS-RMSD failed:', str(ex))
                        reference_results["qs_rmsd"] = {
                            "status": "FAILURE",
                            "error": str(ex)}
            if opts.qs_score:
                ost.LogInfo("-" * 80)
                if skip_score:
                    ost.LogInfo(
                        "Skipping QS-score because consistency check failed")
                    reference_results["qs_score"] = {
                        "status": "FAILURE",
                        "error": "Consistency check failed.",
                        "global_score": 0.0,
                        "best_score": 0.0}
                else:
                    ost.LogInfo("Computing QS-score")
                    try:
                        reference_results["qs_score"] = {
                            "status": "SUCCESS",
                            "error": "",
                            "global_score": qs_scorer.global_score,
                            "best_score": qs_scorer.best_score}
                    except qsscoring.QSscoreError as ex:
                        # default handling: report failure and set score to 0
                        ost.LogError('QSscore failed:', str(ex))
                        reference_results["qs_score"] = {
                            "status": "FAILURE",
                            "error": str(ex),
                            "global_score": 0.0,
                            "best_score": 0.0}
            # Calculate lDDT
            if opts.lddt:
                ost.LogInfo("-" * 80)
                ost.LogInfo("Computing lDDT scores")
                lddt_results = {
                    "single_chain_lddt": list()
                }
                lddt_settings = lDDTSettings(
                    radius=opts.inclusion_radius,
                    sequence_separation=opts.sequence_separation,
                    label="lddt")
                ost.LogInfo("lDDT settings: ")
                ost.LogInfo(str(lddt_settings).rstrip())
                ost.LogInfo("===")
                oligo_lddt_scorer = qs_scorer.GetOligoLDDTScorer(lddt_settings)
                for mapped_lddt_scorer in oligo_lddt_scorer.mapped_lddt_scorers:
                    # Get data
                    lddt_scorer = mapped_lddt_scorer.lddt_scorer
                    model_chain = mapped_lddt_scorer.model_chain_name
                    reference_chain = mapped_lddt_scorer.reference_chain_name
                    if skip_score:
                        ost.LogInfo(
                            " --> Skipping single chain lDDT because "
                            "consistency check failed")
                        lddt_results["single_chain_lddt"].append({
                            "status": "FAILURE",
                            "error": "Consistency check failed.",
                            "model_chain": model_chain,
                            "reference_chain": reference_chain,
                            "global_score": 0.0,
                            "conserved_contacts": 0.0,
                            "total_contacts": 0.0})
                    else:
                        try:
                            ost.LogInfo((" --> Computing lDDT between model "
                                         "chain %s and reference chain %s") % (
                                             model_chain,
                                             reference_chain))
                            ost.LogInfo("Global LDDT score: %.4f" %
                                        lddt_scorer.global_score)
                            ost.LogInfo(
                                "(%i conserved distances out of %i checked, over "
                                "%i thresholds)" % (lddt_scorer.conserved_contacts,
                                                    lddt_scorer.total_contacts,
                                                    len(lddt_settings.cutoffs)))
                            sc_lddt_scores = {
                                "status": "SUCCESS",
                                "error": "",
                                "model_chain": model_chain,
                                "reference_chain": reference_chain,
                                "global_score": lddt_scorer.global_score,
                                "conserved_contacts":
                                    lddt_scorer.conserved_contacts,
                                "total_contacts": lddt_scorer.total_contacts}
                            if opts.save_per_residue_scores:
                                per_residue_sc = \
                                    mapped_lddt_scorer.GetPerResidueScores()
                                ost.LogInfo("Per residue local lDDT (reference):")
                                ost.LogInfo("Chain\tResidue Number\tResidue Name"
                                            "\tlDDT\tConserved Contacts\tTotal "
                                            "Contacts")
                                for prs_scores in per_residue_sc:
                                    ost.LogInfo("%s\t%i\t%s\t%.4f\t%i\t%i" % (
                                        reference_chain,
                                        prs_scores["residue_number"],
                                        prs_scores["residue_name"],
                                        prs_scores["lddt"],
                                        prs_scores["conserved_contacts"],
                                        prs_scores["total_contacts"]))
                                sc_lddt_scores["per_residue_scores"] = \
                                    per_residue_sc
                            lddt_results["single_chain_lddt"].append(
                                sc_lddt_scores)
                        except Exception as ex:
                            ost.LogError('Single chain lDDT failed:', str(ex))
                            lddt_results["single_chain_lddt"].append({
                                "status": "FAILURE",
                                "error": str(ex),
                                "model_chain": model_chain,
                                "reference_chain": reference_chain,
                                "global_score": 0.0,
                                "conserved_contacts": 0.0,
                                "total_contacts": 0.0})
                # perform oligo lddt scoring
                if skip_score:
                    ost.LogInfo(
                        " --> Skipping oligomeric lDDT because consistency "
                        "check failed")
                    lddt_results["oligo_lddt"] = {
                        "status": "FAILURE",
                        "error": "Consistency check failed.",
                        "global_score": 0.0}
                else:
                    try:
                        ost.LogInfo(' --> Computing oligomeric lDDT score')
                        lddt_results["oligo_lddt"] = {
                            "status": "SUCCESS",
                            "error": "",
                            "global_score": oligo_lddt_scorer.oligo_lddt}
                        ost.LogInfo(
                            "Oligo lDDT score: %.4f" %
                            oligo_lddt_scorer.oligo_lddt)
                    except Exception as ex:
                        ost.LogError('Oligo lDDT failed:', str(ex))
                        lddt_results["oligo_lddt"] = {
                            "status": "FAILURE",
                            "error": str(ex),
                            "global_score": 0.0}
                if skip_score:
                    ost.LogInfo(
                        " --> Skipping weighted lDDT because consistency "
                        "check failed")
                    lddt_results["weighted_lddt"] = {
                        "status": "FAILURE",
                        "error": "Consistency check failed.",
                        "global_score": 0.0}
                else:
                    try:
                        ost.LogInfo(' --> Computing weighted lDDT score')
                        lddt_results["weighted_lddt"] = {
                            "status": "SUCCESS",
                            "error": "",
                            "global_score": oligo_lddt_scorer.weighted_lddt}
                        ost.LogInfo(
                            "Weighted lDDT score: %.4f" %
                            oligo_lddt_scorer.weighted_lddt)
                    except Exception as ex:
                        ost.LogError('Weighted lDDT failed:', str(ex))
                        lddt_results["weighted_lddt"] = {
                            "status": "FAILURE",
                            "error": str(ex),
                            "global_score": 0.0}
                reference_results["lddt"] = lddt_results
            model_results[reference_name] = reference_results
            if opts.dump_structures:
                ost.LogInfo("-" * 80)
                ref_output_path = os.path.join(
                    os.path.dirname(opts.reference),
                    reference_name + opts.dump_suffix)
                ost.LogInfo("Saving cleaned up reference to %s" %
                            ref_output_path)
                try:
                    SavePDB(qs_scorer.qs_ent_1.ent,
                            ref_output_path)
                except Exception as ex:
                    ost.LogError("Cannot save reference: %s" % str(ex))
                mdl_output_path = os.path.join(
                    os.path.dirname(opts.model),
                    model_name + opts.dump_suffix)
                ost.LogInfo("Saving cleaned up model to %s" %
                            mdl_output_path)
                try:
                    SavePDB(qs_scorer.qs_ent_2.ent,
                            mdl_output_path)
                except Exception as ex:
                    ost.LogError("Cannot save model: %s" % str(ex))
        result["result"][model_name] = model_results

    if all(skipped) and len(skipped) > 0:
        ost.LogError("Consistency check failed for all model-reference pairs.")
    if opts.output is not None:
        ost.LogInfo("#" * 80)
        ost.LogInfo("Saving output into %s" % opts.output)
        with open(opts.output, "w") as outfile:
            json.dump(result, outfile, indent=4, sort_keys=True)

if __name__ == '__main__':
    _Main()