#!/usr/local/bin/python2.6 -O

import copy
import optparse
import sys
import time

if "/usr/local/lib/python2.6/site-packages" not in sys.path:
    sys.path.insert(0, "/usr/local/lib/python2.6/site-packages")
import Crux

def check_uint(option, opt, value):
    try:
        ret = int(value)
        if ret < 0:
            raise ValueError()
    except ValueError:
        raise optparse.OptionValueError("Invalid value for %s: %r" % \
          (opt, value))

    return ret

def check_bool(option, opt, value):
    try:
        if value.lower() in ("false", "no"):
            ret = False
        elif value.lower() in ("true", "yes"):
            ret = True
        else:
            raise ValueError()
    except ValueError:
        raise optparse.OptionValueError("Invalid value for %s: %r" % \
          (opt, value))

    return ret

class RedpointOption(optparse.Option):
    TYPES = optparse.Option.TYPES + ("uint", "bool")
    TYPE_CHECKER = copy.copy(optparse.Option.TYPE_CHECKER)
    TYPE_CHECKER["uint"] = check_uint
    TYPE_CHECKER["bool"] = check_bool

class RedpointHelpFormatter(optparse.IndentedHelpFormatter):
    def format_usage(self, usage):
        return """\
redpoint is a script based on Crux (http://www.canonware.com/Crux/) that
analyzes a DNA alignment using Metropolis-coupled Markov chain Monte Carlo
(Mc3) to create a representative sample (assuming convergence) of the
stationary posterior distribution, based on mixtures of the GTR+I+G family of
sequence evolution models.  For more information on methods and configuration
parameters, see the relevant Crux documentation (crux --doc=Crux.Mc3.Mc3).

redpoint usage:
  redpoint -h | --help
  redpoint -V | --version
  redpoint [<options>]

  General options:
    -h, --help                   Print usage and exit.
    -V, --version                Print version information and exit.
    -v, --verbose                Enable verbose output (*disabled).
    -q, --quiet                  Disable verbose output (*enabled).
    -t, --threaded               Enable thread parallelism (*enabled).
    -u, --unthreaded             Disable thread parallelism.
    -s <uint>, --seed=<uint>     Set pseudo-random number generator seed
                                   (*based on system time, microsecond
                                   resolution).
    -G <float>, --gcmult=<float> Factor by which to modify Python's garbage
                                   collection thresholds (*1.0).

  Input/output options:
    -g <stages>,
       --stages=<stages>         Comma-separated list composed of:
                                   mc3: Run Mc3 analysis (*enabled).
                                   post: Compute statistics on posterior
                                         distribution (*enabled).
    -i <file>,                   Input DNA alignment, in FASTA format.
       --input-file=<file>
    -p <prefix>,
       --prefix=<prefix>     Input/output path prefix, used as the base name
                               for Mc3 output files.

  Mc3 options (see Crux.Mc3.Mc3.optName for documentation):

    General:                            Model parameter priors:
      --graphDelay=<float>                --rateShapeInvPrior=<float>
      --cvgSampStride=<uint>              --invarPrior=<float>
      --cvgAlpha=<float>                  --brlenPrior=<float>
      --cvgEpsilon=<float>                --rateJumpPrior=<float>
      --minStep=<uint>                    --polytomyJumpPrior=<float>
      --maxStep=<uint>                    --rateShapeInvJumpPrior=<float>
      --stride=<uint>                     --invarJumpPrior=<float>
      --nruns=<uint>                      --freqJumpPrior=<float>
      --ncoupled=<uint>                   --mixtureJumpPrior=<float>
      --heatDelta=<float>               Proposal probabilities:
      --swapStride=<uint>                 --weightProp=<float>
                                          --freqProp=<float>
    Proposal parameters:                  --rmultProp=<float>
      --ncat=<uint>                       --rateProp=<float>
      --catMedian=<bool>                  --rateShapeInvProp=<float>
      --invar=<bool>                      --invarProp=<float>
      --weightLambda=<float>              --brlenProp=<float>
      --freqLambda=<float>                --etbrProp=<float>
      --rmultLambda=<float>               --rateJumpProp=<float>
      --rateLambda=<float>                --polytomyJumpProp=<float>
      --rateShapeInvLambda=<float>        --rateShapeInvJumpProp=<float>
      --invarLambda=<float>               --invarJumpProp=<float>
      --brlenLambda=<float>               --freqJumpProp=<float>
      --etbrPExt=<float>                  --mixtureJumpProp=<float>
      --etbrLambda=<float>

  Fixed parameter overrides (appropriate proposals are implicitly disabled):
    --fixed-topology=<file>      Use the Newick tree in <file> to fix the tree
                                   topology.
    --fixed-tree=<file>          Use the Newick tree in <file> to fix the tree
                                   topology and branch lengths.
    --fixed-nmodels=<uint>       Fix the number of models in the mixture.
    --fixed-rclass=<rclass>      Fix the comma-separated rclass.  Semantics are
                                   strange unless --fixed-nmodels is also
                                   specified.
    --fixed-rates=<rates>        Fix the comma-separated rates associated with
                                   the fixed rclass.  Semantics are strange
                                   unless --fixed-nmodels is also specified.
    --fixed-shape=<alpha>        Fix the +G shape parameter.  Semantics are
                                   strange unless --fixed-nmodels is also
                                   specified.
    --fixed-pinvar=<pinvar>      Fix the +I proportion of invariable sites.
                                   Semantics are strange unless --fixed-nmodels
                                   is also specified.
    --fixed-freqs=<freqs>        Fix the comma-separated state frequencies.
                                   Semantics are strange unless --fixed-nmodels
                                   is also specified.

  Posterior distribution statistics options (*all enabled by default):
    --burnin=<burnin>  Set number of burn-in samples.  Specify --burnin=half for
                         a burn-in equal to half the total samples (*default).
    --sum=<bool>       Write summary statistics to <prefix>.sum.
    --trprobs=<bool>   Write tree topology frequencies with mean branch lengths
                         to <prefix>.trprobs.
    --parts=<bool>     Write partition frequencies and branch length statistics
                         to <prefix>.parts.
    --con=<bool>       Write consensus tree (first tree with mean branch
                         lengths, second with bipartition support values) to
                         <prefix>.con."""

