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.


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

    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:

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.


        Terminates a group of operations.


        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.

    (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.

Metropolis operations:

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

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

    met-1 [ 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.

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


        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.


        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

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

    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] ] 

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

Slicing operations:

    slice-1 [ 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 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

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

        Do overrelaxed slice sampling, for the range of components 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

Tempering/annealing operations:

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


        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.


        Randomize the tempering direction.


        Negate the tempering direction.


        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.  


        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 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 "firstq".  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.


        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, 1998, 1999 by Radford M. Neal