def find_surface(cube):
    """
    Return the `cube` index for the surface layer of for any model grid
    (rgrid, ugrid, sgrid), and any non-dimensional coordinate.

    TODO: Fold this into `find_layer()`

    """
    z = z_coord(cube)
    if not z:
        msg = "Cannot find the surface for cube {!r}".format
        raise ValueError(msg(cube))
    else:
        if np.argmin(z.shape) == 0 and z.ndim == 2:
            points = z[:, 0].points
        elif np.argmin(z.shape) == 1 and z.ndim == 2:
            points = z[0, :].points
        else:
            points = z.points
        positive = z.attributes.get('positive', None)
        if positive == 'up':
            idx = np.unique(points.argmax(axis=0))[0]
        else:
            idx = np.unique(points.argmin(axis=0))[0]
        return idx

def main():
    try:
        prog = sys.argv[0]
#       a = float(sys.argv[1])
#       b = float(sys.argv[2])
    except IndexError:
#        print '\nusage: '+prog+' a b    (where a & b are numbers)\n'
        sys.exit(0)

    unoptFile = open('i.dat','r')
    gamma_data = []
    i_data = []
    for line in unoptFile:
        gamma_data.append(float(line.split()[0]))
        i_data.append(float(line.split()[1]))
    gamma_data, i_data = np.array(gamma_data), np.array(i_data)
    unopt_gamma =  gamma_data[np.argmin(i_data)]
    unoptFile.close()

    optFile = open('omega_opt_ideriv/i.dat','r')

    gamma_data = []
    i_data = []
    for line in optFile:
        gamma_data.append(float(line.split()[0]))
        i_data.append(float(line.split()[1]))
    gamma_data, i_data = np.array(gamma_data), np.array(i_data)
    opt_gamma = gamma_data[np.argmin(i_data)]
    optFile.close()

    print unopt_gamma - opt_gamma, ' '

    def fitChiSq(self, sampleSpec, returnChiSq=False):

        if np.all(sampleSpec.wave == self.wave):
            newSampleSpec = sampleSpec
            grid = self.values
        else:
            print "error in function"
            newSampleSpec = sampleSpec.interpolate(self.wave)            
            minIDx = self.wave.searchsorted(sampleSpec.wave[0])
            maxIDx = self.wave.searchsorted(sampleSpec.wave[-1])
            grid = self.values[:,minIDx:maxIDx]

        
        if newSampleSpec.var != None:
            var = newSampleSpec.var
        else:
            var = 1.

        if  newSampleSpec.dq != None:
            dqMask = newSampleSpec.dq
        else:
            dqMask = np.ones(grid.shape[1]).astype(bool)

        chiSq = ((grid[:,dqMask]-newSampleSpec.flux[dqMask])/var[dqMask])**2
        nu = (np.ones(grid.shape[0])*grid.shape[1]) - len(self.params) - 1
        redChiSq = np.sum(chiSq, axis=1) / nu

        if returnChiSq:
            return np.min(redChiSq), self.points[np.argmin(redChiSq)]
        else:
            return self.points[np.argmin(redChiSq)]

    def click(self,event):
        """
        What to do, if a click on the figure happens:
            1. Check which axis
            2. Get data coord's.
            3. Plot resulting data.
            4. Update Figure
        """
        if event.inaxes==self.overview:
            #Get nearest data
            xpos=np.argmin(np.abs(event.xdata-self.x))
            ypos=np.argmin(np.abs(event.ydata-self.y))
            
            #Check which mouse button:
            if event.button==1:
                #Plot it                
                c,=self.y_subplot.plot(self.y, self.z[:,xpos],label=str(self.x[xpos]))
                self.overview.axvline(self.x[xpos],color=c.get_color(),lw=2)

            elif event.button==3:
                #Plot it                
                c,=self.x_subplot.plot(self.x, self.z[ypos,:],label=str(self.y[ypos]))
                self.overview.axhline(self.y[ypos],color=c.get_color(),lw=2)

        if event.inaxes==self.y_subplot:
            ypos=np.argmin(np.abs(event.xdata-self.y))
            c=self.x_subplot.plot(self.x, self.z[ypos,:],label=str(self.y[ypos]))
            self.overview.axhline(self.y[ypos],color=c.get_color(),lw=2)

        if event.inaxes==self.x_subplot:
            xpos=np.argmin(np.abs(event.xdata-self.x))
            c,=self.y_subplot.plot(self.y, self.z[:,xpos],label=str(self.x[xpos]))
            self.overview.axvline(self.x[xpos],color=c.get_color(),lw=2)
        #Show it
        plt.draw()

