def raise_res(T, W, c, mu=0, sigma=1):
    """Increase the resolution of a wiener series by a factor of c.
        
        Returns a more reolved Wiener series and its associate time series
        
        T = the given Time series.
        W = the associated Wiener series.
        c = Scaling factor (integer greater than 1).
        mu = Mean of W's underlying normal distribution.
        sigma = Standard deviation of W's underlying normal distribution.
    """
    dT = T[1] - T[0]
    dt = float(T[1] - T[0]) / c
    t_series = []
    w_series = []
    for i in range(len(T) - 1):
        t = T[i]
        w_t = W[i]
        t_next = T[i + 1]
        w_next = W[i + 1]
        t_series.append(t)
        w_series.append(w_t)
        for j in range(c - 1):
            t += dt
            dW = w_next - w_t
            drawfrm_cum = np.sqrt(2) * np.sqrt(t_next - t) * sigma * erfc(random())
            if np.sqrt(2) * np.sqrt(t_next - t) * sigma * erfc(-2 * random()) < abs(dW):
                w_t += abs(gauss(0, np.sqrt(dt) * sigma)) * float(dW) / abs(dW)
            else:
                w_t += gauss(0, np.sqrt(dt) * sigma)
            t_series.append(t)
            w_series.append(w_t)
    t_series.append(T[-1])
    w_series.append(W[-1])
    return t_series, w_series

def truncatedgaussian_sample(xmin, xmax, mu, sigma, x):
    roottwo = 1.414213562373095
    temp = math.erfc((-xmin + mu) / (roottwo * sigma)) + x * (
        math.erfc((-xmax + mu) / (roottwo * sigma)) - math.erfc((-xmin + mu) / (roottwo * sigma))
    )
    y = mu + roottwo * sigma * coolmath.inverf(-1.0 + temp)
    return y

 def test_erfc(self):
     import math
     assert math.erfc(0.0) == 1.0
     assert math.erfc(-0.0) == 1.0
     assert math.erfc(float("inf")) == 0.0
     assert math.erfc(float("-inf")) == 2.0
     assert math.isnan(math.erf(float("nan")))
     assert math.erfc(1e-308) == 1.0

 def RandIntVec(self,ListSize, ListSumValue, Distribution='Normal'):
     """
     Inputs:
     ListSize = the size of the list to return
     ListSumValue = The sum of list values
     Distribution = can be 'uniform' for uniform distribution, 'normal' for a normal distribution ~ N(0,1) with +/- 3 sigma  (default), or a list of size 'ListSize' or 'ListSize - 1' for an empirical (arbitrary) distribution. Probabilities of each of the p different outcomes. These should sum to 1 (however, the last element is always assumed to account for the remaining probability, as long as sum(pvals[:-1]) <= 1).  
     Output:
     A list of random integers of length 'ListSize' whose sum is 'ListSumValue'.
     """
     if type(Distribution) == list:
         DistributionSize = len(Distribution)
         if ListSize == DistributionSize or (ListSize-1) == DistributionSize:
             Values = multinomial(ListSumValue,Distribution,size=1)
             OutputValue = Values[0]
     elif Distribution.lower() == 'uniform': #I do not recommend this!!!! I see that it is not as random (at least on my computer) as I had hoped
         UniformDistro = [1/ListSize for i in range(ListSize)]
         Values = multinomial(ListSumValue,UniformDistro,size=1)
         OutputValue = Values[0]
     elif Distribution.lower() == 'normal':
         """
         Normal Distribution Construction....It's very flexible and hideous
         Assume a +-3 sigma range.  Warning, this may or may not be a suitable range for your implementation!
         If one wishes to explore a different range, then changes the LowSigma and HighSigma values
         """
         LowSigma    = -3#-3 sigma
         HighSigma   = 3#+3 sigma
         if (float(ListSize) - 1) == 0:
             StepSize    = 0
         else:
             StepSize    = 1/(float(ListSize) - 1)
         ZValues     = [(LowSigma * (1-i*StepSize) +(i*StepSize)*HighSigma) for i in range(int(ListSize))]
         #Construction parameters for N(Mean,Variance) - Default is N(0,1)
         Mean        = 0
         Var         = 1
         #NormalDistro= [self.NormalDistributionFunction(Mean, Var, x) for x in ZValues]
         NormalDistro= list()
         for i in range(len(ZValues)):
             if i==0:
                 ERFCVAL = 0.5 * math.erfc(-ZValues[i]/math.sqrt(2))
                 NormalDistro.append(ERFCVAL)
             elif i ==  len(ZValues) - 1:
                 ERFCVAL = NormalDistro[0]
                 NormalDistro.append(ERFCVAL)
             else:
                 ERFCVAL1 = 0.5 * math.erfc(-ZValues[i]/math.sqrt(2))
                 ERFCVAL2 = 0.5 * math.erfc(-ZValues[i-1]/math.sqrt(2))
                 ERFCVAL = ERFCVAL1 - ERFCVAL2
                 NormalDistro.append(ERFCVAL)  
         #print "Normal Distribution sum = %f"%sum(NormalDistro)
         Values = multinomial(ListSumValue,NormalDistro,size=1)
         OutputValue = Values[0]
     else:
         raise ValueError ('Cannot create desired vector')
     return OutputValue

