def get_demodulated_data_from_list(filelist,freq=10,supply_index=True,phase_offset=0):
    filelist.sort() #just in case
    dd=[]
    for f in filelist:
        #only use full size files
        stats=os.stat(f)
        if stats.st_size == 10752000:
            print f
            d=demod.demodulate_dat(f,freq,supply_index=True,phase_offset=phase_offset)
            #filename is start of data taking (I think) and we'll just add 1/samprate seconds per rev
            h=np.float64(f[-12:-10])
            m=np.float64(f[-10:-8])
            s=np.float64(f[-8:-6])
            t=h+m/60.+(s+(d['rev']-d['rev'][0])/samprate)/3600.
            d=recf.append_fields(d,'localtime',t)
            ut=np.mod(t+7.,24.)
            if len(f)>21:
                y=np.zeros(len(d),dtype=np.int)+np.int(f[-21:-17])
                mo=np.zeros(len(d),dtype=np.int)+np.int(f[-17:-15])
                dy=np.zeros(len(d),dtype=np.int)+np.int(f[-15:-13])
                ut=np.mod(t+7.,24.)
                utt=t+7.
                dy[utt>ut]=dy[utt>ut]+1
                d=recf.append_fields(d,['year','month','day'],[y,mo,dy])
            d=recf.append_fields(d,'ut',ut)
            dd.append(d)
    return np.concatenate(dd)

