def callback2(data):
    global face_det
    global filename
    global datei
    #rospy.loginfo(rospy.get_caller_id()+"CDIA data:" + str(len(data.head_detections)))
    if len(data.head_detections)==1 and len(data.head_detections[0].face_detections) == 1:
        face_det+=1
        print face_det
        bridge = CvBridge()
        image = bridge.imgmsg_to_cv2(data.head_detections[0].color_image,"rgb8")
        depth = bridge.imgmsg_to_cv2(data.head_detections[0].depth_image,"32FC3")
        #cv2.imshow("image",image)
        cv2.imshow("depth",depth)
        cv2.waitKey()
        print depth.shape
        set_name = datei[1]
        depth_path="/home/stefan/rgbd_db_heads/"+set_name+"/"+filename+"_d.xml"
        img_path="/home/stefan/rgbd_db_heads/"+set_name+"/"+filename+"_c.bmp"
        #path+"/"+str(dir)+"/"+os.path.splitext(file)[0]+".xml"
        depth_slice1 = depth[:,:,0].astype(numpy.float32)
        depth_slice2 = depth[:,:,1].astype(numpy.float32)
        depth_slice3 = depth[:,:,2].astype(numpy.float32)

        for(r,c),value in numpy.ndenumerate(depth_slice1):
            if depth_slice1[r,c]==0:
                depth_slice1[r,c]=numpy.nan
        for(r,c),value in numpy.ndenumerate(depth_slice2):
            if depth_slice2[r,c]==0:
                depth_slice2[r,c]=numpy.nan
        print filename

def detectCircles(img, r, useGradient):
    grayimg = rgb2gray(img)
    edges = cv2.Canny(img,100,200)
    ax[0].imshow(edges, cmap=plt.cm.gray)
    ax[0].set_title('after canny image operation')
    if useGradient == 0:
        accumulator1 = np.zeros(edges.shape)
        for (i,j),value in np.ndenumerate(edges):
            if value:
                for t_idx in np.arange(0,2*math.pi,math.pi/100):
                    a = int(i - (r * math.cos(t_idx)));
                    b = int(j + (r * math.sin(t_idx)));
                    if a>0 and b>0 and a < accumulator1.shape[0] and b < accumulator1.shape[1]:
                        accumulator1[a, b] += 1
        
        print accumulator1
        ax[1].imshow(accumulator1, cmap=plt.cm.gray)
        ax[1].set_title('Accumulator array without using gradient')
    else:
        dx = ndimage.sobel(grayimg, axis=0, mode='constant')
        dy = ndimage.sobel(grayimg, axis=1, mode='constant')
        accumulator = np.zeros(edges.shape)
        for (i,j),value in np.ndenumerate(edges):
            if value:
                gradient = math.atan(-dx[i,j]/(dy[i,j]+0.00001))
                for theta in np.arange(gradient-math.pi/4,gradient+math.pi/4,math.pi/100):
                    a = int(i - (r * math.cos(theta)));
                    b = int(j + (r * math.sin(theta)));
                    if a < accumulator.shape[0] and b < accumulator.shape[1]:
                        accumulator[a, b] += 1
        ax[1].imshow(accumulator, cmap=plt.cm.gray)
        ax[1].set_title('Accumulator array with gradient')
        print accumulator
    return 

    def traverse_cells(self, visitor):
	""" Call a visitor function on each cell in the Partition. The visitor
	should look like this::

	  def visitor(path):
	      pass

        The path passed in is a list of PathNodes describing the nesting of
        the cell within partitions.  From the path, you can get all the
        containing partitions of the element it points to, the element,
        and both n-dimensional and flat indices of the element within each
        partition.  See PathNode for more details.
        """
        if self.children.size:
            for index, child in np.ndenumerate(self.children):
                child.traverse_cells(visitor)
        else:
            for index, elt in np.ndenumerate(self.box):
                # Build a list of PathNodes containing ancestor partitions
                path = [Partition.PathNode(p) for p in self.ancestors]
                path[-1].index = index
                # assign index of elt within each partition to each PathNode
                i = -2
                while i >= -len(path):
                    child = path[i+1]
                    path[i].index = child.partition.self_to_parent(child.index)
                    i -= 1
                # Now visit the element with its path.
                visitor(path)

