PMCParallel.py

import numpy as np
import math
import os
import time
import multiprocessing as mp

def single_run_process(inputargs):
    # Run single process for multiprocessing mapping
    # Naive parallel without updating information
    tmp = list(inputargs)
    MCMCObject = tmp[0]
    loglikelihood = tmp[1]
    if MCMCObject.include and MCMCObject.exclude:
        raise ValueError("Can't specify both included and excluded parameters")

    # Setup mutable state for this run
    _ntarget = int(MCMCObject._ntarget)
    _nevaluated = 0
    _naccepted = 0
    _starttime = 0
    _lasttime = 0

    if not MCMCObject.datapars:
        # For all parameters
        data = np.zeros(
            (MCMCObject._ntarget, MCMCObject._npars), 
            dtype=float)
    else:
        # For only the given parameters
        data = np.zeros(
            (MCMCObject._ntarget, len(MCMCObject.datapars)), 
            dtype=float)

    # Build a list of parameter indices to exclude (don't shift those)
    if MCMCObject.include:
        _excluded = np.array(
            set(range(MCMCObject._npars)) - set(MCMCObject.include))
    elif MCMCObject.exclude:
        _excluded = np.array(MCMCObject.exclude)

    shifts = MCMCObject.proposal(_ntarget-1, os.getpid())    
    # Convert some values to local
    _values = MCMCObject._values

    # Get a spread in a priori
    _values += shifts[0] / MCMCObject.rescale

    # Evaluate and accept the starting point
    last_prob = loglikelihood(_values)
    _nevaluated += 1
    _naccepted += 1
    if MCMCObject.datapars:
        data[_nevaluated-1] = MCMCObject._values[MCMCObject.datapars]
    else:
        data[_nevaluated-1] = MCMCObject._values
    
    _starttime = time.time()
    
    # Loop ends when the correct number of points have been accepted
    for shift in shifts:
        # Propose a new point
        burn_in_scale = max(1.0, MCMCObject._burninrate ** (-MCMCObject._burninnum + _nevaluated))
        _values += shift * burn_in_scale

        # Evaluate the likelihood at that point
        prob = loglikelihood(_values)
        _nevaluated += 1

        # MCMC: if new likelihood > old, keep the point. Otherwise if
        # new/old > random uniform, keep
        if prob >= last_prob or \
                math.exp(prob-last_prob) > np.random.rand():
            _naccepted += 1
            last_prob = prob

        # Case where MCMC doesn't accept the point, reset to last point
        else:
            _values -= shift * burn_in_scale

        # Store the evaluted point
        if MCMCObject.datapars:
            data[_nevaluated-1] = _values[MCMCObject.datapars]
        else:
            data[_nevaluated-1] = _values
        time_left = (_ntarget - _nevaluated) / _nevaluated * (time.time()-_starttime)
        txt = 'Process {process_now:01d}, time left estimated: {time_left:.2f} hours.'
        print(txt.format(process_now=mp.current_process()._identity[0],time_left=time_left/3600))
    
    return data