def VtuMatchLocationsArbitrary(vtu1, vtu2, tolerance=1.0e-6):
    """
  Check that the locations in the supplied vtus match, returning True if they
  match and False otherwise.
  The locations may be in a different order.
  """

    locations1 = vtu1.GetLocations()
    locations2 = vtu2.GetLocations()
    if not locations1.shape == locations2.shape:
        return False

    epsilon = numpy.ones(locations1.shape[1]) * numpy.finfo(numpy.float).eps
    for j in range(locations1.shape[1]):
        epsilon[j] = epsilon[j] * (locations1[:, j].max() - locations1[:, j].min())

    for i in range(len(locations1)):
        for j in range(len(locations1[i])):
            if abs(locations1[i][j]) < epsilon[j]:
                locations1[i][j] = 0.0
            if abs(locations2[i][j]) < epsilon[j]:
                locations2[i][j] = 0.0

    # lexical sort on x,y and z coordinates resp. of locations1 and locations2
    sort_index1 = numpy.lexsort(locations1.T)
    sort_index2 = numpy.lexsort(locations2.T)

    # should now be in same order, so we can check for its biggest difference
    return abs(locations1[sort_index1] - locations2[sort_index2]).max() < tolerance

    def test_sort_index_multicolumn(self):
        import random
        A = np.arange(5).repeat(20)
        B = np.tile(np.arange(5), 20)
        random.shuffle(A)
        random.shuffle(B)
        frame = DataFrame({'A': A, 'B': B,
                           'C': np.random.randn(100)})

        # use .sort_values #9816
        with tm.assert_produces_warning(FutureWarning):
            frame.sort_index(by=['A', 'B'])
        result = frame.sort_values(by=['A', 'B'])
        indexer = np.lexsort((frame['B'], frame['A']))
        expected = frame.take(indexer)
        assert_frame_equal(result, expected)

        # use .sort_values #9816
        with tm.assert_produces_warning(FutureWarning):
            frame.sort_index(by=['A', 'B'], ascending=False)
        result = frame.sort_values(by=['A', 'B'], ascending=False)
        indexer = np.lexsort((frame['B'].rank(ascending=False),
                              frame['A'].rank(ascending=False)))
        expected = frame.take(indexer)
        assert_frame_equal(result, expected)

        # use .sort_values #9816
        with tm.assert_produces_warning(FutureWarning):
            frame.sort_index(by=['B', 'A'])
        result = frame.sort_values(by=['B', 'A'])
        indexer = np.lexsort((frame['A'], frame['B']))
        expected = frame.take(indexer)
        assert_frame_equal(result, expected)

def test_cutting_plane_selector():
    # generate fake data with a number of non-cubical grids
    ds = fake_random_ds(64, nprocs=51)
    assert all(ds.periodicity)

    # test cutting plane against orthogonal plane
    for i, d in enumerate("xyz"):
        norm = np.zeros(3)
        norm[i] = 1.0

        for coord in np.arange(0, 1.0, 0.1):
            center = np.zeros(3)
            center[i] = coord

            data = ds.slice(i, coord)
            data.get_data()
            data2 = ds.cutting(norm, center)
            data2.get_data()

            assert data.shape[0] == data2.shape[0]

            cells1 = np.lexsort((data["x"], data["y"], data["z"]))
            cells2 = np.lexsort((data2["x"], data2["y"], data2["z"]))
            for d2 in "xyz":
                yield assert_equal, data[d2][cells1], data2[d2][cells2]

    def negative_gradient(self, y_true, y_pred, sample_group=None, **kargs):
        y_pred = y_pred.ravel()
        # the lambda terms
        grad = np.empty_like(y_true, dtype=np.float64)

        # for updating terminal regions
        self.weights = np.empty_like(y_true, dtype=np.float64)

        if sample_group is None:
            ix = np.lexsort((y_true, -y_pred))
            inv_ix = np.empty_like(ix)
            inv_ix[ix] = np.arange(len(ix))
            tmp_grad, tmp_weights = _lambda(y_true[ix], y_pred[ix],
                                            self.max_rank)
            grad = tmp_grad[inv_ix]
            self.weights = tmp_weights[inv_ix]
        else:
            for start, end in self._groupby(sample_group):
                ix = np.lexsort((y_true[start:end], -y_pred[start:end]))
                inv_ix = np.empty_like(ix)
                inv_ix[ix] = np.arange(len(ix))

                # sort by current score before passing
                # and then remap the return values
                tmp_grad, tmp_weights = _lambda(y_true[ix + start],
                                                y_pred[ix + start],
                                                self.max_rank)
                grad[start:end] = tmp_grad[inv_ix]
                self.weights[start:end] = tmp_weights[inv_ix]

        return grad

 def set_sort_col(self, col_index, add=False):
     '''Set the column to sort this table by. If add is true, this column
     will be added to the end of the existing sort order (or removed from the
     sort order if it is already present.)
     '''
     if not add:
         if len(self.sortcols)>0 and col_index in self.sortcols:
             # If this column is already sorted, flip it
             self.row_order = self.row_order[::-1]
             self.sortdir = -self.sortdir
         else:
             self.sortdir = 1
             self.sortcols = [col_index]
             # If this column hasn't been sorted yet, then sort descending
             self.row_order = np.lexsort(self.data[:,self.col_order][:,self.sortcols[::-1]].T.tolist())
     else:
         if len(self.sortcols)>0 and col_index in self.sortcols:
             self.sortcols.remove(col_index)
         else:
             self.sortcols += [col_index]
         if self.sortcols==[]:
             # if all sort columns have been toggled off, reset row_order
             self.row_order = np.arange(self.data.shape[0])
         else:
             self.row_order = np.lexsort(self.data[:,self.sortcols[::-1]].T.tolist())
     self.ordered_data = self.data[self.row_order,:][:,self.col_order]