def getPathData(data, param):
    path_data = []
    p = param

    tree = Grid_adaptiveM(data, p.Eps, p)
    tree.buildIndex()
    tree.adjustConsistency()

    leaf_boxes = []
    for (_, _), l1_child in np.ndenumerate(tree.root.children):
        if not l1_child.n_isLeaf and l1_child.children is not None:
            for (_, _), l2_child in np.ndenumerate(l1_child.children):  # child1 is a first-level cell
                leaf_boxes.append((l2_child.n_box, l2_child.a_count))
        leaf_boxes.append((l1_child.n_box, l1_child.a_count))

    for data in leaf_boxes:
        # [[x_min,y_min],[x_max,y_max]]
        path = []
        box = data[0]
        # (x_min, y_min) --> (x_min, y_max) --> (x_max, y_max) --> (x_max, y_min) --> (x_min, y_min)
        path.append((mpath.Path.MOVETO, (box[0][0], box[0][1])))
        path.append((mpath.Path.LINETO, (box[0][0], box[1][1])))
        path.append((mpath.Path.LINETO, (box[1][0], box[1][1])))
        path.append((mpath.Path.LINETO, (box[1][0], box[0][1])))
        path.append((mpath.Path.CLOSEPOLY, (box[0][0], box[0][1])))

        path_data.append((path, data[1]))

    return path_data

def assert_equal_from_matlab(a, b, options=None):
    # Compares a and b for equality. They are all going to be numpy
    # types. hdf5storage and scipy behave differently when importing
    # arrays as to whether they are 2D or not, so we will make them all
    # at least 2D regardless. For strings, the two packages produce
    # transposed results of each other, so one just needs to be
    # transposed. For object arrays, each element must be iterated over
    # to be compared. For structured ndarrays, their fields need to be
    # compared and then they can be compared element and field
    # wise. Otherwise, they can be directly compared. Note, the type is
    # often converted by scipy (or on route to the file before scipy
    # gets it), so comparisons are done by value, which is not perfect.
    a = np.atleast_2d(a)
    b = np.atleast_2d(b)
    if a.dtype.char == 'U':
        a = a.T
    if b.dtype.name == 'object':
        a = a.flatten()
        b = b.flatten()
        for index, x in np.ndenumerate(a):
            assert_equal_from_matlab(a[index], b[index], options)
    elif b.dtype.names is not None or a.dtype.names is not None:
        assert a.dtype.names is not None
        assert b.dtype.names is not None
        assert set(a.dtype.names) == set(b.dtype.names)
        a = a.flatten()
        b = b.flatten()
        for k in b.dtype.names:
            for index, x in np.ndenumerate(a):
                assert_equal_from_matlab(a[k][index], b[k][index],
                                         options)
    else:
        with warnings.catch_warnings():
            warnings.simplefilter('ignore', RuntimeWarning)
            npt.assert_equal(a, b)

    def gradient(self, x, y):
        h = 1e-9
 
        gradw = [np.zeros(w.shape) for w in self.weights]
        gradb = [np.zeros(w.shape) for w in self.weights]

        for gw, w in zip(gradw, self.weights):
            for i, v in np.ndenumerate(w.flat):
                w.flat[i] = v + h
                fxh1 = self.cost.func(self.feedforward(x), y)
                w.flat[i] = v - h
                fxh2 = self.cost.func(self.feedforward(x), y)
                gw.flat[i] = (fxh1 - fxh2)/(2*h)
                w.flat[i] = v

        for gb, b in zip(gradb, self.biases):
            for i, v in np.ndenumerate(b.flat):
                b.flat[i] = v + h
                fxh1 = self.cost.func(self.feedforward(x), y)
                b.flat[i] = v - h
                fxh2 = self.cost.func(self.feedforward(x), y)
                gb.flat[i] = (fxh1 - fxh2)/(2*h)
                b.flat[i] = v

        return gradw, gradb

