def randgen(key, image_size):
    
    from numpy import zeros, remainder, uint8
    dt = 0.001
    n = image_size
    
    xs = zeros((n + 1,))
    ys = zeros((n + 1,))
    zs = zeros((n + 1,))
    
    xs[0], ys[0], zs[0] = key[0], key[1], key[2]
    
    for i in xrange(n) :

        x_dot, y_dot, z_dot = lorenz(xs[i], ys[i], zs[i])
        xs[i + 1] = xs[i] + (x_dot * dt)
        ys[i + 1] = ys[i] + (y_dot * dt)
        zs[i + 1] = zs[i] + (z_dot * dt)
        
    Xs =  remainder(abs(xs*10**14), 2).astype(uint8)
    Ys =  remainder(abs(ys*10**14), 2).astype(uint8)
    Zs =  remainder(abs(zs*10**14), 2).astype(uint8)
    
    rand_array = Xs ^ Ys ^ Zs
    
    return rand_array

def get_spin_density_file(filename):
# edits a CHGCAR file to have only the spin difference table, "spin_up - spin_down"
	import re
	from StringIO import StringIO
	f = open(filename)
	lines = f.read().splitlines()
	n_atoms = numpy.sum(numpy.genfromtxt(StringIO(lines[6]),dtype=(int)))
	n_points = numpy.product(numpy.genfromtxt(StringIO(lines[9+n_atoms]),dtype=(int)))
	f.close
	# Remove the first spin table
	#get start line of the potential table for majority spin
	start_line = 10+n_atoms
	#get lastline of the potential table for majority spin
	if numpy.remainder(n_points,5) == 0:
		last_line = 9+n_atoms+n_points/5
	elif numpy.remainder(n_points,5) != 0:
		last_line = 9+n_atoms+n_points/5+1
	del lines[start_line:last_line]

	# delete lines until you next match the "number of grid points line"
	finished = 0;
	count = 0;
	while finished != 1:
		l_match = re.match(lines[9+n_atoms],lines[9+n_atoms+1],0)
		if (l_match):
			finished = 1
		else:
			del lines[9+n_atoms+1]

	del lines[9+n_atoms+1]
	outfile = 'CHGCAR-spin_density'
	fout = open(outfile,'w')
	for line in lines:
		print>>fout, line
	fout.close

def nmer_neighbors3d(xs,ys,zs,sigmas, nmersize, tol=1e-8):
    """
    Same as neighbors, but for 3D
    """
    xdiff = np.remainder(np.subtract.outer(xs, xs)+.5, 1)-.5
    ydiff = np.remainder(np.subtract.outer(ys, ys)+.5, 1)-.5
    zdiff = np.remainder(np.subtract.outer(zs, zs)+.5, 1)-.5
    sigmadists = np.add.outer(sigmas, sigmas)/2
    dists = np.sqrt((xdiff**2) + (ydiff**2) + (zdiff**2))
    #~ print(dists, sigmadists)
    #~ exit()
    matr = dists - sigmadists < tol
    bigN = len(matr)
    N = bigN // int(nmersize)
    n = int(nmersize)
    smallmatr = np.zeros((N,N))
    for i in range(N):
        li,hi = i*n, (i+1)*n
        for j in range(N):
            lj,hj = j*n, (j+1)*n
            ijcontacts = np.sum(matr[li:hi, lj:hj])
            #~ print('ijcontacts:', i, j, matr[li:hi, lj:hj], ijcontacts)
            #~ print(np.shape(dists), np.shape(sigmadists))
            #~ print('dists - sigmadists:', dists[li:hi, lj:hj] - sigmadists[li:hi, lj:hj])
            if i == j: ijcontacts = 1
            smallmatr[i,j] = ijcontacts
    return smallmatr

