def linearinterpolation(nstep,ndata,dt):
    """used for lienar interpolation between transportation matrices.
    returns weights alpha, beta and indices of matrices.
    Parameters
    -------
    nstep   : int Number of timesteps
    ndata   : int Number of matrices
    Returns
    -------
    alpha,beta : array 
                    coefficients for interpolation
    jalpha,jbeta : array
                    indices for interpolation
    """

    t = np.zeros(nstep,dtype=np.float_)
    for i in range(nstep):
        t[i] = np.fmod(0 + i*dt, 1.0)


    beta    = np.array(nstep,dtype=np.float_)
    alpha   = np.array(nstep,dtype=np.float_)

    w       = t * ndata+0.5
    beta    = np.float_(np.fmod(w, 1.0))
    alpha   = np.float_(1.0-beta)
    jalpha  = np.fmod(np.floor(w)+ndata-1.0,ndata).astype(int)
    jbeta   = np.fmod(np.floor(w),ndata).astype(int)

    return alpha,beta,jalpha,jbeta

def fundamental_arguments(t):

    """Compute the fundamental arguments (mean elements) of Sun and Moon.

    `t` - TDB time in Julian centuries since J2000.0, as float or NumPy array

    Outputs fundamental arguments, in radians:
          a[0] = l (mean anomaly of the Moon)
          a[1] = l' (mean anomaly of the Sun)
          a[2] = F (mean argument of the latitude of the Moon)
          a[3] = D (mean elongation of the Moon from the Sun)
          a[4] = Omega (mean longitude of the Moon's ascending node);
                 from Simon section 3.4(b.3),
                 precession = 5028.8200 arcsec/cy)

    """
    a = fa4 * t
    a += fa3
    a *= t
    a += fa2
    a *= t
    a += fa1
    a *= t
    a += fa0
    fmod(a, ASEC360, out=a)
    a *= ASEC2RAD
    if getattr(t, 'shape', ()):
        return a
    return a[:,0]

 def _visiting(self, step):
     """
     Assignement of components values based on visiting distribution.
     The way of exploring space depends on the Markov chain stepping
     """
     # It it is the first part of the markov chain
     # Changing all components at the same time
     if step < self._x.size:
         visits = np.array([self._visita() for _ in range(self._x.size)])
         visits[visits > 1.e8] = 1.e8 * self._random_state.random_sample()
         visits[visits < -1e8] = -1.e8 * self._random_state.random_sample()
         self._x = visits + self._xbackup
         a = self._x - self._lower
         b = np.fmod(a, self._xrange) + self._xrange
         self._x = np.fmod(b, self._xrange) + self._lower
         self._x[np.fabs(self._x - self._lower) < 1.e-10] += 1.e-10
     else:
         # Second part of the markov chain
         # Now change only one component at a time
         visit = self._visita()
         if visit > 1.e8:
             visit = 1.e8 * self._random_state.random_sample()
         elif visit < -1e8:
             visit = -1.e8 * self._random_state.random_sample()
         index = step - self._x.size
         self._x[index] = visit + self._xbackup[index]
         a = self._x[index] - self._lower[index]
         b = np.fmod(
             a, self._xrange[index]) + self._xrange[index]
         self._x[index] = np.fmod(
             b, self._xrange[index]) + self._lower[index]
         if np.fabs(self._x[index] - self._lower[
                 index]) < 1.e-10:
             self._x[index] += 1.e-10

def mod180deg(x):
    if x >= 0:
        retval = fmod(x+180.0, 360.0)-180.0
    else:
        retval = -(fmod(-x+180.0, 360.0)-180.0)
    assert -180.0 <= retval <= 180.0
    return retval