def preCompute(rowBased_row_array,rowBased_col_array,S_rowBased_data_array):
    """
    format affinity/similarity matrix
    """
    
    # Get parameters
    data_len=len(S_rowBased_data_array)
    row_indptr=sparseAP_cy.getIndptr(rowBased_row_array)
    if row_indptr[-1]!=data_len: row_indptr=np.concatenate((row_indptr,np.array([data_len])))
    row_to_col_ind_arr=np.lexsort((rowBased_row_array,rowBased_col_array))
    colBased_row_array=sparseAP_cy.npArrRearrange_int_para(rowBased_row_array,row_to_col_ind_arr)
    colBased_col_array=sparseAP_cy.npArrRearrange_int_para(rowBased_col_array,row_to_col_ind_arr)
    col_to_row_ind_arr=np.lexsort((colBased_col_array,colBased_row_array))
    col_indptr=sparseAP_cy.getIndptr(colBased_col_array)
    if col_indptr[-1]!=data_len: col_indptr=np.concatenate((col_indptr,np.array([data_len])))
    kk_col_index=sparseAP_cy.getKKIndex(colBased_row_array,colBased_col_array)
    
    #Initialize matrix A, R
    A_rowbased_data_array=np.array([0.0]*data_len)
    R_rowbased_data_array=np.array([0.0]*data_len)
    
    #Add random samll value to remove degeneracies
    random_state=np.random.RandomState(0)
    S_rowBased_data_array+=1e-12*random_state.randn(data_len)*(np.amax(S_rowBased_data_array)-np.amin(S_rowBased_data_array))
    
    #Convert row_to_col_ind_arr/col_to_row_ind_arr data type to np.int datatype so it is compatible with cython code
    row_to_col_ind_arr=row_to_col_ind_arr.astype(np.int)
    col_to_row_ind_arr=col_to_row_ind_arr.astype(np.int)
    
    return S_rowBased_data_array, A_rowbased_data_array, R_rowbased_data_array,col_indptr,row_indptr,row_to_col_ind_arr,col_to_row_ind_arr,kk_col_index

    def fullCheck(self,a):
        # Check that atoms repeats over a
        bp1, bp2, bp3 = self.bp1, self.bp2, self.bp3 
        atomtypes = N.unique(self.snr)
        passed = True
        for atomtype in atomtypes:
            sameatoms = N.argwhere(self.snr==atomtype)
            samexyz = self.xyz[sameatoms]
            samexyz = samexyz.reshape((-1, 3))
            shifted = samexyz+a.reshape((1,3))
            
            # Move inside PBC
            samexyz = moveIntoCell(samexyz,self.pbc[0,:],self.pbc[1,:],self.pbc[2,:],self.accuracy)
            shifted = moveIntoCell(shifted,self.pbc[0,:],self.pbc[1,:],self.pbc[2,:],self.accuracy)

            # Should be the same if sorted!
            ipiv = N.lexsort(N.round(N.transpose(samexyz)/self.accuracy))
            samexyz = samexyz[ipiv,:]
            ipiv = N.lexsort(N.round(N.transpose(shifted)/self.accuracy))
            shifted = shifted[ipiv,:]
            
            if not N.allclose(samexyz,shifted,atol=self.accuracy):
                passed = False
            
        return passed

