# 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


import argparse, os,sys
import numpy as np
import random as rand
import platform, time
import multiprocessing, _thread
import multiprocessing.pool
import csv
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
		L0,Gl,M0,Gm = np.zeros(n_rates),np.zeros(n_rates),np.zeros(n_rates),np.zeros(n_rates)
		if len(Gl_index)==len(L0_index):
			for i in range(n_rates):
				L0[i] = t[j,L0_index[i]]
				Gl[i] = t[j,Gl_index[i]]
				M0[i] = t[j,M0_index[i]]
				Gm[i] = t[j,Gm_index[i]]
		else: # plot when model is equal_G
			for i in range(n_rates):
				L0[i] = t[j,L0_index[i]]
				Gl[i] = t[j,Gl_index[0]]
				M0[i] = t[j,M0_index[i]]
				Gm[i] = t[j,Gm_index[0]]
		
		Garray = np.array([Gl,Gm])
		age_vec,l_vec,m_vec = get_marginal_rates(args.m,L0,M0,Garray,Temp_at_events,shift_ind,root_age)
		marginal_L.append(l_vec)
		marginal_M.append(m_vec)
	
	marginal_L = np.array(marginal_L)
	marginal_M = np.array(marginal_M)
	
	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]))
	for i in range(np.shape(marginal_L)[1]):
		l_vec[i] = np.median(marginal_L[:,i])
		m_vec[i] = np.median(marginal_M[:,i])
		hpd_array_L[:,i] = calcHPD(marginal_L[:,i])

# 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: PyRateContinuous.py
import argparse, os,sys
import numpy as np
import random as rand
import platform, time
import multiprocessing, _thread
import multiprocessing.pool
import csv
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

		L0,Gl,M0,Gm = np.zeros(n_rates),np.zeros(n_rates),np.zeros(n_rates),np.zeros(n_rates)
		if len(Gl_index)==len(L0_index):
			for i in range(n_rates):
				L0[i] = t[j,L0_index[i]]
				Gl[i] = t[j,Gl_index[i]]
				M0[i] = t[j,M0_index[i]]
				Gm[i] = t[j,Gm_index[i]]
		else: # plot when model is equal_G
			for i in range(n_rates):
				L0[i] = t[j,L0_index[i]]
				Gl[i] = t[j,Gl_index[0]]
				M0[i] = t[j,M0_index[i]]
				Gm[i] = t[j,Gm_index[0]]
		
		Garray = np.array([Gl,Gm])
		age_vec,l_vec,m_vec = get_marginal_rates(args.m,L0,M0,Garray,Temp_at_events,shift_ind,root_age)
		marginal_L.append(l_vec)
		marginal_M.append(m_vec)
	
	marginal_L = np.array(marginal_L)
	marginal_M = np.array(marginal_M)
	
	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]))
	for i in range(np.shape(marginal_L)[1]):
		l_vec[i] = np.median(marginal_L[:,i])
		m_vec[i] = np.median(marginal_M[:,i])
		hpd_array_L[:,i] = calcHPD(marginal_L[:,i])

		hpd_array_M[:,i] = calcHPD(marginal_M[:,i])