def largest_odd_factor(var_arr):
    """
    Function that computes the larges odd factors of an array of integers

    Parameters
    -----------------
    var_arr: numpy array
        Array of integers whose largest odd factors needs to be computed

    Returns
    ------------
    odd_d: numpy array
        Array of largest odd factors of each integer in var_arr
    """
    if var_arr.ndim == 1:
        odd_d = np.empty(np.shape(var_arr))
        odd_d[:] = np.NaN

        ind1 = np.where((np.remainder(var_arr, 2) != 0) | (var_arr == 0))[0]
        if np.size(ind1) != 0:
            odd_d[ind1] = var_arr[ind1]

        ind2 = np.where((np.remainder(var_arr, 2) == 0) & (var_arr != 0))[0]
        if np.size(ind2) != 0:
            odd_d[ind2] = largest_odd_factor(var_arr[ind2] / 2.0)
        return odd_d
    else:
        raise Exception('Wrong Input Type')

def edit_LOCPOT_file(filename):
# Removes a spurious set of lines between potential tables for major and minor spins so ASE will read it correctly
	from StringIO import StringIO
	f = open(filename)
	lines = f.read().splitlines()
	n_atoms = numpy.sum(numpy.genfromtxt(StringIO(lines[6]),dtype=(int)))
	n_points = numpy.product(numpy.genfromtxt(StringIO(lines[9+n_atoms]),dtype=(int)))
	f.close
	#get start line of the potential table for majority spin
	start_line = 10+n_atoms

	#get lastline of the potential table for majority spin
	if numpy.remainder(n_points,5) == 0:
		last_line = 9+n_atoms+n_points/5
	elif numpy.remainder(n_points,5) != 0:
		last_line = 9+n_atoms+n_points/5+1
	#print  "%d %s" % (last_line+1, lines[last_line+1])
	if numpy.remainder(n_atoms,5) == 0:
		n_atom_table_lines = n_atoms/5
	elif numpy.remainder(n_atoms,5) != 0:
		n_atom_table_lines = n_atoms/5+1
	last_atom_table_line = n_atom_table_lines+last_line
	#print  "%d %s" % (last_atom_table_line, lines[last_atom_table_line])
	del lines[last_line+1:last_atom_table_line+1]
	outfile = 'editted_LOCPOT_file'
	fout = open(outfile,'w')
	for line in lines:
		print>>fout, line
	fout.close

    def generic_test(self, iname):

        _log.info("Testing image output sizes for %s " % iname)
        inst = webbpsf.Instrument(iname)
        pxscale = inst.pixelscale
        fov_arcsec = 5.0

        PSF = inst.calcPSF(nlambda=1, fov_pixels = 100, oversample=1)
        self.assertEqual(PSF[0].data.shape[0], 100)

        PSF = inst.calcPSF(nlambda=1, fov_arcsec = fov_arcsec, oversample=1)
        fov_pix = int(np.round(fov_arcsec / pxscale))
        self.assertEqual(PSF[0].data.shape[0], fov_pix)

        inst.options['parity'] = 'odd'
        PSF = inst.calcPSF(nlambda=1, fov_arcsec = fov_arcsec, oversample=1)
        self.assertTrue( np.remainder(PSF[0].data.shape[0],2) == 1)

        inst.options['parity'] = 'even'
        PSF = inst.calcPSF(nlambda=1, fov_arcsec = fov_arcsec, oversample=1)
        self.assertTrue( np.remainder(PSF[0].data.shape[0],2) == 0)

        # odd array, even oversampling = even
        inst.options['parity'] = 'odd'
        PSF = inst.calcPSF(nlambda=1, fov_arcsec = fov_arcsec, oversample=2)
        self.assertTrue( np.remainder(PSF[0].data.shape[0],2) == 0)

        # odd array, odd oversampling = odd
        inst.options['parity'] = 'odd'
        PSF = inst.calcPSF(nlambda=1, fov_arcsec = fov_arcsec, oversample=3)
        self.assertTrue( np.remainder(PSF[0].data.shape[0],2) == 1)