def fold(spike_trains, period):
    """Fold `spike_trains` by `period`."""

    # data = {key:[] for key in spike_trains.dtype.names}

    rows = []
    for i,row in spike_trains.iterrows():
        period_num = int( np.ceil(row['duration'] / period) )
        last_period = np.fmod(row['duration'], period)

        spikes = row['spikes']
        for idx in range(period_num):
            lo = idx * period
            hi = (idx+1) * period
            sec = spikes[(spikes>=lo) & (spikes<hi)]
            sec = np.fmod(sec, period)

            r = row.copy()
            r['spikes'] = sec
            r['duration'] = period

            rows.append(r)

        if last_period > 0:
            rows[-1]['duration'] = last_period


    folded_trains = pd.DataFrame(rows)
    folded_trains = folded_trains.reset_index(drop=True)

    return folded_trains

    def render(self, model, params, frame):

        # Scalar animation parameter, based on height and distance
        d = model.edgeCenters[:,2] + 0.5 * model.edgeDistances
        numpy.multiply(d, 1/self.height, d)

        # Add global offset for Z scrolling over time
        numpy.add(d, params.time * self.speed, d)

        # Add an offset that depends on which tree we're in
        numpy.add(d, numpy.choose(model.edgeTree, self.offsets), d)

        # Periodic animation, stored in our color table. Linearly interpolate.
        numpy.fmod(d, self.period, d)
        color = numpy.empty((model.numLEDs, 3))
        color[:,0] = numpy.interp(d, self.colorX, self.colorY[:,0])
        color[:,1] = numpy.interp(d, self.colorX, self.colorY[:,1])
        color[:,2] = numpy.interp(d, self.colorX, self.colorY[:,2])

        # Random flickering noise
        noise = numpy.random.rand(model.numLEDs).reshape(-1, 1)
        numpy.multiply(noise, 0.25, noise)
        numpy.add(noise, 0.75, noise)

        numpy.multiply(color, noise, color)
        numpy.add(frame, color, frame)

def _space_grid(Sm, spacing):
    '''
    Zero everywhere, except for the points at a certain grid interval.
    '''
    Ix, Iy, Iz = np.indices(Sm.shape)
    return Sm * (np.fmod(Ix, spacing) == 0) * (np.fmod(Iy, spacing) == 0) * \
           (np.fmod(Iz, spacing) == 0)

def rebin_counts_2D_indexing_new(x, y, I, xo, yo, Io):
    # use this one for 2D.
    csx = np.empty((len(x), len(y)-1))
    csx[0, :] = 0.0
    csx[1:, :] = np.cumsum(I, axis=0)

    xindices = np.interp(xo, x, np.arange(len(x), dtype='float'))
    xindices_whole = np.floor(xindices).astype(int)
    xindices_frac = np.fmod(xindices, 1.0)
    # the only way an index will match the highest bin edge is for lookups at or outside the range of the
    # source.  In this case the fractional portion will always == 0.0, so clipping this to the bin edge below
    # is safe and allows us to use the source weights unmodified
    xintegral = csx[xindices_whole, :] + I[xindices_whole.clip(max=len(x)-2), :]*xindices_frac[:, None]
    # rebinned over x
    ix = np.diff(xintegral, axis=0)

    csy = np.empty((len(xo)-1, len(y)))
    csy[:, 0] = 0.0
    csy[:, 1:] = np.cumsum(ix, axis=1)

    yindices = np.interp(yo, y, np.arange(len(y), dtype='float'))
    yindices_whole = np.floor(yindices).astype(int)
    yindices_frac = np.fmod(yindices, 1.0)
    yintegral = csy[:, yindices_whole] + ix[:, yindices_whole.clip(max=len(y)-2)]*yindices_frac[None, :]
    # rebinned over x and y
    ixy = np.diff(yintegral, axis=1)

    Io += ixy

def historyids_to_indices(historyids, historyScores): 
	hids=np.array(historyids, dtype=int)
	itr=np.fmod(hids, Global_BINWIDTH)
	sim=np.round(hids/Global_BINWIDTH)
	runlen=max(np.fmod(historyScores[:,0], Global_BINWIDTH))+1
	newi=itr+sim*runlen
	return newi.astype('int')

