def transition(dist, a, f, logspace=0):
    """
    Compute transition probabilities for a HMM. 
    
    to compute transition probabilities between hidden-states,
    when moving from time t to t+1,
    the genetic distance (cM) between the two markers are required.
    
    Assuming known parameters a and f.
    lf logspace = 1, calculations are log-transformed.
    
    Key in dictionary: 0 = not-IBD, 1 = IBD.
    """
    if logspace == 0:
        qk = exp(-a*dist)

        T = { # 0 = not-IBD, 1 = IBD
            1: {1: (1-qk)*f + qk, 0: (1-qk)*(1-f)},
            0: {1: (1-qk)*f, 0: (1-qk)*(1-f) + qk}
            }
        
    else:
        if dist == 0:
            dist = 1e-06

        ff = 1-f
        ad = a*dist
        A = expm1(ad)
        AA = -expm1(-ad)
        T = { # 0 = not-IBD, 1 = IBD
            1: {1: log1p(A*f)-ad, 0: log(AA*ff)},
            0: {1: log(AA*f), 0: log1p(A*ff)-ad}}
    return T

def f_fix(x, N_diploid, s):
    """
    The limits at s=0 and x=0 and x=1 are not implemented yet.
    Ultimately this function will be used as a term of an integrand
    whose integral will possibly be a well known special function.
    """
    return math.expm1(-4*N_diploid*s*x) / math.expm1(-4*N_diploid*s)

def get_fixation_probability_chen(p, s, h):
    """
    This uses the parameter conventions from Christina Chen et al. 2008.
    @param p: initial allele frequency
    @param s: positive when the mutant allele is fitter
    @param h: dominance
    @return: fixation probability
    """
    s_eff = 2.0 * s
    beta = 2.0 * h - 1.0
    if not s_eff:
        return p
    if not beta:
        return math.expm1(-s_eff*p) / math.expm1(-s_eff)
    alpha = h / beta
    if beta * s_eff > 0:
        # overdominant if 0 < alpha < 1
        f = erfi
    elif beta * s_eff < 0:
        # underdominant if 0 < alpha < 1
        f = special.erf
    else:
        raise ValueError
    L = math.sqrt(abs(beta * s_eff))
    a0 = f(L*(0 - alpha))
    a1 = f(L*(p - alpha))
    b0 = f(L*(0 - alpha))
    b1 = f(L*(1 - alpha))
    pfix = (a1 - a0) / (b1 - b0)
    return pfix

def get_pfix_transformed(p0, x4Ns, h):
    """
    Try to get the same result as the function in the kimura module.
    Change of variables according to eqn 3 of Chen et al.
    When I type
    (integral from 0 to p of exp(b*s*(x-a)**2) ) /
    (integral from 0 to 1 of exp(b*s*(x-a)**2) )
    I get
    ( erfi(sqrt(b)*sqrt(s)*a) - erfi(sqrt(b)*sqrt(s)*(a-p)) ) /
    ( erfi(sqrt(b)*sqrt(s)*a) - erfi(sqrt(b)*sqrt(s)*(a-1)) )
    @param p0: proportion of mutant alleles in the population
    @param x4Ns: 4Ns
    @return: fixation probability
    """
    if not x4Ns:
        # This is the neutral case.
        return p0
    if h == 0.5:
        # This is the genic case.
        # Checking for exact equality of 0.5 is OK.
        return math.expm1(-x4Ns*p0) / math.expm1(-x4Ns)
    b = 2.0 * h - 1.0
    a = h / (2.0 * h - 1.0)
    q = cmath.sqrt(b) * cmath.sqrt(x4Ns)
    #
    top = kimura.erfi(q*a) - kimura.erfi(q*(a-p0))
    bot = kimura.erfi(q*a) - kimura.erfi(q*(a-1))
    return top / bot

def J2_integrand(x, S):
    """
    This is part of equation (17) of Kimura and Ohta.
    """
    a = math.expm1(2*S*x)
    b = math.expm1(-2*S*x)
    c = x / (1 - x)
    return -(a * b) / c

	def SetCommand(self,linear=0,angular=0):
		# TODO: Implement expoenetial control and / or thresholded control strategy in here:
		if exponential_control == True:
			vel_linear  = expm1(abs(4*linear) ) / exponential_max
			vel_angular = expm1(abs(4*angular)) / exponential_max
			self.command.linear.x  = copysign(vel_linear , linear )
			self.command.angular.z = -1 * copysign(vel_angular, angular)
		else:
			self.command.linear.x  = linear
			self.command.angular.z = angular

def kimura_sojourn_helper(a, x):
    """
    Computes (exp(ax) - 1) / (exp(a) - 1) and accepts a=0.
    @param a: a scaled selection, can take any value
    @param x: a proportion between 0 and 1
    @return: a nonnegative value
    """
    if not a:
        return x
    else:
        return math.expm1(a*x) / math.expm1(a)

def get_pfix_approx(N_diploid, s):
    """
    This is an approximation of the fixation probability.
    It is derived by applying equation (3) in Kimura 1962,
    where p is the proportion of the preferred gene in the population,
    which in this case is 1/(2N) because it is a new mutant.
    It is given directly as equation (10).
    """
    N = N_diploid * 2
    if s:
        return math.expm1(-s) / math.expm1(-N*s)
    else:
        return 1.0 / N

def expm1(space, d):
    """ Returns exp(number) - 1, computed in a way that is
    accurate even when the value of number is close to zero"""
    try:
        return space.wrap(math.expm1(d))
    except OverflowError:
        return space.wrap(rfloat.INFINITY)