def loc_vector_labels(x):
    """Identify unique labels from the vector of image labels
    
    x - a vector of one label or dose per image
    
    returns labels, labnum, uniqsortvals
    labels - a vector giving an ordinal per image where that ordinal
             is an index into the vector of unique labels (uniqsortvals)
    labnum - # of unique labels in x
    uniqsortvals - a vector containing the unique labels in x
    """
    #
    # Get the index of each image's label in the sorted array
    #
    order = np.lexsort((x,))
    reverse_order = np.lexsort((order,))
    #
    # Get a sorted view of the labels
    #
    sorted_x = x[order]
    #
    # Find the elements that start a new run of labels in the sorted array
    # ex: 0,0,0,3,3,3,5,5,5
    #     1,0,0,1,0,0,1,0,0
    #
    # Then cumsum - 1 turns into:
    #     0,0,0,1,1,1,2,2,2
    #
    # and sorted_x[first_occurrence] gives the unique labels in order
    first_occurrence = np.ones(len(x), bool)
    first_occurrence[1:] = sorted_x[:-1] != sorted_x[1:]
    sorted_labels = np.cumsum(first_occurrence) - 1
    labels = sorted_labels[reverse_order]
    uniqsortvals = sorted_x[first_occurrence]
    return (labels, len(uniqsortvals), uniqsortvals)

    def test_03_01_graph(self):
        '''Make a simple graph'''
        #
        # The skeleton looks something like this:
        #
        #   .   .
        #    . .
        #     .
        #     .
        i,j = np.mgrid[-10:11,-10:11]
        skel = (i < 0) & (np.abs(i) == np.abs(j))
        skel[(i >= 0) & (j == 0)] = True
        #
        # Put a single label at the bottom
        #
        labels = np.zeros(skel.shape, int)
        labels[(i > 8) & (np.abs(j) < 2)] = 1
        np.random.seed(31)
        intensity = np.random.uniform(size = skel.shape)
        workspace, module = self.make_workspace(
            labels, skel, intensity_image = intensity, wants_graph = True)
        module.run(workspace)
        edge_graph = self.read_graph_file(EDGE_FILE)
        vertex_graph = self.read_graph_file(VERTEX_FILE)
        vidx = np.lexsort((vertex_graph["j"], vertex_graph["i"]))
        #
        # There should be two vertices at the bottom of the array - these
        # are bogus artifacts of the object hitting the edge of the image
        #
        for vidxx in vidx[-2:]:
            self.assertEqual(vertex_graph["i"][vidxx], 20)
        vidx = vidx[:-2]
        
        expected_vertices = ((0,0), (0,20), (10,10), (17,10))
        self.assertEqual(len(vidx), len(expected_vertices))
        for idx, v in enumerate(expected_vertices):
            vv = vertex_graph[vidx[idx]]
            self.assertEqual(vv["i"], v[0])
            self.assertEqual(vv["j"], v[1])

        #
        # Get rid of edges to the bogus vertices
        #
        for v in ("v1","v2"):
            edge_graph = edge_graph[vertex_graph["i"][edge_graph[v]-1] != 20]
            
        eidx = np.lexsort((vertex_graph["j"][edge_graph["v1"]-1],
                           vertex_graph["i"][edge_graph["v1"]-1],
                           vertex_graph["j"][edge_graph["v2"]-1],
                           vertex_graph["i"][edge_graph["v2"]-1]))
        expected_edges = (((0,0),(10,10),11, np.sum(intensity[(i <= 0) & (j<=0) & skel])),
                          ((0,20),(10,10),11, np.sum(intensity[(i <= 0) & (j>=0) & skel])),
                          ((10,10),(17,10),8, np.sum(intensity[(i >= 0) & (i <= 7) & skel])))
        for i, (v1, v2, length, total_intensity) in enumerate(expected_edges):
            ee = edge_graph[eidx[i]]
            for ve, v in ((v1, ee["v1"]), (v2, ee["v2"])):
                self.assertEqual(ve[0], vertex_graph["i"][v-1])
                self.assertEqual(ve[1], vertex_graph["j"][v-1])
            self.assertEqual(length, ee["length"])
            self.assertAlmostEqual(total_intensity, ee["total_intensity"], 4)