def compute_lscore_over_time(events, historyScores, run, outputfn):
    runlen = max(np.fmod(historyScores[:, 0], Global_BINWIDTH)) + 1
    myrun = run
    myeventids = ["total"]
    myruncounts = [1]
    hidmin = myrun * Global_BINWIDTH
    hidmax = hidmin + runlen
    tmpi = np.where((historyScores[:, 0] >= hidmin) & (historyScores[:, 0] <= hidmax))[0]
    costsarray = np.mean(historyScores[tmpi, 1:3], axis=1)
    totalscores = np.cumsum(np.exp(-1 * Global_K * costsarray))
    mydata = totalscores / totalscores[runlen - 1]
    for event in events:
        event.histories = listout_ranges(event.histRanges)
        myeventids.append(event.id)
        sys.stderr.write("working on event: %s\n" % (event.id))
        myruncounts.append(event.numsims)
        lscores = np.zeros(runlen)
        hids = np.array(event.histories)
        myhids = hids[np.where((hids <= hidmax) & (hids >= hidmin))]
        myhis = np.fmod(hids, Global_BINWIDTH)
        for i in xrange(myhis.size):
            totalscore = totalscores[myhis[i]]
            lscore = compute_likelihood_histories(myhids[: (i + 1)], historyScores, totalscore)
            lscores[myhis[i]] = lscore
        for i in xrange(max(myhis), hidmax):
            totalscore = totalscores[i]
            lscore = compute_likelihood_histories(myhids, historyScores, totalscore)
            lscores[i] = lscore
        mydata = np.vstack((mydata, lscores))
    np.savetxt(
        outputfn, mydata.T, delimiter="\t", header="\t".join(myeventids) + "\n" + "\t".join(map(str, myruncounts))
    )

 def rotate(self, angle, mask=None):
     """Rotate the grids (arena centered)
     
     Grids to be rotated can be optionally specified by bool/index array
     *mask*, otherwise population is rotated. Specified *angle* can be a
     scalar value to be applied to the population or a population- or
     mask-sized array depending on whether *mask* is specified.
     """
     rot2D = lambda psi: [[cos(psi), sin(psi)], [-sin(psi),  cos(psi)]]
     if mask is not None and type(mask) is np.ndarray:
         if mask.dtype.kind == 'b':
             mask = mask.nonzero()[0]
         if type(angle) is np.ndarray and angle.size == mask.size:
             for i,ix in enumerate(mask):
                 self._phi[ix] = np.dot(self._phi[ix], rot2D(angle[i]))
         elif type(angle) in (int, float, np.float64):
             angle = float(angle)
             self._phi[mask] = np.dot(self._phi[mask], rot2D(angle))
         else:
             raise TypeError, 'angle must be mask-sized array or float'
         self._psi[mask] = np.fmod(self._psi[mask]+angle, 2*pi)
     elif mask is None:
         if type(angle) is np.ndarray and angle.size == self.num_maps:
             for i in xrange(self.num_maps):
                 self._phi[i] = np.dot(self._phi[i], rot2D(angle[i]))
         elif type(angle) in (int, float, np.float64):
             angle = float(angle)
             self._phi = np.dot(self._phi, rot2D(angle))
         else:
             raise TypeError, 'angle must be num_maps array or float'
         self._psi = np.fmod(self._psi+angle, 2*pi)
     else:
         raise TypeError, 'mask must be bool/index array'

