def _calc_bic_(ll, num_snps, num_par, n):
    bic = -2 * (ll) + num_par * sp.log(n)
    extended_bic = bic + \
        2 * _log_choose_(num_snps, num_par - 2)
    modified_bic = bic + \
        2 * (num_par) * sp.log(num_snps / 2.2 - 1)
    return (bic, extended_bic, modified_bic)

def _Psat(Tdb):
    """
    ASHRAE Fundamentals Handbook pag 1.2 eq. 4
    input:
        Dry bulb temperature, K
    return:
        Saturation pressure, Pa
    """
    if 173.15 <= Tdb < 273.15:
        C1 = -5674.5359
        C2 = 6.3925247
        C3 = -0.009677843
        C4 = 0.00000062215701
        C5 = 2.0747825E-09
        C6 = -9.484024E-13
        C7 = 4.1635019
        pws = exp(C1/Tdb + C2 + C3*Tdb + C4*Tdb**2 + C5*Tdb**3 + C6*Tdb**4 +
                  C7*log(Tdb))
    elif 273.15 <= Tdb <= 473.15:
        C8 = -5800.2206
        C9 = 1.3914993
        C10 = -0.048640239
        C11 = 0.000041764768
        C12 = -0.000000014452093
        C13 = 6.5459673
        pws = exp(C8/Tdb + C9 + C10*Tdb + C11*Tdb**2 + C12*Tdb**3 + C13*log(Tdb))
    else:
        raise NotImplementedError("Incoming out of bound")

    return pws

def llfun(act, pred):
    epsilon = 1e-15
    pred = sp.maximum(epsilon, pred)
    pred = sp.minimum(1 - epsilon, pred)
    ll = sum(act * sp.log(pred) + sp.subtract(1, act) * sp.log(sp.subtract(1, pred)))
    ll = ll * -1.0 / len(act)
    return ll

def print_model_probability(logprob):
	"""
	Gives a nice overview of the model probability, allowing
	the practitioner to compare this model's probability to others
	"""
	prob = scipy.exp(logprob)

	limits = {
		'eq'         :logprob, 
		'barely'     :logprob - scipy.log(3),
		'substantial':logprob - scipy.log(10),
		'strong'     :logprob - scipy.log(30),
		'very strong':logprob - scipy.log(100)
	}
	for i in limits:
		limits[i] = "%5.1f" % limits[i]
		limits[i] = " " * (7 - len(limits[i])) + limits[i]
	
	print("Model probability ln(p(D|M, I)): [about 10^%.0f] %.5f" % (logprob / scipy.log(10), logprob))
	print(("""
	Table to compare support against other models (Jeffrey):

	   other model     |
	   ln(p(D|M,I))    | supporting evidence for this model
	 ------------------+-------------------------------------
		  >%%(eq)%s   Negative (supports other model)
	 %%(eq)%s ..%%(barely)%s   Barely worth mentioning
	 %%(barely)%s ..%%(substantial)%s   Substantial
	 %%(substantial)%s ..%%(strong)%s   Strong
	 %%(strong)%s ..%%(very strong)%s   Very strong
		  <%%(very strong)%s   Decisive

	be careful.
	""" % tuple(['s']*10)) % limits)

def summarize_splits(splits, weighted=True):
    rows = []
    if weighted:
        v = [ (x.likelihood * x.weight, x) for x in splits ]
    else:
        v = [ (x.likelihood, x) for x in splits ]
    ptot = sum([ x[0] for x in v ])
    v.sort(); v.reverse()
    opt = scipy.log(v[0][0])

    rows.append(["split", "lnL", "Rel.Prob"])
    sumprob = 0.0
    for L, split in v:
        lnL = scipy.log(L)
        relprob = (L/ptot)
        if sumprob < 0.95:
        #if (opt - lnL) < 2:
            rows.append([str(split), "%.4g" % lnL, "%.4g" % relprob])
        sumprob += relprob
    widths = []
    for i in range(3):
        w = max([ len(x[i]) for x in rows ])
        for x in rows:
            x[i] = x[i].ljust(w)
    return [ "  ".join(x) for x in rows ]