def coord(x,y,field_of_view=100,window_size=(640,480)): # field of view in m
    "Convert world coordinates to pixel coordinates."
    fov_x = field_of_view
    fov_y = field_of_view/float(window_size[0])*float(window_size[1])
    wrapped_coord_x = np.remainder(x+fov_x/2.,fov_x)
    wrapped_coord_y = np.remainder(y+fov_y/2.,fov_y)
    return int(wrapped_coord_x/fov_x*window_size[0]), int(window_size[1]-wrapped_coord_y/fov_y*window_size[1])

    def positionPolygon(self, theta=None):
        """Calculates (R,Z) position at given theta angle by joining points
        by straight lines rather than a spline. This avoids the
        overshoots which can occur with splines.

        Parameters
        ----------
        theta : array_like, optional
            Theta locations to find R, Z at. If None (default), use the
            values of theta stored in the instance

        Returns
        -------
        R, Z : (ndarray, ndarray)
            Value of R, Z at each input theta point
        """
        if theta is None:
            return self.R, self.Z
        n = len(self.R)
        theta = np.remainder(theta, 2.*pi)
        dtheta = 2.*np.pi/n
        ind = np.trunc(theta/dtheta )
        rem = np.remainder(theta, dtheta)
        indp = (ind+1) % n
        return (rem*self.R[indp] + (1.-rem)*self.R[ind]), (rem*self.Z[indp] + (1.-rem)*self.Z[ind])

def generate_matching_df(template,color_lookup,top=20,sample=15):
    '''generate_matching_df (worst function name ever)
    this will generate a dataframe with x,y, corr, and png file path for images that most highly match each sampled pixel. The df gets parsed to json that plugs into d3 grid 
    '''
    # Now read in our actual image
    base = Image.open(template)
    width, height = base.size
    pixels = base.load()
    data = []

    count=0
    new_image = pandas.DataFrame(columns=["x","y","corr","png"])

    for x in range(width):
        for y in range(height):
            # And take only every [sample]th pixel
            if np.remainder(x,sample)==0 and np.remainder(y,sample)==0:
                cpixel = pixels[x, y]
                tmp = color_lookup.copy()
                tmp = (tmp-cpixel).abs().sum(axis=1)
                tmp.sort()
                png = choice(tmp.loc[tmp.index[0:top]].index.tolist(),1)[0]
                new_image.loc[count] = [x,y,0,png]
                count+=1

    new_image["x"] = [int(x) for x in (new_image["x"] / sample) * 10]
    new_image["y"] = [int(x) for x in (new_image["y"] / sample) * 10]
    
    return new_image

def compute_inp_params(lattice, sig_type):
    """
    tau and kmax necessary for possible integer quadruple combinations
    are computed

    Parameters
    ----------------
    lattice: Lattice class
        Attributes of the underlying lattice

    sig_type: {'common', 'specific'}

    Returns
    -----------
    tau: float
        tau is a rational number :math:`= \\frac{\\nu}{\\mu}`

    kmax: float
        kmax is an integer that depends on :math:`\\mu \\ , \\nu`
    """
    lat_params = lattice.lat_params
    cryst_ptgrp = proper_ptgrp(lattice.cryst_ptgrp)

    if cryst_ptgrp == 'D3':
        c_alpha = np.cos(lat_params['alpha'])
        tau = c_alpha / (1 + 2 * c_alpha)
        if sig_type == 'specific':
            [nu, mu] = int_man.rat(tau)
            rho = mu - 3 * nu
            kmax = 4 * mu * rho
        elif sig_type == 'common':
            kmax = []

    if cryst_ptgrp == 'D4':
        tau = (lat_params['a'] ** 2) / (lat_params['c'] ** 2)
        if sig_type == 'specific':
            [nu, mu] = int_man.rat(tau)
            kmax = 4 * mu * nu
        if sig_type == 'common':
            kmax = []

    if cryst_ptgrp == 'D6':
        tau = (lat_params['a'] ** 2) / (lat_params['c'] ** 2)
        if sig_type == 'specific':
            [nu, mu] = int_man.rat(tau)
            if np.remainder(nu, 2) == 0:
                if np.remainder(nu, 4) == 0:
                    kmax = 3 * mu * nu
                else:
                    kmax = 6 * mu * nu
            else:
                kmax = 12 * mu * nu
        if sig_type == 'common':
            kmax = []

    if cryst_ptgrp == 'O':
        tau = 1
        kmax = []

    return tau, kmax

	def create_tag_noleap(self, time, timeunits):
		''' create a datetime object from reference date and ocean_time (noleap version)'''
		# first we need to figure out how many seconds are ellaped between ref_date
		# and start date of the run
		delta_type  = timeunits.split()[0]
		date_string = timeunits.split()[2]
		time_string = timeunits.split()[3]
		ref_string = date_string + ' ' + time_string
		fmt = '%Y-%m-%d %H:%M:%S'
		dateref_dstart = dt.datetime.strptime(ref_string,fmt)

		if delta_type == 'seconds':
			seconds_from_init = float(time)
		elif delta_type == 'days':
			seconds_from_init = float(time) * 86400.

		nyear  = int(np.floor(seconds_from_init / 365 / 86400))
		rm     = np.remainder(seconds_from_init,365*86400)
		ndays  = int(np.floor(rm / 86400))
		rm2    = np.remainder(rm,86400)
		nhours = int(np.floor(rm2 / 3600))
		rm3    = np.remainder(rm2,3600)
		nmin   = int(np.floor(rm3 / 60))
		nsec   = int(np.remainder(rm3,60))

		# pick a year we are sure is not a leap year
		fakeref  = dt.datetime(1901,1,1,0,0)
		fakedate = fakeref + dt.timedelta(days=ndays)
		month    = fakedate.month
		day      = fakedate.day

		tag=dt.datetime(nyear + dateref_dstart.year,month, day, nhours, nmin, nsec)
		return tag