def optical_flow_HS(im1, im2, s, n):
    N,M = im1.shape
    u = np.zeros(im1.shape)
    v = np.zeros(im1.shape)    
    tu = np.zeros(im1.shape)
    tv = np.zeros(im1.shape)
    fx,fy,ft = conv(im1,im2)
   
    for i in range(n):
        for (y, x), n in np.ndenumerate(im1):
            if x >= 2 and x < (M-1) and y >= 2 and y < (N-1): 
                Ex = fx[y][x]
                Ey = fy[y][x]
                Et = ft[y][x]
                AU = (u[y][x-1] + u[y][x+1] + u[y-1][x] + u[y+1][x])/4
                AV = (v[y][x-1] + v[y][x+1] + v[y-1][x] + v[y+1][x])/4
                
                A = (Ex*AU + Ey*AV +Et)
                B = (1 + s*(Ex*Ex + Ey*Ey))
                tu[y][x] = AU - (Ex*s*A/B)
                tv[y][x] = AV - (Ey*s*A/B)
        
        for (y, x), n in np.ndenumerate(im1):
            if x >= 2 and x < (M-1) and y >= 2 and y < (N-1):     
                u[y][x] = tu[y][x] 
                v[y][x] = tv[y][x]
                
    return  u,v

def plotHammingGraph(X):
    print('X:',X)
    indices = [i for i,x in np.ndenumerate(X)] # TODO arrange in circle
    coords = normalizeCoords([i for i,x in np.ndenumerate(X)])
    print('coords:',coords)
    arity = np.max(indices)+1

    from mpl_toolkits.mplot3d import Axes3D

    fig = plt.figure()
    ax = fig.add_subplot(111,projection='3d')
    s = '\n'.join(str(x) for x in np.ndenumerate(X))
    from inspect import getsource
    print('getsource:',getsource(plt.text))
    ax.set_title(Automaton(X))
    fig.text(0.01,0.99,s,horizontalalignment='left',verticalalignment='top',transform=ax.transAxes,fontsize=20)
    # add edges
    from matplotlib import patheffects
    for i,idx in enumerate(indices):
        cur_coord = coords[i]
        ax.text(cur_coord[0], cur_coord[1], cur_coord[2]-0.04,  '%s' % (str(idx)+'\n'+str(X[idx])), horizontalalignment='center', size=20, zorder=2, color='k', fontdict={'fontweight':'bold'}, path_effects=[patheffects.PathPatchEffect(edgecolor='white', facecolor='black', linewidth=1.2)])
        for j in getHammingNeighbors(idx,arity):
            ax.plot(*zip(*[cur_coord, coords[ltoi(j,arity)]]),c='black') # zorder doesn't work, sorry

    #add nodes
    ax.scatter(*zip(*coords),c=X.flatten(),s=500,cmap='gist_ncar')

def combine(a, b, w):
    matches = {}

    # split dictionaries into keys and values
    al = [x for x in a.items()]
    ak, av = zip(*al)
    bl = [x for x in b.items()]
    bk, bv = zip(*bl)

    # scale the values in the range 0-1
    a_scaled = preprocessing.minmax_scale(av, feature_range=(0,1))
    b_scaled = preprocessing.minmax_scale(bv, feature_range=(0,1))

    # build numpy structured arrays combining scaled values and original keys
    names = ['keys', 'values']
    formats = ['S225', 'f8']
    dtype = dict(names=names, formats=formats)
    anp = np.array(list(zip(ak,a_scaled)), dtype=dtype)
    bnp = np.array(list(zip(bk,b_scaled)), dtype=dtype)

    # iterate over numpy structures creating a weighted average between values with the same key
    for i, t1 in np.ndenumerate(anp):
        for j, t2 in np.ndenumerate(bnp):
            if anp['keys'][i] == bnp['keys'][j]:
                stack = np.vstack((anp['values'][i], bnp['values'][j]))
                matches[anp['keys'][i].decode("utf-8")] = np.average(stack, axis=0, weights=w)[0]   # python dictionary

    return matches

def arrays_to_xls(array_list,sheet_name_list,xls_packet_fout):
    '''Writes arrays to Excel spreadsheet.'''
    
    ibook = xlsxwriter.Workbook(xls_packet_fout,{'nan_inf_to_errors': True}) 
    
    # Write the arrays to a bundle of Excel sheets to facilitate array inspection
    for iarray,isheet_name in zip(array_list,sheet_name_list):
                                  
        if (iarray.ndim > 2):    # 3D array           
            for ilay in range(np.shape(iarray)[0]):
                jarray = iarray[ilay,:,:]
                isheet_name = isheet_name + '_' + str(ilay+1)
                isheet = ibook.add_worksheet(isheet_name)
                print 'Writing sheet: ',isheet_name
                for (irow,icol),ival in np.ndenumerate(jarray):
                    isheet.write(irow,icol,ival)
                
        else:   # For one layer        
            isheet = ibook.add_worksheet(isheet_name)
            print 'Writing sheet: ',isheet_name               
            for (irow,icol),ival in np.ndenumerate(iarray):
                isheet.write(irow,icol,ival)                                
    
    ibook.close()
    
    return