def compute_continuous_prob_value(parameters, distribution, rvs):
    mean = float(parameters[0])
    stddev = float(parameters[1]) / 100 * mean
    A = float(parameters[2])
    B = float(parameters[3])
    result = float("-inf")

    if rvs is None:
        while result <= A or result > B:

            if distribution == "normal":
                rvs = stats.norm.rvs
                result = rvs(mean, stddev)

            elif distribution == "lognormal":
                variance = stddev ** 2.0
                mu = log(mean ** 2.0 / sqrt(variance + mean ** 2.0))
                sigma = sqrt(log((variance / mean ** 2.0) + 1.0))
                rvs = stats.lognorm.rvs
                result = rvs(sigma, scale=scipy.exp(mu))

            elif distribution == "gamma":
                betha = (stddev) ** 2 / mean
                alpha = mean / betha
                rvs = stats.gamma.rvs
                result = rvs(alpha, scale=betha)
    else:
        result = 1

    return result

def binary_logloss(p, y):
    epsilon = 1e-15
    p = sp.maximum(epsilon, p)
    p = sp.minimum(1-epsilon, p)
    res = sum(y * sp.log(p) + sp.subtract(1, y) * sp.log(sp.subtract(1, p)))
    res *= -1.0/len(y)
    return res

def KramersKronigFFT(ImX_A):
	'''	Hilbert transform used to calculate real part of a function from its imaginary part
	uses piecewise cubic interpolated integral kernel of the Hilbert transform
	use only if len(ImX_A)=2**m-1, uses fft from scipy.fftpack  '''
	X_A = sp.copy(ImX_A)
	N = int(len(X_A))
	## be careful with the data type, orherwise it fails for large N
	if N > 3e6: A = sp.arange(3,N+1,dtype='float64')
	else:       A = sp.arange(3,N+1)  
	X1 = 4.0*sp.log(1.5)
	X2 = 10.0*sp.log(4.0/3.0)-6.0*sp.log(1.5)
	## filling the kernel
	if N > 3e6: Kernel_A = sp.zeros(N-2,dtype='float64')
	else:       Kernel_A = sp.zeros(N-2)
	Kernel_A = (1-A**2)*((A-2)*sp.arctanh(1.0/(1-2*A))+(A+2)*sp.arctanh(1.0/(1+2*A)))\
	+((A**3-6*A**2+11*A-6)*sp.arctanh(1.0/(3-2*A))+(A+3)*(A**2+3*A+2)*sp.arctanh(1.0/(2*A+3)))/3.0
	Kernel_A = sp.concatenate([-sp.flipud(Kernel_A),sp.array([-X2,-X1,0.0,X1,X2]),Kernel_A])/sp.pi
	## zero-padding the functions for fft
	ImXExt_A = sp.concatenate([X_A[int((N-1)/2):],sp.zeros(N+2),X_A[:int((N-1)/2)]])
	KernelExt_A = sp.concatenate([Kernel_A[N:],sp.zeros(1),Kernel_A[:N]])
	## performing the fft
	ftReXExt_A = -fft(ImXExt_A)*fft(KernelExt_A)
	ReXExt_A = sp.real(ifft(ftReXExt_A))
	ReX_A = sp.concatenate([ReXExt_A[int((3*N+3)/2+1):],ReXExt_A[:int((N-1)/2+1)]])
	return ReX_A

    def _LML_covar(self,hyperparams,debugging=False):
        """
        log marginal likelihood
        """
        try:
            KV = self.get_covariances(hyperparams,debugging=debugging)
        except LA.LinAlgError:
            LG.error('linalg exception in _LML_covar')
            return 1E6
        except ValueError:
            LG.error('value error in _LML_covar')
            return 1E6
 
        lml_quad = 0.5*(KV['Ytilde']*KV['UYU']).sum()
        lml_det =  0.5 * SP.log(KV['S']).sum()
        lml_const = 0.5*self.n*self.t*(SP.log(2*SP.pi))
        
        if debugging:
            # do calculation without kronecker tricks and compare
            _lml_quad = 0.5 * (KV['alpha']*KV['Yvec']).sum()
            _lml_det =  SP.log(SP.diag(KV['L'])).sum()
            assert SP.allclose(_lml_quad,lml_quad),  'ouch, quadratic form is wrong in _LMLcovar'
            assert SP.allclose(_lml_det, lml_det), 'ouch, ldet is wrong in _LML_covar'
        
        lml = lml_quad + lml_det + lml_const

        return lml

    def _visco0(self, rho, T, fase=None, coef=False):

        Visco0 = lambda T: -0.135311743/log(T) + 1.00347841 + \
            1.20654649*log(T) - 0.149564551*log(T)**2 + 0.0125208416*log(T)**3

        def ViscoE(T, rho):
            x = log(T)
            B = -47.5295259/x+87.6799309-42.0741589*x+8.33128289*x**2-0.589252385*x**3
            C = 547.309267/x-904.870586+431.404928*x-81.4504854*x**2+5.37005733*x**3
            D = -1684.39324/x+3331.08630-1632.19172*x+308.804413*x**2-20.2936367*x**3
            return rho.gcc*B+rho.gcc**2*C+rho.gcc**3*D


        if T < 100:
            # Section 4.2.1 for 3.5 < T < 100
            no = Visco0(T)
            ne = ViscoE(T, rho)
            n = exp(no+ne)
        else:
            # Section 4.2.1 for T > 100
            no = 196*T**0.71938*exp(12.451/T-295.67/T**2-4.1249)
            ne = exp(Visco0(T)+ViscoE(T, rho))-exp(Visco0(T)+ViscoE(T, unidades.Density(0)))
            n = no+ne

        if coef:
            return ne
        else:
            return unidades.Viscosity(n*1e-6, "P")

    def _LML_covar(self, hyperparams):
        #calculate marginal likelihood of kronecker GP

        #1. get covariance structures needed:
        try:
            KV = self.get_covariances(hyperparams)
        except linalg.LinAlgError:
            LG.error("exception caught (%s)" % (str(hyperparams)))
            return 1E6
        #2. build lml 
        LML = 0
        LMLc = 0.5* self.nd * SP.log(2.0 * SP.pi)
        #constant part of negative lml
        #quadratic form
        Si = KV['Si']
       
        LMLq = 0.5 * SP.dot(KV['y_rot'].ravel(),KV['YSi'].ravel() )
        #determinant stuff
        LMLd = -0.5 * SP.log(Si).sum()

        if VERBOSE:
            print "costly verbose debugging on"
            K = SP.kron(KV['Kr'],KV['Kc']) + SP.diag(KV['Knoise'])
            Ki = SP.linalg.inv(K)
            LMLq_ = 0.5* SP.dot(SP.dot(self.y.ravel(),Ki),self.y.ravel())
            LMLd_ = 0.5* 2 * SP.log(SP.linalg.cholesky(K).diagonal()).sum()
            check_dist(LMLq,LMLq_)
            check_dist(LMLd,LMLd_)
            

        return LMLc+LMLq+LMLd

