#!/usr/bin/env python
import sys
import numpy as np
import random 
##########################################################
# function to integrate, this function has three humps that become less
# peaked from left to right.
def f(x):
    a = 0.0
    b = 0.0
    c = 0.0
    for xi in x:
        a += -xi*xi/10.
        xi2 = xi - 20.
        b += -xi2*xi2/20.
        xi3 = xi - 40.
        c += -xi3*xi3/200.
    return np.exp(a) + 0.5 * np.exp(b) + np.exp(c)

############################################################
# proposal function
def proposal(x,delta,domain):
    x2=[]
    for xi,di,doi in zip(x,delta,domain):
        dx = di * random.uniform(-0.5,0.5)
        while not( doi[0] < xi+dx < doi[1]):
            dx = di * random.uniform(-0.5,0.5)
        x2.append(xi + dx)
    return np.array(x2)

############################################################
# acceptance rejection step
def metropolis_accept(newx,oldx,denom, invtemp):

############################################################
# metropolis step for temperatures
def heatswap(x, heat, r):

###########################################################
# swap temperatures and find index of cold chain in heat array
def swap(heat,i,j):
    tmp     = heat[i]
    heat[i] = heat[j]
    heat[j] = tmp
    return heat

def findcold(heat):
    for index in range(len(heat)):
        if heat[index]==1.0:
            return index 
    return -9999999999999;

############################################################
# mcmc
def mcmc(....):
    #loop over n steps
    #   loop over temperatures
    #         propose new value
    #         accept reject new value using temperature i
    #   pick a pair of temperatures i, j
    #   metropolis step for temperatures of i j see picture
    #   swap temperatures and find cold chain
    #   report the parameter

##########################################################
#
# main loop
# 
if __name__ == '__main__':
    #DEFAULTS user tunable
    nsteps   = 10000
    delta    = 10.0
    thinning = 1
    chains   = 4
    # defaults not user tunable
    heat = np.linspace(0.0, 1.0, num=chains)
    domain = [[-30,100],[-30,100]]
    startvalue =[-20,-20]

    # arguments
    if len(sys.argv) > 1 :
        nsteps = sys.argv[1]  # number of steps in the mcmc chain
        if "-h" in nsteps or "--help" in nsteps:
            print "python metropolis_print.py nsteps delta thinning"
            sys.exit(-1)
        else:
            nsteps = int(nsteps)
        if len(sys.argv) > 2 :
            delta  = int(sys.argv[2])  # delta for the proposal
            if len(sys.argv) > 3 :
                thinning = int(sys.argv[3])  # thinning for the proposal
                if len(sys.argv) > 4 :
                    chains = int(sys.argv[4])
                    heat = np.linspace(0.0, 1.0, num=chains)
    
    ##########################################################
    # run
    tstart     = datetime.now()
    mcmc(nsteps,[delta,delta],domain,startvalue,thinning)
    tend = datetime.now()
    print >> sys.stderr, "Elapsed time",tend-tstart
    print >> sys.stderr, "Iterations=",nsteps