def get_rectangular_subset(lon,lat,data,view,xmin,xmax,ymin,ymax,buffer=0):
    imin = np.argmin(np.abs(lon-xmin))
    imax = np.argmin(np.abs(lon-xmax))
    jmin = np.argmin(np.abs(lat-ymin))
    jmax = np.argmin(np.abs(lat-ymax))

    return lon[imin-buffer:imax+1+buffer], lat[jmin-buffer:jmax+1+buffer], grid_convert(grid_convert(data,view,'x+y+')[imin-buffer:imax+1+buffer,jmin-buffer:jmax+1+buffer],'x+y+',view)        

def kmeans(xx, centroids, maxIters = 20, minclust=30, maxDiff = 2):

  # Cluster Assignment step
  ca = np.array([np.argmin([np.dot(x_i-y_k, x_i-y_k) for y_k in centroids]) for x_i in xx])
  # all clusters have at least minclust?
  (unique, counts) = np.unique(ca, return_counts=True)
  for cc in counts:
    if cc < minclust:
      return("error: too few", np.array(centroids), ca)
  # Move centroids step
  centroids = np.array([xx[ca == k].mean(axis = 0) for k in range(centroids.shape[0])])

  iter=1
  while (iter<maxIters):
      # Cluster Assignment step
      canew = np.array([np.argmin([np.dot(x_i-y_k, x_i-y_k) for y_k in centroids]) for x_i in xx])
      # all clusters have at least minclust?
      (unique, counts) = np.unique(canew, return_counts=True)
      for cc in counts:
        if cc < minclust:
          return("error: too few", np.array(centroids), canew)
      numdiff = sum(ca != canew)
      if numdiff < maxDiff:
        return("converged", np.array(centroids), canew)
      ca = canew
      # Move centroids step
      centroids = np.array([xx[ca == k].mean(axis = 0) for k in range(centroids.shape[0])])
      iter += 1

  return("error: not converged", np.array(centroids), ca)

 def draw(self):
     self.spec_widg.on_draw(no_draw=True)
     # Add text?
     for kk,lls in enumerate(self.abssys_widg.all_abssys):
         # Label
         ipos = self.abssys_widg.all_items[kk].rfind('_')
         ilbl = self.abssys_widg.all_items[kk][ipos+1:]
         # Add text
         for wv,lbl in self.plt_wv:
             idx = np.argmin(np.abs(self.continuum.dispersion-wv*(1+lls.zabs)))
             self.spec_widg.ax.text(wv.value*(1+lls.zabs),
                 self.continuum.flux[idx],
                 '{:s}_{:s}'.format(ilbl,lbl), ha='center',
                 color='blue', size='small', rotation=90.)
     # Ticks for selected LLS
     idxl = self.get_sngl_sel_sys()
     if idxl is not None:
         lls = self.abssys_widg.all_abssys[idxl]
         # Label
         ipos = self.abssys_widg.all_items[idxl].rfind('_')
         ilbl = self.abssys_widg.all_items[idxl][ipos+1:]
         for line in lls.lls_lines:
             if line.wrest < 915.*u.AA:
                 continue
             idx = np.argmin(np.abs(self.continuum.dispersion-
                 line.wrest*(1+lls.zabs)))
             self.spec_widg.ax.text(line.wrest.value*(1+lls.zabs),
                 self.continuum.flux[idx],
                 '-{:s}'.format(ilbl), ha='center',
                 color='red', size='small', rotation=90.)
     # Draw
     self.spec_widg.canvas.draw()