# Process the command line.
def parse_options():
    helpformatter = RedpointHelpFormatter()
    parser = optparse.OptionParser(formatter=helpformatter,
      option_class=RedpointOption, add_help_option=0)
    parser.disable_interspersed_args()
    parser.add_option("-h", "--help", dest="help", action="store_true",
      default=False)
    parser.add_option("-V", "--version", dest="version", action="store_true",
      default=False)
    parser.add_option("-v", "--verbose", dest="verbose", action="store_true",
      default=False)
    parser.add_option("-q", "--quiet", dest="verbose", action="store_false")
    parser.add_option("-t", "--threaded", dest="threaded", action="store_true",
      default=True)
    parser.add_option("-u", "--unthreaded", dest="threaded",
      action="store_false")
    parser.add_option("-s", "--seed", dest="seed", type="uint",
      default=int(time.time() * 1000000) & 0xffffffff)
    parser.add_option("-G", "--gcmult", dest="gcmult", type="float",
      default=1.0)

    parser.add_option("-g", "--stages", dest="stages", default="mc3,post")
    parser.add_option("-i", "--input-file", dest="inFilename", default=None)
    parser.add_option("-p", "--prefix", dest="prefix", default=None)

    parser.add_option("--graphDelay", dest="graphDelay", type="float",
      default=None)
    parser.add_option("--cvgSampStride", dest="cvgSampStride", type="uint",
      default=None)
    parser.add_option("--cvgAlpha", dest="cvgAlpha", type="float",
      default=None)
    parser.add_option("--cvgEpsilon", dest="cvgEpsilon", type="float",
      default=None)
    parser.add_option("--minStep", dest="minStep", type="uint", default=None)
    parser.add_option("--maxStep", dest="maxStep", type="uint", default=None)
    parser.add_option("--stride", dest="stride", type="uint", default=None)
    parser.add_option("--nruns", dest="nruns", type="uint", default=None)
    parser.add_option("--ncoupled", dest="ncoupled", type="uint", default=None)
    parser.add_option("--heatDelta", dest="heatDelta", type="float",
      default=None)
    parser.add_option("--swapStride", dest="swapStride", type="uint",
      default=None)
    parser.add_option("--ncat", dest="ncat", type="uint", default=None)
    parser.add_option("--catMedian", dest="catMedian", type="bool",
      default=None)
    parser.add_option("--invar", dest="invar", type="bool", default=None)
    parser.add_option("--weightLambda", dest="weightLambda", type="float",
      default=None)
    parser.add_option("--freqLambda", dest="freqLambda", type="float",
      default=None)
    parser.add_option("--rmultLambda", dest="rmultLambda", type="float",
      default=None)
    parser.add_option("--rateLambda", dest="rateLambda", type="float",
      default=None)
    parser.add_option("--rateShapeInvLambda", dest="rateShapeInvLambda",
      type="float", default=None)
    parser.add_option("--invarLambda", dest="invarLambda", type="float",
      default=None)
    parser.add_option("--brlenLambda", dest="brlenLambda", type="float",
      default=None)
    parser.add_option("--etbrPExt", dest="etbrPExt", type="float",
      default=None)
    parser.add_option("--etbrLambda", dest="etbrLambda", type="float",
      default=None)
    parser.add_option("--rateShapeInvPrior", dest="rateShapeInvPrior",
      type="float", default=None)
    parser.add_option("--invarPrior", dest="invarPrior", type="float",
      default=None)
    parser.add_option("--brlenPrior", dest="brlenPrior", type="float",
      default=None)
    parser.add_option("--rateJumpPrior", dest="rateJumpPrior", type="float",
      default=None)
    parser.add_option("--polytomyJumpPrior", dest="polytomyJumpPrior",
      type="float", default=None)
    parser.add_option("--rateShapeInvJumpPrior", dest="rateShapeInvJumpPrior",
      type="float", default=None)
    parser.add_option("--invarJumpPrior", dest="invarJumpPrior", type="float",
      default=None)
    parser.add_option("--freqJumpPrior", dest="freqJumpPrior", type="float",
      default=None)
    parser.add_option("--mixtureJumpPrior", dest="mixtureJumpPrior",
      type="float", default=None)
    parser.add_option("--weightProp", dest="weightProp", type="float",
      default=None)
    parser.add_option("--freqProp", dest="freqProp", type="float",
      default=None)
    parser.add_option("--rmultProp", dest="rmultProp", type="float",
      default=None)
    parser.add_option("--rateProp", dest="rateProp", type="float",
      default=None)
    parser.add_option("--rateShapeInvProp", dest="rateShapeInvProp",
      type="float", default=None)
    parser.add_option("--invarProp", dest="invarProp", type="float",
      default=None)
    parser.add_option("--brlenProp", dest="brlenProp", type="float",
      default=None)
    parser.add_option("--etbrProp", dest="etbrProp", type="float",
      default=None)
    parser.add_option("--rateJumpProp", dest="rateJumpProp", type="float",
      default=None)
    parser.add_option("--polytomyJumpProp", dest="polytomyJumpProp",
      type="float", default=None)
    parser.add_option("--rateShapeInvJumpProp", dest="rateShapeInvJumpProp",
      type="float", default=None)
    parser.add_option("--invarJumpProp", dest="invarJumpProp", type="float",
      default=None)
    parser.add_option("--freqJumpProp", dest="freqJumpProp", type="float",
      default=None)
    parser.add_option("--mixtureJumpProp", dest="mixtureJumpProp", type="float",
      default=None)

    parser.add_option("--fixed-topology", dest="fixed_topology", default=None)
    parser.add_option("--fixed-tree", dest="fixed_tree", default=None)
    parser.add_option("--fixed-nmodels", dest="fixed_nmodels", type="uint",
      default=None)
    parser.add_option("--fixed-rclass", dest="fixed_rclass", default=None)
    parser.add_option("--fixed-rates", dest="fixed_rates", default=None)
    parser.add_option("--fixed-shape", dest="fixed_shape", type="float",
      default=None)
    parser.add_option("--fixed-pinvar", dest="fixed_pinvar", type="float",
      default=None)
    parser.add_option("--fixed-freqs", dest="fixed_freqs", default=None)

    parser.add_option("--burnin", dest="burnin", default="half")

    parser.add_option("--sum", dest="sum", type="bool", default=True)
    parser.add_option("--trprobs", dest="trprobs", type="bool", default=True)
    parser.add_option("--parts", dest="parts", type="bool", default=True)
    parser.add_option("--con", dest="con", type="bool", default=True)

    (opts, extra) = parser.parse_args()
    if (len(extra) > 0):
        print >> sys.stderr, "Incorrect command line usage"
        print >> sys.stderr, helpformatter.format_usage("")
        sys.exit(1)

    if opts.help:
        print helpformatter.format_usage("")
        sys.exit(0)

    if opts.version:
        print "redpoint, Crux version 1.2.0"
        sys.exit(0)

    opts.stage_mc3 = False
    opts.stage_post = False
    stages = opts.stages.split(",")
    for stage in stages:
        if stage == "mc3":
            opts.stage_mc3 = True
        elif stage == "post":
            opts.stage_post = True
        else:
            print >> sys.stderr, "Invalid stage: %s" % stage
            print >> sys.stderr, helpformatter.format_usage("")
            sys.exit(1)

    if opts.inFilename is None:
        print >> sys.stderr, "Input file must be specified"
        print >> sys.stderr, helpformatter.format_usage("")
        sys.exit(1)

    if opts.prefix is None:
        print >> sys.stderr, "Input/output prefix must be specified"
        print >> sys.stderr, helpformatter.format_usage("")
        sys.exit(1)

    opts.fixed = False
    if opts.fixed_tree is not None:
        opts.fixed = True
        opts.etbrProp = 0.0
        opts.polytomyJumpProp = 0.0
        opts.brlenProp = 0.0
        opts.tree = Crux.Tree.Tree(open(opts.fixed_tree).readline())
        opts.tree.deroot()
    elif opts.fixed_topology is not None:
        opts.fixed = True
        opts.etbrProp = 0.0
        opts.polytomyJumpProp = 0.0
        opts.tree = Crux.Tree.Tree(open(opts.fixed_topology).readline())
        opts.tree.deroot()
    else:
        opts.tree = None

    if opts.fixed_nmodels is not None:
        opts.mixtureJumpProp = 0.0

    if opts.fixed_rclass is not None:
        opts.fixed = True
        try:
            opts.fixed_rclass = [int(s) for s in opts.fixed_rclass.split(",")]
        except:
            print >> sys.stderr, "Malformed argument: --fixed-rclass=%s" % \
              opts.fixed_rclass
            print >> sys.stderr, helpformatter.format_usage("")
            sys.exit(1)
        opts.rateJumpProp = 0.0

    if opts.fixed_rates is not None:
        opts.fixed = True
        if opts.fixed_rclass is None:
            print >> sys.stderr, "--fixed-rates requires --fixed-rclass"
            print >> sys.stderr, helpformatter.format_usage("")
            sys.exit(1)
        try:
            opts.fixed_rates = [float(s) for s in opts.fixed_rates.split(",")]
        except:
            print >> sys.stderr, "Malformed argument: --fixed-rates=%s" % \
              opts.fixed_rates
            print >> sys.stderr, helpformatter.format_usage("")
            sys.exit(1)
        if len(opts.fixed_rates) != max(opts.fixed_rclass)+1:
            print >> sys.stderr, \
              "Incorrect number of fixed rates for rclass (%d vs. %d)" % \
              (len(opts.fixed_rates), max(opts.fixed_rclass)+1)
            print >> sys.stderr, helpformatter.format_usage("")
            sys.exit(1)
        opts.rmultProp = 0.0
        opts.rateProp = 0.0

    if opts.fixed_freqs is not None:
        opts.fixed = True
        # Conversion to a list of floats is done later, once the number of
        # character states is known.
        opts.freqProp = 0.0
        opts.freqJumpProp = 0.0

    if opts.fixed_shape is not None:
        opts.fixed = True
        opts.rateShapeInvProp = 0.0
        opts.rateShapeInvJumpProp = 0.0

    if opts.fixed_pinvar is not None:
        opts.fixed = True
        opts.invarProp = 0.0
        opts.invarJumpProp = 0.0

    if opts.burnin == "half":
        opts.burnin = 0xffffffffffffffff
    else:
        try:
            opts.burnin = int(opts.burnin)
            if opts.burnin < 0:
                raise ValueError()
        except:
            print >> sys.stderr, "Invalid burnin: %r" % opts.burnin
            print >> sys.stderr, helpformatter.format_usage("")
            sys.exit(1)

    return opts