def NTU_fPR(P, R, flujo, **kwargs):
    """Calculo de la factor de correccion
    Flujo vendra definido por su acronimo
        CF: Counter flow
        PF: Parallel flow
        CrFMix: Crossflow, both fluids mixed
        CrFSMix: Crossflow, one fluid mixed, other unmixed
        CrFunMix: Crossflow, both fluids unmixed
        1-2TEMAE: 1-2 pass shell and tube exchanger

    kwargs: Opciones adicionales:
        mixed: corriente mezclada para CrFSMix
            Cmin, Cmax
    """

    if flujo == "1-2TEMAE":
        if R == 1:
            NTU = log((1-P)/2-3*P)
        else:
            E = (1+R**2)**0.5
            NTU = log((2-P*(1+R-E))/(2-P*(1+R+E)))/E

    else:
        if R == 1:
            NTU = P/(1-P)
        else:
            NTU = log((1-R/P)/(1-P))/(1-R)

    return NTU

def tfidf(termFrequency):
	""" The student must code this. """
	gf = sp.sum(termFrequency,axis=1).astype(float)
	p = (termFrequency.T/gf).T
	g = sp.sum(p*sp.log(p+1)/sp.log(len(p[0,:])),axis=1) + 1
	a = (sp.log(termFrequency + 1).T*g).T
	return a

