def genenergies(fnR,fnQ,seqsR,seqsQ,gamma,sQ,sR,R0): #Parses seqs and model type then calculates and returns energies R is transcription factor, Q is RNAP
    ematR = np.genfromtxt(fnR,skiprows=1)
    ematQ = np.genfromtxt(fnQ,skiprows=1)
    fR = open(fnR)
    fQ = open(fnQ)
    mattype = fR.read()[:6] #mattype must be the same
    #mattypeQ = fQ.read()[:6]
    energies = np.zeros(len(seqsQ))
    N = len(seqsQ)
    mut_region_lengthQ = len(seqsQ[0])
    mut_region_lengthR = len(seqsR[0])
    
    if mattype == '1Point':
            for i,s in enumerate(seqsR):
                seq_matR = seq2mat(s)
		seq_matQ = seq2mat(seqsQ[i])
		RNAP = (seq_matQ*ematQ).sum()*sQ
		TF = (seq_matR*ematR).sum()*sR + R0
                energies[i] = -RNAP + mp.log(1 + mp.exp(-TF - gamma)) - mp.log(1 + mp.exp(-TF))
    '''
    elif mattype == '2Point':
            for i,s in enumerate(seqs):
                seq_mat = np.zeros(round(sp.misc.comb(mut_region_length,2))*16)
                seq_mat[seq2mat2(s)] = 1
                energies[i] = (seq_mat*(emat.ravel())).sum()
    elif mattype == '3Point':
            for i,s in enumerate(seqs):
                seq_mat = np.zeros(round(sp.misc.comb(mut_region_length,3))*64)
                seq_mat[seq2mat3(s)] = 1
                energies[i] = (seq_mat*(emat.ravel())).sum()
    '''
    return energies

def mean_from_output(output):
    """Given output of the form [params,ll], return mean parameter
    vector"""
    param_vectors,lls = transpose(output)
    z = sum(mpmath.exp(ll) for ll in lls)
    return map(sum,transpose([[(mpmath.exp(ll)/z)*p for p in param_vector]
                               for param_vector,ll in output]))

 def integrate( self, a, b ):
     halfk = self.halfforceconstant
     x0 = self.reference
     beta = self.beta
     avalue = exp( -beta*halfk*(a-x0)**2 )
     bvalue = exp( -beta*halfk*(b-x0)**2 )
     return 0.5*(b-a)*( avalue + bvalue ) 

def fermihalf(x, sgn):
    """ Series approximation to the F_{1/2}(x) or F_{-1/2}(x) 
        Fermi-Dirac integral """

    f = lambda k: mp.sqrt(x ** 2 + np.pi ** 2 * (2 * k - 1) ** 2)

    # if x < -100:
    #    return 0.0
    if x < -9 or True:
        if sgn > 0:
            return mp.exp(x)
        else:
            return mp.exp(x)

    if sgn > 0:  # F_{1/2}(x)
        a = np.array((1.0 / 770751818298, -1.0 / 3574503105, -13.0 / 184757992,
                      85.0 / 3603084, 3923.0 / 220484, 74141.0 / 8289, -5990294.0 / 7995))
        g = lambda k: mp.sqrt(f(k) - x)

    else:  # F_{-1/2}(x)
        a = np.array((-1.0 / 128458636383, -1.0 / 714900621, -1.0 / 3553038,
                      27.0 / 381503, 3923.0 / 110242, 8220.0 / 919))
        g = lambda k: -0.5 * mp.sqrt(f(k) - x) / f(k)

    F = np.polyval(a, x) + 2 * np.sqrt(2 * np.pi) * sum(map(g, range(1, 21)))
    return F

