def gallery_gray_patches(W,show=False, rescale=False):
    """Create a gallery of image patches from <W>, with
    grayscale patches aligned along columns"""

    n_vis, n_feats = W.shape;

    n_pix = np.sqrt(n_vis)
    n_rows = np.floor(np.sqrt(n_feats))
    n_cols = np.ceil(n_feats / n_rows)    
    border_pix = 1;

    # INITIALIZE GALLERY CONTAINER    
    im_gallery = np.nan*np.ones((border_pix+n_rows*(border_pix+n_pix),
                   border_pix+n_cols*(border_pix+n_pix)))

    for iw in xrange(n_feats):
        # RESCALE EACH IMAGE
        W_tmp = W[:,iw].copy()
        if rescale:
            W_tmp = (W_tmp - W_tmp.mean())/np.max(np.abs(W_tmp)); 

        W_tmp = W_tmp.reshape(n_pix, n_pix)

        # FANCY INDEXING INTO IMAGE GALLERY
        im_gallery[border_pix + np.floor(iw/n_cols)*(border_pix+n_pix): 
                        border_pix + (1 + np.floor(iw/n_cols))*(border_pix+n_pix) - border_pix, 
                  border_pix + np.mod(iw,n_cols)*(border_pix+n_pix): 
                        border_pix + (1 + np.mod(iw,n_cols))*(border_pix+n_pix) - border_pix] = W_tmp
    if show:
        plt.imshow(im_gallery,interpolation='none')
        plt.axis("image")
        plt.axis("off")

    return im_gallery

    def __call__(self, x, pos=None):
        '''Return the time format as pos'''

        _, dmax = self.axis.get_data_interval()
        vmin, vmax = self.axis.get_view_interval()

        # In lag-time axes, anything greater than dmax / 2 is negative time
        if self.lag and x >= dmax * 0.5:
            # In lag mode, don't tick past the limits of the data
            if x > dmax:
                return ''
            value = np.abs(x - dmax)
            # Do we need to tweak vmin/vmax here?
            sign = '-'
        else:
            value = x
            sign = ''

        if vmax - vmin > 3600:
            s = '{:d}:{:02d}:{:02d}'.format(int(value / 3600.0),
                                            int(np.mod(value / 60.0, 60)),
                                            int(np.mod(value, 60)))
        elif vmax - vmin > 60:
            s = '{:d}:{:02d}'.format(int(value / 60.0),
                                     int(np.mod(value, 60)))
        else:
            s = '{:.2g}'.format(value)

        return '{:s}{:s}'.format(sign, s)

    def _find_tails(self, mag, rounding=True, half=5, full=10, flag=50):
        '''
        Find how many of each of the tail pieces is necessary.  Flag
        specifies the increment for a flag, barb for a full barb, and half for
        half a barb. Mag should be the magnitude of a vector (ie. >= 0).

        This returns a tuple of:

            (*number of flags*, *number of barbs*, *half_flag*, *empty_flag*)

        *half_flag* is a boolean whether half of a barb is needed,
        since there should only ever be one half on a given
        barb. *empty_flag* flag is an array of flags to easily tell if
        a barb is empty (too low to plot any barbs/flags.
        '''

        #If rounding, round to the nearest multiple of half, the smallest
        #increment
        if rounding:
            mag = half * (mag / half + 0.5).astype(np.int)

        num_flags = np.floor(mag / flag).astype(np.int)
        mag = np.mod(mag, flag)

        num_barb = np.floor(mag / full).astype(np.int)
        mag = np.mod(mag, full)

        half_flag = mag >= half
        empty_flag = ~(half_flag | (num_flags > 0) | (num_barb > 0))

        return num_flags, num_barb, half_flag, empty_flag