def cdi_info(energy, h, z, pix, del_x_d, verbose = False):
    """
    h - object size\nz - sam-det dist\npix - # of pix\ndel_x_d - pixel size
    """
    x = (pix/2.)*del_x_d
    l = energy_to_wavelength(energy)
    NF = lambda nh, nl, nz : nh**2./(nl*nz)
    del_x_s = lambda l, z, x : (l*z)/(2.*x)
    nNF = NF(h,l,z)
    OS = lambda l,z,x,h,pix : ((pix*del_x_s(l,z,x))**2.)/(h**2.)
    nOS = OS(l,z,x,h,pix)
    if verbose:
        pyl.figure()
        zrange = sp.linspace(0, 2*z, 100)
        pyl.plot(zrange, sp.log(NF(h,l,zrange)))
        pyl.title('NF')
        pyl.xlabel('z [m]')
        pyl.ylabel('log NF')
        pyl.figure()
        pyl.plot(zrange, sp.log(OS(l,zrange, x, h, pix)))
        pyl.title('OS')
        pyl.xlabel('z [m]')
        pyl.ylabel('log OS')
    
    print 'NF: %1.2e\nOS: %1.2e\ndel_x_d: %1.2e\nw_d: %1.2e\ndel_x_s: %1.2e\nw_s: %1.2e' % (nNF, nOS, del_x_d, pix*del_x_d, del_x_s(l,z,x), del_x_s(l,z,x)*pix)
    aperture_stats(energy, z, x)

    def loglike(self, data, paravec, sign = 1): 

        lu, lsig = paravec
        loglike = (-1/(2*np.exp(lsig)**2)) * sum((sp.log(data)-np.exp(lu))**2) - (len(data)/2) * sp.log(2*sp.pi) \
        - len(data) * sp.log(np.exp(lsig)) - sum(sp.log(data))
        loglike = sign*loglike
        return loglike

    def TB_Cv_exceso(self, T, P):
        """Método de cálculo de la capacidad calorífica a volumen constante de exceso mediante la ecuación de estado de Trebble-Bishnoi"""
        a, b, c, d, q1, q2=self.TB_lib(T, P)
        v=self.TB_V(T, P)
        z=P*v/R_atml/T
        t=1+6*c/b+c**2/b**2+4*d**2/b**2
        tita=abs(t)**0.5
        A=a*P/R_atml**2/T**2
        B=b*P/R_atml/T
        u=1+c/b
        delta=v**2+(b+c)*v-b*c-d**2
        beta=1.+q2*(1-self.tr(T)+log(self.tr(T)))
        da=-q1*a/self.Tc
        dda=q1**2*a/self.Tc**2
        if self.tr(T)<=1.0:
            db=b/beta*(1/T-1/self.Tc)
            ddb=-q2*b/beta/T**2
        else:
            db=0
            ddb=0

        dt=-db/b**2*(6*c+2*c**2/b+8*d**2/b)
        dtita=abs(dt)/20
        if t>=0:
            lamda=log((2*z+B*(u-tita))/(2*z+B*(u+tita)))
            dlamda=(db-db*tita-b*dtita)/(2*v+b+c-b*tita)-(db+db*tita+b*dtita)/((2*v+b+c+b*tita))
        else:
            lamda=2*arctan((2*z+u*B)/B/tita)-pi
            dlamda=2/(1+((2*v+b+c)/b/tita)**2)*(db/b/tita-(2*v+b+c)*(db/b**2/tita+dtita/b/tita**2))

        Cv=1/b/tita*(dlamda*(a-da*T)-lamda*dda*T-lamda*(a-da*T)*(db/b+dtita/tita))+(ddb*T+db)*(-R_atml*T/(v-b)+a/b**2/t*((v*(3*c+b)-b*c+c**2-2*d**2)/delta+(3*c+b)*lamda/b/tita))+db*T*(-R_atml/(v-b)-R_atml*T*db/(v-b)**2+1/b**2/t*(da-2*a*db/b-a*dt/t)*((v*(3*c+b)-b*c+c**2-2*d**2)/delta+(3*c+b)*lamda/b/tita)+a/b**2/t*(db*(v-c)*(v**2-2*c*v-c**2+d**2)/delta**2+db*lamda/b/tita+(3*c+b)/b/tita*(dlamda-lamda*(db/b+dtita/tita))))
        return unidades.SpecificHeat(Cv*101325/1000/self.peso_molecular, "JkgK")