def dedekind(tau, floatpre):
    """
    Algorithm 22 (Dedekind eta)
    Input : tau in the upper half-plane, k in N
    Output : eta(tau)
    """
    a = 2 * mpmath.pi / mpmath.mpf(24)
    b = mpmath.exp(mpmath.mpc(0, a))
    p = 1
    m = 0
    while m <= 0.999:
        n = nearest_integer(tau.real)
        if n != 0:
            tau -= n
            p *= b ** n
        m = tau.real * tau.real + tau.imag * tau.imag
        if m <= 0.999:
            ro = mpmath.sqrt(mpmath.power(tau, -1) * 1j)
            if ro.real < 0:
                ro = -ro
            p = p * ro
            tau = (-p.real + p.imag * 1j) / m
    q1 = mpmath.exp(a * tau * 1j)
    q = q1 ** 24
    s = 1
    qs = mpmath.mpc(1, 0)
    qn = 1
    des = mpmath.mpf(10) ** (-floatpre)
    while abs(qs) > des:
        t = -q * qn * qn * qs
        qn = qn * q
        qs = qn * t
        s += t + qs
    return p * q1 * s

	def sigma(self,z):
		# A+S 18.10.
		from mpmath import pi, jtheta, exp, mpc, sqrt, sin
		Delta = self.Delta
		e1, _, _ = self.__roots
		om = self.__periods[0] / 2
		omp = self.__periods[1] / 2
		if self.__ng3:
			z = mpc(0,1) * z
		if Delta > 0:
			tau = omp / om
			q = (exp(mpc(0,1) * pi() * tau)).real
			eta = -(pi()**2 * jtheta(n=1,z=0,q=q,derivative=3)) / (12 * om * jtheta(n=1,z=0,q=q,derivative=1))
			v = (pi() * z) / (2 * om)
			retval = (2 * om) / pi() * exp((eta * z**2)/(2 * om)) * jtheta(n=1,z=v,q=q)/jtheta(n=1,z=0,q=q,derivative=1)
		elif Delta < 0:
			om2 = om + omp
			om2p = omp - om
			tau2 = om2p / (2 * om2)
			q = mpc(0,(mpc(0,1) * exp(mpc(0,1) * pi() * tau2)).imag)
			eta2 = -(pi()**2 * jtheta(n=1,z=0,q=q,derivative=3)) / (12 * om2 * jtheta(n=1,z=0,q=q,derivative=1))
			v = (pi() * z) / (2 * om2)
			retval = (2 * om2) / pi() * exp((eta2 * z**2)/(2 * om2)) * jtheta(n=1,z=v,q=q)/jtheta(n=1,z=0,q=q,derivative=1)
		else:
			g2, g3 = self.__invariants
			if g2 == 0 and g3 == 0:
				retval = z
			else:
				c = e1 / 2
				A = sqrt(3 * c)
				retval = (1 / A) * sin(A*z) * exp((c*z**2) / 2)
		if self.__ng3:
			return mpc(0,-1) * retval
		else:
			return retval

def test_stoch_eig_high_prec():
    n = 1e-100
    with mp.workdps(100):
        P = mp.matrix([[1-3*(mp.exp(n)-1), 3*(mp.exp(n)-1)],
                       [mp.exp(n)-1      , 1-(mp.exp(n)-1)]])

    run_stoch_eig(P, verbose=VERBOSE)

def test_gth_solve_high_prec():
    n = 1e-100
    with mp.workdps(100):
        P = mp.matrix([[-3*(mp.exp(n)-1), 3*(mp.exp(n)-1)],
                       [mp.exp(n)-1     , -(mp.exp(n)-1) ]])

    run_gth_solve(P, verbose=VERBOSE)

def dfdy (y,x,b,c):

    global FT,XO,XT,XF

    ft=FT

    v=x[0]
    i=y[0]

    iss=IS=b[0]
    n=N=b[1]
    ikf=IKF=b[2]
    isr=ISR=b[3]
    nr=NR=b[4]
    vj=VJ=b[5]
    m=M=b[6]
    rs=RS=b[7]

    #fh = iss**mpm.mpf(2)*rs*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - mpm.mpf(1))))*(mpm.exp((-i*rs + v)/(ft*n)) - mpm.mpf(1))*mpm.exp((-i*rs + v)/(ft*n))/(mpm.mpf(2)*ft*n*(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - mpm.mpf(1)))) - iss*rs*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - mpm.mpf(1))))*mpm.exp((-i*rs + v)/(ft*n))/(ft*n)
    #sh = isr*m*rs*(mpm.mpf(1) - (-i*rs + v)/vj)*((mpm.mpf(1) - (-i*rs + v)/vj)**mpm.mpf(2) + mpm.mpf('0.005'))**(m/mpm.mpf(2))*(mpm.exp((-i*rs + v)/(ft*nr)) - mpm.mpf(1))/(vj*((mpm.mpf(1) - (-i*rs + v)/vj)**mpm.mpf(2) + mpm.mpf('0.005'))) - isr*rs*((mpm.mpf(1) - (-i*rs + v)/vj)**mpm.mpf(2) + mpm.mpf('0.005'))**(m/mpm.mpf(2))*mpm.exp((-i*rs + v)/(ft*nr))/(ft*nr)

    fh = iss**XT*rs*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*(mpm.exp((-i*rs + v)/(ft*n)) - XO)*mpm.exp((-i*rs + v)/(ft*n))/(XT*ft*n*(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO))) - iss*rs*mpm.sqrt(ikf/(ikf + iss*(mpm.exp((-i*rs + v)/(ft*n)) - XO)))*mpm.exp((-i*rs + v)/(ft*n))/(ft*n)
    sh = isr*m*rs*(XO - (-i*rs + v)/vj)*((XO - (-i*rs + v)/vj)**XT + XF)**(m/XT)*(mpm.exp((-i*rs + v)/(ft*nr)) - XO)/(vj*((XO - (-i*rs + v)/vj)**XT + XF)) - isr*rs*((XO - (-i*rs + v)/vj)**XT + XF)**(m/XT)*mpm.exp((-i*rs + v)/(ft*nr))/(ft*nr)

    return mpm.matrix ([[fh+sh]])

 def bimax_integrand(self, z, wc, l, n, t):
     """
     Integrand of electron-noise integral.
     
     """
     return f1(wc*l/z/mp.sqrt(2)) * z * \
         (mp.exp(-z**2) + n/mp.sqrt(t)*mp.exp(-z**2 / t)) / \
         (mp.fabs(BiMax.d_l(z, wc, n, t))**2 * wc**2)