def getnormcdf(x,lowertail=True):
  """
  Get the normal CDF function. used to calculate p-value
  """
  # ax=math.fabs(x)
  #axv=math.erfc(x/(2**0.5))/2; # higher tail
  if lowertail==False:
    #return axv
    return math.erfc(x/(2**0.5))/2
  else:
    #return 1-axv
    return math.erfc(-x/(2**0.5))/2

 def de_marsily_no_reaction(self,x,t):
   # Based on Equation 10.3.2 in Quantitative Hydrogeology, Ghislain de
   # Marsily, 1986.
   D_ = self.D
   U_over_porR = self.U/(self.porosity*self.R)
   two_sqrt_Dt_over_porR = 2.*math.sqrt(D_*t/(self.porosity*self.R))
   temp = 0.5* \
          (math.erfc((x-t*U_over_porR)/two_sqrt_Dt_over_porR) + \
           math.exp(self.U*x/D_) * \
           math.erfc((x+t*U_over_porR)/two_sqrt_Dt_over_porR))
   value = temp*(self.c1-self.c0) + self.c0
   return value

 def ogata_banks(self,x,t):
   # Based on "Solution of the Differential Equation of Longitudinal
   # Dispersion in Porous Media", USGS Professional Paper 411-A by
   # Akio Ogata and R.B. Banks, 1961.
   D_ = self.D/self.porosity
   v = self.U/self.porosity
   temp = 0.5* \
          (math.erfc((x-v/self.R*t)/(2.*math.sqrt(D_/self.R*t))) + \
           math.exp(v*x/D_) * \
           math.erfc((x+v/self.R*t)/(2.*math.sqrt(D_/self.R*t))))
   value = temp*(self.c1-self.c0) + self.c0
   return value

def blackscholes(stockprice, strike, riskfree, time, volatility):
	d1 = (log(stockprice/strike)+(riskfree+.5*volatility**2)*time)/(volatility*sqrt(time))
	d2 = d1 - volatility*sqrt(time)

	callprice = stockprice*.5*erfc(-d1/sqrt(2))-strike*exp(-riskfree*time)*.5*erfc(-d2/sqrt(2))

	putprice = strike*exp(-riskfree*time)*.5*erfc(d2/sqrt(2))-stockprice*.5*erfc(d1/sqrt(2))

	calldelta = .5*erfc(-d1/sqrt(2))

	putdelta = calldelta - 1

	return [callprice, putprice, calldelta, putdelta]