def get_cmd_shape(color, mag):
    '''
    gets the outline of a cmd. Guesses at a large polygon, and then add points
    that are outside of the polygon, ignores points within.

    then polar sorts the result.
    returns: N,2 array.
    '''
    # make a guess at the polygon.
    left = (np.min(color), mag[np.argmin(color)])
    right = (np.max(color), mag[np.argmax(color)])
    up = (color[np.argmin(mag)], np.min(mag))
    down = (color[np.argmax(mag)], np.max(mag))
    verts = np.vstack((left, right, up, down))

    points = np.column_stack((color, mag))
    for point in points:
        if nxutils.pnpoly(point[0], point[1], verts) == 0.:
            # add point to verts
            col = verts[:, 0]
            m = verts[:, 1]
            col = np.append(col, point[0])
            m = np.append(m, point[1])
            # order the new points in a circle
            verts = polar_sort(zip(col, m))

    verts = np.append(verts, [verts[0]], axis=0)
    # plt.plot(verts[:, 0], verts[:, 1], lw = 2)
    return verts

    def _augmenting_row_reduction(self):
        """
        Augmenting row reduction step from LAPJV algorithm
        """
        unassigned = np.where(self._x == -1)[0]
        for i in unassigned:
            while True:
                #find smallest 2 values and indices
                temp = self.c[i] - self._v
                j1 = np.argmin(temp)
                u1 = temp[j1]
                temp[j1] = np.max(temp) + 1
                j2 = np.argmin(temp)
                u2 = temp[j2]

                if u1 < u2:
                    self._v[j1] -= u2 - u1
                elif self._y[j1] != -1:
                    j1 = j2
                k = self._y[j1]
                if k != -1:
                    self._x[k] = -1
                    self._x[i] = j1
                    self._y[j1] = i
                    i = k
                if np.allclose(u1, u2) or k == -1:
                    break

    def _make_tuples(self, key):
        print('Populating', key)
        spikes = np.vstack([s.squeeze() for s in (preprocess.Spikes.RateTrace() & key).fetch('rate_trace')])
        s = spikes.sum(axis=0)
        nans = np.isnan(s)

        key['leading_nans'] = int(nans[0])
        key['trailing_nans'] = int(nans[1])

        t = (preprocess.Sync() & key).fetch1('frame_times')  # does not need to be unique

        flip_first = (vis.Trial() * preprocess.Sync().proj('psy_id', trial_idx='first_trial') & key).fetch1('flip_times')
        flip_last = (vis.Trial() * preprocess.Sync().proj('psy_id', trial_idx='last_trial') & key).fetch1('flip_times')

        # (vis.Trial() * preprocess.Sync() & 'trial_idx between first_trial and last_trial')
        fro = np.atleast_1d(flip_first.squeeze())[0]
        to = np.atleast_1d(flip_last.squeeze())[
            -1]  # not necessarily where the stimulus stopped, just last presentation
        idx_fro = np.argmin(np.abs(t - fro))
        idx_to = np.argmin(np.abs(t - to)) + 1
        key['stimulus_nans'] = int(np.any(nans[idx_fro:idx_to]))
        if np.any(nans):
            key['nan_idx'] = nans
        key['stimulus_start'] = idx_fro + 1
        key['stimulus_end'] = idx_to

        self.insert1(key)

    def find_boundaries(self, all=False):
        """Find the local minima on either side of each peak

        Arguments:
            (none)
        """
        try:
            prevb = np.argmin(self.y[0:self._idx[0]])
        except IndexError:
            prevb = 0

        bounds = []

        if not all:
            pos = self._idx[self._keep]
        else:
            pos = self._idx

        npks = len(pos)

        for i in range(npks):
            thismax = pos[i]
            if i < npks-1:
                nextmax = pos[i + 1]
                relb = np.argmin(self.y[thismax:nextmax])
                nextb = relb + thismax
            else:
                nextmax = len(self.y)-1
                nextb = len(self.y)-1

            bounds.append([prevb, nextb])
            prevb = nextb

        self._bounds = np.array(bounds)

 def _transform_indices(self, key):
     """
     Transforms indices by snapping to the closest value if
     values are numeric, otherwise applies no transformation.
     """
     ndims = self.ndims
     if all(not isinstance(el, slice) for el in key):
         dim_inds = []
         for dim in self._cached_index_names:
             dim_type = self.get_dimension_type(dim)
             if isinstance(dim_type, type) and issubclass(dim_type, Number):
                 dim_inds.append(self.get_dimension_index(dim))
         str_keys = iter(key[i] for i in range(self.ndims)
                         if i not in dim_inds)
         num_keys = []
         if len(dim_inds):
             keys = list({tuple(k[i] if ndims > 1 else k for i in dim_inds)
                          for k in self.keys()})
             q = np.array([tuple(key[i] if ndims > 1 else key for i in dim_inds)])
             idx = np.argmin([np.inner(q - np.array(x), q - np.array(x))
                              if len(dim_inds) == 2 else np.abs(q-x)
                                  for x in keys])
             num_keys = iter(keys[idx])
         key = tuple(next(num_keys) if i in dim_inds else next(str_keys)
                     for i in range(self.ndims))
     elif any(not isinstance(el, slice) for el in key):
         index_ind = [idx for idx, el in enumerate(key)
                      if not isinstance(el, (slice, str))][0]
         dim_keys = np.array([k[index_ind] for k in self.keys()])
         snapped_val = dim_keys[np.argmin(np.abs(dim_keys-key[index_ind]))]
         key = list(key)
         key[index_ind] = snapped_val
         key = tuple(key)
     return key