def boundaryGaussian2(N,P,
                      A00,A01,alphaX00,alphaX01,x00,x01,
                      A10,A11,alphaX10,alphaX11,x10,x11):
    #
    # f0(x) = A00exp(-alphaX00(x-x00)^2) + A01exp(-alphaX01(x-x01)^2)
    # f1(x) = A10exp(-alphaX10(x-x10)^2) + A11exp(-alphaX11(x-x11)^2)
    #

    x00 = np.mod(x00,1.)
    x01 = np.mod(x01,1.)
    x10 = np.mod(x10,1.)
    x11 = np.mod(x11,1.)

    # Defines f0 and f1
    X  = np.linspace( 0.0 , 1.0 , N + 1 )

    f0 = ( A00 * np.exp( -alphaX00 * np.power( X - x00 , 2 ) ) +
           A01 * np.exp( -alphaX01 * np.power( X - x01 , 2 ) ) )

    f1 = ( A10 * np.exp( -alphaX10 * np.power( X - x10 , 2 ) ) +
           A11 * np.exp( -alphaX11 * np.power( X - x11 , 2 ) ) )

    temporalBoundaries = grid.TemporalBoundaries( N , P , f0 , f1 )
    spatialBoundaries  = grid.SpatialBoundaries( N , P )

    return grid.Boundaries( N , P ,
                            temporalBoundaries, spatialBoundaries )

def j2000tob1950(ra, dec):
    """
    Convert J2000 to B1950 coordinates.

    This routine was derived by taking the inverse of the b1950toj2000 routine
    """

    # Convert to radians
    ra = np.radians(ra)
    dec = np.radians(dec)

    # Convert RA, Dec to rectangular coordinates
    x = np.cos(ra) * np.cos(dec)
    y = np.sin(ra) * np.cos(dec)
    z = np.sin(dec)

    # Apply the precession matrix
    x2 = P2[0, 0] * x + P2[1, 0] * y + P2[2, 0] * z
    y2 = P2[0, 1] * x + P2[1, 1] * y + P2[2, 1] * z
    z2 = P2[0, 2] * x + P2[1, 2] * y + P2[2, 2] * z

    # Convert the new rectangular coordinates back to RA, Dec
    ra = np.arctan2(y2, x2)
    dec = np.arcsin(z2)

    # Convert to degrees
    ra = np.degrees(ra)
    dec = np.degrees(dec)

    # Make sure ra is between 0. and 360.
    ra = np.mod(ra, 360.0)
    dec = np.mod(dec + 90.0, 180.0) - 90.0

    return ra, dec