def log_d_pois_like_trunc_5(d,s1,s2,a,p):
    """double poisson w max 5 goals"""
    #dp = np.sign(d)*np.power(np.abs(d),p)
    dp = 1.5*np.arctan(d)    #print(dp)
    return ( log(a)*(s1+s2)+dp*(s1-s2) - 2*a*cosh(dp)
         -gammaln(s1+1) - gammaln(s2+1) 
        -log(gammaincc(6,a*exp(-dp))*gammaincc(6,a*exp(dp)) ) ) 

    def LML(self,params=None,*kw_args):
        """
        calculate LML
        """
        if params is not None:
            self.setParams(params)

        self._update_cache()
        
        start = TIME.time()

        #1. const term
        lml  = self.N*self.P*SP.log(2*SP.pi)

        #2. logdet term
        lml += SP.sum(SP.log(self.cache['Sc2']))*self.N
        lml += 2*SP.log(SP.diag(self.cache['cholB'])).sum()

        #3. quatratic term
        lml += SP.sum(self.cache['LY']*self.cache['LY'])
        lml -= SP.sum(self.cache['WLY']*self.cache['BiWLY'])

        lml *= 0.5

        smartSum(self.time,'lml',TIME.time()-start)
        smartSum(self.count,'lml',1)

        return lml

def evaluate_ll(y, yhat):
    epsilon = 1e-15
    yhat = sp.maximum(epsilon, yhat)
    yhat = sp.minimum(1-epsilon, yhat)
    ll = sum(y*sp.log(yhat) + sp.subtract(1,y)*sp.log(sp.subtract(1,yhat)))
    ll = ll * -1.0/len(y)
    return ll

 def fit(self, kk=None):
     """
     Fit Fourier spectrum with the function set at class instantination
     ==> NB: fitting is done in logarithmic coordinates
     and fills plotting arrays with data
     --------
     Options:
     --------
     kk
        (k1,k2) <None> spectral interval for function fitting
        by default interval [ kk[1], kk[imax__kk] ] will be fitted
        ==> i.e. k=0 is excluded
     """
     # fitting interval
     if kk:
         ik_min=(self.fft_data.kk[1:self.fft_data.imax__kk]<=kk[0]).nonzero()[0][-1]
         ik_max=(self.fft_data.kk[1:self.fft_data.imax__kk]<=kk[1]).nonzero()[0][-1]
     else:
         ik_min=1;
         ik_max=self.fft_data.imax__kk
     # do fitting
     self.__popt,self.__pcov = scipy.optimize.curve_fit(self.__func_fit,
                                                        scipy.log(self.fft_data.kk[ik_min:ik_max]),
                                                        scipy.log(self.fft_data.Ik[ik_min:ik_max]) )
     # boundaries of fitted interval
     self.kmin = self.fft_data.kk[ik_min]
     self.kmax = self.fft_data.kk[ik_max]
     # fill plot arrays <===============
     self.kk_plot=scipy.logspace( scipy.log10(self.kmin),
                                  scipy.log10(self.kmax),
                                  self.nk_plot )
     self.Ik_plot=self.fitting_function(self.kk_plot)