#!/usr/bin/env python
"""
A python module with logging, doctest to implement simulation of
reversible M/M/1 queueing processes.

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.1 $ $Date: 2008/07/03 20:22:32 $'


# \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, bdi


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


# \section{Classes}
class MM1(bdi.BDI):
    """
    This class implements a reversible M/M/1 queueing process
    with rates as follows: immigration 'lambd', service 'mu'::

     >>> chain = MM1(seed=2) ; print chain
     []
     >>> chain.equilibrium_draw()
     (0, 0)
     >>> chain.evolution(1.5) ; print ['%4.2f %d' % x for x in chain]
     ['0.00 0', '0.05 1', '0.88 0', '1.06 1', '1.19 2', '1.36 1', '1.50 1']
     >>> print ['%4.2f %d %d' % x for x in chain.reversed_statistics()]
     ['0.17 1 2', '0.13 2 1', '0.18 1 0', '0.82 0 1', '0.05 1 0']
    """
    def __init__(self, lambd=1, mu=2, rng=None, seed=None):
        self.lambd, self.mu = lambd, mu
        super(bdi.BDI, self).__init__(rng=rng, seed=seed)

    def equilibrium_draw(self):
        """
        Draws from equilibrium of chain, result used to initialize
        chain. The equilibrium distribution is geometric: number of
        failures before first success for independent success
        probability 'lambd/mu'. Note that state is tuple of time and
        place::

         >>> chain = MM1(seed=1)
         >>> N = 10000
         >>> result = {}
         >>> for t, n in [chain.equilibrium_draw() for i in range(N)]:
         ...   if result.has_key(n):
         ...      result[n] += 1
         ...   else:
         ...      result[n] = 1.0
         >>> res = sorted([(i, '%4.2f' % (f / N)) for i, f in result.items()])
         >>> [x[1] for x in res][:7]
         ['0.50', '0.26', '0.12', '0.06', '0.03', '0.02', '0.01']
        """
        n, rho = 0, float(self.lambd) / self.mu
        while self.rng.random() < rho:
            n += 1
        self[:]=[(0, n)]
        return (0, n)

    def update(self):
        """
        Updates chain by single step, which is proposed to be '+1'
        with probability 'lambd', '-1' with probability 'mu' if 'x' is
        non-zero, with time-increment being exponential of rate 'lambd
        + mu' if 'x' is non-zero, 'lambd' if 'x' is zero. The
        accumulated path is stored in 'self'.
        """
        t, x = self[-1]
        if self[-1][1] == 0:
            rate = self.lambd
            x += 1
        else:
            rate = self.lambd + self.mu
            if rate * self.rng.random() < self.mu:
                x -= 1
            else:
                x += 1

        t += self.rng.expovariate(rate)
        self.append((t, x))
       

# \section{Functions}
def main(opt, arguments):
    """
    Make a single run.
    """
    logging.info('Entering main function')
    T = 2
    chain = MM1()
    chain.equilibrium_draw()
    chain.evolution(T)
    print ['(%5.3f %d %d)' % x for x in chain.reversed_statistics()]


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: mm1.py,v $
#\\ Revision 1.1  2008/07/03 20:22:32  wskendall
#\\ Initial revision
#\\
#\\ 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: mm1.py,v 1.1 2008/07/03 20:22:32 wskendall Exp wskendall $'