def sortrowsByMultiCol(data, list_sortCol, isAscending=1):
    """
    Sort 2D numpy array by multiple columns
    :param data:
    :param list_sortCol: e.g. [0, 2, 3], 1st element is most important column to sort, 2nd element is 2nd most important, etc.
    :return: sorted data
    """
    # data_sort = data # data_sort has to be number
    # k_col = len(list_sortCol)-1
    # while k_col >= 0:
    #     idx_col = list_sortCol[k_col]
    #     t_data = np.float64(data_sort[:,idx_col])
    #     if isAscending==1:
    #         data_sort = data_sort[t_data.argsort(),:]
    #     elif isAscending==0:
    #         data_sort = data_sort[(-t_data).argsort(),:] # descending order
    #     k_col -= 1

    dataForSort = np.transpose(data[:, list_sortCol[::-1]])  # [start:stop:step] -1 is meant to reverse the order
    if isAscending == 1:
        idx_sort = np.lexsort(dataForSort)  # last row in dataForSort is most important to sort
    elif isAscending == 0:
        idx_sort = np.lexsort(-dataForSort)  # last row in dataForSort is most important to sort

    data_sorted = data[idx_sort, :]

    return data_sorted

def pearson_correlation(movie_id_1, movie_id_2):
    rated_movie1 = ratings[ratings[:, 1] == movie_id_1]
    rated_movie2 = ratings[ratings[:, 1] == movie_id_2]
    rated_both = np.intersect1d(rated_movie1[:, 0], rated_movie2[:, 0], True)

    if len(rated_both) < 15:
        return 0

    ratings_movie1 = rated_movie1[np.in1d(rated_movie1[:, 0], rated_both), :]
    ratings_movie2 = rated_movie2[np.in1d(rated_movie2[:, 0], rated_both), :]
    sorted_movie1 = ratings_movie1[np.lexsort((ratings_movie1[:, 0], ))][:, [0, 2]]
    sorted_movie2 = ratings_movie2[np.lexsort((ratings_movie2[:, 0], ))][:, [0, 2]]

    mean1 = np.mean(ratings_movie1[:, 2])
    mean2 = np.mean(ratings_movie2[:, 2])

    numerator = 0
    denomX = 0
    denomY = 0
    for i in range(len(sorted_movie1)):
        x = sorted_movie1[i][1] - mean1
        y = sorted_movie2[i][1] - mean2
        numerator += x * y
        denomX += x * x
        denomY += y * y

    if (denomX == 0 or denomY == 0):
        return 0

    return round(numerator / m.sqrt(denomX * denomY), 3)

def test_hemisphere_subdivide():

    def flip(vertices):
        x, y, z = vertices.T
        f = (z < 0) | ((z == 0) & (y < 0)) | ((z == 0) & (y == 0) & (x < 0))
        return 1 - 2*f[:, None]

    decimals = 6
    # Test HemiSphere.subdivide
    # Create a hemisphere by dividing a hemi-icosahedron
    hemi1 = HemiSphere.from_sphere(unit_icosahedron).subdivide(4)
    vertices1 = np.round(hemi1.vertices, decimals)
    vertices1 *= flip(vertices1)
    order = np.lexsort(vertices1.T)
    vertices1 = vertices1[order]

    # Create a hemisphere from a subdivided sphere
    sphere = unit_icosahedron.subdivide(4)
    hemi2 = HemiSphere.from_sphere(sphere)
    vertices2 = np.round(hemi2.vertices, decimals)
    vertices2 *= flip(vertices2)
    order = np.lexsort(vertices2.T)
    vertices2 = vertices2[order]

    # The two hemispheres should have the same vertices up to their order
    nt.assert_array_equal(vertices1, vertices2)

    # Create a hemisphere from vertices
    hemi3 = HemiSphere(xyz=hemi1.vertices)
    nt.assert_array_equal(hemi1.faces, hemi3.faces)
    nt.assert_array_equal(hemi1.edges, hemi3.edges)