def get_cat(cat_filename, cut=0):
  '''  
  Opens the target file using numpy genfromtxt.  Assumes the format is given by thehead.  If coordinates are sexagesimal then converts them to degrees and appends a two new columns ('degra' and 'degdec') to the record array.  Even if coordinates are already in degrees, the 'degra' and 'degdec' columns are appended but will be identical to the original 'ra' and 'dec'.
  '''

  f = open(cat_filename,'r')
  header = f.readline().lower().split()
  f.close()

  thehead = ['hstid', 'field', 'ra', 'dec', 'v', 'verr', 'bvcol', 'bvcolerr', 'vicol', 'vicolerr']
  dtype = ['S30', 'S30', 'S30', 'S30', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8']
  
  if header == thehead:
    data = np.genfromtxt(cat_filename, names=header, dtype=dtype, skip_header=True, skip_footer=cut)
  
    if ':' in data['ra'][0]:
      ra_degrees= np.array([sex2deg(x) for x in data['ra']])
      dec_degrees = np.array([sex2deg(x, RA=False) for x in data['dec']])      
      
      cat = rfn.append_fields(data, names=['degra','degdec'], data=[ra_degrees, dec_degrees], 
                              dtypes=['f8','f8'], usemask=False)
    
    else:
      cat = rfn.append_fields(data, names=['degra','degdec'], data=[data['ra'], data['dec']], 
                              dtypes=['f8','f8'], usemask=False)
    
    return cat
  
  else:
    sys.exit('Columns need to be title: %s' % thehead)

 def readH5(self, fname):
     """ 
     Reads Catalog from H5 file, specified as argument
     """
     of=h5py.File(fname, "r")
     self.data=of["objects"].value
     self.meta=of["meta"].attrs
     if "dNdz" in of.keys():
         self.dNdz=of['dNdz'].value 
     if "bz" in of.keys():
         self.bz=of['bz'].value
     self.window=window.readWindowH5(of['window'])
     self.photoz=photoz.readPhotoZH5(of['photoz'])
     cversion=float(self.meta['version'])
     if cversion==0.1:
         print("updating 0.1 to version ", self.version)
         self.data=recfunctions.append_fields(self.data,'sigma_pz',(1+self.data["z_real_t"])*self.data["z_error"],
                                              usemask=False)
         self.data=recfunctions.append_fields(self.data,'z',self.data["z_real_t"]+(1+self.data["z_real_t"])*self.data["z_error"],
                                              usemask=False)
         self.data=self.data[ [ name for name in self.data.dtype.names if name not in ["z_real_t", "z_rsd_t","z_error"] ] ]
     if cversion==0.2:
         print("WARNING: upgrading from 0.2 to 0.3, photozs internally slightly inconsistent.")
         self.data=recfunctions.append_fields(self.data,'sigma_pz',(1+self.data["z"])*self.photoz.sigma,
                                              usemask=False)

    def __init__(self,d={}):
# time: timestep number
        self.time=d.get('time',0)
# box: numpy array of lattice vectors
        box = d.get('box',N.zeros((3,3)))
        self.box = box
        if box.shape == (3,3): self.vc = box
        elif box.shape == (3): self.vc = N.diag(box)
        else:                  raise ValueError ('Box should be (3,3) or (3) array')
# atoms: atoms numpy array
        atoms = d.get('atoms',[])
        
        leg_list = atoms.dtype.names

        if not 'id' in leg_list:       # add 'id' column
            atoms = nlrf.append_fields(atoms, 'id', N.arange(len(atoms))+1, asrecarray=True, usemask=False)
   
        if not 'itype' in leg_list:    # add 'itype' column
            if 'label' in leg_list:
                labels = list(N.unique(atoms['label']))
                labels = dict(zip(labels, range(1,len(labels)+1)))
                ityp = N.array([labels[atom['label']] for atom in atoms])
                atoms = nlrf.append_fields(atoms, 'itype', ityp, asrecarray=True, usemask=False)
            else:   
                atoms = nlrf.append_fields(atoms, 'itype', N.ones(len(atoms)), asrecarray=True, usemask=False)

        self.atoms = atoms

    def analyze_chamber_data(self,raw_chamber_data):
        ethanol_data = raw_chamber_data[raw_chamber_data['status']=='Ethanol']
        analyzed_ethanol_data = self.analyze_data(ethanol_data)
        status_array = numpy.array(['Ethanol']*len(analyzed_ethanol_data),dtype='|S25')
        analyzed_chamber_data = recfunctions.append_fields(analyzed_ethanol_data,
                                                           'status',
                                                           status_array,
                                                           dtypes='|S25',
                                                           usemask=False)

        air_before_data = raw_chamber_data[raw_chamber_data['status']=='AirBefore']
        if air_before_data.size != 0:
            analyzed_air_before_data = self.analyze_data(air_before_data)
            status_array = numpy.array(['AirBefore']*len(analyzed_air_before_data),dtype='|S25')
            analyzed_air_before_data = recfunctions.append_fields(analyzed_air_before_data,
                                                                  'status',
                                                                  status_array,
                                                                  dtypes='|S25',
                                                                  usemask=False)
            analyzed_chamber_data = recfunctions.stack_arrays((analyzed_air_before_data,analyzed_chamber_data),usemask=False)


        air_after_data = raw_chamber_data[raw_chamber_data['status']=='AirAfter']
        if air_after_data.size != 0:
            analyzed_air_after_data = self.analyze_data(air_after_data)
            status_array = numpy.array(['AirAfter']*len(analyzed_air_after_data),dtype='|S25')
            analyzed_air_after_data = recfunctions.append_fields(analyzed_air_after_data,
                                                                  'status',
                                                                  status_array,
                                                                  dtypes='|S25',
                                                                  usemask=False)
            analyzed_chamber_data = recfunctions.stack_arrays((analyzed_chamber_data,analyzed_air_after_data),usemask=False)

        return analyzed_chamber_data

def whichGalaxyProfile(sdss):

    exp_L = np.exp(np.array([sdss['LNLEXP_G'],sdss['LNLEXP_R'],sdss['LNLEXP_I'],sdss['LNLEXP_Z']])).T
    dev_L = np.exp(np.array([sdss['LNLDEV_G'],sdss['LNLDEV_R'],sdss['LNLDEV_I'],sdss['LNLDEV_Z']])).T
    star_L = np.exp(np.array([sdss['LNLSTAR_G'],sdss['LNLSTAR_R'],sdss['LNLSTAR_I'],sdss['LNLSTAR_Z']])).T

    expfracL = exp_L /(exp_L + dev_L + star_L)
    devfracL = dev_L /(exp_L + dev_L + star_L)
    
    modelmode = np.zeros((len(sdss), 4), dtype=np.int32)

    expmodel = (expfracL > 0.5)
    modelmode[expmodel] = 0
    devmodel = (devfracL > 0.5)
    modelmode[devmodel] = 1
    neither = - (expmodel | devmodel)
    modelmode[neither] = 2
    
    sdss = rf.append_fields(sdss, 'BESTPROF_G', modelmode[:,0])
    sdss = rf.append_fields(sdss, 'BESTPROF_R', modelmode[:,1])
    sdss = rf.append_fields(sdss, 'BESTPROF_I', modelmode[:,2])
    sdss = rf.append_fields(sdss, 'BESTPROF_Z', modelmode[:,3])
    
    #print ' exp :', np.sum(expmodel),' dev :', np.sum(devmodel), 'neither :', np.sum(neither)
    return sdss

 def training_split(self, holdout_unit='none', holdout_prop=.2):
     ''' Splits the data up into test and train subsets '''
     if holdout_prop > .99 or holdout_prop < .01:
         raise ValueError('The holdout proportion must be between .1 and .99.')
     if holdout_unit == 'none':
         self.training_data = self.observation_matrix
         self.test_data = self.prediction_matrix
         self.training_type = 'make predictions'
         print 'Fitting model to all data'
     elif holdout_unit == 'datapoint':
         holdouts = np.random.binomial(1, holdout_prop, self.data_rows)
         self.training_data = np.delete(self.observation_matrix, np.where(holdouts==1)[0], axis=0)
         self.test_data = np.delete(self.observation_matrix, np.where(holdouts==0)[0], axis=0)
         self.training_type = 'datapoint'
         print 'Fitting model to ' + str((1-holdout_prop)*100) + '% of datapoints'
     elif holdout_unit == 'country-year':
         country_years = [self.observation_matrix.country[i] + '_' + str(self.observation_matrix.year[i]) for i in range(self.data_rows)]
         data_flagged = recfunctions.append_fields(self.observation_matrix, 'holdout', np.zeros(self.data_rows)).view(np.recarray)
         for i in np.unique(country_years):
             data_flagged.holdout[np.where(data_flagged.country + '_' + data_flagged.year.astype('|S4')==i)[0]] = np.random.binomial(1, holdout_prop)
         self.training_data = np.delete(data_flagged, np.where(data_flagged.holdout==1)[0], axis=0)
         self.test_data = np.delete(data_flagged, np.where(data_flagged.holdout==0)[0], axis=0)
         self.training_type = 'country-year'
         print 'Fitting model to ' + str((1-holdout_prop)*100) + '% of country-years'
     elif holdout_unit == 'country':
         data_flagged = recfunctions.append_fields(self.observation_matrix, 'holdout', np.zeros(self.data_rows)).view(np.recarray)
         for i in self.country_list:
             data_flagged.holdout[np.where(data_flagged.country==i)[0]] = np.random.binomial(1, holdout_prop)
         self.training_data = np.delete(data_flagged, np.where(data_flagged.holdout==1)[0], axis=0)
         self.test_data = np.delete(data_flagged, np.where(data_flagged.holdout==0)[0], axis=0)
         self.training_type = 'country'
         print 'Fitting model to ' + str((1-holdout_prop)*100) + '% of countries'
     else:
         raise ValueError("The holdout unit must be either 'datapoint', 'country-year', or 'country'.")

    def getClassifiers(self):
        if not os.path.exists(self.outDir):
            os.mkdir(self.outDir)
        outDir = self.outDir + os.sep + "classPickle"
        if not os.path.exists(outDir):
            os.mkdir(outDir)
        class1Save = outDir + os.sep + "classifier1.pkl"
        class2Save = outDir + os.sep + "classifier2.pkl"
        
        class1Exists = os.path.exists(class1Save)
        class2Exists = os.path.exists(class2Save)

        if not (class1Exists and class2Exists):
            self._setupTempDir()
            self.fitsFiles = [f[:-5] for f in os.listdir(self.fitsFolder) if ".fits" in f]
            self.fitsFilesLoc = [os.path.abspath(self.fitsFolder + os.sep + f) for f in os.listdir(self.fitsFolder) if ".fits" in f]
            
            for f in self.fitsFiles:
                self.mainCatalog[f] = self.getCatalog(self.fitsFolder + os.sep + f + ".fits", ishape=True)
                self.candidateMask[f] = self._getCandidateMask(self.mainCatalog[f], np.loadtxt(self.fitsFolder + os.sep + f + ".txt"))
                self.mainCatalog[f] = append_fields(self.mainCatalog[f], 'WEIGHT', self.candidateMask[f] * 1.0, usemask=False)    
                self.mainCatalog[f] = append_fields(self.mainCatalog[f], 'EXTENDED', self.candidateMask[f], usemask=False)    
                self.mainCatalog[f] = append_fields(self.mainCatalog[f], 'HLR', np.zeros(self.mainCatalog[f].shape), usemask=False)    
                self.mainCatalog[f] = append_fields(self.mainCatalog[f], 'MAG', np.zeros(self.mainCatalog[f].shape), usemask=False)
            self._trainClassifier()
            joblib.dump(self.sc, class1Save) 
            joblib.dump(self.sc2, class2Save) 
        else:
            self.sc = joblib.load(class1Save)
            self.sc2 = joblib.load(class2Save)
            

        #self._testClassifier(catalog, candidateMask)
        #self._cleanTempDir()
        self._debug("Classifier generated. Now you can invoke .clasify(catalog)")

def append_fields(base, names, data, dtypes=None, fill_value=-1,
                  usemask=False,   # Different from recfunctions default
                  asrecarray=False):
    """Append fields to numpy structured array
    If fields already exists in data, will overwrite
    """
    if isinstance(names, (tuple, list)):
        # Add multiple fields at once
        if dtypes is None:
            dtypes = [d.dtype for d in data]
        # Convert to numpy arrays so we can use boolean index arrays
        names = np.array(names)
        data = np.array(data)
        dtypes = np.array(dtypes)
        not_yet_in_data = True ^ np.in1d(names, base.dtype.names)
        # Append the fields that were not in the data
        base = recfunctions.append_fields(base,
                                          names[not_yet_in_data].tolist(),
                                          data[not_yet_in_data].tolist(),
                                          dtypes[not_yet_in_data].tolist(),
                                          fill_value, usemask, asrecarray)
        # Overwrite the fields that are already in the data
        for i in np.where(True ^ not_yet_in_data)[0]:
            base[names[i]] = data[i]
        return base
    else:
        # Add single field
        if names in base.dtype.names:
            # Field already exists: overwrite data
            base[names] = data
            return base
        else:
            return recfunctions.append_fields(base, names, data, dtypes,
                                              fill_value, usemask, asrecarray)

 def load_training_data(self, training_data='', **extras):
     """Read an HDF5 file with `parameters` a structured ndarray and
     `spectra` an ndarray.  Convert to a structured array of labels of
     length `ntrain` with `nlabel` fields. and an ndarray of training
     spectra of shape (nwave, ntrain).
     """
     self.has_errors = False
     with h5py.File(training_data, "r") as f:
         self.library_spectra = f['spectra'][:]
         self.library_labels = f['parameters'][:]
         self.wavelengths = f['wavelengths'][:]
         try:
             self.library_snr = self.library_spectra / f['uncertainty'][:]
             self.has_errors = True
         except:
             pass
         ancillary = f['ancillary'][:]
     # add and rename labels here.  Note that not all labels need to be or
     # will be used in the feature generation
     newfield = ['logt', 'miles_id']
     newdata = [np.log10(self.library_labels['teff']), ancillary['miles_id']]
     self.library_labels = rfn.append_fields(self.library_labels,
                                              newfield, newdata, usemask=False)
     try:
         # assuming f_nu
         fbol = np.trapz(self.library_spectra / self.wavelengths**2, self.wavelengths)
         newfield = ['logl', 'luminosity', 'fbol']
         newdata = [ancillary['logl'], 10**ancillary['logl'], fbol]
         self.library_labels = rfn.append_fields(self.library_labels,
                                                  newfield, newdata, usemask=False)            
     except:
         pass
     self.reset_mask()

def merge_paradigm(trial_info, paradigm, behavioural=None, **conf):
    
    
    baseline_condition = ''
    for arg in conf:
        if arg == 'baseline':
            baseline_condition = conf[arg]

    mask_blink_outlier = np.in1d(paradigm['Trial'], trial_info['Trial'])
        
    trial_info = nprec.append_fields(trial_info, 
                                     'Label', 
                                     paradigm['Label'][mask_blink_outlier]).data

    
    mask_task = paradigm['Label'] != baseline_condition
        
    
    print 'Trials no.' + str(len(trial_info))

    if behavioural != None:
        m = mask_task * mask_blink_outlier
        m = m[1::2]
    
        trial_task_info = trial_info[trial_info['Label'] != baseline_condition]
    
        trial_cond = nprec.append_fields(trial_task_info,
                                     behavioural.dtype.names, 
                                     [behavioural[b][m] for b in behavioural.dtype.names]).data
        
        return trial_cond, trial_info
    
    else:
        
        return trial_info

def Tracks2Cells(tracks, falarms=None) :
    """
    Convert lists of tracks (and falarms) into a single recarray of storm cells
    with track IDs.

    This can be reversed with Cells2Tracks().
    """
    if falarms is None :
        falarms = []

    # NOTE: This function can not handle arrays of tracks that do/do not have
    #       a trackID field in a mix.  Either they all have it, or not.
    if not any('trackID' in aTrack.dtype.names for aTrack in tracks) :
        tracks = [nprf.append_fields(aTrack, 'trackID',
                                     [trackIndex] * len(aTrack),
                                     usemask=False)
                  for trackIndex, aTrack in enumerate(tracks)]

    if not any('trackID' in aTrack.dtype.names for aTrack in falarms) :
        falarms = [nprf.append_fields(aTrack, 'trackID',
                                      [-trackIndex - 1] * len(aTrack),
                                      usemask=False)
                   for trackIndex, aTrack in enumerate(falarms)]

    # If both are empty, then create an array without hstack()
    if len(tracks) != 0 or len(falarms) != 0 :
        allCells = np.hstack(tracks + falarms)
    else :
        allCells = np.array([], dtype=volume_dtype)

    return allCells

 def getRAandDec(self, fitsFile, catalog):
     self._debug("\tRunning xy2sky")
     tempDir = self.tempDir + os.sep + "xy2sky"
     self._debug("\tGenerating temp directory at %s" % tempDir)
     
     try:
         shutil.rmtree(tempDir)  # delete directory
     except OSError as exc:
         if exc.errno != errno.ENOENT:
             raise  # re-raise exception
     os.mkdir(tempDir)
     
     imfile = tempDir + os.sep + "imfile.txt"
     np.savetxt(imfile, catalog[['X_IMAGE','Y_IMAGE']], fmt="%0.3f")
     outfile = tempDir + os.sep + "skys.txt"
     # | cut -d \" \" -f 1-2 > %s
     commandline = wcsToolsPath + "/xy2sky -d %s @%s | awk '{print $1,$2}'> %s" % (fitsFile, imfile, outfile)
     p = subprocess.Popen(["/bin/bash", "-i", "-c", commandline], stderr=subprocess.PIPE, stdout=subprocess.PIPE)
     output = p.communicate() #now wait
     
     res = np.loadtxt(outfile)
     catalog = append_fields(catalog, 'RA', res[:,0], usemask=False)    
     catalog = append_fields(catalog, 'DEC', res[:,1], usemask=False)    
     
     return catalog

    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 __init__(self, pathToFile=""):
        self.path = pathToFile
        self._points = np.genfromtxt(pathToFile, delimiter=' ', names='birth, death', dtype='f8, f8')
        self.lifespan_coords = []
        self.birth_coords = []

        # Generate commonly-used derived fields: lifespan, avg_coord
        self._points = rf.append_fields(self._points, 'lifespan', self._points['death'] - self._points['birth'], dtypes='f8')
        self._points = rf.append_fields(self._points, 'avg_coord', (self._points['death'] + self._points['birth'])/2, dtypes='f8')

def add_constant(data, prepend=False):
    '''
    This appends a column of ones to an array if prepend==False.

    For ndarrays it checks to make sure a constant is not already included.
    If there is at least one column of ones then the original array is
    returned.  Does not check for a constant if a structured or recarray is
    given.

    Parameters
    ----------
    data : array-like
        `data` is the column-ordered design matrix
    prepend : bool
        True and the constant is prepended rather than appended.

    Returns
    -------
    data : array
        The original array with a constant (column of ones) as the first or
        last column.

    Notes
    -----

    .. WARNING::
       The default of prepend will be changed to True in the next release of
       statsmodels. We recommend to use an explicit prepend in any permanent
       code.

    '''
    import warnings
    warnings.warn("The default of `prepend` will be changed to True in the "
                  "next release, use explicit prepend", FutureWarning)
    if not data.dtype.names:
        data = np.asarray(data)
        if np.any(data[0]==1):
            ind = np.squeeze(np.where(data[0]==1))
            if ind.size == 1 and np.all(data[:,ind] == 1):
                return data
            elif ind.size > 1:
                for col in ind:
                    if np.all(data[:,col] == 1):
                        return data
        data = np.column_stack((data, np.ones((data.shape[0], 1))))
        if prepend:
            return np.roll(data, 1, 1)
    else:
        return_rec = data.__class__ is np.recarray
        if prepend:
            ones = np.ones((data.shape[0], 1), dtype=[('const', float)])
            data = nprf.append_fields(ones, data.dtype.names, [data[i] for
                i in data.dtype.names], usemask=False, asrecarray=return_rec)
        else:
            data = nprf.append_fields(data, 'const', np.ones(data.shape[0]),
                    usemask=False, asrecarray = return_rec)
    return data

def DESdperp_to_SDSSdperp(fullsdss, fulldes):
    
    modelmag_g = fullsdss['MODELMAG_G'] - fullsdss['EXTINCTION_G']
    modelmag_r = fullsdss['MODELMAG_R'] - fullsdss['EXTINCTION_R']
    modelmag_i = fullsdss['MODELMAG_I'] - fullsdss['EXTINCTION_I']
    modelmag_z = fullsdss['MODELMAG_Z'] - fullsdss['EXTINCTION_Z']
    dperp_sdss = (modelmag_r - modelmag_i) - (modelmag_g - modelmag_r)/8.0

    des, sdss = match(fulldes, fullsdss)
    
    modelmag_g_des = des['MODELMAG_G_DES'] - des['XCORR_SFD98_G']
    modelmag_r_des = des['MODELMAG_R_DES'] - des['XCORR_SFD98_R']
    modelmag_i_des = des['MODELMAG_I_DES'] - des['XCORR_SFD98_I']
    modelmag_z_des = des['MODELMAG_Z_DES'] - des['XCORR_SFD98_Z']
    dperp_des = (modelmag_r_des - modelmag_i_des) - (modelmag_g_des - modelmag_r_des)/8.0

    expcut = (des['IM3_GALPROF'] == 1)
    devcut = (des['IM3_GALPROF'] == 2)
    
    #magcut = ((des['MODELMAG_G_DES'] < 22.0) & (des['MODELMAG_R_DES'] < 22.0) &
    #(des['MODELMAG_I_DES'] < 22.0) & (des['MODELMAG_Z_DES'] < 22.0) )
    #magcut = ((sdss['MODELMAG_R'] < 22.0) & (sdss['MODELMAG_I'] < 22.0)
    #          &(sdss['MODELMAG_G'] < 22.0) & (sdss['MODELMAG_Z'] < 22.0))
    
    use =  (#(18.0 < des['CMODELMAG_I_DES']) &
            #(20.9 > des['CMODELMAG_I_DES']) &
            #(des['MODELMAG_R_DES'] - des['MODELMAG_I_DES'] < 2.) &
            (des['FIBER2MAG_I_DES'] < 21.5 )
            )
            
    #use = use  & magcut
            
    des_exp = des[expcut & use]
    des_dev = des[devcut & use]
    sdss_exp = sdss[expcut & use]
    sdss_dev = sdss[devcut & use]
    
    #des_exp_dperp = dperp_fitting(des_exp, sdss_exp)
    #des_dev_dperp = dperp_fitting(des_dev, sdss_dev)
    
    #SDSSlike_dperp = np.zeros(len(fulldes), dtype=np.float32)
    #SDSSlike_dperp[expcut & use] = des_exp_dperp
    #SDSSlike_dperp[devcut & use] = des_dev_dperp
    
    try :
        #fulldes = rf.append_fields(fulldes, 'DPERP_DES', SDSSlike_dperp)
        fulldes = rf.append_fields(fulldes, 'DPERP', dperp_des)
        fullsdss = rf.append_fields(fullsdss, 'DPERP', dperp_sdss)

    except ValueError:
        #fulldes['DPERP_DES'] = SDSSlike_dperp
        fulldes['DPERP'] = dperp_des
        fullsdss['DPERP'] = dperp_sdss

    return fullsdss, fulldes

def count_good_trials():
    count_file = open(path_b_2_8+'/count_trials_blink.txt', 'w')      
    count_file.write('Subj C_inc C_tot  NC_inc NC_tot 1_inc 1_tot 2_inc 2_tot 3_inc 3_tot 4_inc 4_tot\r\n')
    for file in file_list:
        d_data = load_data_eye(path_blink, file)
        trial_info = extract_trials_info(d_data)
        mask_blink_outlier = np.in1d(paradigm['Trial'], trial_info['Trial'])
        trial_info = nprec.append_fields(trial_info, 
                                         'Condition', 
                                         paradigm['Condition'][mask_blink_outlier]).data
        task_trial = trial_info[trial_info['Condition'] != 'FIX']
        name = file.split('.')[0]
        try:
            behavioural = open_behavioural(path_b, name+'.xlsx')
        except IOError, err:
            print err
            continue
        m = mask_blink_outlier[1::2]
        trial_cond = nprec.append_fields(task_trial,
                                     ['Accuracy', 'Combination'], 
                                     [behavioural['Accuracy'][m], 
                                     behavioural['Combination'][m]]).data
        par = nprec.append_fields(paradigm[1::2],
                                     'Accuracy', 
                                     behavioural['Accuracy']).data
        '''
        trial_cond = trial_cond[trial_cond['Accuracy'] == 1]
        par = par[par['Accuracy'] == 1]
        behavioural = behavioural[behavioural['Accuracy'] == 1]
        '''
        count_file.write(file+' ')
        count_file.write(str(np.count_nonzero(trial_cond['Condition'] == 'C')))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(par['Condition'] == 'C')))
        count_file.write(' ')     
        count_file.write(str(np.count_nonzero(trial_cond['Condition'] == 'NC')))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(par['Condition'] == 'NC'))) 
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(trial_cond['Combination'] == 1)))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(behavioural['Combination'] == 1)))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(trial_cond['Combination'] == 2)))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(behavioural['Combination'] == 2)))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(trial_cond['Combination'] == 3)))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(behavioural['Combination'] == 3)))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(trial_cond['Combination'] == 4)))
        count_file.write(' ')  
        count_file.write(str(np.count_nonzero(behavioural['Combination'] == 4)))
        count_file.write(' \r\n')

    def readascii(self, infile, form='short'):

        if form=='short':
            dnames=('name', 'targ_ra', 'targ_dec', 'equinox', 'mag', 'band', 'priority')
            dformat=('S30', 'f4', 'f4', 'i4', 'f4', 'S1', 'f4')
            object_arr=np.loadtxt(infile, dtype={'names': dnames, 'formats': dformat},
                         converters={1:ra_read, 2:dec_read})
            #determine the missing values
            mnames=[]
            mtypes=[]
            for i in range(len(self.dnames)):
                if self.dnames[i] not in dnames:
                   mnames.append(self.dnames[i])
                   mtypes.append(self.dformat[i])
            #set up the default values
            default_list=[np.zeros(len(object_arr))]*len(mnames)
            default_list[0]=default_list[0]+self.default_width
            default_list[1]=default_list[1]+0.5*self.default_length
            default_list[2]=default_list[2]+0.5*self.default_length
            object_arr=rfn.append_fields(object_arr, names=mnames, data=default_list, dtypes=mtypes,
                     fill_value=0, usemask=False)
        elif form=='long':
            dnames=('name', 'targ_ra', 'targ_dec', 'equinox', 'mag', 'band', 'priority', 'width', 'length', 'tilt')
            dformat=('S30', 'f4', 'f4', 'i4', 'f4', 'S1', 'f4', 'f4', 'f4', 'f4')
            object_arr=np.loadtxt(infile, dtype={'names': dnames, 'formats': dformat},
                         converters={1:ra_read, 2:dec_read})
            #determine the missing values
            mnames=[]
            mtypes=[]
            for i in range(len(self.dnames)):
                if self.dnames[i] not in dnames:
                   mnames.append(self.dnames[i])
                   mtypes.append(self.dformat[i])
            #set up the default values
            default_list=[np.zeros(len(object_arr))]*len(mnames)
            object_arr=rfn.append_fields(object_arr, names=mnames, data=default_list, dtypes=mtypes,
                     fill_value=0, usemask=False)
        else:
            message='This format is not supported'
            raise SlitError(message)

        #set objects that are preselected
        object_arr['inmask_flag'] = 1.0*(object_arr['priority'] >= 1.0)
        #set reference stars
        object_arr['refstar_flag'] = 1.0*(object_arr['priority'] == -1.0)
     
        #stack the data if it already exists
        if self.data is None:  
           self.data=object_arr
        else:
           self.data=self.add_arrays(self.data, object_arr)
        # total number of objects:
        self.nobjects=len(self.data)
        self.update_flags()

