# Path: pyrate_lib/lib_utilities.py
# def calcHPD(data, level=0.95) :
#     assert (0 < level < 1)    
#     d = list(data)
#     d.sort()    
#     nData = len(data)
#     nIn = int(round(level * nData))
#     if nIn < 2 :
#         raise RuntimeError("not enough data")    
#     i = 0
#     r = d[i+nIn-1] - d[i]
#     for k in range(len(d) - (nIn - 1)) :
#         rk = d[k+nIn-1] - d[k]
#         if rk < r :
#             r = rk
#             i = k
#     assert 0 <= i <= i+nIn-1 < len(d)    
#     return np.array([d[i], d[i+nIn-1]])
# 
# Path: pyrate_lib/lib_utilities.py
# def print_R_vec(name,v):
#     new_v=[]
#     if len(v)==0: vec= "%s=c()" % (name)
#     elif len(v)==1: vec= "%s=c(%s)" % (name,v[0])
#     elif len(v)==2: vec= "%s=c(%s,%s)" % (name,v[0],v[1])
#     else:
#         for j in range(0,len(v)): 
#             value=v[j]
#             if isnan(v[j]): value="NA"
#             new_v.append(value)
#         
#         vec="%s=c(%s, " % (name,new_v[0])
#         for j in range(1,len(v)-1): vec += "%s," % (new_v[j])
#         vec += "%s)"  % (new_v[j+1])
#     return vec
# 
# Path: pyrate_lib/lib_utilities.py
# def get_mode(data):
#     # determine bins Freedman-Diaconis rule
#     iqr = np.subtract(*np.percentile(data,[75,25])) # interquantile range
#     h_temp = 2 * iqr * len(data) **(-1./3) # width
#     h = max(h_temp,0.0001)
#     n_bins= int((max(data)-min(data))/h)
#     hist=np.histogram(data,bins=n_bins)
#     return hist[1][np.argmax(hist[0])] # modal value


import argparse, os,sys
import numpy as np
import scipy.special
import random as rand
import platform, time
import multiprocessing, _thread
import multiprocessing.pool
import os, csv, glob
import scipy.stats
import _thread
import pyrate_lib.lib_utilities as lib_utilities
from numpy import *
from scipy.special import gamma
from scipy.special import beta as f_beta
from scipy.special import gdtr, gdtrix
from scipy.special import betainc
from multiprocessing import Pool, freeze_support
from pyrate_lib.lib_updates_priors import *
from pyrate_lib.lib_DD_likelihood  import *
from pyrate_lib.lib_utilities import calcHPD as calcHPD
from pyrate_lib.lib_utilities import print_R_vec as print_R_vec
from pyrate_lib.lib_utilities import get_mode as get_mode
    
            marginal_L.append(trasfRate_general(baseline_L, G_temp[fixed_focal_clade,0,:],Dtraj))
            marginal_M.append(trasfRate_general(baseline_M, G_temp[fixed_focal_clade,1,:],Dtraj))
            if j % 100 ==0: 
                sys.stdout.write(".")
                sys.stdout.flush()

        if i== -1: print("\nCalculating mean rates and HPDs...")            
        else: print("\nProcessing variable:", variable_names[i])
        
        marginal_L = np.array(marginal_L)
        marginal_M = np.array(marginal_M)
        #print np.shape(marginal_L)

        l_vec= np.zeros(np.shape(marginal_L)[1])
        m_vec= np.zeros(np.shape(marginal_L)[1])
        hpd_array_L= np.zeros((2,np.shape(marginal_L)[1]))
        hpd_array_M= np.zeros((2,np.shape(marginal_L)[1]))
        hpd_array_L50= np.zeros((2,np.shape(marginal_L)[1]))
        hpd_array_M50= np.zeros((2,np.shape(marginal_L)[1]))
        
        if i>=0:
            l_vec = np.mean(marginal_L, axis=0) # get_mode
            m_vec = np.mean(marginal_M, axis=0) # get_mode
        else:        
            for ii in range(np.shape(marginal_L)[1]): # loop over marginal rates
                l_vec[ii] = np.mean(marginal_L[:,ii]) # get_mode
                m_vec[ii] = np.mean(marginal_M[:,ii]) # get_mode
                hpd_array_L[:,ii] = calcHPD(marginal_L[:,ii])
                hpd_array_M[:,ii] = calcHPD(marginal_M[:,ii])