def calculate_results(representedPalindromes, kmer_counts, 
	    	  totalCharCount, hypothesesCount):
    """
    INPUT:
	representedPalindromes - A collections.Counter() 
		object
	kmer_counts - A collections.Counter() object. 
	totalCharCount - An integer 
	hypothesesCount - An integer
    OUTPUT:
	tupleList - A list of tuples, each tuple
        contains the palindrome name, observed count,
        expected count, z-score and E-value. 

    Takes in a collections.Counter representedPalindromes,
    a collections.Counter kmer_counts,
    an int totalCharCount, and an int hypothesesCount.  For 
    each palindrome in representedPalindromes the expected
    value, z-Score, and e-value is calculated. These three
    values are coupled with the palindrome name and observed
    count, to create a tuple with five elements.  These 
    tuples are loaded into a list, which is returned after 
    all palindromes in representedPalindromes has been 
    iterated over. 
    """
    tupleList = []
    # The list of tuples. Each element in this list is a 
    # tuple, with five elements. The tuple elements
    # are palindrome name, observed count, expected count,
    # z-score and E-value. 
    for key, value in representedPalindromes.items():
        if (len(key)%2) == 1:
            keyList = list(key)
            if(keyList[len(key)/2]) in "ATCG": continue
	    # Throw out odd palindromes with middle characters
	    # ATCG, they are not desired.  
        expectedValue = expected_value(key,kmer_counts)
        standardDeviation = (math.sqrt(expectedValue *
		 (1 - (expectedValue/ totalCharCount))))
        zScore = ((kmer_counts[key] - expectedValue)
		/standardDeviation)
        if zScore <= 0:
           e_Value = ((math.erfc(-zScore/math.sqrt(2))/2)*
		2*hypothesesCount)
        else:
           e_Value = ((math.erfc(zScore/math.sqrt(2))/2)
		*2*hypothesesCount)
        tupleList.append((key, str(value), str(expectedValue),
		 zScore, e_Value))
    return tupleList

def rspace_sum(rs, qs, basis, kappa, rlim):
    """
    Real space part of the sum
    Parameters:
       rs -- list of particle positions
       qs -- list of particle charges
       basis -- real space basis
       kappa -- splitting parameter
       rlim -- size of lattice (one side of a cube of points)
    """
    rspace_sum = 0.0
    for i in range(len(rs)):
        for j in range(len(rs)):
            q = qs[i]*qs[j]
            r = rs[j] - rs[i]
            for n1 in range(-rlim+1, rlim):
                for n2 in range(-rlim+1, rlim):
                    for n3 in range(-rlim+1, rlim):
                        if i == j and n1 == 0 and n2 == 0 and n3 == 0:
                            continue
                        lat = n1*basis[0] + n2*basis[1] + n3*basis[2]
                        d = r + lat
                        rd = math.sqrt(d.dot(d))
                        rspace_sum += q * math.erfc(kappa*rd)/rd
    return rspace_sum