def plot_density(count_trips,count,title):
    grid = np.zeros((config.bins,config.bins))
    for (i,j),z in np.ndenumerate(grid):
        try:
            grid[j,i] = float(count[(i,j)]) / float(count_trips[(i,j)])
        except:
            grid[j,i] = 0
        #print "----"
        #print grid[i,j], i, j
        #print count[(i,j)]
        #print count_trips[(i,j)]
    grid = np.flipud(grid) #to counter matshow vertical flip
    fig, ax = plt.subplots(figsize=(10, 10))
    ax.matshow(grid, cmap='spectral')
    ax.xaxis.set_ticks_position('bottom')
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')
    xticks = np.linspace(config.minlong,config.maxlong,num=round(config.bins/2))
    yticks = np.linspace(config.minlat,config.maxlat,num=round(config.bins/2))
    yticks = yticks[::-1]
    xticks = np.around(xticks,decimals=1)
    yticks = np.around(yticks,decimals=1)
    xspace = np.linspace(0,config.bins-1,config.bins/2)
    yspace = np.linspace(0,config.bins-1,config.bins/2)
    plt.xticks(xspace,xticks)
    plt.yticks(yspace,yticks)
    for (i,j),z in np.ndenumerate(grid):
        ax.text(j, i, '{:0.2f}'.format(z), ha='center', va='center')
    plt.title(title)
    plt.show()

def make_G_hat(M, alpha=1, beta=1):
    '''G hat is Markov chain of length 2
    Gcp is a matrix to go from  contries to product and then 
    Gpc is a matrix to go from products to ccountries'''
    
    k_c  = M.sum(axis=1) #aka k_c summing over the rows
    k_p = M.sum(axis=0) #aka k_p summering over the columns
    
    G_cp = np.zeros(shape=M.shape)
    #Gcp_beta
    for [c, p], val in np.ndenumerate(M):
        numerateur = (M[c,p]) * (k_c[c] ** ((-1) * beta))
        denominateur = Gcp_denominateur(M, p, k_c, beta)
        G_cp[c,p] = numerateur / float(denominateur)
    
    
    G_pc = np.zeros(shape=M.T.shape)
    #Gpc_alpha
    for [p, c], val in np.ndenumerate(M.T):
        numerateur = (M.T[p,c]) * (k_p[p] ** ((-1) * alpha))
        denominateur = Gpc_denominateur(M, c, k_p, alpha)
        G_pc[p,c] = numerateur / float(denominateur)
    
    
    return {'G_cp': G_cp, "G_pc" : G_pc}

def trend_significance(residuals, sigma=0.05):
    nt = len(residuals)
    count = 0
    x = len(residuals[0, :, 0])
    y = len(residuals[0, 0, :])
    rcorrs = np.empty(shape=[x, y])
    for (i,j), value in np.ndenumerate(rcorrs):
        count += 1
        r_corr,_ = sp.stats.pearsonr(residuals[: -1, i, j], residuals[1:, i, j])
        if r_corr < 0:
            r_corr = 0
        rcorrs[i][j] = r_corr
    
    cs = np.empty(shape=[x, y])    
    for (i,j), rcor in np.ndenumerate(rcorrs):
        neff = float(nt * (1-rcor) / (1 + rcor))
        #neff = nt
        a = residuals[:,i,j]
        b = a * a
        d = sum(b)
        se = np.sqrt( d / ( neff - 2 ) )
        sb = se / np.sqrt( sum( ( np.arange(nt) - np.mean( np.arange(nt) ) )**2 ) )

        tcrit = sp.stats.t.isf(sigma/2.0, nt - 2 )

        c = tcrit * sb

        cs[i][j] = c
    return cs

    def OnPaint(self, event):
        dc = wx.PaintDC(self.frame.drawPanel)

        # calculate some basic variables for drawing
        boardPixelWidth, boardPixelHeight = dc.GetSize()
        pieceWidth = boardPixelWidth // self.tetris.BOARD_WIDTH
        pieceHeight = boardPixelHeight // self.tetris.BOARD_HEIGHT
        offsetX = (boardPixelWidth - pieceWidth * self.tetris.BOARD_WIDTH) // 2
        offsetY = boardPixelHeight - pieceHeight * self.tetris.BOARD_HEIGHT

        # draw board
        for (y, x), value in numpy.ndenumerate(self.tetris.board):
            if value != 0:
                self.DrawRect(dc, x, y, pieceWidth, pieceHeight, pieces.Pieces.colors[value], offsetX, offsetY)

        # draw current tile
        if self.tetris.cPiece is not None:
            for (y, x), value in numpy.ndenumerate(self.tetris.cPiece.shape):
                if value != 0:
                    self.DrawRect(
                        dc,
                        x + self.tetris.cPiece.x,
                        y + self.tetris.cPiece.y,
                        pieceWidth,
                        pieceHeight,
                        self.tetris.cPiece.color,
                        offsetX,
                        offsetY,
                    )