def toKepler(u, which = 'Pueyo', mass = 1, referenceTime = None):
    """
    """
    if which == 'Pueyo':
        res = np.zeros(6)
        res[1] = u[1]
        res[5] = u[5]
        
        res[0] = semimajoraxis(math.exp(u[0]), starMass = mass)
        res[2] = math.degrees(math.acos(u[2]))
        res[3] = np.mod((u[3]-u[4])*0.5,360)
        res[4] = np.mod((u[3]+u[4])*0.5,360)
        return res
    elif which == 'alternative':
        res = np.zeros(6)
        res[1] = u[1]
        res[5] = u[5]
        
        res[0] = semimajoraxis(math.exp(u[0]), starMass = mass)
        res[2] = math.degrees(math.acos(u[2]))
        res[3] = u[3]
        res[4] = u[4]
        return res        
    elif which == 'Chauvin':
        stat = StatisticsMCMC()
        res = stat.xFROMu(u,referenceTime,mass)    
        return res
    
    return None

 def getBits(self,cell):
         zero=[-self.markerArea[i]/2. for i in [0,1]]
         bitx=[int(i) for i in bin(int(cell[0]))[::-1][:-2]]
         bity=[int(i) for i in bin(int(cell[1]))[::-1][:-2]]
         s0=int(np.log2(self.cellsPerBlock[0]*self.noBlocks[0]))
         s1=int(np.log2(self.cellsPerBlock[1]*self.noBlocks[1]))
         for i in range(s0-len(bitx)):
             bitx.append(0)
         for i in range(s1-len(bity)):
             bity.append(0)
         tx=np.zeros(s0,dtype=np.bool)
         ty=np.zeros(s1,dtype=np.bool)
         px=np.empty((s0,2))
         py=np.empty((s1,2))
         for i,b in enumerate(bitx):
             x=zero[0]+mod(i+1,self.noBitsX)*self.bitDistance
             y=zero[1]+((i+1)/self.noBitsY)*self.bitDistance
             px[i]=(x,y)
             tx[i]=b
         for i,b in enumerate(bity):
             x=zero[0]+(self.noBitsX-mod(i+1,self.noBitsX)-1)*self.bitDistance
             y=zero[1]+(self.noBitsY-(i+1)/self.noBitsY-1)*self.bitDistance
             py[i]=(x,y)
             ty[i]=b
         return px,py,tx,ty

 def getBoxFilter(self,
         time, point_coords,
         data_set = 'isotropic1024coarse',
         make_modulo = False,
         field = 'velocity',
         filter_width = 7*2*np.pi / 1024):
     if not self.connection_on:
         print('you didn\'t connect to the database')
         sys.exit()
     if not (point_coords.shape[-1] == 3):
         print ('wrong number of values for coordinates in getBoxFilter')
         sys.exit()
         return None
     if not (point_coords.dtype == np.float32):
         print 'point coordinates in getBoxFilter must be floats. stopping.'
         sys.exit()
         return None
     npoints = point_coords.shape[0]
     for i in range(1, len(point_coords.shape)-1):
         npoints *= point_coords.shape[i]
     if make_modulo:
         pcoords = np.zeros(point_coords.shape, np.float64)
         pcoords[:] = point_coords
         np.mod(pcoords, 2*np.pi, point_coords)
     result_array = point_coords.copy()
     self.lib.getBoxFilter(self.authToken,
              ctypes.c_char_p(data_set),
              ctypes.c_char_p(field),
              ctypes.c_float(time),
              ctypes.c_float(filter_width),
              ctypes.c_int(npoints),
              point_coords.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))),
              result_array.ctypes.data_as(ctypes.POINTER(ctypes.POINTER(ctypes.c_float))))
     return result_array

    def draw(self, dt):
        if self._mixer.is_onset():
            self.onset_speed_boost = self.parameter('onset-speed-boost').get()

        self.center_offset_angle += dt * self.parameter('center-speed').get() * self.onset_speed_boost
        self.hue_inner += dt * self.parameter('hue-speed').get() * self.onset_speed_boost
        self.wave_offset += dt * self.parameter('wave-speed').get() * self.onset_speed_boost
        self.color_offset += dt * self.parameter('speed').get() * self.onset_speed_boost

        self.onset_speed_boost = max(1, self.onset_speed_boost - self.parameter('onset-speed-decay').get())

        wave_hue_period = 2 * math.pi * self.parameter('wave-hue-period').get()
        wave_hue_width = self.parameter('wave-hue-width').get()
        radius_hue_width = self.parameter('radius-hue-width').get()
        angle_hue_width = self.parameter('angle-hue-width').get()

        cx, cy = self.scene().center_point()
        self.locations = np.asarray(self.scene().get_all_pixel_locations())
        x,y = self.locations.T
        x -= cx + math.cos(self.center_offset_angle) * self.parameter('center-distance').get()
        y -= cy + math.sin(self.center_offset_angle) * self.parameter('center-distance').get()
        self.pixel_distances = np.sqrt(np.square(x) + np.square(y))
        self.pixel_angles = np.arctan2(y, x) / (2.0 * math.pi)
        self.pixel_distances /= max(self.pixel_distances)

        angles = np.mod(1.0 + self.pixel_angles + np.sin(self.wave_offset + self.pixel_distances * wave_hue_period) * wave_hue_width, 1.0)
        hues = self.color_offset + (radius_hue_width * self.pixel_distances) + (2 * np.abs(angles - 0.5) * angle_hue_width)
        hues = np.int_(np.mod(hues, 1.0) * self._fader_steps)
        colors = self._fader.color_cache[hues]
        colors = colors.T
        colors[0] = np.mod(colors[0] + self.hue_inner, 1.0)
        colors = colors.T

        self._pixel_buffer = colors

