#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# MEDELLER : template-based membrane protein structure prediction
#
# Author: Sebastian Kelm
# Created: 10/09/2009
#
#



import sys
import os
import tempfile
import shutil

if __name__ == "__main__":
  # Make sure we can import stuff
  sys.path.append(os.path.abspath(os.path.dirname(sys.argv[0])+"/../lib-python"))

from prosci.common import *
from prosci.util.pdb import Pdb, diffPdb
from prosci.util.residue import ResidueList
from prosci.util.protein import Protein
from prosci.util.ali import Ali
from prosci.util.gaps import count_nongaps, deGappify
from prosci.medeller.corebuilder import CoreBuilder, SubstTableSelector, calculateAccuracy
from prosci.medeller.fread import getFreadPath, modelAllGaps, calculateCoreCoverage, countMaxModellableResidues, mapCore2TargetSeq, DatabaseDependencyError
from prosci.loops.fread import FREAD_RANKING, FREAD_RANKING_BY_ESSS

def write_corresponding_native_structure(native, target_residue_indeces, model, outfname_native, outfname_model, model_fname):
    if native is not None:
      native = native.to_residuelist()
    if model is not None:
      model = model.to_residuelist()
    "BENCHMARKING: Write a file containing the native structure residues corresponding to the core model"
    if native is not None and os.path.exists(model_fname):
      out = ResidueList([])
      for i, res in zip(target_residue_indeces, model):
        native_res = native[i].copy()
        assert native_res.res == res.res
        native_res.ires = res.ires
        native_res.inscode = res.inscode
        native_res.chain = res.chain
        out.append(native_res)
      write_file(outfname_native, str(out))
      
      model = model.deep_copy()
      for i, res in zip(target_residue_indeces, model):
        native_res = native[i]
        assert native_res.res == res.res
        res.ires = native_res.ires
        res.inscode = native_res.inscode
        res.chain = native_res.chain
      write_file(outfname_model, str(model))



COMMON_HICOV_LOOP_STRATEGY = ("mem:sol", "mem|sol", "mem", "sol")