def formatVecStr(s):
    strs = []
    for i in xrange(0, len(s), 5):
        strs.append(s[i:(i+5 if i+5 < len(s) else len(s))])

    return " ".join(strs)

def writePostStats(post, prefix, opts):
    if opts.sum:
        post.parseP()
        fSum = open("%s.sum" % prefix, "w")
        files = [fSum]
        if opts.verbose:
            files.append(sys.stdout)

        char_ = alignment.charType.get()
        states = char_.pcodes()

        for f in files:
            if post.burnin == 0xffffffffffffffff:
                f.write("burnin: <half>\n")
            else:
                f.write("burnin: %d\n" % post.burnin)
            f.write("nsamples: %d\n" % (post.nsamples * post.mc3.nruns))
            f.write("stepFirst: %d\n" % post.stepFirst)
            f.write("stepLast: %d\n" % post.stepLast)

            f.write("\n")

            f.write("Parameter\tRcov\t\t1.0-alpha\t1.0-alpha-epsilon\n")
            f.write("--------------------------------------------------" \
              "---------------\n")
            if post.mc3.nruns > 1:
                f.write("lnL\t\t%.6f\t%.6f\t%.6f\n" % \
                  (post.rcovLnL, 1.0 - post.mc3.cvgAlpha, \
                  1.0 - post.mc3.cvgAlpha - post.mc3.cvgEpsilon))
            else:
                f.write("lnL\t\tN/A\n")
            if post.mc3.nruns > 1 and post.maxModels == 1:
                f.write("rclass\t\t%.6f\t%.6f\n" % \
                  (post.rcovRclass, 1.0 - post.mc3.cvgAlpha))
            else:
                f.write("rclass\t\tN/A\n")

            f.write("\n")

            f.write("\t\t\t\t\t\t%.2f%% Cred. Interval\n" % \
              ((1.0 - post.mc3.cvgAlpha) * 100.0))
            f.write("\t\t\t\t\t\t-------------------------------\n")
            f.write("Parameter\tMean\t\tVariance\tLower\t\tUpper\t\tMedian\n")
            f.write("-----------------------------------------------------" \
              "-----------------------------------------\n")
            f.write("nmodels\t\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n" % \
              (post.nmodelsMean, post.nmodelsVar, post.nmodelsMinCred, \
              post.nmodelsMaxCred, post.nmodelsMed))
            f.write("lnL\t\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n" % \
              (post.lnLMean, post.lnLVar, post.lnLMinCred, post.lnLMaxCred, \
              post.lnLMed))

            # Relative mutation rates.
            ratesLabels = []
            for i in xrange(len(states)):
                for j in xrange(i+1, len(states)):
                    ratesLabels.append("r(%s<->%s)" % (states[i], states[j]))
            if post.maxModels == 1 and post.mc3.rateJumpProp == 0.0:
                ratesMinCred = post.ratesMinCred
                ratesMaxCred = post.ratesMaxCred
                ratesMean = post.ratesMean
                ratesVar = post.ratesVar
                ratesMed = post.ratesMed
                for i in xrange(len(ratesLabels)):
                    f.write("%s\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n" % \
                      (ratesLabels[i], ratesMean[i], ratesVar[i], \
                      ratesMinCred[i], ratesMaxCred[i], ratesMed[i]))
            else:
                for i in xrange(len(ratesLabels)):
                    f.write("%s\tN/A\n" % (ratesLabels[i]))

            # State frequencies.
            freqsLabels = []
            for i in xrange(len(states)):
                freqsLabels.append("pi(%s)" % states[i])
            if post.mc3.freqJumpProp == 0.0 and post.maxModels == 1:
                freqsMinCred = post.freqsMinCred
                freqsMaxCred = post.freqsMaxCred
                freqsMean = post.freqsMean
                freqsVar = post.freqsVar
                freqsMed = post.freqsMed
                for i in xrange(len(freqsLabels)):
                    f.write("%s\t\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n" % \
                      (freqsLabels[i], freqsMean[i], freqsVar[i], \
                      freqsMinCred[i], freqsMaxCred[i], freqsMed[i]))
            else:
                for i in xrange(len(states)):
                    f.write("%s\t\tN/A\n" % freqsLabels[i])

            # +G shape.
            if post.maxModels == 1 and post.mc3.ncat > 1 \
              and post.mc3.rateShapeInvJumpProp == 0.0:
                f.write("+G alpha\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n" % \
                  (post.alphaMean, post.alphaVar, post.alphaMinCred, \
                  post.alphaMaxCred, post.alphaMed))
            else:
                f.write("+G alpha\tN/A\n")

            # +I pinvar.
            if post.maxModels == 1 and post.mc3.invar \
              and post.mc3.invarJumpProp == 0.0:
                f.write("+I pinvar\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n" % \
                  (post.pinvarMean, post.pinvarVar, post.pinvarMinCred, \
                  post.pinvarMaxCred, post.pinvarMed))
            else:
                f.write("+I pinvar\tN/A\n")

        fSum.close()

    if opts.trprobs:
        trprobs = post.sumt.trprobs
        if len(trprobs) > 0:
            f = open("%s.trprobs" % prefix, "w")
            f.write("[[#obs p P MSE #runs] tree]\n")

            P = 0.0
            for trprob in trprobs:
                p = float(trprob.nobs) / float(post.mc3.nruns * post.nsamples)
                P += p
                f.write("[%d %f %f %f %d] " % (trprob.nobs, p, P, trprob.mse, \
                  trprob.nruns))
                trprob.tree.render(lengths=True, lengthFormat="%.12e", \
                  outFile=f)

            f.close()

    if opts.parts:
        parts = post.sumt.parts
        if len(parts) > 0:
            f = open("%s.parts" % prefix, "w")
            part = parts[0]
            plen = max(len(formatVecStr("%r" % part.vec)), len("Partition"))
            olen = max(len("%d" % (post.mc3.nruns * post.nsamples)), \
              len("#obs"))

            # Print taxa key.
            taxa = post.runs[0][0].tree.taxa
            for i in xrange(len(taxa)):
                taxon = taxa[i]
                if i % 5 == 0 and i != 0:
                    f.write("%s\n" % \
                      (formatVecStr(("|    "*((i+(5-(i%5)))/5))[:i+1])[:-1]) \
                      .rstrip())
                f.write("%s%s\n" % \
                  (formatVecStr(("|    "*((i+(5-(i%5)))/5))[:i+1])[:-1], \
                  taxon.label))
            i = len(taxa) - 1
            pstr = \
              (formatVecStr(("|    "*((i+(5-(i%5)))/5))[:i+1])[:-1]).rstrip()
            f.write("%s\n" % pstr)

            f.write( \
              "%-*s  %-*s\tProbability\tMSE\t\tmean(v)\t\tvar(v)\t\t#runs\n" % \
              (plen, pstr, olen, "#obs"))

            for part in parts:
                f.write("%-*s  %*d\t%f\t%f\t%.5e\t%.5e\t%d\n" % \
                  (plen, formatVecStr("%r" % part.vec), olen, part.nobs, \
                  float(part.nobs) / float(post.mc3.nruns * post.nsamples), \
                  part.mse, part.mean, part.var, part.nruns))

            f.close()

    if opts.con:
        conTree = post.sumt.getConTree()
        f = open("%s.con" % prefix, "w")
        # Render with mean branch lengths.
        conTree.render(lengths=True, lengthFormat="%.12e", outFile=f)

        # Render with bipartition support values.
        for edge in conTree.edges:
            edge.length = edge.aux
        conTree.render(lengths=True, lengthFormat="%.6f", outFile=f)
        f.close()