def go_back(x,y,a=40,b=50,R_0 = 90,l=10,m=3,aa=0.0,q0=3.0,eps=.07):

    hit_divert = (x>b)
    inCORE = x<b
    stopevolve = hit_divert
    
    #eps = .2
    C = ((2*m*l*a**2)/(R_0*q0*b**2))*eps
    
    def func(y_out):
        return (-y + y_out - C*(x/b)**(m-2) *np.cos(m*y_out))**2
    
    def func2(y_out):
        return (-y_old + y_out + (2.0*np.pi/q) + aa*np.cos(y_out))**2
    
    y_old = copy(y)
    y_old = (newton_krylov(func,y))
    y_old = np.mod(y_old,2.0*np.pi)

    x_old = x + (m*b*C)/(m-1)*(x/b)**(m-1) *np.sin(m*y_old)

    q = q0*(x_old/a)**2

    y_old2 = copy(y_old)
    y_old2 = (newton_krylov(func2,y_old))

    #y_old2 = y_old - 2*np.pi/q #- aa*np.cos(
    y_old2 = np.mod(y_old2,2.0*np.pi)
    x_old2 = x_old*(1.0 -aa*np.sin(y_old2))



    return x_old2,y_old2

def L_sys_curve(function, level):
    axiom, rules, angleL, angleR, angle0 = function()
    angleL, angleR, angle0 = np.array([angleL, angleR, angle0]) * np.pi / 180.
    gen = generation(axiom, rules, level)
    print(gen)
    deux_pi, radius = 2 * np.pi, 1.0
    x0, y0, angle, stack = 0., 0., angle0, []
    x, y = [x0], [y0]
    plt.clf()
    for c in gen:
        if c == '[':
            stack.append((x0, y0, angle))
        elif c == ']':
            plt.plot(x, y, 'g')
            x0, y0, angle = stack.pop()
            x, y = [x0], [y0]
        elif c == '+':
            angle = np.mod(angle + angleR, deux_pi)
        elif c == '-':
            angle = np.mod(angle + angleL, deux_pi)
        else:
            if c == 'f':  # jump
                plt.plot(x, y, 'b')
                x, y = [], []
            x0 = x0 + radius * np.cos(angle)
            y0 = y0 + radius * np.sin(angle)
            x.append(x0)
            y.append(y0)
    plt.axis('off')  # ,axisbg=(1, 1, 1))
    plt.plot(x, y, 'g')
    plt.show()
    return

def go_forward(x,y,a=40,b=50,R_0 = 90,l=10,m=3,aa=0.0,q0=3.0):

    hit_divert = (x>b)
    inCORE = x<b
    stopevolve = hit_divert
    
    eps = .2
    C = ((2*m*l*a**2)/(R_0*q0*b**2))*eps

    x_new = x/(1-aa*np.sin(y))
    q = q0*(x_new/a)**2
    y_new =  (y+ 2*np.pi/q + aa*np.cos(y))
    y_new = np.mod(y_new,2*np.pi)

    def func(x_out):
        return (-x_new + x_out +(m*b*C)/(m-1)*(x_out/b)**(m-1) *np.sin(m*y_new))**2
    
    x_new2 = (newton_krylov(func,x_new,method='gmres',maxiter=50))
    y_new2 = (y_new - C*(x_new2/b)**(m-2) * np.cos(m*y_new))
                                
    #print 'xchange:', x_new2/x

    x_new = x_new2
    y_new = np.mod(y_new2,2*np.pi)

    
    return x_new,y_new

