def agg_lut(lathi,lonhi,latlo,lonlo):
    # outputlut=np.zeros((lathi.shape[0],lathi.shape[1],2))-9999
    print("Computing Geographic Look Up Table")
    outputlut=np.empty((latlo.shape[0],latlo.shape[1],2),dtype=list)
    nx=lathi.shape[1]
    ny=lathi.shape[0]

    maxdist=((latlo[0,0]-latlo[1,1])**2 + (lonlo[0,0]-lonlo[1,1])**2)/1.5
    
    for i in range(ny):
        dists=(latlo-lathi[i,0])**2+(lonlo-lonhi[i,0])**2
        y,x=np.unravel_index(dists.argmin(),dists.shape)
        # print(i,ny,y,x,lathi[i,0],lonhi[i,0],latlo[y,x],lonlo[y,x])
        for j in range(nx):
            xmax=min(nx,x+3)
            xmin=max(0,x-3)
            ymax=min(ny,y+3)
            ymin=max(0,y-3)
            
            windists=((latlo[ymin:ymax,xmin:xmax]-lathi[i,j])**2+
                      (lonlo[ymin:ymax,xmin:xmax]-lonhi[i,j])**2)
            yoff,xoff=np.unravel_index(windists.argmin(),windists.shape)
            x=xoff+xmin
            y=yoff+ymin
            if not outputlut[y,x,0]:
                outputlut[y,x,0]=list()
                outputlut[y,x,1]=list()
            
            if windists[yoff,xoff]<maxdist:
                outputlut[y,x,0].append(i)
                outputlut[y,x,1].append(j)
                
    return outputlut

    def test_grad(self):
        eps = 1e-7
        f, args, vals = self.get_args()
        output0 = f(*vals)

        # Go through and backpropagate all of the gradients from the outputs
        grad0 = []
        for i in range(len(output0) - 2):
            grad0.append([])
            for j in range(output0[i].size):
                ind = np.unravel_index(j, output0[i].shape)

                g = theano.function(
                    args, theano.grad(self.op(*args)[i][ind], args))
                grad0[-1].append(g(*vals))

        # Loop over each input and numerically compute the gradient
        for k in range(len(vals)):
            for l in range(vals[k].size):
                inner = np.unravel_index(l, vals[k].shape)
                vals[k][inner] += eps
                plus = f(*vals)
                vals[k][inner] -= 2*eps
                minus = f(*vals)
                vals[k][inner] += eps

                # Compare to the backpropagated gradients
                for i in range(len(output0) - 2):
                    for j in range(output0[i].size):
                        ind = np.unravel_index(j, output0[i].shape)
                        delta = 0.5 * (plus[i][ind] - minus[i][ind]) / eps
                        ref = grad0[i][j][k][inner]
                        assert np.abs(delta - ref) < 2*eps, \
                            "{0}".format((k, l, i, j, delta, ref, delta-ref))