def run(alignment, prefix, opts):
    mc3 = Crux.Mc3.Mc3(alignment, prefix)

    # Set optional Mc3 configuration parameters.
    if opts.graphDelay is not None: mc3.graphDelay = opts.graphDelay
    if opts.cvgSampStride is not None: mc3.cvgSampStride = opts.cvgSampStride
    if opts.cvgAlpha is not None: mc3.cvgAlpha = opts.cvgAlpha
    if opts.cvgEpsilon is not None: mc3.cvgEpsilon = opts.cvgEpsilon
    if opts.minStep is not None: mc3.minStep = opts.minStep
    if opts.maxStep is not None:
        # Automatically reduce the default minStep if necessary to avoid a
        # validation error.
        if opts.minStep is None and mc3.minStep > opts.maxStep:
            mc3.minStep = opts.maxStep
        mc3.maxStep = opts.maxStep
    if opts.stride is not None: mc3.stride = opts.stride
    if opts.nruns is not None: mc3.nruns = opts.nruns
    if opts.ncoupled is not None: mc3.ncoupled = opts.ncoupled
    if opts.heatDelta is not None: mc3.heatDelta = opts.heatDelta
    if opts.swapStride is not None: mc3.swapStride = opts.swapStride
    if opts.fixed_nmodels is not None: mc3.nmodels = opts.fixed_nmodels
    if opts.ncat is not None: mc3.ncat = opts.ncat
    if opts.catMedian is not None: mc3.catMedian = opts.catMedian
    if opts.invar is not None: mc3.invar = opts.invar
    if opts.weightLambda is not None: mc3.weightLambda = opts.weightLambda
    if opts.freqLambda is not None: mc3.freqLambda = opts.freqLambda
    if opts.rmultLambda is not None: mc3.rmultLambda = opts.rmultLambda
    if opts.rateLambda is not None: mc3.rateLambda = opts.rateLambda
    if opts.rateShapeInvLambda is not None:
        mc3.rateShapeInvLambda = opts.rateShapeInvLambda
    if opts.invarLambda is not None:
        mc3.invarLambda = opts.invarLambda
    if opts.brlenLambda is not None: mc3.brlenLambda = opts.brlenLambda
    if opts.etbrPExt is not None: mc3.etbrPExt = opts.etbrPExt
    if opts.etbrLambda is not None: mc3.etbrLambda = opts.etbrLambda
    if opts.rateShapeInvPrior is not None:
        mc3.rateShapeInvPrior = opts.rateShapeInvPrior
    if opts.invarPrior is not None:
        mc3.invarPrior = opts.invarPrior
    if opts.brlenPrior is not None: mc3.brlenPrior = opts.brlenPrior
    if opts.rateJumpPrior is not None: mc3.rateJumpPrior = opts.rateJumpPrior
    if opts.polytomyJumpPrior is not None:
        mc3.polytomyJumpPrior = opts.polytomyJumpPrior
    if opts.rateShapeInvJumpPrior is not None:
        mc3.rateShapeInvJumpPrior = opts.rateShapeInvJumpPrior
    if opts.invarJumpPrior is not None:
        mc3.invarJumpPrior = opts.invarJumpPrior
    if opts.freqJumpPrior is not None: mc3.freqJumpPrior = opts.freqJumpPrior
    if opts.mixtureJumpPrior is not None:
        mc3.mixtureJumpPrior = opts.mixtureJumpPrior
    if opts.weightProp is not None: mc3.weightProp = opts.weightProp
    if opts.freqProp is not None: mc3.freqProp = opts.freqProp
    if opts.rmultProp is not None: mc3.rmultProp = opts.rmultProp
    if opts.rateProp is not None: mc3.rateProp = opts.rateProp
    if opts.rateShapeInvProp is not None:
        mc3.rateShapeInvProp = opts.rateShapeInvProp
    if opts.invarProp is not None:
        mc3.invarProp = opts.invarProp
    if opts.brlenProp is not None: mc3.brlenProp = opts.brlenProp
    if opts.etbrProp is not None: mc3.etbrProp = opts.etbrProp
    if opts.rateJumpProp is not None: mc3.rateJumpProp = opts.rateJumpProp
    if opts.polytomyJumpProp is not None:
        mc3.polytomyJumpProp = opts.polytomyJumpProp
    if opts.rateShapeInvJumpProp is not None:
        mc3.rateShapeInvJumpProp = opts.rateShapeInvJumpProp
    if opts.invarJumpProp is not None:
        mc3.invarJumpProp = opts.invarJumpProp
    if opts.freqJumpProp is not None: mc3.freqJumpProp = opts.freqJumpProp
    if opts.mixtureJumpProp is not None:
        mc3.mixtureJumpProp = opts.mixtureJumpProp

    if opts.fixed:
        # --fixed-freqs has to be handled here since the number of character
        # states wasn't yet known in parse_options().
        if opts.fixed_freqs is not None:
            try:
                opts.fixed_freqs = [float(s) for s in \
                  opts.fixed_freqs.split(",")]
            except:
                print >> sys.stderr, "Malformed argument: --fixed-freqs=%s" % \
                  opts.fixed_freqs
                helpformatter = RedpointHelpFormatter()
                print >> sys.stderr, helpformatter.format_usage("")
                sys.exit(1)
            char_ = alignment.charType.get()
            states = char_.pcodes()
            if len(opts.fixed_freqs) != len(states):
                print >> sys.stderr, \
                  "Incorrect number of fixed rates for rclass (%d vs. %d)" % \
                  (len(opts.fixed_freqs), len(states))
                helpformatter = RedpointHelpFormatter()
                print >> sys.stderr, helpformatter.format_usage("")
                sys.exit(1)

        liks = []
        for i in xrange(mc3.nruns * mc3.ncoupled):
            lik = mc3.randomLik(opts.tree)
            for m in xrange(mc3.nmodels):
                if opts.fixed_rates is not None:
                    lik.setRclass(m, opts.fixed_rclass, opts.fixed_rates)
                elif opts.fixed_rclass is not None:
                    lik.setRclass(m, opts.fixed_rclass)
                if opts.fixed_shape is not None:
                    lik.setAlpha(m, opts.fixed_shape)
                if opts.fixed_pinvar is not None:
                    lik.setWVar(m, 1.0-opts.fixed_pinvar)
                    lik.setWInvar(m, opts.fixed_pinvar)
                if opts.fixed_freqs is not None:
                    for f in xrange(len(opts.fixed_freqs)):
                        lik.setFreq(m, f, opts.fixed_freqs[f])
            liks.append(lik)
    else:
        liks = None

    mc3.run(verbose=opts.verbose, liks=liks)