def label_components(m):
    linked = [{0}]
    label = 0
    b = np.zeros(m.shape,dtype=int)
    for (i,j),v in np.ndenumerate(m):
        if v:
            N = b[i-1,j]
            W = b[i,j-1]
            if N > 0:
                b[i,j] = N
                if W > 0 and N != W:
                    linked[N].add(W)
                    linked[W].add(N)
            elif W > 0:
                b[i,j] = W
            else:
                label += 1
                b[i,j] = label
                linked += [{label}]
        else:
            b[i,j] = 0

    for (i,j),v in np.ndenumerate(b):
        if v > 0 and len(linked[v])> 1:
            b[i,j] = min(linked[v])

    labels = list({min(linked[v]) for v in range(label+1)})
    for i in range(1,len(labels)):
        b[b==labels[i]] = i

    return b,len(labels)-1

def limpa_imagem(img_cinza):
    #binariza a imagem em escala de cinza
    img_bin_cinza = np.where(img_cinza < np.mean(img_cinza), 0, 255)
    
    # aplica lbp sobre a imagem em escala de cinza
    # lbp foi aplicado para evitar perda de informacao em regioes
    # proximas as regioes escuras (provaveis celulas)
    lbp_img = local_binary_pattern(img_cinza, 24, 3, method='uniform')
    
    # aplica efeito de blurring sobre a imagem resultante do lbp 
    blur_img = gaussian(lbp_img,sigma=6)    
    img_bin_blur = np.where(blur_img < np.mean(blur_img), 0, 255)
     
    # junta as duas regiões definidas pela binarizacao da imagem em escala
    # de cinza e a binarizacao do blurring    
    mascara = np.copy(img_bin_cinza)    
    for (a,b), valor in np.ndenumerate(img_bin_blur):
        if valor == 0:        
            mascara[a][b] = 0 
            
    # aplica a mascara obtida sobre a imagem original (em escala de cinza)
    # para delimitar melhor as regiões que não fornecerao informacoes (regioes
    # totalmente brancas)
    img_limpa = np.copy(img_cinza)
    for (a,b), valor in np.ndenumerate(mascara):
        if valor == 255:
            img_limpa[a][b] = 255

    return (img_limpa)

    def axis_radio_event(self):
        matrix = numpy.array([[button.isChecked() for button in row] for row in self.button_matrix])

        if (matrix.sum() < 3) :
            empty_col = matrix.sum(0).tolist().index(0)
            empty_row = matrix.sum(1).tolist().index(0)
            matrix[empty_row, empty_col] = True

            for (row, col), value in numpy.ndenumerate(matrix):
                self.button_matrix[row][col].setChecked(value)

        elif (matrix.sum() > 3) :

            double_col = matrix.sum(0).tolist().index(2)
            double_row = matrix.sum(1).tolist().index(2)
            matrix[double_row, :] = False
            matrix[:, double_col] = False
            matrix[double_row, double_col] = True

            empty_col = matrix.sum(0).tolist().index(0)
            empty_row = matrix.sum(1).tolist().index(0)
            matrix[empty_row, empty_col] = True

            for (row, col), value in numpy.ndenumerate(matrix):
                self.button_matrix[row][col].setChecked(value)