def plot_fft_results(freq, pow, cum, angles, x, y, freqi, powi, angi):
	plt.figure(figsize=(10,15))
	plt.subplot(5,1,1)
	plt.plot(freq, pow, 'x')
	plt.gca().set_xscale('log')
	plt.xlabel('frequency')
	plt.ylabel('power')
	plt.hlines(powi, plt.xlim()[0], plt.xlim()[1], alpha=0.3)
	plt.subplot(5,1,2)
	plt.plot(freq, cum, 'x')
	plt.gca().set_xscale('log')
	plt.ylabel('cumulative probability')
	plt.xlabel('frequency')
	plt.vlines(freqi, 0, 1, alpha=0.3)
	plt.subplot(5,1,3)
	plt.plot(cum, pow**0.5, 'x')
	plt.xlabel('cumulative probability')
	plt.ylabel('amplitude')
	plt.hlines(powi**0.5, plt.xlim()[0], plt.xlim()[1], alpha=0.3)
	plt.subplot(5,1,4)
	plt.plot(cum, numpy.fmod(angles+10*numpy.pi, 2*numpy.pi), 'x')
	plt.xlabel('cumulative probability')
	plt.ylabel('phase')
	plt.hlines(numpy.fmod(angi+10*numpy.pi, 2*numpy.pi), plt.xlim()[0], plt.xlim()[1], alpha=0.3)
	plt.subplot(5,1,5)
	plt.plot(x, y, 'x')
	plt.xlabel('time')
	plt.ylabel('value')
	xnew = numpy.linspace(x.min(), x.max(), 200)
	ynew = powi**0.5 * numpy.cos(xnew * freqi * pi * 2 + angi)
	plt.plot(xnew, ynew, '-', alpha=0.5, lw=3, 
		label="freq=%.5f\npow=%.3f\nang=%.3f" % (freqi, powi, angi))
	plt.legend(loc='best', prop=dict(size=10))

 def test_impulse_negative(self):
     """ check the transform of a negative impulse at a random place """
     win_s = 256
     i = int(floor(random()*win_s))
     impulse = -.1
     f = fft(win_s)
     timegrain = fvec(win_s)
     timegrain[0] = 0
     timegrain[i] = impulse
     fftgrain = f ( timegrain )
     #self.plot_this ( fftgrain.phas )
     assert_almost_equal ( fftgrain.norm, abs(impulse), decimal = 5 )
     if impulse < 0:
         # phase can be pi or -pi, as it is not unwrapped
         #assert_almost_equal ( abs(fftgrain.phas[1:-1]) , pi, decimal = 6 )
         assert_almost_equal ( fftgrain.phas[0], pi, decimal = 6)
         assert_almost_equal ( np.fmod(fftgrain.phas[-1], pi), 0, decimal = 6)
     else:
         #assert_equal ( fftgrain.phas[1:-1] == 0, True)
         assert_equal ( fftgrain.phas[0], 0)
         assert_almost_equal ( np.fmod(fftgrain.phas[-1], pi), 0, decimal = 6)
     # now check the resynthesis
     synthgrain = f.rdo ( fftgrain )
     #self.plot_this ( fftgrain.phas.T )
     assert_equal ( fftgrain.phas <= pi, True)
     assert_equal ( fftgrain.phas >= -pi, True)
     #self.plot_this ( synthgrain - timegrain )
     assert_almost_equal ( synthgrain, timegrain, decimal = 6 )

def analytic_phase_reset_old(phi, dx=0., dy=0.,
        a=0., l_ccw=0.2, l_cw=-1., X=1., Y=1.):

    # this is an older expression derived with several assumptions
    assert(X == 1.)
    assert(Y == 1.)
    assert(l_cw == -1.)

    r0 = iris_fixedpoint(a, l_ccw, l_cw, X, Y)
    T = iris_period(a, l_ccw, l_cw, X, Y)
    if r0 == None:
        raise RuntimeError("No limit cycle found")
    else:
        quad1or2 = np.fmod(phi, 2 * math.pi) < math.pi
        quad1or3 = np.fmod(phi, math.pi) < math.pi/2
        du = (quad1or3 * dy + (1 - quad1or3) * dx) * (1 - 2*quad1or2)
        ds = (quad1or3 * dx + (1 - quad1or3) * dy) * (
                1 - 2*quad1or2*quad1or3 - 2*(1 - quad1or2)*(1 - quad1or3))
        t = np.fmod(phi, math.pi/2)/(math.pi/2) * T/4
        Q = 1/(l_ccw * r0) * np.exp(-l_ccw * t)
        dt0 = -Q * du
        dr = np.exp(-T/4) * (X * -dt0 + np.exp(t) * ds)
        return (
                (dt0
                + -1./(l_ccw * r0)
                * 1./(1 - 1./(l_ccw * r0) * (r0/Y)**(1/l_ccw)) * dr)
                / T * 2*math.pi
                )