def estimate_pk_parms_1d(x,f,pktype):
    """
    Gives initial guess of parameters for analytic fit of one dimensional peak
    data.

    Required Arguments:
    x -- (n) ndarray of coordinate positions
    f -- (n) ndarray of intensity measurements at coordinate positions x
    pktype -- string, type of analytic function that will be used to fit the data,
    current options are "gaussian","lorentzian","pvoigt" (psuedo voigt), and
    "split_pvoigt" (split psuedo voigt)
    

    Outputs:
    p -- (m) ndarray containing initial guesses for parameters for the input peaktype 
    (see peak function help for what each parameters corresponds to)
    """

    
    data_max=np.max(f)
#    lbg=np.mean(f[:2])
#    rbg=np.mean(f[:2]) 
    if((f[0]> (0.25*data_max)) and (f[-1]> (0.25*data_max))):#heuristic for wide peaks  
        bg0=0.
    elif (f[0]> (0.25*data_max)): #peak cut off on the left 
        bg0=f[-1]
    elif (f[-1]> (0.25*data_max)): #peak cut off on the right 
        bg0=f[0]
    else:
        bg0=(f[0]+f[-1])/2.    
    #bg1=(rbg-lbg)/(x[-1]-x[0])
    
    cen_index=np.argmax(f) 
    x0=x[cen_index]    
    A=data_max-bg0#-(bg0+bg1*x0)  
    
    num_pts=len(f)
   
    #checks for peaks that are cut off
    if cen_index == (num_pts-1):
        FWHM=x[cen_index]-x[np.argmin(np.abs(f[:cen_index]-A/2.))]#peak cut off on the left 
    elif cen_index == 0:
        FWHM=x[cen_index+np.argmin(np.abs(f[cen_index+1:]-A/2.))]-x[0] #peak cut off on the right 
    else:    
        FWHM=x[cen_index+np.argmin(np.abs(f[cen_index+1:]-A/2.))]-x[np.argmin(np.abs(f[:cen_index]-A/2.))]
        
    if FWHM <=0:##uh,oh something went bad
        FWHM=(x[-1]-x[0])/4. #completely arbitrary, set peak width to 1/4 window size

    
    
    if pktype=='gaussian' or pktype=='lorentzian':
        p=[A,x0,FWHM,bg0,0.]
    elif pktype=='pvoigt':
        p=[A,x0,FWHM,0.5,bg0,0.]
    elif pktype=='split_pvoigt':
        p=[A,x0,FWHM,FWHM,0.5,0.5,bg0,0.]
        
    p=np.array(p)
    return p