def check_gradient(x_ad, x_np, *vars):
    grad1 = x_ad.gradient(vars)

    grad2 = []
    for v in vars:
        curr = np.zeros(shape=x_np.shape)
        for i,xi in np.ndenumerate(x_np):
            curr[i] = xi.d(v)
        grad2.append(curr)

    grad1,grad2 = [np.array(g) for g in grad1,grad2]
    assert np.max(np.abs(np.log(grad1 / grad2))) < 1e-12

    if get_order() == 1:
        return

    hess1 = x_ad.hessian(vars)
    hess2 = []
    for v in vars:
        currRow1,currRow2 = [],[]
        #hess1.append(currRow1)
        hess2.append(currRow2)

        for u in vars:
            # make d2 entry for hess2
            curr = np.zeros(shape=x_np.shape)
            for i,xi in np.ndenumerate(x_np):
                if u is v:
                    curr[i] = xi.d2(v)
                else:
                    curr[i] = xi.d2c(u,v)
            currRow2.append(curr)

    hess1,hess2 = [np.array(h) for h in hess1,hess2]
    assert np.max(np.abs(np.log(hess1 / hess2))) < 1e-12

    def test_convenient_facetgrid(self):
        a = easy_array((10, 15, 4))
        d = DataArray(a, dims=["y", "x", "z"])
        g = self.plotfunc(d, x="x", y="y", col="z", col_wrap=2)

        self.assertArrayEqual(g.axes.shape, [2, 2])
        for (y, x), ax in np.ndenumerate(g.axes):
            self.assertTrue(ax.has_data())
            if x == 0:
                self.assertEqual("y", ax.get_ylabel())
            else:
                self.assertEqual("", ax.get_ylabel())
            if y == 1:
                self.assertEqual("x", ax.get_xlabel())
            else:
                self.assertEqual("", ax.get_xlabel())

        # Infering labels
        g = self.plotfunc(d, col="z", col_wrap=2)
        self.assertArrayEqual(g.axes.shape, [2, 2])
        for (y, x), ax in np.ndenumerate(g.axes):
            self.assertTrue(ax.has_data())
            if x == 0:
                self.assertEqual("y", ax.get_ylabel())
            else:
                self.assertEqual("", ax.get_ylabel())
            if y == 1:
                self.assertEqual("x", ax.get_xlabel())
            else:
                self.assertEqual("", ax.get_xlabel())

def inverse_relation(R, size_rhs=None, with_indices=False):
    """Computes the inverse relation of a relation.

    If `r` is a relation, then the inverse relation `ri` is defined by

        x ri y  <=>  y r x

    Parameters
    ----------
    R
        2D |NumPy array| of integers representing a relation r on the
        natural numbers via ::

            x r y <=> (x < R.size[0] and y in R[x]).

        Rows of `R` which are to short are padded with -1.
    size_rhs
        Can be provided for speedup. Has to be greater than `R.max()`.
    with_indices
        If `True`, also return the matrix `RINVI`.

    Returns
    -------
    RINV
        2D |NumPy array| representation of the inverse relation.
    RINVI
        |NumPy array| such that for `RINV[i, j] != -1`::

            R[RINV[i, j], RINVI[i, j]] = i.

        Only returned if `with_indices` is `True`.
    """

    assert R.ndim == 2
    logger.warn('Call to unoptimized function inverse_relation')

    num_columns_RINV = np.bincount(R.ravel()).max()
    size_rhs = size_rhs or (R.max() + 1)
    RINV = np.empty((size_rhs, num_columns_RINV), dtype=R.dtype)
    RINV.fill(-1)
    if with_indices:
        RINVI = np.empty_like(RINV)
        RINVI.fill(-1)

    RINV_COL_COUNTS = np.zeros(size_rhs, dtype=np.int32)

    if not with_indices:
        for index, x in np.ndenumerate(R):
            if x >= 0:
                RINV[x, RINV_COL_COUNTS[x]] = index[0]
                RINV_COL_COUNTS[x] += 1
        return RINV
    else:
        for index, x in np.ndenumerate(R):
            if x >= 0:
                RINV[x, RINV_COL_COUNTS[x]] = index[0]
                RINVI[x, RINV_COL_COUNTS[x]] = index[1]
                RINV_COL_COUNTS[x] += 1
        return RINV, RINVI