    MC-SPEC:  Specify how to do the Markov chain simulation.

    MC-spec writes records to a log file that specify the operations
    making up a single Markov chain iteration, and that also specify how
    dynamic trajectories are to be computed.  When invoked with just a log
    file as argument, it displays the specifications stored at the given
    index, or the last specifications stored if no index is given.

    MC-spec is often invoked several times during a single simulation run.
    The last specification in the log file is used for further iterations.

    Usage:

        mc-spec log-file { operation-spec } [ / trajectory-spec ]

    or:

        mc-spec log-file [ index | "all" ]

    The log file must already exist.  For the first form, specifications
    are appended to the log file with index equal that of the last record
    in the file, or 0 if the log file has no non-negative records.

    For the second form, with no index specified, the last specifications
    stored in the log file are displayed.  If an index is specified in the
    second form, the specifications in the record with the given index are
    displayed.  If the string "all" (no quotes) is given after the log
    file, specifications at all indexes are displayed.

    The allowed operations can be grouped in the categories general,
    metropolis, dynamical, tempering, and slicing.  Note that the slicing
    operations should not be used when coupling chains, since they can
    cause the random number sequence for different chains to get out of
    synchronization.

    General operations:

        repeat times

            Repeat the following operations the specified number of times.
            The repeated group of operations is terminated by 'end', or
            the end of the list.

        end

            Terminates a group of operations.

        plot

            Prints the current values of those quantities specified in the
            "xxx-mc" command on standard output, preceded by an index that
            is reset to zero every iteration.

        multiply-stepsizes factor [ first[:last] ] ]

            Multiplies the stepsizes for coordinates in the range first:last
            (default is all if first:last is omitted, just first if last is
            omitted) by the given factor.  This modification persists for
            the remainder of the operations in the current iteraton, unless
            modified again, or unless the automatic stepsizes are overridden
            by a stepsize specification starting with "-".  (This last provision
            is just to make implementation easier.)

            If a multiply-stepsizes operation is inside a "repeat", it will
            be done repeatedly.  In particular, if factor is less than one,
            the stepsizes will be successively smaller for each repetition.

        (any operation not otherwise defined here) [ number [ number ] ]

            Invoke the application-specific update procedure, passing the
            name of the operation, and the numerical parameters following
            it, which default to zero.

            These operations may or may not be allowed when coupling.

    Metropolis/Gibbs operations:

    For descriptionsof the stepsize and the "-b" argument, see below.

        metropolis [ -b ] [ stepsize-adjust[:stepsize-alpha] ]

            Do a simple Metropolis update, to all components at once, using
            a Gaussian proposal distribution.

        met [ -b ]  [ stepsize-adjust[:stepsize-alpha] ]

            Same as "metropolis".

        met-1 [ -b ] [ -r ] [ stepsize-adjust[:stepsize-alpha] [ first[:last] ] ]

            Do single-component Metropolis updates, for the range of components
            given.  The default if no components are given is all components.
            The default for last if first alone is given is last=first.  If no
            -r argument is given, the components are all updated in order; with
            the -r option, just one randomly-chosen component is updated.

        rgrid-met [ -b ] [ stepsize-adjust[:stepsize-alpha] ]

            Do a random-grid Metropolis update, to all components at once, using
            a uniform proposal distribution centred on the current state and
            extending to plus and minus the stepsize specified.  This uniform
            proposal is found by randomly positioning a grid of hypercubes and
            then proposing a move to the centre of the hypercube containing the
            current state.  This leads to exact coalescence of coupled chains
            whenever both chains' states are in the same hypercube.

        rgrid-met-1 [ -b ] [ -r ] [ stepsize-adjust[:stepsize-alpha]
                                    [ first[:last] ] ]

            Do single-component random-grid Metropolis updates, for the range
            of components given.  The default if no components are given is all
            components. The default for last if first alone is given is last=first.
            If no -r argument is given, the components are all updated in order;
            with the -r option, just one randomly-chosen component is updated.

        gaussian-gibbs [ -r ] [ first[:last] ] ]

            Do Gibbs sampling updates for the range of corrdinates given
            (default, is all coordinates; last defaults to first), or do a
            Gibbs sampling update for one of these coordinates chosen at
            random, if -r is specified.  The Gibbs sampling is done by just
            ASSUMING that the conditional distribution is Gaussian, with
            its mean and variance being found by evaluating the energy at
            values of -1, 0, and +1 for the coordinate being updated.  This
            WILL NOT BE VALID if the conditional distribution is not Gaussian.

    Dynamical operations:

        heatbath [ decay ]

            Do a heatbath update of the momentum variables.  If decay is
            zero (the default), the current momentum is forgotten, and
            new values are picked randomly from their distribution.  If
            decay is non-zero, the momentum variables are multiplied by
            decay, and Gaussian noise with variance 1-decay^2 is then
            added.

        radial-heatbath

            Do a radial heatbath update of the momentum variables, in
            which the squared magnitude of the momentum is sampled
            from a chi-squared distribution, while the direction is
            left unchanged.

        mix-momentum factor

            Does a random rotation of all momentum variables, by adding
            mean zero Gaussian noise to each momentum variable, and then
            multiplying all by the factor needed to produce the same
            kinetic energy as before.  The variance of the noise added is
            the average squared magnitude of the current momentum variables
            times the specified factor squared.

        negate

            Negate the momentum variables.

        dynamic steps [ stepsize-adjust[:stepsize-alpha] ]

            Follow a dynamical trajectory for the given number of steps,
            always accepting the result (hence does not leave distribution
            exactly invariant).

        permuted-dynamic steps [ stepsize-adjust[:stepsize-alpha] ]

            Like 'dynamic', but the order of approximations is randomly
            permuted.

        hybrid steps[:window[:jump]] [ stepsize-adjust[:stepsize-alpha] ]

            OR

        hybrid max-steps/max-ok[:jump] [ stepsize-adjust[:stepsize-alpha] ]

            Use the results of following a dynamical trajectory as a
            candidate state for a Metropolis update.  There are two forms,
            differing in the way length of a trajectory and its acceptance
            are determined.  In both forms, states along a trajectory are
            looked at only at every 'jump' steps (and either 'steps' or
            'max-steps' must be multiples of 'jump').  The default is a jump
            of one.  The stepsize is determined as for the 'dynamic' operation.
            In both forms, the momentum is negated if the proposal is accepted
            in order to make the step reversible.

            In the first form, acceptance is based on a 'window' of states
            at the beginning and end of the trajectory.  In the second form,
            acceptance is based on a single state, but the number of steps
            in the trajectory is not fixed - instead, the trajectory ends
            after max-ok states that would be accepted have been found (looking
            only every 'jump' states), or when max-steps states of any sort have
            been produced; if the trajectory ends for the first reason, acceptance
            is guaranteed.  The default if neither 'max-ok' nor 'window' is
            specified is standard hybrid Monte Carlo (ie, a window of one).

        tempered-hybrid temper-factor steps[:window[:jump]]
                      [ stepsize-adjust[:stepsize-alpha] ]

            Like hybrid, but the trajectory is "tempered" by multiplying
            the momenta by temper-factor in the first half of the portion
            of the trajectory outside any windows, and by 1/temper-factor
            in the second half.

        spiral temper-factor steps [ stepsize-adjust[:stepsize-alpha] ]

            Performs spiral dynamics, in which steps of Hamiltonian dynamics
            alternate with multiplication/division of the momentum by the
            specified tempering factor.  The specified number of steps are
            randomly divided into steps "before" and "after" the current
            state.

        double-spiral temper-factor steps [ stepsize-adjust[:stepsize-alpha] ]

            Like spiral, but with a reversal of direction at a randomly
            chosen point, producing a double spiral.

    Slicing operations:

        slice-1 [ -r ] [ -s [-]factor[/threshold] ]
                       [ stepsize-adjust[:stepsize-alpha]
                       [ max-steps [ first[:last] ] ] ]

            Do single-component slice sampling, for the range of components
            given.  The default for the range of components (first:last) is
            all components.  The default for last if first alone is given is
            last=first.  Components are numbered starting with zero.  If no -r
            argument is given, the components are all updated in order; with
            the -r option, just one randomly-chosen component is updated.

            If the max-steps argument is positive, it gives the maximum number
            of intervals to create when using the stepping out procedure to
            find an interval around the current point.  A value of zero
            indicates that stepping out should be done with no limit (this
            is the default).  If max-steps is negative, the interval is found
            by doubling, with minus max-steps again being the limit on the
            number of intervals (so setting max-steps to -1 gives the same
            results as setting it to 1).  (There is no way to double without
            limit, but setting max-steps to -1000 is the same for all practical
            purposes.)

            With no -s option, when a point is rejected, one end of the interval
            shrinks to the rejected point.  If the -s option is given, with no
            "-" before "factor", the interval may then be shrunk further by
            dividing it up into "factor" subintervals, and retaining only the
            subinterval that contains the current point.  This additional
            shrinkage is done if the energy is more than "threshold" above the
            slice level (the default is a threshold of zero).  If "-" is given
            before "factor", the initial shrinkage to the rejected point is not
            done - only the shrinkage by "factor" is done.  As a consequence,
            shrinkage of the interval can be disabled entirely using "-s -1".

        slice [ -g | -G ] [ stepsize-adjust[:stepsize-alpha] ]

            Do multi-dimensional slice sampling with the given stepsize.  At
            present, the initial hyperrectangle found cannot be expanded (either
            by stepping out or by doubling).  After a point is randomly chosen
            from the hyperrectangle, it is shrunk in certain directions if the
            chosen point is not in the slice, and another random point chosen,
            until eventually a point in the slice is found.  If neither the -g
            nor the -G option is specified, all coordinate directions are shrunk.
            If -g is specified, only one coordinate direction is shrunk - that
            for which the product of the absolute value of the energy gradient
            and the current dimension of the hyperrectange in that coordinate
            is the maximum.  If -G is specified, this coordinate direction
            is shrunk (as for -g), along with all other coordinate directions
            for which the corresponding product is at least half this maximum.

        slice-gaussian [ -e ] [ stepsize-adjust[:stepsize-alpha] ]

            Do multi-dimensional slice sampling with Gaussian-distributed
            "crumbs".  The distribution of the first crumb has diagonal
            covariance, with standard deviations given by the stepsizes.
            If no option is given, this is also the distribution of later
            crumbs.  When -e is specified, later crumbs have their standard
            deviations rescaled based on the energy of the most recent trial
            point, so as to (hopefully) lead to a good chance of the next
            trial point lying in the slice.

        slice-over [ -r ] [ refinements [ refresh-prob
          [ stepsize-adjust[:stepsize-alpha] [ max-steps [ first[:last] ] ] ] ] ]

            Do overrelaxed slice sampling, for the range of components given,
            or for one such component chosen randomly, if the -r option is given.
            The endpoints are computed using the given number of refinements
            (default zero).  The refresh probability is the probability of doing
            an ordinary slice sampling update rather than an overrelaxed one;
            it defaults to zero.  The default range of components if no range is
            given is all components. The default for last if first alone is given
            is last=first.  The meaning of max-steps is the same as for slice-1.

        slice-inside steps [ stepsize-adjust[:stepsize-alpha] ]

            Performs multivariate slice sampling by reflection from inside
            points, proceeding for the indicated number of steps, with the
            indicated stepsize adjustment factor.  The momentum is negated
            in such a way as to make the operation reversible.

        slice-outside steps[/in-steps] [ stepsize-adjust[:stepsize-alpha] ]

            Performs multivariate slice sampling by reflection from outside
            points, proceeding for the indicated number of steps, with the
            indicated stepsize adjustment factor. The momentum is negated
            at the end in such a way as to make the operation reversible.
            The in-steps argument gives a limit on the number of steps that
            can be inside the slice; it defaults to the same as steps (thereby
            having no effect).  Setting in-steps to less than steps may decrease
            the chances of the trajectory ending on an outside point, and being
            rejected.

    Tempering/annealing operations:

    See also the descriptions under "dynamical operations" of the spiral,
    double-spiral, and tempered-hybrid operations.

        sim-temp

            Do a metropolis update of the simulated tempering inverse
            temperature, with a proposal of changing the temperature
            index in accord with the current tempering direction.  The
            direction is negated if the proposal is accepted to make the
            step reversible.

        rand-dir

            Randomize the tempering direction.

        neg-dir

            Negate the tempering direction.

        temp-trans

            Perform a tempered transition, with components given by the
            following operations (terminated by 'end').  The components
            are done in forward order for the first half of the tempered
            transition, in reverse order for the second.  For this to
            work correctly, the components must be reversible in themselves.

        AIS

            Proceed to the next step of an annealed importance sampling
            run, with the tempering index increased by one, adjusting the
            importance weight appropriately.  If the tempering index is
            already at the final value (inverse temperature of one), a new
            state is randomly generated from the distribution at inverse
            temperature zero, and the tempering index is set to the beginning
            of the tempering schedule.

    Test operations:

        The following operations are used for testing and research only.
        They do not leave the desired distribution invariant, and hence
        should not be used as part of a sampling scheme.

        multiply-momentum factor

            All momentum variables are multiplied by the given factor.

        set-momentum value

            All momentum variables are set to the given value.

    Certain of these operations are normally used in standard combinations,
    in particular the following:

     heatbath hybrid <steps>                 Standard hybrid Monte Carlo
     heatbath <decay> hybrid <steps> negate  Persistent Hybrid Monte Carlo
     rand-dir sim-temp                       Standard simulated tempering
     sim-temp neg-dir                        Persistent Simulated tempering

    All operations are reversible (other than 'ais', 'repeat', and 'end'
    for which the concept is not applicable), except for 'dynamic',
    'permuted-dynamic', 'multiply-momentum', 'set-momentum', and perhaps
    the application-specific operations.  However, note that in general
    sequential combinations of reversible operations are not reversible.

    Defaults are decay of zero, stepsize-adjust of one, stepsize-alpha of
    infinity, window of one, and jump of one.  A minus sign before a
    stepsize-adjust value results in the adjustment being applied to
    uniform stepsizes of one; otherwise, the adjustment is applied to the
    application-specific stepsizes (which may be non-uniform).  A minus
    sign in front of stepsize alpha means that instead of the usual Gamma
    distribution for stepsizes, a distribution over 'alpha' orders of
    magnitude, uniform in the log domain, centred at the value found using
    stepsize-adjust, is used.

    The "-b" option for Metropolis operations causes the "Barker" or
    "Boltzmann" acceptance probability to be used, rather than the usual
    "Metropolis" form.  If R is the ratio of probability densities for the
    new and old states, the "-b" option acceptance probability is R/(1+R),
    whereas the Metropolis form is min(1,R).

    The trajectory specification can have one of the following forms:

        leapfrog [ halfp | halfq ] [ N-approx ]

            Use the leapfrog method, repeated N-approx times using a randomly
            selected ordering of that number of energy approximations (whose
            average must be the true value).  The default for N-approx is one.
            The "halfp" or "halfq" option specifies whether the half-steps
            at the beginning and end of the trajectory are for p (momentum)
            or q (position).  The default is "halfp".

        opt2 [ rev | sym ] [ firstq | firstp ]

            Apply McLachlan and Atela's "optimal" two stage method, with the
            first step being done for the position (firstq) or momentum
            (firstp) variables.  The default is "firstp".  This method is
            not symmetrical.  The "rev" option gives the reversed version.
            The "sym" option gives the ordinary followed by the reversed
            version, with half the stepsize for each, thereby producing
            a symmetrical method.  Only the symmetrical method should be
            used in conjuntion with "hybrid" and "tempered hybrid" operations.

        gen2 [ rev | sym ] [ firstq | firstp ] a2

            Apply a two stage second-order method with parameter a2, as
            described by McLachlan and Atela.  Firstq and firstp are as for
            opt2.  This is probably only interesting for research, since the
            most useful values for a2 are covered by "leapfrog" and "opt2".
            The "rev" and "sym" options are as for opt2.

        opt4 [ rev | sym ]

            Apply McLachlan and Atela's "optimal" four stage method for
            quadratic kinetic energy. The "rev" and "sym" options are as for opt2.

    The default if no specification is given is "leapfrog halfp 1".  The
    opt2, gen2, and opt4 methods are taken from a paper by Robert
    I. McLachlan and Pau Atela, "The accuracy of symplectic integrators",
    Nonlinearity, vol. 5, pp. 541-562, 1992.  (Note: there is a typo in
    their description of opt2 in table 2. b1=1-b2 should be 1-1/sqrt(2).)

    The list of operations is stored in a record of type 'o'; the
    trajectory specification in a record of type 't'.

                Copyright (c) 1995-2004 by Radford M. Neal