def gpu_nnc_predict(trX, trY, teX, metric='cosine', batch_size=4096):
    if metric == 'cosine':
        metric_fn = cosine_dist
    else:
        metric_fn = euclid_dist
    idxs = []
    for i in range(0, len(teX), batch_size):
        mb_dists = []
        mb_idxs = []
        for j in range(0, len(trX), batch_size):
            dist = metric_fn(floatX(teX[i:i+batch_size]), floatX(trX[j:j+batch_size]))
            if metric == 'cosine':
                mb_dists.append(np.max(dist, axis=1))
                mb_idxs.append(j+np.argmax(dist, axis=1))
            else:
                mb_dists.append(np.min(dist, axis=1))
                mb_idxs.append(j+np.argmin(dist, axis=1))                
        mb_idxs = np.asarray(mb_idxs)
        mb_dists = np.asarray(mb_dists)
        if metric == 'cosine':
            i = mb_idxs[np.argmax(mb_dists, axis=0), np.arange(mb_idxs.shape[1])]
        else:
            i = mb_idxs[np.argmin(mb_dists, axis=0), np.arange(mb_idxs.shape[1])]
        idxs.append(i)
    idxs = np.concatenate(idxs, axis=0)
    nearest = trY[idxs]
    return nearest

	def optimize_pa(self,fixed_pa=False,step=10,pa_init=0,pa_max=180):
		self.fixed_pa_best=False
		if fixed_pa:
			self.fixed_pa_best=True
			self.pa_best=fixed_pa
		else:
			self.pas=[]
			self.chi2=[]

			pa=pa_init
			while True:
				##
				## this could be made more efficient by only calling the sub-routine rotate_and_fft
				## but then need to watch out for differential vs. absolute rotations
				i=img2vis(self.f_model, self.pxscale, self.lam, oifits=self.oifits, pa=pa, phot=self.phot)

				self.chi2.append(i.vis_chi2())
				self.pas.append(pa)
				pa+=step
				if pa > pa_max:
					break
			
			## chose global chi2 minimum here for the moment, and plot PA vs. chi2 
			##    so that we see if global minimum is bad.
			self.pa_best = self.pas[np.argmin(self.chi2)]
			##
			## write out best chi**2 and name of model
			with open("chi2_min.txt","a") as f:
				txt="{0:06.0f} -- {1:5.3f} -- {2:5.2f} -- {3}\n".format(self.chi2[np.argmin(self.chi2)], self.pxscale, self.f_p, self.f_model)
				f.write(txt)
		self.rotate_and_fft(self.pa_best)

 def correct_shift_at_349nm(self):
     shift = self.values[closest_wavelength_index(self.nm,349)] - self.values[closest_wavelength_index(self.nm,348)]
     min_ = min(np.argmin(self.nm), closest_wavelength_index(self.nm,348))
     max_ = min(np.argmin(self.nm), closest_wavelength_index(self.nm,348))
     for i in range(min_,max_):
         self.values[i] -= shift            
     self.state += ' 349nm_shift_corrected'     

def online_k_means(k,b,t,X_in):
    random_number = 11232015
    random_num = np.random.randint(X_in.shape[0], size =300 )
    rng = np.random.RandomState(random_number)
    permutation1 = rng.permutation(len(random_num))
    random_num = random_num[permutation1]
    x_input = X_in[random_num]
    c,l = mykmeansplusplus(x_input,k,t)
    v = np.zeros((k))
    for i in range(t):
        random_num = np.random.randint(X_in.shape[0], size = b)
        rng = np.random.RandomState(random_number)
        permutation1 = rng.permutation(len(random_num))
        random_num = random_num[permutation1]
        M = X_in[random_num]
        Y = cdist(M,c,metric='euclidean', p=2, V=None, VI=None, w=None)
        clust_index = np.argmin(Y,axis = 1)
        for i in range(M.shape[0]):
            c_in = clust_index[i]
            v[c_in] += 1
            ita = 1 / v[c_in]
            c[c_in] = np.add(np.multiply((1 - ita),c[c_in]),np.multiply(ita,M[i]))
    Y_l = cdist(X_in,c,metric='euclidean', p=2, V=None, VI=None, w=None)
    l = np.argmin(Y_l,axis = 1)        
    return c,l

