def test_drop_fields(self):
        # Test drop_fields
        a = np.array([(1, (2, 3.0)), (4, (5, 6.0))], dtype=[("a", int), ("b", [("ba", float), ("bb", int)])])

        # A basic field
        test = drop_fields(a, "a")
        control = np.array([((2, 3.0),), ((5, 6.0),)], dtype=[("b", [("ba", float), ("bb", int)])])
        assert_equal(test, control)

        # Another basic field (but nesting two fields)
        test = drop_fields(a, "b")
        control = np.array([(1,), (4,)], dtype=[("a", int)])
        assert_equal(test, control)

        # A nested sub-field
        test = drop_fields(a, ["ba"])
        control = np.array([(1, (3.0,)), (4, (6.0,))], dtype=[("a", int), ("b", [("bb", int)])])
        assert_equal(test, control)

        # All the nested sub-field from a field: zap that field
        test = drop_fields(a, ["ba", "bb"])
        control = np.array([(1,), (4,)], dtype=[("a", int)])
        assert_equal(test, control)

        test = drop_fields(a, ["a", "b"])
        assert_(test is None)

def _build_trajectories(data):
# """
# build_trajectories(data) is responsible for the book keeping of the trajectories, 
# using the prev,next fields in the frames,creating new set of data, in the form of trajectories
# """
	# first frame, initialization
	trajid = 0 # running trajid counter
	frame = data[0]
	frame.trajid = n.nan*n.empty_like(frame.x)
	ind = frame.next > -2
	frame.trajid[ind] = range(trajid,trajid+ind.size)
	trajid = trajid+ind.size


	for i,frame in data[1:]:
		frame.trajid = n.nan*n.empty_like(frame.x)
		old = frame.prev > -1
		frame.trajid[old] = data[i-1].trajid[frame.prev[old]]
		ind = frame.prev < 0 and frame.next > -2
		frame.trajid[ind] = range(trajid,trajid+ind.size)
		trajid = trajid+ind.size
		drop_fields(frame,['prev','next'])
		
	for frame in data:
		frame = frame[~n.isnan(frame)]
		frame.t = frame.t*n.ones_like(frame.x)
	
	return data

    def test_drop_fields(self):
        # Test drop_fields
        a = np.array([(1, (2, 3.0)), (4, (5, 6.0))],
                     dtype=[('a', int), ('b', [('ba', float), ('bb', int)])])

        # A basic field
        test = drop_fields(a, 'a')
        control = np.array([((2, 3.0),), ((5, 6.0),)],
                           dtype=[('b', [('ba', float), ('bb', int)])])
        assert_equal(test, control)

        # Another basic field (but nesting two fields)
        test = drop_fields(a, 'b')
        control = np.array([(1,), (4,)], dtype=[('a', int)])
        assert_equal(test, control)

        # A nested sub-field
        test = drop_fields(a, ['ba', ])
        control = np.array([(1, (3.0,)), (4, (6.0,))],
                           dtype=[('a', int), ('b', [('bb', int)])])
        assert_equal(test, control)

        # All the nested sub-field from a field: zap that field
        test = drop_fields(a, ['ba', 'bb'])
        control = np.array([(1,), (4,)], dtype=[('a', int)])
        assert_equal(test, control)

        test = drop_fields(a, ['a', 'b'])
        assert_(test is None)

    def get_raw_chamber_data(self,filtered_data):
        # chamber_dtype = numpy.dtype([('time_secs', '<u4'),
        #                              ('time_nsecs', '<u4'),
        #                              ('time_rel', '<f4'),
        #                              ('status', '|S25'),
        #                              ('tunnel', '<u2'),
        #                              ('fly_x', '<f4'),
        #                              ('fly_y', '<f4'),
        #                              ('fly_angle', '<f4'),
        #                              ])
        header = list(FILE_TOOLS.chamber_dtype.names)
        tracking_chamber_data = filtered_data[filtered_data['status'] != 'Walk To End']
        tracking_chamber_data = tracking_chamber_data[header]
        tracking_chamber_data = tracking_chamber_data.astype(FILE_TOOLS.chamber_dtype)
        tracking_chamber_data['tunnel'] = tracking_chamber_data['tunnel']+1
        indicies = tracking_chamber_data['status'] == 'End Chamber Ethanol'
        raw_chamber_data_ethanol = tracking_chamber_data[indicies]
        raw_chamber_data_ethanol = recfunctions.drop_fields(raw_chamber_data_ethanol,
                                                            'status',
                                                            usemask=False)
        status_array = numpy.array(['Ethanol']*len(raw_chamber_data_ethanol),dtype='|S25')
        raw_chamber_data_ethanol = recfunctions.append_fields(raw_chamber_data_ethanol,
                                                              'status',
                                                              status_array,
                                                              dtypes='|S25',
                                                              usemask=False)
        raw_chamber_data = raw_chamber_data_ethanol

        ethanol_start_time = raw_chamber_data_ethanol['time_rel'][0]
        indicies = tracking_chamber_data['status'] == 'End Chamber Air'
        indicies &= tracking_chamber_data['time_rel'] < ethanol_start_time
        raw_chamber_data_air_before = tracking_chamber_data[indicies]
        raw_chamber_data_air_before = recfunctions.drop_fields(raw_chamber_data_air_before,
                                                               'status',
                                                               usemask=False)
        status_array = numpy.array(['AirBefore']*len(raw_chamber_data_air_before),dtype='|S25')
        raw_chamber_data_air_before = recfunctions.append_fields(raw_chamber_data_air_before,
                                                                 'status',
                                                                 status_array,
                                                                 dtypes='|S25',
                                                                 usemask=False)
        raw_chamber_data = recfunctions.stack_arrays((raw_chamber_data_air_before,raw_chamber_data),usemask=False)

        indicies = tracking_chamber_data['status'] == 'End Chamber Air'
        indicies &= tracking_chamber_data['time_rel'] > ethanol_start_time
        raw_chamber_data_air_after = tracking_chamber_data[indicies]
        raw_chamber_data_air_after = recfunctions.drop_fields(raw_chamber_data_air_after,
                                                               'status',
                                                               usemask=False)
        status_array = numpy.array(['AirAfter']*len(raw_chamber_data_air_after),dtype='|S25')
        raw_chamber_data_air_after = recfunctions.append_fields(raw_chamber_data_air_after,
                                                                 'status',
                                                                 status_array,
                                                                 dtypes='|S25',
                                                                 usemask=False)
        raw_chamber_data = recfunctions.stack_arrays((raw_chamber_data,raw_chamber_data_air_after),usemask=False)

        return raw_chamber_data