def nearest_odd(N):
    """
    Get the nearest odd number for each value of N.

    :param N: int / sequence of ints
    :return: int / sequence of ints
    :Example:
    >>> from __future__ import print_function
    >>> print(nearest_odd(range(1, 11)))
    [  1.   3.   3.   5.   5.   7.   7.   9.   9.  11.]
    >>> nearest_odd(0)
    1
    >>> nearest_odd(3)
    3.0
    """
    if hasattr(N, "__iter__"):
        N = np.array(N)
        y = np.floor(N)
        y[np.remainder(y, 2) == 0] = np.ceil(N[np.remainder(y, 2) == 0])
        y[np.remainder(y, 2) == 0] += 1
        return y
    if N % 2 == 0:
        return N + 1
    elif np.floor(N) % 2 == 0:
        return np.ceil(N)
    elif np.floor(N) % 2 != 0:
        return np.floor(N)
    return N

def plot_mwd(RA,Dec,org=0,title='Mollweide projection', projection='mollweide'):
    ''' RA, Dec are arrays of the same length.
    RA takes values in [0,360), Dec in [-90,90],
    which represent angles in degrees.
    org is the origin of the plot, 0 or a multiple of 30 degrees in [0,360).
    title is the title of the figure.
    projection is the kind of projection: 'mollweide', 'aitoff', 'hammer', 'lambert'
    '''
    x = N.remainder(RA+360-org,360) # shift RA values
    ind = x>180
    x[ind] -=360    # scale conversion to [-180, 180]
    x=-x    # reverse the scale: East to the left
    tick_labels = N.array([150, 120, 90, 60, 30, 0, 330, 300, 270, 240, 210])
    tick_labels = N.remainder(tick_labels+360+org,360)
    fig = P.figure(figsize=(10, 5))
    ax = fig.add_subplot(111, projection=projection,)# axisbg ='LightCyan')
    ax.scatter(N.radians(x),N.radians(Dec))  # convert degrees to radians
    ax.set_xticklabels(tick_labels)     # we add the scale on the x axis
    ax.set_title(title)
    ax.title.set_fontsize(15)
    ax.set_xlabel("RA")
    ax.xaxis.label.set_fontsize(12)
    ax.set_ylabel("Dec")
    ax.yaxis.label.set_fontsize(12)
    ax.grid(True)