def rescale_frontoparallel(p_fp,box_fp,p_im):
    """
    The fronto-parallel image region is rescaled to bring it in 
    the same approx. size as the target region size.

    p_fp : nx2 coordinates of countour points in the fronto-parallel plane
    box  : 4x2 coordinates of bounding box of p_fp
    p_im : nx2 coordinates of countour in the image

    NOTE : p_fp and p are corresponding, i.e. : p_fp[i] ~ p[i]

    Returns the scale 's' to scale the fronto-parallel points by.
    """
    l1 = np.linalg.norm(box_fp[1,:]-box_fp[0,:])
    l2 = np.linalg.norm(box_fp[1,:]-box_fp[2,:])

    n0 = np.argmin(np.linalg.norm(p_fp-box_fp[0,:][None,:],axis=1))
    n1 = np.argmin(np.linalg.norm(p_fp-box_fp[1,:][None,:],axis=1))
    n2 = np.argmin(np.linalg.norm(p_fp-box_fp[2,:][None,:],axis=1))

    lt1 = np.linalg.norm(p_im[n1,:]-p_im[n0,:])
    lt2 = np.linalg.norm(p_im[n1,:]-p_im[n2,:])

    s =  max(lt1/l1,lt2/l2)
    if not np.isfinite(s):
        s = 1.0
    return s

def lattice(unit_cell_size, shape):
    """
    We should just have delta functions at 1/d
    however, if the unit cell does not divide the 
    detector shape evenly then these fall between 
    pixels. 
    """
    # generate the q-space coordinates
    qi = np.fft.fftfreq(shape[0])
    qj = np.fft.fftfreq(shape[1])
    qk = np.fft.fftfreq(shape[2])

    # generate the recirocal lattice points from the unit cell size
    qs_unit = np.meshgrid(
        np.fft.fftfreq(unit_cell_size[0]),
        np.fft.fftfreq(unit_cell_size[1]),
        np.fft.fftfreq(unit_cell_size[2]),
        indexing="ij",
    )

    # now we want qs[qs_unit] = 1.
    lattice = np.zeros(shape, dtype=np.float)
    for ii, jj, kk in zip(qs_unit[0].ravel(), qs_unit[1].ravel(), qs_unit[2].ravel()):
        i = np.argmin(np.abs(ii - qi))
        j = np.argmin(np.abs(jj - qj))
        k = np.argmin(np.abs(kk - qk))
        lattice[i, j, k] = 1.0

    return lattice

  def _readDem(self):
    """
    Read coordinates defining DEM and create vectors of x, y, and z values.
    """

    # Load each coordinate as a numpy array.
    x, y, z = numpy.loadtxt(self.inputDem, dtype=numpy.float64, unpack=True)

    self.numZIn = len(z)
    if (y[0] == y[1]):
      # Ordered by rows.
      self.numXIn = max(numpy.argmax(x) + 1, numpy.argmin(x) + 1)
      self.xIn = x[0:self.numXIn]
      self.numYIn = self.numZIn/self.numXIn
      self.yIn = y[0:self.numZIn:self.numXIn]
      self.zIn = numpy.reshape(z, (self.numYIn, self.numXIn))
    else:
      # Ordered by columns.
      self.numYIn = max(numpy.argmax(y) + 1, numpy.argmin(y) + 1)
      self.yIn = y[0:self.numYIn]
      self.numXIn = self.numZIn/self.numYIn
      self.xIn = x[0:self.numZIn:self.numYIn]
      self.ZIn = numpy.transpose(numpy.reshape(z, (self.numXIn, self.numYIn)))

    if (self.xIn[0] > self.xIn[1]):
      self.xIn = numpy.flipud(self.xIn)
      self.zIn = numpy.fliplr(self.zIn)
    if (self.yIn[0] > self.yIn[1]):
      self.yIn = numpy.flipud(self.yIn)
      self.zIn = numpy.flipud(self.zIn)

    return