if __name__ == "__main__":
    
    
    FREAD_IS_INSTALLED = bool(getFreadPath())
    
    
    import prosci.shell
    
    ################################
    # Command line option handling #
    ################################
    
    # DELETED OPTIONS
    # Without options: 'history', 'loopaccuracy', 'readonly_loops', 'readonly_models',
    # With options: 'model', 'target', 
    
    params = prosci.shell.Params(
       allowed=['all', 'pairs', 'help', 'S', 'W', 'F', 'M', '1', '2', '3', '4', '5', '6', 'v', 'b', 'heuristic', 'exhaustive', 'nomerge', 'atm', 'verbose', 'naive', 'simple', 'refine_sidechains', 'scwrl', 'clash', 'debug'],
       required=['a', 'i'],
       withargument=['a', 'i', 'list', 'd', 'prefix', 'substtables', 'minimal_core', 'exclude_pdbs', 'native_structure', 'ligands', 'contact_identity', 'hiacc_loop_strategy', 'hicov_loop_strategy', 'hiacc_anchor_rmsd', 'hicov_anchor_rmsd', 'hiacc_score', 'hicov_score', 'pool_size', 'chaincode'],
       helpoption='help',
       usage="""
USAGE:
  
  Use all templates in alignment file
    $0 [OPTIONS] -a <alignment> -i <target_id> --all
  
  Use only specified templates.
    $0 [OPTIONS] -a <alignment> -i <target_id>
        <template_id>...
  
  Specify templates and their associated structures.
    $0 [OPTIONS] -a <alignment> -i <target_id>
        --pairs {<template_id> <template_struc>}...
  
  Specify a list of templates and their structures, read from a separate file:
    $0 [OPTIONS] -a <alignment> -i <target_id>
        --list <template_list>
      \n""",
       
       help="""
OPTIONS:
  
  REQUIRED INPUT:

    -a FILE   FILE is an alignment file in ALI format.
              Templates should be aligned and pre-annotated with
              JOY: "secondary structure and phi angle", and
              iMembrane: "membrane layer (extended)".
    
    -i ID     ID is the entry code, in the above ALI file, for
              the target protein.
    
  
  ALTERNATIVE INPUT:
  
    --list FILE    A whitespace-separated file listing the templates to
                   use, one per line, in the format:
                      ID   /path/to/file.pdb
                      ID   /path/to/file.pdb
                      ...
  
  
  MISC:
    
    --help
        Print this message and exit
    
    --substtables DIRECTORY
        Specifies the path to the directory containing the substitution
        tables.
        (Optional. Defaults to a directory called "substtables"
        parallel to the directory where this script is located.)

    -d DIR
        Search for input files in this directory, before searching
        the current working directory. (This only applies to filenames
        provided as relative paths.)
        Separate multiple directory names using the colon ':' character.
    
    --prefix PREFIX
        Sets the prefix for output files to PREFIX.
        Default for a single template:
          TARGETID_TEMPLATEID
        Default for multiple templates (N=number of templates):
          TARGETID_Ntemplates
    
    --chaincode X
        Sets the chain code of the output structure to X (default: 'A').
    
    --atm
          Output files with the .atm extension, instead of .pdb
    
    --verbose
          Print status messages
  
  
  BENCHMARKING:
    
    --exclude_pdbs PDB1,PDB2,PDB3,...
        Exclude loop fragments coming from the given PDBs. Argument is a comma-
        separated list of PDB codes (only the first 4 characters of each ID are
        used).
    
    --native_structure FILE
        Specify the PDB structure of target protein chain. This will result in
        files ending in .native.pdb which contain the equivalent target residues
        to each MEDELLER model. Files ending .renumbered.pdb will also be
        created, containing the model residues renumbered to correspond to the
        native structure's numbering.
        NOTE: The given native structure must correspond exactly to the target
        sequence given in the input alignment.
    
  
  ALGORITHM BEHAVIOUR:
    --clash
        Don't check for clashes in FREAD loop predictions.
    
    --pool_size INT
        Number of processes to spawn for FREAD loop modelling. Default is 1. Do
        not set this higher than the number of processor cores of your machine.
    
    --ligands RESTYPE[,RESTYPE...]
        Copy ligands of the specified types from the template(s) into the model.
        If the empty string "" is specified, all ligands will be copied.
    
    -contact_identity FLOAT
        Calculate contacts during loop modelling and match contacts of database
        fragments with the contacts in the model being built. Only loops with
        at least FLOAT (fraction [0,1]) of its contacts being conserved are
        considered. Specify 0 if you only wish to see what contacts exist for
        each decoy, without letting it affect the model. A negative number
        disables contact calculation (which is the default). Use of this option
        requires the FREAD database bein searched to contain contact
        information.
    
    --hiacc_loop_strategy MODE
        Strategy to follow for filling gaps in the High Accuracy stage.
        By default the "mem" database is searched. (Default: mem)
    
    --hiacc_anchor_rmsd FLOAT
        Set the anchor RMSD cut-off for loop modlling in the High Accuracy
        model. Only smaller values <= 1.0 are recommended. (Default: 1.0)
    
    --hiacc_score INT
        Set the substitution score cut-off for loop modlling in the
        High Accuracy model. Increasing this cut-off will reduce both runtime
        and coverage, while decreasing it has the opposit effect. (Default: 25)
    
    --hicov_loop_strategy MODE
        (e.g. """+",".join(COMMON_HICOV_LOOP_STRATEGY)+"""; default: """+COMMON_HICOV_LOOP_STRATEGY[0]+""")
        Strategy to follow for filling gaps in the High Coverage stage.
        By default the "mem" database is searched and, if no hits were found,
        the "sol" is searched. This increases the chances of filling all gaps.
    
    --hicov_anchor_rmsd FLOAT
        Set the anchor RMSD cut-off for loop modlling in the high-coverage
        model. Increasing this above the default 1.0 will require a compiled
        version of CCD (source code shipped with MEDELLER). Higher cut-offs
        generally increase the amount of low quality loop decoys found. Good
        loops should still rank higher than bad ones, meaning the main effect
        of increasing the cut-off will be longer run times and higher coverage
        "high coverage" models.
      
    --hicov_score INT
        Set the substitution score cut-off for loop modlling in the
        high-coverage model. Increasing this cut-off will reduce both runtime
        and coverage, while decreasing it has the opposit effect. Default is 1.
      
    --scwrl
        Model missing sidechains using SCWRL4 or SCWRL3. The executable 'Scwrl4'
        or 'scwrl3' must be in your medeller/bin folder or in your PATH
        environment variable. Models without SCWRL sidechains are still output
        and contain the word "noscwrl" in the filename.
    
    --refine_sidechains
        Refine sidechains using Modeller. This makes the final modelling phase
        (phase 6) much slower, but will reduce sidechain clashes. Use this if
        you want decent sidechains but don't have SCWRL installed. This only
        has an effect if you are creating a complete model. When using the
        --scwrl option, refined Modeller sidechains will be overwritten with
        SCWRL sidechains, so the effect on the model may be negligible.
    
    --naive
        Build the core as a naive model (accepting the backbone coordinates of
        all aligned residues), without using the substitution score or column
        masking rules. This is useful for benchmarking.
        You may want to combine this option with the options -4, -5 and -6 to
        switch off any loop modelling.
    
    --simple
        Like --naive, but restrict copying to residues fulfilling the
        "minimal core" criteria. This is the same as -123456M combined with
        setting the minimal core word cut-off to 0 (by default 2 residues
        are cut off each end of each minimal core fragment).
    
    --minimal_core ANNOTTYPE:VALUES
        Specify the annotation type and values used to define the minimal core.
        These must be present in standard TEM format in the input alignment.
        For example:
        
        "membrane layer:T" would select all residues in the membrane tail layer
        This is the default behaviour.
        
        "membrane layer:TH" would select all residues in the membrane tail and
        head layers. This is the behaviour produced by the "-b" option.
    
    --nomerge
        Never merge overlapping fragments. This will prevent any errors caused
        by averaging coordinates, but may result in a smaller core.
        (Only applicable when modelling with multiple templates.)
    
    --heuristic    (enabled by default)
        Use a fast heuristic to choose fragments included in the core
        model. This will not necessarily give the best possible core,
        but is very fast.
        (Only applicable when modelling with multiple templates.)
    
    --exhaustive
        Use an exhaustive search to find the best-scoring combination
        of fragments to include in the core model. This is slower than
        the default heuristic method and the speed and memory usage
        depend on the number of templates and conserved fragmets found.
        Most times, the result is the same as for the heuristic method,
        but it occasionally makes a difference.
        (Only applicable when modelling with multiple templates.)
    
    -b    
        Start from a bigger minimal core, which includes residues in 
        the polar head group layer (H) of the membrane.
        This may result in higher coverage, but is more likely to
        include regions that do not span the entire membrane and are
        less structurally conserved.
    
    -S    
        Disable use of secondary structure information.
        Use membrane layer annotation only.
        This disables loop masking as well.
    
    -W  Disable window-based smoothing
    -F  Disable fragment score. Score only the added residue
    -M  Disable masking of long loops and loops containing gaps
    
    -1  Disable phase 1 : copying secondary structures
    -2  Disable phase 2 : copying long loops
    -3  Disable phase 3 : copying gappy loops + terminal loops
    -4  Disable phase 4 : high accuracy database search loop modelling
    -5  Disable phase 5 : high coverage database search loop modelling
    -6  Disable phase 6 : ab initio loop modelling
  
  
  DEBUGGING:
  
    -v  
        Be very verbose. Print lots of debug info to STDOUT
        (Errors and warnings are still printed to STDERR as normal)

    --debug
        Enable debugging mode. Only developers should use this option.

  
  DEFAULT OUTPUT FILES:
  
    Output filenames will have the format PREFIX.SUFFIX
    For details on PREFIX, see --prefix
    SUFFIX is one of the following:
      
    
      SUFFIX                  FILE CONTENT
      
      core.pdb                Core model structure
                              (no loop modelling)
      
      hiacc.messages          Log output of high accuracy FREAD
      hiacc.summary           Summary of loop models found by hiacc FREAD
      hiacc.pdb               High accuracy model structure
                              (core + high accuracy database loops)
      
      hicov.messages          Log output of high coverage FREAD
      hicov.summary           Summary of loop models found by hicov FREAD
      hicov.pdb               High coverage model structure
                              (hiacc + low accuracy db loops)
      
      complete.pdb            Complete model structure
                              (hicov + ab initio loops)
      \n""")
    
    
    if (not params.args and not params.isOpt("all")) or \
       (params.args and params.isOpt("all")) or \
       (len(params.args)%2 and params.isOpt("pairs")) or \
       (params.args and params.isOpt("list")):
      sys.stderr.write("""
Incorrect usage. Must specify either
  * any number of IDs as arguments, or
  * '--all' and no arguments, or
  * '--list FILENAME' and no arguments, or
  * '--pairs' and an even number of arguments, i.e.: ID filename ID filename ...
      """.strip())
      sys.stderr.write("\n")
      params.writeUsage()
      sys.exit(0)
    
    params.setDefaultOptions(d=None, hiacc_loop_strategy="mem", hicov_loop_strategy=COMMON_HICOV_LOOP_STRATEGY[0], contact_identity=-1, hiacc_anchor_rmsd=1.0, hicov_anchor_rmsd=1.0, hiacc_score=25, hicov_score=1, pool_size=1, minimal_core="membrane layer:T", chaincode="A")
    
    DEBUG = params.isOpt("debug")

    if params.isOpt('verbose'):
      print "Options:"
      for k in sorted(params.opts):
        print "   ", k, "=", params.opts[k]
      print "Positional arguments:"
      for arg in params.args:
        print "  ", arg
    
    pool_size = int(params.getOpt("pool_size"))
    contact_identity = float(params.getOpt("contact_identity"))
    hicov_loop_strategy = params.getOpt("hicov_loop_strategy")
    hicov_anchor_rmsd = float(params.getOpt("hicov_anchor_rmsd"))
    hicov_score = int(params.getOpt("hicov_score"))
    hiacc_loop_strategy = params.getOpt("hiacc_loop_strategy")
    hiacc_anchor_rmsd = float(params.getOpt("hiacc_anchor_rmsd"))
    hiacc_score = int(params.getOpt("hiacc_score"))
    
    assert len(params.getOpt("chaincode")) == 1
    
    annotation_type, annotation_values = params.getOpt("minimal_core").split(":")
    if params.isOpt('b'):
      annotation_values += "H"
    
    
    template_searchdirs = params.getOpt('d')
    if template_searchdirs:
      template_searchdirs = template_searchdirs.strip(':').split(':')
      for i,dir in enumerate(template_searchdirs):
        template_searchdirs[i] = os.path.normpath(dir)
    else:
      template_searchdirs = []
      
    
    def resolve(fname):
      if template_searchdirs and not os.path.isabs(fname):
        for dir in template_searchdirs:
          if os.path.exists("%s/%s"%(dir, fname)):
            return os.path.abspath("%s/%s"%(dir, fname))
      return os.path.abspath(fname)
    
    
    alignment_filename = params.getOpt('a')
    target_id          = params.getOpt('i')
    subst_table_dir    = params.getOpt("substtables", os.path.dirname(params.scriptdir) + "/substtables")
    
    assert os.path.exists(subst_table_dir)
    
    
    if params.isOpt('verbose'):
      print "Parsing input files"

    
    
    # Read in alignment file
    alignment = Ali(file(resolve(alignment_filename)))
    assert len(alignment), "Error parsing ALI (or TEM) alignment file. No data?"
    
    # Template selection
    #
    template_id_list    = []
    template_struc_list = []
    if params.isOpt("pairs"):
      for i,x in enumerate(params.args):
        if i%2:
          template_struc_list.append(x)
        else:
          template_id_list.append(x)
    elif params.isOpt("list"):
      for line in file(resolve(params.opts["list"])):
        line = line.strip()
        fields = line.split()
        template_id_list.append(fields[0])
        template_struc_list.append(fields[1])
    else:
      for x in alignment.get_entries_by_type("structure"):
        if x.code == target_id:
          continue
        if not params.isOpt("all") and x.code not in params.args:
          continue
        template_id_list.append(x.code)
        template_struc_list.append(x.getStructureFilename(template_searchdirs))
    
    # Remove annoying MODELLER-style "P1;" headings
    target_id.replace("P1;", "")
    for i,x in enumerate(template_id_list):
      template_id_list[i] = x.replace("P1;", "")
    for eg in alignment:
      for entry in eg:
        entry.code = entry.code.replace("P1;", "")
        
    
    
    outsuffix = "pdb"
    if params.isOpt("atm"):
      outsuffix = "atm"
    
    # Output filenames
    if params.isOpt("prefix"):
      outprefix = params.opts["prefix"]
    else:
      if len(template_id_list) == 1:
        outprefix = "%s_%s" %(target_id, template_id_list[0])
      else:
        outprefix = "%s_%dtemplates" %(target_id, len(template_id_list))
    
    outdir = os.path.dirname(os.path.abspath(outprefix))
    if not os.path.exists(outdir):
      os.makedirs(outdir)
    
    
    def make_filename(shortname, suffix=None):
      if suffix is None:
        suffix = outsuffix
      return "%s.%s.%s" %(outprefix, shortname, suffix)
    
    
    
    ############
    # SETTINGS #
    ############
    
    exclude_pdbs = params.getOpt("exclude_pdbs", "").split(',')
    if not exclude_pdbs[0]:
      exclude_pdbs = None
    else:
      pool_size = 1
      print "DEBUG: Setting pool_size = 1. Multiprocessing is not compatible with the PDB exclusion filter"
    
    wordcutoffs=2
    minwordlength=2
    
    if not params.isOpt('W'):
      windowLeft=1
      windowRight=1
    else:
      windowLeft=0
      windowRight=0
    
    maskLoops=(not params.isOpt('M'))
    useFragmentScore=(not params.isOpt('F'))
    
    phasesToPerform = [1, 2, 3]
    if params.isOpt('1'):
      phasesToPerform.remove(1)
    if params.isOpt('2'):
      phasesToPerform.remove(2)
    if params.isOpt('3'):
      phasesToPerform.remove(3)
    
    if params.isOpt("naive") or params.isOpt("simple"):
      wordcutoffs = 0
      maskLoops = False
      phasesToPerform = []
    
    
    # Resolve template filenames
    #
    for i,fname in enumerate(template_struc_list):
      template_struc_list[i] = resolve(fname)
    
    # Add directories where templates reside to search path (for Modeller)
    #
    for i,fname in enumerate(template_struc_list):
      template_searchdirs.append(os.path.dirname(fname))
    
    # Read template structures
    #
    for i,fname in enumerate(template_struc_list):
      template_struc_list[i] = Protein(fname, code=template_id_list[i])
      assert len(template_struc_list[i]), "Error parsing template PDB file '%s'. No data?"%(fname)
    
    
    if params.isOpt('verbose'):
      print "Initialising CoreBuilder"
    
    # Initialise MEDELLER's main datastructure
    #
    try:
      coreBuilder = CoreBuilder(alignment[target_id], [alignment[x] for x in template_id_list], template_struc_list, template_directory=template_searchdirs, score_in_occupancy=True, useFragmentScore=useFragmentScore, chaincode=params.getOpt("chaincode"))
    except IndexError, error:
      sys.stderr.write("IndexError: %s\n"%error)
      sys.exit(1)
    
    
    # Length of target
    # Used to write sequence files with the appropriate number of terminal gaps
    #
    UNGAPPED_TARGET_SEQ = deGappify(coreBuilder.getTargetSeq())
    TARGET_LENGTH = len(UNGAPPED_TARGET_SEQ)
    
    
    # BENCHMARKING: Load the native structure, if given
    #
    native_pdb=None
    if params.isOpt("native_structure"):
      native_pdb = Protein(resolve(params.getOpt("native_structure")))
      if native_pdb.get_seq() != UNGAPPED_TARGET_SEQ:
        sys.stderr.write("Error: Sequence of native target structure must be identical to target's sequence in the alignment!:\nStruc:%s\nAlign:%s\n"%(native_pdb.get_seq(), UNGAPPED_TARGET_SEQ))
        sys.exit(2)
    
    
    
    coreBuilder.setVerbose(params.isOpt('v'))
    coreBuilder.setWindowSize(windowLeft, windowRight)
    #coreBuilder.setVerbose(False)
    
    
    ###################################
    # Read in the Substitution Tables #
    ###################################

    if params.isOpt('verbose'):
      print "Loading substitution tables"
    
    if params.isOpt('S'):
        # Membrane Layer annotation only
        #
        
        fileprefix = "%s/jsubst.layer.ENV_"%(subst_table_dir)
        filesuffix = "-LOG_ODDS.mat"
        envs   = ["membrane layer"]
        values = ["NHT?"]
        env2esst = dict()
        for i in xrange(len(values[0])):
          val = values[0][i]
          env2esst[val] = fileprefix + val.replace("?", "_") + filesuffix
        
        # We have no secondary structure information.
        # Thus, we cannot mask loops - we don't know where they are!
        #
        maskLoops = False
    else:
        # Membrane Layer + Secondary Structure annotation
        #
        fileprefix = "%s/jsubst.2struc_layer.ENV_"%(subst_table_dir)
        filesuffix = "-LOG_ODDS.mat"
        envs   = ["membrane layer", "secondary structure and phi angle"]
        values = ["NHT?", "CHEP"]
        env2esst = dict()
        for i in xrange(len(values[0])):
          for j in xrange(len(values[1])):
            val = values[0][i] + values[1][j]
            env2esst[val] = fileprefix + val.replace("?", "_") + filesuffix
    
    ESSTs = SubstTableSelector(env2esst, envs)
    coreBuilder.setSubstTables(ESSTs, maskLoops)
    
    
    
    ##################
    # Model the core #
    ##################
    
    if params.isOpt('verbose'):
      print "Building model core"
    
    # Allow absence of substitution tables, if we're not doing any gap extension
    if not phasesToPerform:
      try:
        for c in coreBuilder.cc:
          c.calculateSubstitutionScores()
      except Ali.InvalidEntryError:
        for c in coreBuilder.cc:
          c.setDummySubstitutionScores()
    
    if params.isOpt("naive"):
      # Select all aligned residues
      coreBuilder.selectAll()
    else:
      # Select minimal core using the annotation specified by annotation_type and annotation_values
      coreBuilder.selectMinimalCore(wordcutoffs, minwordlength, annotation_type, annotation_values)
      
      # Extend core model using MEDELLER algorithm
      coreBuilder.extendCore(phasesToPerform)
      
    # Allow overlapping fragments
    #
    if not params.isOpt("nomerge"):
      coreBuilder.createMergedFragments() # use default merging settings
    
    if params.isOpt("exhaustive"):
      coreBuilder.chooseFragmentsExhaustive()
    else: #if params.isOpt("heuristic"):
      coreBuilder.chooseFragmentsHeuristic()
    
    #else:
      #coreBuilder.chooseFragmentsAuto(cutoff=30)
    
    #coreBuilder.setVerbose(True)
    coreModel, target_residue_indeces = coreBuilder.buildCoreModel()
    coreModel = Protein(coreModel)
    #coreBuilder.setVerbose(False)
    
    
    coreModel.code = target_id + "_core"
    
    if params.isOpt("ligands"):
      ligand_types = params.getOpt("ligands").upper().split(",")
      if len(ligand_types) == 1 and len(ligand_types[0]) == 0:
        ligand_types = []
      if params.isOpt('verbose'):
        print "Copying ligands"
      ligs = ResidueList([])
      for template in template_struc_list:
        for chain in template.ligands:
          for res in chain:
            if (not ligand_types) or (res.res in ligand_types):
              newres = res.copy()
              newres.chain = "Z"
              ligs.append(res)
              if params.isOpt('verbose'):
                print "    Adding ligand", res.res
            else:
              if params.isOpt('verbose'):
                print "    Not adding ligand", res.res
      coreModel.add_ligands(ligs)
    
    # Write core structure file
    #
    #if not params.isOpt("readonly_models"):
    write_file(make_filename("core"), str(coreModel))
    
    
    target_residue_indices, model_alignment = mapCore2TargetSeq(UNGAPPED_TARGET_SEQ, coreModel, renumber=True, targetid=target_id)
    
    
    
    # BENCHMARKING: Write a file containing the native structure residues corresponding to the core model
    #
    write_corresponding_native_structure(native_pdb, target_residue_indices, coreModel, make_filename("core.native"), make_filename("core.renumbered"), make_filename("core"))
    
    
    latestModel = coreModel
    
    
    
    #####################################
    # Fill gaps in the model with FREAD #
    #####################################
    
    
    def fill_gaps_with_fread(inshortname, outshortname, renumber=False, **kwargs):
      "Run FREAD with specified parameters, write corresponding native structure."
      
      outfilename = make_filename(outshortname)
      
      if DEBUG and os.path.exists(outfilename):
        model = Protein(outfilename, code=outshortname)
        target_residue_indices, model_alignment = mapCore2TargetSeq(UNGAPPED_TARGET_SEQ, model, renumber=renumber, targetid=target_id)
        return model, target_residue_indices

      infilename = make_filename(inshortname)
      
      model = Protein(infilename)
      
      target_residue_indeces, model_alignment = mapCore2TargetSeq(UNGAPPED_TARGET_SEQ, model, renumber=renumber, targetid=target_id)
      
      vdw_factor=0.7
      if params.isOpt("clash"):
        vdw_factor = 0.0
      
      try:
        modelAllGaps(infilename, UNGAPPED_TARGET_SEQ, target_residue_indeces, outfilename, make_filename(outshortname, "decoys"), make_filename(outshortname, "messages"), summary_file_prefix=make_filename(outshortname, "summary"), exclude_pdbs=exclude_pdbs, pool_size=pool_size, vdw_factor=vdw_factor, **kwargs)
      except DatabaseDependencyError, e:
        print e
      
      model = Protein(outfilename, code=outshortname)
      
      target_residue_indeces, model_alignment = mapCore2TargetSeq(UNGAPPED_TARGET_SEQ, model, renumber=renumber, targetid=target_id)
      
      # BENCHMARKING: Write a file containing the native structure residues corresponding to the model
      #
      write_corresponding_native_structure(native_pdb, target_residue_indeces, model, make_filename(outshortname+".native"), make_filename(outshortname+".renumbered"), outfilename)
      
      return model, target_residue_indeces
    
    
    
    if params.isOpt('verbose'):
      if not FREAD_IS_INSTALLED and not (params.isOpt('4') and params.isOpt('5')):
        print "NOTICE: FREAD is not installed. Cannot produce high-accuracy and high-coverage models."
    
    
    if not params.isOpt('4') and FREAD_IS_INSTALLED:
      if params.isOpt('verbose'):
        print "Modelling loops (high-accuracy)"
      # Model loops with high accuracy
      latestModel, target_residue_indices = fill_gaps_with_fread("core", "hiacc", renumber=True, score_cutoff=hiacc_score, anchor_rmsd=hiacc_anchor_rmsd, max_extension=1000000, minb=4, dbtype=hiacc_loop_strategy, ranking=FREAD_RANKING, contact_identity=contact_identity)
    
    
    if not params.isOpt('5') and FREAD_IS_INSTALLED:
      if params.isOpt('verbose'):
        print "Modelling loops (high-coverage)"
      # Model loops with high coverage
      latestModel, target_residue_indices = fill_gaps_with_fread("hiacc", "hicov", renumber=False, score_cutoff=hicov_score, anchor_rmsd=hicov_anchor_rmsd, max_extension=1000000, minb=5, dbtype=hicov_loop_strategy, ranking=FREAD_RANKING_BY_ESSS, contact_identity=contact_identity)
    
    
    
    ###################################
    # Completing model using MODELLER #
    ###################################
    
    
    if not params.isOpt('6'):
      try:
        from prosci.medeller.modellerif import completeModel as completeAllAtom
      except ImportError:
        if params.isOpt('verbose'):
          print "NOTICE: Modeller is not installed. Cannot produce complete model."
      else:
        if params.isOpt('verbose'):
          print "Completing model"
          if not FREAD_IS_INSTALLED:
            print "NOTICE: FREAD is not installed. All loops are modelled using Modeller."
        
        outdir = tempfile.mkdtemp()
        try:
          # Write all pre-parsed templates to a temporary directory
          #
          # Modify alignment to delete any residues in the template sequence
          # that doesn't have corresponding co-ordinates in the PDB file
          #
          for t, tid in zip(template_struc_list, template_id_list):
            write_file(os.path.join(outdir, t.code + ".pdb"),str(t))
            annotation = alignment[tid]
            residues_with_structure = set(t.to_residuelist().map_to_seq(annotation.master.seq))
            seq = list(annotation.master.seq)
            for i, c in enumerate(seq):
              if i not in residues_with_structure:
                seq[i] = "-"
            annotation.master.seq = "".join(seq)

          template_searchdirs = [outdir]
          
          # Complete the model using MODELLER
          completeModel = completeAllAtom(latestModel, alignment, target_id, template_ids=template_id_list, template_directory=template_searchdirs, outFileName=make_filename("complete"), bfactor=15, refine_sidechains=params.isOpt("refine_sidechains"), debug=DEBUG)
          latestModel = completeModel
          latestModel.code = "complete"
        finally:
          shutil.rmtree(outdir)
        
        target_residue_indeces, model_alignment = mapCore2TargetSeq(UNGAPPED_TARGET_SEQ, completeModel, targetid=target_id)
        
        
        # BENCHMARKING: Write a file containing the native structure residues corresponding to the hicov model
        #
        write_corresponding_native_structure(native_pdb, target_residue_indeces, completeModel, make_filename("complete.native"), make_filename("complete.renumbered"), make_filename("complete"))
        
        
        latestModel = completeModel



    #####################################
    # Modelling side chains using SCWRL #
    #####################################
    if params.isOpt("scwrl"):
      try:
        from prosci.medeller.scwrl import modelMissingSidechains
      except ImportError, e:
        print "NOTICE:", e
      else:
        for shortname in ("core", "hiacc", "hicov", "complete"):
          if os.path.exists(make_filename(shortname)):
            model = Protein(make_filename(shortname))
            if params.isOpt('verbose'):
              print "Modelling sidechains (%s)" % (shortname)
            modelWithSidechains = modelMissingSidechains(model)
            if len(modelWithSidechains) < len(model):
              sys.stderr.write("Warning: Something went wrong during sidechain modelling for %s model\n" % (shortname))
            else:
              shutil.move(make_filename(shortname), make_filename(shortname+".noscwrl"))
              write_file(make_filename(shortname), str(modelWithSidechains))



if native_pdb is not None:
  print "Target length: %d" % (TARGET_LENGTH)