def test_normalization():
    """Test that `match_template` gives the correct normalization.

    Normalization gives 1 for a perfect match and -1 for an inverted-match.
    This test adds positive and negative squares to a zero-array and matches
    the array with a positive template.
    """
    n = 5
    N = 20
    ipos, jpos = (2, 3)
    ineg, jneg = (12, 11)
    image = np.full((N, N), 0.5)
    image[ipos:ipos + n, jpos:jpos + n] = 1
    image[ineg:ineg + n, jneg:jneg + n] = 0

    # white square with a black border
    template = np.zeros((n + 2, n + 2))
    template[1:1 + n, 1:1 + n] = 1

    result = match_template(image, template)

    # get the max and min results.
    sorted_result = np.argsort(result.flat)
    iflat_min = sorted_result[0]
    iflat_max = sorted_result[-1]
    min_result = np.unravel_index(iflat_min, result.shape)
    max_result = np.unravel_index(iflat_max, result.shape)

    # shift result by 1 because of template border
    assert np.all((np.array(min_result) + 1) == (ineg, jneg))
    assert np.all((np.array(max_result) + 1) == (ipos, jpos))

    assert np.allclose(result.flat[iflat_min], -1)
    assert np.allclose(result.flat[iflat_max], 1)

 def classify_obj(self, zchi2arr1, zfit1, flags1, zchi2arr2, zfit2, flags2):
     flag_val = int('0b100',2) # From BOSS zwarning flag definitions
     for ifiber in xrange(zchi2arr1.shape[0]):
         self.minrchi2[ifiber,0] = n.min(zchi2arr1[ifiber]) / (self.npixflux) # Calculate reduced chi**2 values to compare templates of differing lengths
         if zchi2arr2 != None: self.minrchi2[ifiber,1] = n.min(zchi2arr2[ifiber]) / (self.npixflux)
         minpos = self.minrchi2[ifiber].argmin() # Location of best chi2 of array of best (individual template) chi2s
         
         if minpos == 0: # Means overall chi2 minimum came from template 1
             self.type.append('GALAXY')
             self.minvector.append(zfit1.minvector[ifiber])
             self.z[ifiber] = zfit1.z[ifiber]
             self.z_err[ifiber] = zfit1.z_err[ifiber]
             minloc = n.unravel_index(zchi2arr1[ifiber].argmin(), zchi2arr1[ifiber].shape)[:-1]
             self.zwarning = n.append(self.zwarning, flags1[ifiber])
             argsort = self.minrchi2[ifiber].argsort()
             if len(argsort) > 1:
                 if argsort[1] == 1:
                     if ( n.min(zchi2arr2[ifiber]) - n.min(zchi2arr1[ifiber]) ) < zfit1.threshold: self.zwarning[ifiber] = int(self.zwarning[ifiber]) | flag_val #THIS THRESHOLD PROBABLY ISN'T RIGHT AND NEEDS TO BE CHANGED
         
         elif minpos == 1: # Means overall chi2 minimum came from template 2
             self.type.append('STAR')
             self.minvector.append(zfit2.minvector[ifiber])
             self.z[ifiber] = zfit2.z[ifiber]
             self.z_err[ifiber] = zfit2.z_err[ifiber]
             minloc = n.unravel_index(zchi2arr2[ifiber].argmin(), zchi2arr2[ifiber].shape)[:-1]
             self.zwarning = n.append(self.zwarning, flags2[ifiber])
             argsort = self.minrchi2[ifiber].argsort()
             if argsort[1] == 0:
                 if ( n.min(zchi2arr1[ifiber]) - n.min(zchi2arr2[ifiber]) ) < zfit2.threshold: self.zwarning[ifiber] = int(self.zwarning[ifiber]) | flag_val