def sortContours2( contours, direction = "x" ):#TODO
    contourPoints = np.zeros((len(contours),2), dtype = int)
    if direction == "x":
        a = 1
        b = 0
    elif direction == "y":
        a = 0
        b = 1
    
    counter = 0
    for cnt in contours:
        conResh = np.reshape(cnt,(-1,2))
        idx = np.lexsort( (conResh[:,a],conResh[:,b]) )
        sortedContours = conResh[idx,:]
        contourPoints[counter,:] = sortedContours[0,:] # The coordinate of reference point.
        counter = counter + 1
        
    sortedIdx = np.lexsort((contourPoints[:,a], contourPoints[:,b]))
    sortedContours = []
    referencePoints = []
    for idx in sortedIdx:
        sortedContours.append(contours[idx])
        referencePoints.append(contourPoints[idx])
        
    return sortedContours, referencePoints

def doubleParetoSorting(x0, x1):
    fronts = [[]]
    left = [[]]
    right = [[]]
    idx = np.lexsort((x1, x0))
    
    idxEdge = np.lexsort((-np.square(x0-0.5), x1))
    
    fronts[-1].append(idxEdge[0])
    left[-1].append(x0[idxEdge[0]])
    right[-1].append(x0[idxEdge[0]])
    for i0 in idxEdge[1:]:
        if x0[i0]>=left[-1] and x0[i0]<=right[-1]:
            #add a new front
            fronts.append([])
            left.append([])
            right.append([])
            fronts[-1].append(i0)
            left[-1].append(x0[i0])
            right[-1].append(x0[i0])
        else:
            #check existing fonts
            for i1 in range(len(fronts)):
                if x0[i0]<left[i1] or x0[i0]>right[i1]:
                    if x0[i0]<left[i1]:
                        left[i1] = x0[i0]
                        fronts[i1].insert(0, i0)
                    else:
                        right[i1] = x0[i0]
                        fronts[i1].append(i0)
                    break    
    return (fronts, idx)

def ssea_ranker (ssea_list, q_val_cutoff, peak_type):
    if peak_type == "Amplification":
        sign = 1
    elif peak_type == "Deletion":
        sign = -1
    
    thresh=(ssea_list[:,1]<q_val_cutoff).astype(int)
    sign_adj = ssea_list*sign
    rank_big=np.lexsort((thresh, -sign_adj[:,0]))
    ranks_ranked = []
    for x in enumerate(rank_big):
        ranks_ranked.append(x)
    ranked = np.array(ranks_ranked)
    
    ranked_sorted = ranked[np.lexsort([ranked[:,0], ranked[:,1]])]
    ranks = []
    for x in ranked_sorted:
        ranks.append(x[0])
    
    
    ranks_out = []
    for x in xrange(len(ssea_list[:,0])):
        if thresh[x]==1: 
            ranks_out.append(ranks[x])
        else:
            ranks_out.append('not significant')
    return ranks_out

def compute_orderings(path):

    logging.info("Computing orderings of features")
    job = load_job(path)
    original = np.arange(len(job.input.feature_ids))
    stats = job.results.feature_to_score[...]
    rev_stats = 0.0 - stats

    logging.info("  Computing ordering by score for each tuning param")
    by_score_original = np.zeros(np.shape(job.results.raw_stats), int)
    for i in range(len(job.settings.tuning_params)):
        by_score_original[i] = np.lexsort(
            (original, rev_stats[i]))

    order_by_score_original = by_score_original

    logging.info("  Computing ordering by fold change")
    by_foldchange_original = np.zeros(np.shape(job.results.fold_change.table), int)
    foldchange = job.results.fold_change.table[...]
    rev_foldchange = 0.0 - foldchange
    for i in range(len(job.results.fold_change.header)):
        keys = (original, rev_foldchange[..., i])

        by_foldchange_original[..., i] = np.lexsort(keys)

    order_by_foldchange_original = by_foldchange_original

    with h5py.File(path, 'r+') as db:
        orderings = db.create_group('orderings')
        orderings['by_score_original'] = order_by_score_original
        orderings['by_foldchange_original'] = order_by_foldchange_original        