#===============================================================================
# Beginning of main execution.

# Parse command line options.
opts = parse_options()

if opts.gcmult != 1.0:
    import gc
    threshold = gc.get_threshold()
    gc.set_threshold(int(threshold[0] * opts.gcmult),
      int(threshold[1] * opts.gcmult), int(threshold[2] * opts.gcmult))

Crux.Config.verbose = opts.verbose

# Set threading mode.
if opts.threaded:
    Crux.threaded()

# Seed the PRNG.
Crux.seed(opts.seed)

cTMatrix = Crux.CTMatrix.CTMatrix(open(opts.inFilename, "r"))
alignment = Crux.CTMatrix.Alignment(matrix=cTMatrix)
alignment.compact()

if opts.stage_mc3:
    run(alignment, opts.prefix, opts)

if 0:
    # Only one node needs to do post-processing.
    from mpi4py import MPI
    rank = MPI.COMM_WORLD.Get_rank()
else:
    rank = 0

if opts.stage_post and rank == 0:
    post = Crux.Mc3.Post.Post(alignment, opts.prefix, opts.burnin, \
      opts.verbose, lazy=True)
    post.parseS()
    if post.nsamples == 0:
        print >> sys.stderr, "No samples"
        sys.exit(1)
    writePostStats(post, opts.prefix, opts)