def temp_t(t, x=mpf(1), Q=mpf(1), A=mpf(1), Ti=mpf(10)):
    u'''
    Definição da derivada da temperatura em relação ao tempo, temp_t = T_t(t, x)
    '''
    t = mpf(t)
    termo1 = sqrt(A/(t*pi)) * exp(-x**2/(t*4*A))
    termo2 = sqrt(A/(t*pi)) * (x**2/(t**2 * 2 * A)) * exp(-x**2/(t*4*A))
    termo3 = (x**2/(sqrt(t**3 * A) * 4)) * erfc_z(x / (mpf(2) * sqrt(t * A)))
    
    return Q * (termo1 + termo2 + termo3)

def fl2(x0, seta2, seff2):

    coeff = 1./np.sqrt(np.pi)
    xm = x0/np.sqrt(2.*seta2)
    tau = lamda1*seff2/np.sqrt(2.*seta2)

    intgr1 = -coeff*(tau-xm)*exp(-(tau-xm)**2) +(.5+(tau-xm)**2)*erfc(tau-xm)
    intgr2 = -coeff*(tau+xm)*exp(-(tau+xm)**2) +(.5+(tau+xm)**2)*erfc(tau+xm)

    return (seta2/(1.+lamda2*seff2)**2)*(intgr1+intgr2)  

def sph_i2n_exact(n, z):
    """Return the value of i^{(2)}_n computed using the exact formula.

    The expression used is http://dlmf.nist.gov/10.49.E10 .

    """
    zm = mpmathify(z)
    s1 = sum(mpc(-1,0)**k * _a(k, n)/zm**(k+1) for k in xrange(n+1))
    s2 = sum(_a(k, n)/zm**(k+1) for k in xrange(n+1))
    return exp(zm)/2 * s1 + mpc(-1,0)**n*exp(-zm)/2 * s2

def DedekindEtaA4(tau):
    ''' Compute the derivative of the Dedekind Eta function for imaginary argument tau. 
        Numerically. '''
    try:
        import mpmath as mp
        mpmath_loaded = True
    except ImportError:
        mpmath_loaded = False 
    
    return mp.cbrt(0.5*mp.jtheta(2,0,mp.exp(-mp.pi*tau))*mp.jtheta(3,0,mp.exp(-mp.pi*tau))*mp.jtheta(4,0,mp.exp(-mp.pi*tau)))

def DedekindEtaA2(tau):
    ''' Compute the derivative of the Dedekind Eta function for imaginary argument tau. 
        Numerically. '''
    try:
        import mpmath as mp
        mpmath_loaded = True
    except ImportError:
        mpmath_loaded = False 
    
    return mp.exp(-mp.pi/12.0)*mp.jtheta(3,mp.pi*(mp.j*tau+1.0)/2.0,mp.exp(-3.0*mp.pi))

def fl1(x0, seta2, seff2):

    coeff = 1./np.sqrt(np.pi)
    xm = x0/np.sqrt(2.*seta2)
    tau = lamda1*seff2/np.sqrt(2.*seta2)

    intgr1 = coeff*exp(-(tau-xm)**2) -(tau-xm)*erfc(tau-xm)
    intgr2 = coeff*exp(-(tau+xm)**2) -(tau+xm)*erfc(tau+xm)

    return .5*(np.sqrt(2.*seta2)/(1.+lamda2*seff2))*(intgr1+intgr2)