def rotate_struct(ev, ra, dec):
    r"""Wrapper around the rotate-method in skylab.utils for structured
    arrays.

    Parameters
    ----------
    ev : structured array
        Event information with ra, sinDec, plus true information

    ra, dec : float
        Coordinates to rotate the true direction onto

    Returns
    --------
    ev : structured array
        Array with rotated value, true information is deleted

    """
    names = ev.dtype.names

    rot = np.copy(ev)

    # Function call
    rot["ra"], rot_dec = rotate(ev["trueRa"], ev["trueDec"],
                                ra * np.ones(len(ev)), dec * np.ones(len(ev)),
                                ev["ra"], np.arcsin(ev["sinDec"]))

    if "dec" in names:
        rot["dec"] = rot_dec
    rot["sinDec"] = np.sin(rot_dec)

    # "delete" Monte Carlo information from sampled events
    mc = ["trueRa", "trueDec", "trueE", "ow"]

    return drop_fields(rot, mc)

def DoJoin(balrog, row, size, odir, zz, names, end=None, cols=False, field=None):
    if not os.path.exists(odir):
        os.makedirs(odir)

    if end is None:
        end = row + len(zz)

    if end > size:
        end = size
    ee = end - row

    b = balrog[-1].read(rows=np.arange(row,end))
    d = []
    for name in names:
        d.append(zz[name][:ee])

    n = list(names)
    n.append('field')
    d.append(np.array([field]*len(b)))

    c = rec.append_fields(b, n, d)
    if 'table' in c.dtype.names:
        c = rec.drop_fields(c, 'table')

    ofile = os.path.join(odir, '%i-%i.fits'%(row,end))
    esutil.io.write(ofile, c, clobber=True)
    
    if cols:
        return end, c.dtype.names
    else:
        return end

def _prepare_experiments(experiments):
    '''
    transform the experiments structured array into a numpy array.

    Parameters
    ----------
    experiments : structured array
    
    Returns
    -------
    ndarray
    
    '''
    experiments = recfunctions.drop_fields(experiments, "scenario_id", 
                                           asrecarray=True)
    uncs = recfunctions.get_names(experiments.dtype)

    temp_experiments = np.zeros((experiments.shape[0], len(uncs)))
    
    for i, u in enumerate(uncs):
        try: 
            temp_experiments[:,i] = experiments[u].astype(np.float)
        except ValueError:
            
            data = experiments[u]
            entries = sorted(list(set(data)))
            
            for j, entry in enumerate(entries):
                temp_experiments[data==entry,i] = j
    
    return temp_experiments, uncs