def add_constant(data, prepend=True, has_constant='skip'):
    '''
    This appends a column of ones to an array if prepend==False.

    Parameters
    ----------
    data : array-like
        `data` is the column-ordered design matrix
    prepend : bool
        True and the constant is prepended rather than appended.
    has_constant : str {'raise', 'add', 'skip'}
        Behavior if ``data'' already has a constant. The default will return
        data without adding another constant. If 'raise', will raise an
        error if a constant is present. Using 'add' will duplicate the
        constant, if one is present. Has no effect for structured or
        recarrays. There is no checking for a constant in this case.

    Returns
    -------
    data : array
        The original array with a constant (column of ones) as the first or
        last column.
    '''
    if _is_using_pandas(data, None):
        # work on a copy
        return _pandas_add_constant(data.copy(), prepend, has_constant)
    else:
        data = np.asarray(data)
    if not data.dtype.names:
        var0 = data.var(0) == 0
        if np.any(var0):
            if has_constant == 'raise':
                raise ValueError("data already contains a constant.")
            elif has_constant == 'skip':
                return data
            elif has_constant == 'add':
                pass
            else:
                raise ValueError("Option {0} not understood for "
                                 "has_constant.".format(has_constant))
        data = np.column_stack((data, np.ones((data.shape[0], 1))))
        if prepend:
            return np.roll(data, 1, 1)
    else:
        return_rec = data.__class__ is np.recarray
        if prepend:
            ones = np.ones((data.shape[0], 1), dtype=[('const', float)])
            data = nprf.append_fields(ones, data.dtype.names,
                                      [data[i] for i in dat