def _level1_xfm_no_highpass(X, h0o, h1o, ext_mode):
    """Perform level 1 of the 3d transform discarding highpass subbands.

    """
    # Check shape of input according to ext_mode. Note that shape of X is
    # double original input in each direction.
    if ext_mode == 4 and np.any(np.fmod(X.shape, 2) != 0):
        raise ValueError('Input shape should be a multiple of 2 in each direction when self.ext_mode == 4')
    elif ext_mode == 8 and np.any(np.fmod(X.shape, 4) != 0):
        raise ValueError('Input shape should be a multiple of 4 in each direction when self.ext_mode == 8')

    out = np.zeros_like(X)

    # Loop over 2nd dimension extracting 2D slice from first and 3rd dimensions
    for f in xrange(X.shape[1]):
        # extract slice
        y = X[:, f, :].T
        out[:, f, :] = colfilter(y, h0o).T

  # Loop over 3rd dimension extracting 2D slice from first and 2nd dimensions
    for f in xrange(X.shape[2]):
        y = colfilter(out[:, :, f].T, h0o).T
        out[:, :, f] = colfilter(y, h0o)

    return out

 def visiting(self, x, step, temperature):
     dim = x.size
     if step < dim:
         # Changing all coordinates with a new visting value
         visits = np.array([self.visit_fn(
             temperature) for _ in range(dim)])
         upper_sample = self.rs.random_sample()
         lower_sample = self.rs.random_sample()
         visits[visits > self.tail_limit] = self.tail_limit * upper_sample
         visits[visits < -self.tail_limit] = -self.tail_limit * lower_sample
         x_visit = visits + x
         a = x_visit - self.lower
         b = np.fmod(a, self.b_range) + self.b_range
         x_visit = np.fmod(b, self.b_range) + self.lower
         x_visit[np.fabs(
             x_visit - self.lower) < self.min_visit_bound] += 1.e-10
     else:
         # Changing only one coordinate at a time based on Markov chain step
         x_visit = np.copy(x)
         visit = self.visit_fn(temperature)
         if visit > self.tail_limit:
             visit = self.tail_limit * self.rs.random_sample()
         elif visit < -self.tail_limit:
             visit = -self.tail_limit * self.rs.random_sample()
         index = step - dim
         x_visit[index] = visit + x[index]
         a = x_visit[index] - self.lower[index]
         b = np.fmod(a, self.b_range[index]) + self.b_range[index]
         x_visit[index] = np.fmod(b, self.b_range[
             index]) + self.lower[index]
         if np.fabs(x_visit[index] - self.lower[
                 index]) < self.min_visit_bound:
             x_visit[index] += self.min_visit_bound
     return x_visit

    def transit_depths(self, p):
        p = self.to_params(p)

        dt = self.times[1] - self.times[0]

        depths = np.zeros(self.times.shape)

        left_bounds = self.times - dt/2.0
        right_bounds = self.times + dt/2.0

        left_time_since_transit = np.fmod(left_bounds - p['t0'], p['P'])
        right_time_since_transit = np.fmod(right_bounds - p['t0'], p['P'])

        left_in_transit = (left_time_since_transit > 0) & (left_time_since_transit < p['T'])
        right_in_transit = (right_time_since_transit > 0) & (right_time_since_transit < p['T'])

        depths[left_in_transit & right_in_transit] = 1.0

        entries = right_in_transit & (~left_in_transit)
        exits = left_in_transit & (~right_in_transit)

        depths[entries] = right_time_since_transit[entries] / dt
        depths[exits] = (p['T']-left_time_since_transit[exits])/dt        

        return depths

	def FindPosition(self,pos,dx,dy):
		"""Finds which index the walker belongs to. 
		Implements periodic boundary conditions on the walk-area
		Should work in 1d as well now"""
		indx = [-1,-1]
		if self.d==1:
			for i in xrange(len(self.x)+1):
				# print 'i = %d, x[i] = '%i,self.x[i],pos-self.x[i]
				if np.abs(pos-self.x[i])<dx/2.0:
					"This test must be implemented in 2D, and in self.checkpos()!"
					return i
		elif self.d==2:
			if pos[0]>self.x1_ or pos[0]<self.x0_:
				print 'this should not happen now! pos[x] = ',pos[0]
				pos[0] = np.fmod(pos[0],(self.X[0,-1]-self.X[0,0])+dx)+self.X[0,0]
				pos[0] *= ((self.X[0,-1]-self.X[0,0])+dx) if pos[0]<0 else 1
			for i in xrange(len(self.X)):
				if np.abs(pos[0]-self.X[i,i])<dx/2.0:
					indx[0] = i
					break
			if pos[1]>self.y1_ or pos[1]<self.y0_:
				print 'this should not happen now, pos[y] = ',pos[1]
				pos[1] = np.fmod(pos[1],(self.Y[-1,0]-self.Y[0,0])+dy)+self.Y[0,0]
				pos[1] *= ((self.Y[-1,0]-self.Y[0,0])+dx) if pos[1]<0 else 1
				# print pos[1]
			for j in xrange(len(self.Y)):
				if np.abs(pos[1]-self.Y[j,j])<dy/2.0:
					indx[1] = j
					break
			if indx[0] == -1 or indx[-1] == -1:
				print 'økadjs g'
			# print indx
			return indx