def VtuMatchLocationsArbitrary(vtu1, vtu2, tolerance = 1.0e-6):
  """
  Check that the locations in the supplied vtus match, returning True if they
  match and False otherwise.
  The locations may be in a different order.
  """
   
  locations1 = vtu1.GetLocations()
  locations2 = vtu2.GetLocations()
  if not locations1.shape == locations2.shape:
    return False   
    
  for j in range(locations1.shape[1]):
    # compute the smallest possible precision given the range of this coordinate
    epsilon = numpy.finfo(numpy.float).eps * numpy.abs(locations1[:,j]).max()
    if tolerance<epsilon:
      # the specified tolerance is smaller than possible machine precision
      # (or something else went wrong)
      raise Exception("ERROR: specified tolerance is smaller than machine precision of given locations")
    # ensure epsilon doesn't get too small (might be for zero for instance)
    epsilon=max(epsilon,tolerance/100.0)

    # round to that many decimal places (-2 to be sure) so that
    # we don't get rounding issues with lexsort
    locations1[:,j]=numpy.around(locations1[:,j], int(-numpy.log10(epsilon))-2)
    locations2[:,j]=numpy.around(locations2[:,j], int(-numpy.log10(epsilon))-2)

  # lexical sort on x,y and z coordinates resp. of locations1 and locations2
  sort_index1=numpy.lexsort(locations1.T)
  sort_index2=numpy.lexsort(locations2.T)
  
  # should now be in same order, so we can check for its biggest difference
  return numpy.allclose(locations1[sort_index1],locations2[sort_index2], atol=tolerance)

def test_set_synaptic_parameters_fully_connected(sim):
    sim.setup()
    mpi_rank = sim.rank()
    p1 = sim.Population(4, sim.IF_cond_exp())
    p2 = sim.Population(2, sim.IF_cond_exp())
    syn = sim.TsodyksMarkramSynapse(U=0.5, weight=0.123, delay=0.1)
    prj = sim.Projection(p1, p2, sim.AllToAllConnector(), syn)

    expected = numpy.array([
        (0.0, 0.0, 0.123, 0.1, 0.5),
        (0.0, 1.0, 0.123, 0.1, 0.5),
        (1.0, 0.0, 0.123, 0.1, 0.5),
        (1.0, 1.0, 0.123, 0.1, 0.5),
        (2.0, 0.0, 0.123, 0.1, 0.5),
        (2.0, 1.0, 0.123, 0.1, 0.5),
        (3.0, 0.0, 0.123, 0.1, 0.5),
        (3.0, 1.0, 0.123, 0.1, 0.5),
    ])
    actual = numpy.array(prj.get(['weight', 'delay', 'U'], format='list'))
    if mpi_rank == 0:
        ind = numpy.lexsort((actual[:, 1], actual[:, 0]))
        assert_arrays_almost_equal(actual[ind], expected, 1e-16)

    positional_weights = numpy.array([[0, 1], [2, 3], [4, 5], [6, 7]], dtype=float)
    prj.set(weight=positional_weights)
    expected = positional_weights
    actual = prj.get('weight', format='array')
    if mpi_rank == 0:
        assert_arrays_equal(actual, expected)

    u_list = [0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2]
    prj.set(U=u_list)
    expected = numpy.array([[0.9, 0.8], [0.7, 0.6], [0.5, 0.4], [0.3, 0.2]])
    actual = prj.get('U', format='array')
    if mpi_rank == 0:
        assert_arrays_equal(actual, expected)

    f_delay = lambda d: 0.5+d
    prj.set(delay=f_delay)
    expected = numpy.array([[0.5, 1.5], [1.5, 0.5], [2.5, 1.5], [3.5, 2.5]])
    actual = prj.get('delay', format='array')
    if mpi_rank == 0:
        assert_arrays_equal(actual, expected)

    # final sanity check
    expected = numpy.array([
        (0.0, 0.0, 0.0, 0.5, 0.9),
        (0.0, 1.0, 1.0, 1.5, 0.8),
        (1.0, 0.0, 2.0, 1.5, 0.7),
        (1.0, 1.0, 3.0, 0.5, 0.6),
        (2.0, 0.0, 4.0, 2.5, 0.5),
        (2.0, 1.0, 5.0, 1.5, 0.4),
        (3.0, 0.0, 6.0, 3.5, 0.3),
        (3.0, 1.0, 7.0, 2.5, 0.2),
    ])
    actual = numpy.array(prj.get(['weight', 'delay', 'U'], format='list'))
    if mpi_rank == 0:
        ind = numpy.lexsort((actual[:, 1], actual[:, 0]))
        assert_arrays_equal(actual[ind], expected)