def sum_gaussian_theta(variance, offset, interval):
    factor = mpmath.sqrt(1.0/(variance*tau))
    factor_exponent=mpmath.mpf(-1.0)/(2*variance)
    exponent_factor = mpmath.mpf(offset)*offset
    factor_full = factor*mpmath.exp(exponent_factor*factor_exponent)
    exponent_interval = interval*interval
    exponent_offset = 2*interval*offset
    q = mpmath.exp(factor_exponent*exponent_interval)
    z = factor_exponent*exponent_offset/(2*mpmath.j)
    theta = mpmath.jtheta(3,z,q)
    return factor_full*theta

 def gaussian_total(self,offset):
     factor = mpmath.sqrt(1.0/(self._period*self._tau))
     factor_exponent = -1.0/(2*self._period)
     exponent_factor = mpmath.mpf(offset)*offset
     exponent_interval = self._tau*self._tau
     exponent_offset = 2*self._tau*offset
     factor_full = factor*mpmath.exp(exponent_factor*factor_exponent)
     q = mpmath.exp(factor_exponent*exponent_interval)
     z = factor_exponent*exponent_offset/(2*mpmath.j)
     theta = mpmath.jtheta(3,z,q).real
     return factor_full*theta

def fq(x0, seta2, seff2):

    coeff = 1./np.sqrt(np.pi)
    xm = x0/np.sqrt(2.*seta2)
    tau = lamda1*seff2/np.sqrt(2.*seta2)
    psi = lamda2*seff2*x0/np.sqrt(2.*seta2)
    
    intgr1 = coeff*(-tau+2*psi+xm)*exp(-(tau+xm)**2) +.5*(1.+2.*(tau-psi)**2)*erfc(tau+xm)
    intgr2 = coeff*(-tau-2*psi-xm)*exp(-(tau-xm)**2) +.5*(1.+2.*(tau+psi)**2)*erfc(tau-xm)
    intgr3 = .5*x0**2*(erf(tau+xm)+erf(tau-xm))

    return (seta2/(1.+lamda2*seff2)**2)*(intgr1+intgr2)+intgr3

    def energy(self, clustering):
        energy = mpmath.mpf(0.0)
        new_vertex_distributions = _combine_vertex_distributions_given_clustering(
            self.vertex_distributions, clustering)

        # likelihood
        likelihood_energy = -self._log_likelihood(clustering, new_vertex_distributions)

        # prior on similarity:
        # We prefer the cluster whose minimum similarity is large.
        # - the similarity of a pair of vertexes is measured by the similarity
        #   of top 10 words in the distribution. (measure each word type
        #   respectively and take average)
        intra_cluster_energy = mpmath.mpf(0.0)
        for cluster_id, cluster_vertex_set in enumerate(clustering):
            min_similarity_within_cluster = self._min_similarity_within_cluster(cluster_vertex_set, new_vertex_distributions[cluster_id])
            intra_cluster_energy += -mpmath.log(mpmath.exp(min_similarity_within_cluster - 1))

        # Between cluster similarity:
        #  - For each pair of clusters, we want to find the pair of words with maximum similarity
        #    and prefer this similarity value to be small.
        inter_cluster_energy = mpmath.mpf(0.0)
        if len(clustering) > 1:
            for i in range(0, len(clustering)-1):
                for j in range(i+1, len(clustering)):
                    max_similarity_between_clusters = self._max_similarity_between_clusters(clustering[i], clustering[j])
                    inter_cluster_energy += -mpmath.log(mpmath.exp(-max_similarity_between_clusters))

        # prior on clustering complexity: prefer small number of clusters.
        length_energy = -mpmath.log(mpmath.exp(-len(clustering)))

        # classification: prefer small number of categories.
        class_energy = 0.0
        if self._classifier is not None:
            num_classes = self._calculate_num_of_categories(clustering, new_vertex_distributions)
            class_energy = -mpmath.log(mpmath.exp(-(abs(num_classes-len(clustering)))))

        # classification confidence: maximize the classification confidence
        confidence_energy = 0.0
        for cluster_id, cluster_vertex_set in enumerate(clustering):
            (category, confidence) = self._predict_label(new_vertex_distributions[cluster_id])
            confidence_energy += -mpmath.log(confidence)

        energy += (0.5)*likelihood_energy + intra_cluster_energy + inter_cluster_energy + 30.0*length_energy + 20.0*class_energy + confidence_energy
        logging.debug('ENERGY: {0:12.6f}\t{1:12.6f}\t{2:12.6f}\t{3:12.6f}\t{4:12.6f}\t{5:12.6f}'.format(
            likelihood_energy.__float__(),
            intra_cluster_energy.__float__(),
            inter_cluster_energy.__float__(),
            length_energy.__float__(),
            class_energy.__float__(),
            confidence_energy.__float__()))
        return energy