def alpha_merge_eqn(x, alpha, beta, x0, opthin=False):
    """Equation we need the root for to merge power law to modified
    blackbody

    Parameters
    ----------
    x : float
      h nu / k T to evaluate at

    alpha : float
      blue side power law index

    beta : float
      Dust attenuation power law index

    x0 : float
      h nu_0 / k T

    opthin : bool
      Assume optically thin case
    """

    try:
        # This can overflow badly
        xox0beta = (x / x0)**beta
        bterm = xox0beta / math.expm1(xox0beta)
    except OverflowError:
        # If xox0beta is very large, then the bterm is zero
        bterm = 0.0
    return x - (1.0 - math.exp(-x)) * (3.0 + alpha + beta * bterm)

    def __init__(self, T, fnorm, wavenorm=500.0):
        """Initializer

        Parameters:
        -----------
        T : float
          Temperature/(1+z) in K

        fnorm : float
          Normalization flux, in mJy

        wavenorm : float
          Wavelength of normalization flux, in microns (def: 500)
        """

        self._T = float(T)
        self._fnorm = float(fnorm)
        self._wavenorm = float(wavenorm)

        # Some constants -- eventually, replace these with
        # astropy.constants, but that is in development, so hardwire
        # for now
        self._hcokt = h * c / (k * self._T)
        self._xnorm = self._hcokt / self._wavenorm
        self._normfac = self._fnorm * math.expm1(self._xnorm) / \
            self._xnorm**3

def gen_modified_branch_history_sample(
        initial_state, final_state, blen_in, rates, P, t0=0.0):
    """
    This is a helper function for Nielsen modified rejection sampling.
    Yield (transition time, new state) pairs.
    The idea is to sample a path which may need to be rejected,
    and it is slightly clever in the sense that the path does not
    need to be rejected as often as do naive forward path samples.
    In more detail, this path sampler will generate paths
    conditional on at least one change occurring on the path,
    when appropriate.
    @param initial_state: initial state
    @param final_state: initial state
    @param blen_in: length of the branch
    @param rates: the rate away from each state
    @param P: transition matrix conditional on leaving a state
    @param t0: initial time
    """
    t = t0
    state = initial_state
    if state != final_state:
        rate = rates[initial_state]
        u = random.random()
        delta_t = -math.log1p(u*math.expm1(-blen_in*rate)) / rate
        t += delta_t
        if t >= blen_in:
            return
        distn = P[state]
        state = cmedbutil.random_category(distn)
        yield t, state
    for t, state in gen_branch_history_sample(state, blen_in, rates, P, t0=t):
        yield t, state

 def __init__(self, arms, rounds):
     self._arms = arms
     self._dist = np.ones(arms) / arms
     sim_learning_rate = min(1.0,
                             sqrt((arms * log(arms)) /
                                  (expm1(1.0) * rounds)))
     self._alpha = 3 * sim_learning_rate

def J1_integrand(x, S):
    """
    This is part of equation (17) of Kimura and Ohta.
    """
    a = math.expm1(2*S*x)
    b = math.exp(-2*S*x) - math.exp(-2*S)
    c = x * (1 - x)
    return (a * b) / c

 def fprime(curves):
      if curves == 0:
           return work_per_curve / (2*odds_factor_exists)
      arg = curves/median_curves
      l = -expm1(-arg) # 1-exp(-arg)
      # d/dx (x/l) = (l - x*l')/l^2
      # x*l' = arg*(1-l)
      return work_per_curve / odds_factor_exists * (l+arg*(l-1)) / (l*l)

    def calculateFalsePositiveProbabilityForNumberOfHashes(self, aNumberOfHashes):
        theInnerPower = -1.0 * (float( aNumberOfHashes ) * float(self.myNumWordsAdded)) / float(self.getBitArraySize())
        theInnerExpression = -1.0 * math.expm1(theInnerPower)

        theOuterPower = float( aNumberOfHashes )
        theOuterExpression = math.pow(theInnerExpression, theOuterPower)

        return theOuterExpression * 100.0

def log1mexp(a):
    """
            Computes log(1-exp(a)) according to Machler, "Accurately computing ..."
            Note: a should be a large negative value!
    """
    if a > 0: print >>sys.stderr, "# Warning, log1mexp with a=", a, " > 0"
    if a < -log(2.0): return log1p(-exp(a))
    else:             return log(-expm1(a))

 def __init__(self, options, rounds, highprob=False, lr=None):
     if lr == None:
         self._learning_rate = min(1.0, sqrt((options * log(options)) / (expm1(1) * rounds)))
     else:
         self._learning_rate = lr
     self._weights = np.ones(options, dtype=float)
     self.prob = self._weights / options
     self.gains = np.zeros(options)

def log1mexp_numba(a):
    if a > 0:
        print "LOGEXP"
        return

    if a < -log(2.0):
        return log1p(-exp(a))
    else:
        return log(-expm1(a))

def ciclo_none(rtot_species, period, Nindivs, invK, L_abs):
    rcal = rtot_species
    rspneq = calc_r_periodo_vh(rcal,invperiod) if (rcal>=0) else calc_r_periodo_vh(-rcal,invperiod)
    incNmalth= np.random.binomial(Nindivs,-math.expm1(-rspneq))
    if (rcal>0):
        inc_pop = Nindivs + incNmalth 
    else:
        inc_pop = Nindivs - incNmalth
    return([inc_pop,rcal])