def read_positions():
    head,points1 = csv_parse.read('../Baltay-fibers_random.csv',delimiter=' ')
    head,points2 = csv_parse.read('../Baltay-fibers_residual.csv',delimiter=' ')
    points2 = rec.drop_fields(points2,('r','theta'))
    points2['Number'] += 10000 # to distinguish them from the "randoms"
    points = np.hstack((points1,points2))
    return points

    def remove_columns(self, col_names=None):
        '''
        This function will remove the all the columns within with names in
        col_names from all the datasets in self.columnar_data.

        Parameters
        ----------
        col_names : string or list
            The name or names of columns to be removed

        '''
        
        if col_names != None:

            if type(col_names) == str:
                col_names = [col_names]
            else:
                col_names = list(col_names)

            # Format column names
            col_names = ff.format_headers(col_names)

            removed_data = []
            for data in self.columnar_data:
                removed_data.append(drop_fields(data, col_names))
            self.columnar_data = removed_data

    def __init__(self, x, y, mass_min=0.05, mode=sdutil.BINARY):
        ''' init

        '''
        x = recfunctions.drop_fields(x, "scenario_id", asrecarray=True)

        self.x = x
        self.y = y
        self.mass_min = mass_min
        self.mode = mode

        # we need to transform the structured array to a ndarray
        # we use dummy variables for each category in case of categorical
        # variables. Integers are treated as floats
        self.feature_names = []
        columns = []
        for unc, dtype in x.dtype.descr:
            dtype = x.dtype.fields[unc][0]
            if dtype == np.object:
                categories = sorted(list(set(x[unc])))
                for cat in categories:
                    label = '{}{}{}'.format(unc, self.sep, cat)
                    self.feature_names.append(label)
                    columns.append(x[unc] == cat)
            else:
                self.feature_names.append(unc)
                columns.append(x[unc])

        self._x = np.column_stack(columns)
        self._boxes = None
        self._stats = None

def data_save(data, output_filename):
    # This isn't too hard, except we're going to put a copy of the
    # measures we actually care about at the beginning!
    names = list(data.dtype.names)
    
    # Find all the columns that have 'av' in their title and not
    # and not '_mask'
    drop_names = [ name for name in names if (name.find('_av_') == -1) | (name.find('_mask') > 0) ]
    drop_names.pop(0)

    important_data = rec.drop_fields(data, drop_names, usemask=False, asrecarray=True)
    
    names = list(important_data.dtype.names)
    
    # Strip the beginning part to get shorter and easy to manage variable names
    names[1:] = [ name[6:] for name in names[1:] ]
    names[1:] = [ name[:(-8)] for name in names[1:] ]
    names[1:] = [ name[0].upper() + name[1:] + 'Cort' for name in names[1:] ]
    names[0] = 'SubID'
    important_data.dtype.names = names

    # Create two temporaray output_filenames:
    temp_filename1 = output_filename + '_temp1'
    temp_filename2 = output_filename + '_temp2'
    
    plt.rec2csv(data, temp_filename1, delimiter='\t', formatd=None, withheader=True)
    plt.rec2csv(important_data, temp_filename2, delimiter='\t', formatd=None, withheader=True)
    
    mcf.KW_paste(temp_filename2, temp_filename1, output_filename)
    mcf.KW_rmforce(temp_filename1)
    mcf.KW_rmforce(temp_filename2)

def merge_cort(data, cortisol_filename):
    
    cort_data = np.genfromtxt(cortisol_filename, dtype=None, names=True, delimiter='\t')
    
    names = list(cort_data.dtype.names)
    
    # Find all the columns in cort_data that have 'av' in their title
    # and not '_mask'
    drop_names = names[8:]

    cort_data = nprf.drop_fields(cort_data, drop_names, usemask=False, asrecarray=True)
    
    data = nprf.join_by('SubID', data, cort_data, jointype='leftouter',
                            r1postfix='KW', r2postfix='KW2', usemask=False,asrecarray=True)
    
    # Bizzarely, the join_by function pads with the biggest numbers it can think of!
    # So we're going to replace everything over 999 with 999
    for name in names[1:8]:
        data[name][data[name]>999] = 999
    
    # Define a UsableCort field: 1 if ANY of the cortisol values are not 999
    cort_array = np.vstack( [ data[name] for name in names[1:8]])
    usable_cort_array = np.ones(cort_array.shape[1])
    usable_cort_array[np.any(cort_array<>999, axis=0)] = 1
    
    data = nprf.append_fields(base = data, names='UsableCort', data = usable_cort_array, usemask=False)

    return data

    def _rotate_subset(self, value, orig_experiments, logical): 
        '''
        rotate a subset
        
        Parameters
        ----------
        value : list of strings
        orig_experiment : numpy structured array
        logical : boolean array
        
        '''
        list_dtypes = [(name, "<f8") for name in value]
        
        #cast everything to float
        drop_names = set(rf.get_names(orig_experiments.dtype)) - set(value)
        orig_subset = rf.drop_fields(orig_experiments, drop_names, 
                                               asrecarray=True)
        subset_experiments = orig_subset.astype(list_dtypes).view('<f8').reshape(orig_experiments.shape[0], len(value))
 
        #normalize the data
        mean = np.mean(subset_experiments,axis=0)
        std = np.std(subset_experiments, axis=0)
        std[std==0] = 1 #in order to avoid a devision by zero
        subset_experiments = (subset_experiments - mean)/std
        
        #get the experiments of interest
        experiments_of_interest = subset_experiments[logical]
        
        #determine the rotation
        rotation_matrix =  self._determine_rotation(experiments_of_interest)
        
        #apply the rotation
        subset_experiments = np.dot(subset_experiments,rotation_matrix)
        return rotation_matrix, subset_experiments