def PsplineXY(x,a):
    """
    periodic spline (P=1)

    a = [x0,y0, x1,y1, x2,y2...]
    x in [0,1]
    """
    coef = numpy.array(a)
    xp = numpy.zeros(3*len(coef)/2)
    yp = numpy.zeros(3*len(coef)/2)
    x0 = numpy.mod(coef[::2], 1.0)
    s = numpy.array(x0).argsort()
    xp[0:len(coef)//2]      = numpy.mod(x0[s], 1.0)-1
    xp[len(coef)//2:len(coef)] = xp[0:len(coef)//2]+1
    xp[-len(coef)//2:]      = xp[0:len(coef)//2]+2
    yp[0:len(coef)//2]      = coef[1::2][s]
    yp[len(coef)//2:len(coef)] = yp[0:len(coef)//2]
    yp[-len(coef)//2:]      = yp[0:len(coef)//2]
        
    xx = numpy.array(x).flatten()
    res = interp1d(xp, yp, kind='quadratic', \
                   bounds_error=False, fill_value=0.0)\
                   (numpy.mod(xx, 1))
    res = res.reshape(numpy.array(x).shape)
    return res

def calc_geoinc(trace,metric=True):
    """docstring for calc_geoinc"""
    stla=trace.stats.sac.stla
    stlo=trace.stats.sac.stlo
    stdp=trace.stats.sac.stdp
    evla=trace.stats.sac.evla
    evlo=trace.stats.sac.evlo
    evdp=trace.stats.sac.evdp
    
    
    if metric:
        baz=np.rad2deg(np.arctan2((evlo-stlo),(evla-stla)))
        EpiDist = np.sqrt((evlo-stlo)**2. +  (evla-stla)**2.)
        inc = np.rad2deg(np.arctan(EpiDist/ (evdp-stdp)))
        
        HypoDist = np.sqrt((evdp-stdp)**2. + EpiDist**2)
         
    if baz<0.:
        baz=baz+360.
        
    
    azi=np.mod(baz+180.0,360.0)
    inc=np.mod(inc+180.0,180.)

    
    return azi,inc,baz,HypoDist,EpiDist,stdp

    def compute_LDPs(self,ln,RAT):
        """compute edge LDP

        Parameters
        ----------

        n1      : float/string
            node ID
        n2      : float/string
            node ID
        RAT     : string
            A specific RAT which exist in the network ( if not , raises an error)
        value    : list : [LDP value , LDP standard deviation] 
        method    : ElectroMagnetic Solver method ( 'direct', 'Multiwall', 'PyRay'


        """
        p=nx.get_node_attributes(self.SubNet[RAT],'p')
        epwr=nx.get_node_attributes(self.SubNet[RAT],'epwr')
        sens=nx.get_node_attributes(self.SubNet[RAT],'sens')
        e=self.link[RAT]#self.SubNet[RAT].edges()
        re=self.relink[RAT] # reverse link aka other direction of link
        lp,lt, d, v= self.EMS.solve(p,e,'all',RAT,epwr,sens)
        lD=[{'Pr':lp[i],'TOA':lt[np.mod(i,len(e))] ,'d':d[np.mod(i,len(e))],'vis':v[i]} for i in range(len(d))]
        self.update_LDPs(iter(e+re),RAT,lD)

    def champ(self):
        if self.structure: N_lame = self.N_lame-self.struct_N
        else: N_lame = self.N_lame

        force = np.zeros_like(self.lames[2, :N_lame])
        damp_min = 0.8
        damp_tau = 15.
        damp = lambda t: damp_min #+ (1.-damp_min)*np.exp(-np.abs(np.mod(t+self.period/2, self.period)-self.period/2)/damp_tau)

        smooth = lambda t: 1.-np.exp(-np.abs(np.mod(t+self.period/2, self.period)-self.period/2)**2/damp_tau**2)
        on_off = lambda t, freq: (np.sin(2*np.pi*t/self.period*freq) > 0.)

        noise = lambda t: .4 * smooth(t)
        np.random.seed(12345)
        random_timing = np.random.rand(N_lame)
        #print(np.mod(random_timing + self.t/self.period, 1))
        struct_angles = np.random.permutation( np.hstack((self.struct_angles, -np.array(self.struct_angles)))) * np.pi / 180
        #print(struct_angles, (np.mod(random_timing + self.t/self.period, 1)*len(struct_angles)).astype(np.int)) 
        angle_desired = np.zeros(N_lame)
        for i, idx in enumerate((np.mod(random_timing + self.t/self.period, 1)*len(struct_angles)).astype(np.int)):
            angle_desired[i] = struct_angles[idx]
        
        force -= 20 * (np.mod(self.lames[2, :N_lame] - angle_desired +np.pi/2, np.pi) - np.pi/2 ) * smooth(self.t)
        force -= 80 * (np.mod(self.lames[2, :N_lame] + np.pi/2, np.pi) - np.pi/2) * (1- smooth(self.t) )
        force += noise(self.t)*np.pi*np.random.randn(N_lame)
        force -= damp(self.t) * self.lames[3, :N_lame]/self.dt
        force = .02 * 100 * np.tanh(force/100)
        return force    

def toSynthetic(theta,which = 'Pueyo', mass = 1, referenceTime = None):
    """
    """
    if which == 'Pueyo':
        res = np.zeros(6)
        temp = theta
        res[0] = math.log(period(temp[0],starMass = mass))
        res[2] = math.cos(math.radians(temp[2]))
        res[3] = np.mod(temp[4]+temp[3],360)
        res[4] = np.mod(temp[4]-temp[3],360)        
        res[1] = temp[1]
        res[5] = temp[5]
    elif which == 'alternative':
        res = np.zeros(6)
        temp = theta
        res[0] = math.log(period(temp[0],starMass = mass))
        res[2] = math.cos(math.radians(temp[2]))
        res[3] = temp[3]
        res[4] = temp[4]        
        res[1] = temp[1]
        res[5] = temp[5]        
    elif which == 'Chauvin':
        stat = StatisticsMCMC()
        res = stat.uFROMx(theta,referenceTime,mass)

    return res

def hist_lastAxSub(data, row, col, plot_title = '', chBool = None, tag = 'CH'):
    """imshow basically. takes 3D data, plots first 2 dim, separates by 3rd dim"""
    numCh = data.shape[-1]
    print "numCh:", numCh
    if numCh > row*col:
        numCh = row*col
    f, axarr = plt.subplots(row, col, sharex='col', sharey='row')
    if row == 1 or col == 1:
        for v in range(numCh):
            axarr[v].hist(data[:,v])
            axarr[v].set_title(tag + str(v))
            try:
                if chBool[v] == False:
                    axarr[v].set_title('(' + tag.lower() + + str(v) + ')')
            except:
                pass
        plt.suptitle(plot_title)
    else:
        for v in range(numCh):
            axarr[v/col, np.mod(v,col)].hist(data[:,v])
            axarr[v/col, np.mod(v,col)].set_title(tag + str(v))
            axarr[v/col, np.mod(v,col)].locator_params(nbins=2)
            try:
                if chBool[v] == False:
                    axarr[v/col, np.mod(v,col)].set_title('(' + tag.lower() + + str(v) + ')')
            except:
                pass                         
        plt.suptitle(plot_title)

  def movementCompute(self, displacement, noiseFactor = 0):
    """
    Shift the current active cells by a vector.

    @param displacement (pair of floats)
    A translation vector [di, dj].
    """

    if noiseFactor != 0:
      displacement = copy.deepcopy(displacement)
      xnoise = np.random.normal(0, noiseFactor)
      ynoise = np.random.normal(0, noiseFactor)
      displacement[0] += xnoise
      displacement[1] += ynoise


    # Calculate delta in the module's coordinates.
    phaseDisplacement = (np.matmul(self.rotationMatrix, displacement) *
                         self.phasesPerUnitDistance)

    # Shift the active coordinates.
    np.add(self.activePhases, phaseDisplacement, out=self.activePhases)

    # In Python, (x % 1.0) can return 1.0 because of floating point goofiness.
    # Generally this doesn't cause problems, it's just confusing when you're
    # debugging.
    np.round(self.activePhases, decimals=9, out=self.activePhases)
    np.mod(self.activePhases, 1.0, out=self.activePhases)

    self._computeActiveCells()
    self.phaseDisplacement = phaseDisplacement

def magic(n):
    ix = np.arange(n)+1
    J, I = np.meshgrid(ix,ix)
    A = np.mod(I+J-(n+3)/2,n)
    B = np.mod(I+2*J-2,n)
    M = n*A + B + 1
    return M