#!/usr/bin/env python
"""
A python module with logging, doctest to implement simulation of
reversible reflected simple random walk.

Metadata:

  - Author:    W.S.Kendall
  - Copyright: opensource (c) W.S.Kendall 2008
  - Date:      2008-07-03
  - URL:       http://www.wilfridkendall.co.uk
  - Generator: 'E:/HOME/WSK/lib/python/generate.py' version 1.18.


Command-line options::

  usage: %prog [options] arguments
    --doctest:		   perform doctests;
    --logging=LOGGING:	   logging verbosity level.

Notes:
 - values for 'logging' option:
 'DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'
 (decreasing level of verbosity), with 'WARNING' as default.

"""
__RCSID__ = '$Revision: 1.5 $ $Date: 2008/07/03 12:35:28 $'


# \section{Definitions}
# We need the logging module regardless of whether the #logging# option
# is deployed, since #logging.info# calls are scattered through the text.
import logging, math, reversible


#\section{Exceptions}
class Error(Exception):
    """
    Generic exception to be raised by this module.
    """


# \section{Classes}
class RandomWalk(reversible.Reversible):
    """
    This class implements a reversible reflected simple random walk
    with probability 'p' of a '+1' jump, probability '1-p' of a '-1'
    jump, reflected by prohibition at lower boundary 'x0' and upper
    boundary 'x1'::

     >>> chain = RandomWalk(seed=3, p=0.4) ; print chain
     []
     >>> chain.equilibrium_draw() ; print chain
     [6]
     >>> chain.evolution(10) ; print chain
     [6, 7, 6, 7, 8, 7, 6, 5, 4, 3]
     >>> print chain.reversed_statistics()
     [(3, 4), (4, 5), (5, 6), (6, 7), (7, 8), (8, 7), (7, 6), (6, 7), (7, 6)]
    """
    def __init__(self, p=0.5, x0=0, x1=10, rng=None, seed=None):
	self.p, self.x0, self.x1 = p, x0, x1
	super(RandomWalk, self).__init__(rng=rng, seed=seed)

    def equilibrium_draw(self):
	"""
	Draws from equilibrium of chain (use detailed balance to
	compute this over range of integers from 'x0' to 'x1'
	inclusive), result used to initialize chain. It makes
	numerical sense to use rejection sampling!
	"""
	# Compute \(\rho = p/(1-p)\).
	rho = self.p / (1 - self.p)
	# Separate treatment according to whether or not \(\rho \leq 1\).
	accept = False
	if rho <= 1:
	    while not accept:
		x = self.rng.choice(range(self.x0, self.x1+1))
		if self.rng.random() < math.pow(rho, (x - self.x0)):
		    accept = True
	else:
	    invrho = 1 / rho
	    while not accept:
		x = self.rng.choice(range(self.x0, self.x1+1))
		if self.rng.random() < math.pow(invrho, (self.x1 - x)):
		    accept = True
	self[:] = [x]

    def update(self):
	"""
	Updates chain by single step, which is proposed to be '+1'
	with probability 'p', '-1' with probability '1-p', with
	out-of-range transitions prohibited. The accumulated path is
	stored in 'self'.
	"""
	proposal = self[-1] + [-1, +1][self.rng.random()<self.p]
	if proposal >= self.x0 and proposal <= self.x1:
	    self.append(proposal)
	else:
	    self.append(self[-1])


# \section{Functions}
def main(opt, arguments):
    """
    Produce crude statistical analysis of a single run.
    """
    logging.info('Entering main function')
    p = 0.4
    n = 200
    chain = RandomWalk(p=p)
    chain.equilibrium_draw()
    chain.evolution(n)
    path = [a for a in chain.reversed_statistics() \
		if not a[0] == 0 and not a[0] == 10]
    print 'Crude statistical analysis: mean and 95% confidence interval.'
    print sum(map(lambda a: a[1] - a[0], path))/float(len(path))
    from math import sqrt
    mean, se = p - (1 - p), 1/sqrt(len(path))
    print '[%f, %f]' % (mean - 2 * se, mean + 2 * se)


def _doctest():
    """
    Doctests.
    """
    logging.info('Commencing doctests')
    import doctest
    doctest.testmod()
    print __RCSID__.replace('$', '').strip()


# \section{Main script}
if __name__ == '__main__':
    import optionparse
    option, args = optionparse.parse(__doc__, version=__RCSID__)

    if option.logging is not None:
	logging.basicConfig(level=logging.__dict__[option.logging])
	logging.info('Logging level %s' % option.logging)

    if option.doctest:
	_doctest()


    # The action of the module as script should be concentrated
    # in the following function, to facilitate profiling.
    main(option, args)


    logging.info('Finishing')
    raise SystemExit


# \newpage \scriptsize
# \begin{multicols}{2}
# \section{History}
# \text{ }
#\\ $Log: randomwalk.py $
#\\ Revision 1.5  2008/07/03 12:35:28  wskendall
#\\ Display confidence interval in function #main#.
#\\
#\\ Revision 1.2  2008/07/03 12:22:32  wskendall
#\\ Added doctest, also statistical test of a longer run.
#\\
#\\ Revision 1.1  2008/07/03 11:36:12  wskendall
#\\ Initial revision
#\\
#
# \end{multicols}

'$Id: randomwalk.py 1.5 2008/07/03 12:35:28 wskendall Development wskendall $'