def compute_kullback_leibler_check_statistic(n=100, prngstate=None):
    """Compute the lowest of the survival function and the CDF of the exact KL
    divergence KL(N(mu1,s1)||N(mu2,s2)) w.r.t. the sample distribution of the
    KL divergence drawn by computing log(P(x|N(mu1,s1)))-log(P(x|N(mu2,s2)))
    over a sample x~N(mu1,s1). If we are computing the KL divergence
    accurately, the exact value should fall squarely in the sample, and the
    tail probabilities should be relatively large.

    """
    if prngstate is None:
        raise TypeError('Must explicitly specify numpy.random.RandomState')
    mu1 = mu2 = 0
    s1 = 1
    s2 = 2
    exact = gaussian_kl_divergence(mu1, s1, mu2, s2)
    sample = prngstate.normal(mu1, s1, n)
    lpdf1 = gaussian_log_pdf(mu1, s1)
    lpdf2 = gaussian_log_pdf(mu2, s2)
    estimate, std = kl.kullback_leibler(sample, lpdf1, lpdf2)
    # This computes the minimum of the left and right tail probabilities of the
    # exact KL divergence vs a gaussian fit to the sample estimate. There is a
    # distinct negative skew to the samples used to compute `estimate`, so this
    # statistic is not uniform. Nonetheless, we do not expect it to get too
    # small.
    return erfc(abs(exact - estimate) / std) / 2

 def F8_f(self):
     self.F8 = 0.5 * math.erfc(-self.z_score/math.sqrt(2))
     if self.F8 == 0:
         self.F8 = 10^-16
     else:
         self.F8 = self.F8
     return self.F8

  def Ea(self,n_1,r):
      '''
      This function computes Easiness of task dependent on Reaction Time and Accuracy
      Arguments:
      n_1 - One of the numerical numbers in the task
      r - Ratio of numbers in the task

      In order to compute the P_Error an erfc formula consisting of the two numbers and the weber fraction is used
      '''
      m = (self.m)
      intercept = self.i
      w = self.w

      dist = int(round(n_1 / r)) - n_1
      rt = (dist * m) + intercept # Computing Reaction_Time from the values obtained from the RT graph
      n_2 = int(round(n_1/r))
      numer = abs(n_1-n_2)
      denom = math.sqrt(2)*w*(((n_1**2)+(n_2**2))**0.5)
      P_Err = 0.5*math.erfc(numer/denom)
      P_A = 1 - P_Err
      #print(P_A)
      array_poss = np.random.choice([0,1],size=(10),p=[1-P_A, P_A]) # Generating bin values array (consisting of 0's and 1's)using P_Acc probability
      val = np.random.choice(array_poss) # Randomly choosing samples from array_poss
      rt = np.random.normal(rt,intercept) # Normal Distribution sampling of RT vals
      if rt<500:
          rt = 500 # re-evaluate RT's < 500
          rt = np.random.normal(rt,130)
      E = val - (rt/2000.) # Computing discrete Easiness value of a question

      return E, val, rt

    def _calc_real_and_point(self):
        """
        Determines the self energy -(eta/pi)**(1/2) * sum_{i=1}^{N} q_i**2
        
        If cell is charged a compensating background is added (i.e. a G=0 term)
        """
        all_nn = self._s.get_all_neighbors(self._rmax, True)

        forcepf = 2.0 * self._sqrt_eta / sqrt(pi)
        coords = self._coords
        numsites = self._s.num_sites
        ereal = np.zeros((numsites, numsites))
        epoint = np.zeros((numsites))
        forces = np.zeros((numsites, 3))
        for i in xrange(numsites):
            nn = all_nn[i] #self._s.get_neighbors(site, self._rmax)
            qi = self._oxi_states[i]
            epoint[i] = qi * qi
            epoint[i] *= -1.0 * sqrt(self._eta / pi)
            epoint[i] += qi * pi / (2.0 * self._vol * self._eta)  #add jellium term
            for j in range(len(nn)):  #for (nsite, rij)  in nn:
                nsite = nn[j][0]
                rij = nn[j][1]
                qj = compute_average_oxidation_state(nsite)
                erfcval = erfc(self._sqrt_eta * rij)
                ereal[nn[j][2], i] += erfcval * qi * qj / rij
                fijpf = qj / pow(rij, 3) * (erfcval + forcepf * rij * exp(-self._eta * pow(rij, 2)))
                forces[i] += fijpf * (coords[i] - nsite.coords) * qi * EwaldSummation.CONV_FACT

        ereal = ereal * 0.5 * EwaldSummation.CONV_FACT
        epoint = epoint * EwaldSummation.CONV_FACT
        return (ereal, epoint, forces)

def graph_FWHM_data_range(start_date=datetime.datetime(2015,3,6),
                          end_date=datetime.datetime(2015,4,15),tenmin=True,
                          path='/home/douglas/Dropbox (Thacher)/Observatory/Seeing/Data/',
                          write=True,outpath='./'):
    
    
    plot_params()
    fwhm = get_FWHM_data_range(start_date = start_date, end_date=end_date, path=path, tenmin=tenmin)

    # Basic stats
    med = np.median(fwhm)
    mean = np.mean(fwhm)
    fwhm_clip, low, high = sigmaclip(fwhm,low=3,high=3)
    meanclip = np.mean(fwhm_clip)

    # Get mode using kernel density estimation (KDE)
    vals = np.linspace(0,30,1000)
    fkde = gaussian_kde(fwhm)
    fpdf = fkde(vals)
    mode = vals[np.argmax(fpdf)]
    std = np.std(fwhm)


    plt.ion()
    plt.figure(99)
    plt.clf()
    plt.hist(fwhm, color='darkgoldenrod',bins=35)
    plt.xlabel('FWHM (arcsec)',fontsize=16)
    plt.ylabel('Frequency',fontsize=16)
    plt.annotate('mode $=$ %.2f" ' % mode, [0.87,0.85],horizontalalignment='right',
                 xycoords='figure fraction',fontsize='large')
    plt.annotate('median $=$ %.2f" ' % med, [0.87,0.8],horizontalalignment='right',
                 xycoords='figure fraction',fontsize='large')
    plt.annotate('mean $=$ %.2f" ' % mean, [0.87,0.75],horizontalalignment='right',
                 xycoords='figure fraction',fontsize='large')

    xvals = np.linspace(0,30,1000)
    kde = gaussian_kde(fwhm)
    pdf = kde(xvals)
    dist_c = np.cumsum(pdf)/np.sum(pdf)
    func = interp1d(dist_c,vals,kind='linear')
    lo = np.float(func(math.erfc(1./np.sqrt(2))))
    hi = np.float(func(math.erf(1./np.sqrt(2))))

    disthi = np.linspace(.684,.999,100)
    distlo = disthi-0.6827
    disthis = func(disthi)
    distlos = func(distlo)

    interval = np.min(disthis-distlos)

    plt.annotate('1 $\sigma$ int. $=$ %.2f" ' % interval, [0.87,0.70],horizontalalignment='right',
                 xycoords='figure fraction',fontsize='large')
    
    
    plt.rcdefaults()

    plt.savefig(outpath+'Seeing_Cumulative.png',dpi=300)

    return

 def bear(self,x,t):
   # Based on Equation 10.6.22 in "Dynamics of Fluids in Porous Media", 
   # Jacob Bear, 1988.
   D_ = self.D/self.porosity
   v = self.U/self.porosity
   beta = math.sqrt(v*v/(4.*D_*D_)+ self.lam/D_)
   temp = 0.5 * \
          math.exp(v*x/(2.*D_)) * \
          (math.exp(-1.*beta*x) * \
           math.erfc((x-math.sqrt(v*v+4.*self.lam*D_)*t) / \
                     (2.*math.sqrt(D_*t))) + \
           math.exp(beta*x) * \
           math.erfc((x+math.sqrt(v*v+4.*self.lam*D_)*t) / \
                     (2.*math.sqrt(D_*t))))
   value = temp*(self.c1-self.c0) + self.c0
   return value