def idx2xyz(idx, xl, yl, zl):
    """Transform a list of indices of 1D array into coordinates of a 3D volume of certain sizes.
    
    Parameters
    ----------
    idx : np.ndarray
        1D array to be converted. An increment of ``idx``
        corresponds to a an increment of x. When reaching ``xl``, x is reset and 
        y is incremented of one. When reaching ``yl``, x and y are reset and z is
        incremented.
    
    xl, yl, zl : int 
        Sizes for 3D volume.
    
    Returns
    -------
    list 
        List of 3 ``np.ndarray`` objects (for x, y and z), containing coordinate value.

    """
    
    z = np.floor(idx / (xl*yl))
    r = np.remainder(idx, xl*yl)
    y = np.floor(r / xl)
    x = np.remainder(r, xl)
    return x, y, z

    def __init__(self, x,y,z,sig):
        """Requires periodic BC, box length 1"""
        self.points = np.remainder(array((x,y,z), dtype=float), 1)
        self.sigmas = array(sig, dtype=float)
        self.N = N = len(sig)
        pts = array(list(self.points) + [sig])
        pts1 = np.concatenate([pts.T +(n,m,p,0) for n in [0,-1] for m in [0,-1] for p in [0,-1]], axis=0)
        self.allpoints = pts2 = np.remainder(pts1[:,:3] + 0.5,2)-1
        self.allsigmas = s2 = pts1[:,3]
        d = self.delaunay = Delaunay(pts2)
        
        d.simplices
        triangs = [(d.simplices[:,0], d.simplices[:,1], d.simplices[:,2]),
                   (d.simplices[:,0], d.simplices[:,1], d.simplices[:,3]),
                   (d.simplices[:,0], d.simplices[:,2], d.simplices[:,3]),
                   (d.simplices[:,1], d.simplices[:,2], d.simplices[:,3])]

        triangs = np.concatenate(triangs,axis=1).T
        #print(shape(array(triangs)))

        triangs.sort(1)
        triangs2 = triangs[triangs[:,0] < self.N]
        #print(shape(array(triangs2)))

        trirem = np.remainder(triangs2,N)
        #trirem.sort(1)
        self.triangles = triangs2[unique_rows(trirem)]

def _detect(alpha, beta_1, beta_2, phi_range, tau, As, Au, tau_shift):
	omega_p = np.pi/(2*T)
	phi_p   = omega_p * tau

	tau_1_  = np.remainder(tau, 2*T)
	tau_2_  = np.remainder(tau+tau_shift, 2*T)
	k_tau_1 = int(np.floor(tau / (2*T)))
	k_tau_2 = int(np.floor((tau+tau_shift) / (2*T)))

	bkn_1, bk_1 = shift_indices(beta_1, k_tau_1)
	bkn_2, bk_2 = shift_indices(beta_2, k_tau_2)

	arg1 = np.cos(phi_p) * (tau_1_ * bkn_1 + (2*T - tau_1_) * bk_1)
	arg2 = ((2*T) / np.pi) * np.sin(phi_p) * (bkn_1 - bk_1)
	arg3 = np.sin(phi_p) * (tau_2_ * bkn_2 + (2*T - tau_2_) * bk_2)
	arg4 = ((2*T) / np.pi) * np.cos(phi_p) * (bkn_2 - bk_2)

	PHI_C, ALPHA = np.meshgrid(phi_range, alpha)
	ARG12 = np.meshgrid(phi_range, (arg1 - arg2))[1]
	ARG34 = np.meshgrid(phi_range, (arg3 + arg4))[1]

	# TODO: fix Au in partial overlap! One strong chip can dominate the complete correlation of a symbol ...
	result = (T/2) * ALPHA * As + (Au/4) * (np.cos(PHI_C) * ARG12 - np.sin(PHI_C) * ARG34)

	# modified for Au only transmissions
	if As > 0:
		return result / ((T/2) * As * Au)
	else:
		return result / ((T/2) * Au)

def get_indices_peak(sequence, ind_max, threshold=0.8):
    """returns the indices of the peak around ind_max,
    with values down to  ``threshold * sequence[ind_max]``
    """
    thres_value = sequence[ind_max] * threshold
    # considering sequence is circular, we put ind_max to the center,
    # and find the bounds:
    lenSeq = sequence.size
    ##midSeq = lenSeq / 2
    ##indices = np.remainder(np.arange(lenSeq) - midSeq + ind_max, lenSeq)
    ##
    #newseq = sequence[]
    indSup = ind_max
    indInf = ind_max
    while sequence[indSup]>thres_value:
        indSup += 1
        indSup = np.remainder(indSup, lenSeq)
    while sequence[indInf]>thres_value:
        indInf -= 1
        indInf = np.remainder(indInf, lenSeq)
        
    indices = np.zeros(lenSeq, dtype=bool)
    if indInf < indSup:
        indices[(indInf+1):indSup] = True
    else:
        indices[:indSup] = True
        indices[(indInf+1):] = True
    
    return indices