def setup_cart(results, classify, incl_unc=[], mass_min=0.05):
    """helper function for performing cart
    
    Parameters
    ----------
    results : tuple of structured array and dict with numpy arrays
              the return from :meth:`perform_experiments`.
    classify : string, function or callable
               either a string denoting the outcome of interest to 
               use or a function. 
    incl_unc : list of strings
    mass_min : float
    
    
    Raises
    ------
    TypeError 
        if classify is not a string or a callable.
    
    """
    
    if not incl_unc:
        x = np.ma.array(results[0])
    else:
        drop_names = set(recfunctions.get_names(results[0].dtype))-set(incl_unc)
        x = recfunctions.drop_fields(results[0], drop_names, asrecarray = True)
    if type(classify)==types.StringType:
        y = results[1][classify]
    elif callable(classify):
        y = classify(results[1])
    else:
        raise TypeError("unknown type for classify")
    
    return CART(x, y, mass_min)

 def __init__(self, filename, date_sep='-', time_sep=':', format='stroke_DC3'):
     """ Load NLDN data from a file, into a numpy named array stored in the
         *data* attribute. *data*['time'] is relative to the *basedate* datetime
         attribute
         """
     self.format=format
     
     dtype_specs = getattr(self, format)
     
     
     nldn_initial = np.genfromtxt(filename, dtype=dtype_specs['columns'])
     date_part = np.genfromtxt(nldn_initial['date'],
                     delimiter=date_sep, dtype=dtype_specs['date_dtype'])
     time_part = np.genfromtxt(nldn_initial['time'],
                     delimiter=time_sep, dtype=dtype_specs['time_dtype'])
     dates = [datetime(a['year'], a['month'], a['day'], b['hour'], b['minute']) 
                 for a, b in zip(date_part, time_part)]
     min_date = min(dates)
     min_date = datetime(min_date.year, min_date.month, min_date.day)
     t = np.fromiter( ((d-min_date).total_seconds() for d in dates), dtype='float64')
     t += time_part['second']
     
     self.basedate = min_date
     data = drop_fields(nldn_initial, ('date', 'time'))
     data = append_fields(data, 'time', t)
     
     self.data = data
     

def read_originals():
    """Return the originally defined fiber positions, sorted by x and y."""
    head, points1 = csv_parse.read("../Baltay-fibers_random.csv", delimiter=" ")
    head, points2 = csv_parse.read("../Baltay-fibers_residual.csv", delimiter=" ")
    points2 = recfunc.drop_fields(points2, ("r", "theta"))
    points2["Number"] += 10000  # to distinguish them from the "randoms"
    return np.hstack((points1, points2))

  def __update(s):

    # Remove inactive channels
    names = s.records.dtype.names
    s.records = rcf.drop_fields(s.records, drop_names=s.inactive)
    s.chans = [s.chans[i] for i in xrange(len(names)) if names[i] not in s.inactive]

    s.__refresh_active()