def get_x_y_wind(lat, lon, dataset, model):
    indices = []
    # m = Basemap(rsphere=(6371229.00, 6356752.3142), projection='merc',
    #             llcrnrlat=40.5833284543, urcrnrlat=47.4999927992,
    #                    llcrnrlon=-129, urcrnrlon=-123.7265625)
    # xpt, ypt = m(lon,lat)
    # print "x ", xpt/100
    # print "y ", ypt/100
    print "lat ", lat
    print "lon ", lon
    lat_name = 'lat'
    lon_name = 'lon'
    file_lats = dataset.variables[lat_name][:]
    file_lons = dataset.variables[lon_name][:]
    file_lat = np.abs(file_lats - lat).argmin()
    file_lon = np.abs(file_lons - lon).argmin()
    print "Argmin lat", file_lat
    print "Argmin lon", file_lon
    file_lat = np.unravel_index(file_lat, file_lats.shape)
    file_lon = np.unravel_index(file_lon, file_lons.shape)
    print "Unravel lat", file_lat
    print "Unravel lon", file_lon
    file_lat_y = file_lat[0]
    file_lat_x = file_lat[1]
    print"lat y x", file_lat_y, file_lat_x
    file_lon_y = file_lon[0]
    file_lon_x = file_lon[1]
    print"lon y x", file_lon_y, file_lon_x
    print "lat ", file_lats[file_lat_y][file_lat_x]
    print "lon ", file_lons[file_lon_y][file_lon_x]
    indices.append(file_lat_y)
    indices.append(file_lat_x)
    return indices

    def _fixEvenKernel(kernel):
        """Take a kernel with even dimensions and make them odd, centered correctly.

        Parameters
        ----------
        kernel : `numpy.array`
            a numpy.array

        Returns
        -------
        out : `numpy.array`
            a fixed kernel numpy.array
        """
        # Make sure the peak (close to a delta-function) is in the center!
        maxloc = np.unravel_index(np.argmax(kernel), kernel.shape)
        out = np.roll(kernel, kernel.shape[0]//2 - maxloc[0], axis=0)
        out = np.roll(out, out.shape[1]//2 - maxloc[1], axis=1)
        # Make sure it is odd-dimensioned by trimming it.
        if (out.shape[0] % 2) == 0:
            maxloc = np.unravel_index(np.argmax(out), out.shape)
            if out.shape[0] - maxloc[0] > maxloc[0]:
                out = out[:-1, :]
            else:
                out = out[1:, :]
            if out.shape[1] - maxloc[1] > maxloc[1]:
                out = out[:, :-1]
            else:
                out = out[:, 1:]
        return out

def nearest(array, value):
    """
    Find nearest position in array to value.
    Parameters
    ----------
    array : 'array_like'
    value: 'float', 'list'

    Returns
    ----------
    Searches n x 1 and n x m x 1  arrays for floats.
    Searches n x m and n x m x p arrays for lists of size m or p.
    """

    if isinstance(array, (list, tuple)):
        array = np.asarray(array)

    if isinstance(value, (float, int)):
        pos = (np.abs(array - value)).argmin()
        if len(array.shape) == 2:
            return np.unravel_index(pos, array.shape)
        else:
            return pos
    else:
        pos = (np.sum((array - value)**2, axis=1)**(1 / 2)).argmin()
        if len(array.shape) == 3:
            return np.unravel_index(pos, array.shape)
        else:
            return pos

    def write_cv_results_toFile(self, scores, precision_scores, recall_scores, param_grid, result_path):
        
        final_results_dict = {}
        all_scores_dict = {'error' : scores, 'precision' : precision_scores, 'recall' : recall_scores}
        for score in all_scores_dict:
            avg_score = np.mean(all_scores_dict[score], axis=1)
            std_score = np.std(all_scores_dict[score], axis=1)

            if score == 'error':
                opt_ind = np.unravel_index(np.argmin(avg_score), avg_score.shape)
            else:
                opt_ind = np.unravel_index(np.argmax(avg_score), avg_score.shape)
            
            final_results_dict[score] = avg_score[opt_ind]
            final_results_dict[score + '_std'] = std_score[opt_ind]
            params_dict = {}
            for element in self.grid_dictionary:
                params_dict[element] =  param_grid[opt_ind[0]][self.grid_dictionary[element]]
            final_results_dict[score + '_params'] = params_dict
            n_estimators = opt_ind[1] + 1
            final_results_dict[score + '_params'].update(dict(n_trees = str(n_estimators)))
            if not 'fe_params' in params_dict:
                final_results_dict[score + '_params'].update(dict(fe_params = None))
        
        final_results_dict['channel_type'] = self.config.channel_type

        json.dump(final_results_dict, open(result_path,'w'))  

def nonMaximumSuppression(I, dims, pos, c, maxima):
    # if pos== true, find local maximum, else - find local minimum
    #copy image
    I_temp = I
    w = dims[0]
    h = dims[1]
    n = 5
    margin = 5
    tau = 50
#try np here too
    for i in range(n + margin, w - n - margin):
        #    print "1 entered 1st cycle i=", i
        for j in range(n + margin, h - n - margin):
            window = I_temp[i:i+n, j:j+n]

            if np.sum(window) < 0:
                mval = 0
            elif pos:
                mval = np.amax(window)
                mcoords = np.unravel_index(np.argmax(window), window.shape) + np.array((i,j))
            else:
                mval = np.amin(window)
                mcoords = np.unravel_index(np.argmax(window), window.shape) + np.array((i,j))
            I_temp [i:i+n, j:j+n] = -1
            print mcoords, "mval= ", mval
            if pos:
                if mval >= tau:
                    maxima.append(maximum(mcoords[0],mcoords[1], mval, c))

            else:
                if mval <= tau and mval > 0:
                    maxima.append(maximum(mcoords[0],mcoords[1], mval, c))

  def reorder_in_x(self,xmat):

    xfactors=self.xfactors
    tdim=self.xmatrix.shape
    ndim=len(tdim)
    mf=xfactors.shape[0]
    xmatredim=np.zeros(mf,dtype=int)
    for i in range(mf):
      xmatredim[i]=np.prod(xfactors[i,:])
    xmatre=np.zeros(xmatredim)

    pdim=np.prod(tdim)
    i_in=np.arange(pdim)
    i_inx=np.array(np.unravel_index(i_in,tdim))
    i_outx=[None]*ndim
    for i in range(ndim):
      i_outx[i]=np.array(np.unravel_index(i_inx[i],xfactors[:,i]))
    i_outx=np.array(i_outx)
    
    i_out=[None]*mf
    for i in range(mf):
      i_out[i]=np.array(np.ravel_multi_index(i_outx[:,i],xfactors[i]))
    i_out=np.array(i_out)

    xmatre[tuple(i_out)]=self.xmatrix[tuple(i_inx)]

    return(xmatre)

  def invert_reorder_in_x(self,xmatre,xfactors):

    (mf,nf)=xfactors.shape
    ## tdim2=xmatre.shape
    ## ndim2=len(tdim2)

    ndim=nf
    tdim=np.zeros(nf,dtype=int)
    for i in range(nf):
      tdim[i]=np.prod(xfactors[:,i])
    xmatrix=np.zeros(tdim)

    pdim=np.prod(tdim)
    i_in=np.arange(pdim)
    i_inx=np.array(np.unravel_index(i_in,tdim))

    i_outx=[None]*ndim
    for i in range(ndim):
      i_outx[i]=np.array(np.unravel_index(i_inx[i],xfactors[:,i]))
    i_outx=np.array(i_outx)
    
    i_out=[None]*mf
    for i in range(mf):
      i_out[i]=np.array(np.ravel_multi_index(i_outx[:,i],xfactors[i]))
    i_out=np.array(i_out)

    xmatrix[tuple(i_inx)]=xmatre[tuple(i_out)]
    xmatrix=xmatrix.astype(np.uint8)

    return(xmatrix)

def seuillage_80(masque,image,handle=False):
	"""
		Permet de retirer 20 pourcent des pixels les plus clairs.
	"""
	nb_pixel=sum(sum(masque))
	seuil=nb_pixel*20/100
	image=humuscle(image)
	if handle:
		print('nombre de pixel avant seuillage')
		print(nb_pixel)
		print('nombre de pixel a supprimer')
		print(seuil)

	if handle:
		imshow('masque initial',masque)
	k=0
	while(k<seuil):
		masque[np.unravel_index(np.argmax(image),image.shape)]=0
		image[np.unravel_index(np.argmax(image),image.shape)]=image.min()
		k+=1

	if handle:
		imshow('masque apres seuillage a 80%',masque)
		print('nombre de pixel apres seuillage')
		print(sum(sum(masque)))
	return masque

def klst(X):
    history=[]
    add_vec =[]    
    print("データ行列の大きさ ", X.shape)
    count = X.shape[0]
    dist = calc(X)
    print(dist.shape)
    exit()
    i,j = np.unravel_index(dist.argmax(), dist.shape)
    max = dist[i,j]
    ind=[x for x in range(0,count)]
    
    while dist.shape[0] > 1:
        dist = pd.DataFrame(dist, index = ind ,columns = ind)

        p =slm(dist.values)
        print("最短距離ベクトル", p)
    #     print(dist)
        dist, add_vec = cut_n_comb(dist, p, add_vec)
        ind= dist.index
        dist = dist.values
        i,j = np.unravel_index(dist.argmax(), dist.shape)
        max = dist[i,j]
        if max == 0:#の例外処理が必要
            break
        print("合成結果：", add_vec)
        count -= 1
        print(count-1, "残り処理回数")
        print("-----------------------------")

def plot_drain_pit(self, pit, drain, prop, s, elev, area):
    from matplotlib import pyplot
    
    cmap = 'Greens'

    ipit, jpit = np.unravel_index(pit, elev.shape)
    Idrain, Jdrain = np.unravel_index(drain, elev.shape)
    Iarea, Jarea = np.unravel_index(area, elev.shape)

    imin = max(0, min(ipit, Idrain.min(), Iarea.min())-1)
    imax = min(elev.shape[0], max(ipit, Idrain.max(), Iarea.max()) + 2)
    jmin = max(0, min(jpit, Jdrain.min(), Jarea.min())-1)
    jmax = min(elev.shape[1], max(jpit, Jdrain.max(), Jarea.max()) + 2)
    roi = (slice(imin, imax), slice(jmin, jmax))

    pyplot.figure()
    
    ax = pyplot.subplot(221);
    pyplot.axis('off')
    im = pyplot.imshow(elev[roi], interpolation='none'); im.set_cmap(cmap)
    pyplot.plot(jpit-jmin, ipit-imin, lw=0, color='k', marker='$P$', ms=10)
    pyplot.plot(Jarea-jmin, Iarea-imin, 'k.')
    pyplot.plot(Jdrain-jmin, Idrain-imin, lw=0, color='k', marker='$D$', ms=10)

    pyplot.subplot(223, sharex=ax, sharey=ax); pyplot.axis('off')
    im = pyplot.imshow(elev[roi], interpolation='none'); im.set_cmap(cmap)
    for i, j, val in zip(Idrain, Jdrain, prop):
        pyplot.plot(j-jmin, i-imin, lw=0, color='k', marker='$%.2f$' % val, ms=16)

    pyplot.subplot(224, sharex=ax, sharey=ax); pyplot.axis('off')
    im = pyplot.imshow(elev[roi], interpolation='none'); im.set_cmap(cmap)
    for i, j, val in zip(Idrain, Jdrain, s):
        pyplot.plot(j-jmin, i-imin, lw=0, color='k', marker='$%.2f$' % val, ms=16)

    pyplot.tight_layout()

def calculateArea(flowDirectionGrid,flowAcc,areas,noDataValue,mask = None):

    # Get the sorted indices of the array in reverse order (e.g. largest first)
    idcs = flowAcc.argsort(axis= None)

    #Mask out the indexes outside of the mask
    if not mask is None:
        allRows,allCols = np.unravel_index(idcs,flowAcc.shape)
        gdIdcs = mask[allRows,allCols]
        idcs = idcs[gdIdcs]

    #Loop through all the indices in sorted order
    for idx in idcs:
        #Get the currents row column index
        [i, j] = np.unravel_index(idx, flowAcc.shape)  # Get the row/column indices

        #Get the index of the point downstream of this and the distance to that point
        [downI,downJ,inBounds] = getFlowToCell(i,j,flowDirectionGrid[i,j],flowAcc.shape[0],flowAcc.shape[1])

        #So long as we are not draining to the outskirts, proceed
        if inBounds and not flowAcc[downI,downJ] == noDataValue:

            #Accumulate area
            areas[downI,downJ] += areas[i,j]

    return areas

def grid_maximum(matrix):
    """ Find where in grid highest probability lies """
    
    arr = np.copy(matrix)
    # 1st maximum
    i, j = np.unravel_index(arr.argmax(), arr.shape)
        
    lon_bins = np.linspace(llcrnrlon, urcrnrlon+2, xBins)
    lat_bins = np.linspace(llcrnrlat, urcrnrlat+1, xBins*xyRatio)

    lon_corr = lon_bins[1]-lon_bins[0]
    lat_corr = lat_bins[1]-lat_bins[0]
    
    #for lon in lon_bins:
    #    for lat in lat_bins:
    #        coord = lat+lat_corr*0.5, lon+lon_corr*0.5
    #        print rg.get(coord)['admin1']

    maximum = lat_bins[i]+lat_corr*0.5, lon_bins[j]+lon_corr*0.5

    # 2nd maximum
    arr[i, j] = 0
    i, j = np.unravel_index(arr.argmax(), arr.shape)
    second_maximum = lat_bins[i]+lat_corr*0.5, lon_bins[j]+lon_corr*0.5

    # 3rd maximum
    arr[i, j] = 0
    i, j = np.unravel_index(arr.argmax(), arr.shape)
    third_maximum = lat_bins[i]+lat_corr*0.5, lon_bins[j]+lon_corr*0.5

    return maximum, second_maximum, third_maximum

def split_quater(img1,img2,n):
    iq1 = img1
    iq2 = img2
    if n == 1:
        iq1 = np.copy(iq1[0:239,0:239])
        iq2 = np.copy(iq2[0:239,0:239])
    elif n == 2:
        iq1 = np.copy(iq1[0:239,240:479])
        iq2 = np.copy(iq2[0:239,240:479])
    elif n == 3:
        iq1 = np.copy(iq1[240:479,0:239])
        iq2 = np.copy(iq2[240:479,0:239])
    elif n == 4:
        iq1 = np.copy(iq1[240:479,240:479])
        iq2 = np.copy(iq2[240:479,240:479])
    elif n == 5:
        iq1 = img1
        iq2 = img2
    iq1 = cv2.GaussianBlur(iq1,(5,5),0)
    iq2 = cv2.GaussianBlur(iq2,(5,5),0)

    corr_img11 = cross_image(iq1,iq1)
    corr_img12 = cross_image(iq1,iq2)

    #cv2.imwrite('corr_img11.jpg',corr_img11)
    #cv2.imwrite('corr_img12.jpg',corr_img22)


    y_self,x_self = np.unravel_index(np.argmax(corr_img11), corr_img11.shape)
    y,x = np.unravel_index(np.argmax(corr_img12), corr_img12.shape)

    y_diff = y - y_self
    x_diff = x - x_self
    print(y_diff,x_diff)
    return (y_diff,x_diff)

    def from_hdu(cls, hdu, hdu_bands=None):
        """Make a WcsNDMap object from a FITS HDU.

        Parameters
        ----------
        hdu : `~astropy.io.fits.BinTableHDU` or `~astropy.io.fits.ImageHDU`
            The map FITS HDU.
        hdu_bands : `~astropy.io.fits.BinTableHDU`
            The BANDS table HDU.
        """
        geom = WcsGeom.from_header(hdu.header, hdu_bands)
        shape = tuple([ax.nbin for ax in geom.axes])
        shape_wcs = tuple([np.max(geom.npix[0]),
                           np.max(geom.npix[1])])
        meta = cls._get_meta_from_header(hdu.header)
        map_out = cls(geom, meta=meta)

        # TODO: Should we support extracting slices?
        if isinstance(hdu, fits.BinTableHDU):
            pix = hdu.data.field('PIX')
            pix = np.unravel_index(pix, shape_wcs[::-1])
            vals = hdu.data.field('VALUE')
            if 'CHANNEL' in hdu.data.columns.names and shape:
                chan = hdu.data.field('CHANNEL')
                chan = np.unravel_index(chan, shape[::-1])
                idx = chan + pix
            else:
                idx = pix

            map_out.set_by_idx(idx[::-1], vals)
        else:
            map_out.data = hdu.data

        return map_out

    def maxmindam(self,probsize,meanvasdam):

        m = probsize[0]
        n = probsize[1]

        damage = vasdam.damcalc(self,probsize,meanvasdam)

        maxdam = np.argmax(damage)
        mindam = np.argmin(damage)

        maxdamval = np.max(damage)
        mindamval = np.min(damage)
        
        dimsdam = damage.shape
        
        maxidx = np.unravel_index(maxdam,dimsdam)
        minidx = np.unravel_index(mindam,dimsdam)

        maxvasdam = np.array([maxidx[0],m+maxidx[0],maxidx[1],n+maxidx[1]]).reshape([2,2])
        minvasdam = np.array([minidx[0],m+minidx[0],minidx[1],n+minidx[1]]).reshape([2,2])

        location = np.array([maxvasdam,int(maxdamval),minvasdam,int(mindamval)])
        

        return  location

 def iterkeys(self,index):
     #import pdb; pdb.set_trace()
     if not isinstance(index,tuple) and self.shape[0] == 1:
         index = (1,index)
     if isinstance(index, int):
         key = np.unravel_index(index-1, self.shape, order='F')
         yield tuple(k+1 for k in key)
     elif isinstance(index,slice):
         index = range((index.start or 1)-1,
                       index.stop or np.prod(self.shape),
                       index.step or 1)
         for key in np.transpose(np.unravel_index(index, self.shape, order='F')): # 0-based
             yield tuple(k+1 for k in key)
     elif isinstance(index,(list,np.ndarray)):
         index = np.asarray(index)-1
         for key in np.transpose(np.unravel_index(index, self.shape, order='F')):
             yield tuple(k+1 for k in key)
     else:
         assert isinstance(index,tuple),index.__class__
         indices = []  # 1-based
         for i,ix in enumerate(index):
             if isinstance(ix,slice):
                 indices.append(np.arange((ix.start or 1),
                                          (ix.stop  or self.shape[i]) + 1,
                                          ix.step   or 1,
                                          dtype=int))
             else:
                 indices.append(np.asarray(ix))
         assert len(index) == 2
         indices[0].shape = (-1,1)
         for key in np.broadcast(*indices):
             yield tuple(map(int,key))