def createDistanceList(mergedDict, direction='ra'):

    cov = []
    unit = []
    X = 4
    indexDict=mergedDict[direction]
    if direction=='ra':
        for RA_grid, starList in indexDict.items():
            if np.remainder(RA_grid, X)==0:
                for i in range(len(starList)-1):
                    for j in range(i+1, len(starList)):
                        A = starList[i]
                        B = starList[j]
                        cov.append(A.e1*B.e1 + A.e2*B.e2)
                        unit.append(int(abs(B.DEC_grid - A.DEC_grid)))
    if direction=='dec':
        for DEC_grid, starList in indexDict.items():
            if np.remainder(DEC_grid, X)==0:
                for i in range(len(starList)-1):
                    for j in range(i+1, len(starList)):
                        A = starList[i]
                        B = starList[j]
                        cov.append(A.e1*B.e1 + A.e2*B.e2)
                        unit.append(int(abs(B.RA_grid - A.RA_grid)))

    x=list(set(unit))
    mean = [0]*len(x)
    d = [[] for i in range(len(x))]
    for i in range(len(unit)):
        d[unit[i]-1].append(cov[i])

    for j in range(len(d)):
        mean[j] = np.mean(d[j])

    return x, mean

def TimestampMerge(min_node,tot,stamp_ID):
    """Merge timestamps files across nodes and write merged list to file"""
    masterlist=[]
    diff=[]
    n=min_node
    max_node = min_node + tot
    while n < max_node:
        mastertime=[]
        j=0
        while j<5:
            fname='node{0}/timestamp{1}.{2}.dat'.format(n,stamp_ID,j)
            try:
                times=man.LoadData(fname)
                year=man.IterativeStrAppend(times,0)
                month=man.IterativeStrAppend(times,1)
                day=man.IterativeStrAppend(times,2)
                hour=man.IterativeStrAppend(times,3)
                minute=man.IterativeStrAppend(times,4)

                seconds=[]
                for i in range(len(times)):
                    a=man.IterativeIntAppend(times,5)
                    b=man.IterativeFloatAppend(times,6)
                    point=a[i]+b[i]
                    if np.remainder(i,100) == 0:
                        print 'Done {0} of {1} seconds'.format(i,len(times))
                    seconds.append(point)

                time=[]
                for i in range(len(times)):
                    point="{0}-{1}-{2} {3}:{4}:{5}".format(year[i],month[i],day[i],hour[i],minute[i],seconds[i])
                    if np.remainder(i,100) == 0:
                        print 'Done {0} of {1} stamps'.format(i,len(times))
                    time.append(point)

                t=Time(time, format='iso',scale='utc')
                mjd=t.mjd
                mjd.sort()
                diff.append(mjd)
                r=Time(mjd,format='mjd',scale='utc')
                iso=r.iso
                for i in range(len(iso)):
                    point=[iso[i],j-1]
                    mastertime.append(point)
                    masterlist.append(point)
                print "Done node{0}/timestamp{1}.{2}.dat".format(n,stamp_ID,j)
                j+=1 
            except IOError:
                print "Missing node{0}/timestamp{1}.{2}.dat".format(n,stamp_ID,j)
                j+=1
                pass

        mastertime.sort()
        name='node{0}/MergedTimeStamp{1}.dat'.format(n,stamp_ID)
        man.WriteFileCols(mastertime,name)
        print 'Done master time stamp node{0}'.format(n)
        n+=1
    masterlist.sort()
    name='MasterTimeStamp{0}.dat'.format(stamp_ID)

def matrix_mod_exp(n,T,I,mod):
  Accum = I
  while n:
    if n % 2:
      Accum = np.remainder(T.dot(Accum),mod)
    T = np.remainder(T.dot(T),mod)
    n /= 2
  return Accum