def monobit(bits):
	ceros=	bits.count("0")
	unos =  bits.count("1")
	frequencia = len(bits)
	SN = ceros-unos
	test = abs(SN)/sqrt(frequencia)
	return erfc(test/sqrt(2))

def distparams(dist):
    """
    Description:
    ------------
    Return robust statistics of a distribution of data values

    Example:
    --------
    med,mode,interval,lo,hi = distparams(dist)
    """

    from scipy.stats.kde import gaussian_kde
    from scipy.interpolate import interp1d
    vals = np.linspace(np.min(dist)*0.5,np.max(dist)*1.5,1000)
    kde = gaussian_kde(dist)
    pdf = kde(vals)
    dist_c = np.cumsum(pdf)/np.sum(pdf)
    func = interp1d(dist_c,vals,kind='linear')
    lo = np.float(func(math.erfc(1./np.sqrt(2))))
    hi = np.float(func(math.erf(1./np.sqrt(2))))
    med = np.float(func(0.5))
    mode = vals[np.argmax(pdf)]

    disthi = np.linspace(.684,.999,100)
    distlo = disthi-0.6827
    disthis = func(disthi)
    distlos = func(distlo)
    
    interval = np.min(disthis-distlos)

    return med,mode,interval,lo,hi

def TraditionalEwald(L,D,rs,nMax):
    alpha = 7./(L)

    VShort = lambda r,a: math.erfc(a*r)/r
    fVLong = lambda k2,a: 4*pi*math.exp(-k2/(4*a*a))/(k2*(L**D))

    # Short-ranged r-space
    Vss = []
    for r in rs:
      Vs = VShort(r,alpha)
      Vss.append(Vs)

    # Long-ranged k-space
    ns = Genns(nMax,D)
    ks = []
    fVls = []
    for n in ns:
      magn2 = dot(n,n)
      k2 = 4.*pi*pi*magn2/(L*L)
      ks.append(math.sqrt(k2))
      if magn2 != 0:
        fVl = fVLong(k2,alpha)
        fVls.append(fVl)
      else:
        fVls.append(0.)

    ks, fVls = GetUnique(ks,fVls)
    for (k,fVl) in sorted(zip(ks,fVls)):
      print k, fVl

    return [Vss,fVls]

def monobit (bits, p_value=0.01):
    """Monobit Test
    returns a tuple representing (metric, pass)"""
    S = bits.count(1) - bits.count(0) 
    S_obs = abs(S) / math.sqrt(len(bits))
    P = math.erfc(S_obs / math.sqrt(2))
    return P, P >= p_value