def convert(ifile):
    folder = "/lustre/scratch/astro/cs390/LGalaxies_Hen15_PublicRelease/MergerTrees/MR/treedata/"
    lastsnap = 63
    alistfile = "/lustre/scratch/astro/cs390/LGalaxies_Hen15_PublicRelease/input/zlists/zlist_MR.txt"
    f = h5py.File(folder+'/trees_'+str(ifile)+".hdf5", 'w')
    # Version
    f.attrs.create('Version', 0, dtype=numpy.int32)
    # Subversion
    f.attrs.create('Subversion', 1, dtype=numpy.int32)
    # Title
    f.attrs.create('Title', "The Mighty Peter")
    # Description
    f.attrs.create('Description', "This is for testing")
    # BoxsizeMpc -- I'm not convinced that we should use Mpc instead Mpc/h (It's quite difficult to remember)
    # so I will use Mpc/h to avoid the errors from myself
    f.attrs.create('BoxsizeMpc_h', 62.5, dtype=numpy.float32)
    # OmegaBaryon
    f.attrs.create('OmegaBaryon', 0.044, dtype=numpy.float32)
    # OmegaCDM
    f.attrs.create('OmegaCDM', 0.27-0.044, dtype=numpy.float32)
    # H100
    f.attrs.create('H100', 0.704, dtype=numpy.float32)
    # Sigma8
    f.attrs.create('Sigma8', 0.807, dtype=numpy.float32)
    
    #Group -- Snapshot
    snapshot_grp = f.create_group("Snapshots")
    (nsnaps,snapshot_data) = load_snapshot(alistfile)
    #NSnap
    print numpy.int32(nsnaps)
    snapshot_grp.attrs['NSnap'] = numpy.int32(nsnaps)
    #Snap
    snapshot_snap = snapshot_grp.create_dataset('Snap', data=snapshot_data)

    #Group -- MergerTrees
    mergertree_grp = f.create_group("MergerTrees")
    verbose = 1
    print "Reading tree",ifile
    (nTrees,nHalos,nTreeHalos,output_Halos,output_HaloIDs) = read_lgal_input_fulltrees_withids(folder,lastsnap,ifile,verbose)
    print "Done reading tree",ifile
    #TableFlag
    mergertree_grp.attrs['TableFlag'] = numpy.int32(1)
    #NTree
    mergertree_grp.attrs['NTrees'] = numpy.int32(nTrees)
    #NHalo
    mergertree_grp.attrs['NHalos'] = numpy.int32(nHalos)
    #NHalosInTree

    nhalosintree_data = mergertree_grp.create_dataset('NHalosInTree', data=nTreeHalos.astype(numpy.int32))
    #Halo
    print "Merging arrays"
    #halo = rfn.merge_arrays((output_Halos,output_HaloIDs), flatten = True, usemask = False)
    halo = join_struct_arrays((output_Halos,output_HaloIDs))
    print "Done merging arrays"
    halo = rfn.drop_fields(halo,['dummy','PeanoKey'])
    print "Outputting merger trees"
    nhalosintree_data = mergertree_grp.create_dataset('Halo', data=halo)
    print "Done"

    def drop_extra_columns(self):
        """Remove any optional columns from this CopyNumArray.

        Returns a new copy with only the core columns retained:
            log2 value, chromosome, start, end, bin name.
        """
        result = self.__class__(self.sample_id)
        result.data = rfn.drop_fields(self.data, self._xtra)
        return result

def plot_cdfs(x, y, ccdf=False):
    '''plot cumulative density functions for each column in x, based on the 
    classification specified in y.
    
    Parameters
    ----------
    x : recarray
        the experiments to use in the cdfs
    y : ndaray 
        the categorization for the data
    ccdf : bool, optional
           if true, plot a complementary cdf instead of a normal cdf.
    
    '''
    x = rf.drop_fields(x, "scenario_id", asrecarray=True)
    uncs = rf.get_names(x.dtype)
    cp = sns.color_palette()
    
    n_col = 4
    n_row = len(uncs)//n_col +1
    size = 3 
    aspect = 1
    figsize = n_col * size * aspect, n_row * size
    fig, axes = plt.subplots(n_row, n_col,
                             figsize=figsize,
                             squeeze=False)

    for i, unc in enumerate(uncs):
        discrete = False
        
        i_col = i % n_col
        i_row = i // n_col
        ax = axes[i_row, i_col]
        
        data = x[unc]
        if x.dtype[unc] == np.dtype('O'):
            discrete = True
        plot_cdf(ax, unc, data, y, discrete, ccdf=ccdf)
    
    # last row might contain empty axis, 
    # let's make them disappear
    i_row = len(uncs) // n_col
    i_col = len(uncs) % n_col
    for i_col in range(i_col, n_col):
        ax = axes[i_row, i_col]
        ax.set_xticklabels([])
        ax.set_xticks([])
        ax.set_yticklabels([])
        ax.set_yticks([])
        
        sns.despine(ax=ax, top=True, right=True, left=True, bottom=True)
    
    proxies, labels = build_legend(x, y)
    
    fig.legend(proxies, labels, "upper center")

    return fig