def fix_ps1_coord_bugs(cols, file, hduname):
	# FIXME: Work around PS1 bugs
	ra, dec = cols['ra'], cols['dec']

	if np.any(ra < 0):
		ra[ra < 0] = np.fmod(np.fmod(ra[ra < 0], 360) + 360, 360)

	if np.any(ra >= 360):
		ra[ra >= 360] = np.fmod(np.fmod(ra[ra >= 360], 360) + 360, 360)

	if np.any(np.abs(dec) > 90):
		logger.warning("Encountered %d instances of dec > +/- 90 in file %s, %s HDU. Truncating to +/-90." % (sum(np.abs(dec) > 90), file, hduname))
		dec[dec > 90] = 90
		dec[dec < -90] = -90

	cols['ra'], cols['dec'] = ra, dec

	# Remove any NaN rows
	if np.isnan(cols['ra'].sum()):
		logger.warning("Encountered %d instances of ra == NaN in file %s, %s HDU. Removing those rows." % (sum(np.isnan(ra)), file, hduname))
		keep = ~np.isnan(cols['ra'])
		for name in cols: cols[name] = cols[name][keep]

	if np.isnan(cols['dec'].sum()):
		logger.warning("Encountered %d instances of dec == NaN in file %s, %s HDU. Removing those rows." % (sum(np.isnan(dec)), file, hduname))
		keep = ~np.isnan(cols['dec'])
		for name in cols: cols[name] = cols[name][keep]

        def test_infix(lhs, rhs):
            ovl0 = no_op(lhs)

            def test_np(fcn):
                ovl_l = fcn(ovl0, rhs)
                ovl_r = fcn(rhs, ovl0)
                ovl_l, ovl_r = evaluate([ovl_l, ovl_r])
                np_l = fcn(lhs, rhs)
                np_r = fcn(rhs, lhs)

                assert np.all(np.equal(ovl_l, np_l))
                assert np.all(np.equal(ovl_r, np_r))

            test_np(lambda x, y: x + y)
            test_np(lambda x, y: x - y)
            test_np(lambda x, y: x * y)
            test_np(lambda x, y: x / y)
            test_np(lambda x, y: x == y)
            test_np(lambda x, y: x != y)
            test_np(lambda x, y: x < y)
            test_np(lambda x, y: x <= y)
            test_np(lambda x, y: x > y)
            test_np(lambda x, y: x >= y)

            # OVL uses c-style fmod, not python style mod, so use numpy fmod function for test
            # see : http://docs.scipy.org/doc/numpy/reference/generated/numpy.fmod.html
            ovl_left, ovl_right = evaluate([ovl0 % rhs, rhs % ovl0])
            np_left = np.fmod(lhs, rhs)
            np_right = np.fmod(rhs, lhs)

            assert np.all(np.equal(ovl_left, np_left))
            assert np.all(np.equal(ovl_right, np_right))

            ovl_neg = evaluate([-ovl0])
            assert np.all(np.equal(ovl_neg, -lhs))