# Path: pyrate_lib/lib_utilities.py
# def calcHPD(data, level=0.95) :
#     assert (0 < level < 1)    
#     d = list(data)
#     d.sort()    
#     nData = len(data)
#     nIn = int(round(level * nData))
#     if nIn < 2 :
#         raise RuntimeError("not enough data")    
#     i = 0
#     r = d[i+nIn-1] - d[i]
#     for k in range(len(d) - (nIn - 1)) :
#         rk = d[k+nIn-1] - d[k]
#         if rk < r :
#             r = rk
#             i = k
#     assert 0 <= i <= i+nIn-1 < len(d)    
#     return np.array([d[i], d[i+nIn-1]])
# 
# Path: pyrate_lib/lib_utilities.py
# def print_R_vec(name,v):
#     new_v=[]
#     if len(v)==0: vec= "%s=c()" % (name)
#     elif len(v)==1: vec= "%s=c(%s)" % (name,v[0])
#     elif len(v)==2: vec= "%s=c(%s,%s)" % (name,v[0],v[1])
#     else:
#         for j in range(0,len(v)): 
#             value=v[j]
#             if isnan(v[j]): value="NA"
#             new_v.append(value)
#         
#         vec="%s=c(%s, " % (name,new_v[0])
#         for j in range(1,len(v)-1): vec += "%s," % (new_v[j])
#         vec += "%s)"  % (new_v[j+1])
#     return vec
# 
# Path: pyrate_lib/lib_utilities.py
# def get_mode(data):
#     # determine bins Freedman-Diaconis rule
#     iqr = np.subtract(*np.percentile(data,[75,25])) # interquantile range
#     h_temp = 2 * iqr * len(data) **(-1./3) # width
#     h = max(h_temp,0.0001)
#     n_bins= int((max(data)-min(data))/h)
#     hist=np.histogram(data,bins=n_bins)
#     return hist[1][np.argmax(hist[0])] # modal value


# Path: PyRateMBD.py
import argparse, os,sys
import numpy as np
import scipy.special
import random as rand
import platform, time
import multiprocessing, _thread
import multiprocessing.pool
import os, csv, glob
import scipy.stats
import _thread
import pyrate_lib.lib_utilities as lib_utilities
from numpy import *
from scipy.special import gamma
from scipy.special import beta as f_beta
from scipy.special import gdtr, gdtrix
from scipy.special import betainc
from multiprocessing import Pool, freeze_support
from pyrate_lib.lib_updates_priors import *
from pyrate_lib.lib_DD_likelihood  import *
from pyrate_lib.lib_utilities import calcHPD as calcHPD
from pyrate_lib.lib_utilities import print_R_vec as print_R_vec
from pyrate_lib.lib_utilities import get_mode as get_mode

    
            marginal_L.append(trasfRate_general(baseline_L, G_temp[fixed_focal_clade,0,:],Dtraj))
            marginal_M.append(trasfRate_general(baseline_M, G_temp[fixed_focal_clade,1,:],Dtraj))
            if j % 100 ==0: 
                sys.stdout.write(".")
                sys.stdout.flush()

        if i== -1: print("\nCalculating mean rates and HPDs...")            
        else: print("\nProcessing variable:", variable_names[i])
        
        marginal_L = np.array(marginal_L)
        marginal_M = np.array(marginal_M)
        #print np.shape(marginal_L)

        l_vec= np.zeros(np.shape(marginal_L)[1])
        m_vec= np.zeros(np.shape(marginal_L)[1])
        hpd_array_L= np.zeros((2,np.shape(marginal_L)[1]))
        hpd_array_M= np.zeros((2,np.shape(marginal_L)[1]))
        hpd_array_L50= np.zeros((2,np.shape(marginal_L)[1]))
        hpd_array_M50= np.zeros((2,np.shape(marginal_L)[1]))
        
        if i>=0:
            l_vec = np.mean(marginal_L, axis=0) # get_mode
            m_vec = np.mean(marginal_M, axis=0) # get_mode
        else:        
            for ii in range(np.shape(marginal_L)[1]): # loop over marginal rates
                l_vec[ii] = np.mean(marginal_L[:,ii]) # get_mode
                m_vec[ii] = np.mean(marginal_M[:,ii]) # get_mode
                hpd_array_L[:,ii] = calcHPD(marginal_L[:,ii])
                hpd_array_M[:,ii] = calcHPD(marginal_M[:,ii])

                hpd_array_L50[:,ii] = calcHPD(marginal_L[:,ii],0.75)