def _cplxreal(z, tol=None):

    import numpy as np
    from numpy import atleast_1d, atleast_2d, array

    z = atleast_1d(z)
    if z.size == 0:
        return z, z
    elif z.ndim != 1:
        raise ValueError('_cplxreal only accepts 1D input')

    if tol is None:
        # Get tolerance from dtype of input
        tol = 100 * np.finfo((1.0 * z).dtype).eps

    # Sort by real part, magnitude of imaginary part (speed up further sorting)
    z = z[np.lexsort((abs(z.imag), z.real))]

    # Split reals from conjugate pairs
    real_indices = abs(z.imag) <= tol * abs(z)
    zr = z[real_indices].real

    if len(zr) == len(z):
        # Input is entirely real
        return array([]), zr

    # Split positive and negative halves of conjugates
    z = z[~real_indices]
    zp = z[z.imag > 0]
    zn = z[z.imag < 0]

    if len(zp) != len(zn):
        raise ValueError('Array contains complex value with no matching '
                         'conjugate.')

    # Find runs of (approximately) the same real part
    same_real = np.diff(zp.real) <= tol * abs(zp[:-1])
    diffs = numpy.diff(concatenate(([0], same_real, [0])))
    run_starts = numpy.where(diffs > 0)[0]
    run_stops = numpy.where(diffs < 0)[0]

    # Sort each run by their imaginary parts
    for i in range(len(run_starts)):
        start = run_starts[i]
        stop = run_stops[i] + 1
        for chunk in (zp[start:stop], zn[start:stop]):
            chunk[...] = chunk[np.lexsort([abs(chunk.imag)])]

    # Check that negatives match positives
    if any(abs(zp - zn.conj()) > tol * abs(zn)):
        raise ValueError('Array contains complex value with no matching '
                         'conjugate.')

    # Average out numerical inaccuracy in real vs imag parts of pairs
    zc = (zp + zn.conj()) / 2

    return zc, zr

def convert(cst_file_in, scalar_file_out):
    """
    Calculates a scalar element pattern file from a CST element pattern file

    Parameters
    ----------
    cst_file_in : string
        Input CST format element pattern file
    scalar_file_out : string
        Output scalar format element pattern file

    Notes
    -----
    This function is designed to be used to create scalar element input files
    for the oskar_fit_element_data application.
    """

    import numpy as np

    # Load the CST element pattern data for X. (Ignore lines that don't consist
    # of 8 floats)
    X = load_cst_file(cst_file_in)
    # Only require a columns for:
    # Theta, Phi, Abs(Theta), Phase(Theta), Abs(Phi), Phase(Phi)
    X = np.copy(X[:, [0, 1, 3, 4, 5, 6]])

    # Generate the rotated data for Y from X by adding 90 degrees to the phi
    # values
    Y = np.copy(X)
    Y[:, 1] += 90.0
    Y[Y[:, 1] >= 360.0, 1] -= 360.0

    # Linked column sort by phi and then theta for both X and Y.
    X = X[np.lexsort((X[:, 1], X[:, 0])), :]
    Y = Y[np.lexsort((Y[:, 1], Y[:, 0])), :]
    assert(np.sum(X[:, 0] == Y[:, 0]) == len(X[:, 0]))
    assert(np.sum(X[:, 1] == Y[:, 1]) == len(X[:, 1]))

    # Generate scalar values from sorted data.
    X_theta = X[:, 2] * np.exp(1j * X[:, 3] * np.pi / 180.0)
    X_phi = X[:, 4] * np.exp(1j * X[:, 5] * np.pi / 180.0)
    Y_theta = Y[:, 2] * np.exp(1j * Y[:, 3] * np.pi / 180.0)
    Y_phi = Y[:, 4] * np.exp(1j * Y[:, 5] * np.pi / 180.0)
    s = X_theta * np.conj(X_theta) + X_phi * np.conj(X_phi) + \
        Y_theta * np.conj(Y_theta) + Y_phi * np.conj(Y_phi)

    # Take the sqrt to convert to a 'voltage'
    s = np.sqrt(0.5 * s)
    s_amp = np.absolute(s)
    s_phase = np.angle(s, deg=True)

    # Write scalar values to file Columns = (theta, phi, amp, phase).
    o = np.column_stack((X[:, 0], X[:, 1], s_amp, s_phase))
    np.savetxt(scalar_file_out, o, fmt=['%12.4f', '%12.4f', '%20.6e',
                                        '%12.4f'])