def loadlocal_mnist(images_path, labels_path):
    """ Read MNIST from ubyte files.

    Parameters
    ----------
    images_path : str
        path to the test or train MNIST ubyte file
    labels_path : str
        path to the test or train MNIST class labels file

    Returns
    --------
    images : [n_samples, n_pixels] numpy.array
        Pixel values of the images.
    labels : [n_samples] numpy array
        Target class labels

    """
    with open(labels_path, 'rb') as lbpath:
        magic, n = struct.unpack('>II',
                                 lbpath.read(8))
        labels = np.fromfile(lbpath,
                             dtype=np.uint8)
    with open(images_path, 'rb') as imgpath:
        magic, num, rows, cols = struct.unpack(">IIII",
                                               imgpath.read(16))
        images = np.fromfile(imgpath,
                             dtype=np.uint8).reshape(len(labels), 784)

    return images, labels

    def Read(self):

        #return numpy.ones((256, 819)).astype('float32'), numpy.ones(256).astype('int32')

        with open(self.featureFile,"rb") as f:

            dt = numpy.dtype([('numSamples',(numpy.int32,1)),('sampPeriod',(numpy.int32,1)),('sampSize',(numpy.int16,1)),('sampKind',(numpy.int16,1))])
            header =  numpy.fromfile(f,dt.newbyteorder('>' if self.byteOrder==ByteOrder.BigEndian else '<'),count=1)

            numSamples = header[0]['numSamples']
            sampPeriod = header[0]['sampPeriod']
            sampSize   = header[0]['sampSize']
            sampKind   = header[0]['sampKind']

            # print 'Num samples = {}'.format(numSamples)
            # print 'Sample period = {}'.format(sampPeriod)
            # print 'Sample size = {}'.format(sampSize)
            # print 'Sample kind = {}'.format(sampKind)
            dt = numpy.dtype([('sample',(numpy.float32,sampSize/4))])
            samples = numpy.fromfile(f,dt.newbyteorder('>' if self.byteOrder==ByteOrder.BigEndian else '<'),count=numSamples)

        self._markDone()

        if self.labelFile is None:
            labels = None
        else:
            labels = ReadLabel(self.labelFile)

        return samples[:]['sample'], labels

def getCoincidences(fTimes, fChans, gate, radius, heraldChan):
  bufRes = 156.25e-12
  gate = int(gate/bufRes)
  radius = int(radius/bufRes)
  coin = np.zeros([8,8], dtype = np.uint64)
  times = np.fromfile(fTimes, dtype = np.uint64)
  chans = np.fromfile(fChans, dtype = np.uint8)
  #print "len(times), len(chans)", len(times), len(chans)
  for chan in range(8,16):
    colIdx = np.where(chans==chan)[0]
    for idx in colIdx:
      #print "chans[idx]: %d"%chans[idx]
      #print "few chans: ",chans[idx-3:idx+3]
      #print "few times: ",times[idx-3:idx+3]
      j = idx + 1
      while (j < len(times)) and (chans[j]==heraldChan) and (times[j] - gate <= times[idx]):
        i = idx - 1
        while (i >= 0):
          if (times[i] + radius >= times[idx]) and (chans[idx] != chans[i]) and chans[i] < 8:
            row = chans[i]
            col = chans[idx] % 8
            coin[row, col] += 1
            break
          elif (times[i] + radius <= times[idx]):
            row = heraldChan % 8 #works even if for some reason we had the rows plugged into channels 8-15 of the tagger
            col = chans[idx] % 8
            coin[row, col] += 1
            break
          i -= 1
        j += 1
  return coin
        
  """

def read_lgal_input_fulltrees_withids(folder,lastsnap,file,verbose):
    firstfile = file
    lastfile = file 
    nTrees = 0
    nHalos = 0
    nTreeHalos = numpy.array([],dtype=numpy.int32)
    output_Halos = numpy.array([],dtype=struct_lgalinput)
    output_HaloIDs = numpy.array([],dtype=struct_lgaldbidsinput)
    ifile = file
    filename = folder+'/trees_'+"%03d"%(lastsnap)+'.'+"%d"%(ifile)
    f = open(filename,"rb")
    this_nTrees = numpy.fromfile(f,numpy.int32,1)[0]
    nTrees += this_nTrees
    this_nHalos = numpy.fromfile(f,numpy.int32,1)[0]
    nHalos += this_nHalos
    if(verbose):
        print "File ", ifile," nHalos = ",this_nHalos
    nTreeHalos = numpy.fromfile(f,numpy.int32,this_nTrees)
    output_Halos = numpy.fromfile(f,struct_lgalinput,this_nHalos)
    f.close()
    filename = folder+'/tree_dbids_'+"%03d"%(lastsnap)+'.'+"%d"%(ifile)
    f = open(filename,"rb")
    output_HaloIDs = numpy.fromfile(f,struct_lgaldbidsinput,this_nHalos)
    f.close()
    return (nTrees,nHalos,nTreeHalos,output_Halos,output_HaloIDs)

    def _load_ahf_particle_block(self, f):
        """Load the particles for the next halo described in particle file f"""
        ng = len(self.base.gas)
        nds = len(self.base.dark) + len(self.base.star)
        nparts = int(f.readline().split()[0])

        if self.isnew:
            if isinstance(f, file):
                data = (np.fromfile(
                    f, dtype=int, sep=" ", count=nparts*2).reshape(nparts, 2))[:, 0]
            else:
                # unfortunately with gzipped files there does not
                # seem to be an efficient way to load nparts lines
                data = np.zeros(nparts, dtype=int)
                for i in xrange(nparts):
                    data[i] = int(f.readline().split()[0])

            if self._use_iord :
                data = self._iord_to_fpos[data]
            else :
                hi_mask = data >= nds
                data[np.where(hi_mask)] -= nds
                data[np.where(~hi_mask)] += ng
        else:
            if isinstance(f, file):
                data = np.fromfile(f, dtype=int, sep=" ", count=nparts)
            else:
                # see comment above on gzipped files
                data = np.zeros(nparts, dtype=int)
                for i in xrange(nparts):
                    data[i] = int(f.readline())
        data.sort()
        return data

	def read_from_file(self, filename):
		'''
		Read data from file. Sets the instance variables
		self.raw_velocity and self.kmsrho8
		
		Parameters:
			* filename (string): the file to read from.
		Returns:
			Nothing
		'''
		print_msg('Reading velocity file: %s...' % filename)
		self.filename = filename

		#Read raw data from velocity file
		f = open(filename, 'rb')
		temp_mesh = np.fromfile(f, count=3, dtype='int32')
		self.mesh_x, self.mesh_y, self.mesh_z = temp_mesh
		self.raw_velocity = np.fromfile(f, dtype='float32').astype('float64')
		f.close()
		self.raw_velocity = self.raw_velocity.reshape((3, self.mesh_x, self.mesh_y, self.mesh_z), order='F')

		#Store the redshift from the filename
		try:
			import os.path
			name = os.path.split(filename)[1]
			self.z = float(name.split('v_')[0])
		except:
			print_msg('Could not determine redshift from file name')
			self.z = -1

		#Convert to kms/s*(rho/8)
		self.kmsrho8 = self.raw_velocity*conv.velconvert(z = self.z)


		print_msg('...done')

def openDATfile(filename,ftype,srate=25000):
    fh = open(filename,'r')
    fh.seek(0)
    if ftype == 'amp':
        data = np.fromfile(fh, dtype=np.int16)
        fh.close()
        data = np.double(data)
        data *= 0.195 # according the Intan, the output should be multiplied by 0.195 to be converted to micro-volts
    elif ftype == 'adc':
        data = np.fromfile(fh, dtype=np.uint16)
        fh.close()
        data = np.double(data)
        data *= 0.000050354 # according the Intan, the output should be multiplied by 0.195 to be converted to micro-volts
        data -= np.mean(data)
    
    elif ftype == 'aux':
        data = np.fromfile(fh, dtype=np.uint16)
        fh.close()
        data = np.double(data)
        data *= 0.0000748 # according the Intan, the output should be multiplied by 0.195 to be converted to micro-volts
        
    elif ftype == 'time':
        data = np.fromfile(fh, dtype=np.int32)
        fh.close()
        data = np.double(data)
        data /= srate # according the Intan, the output should be multiplied by 0.195 to be converted to micro-volts
    return data

    def __init__(self, filename, verbose = False):
        super(SpeReader, self).__init__(filename, verbose = verbose)

        # open the file & read the header
        self.header_size = 4100
        self.fileptr = open(filename, "rb")

        self.fileptr.seek(42)
        self.image_width = int(numpy.fromfile(self.fileptr, numpy.uint16, 1)[0])
        self.fileptr.seek(656)
        self.image_height = int(numpy.fromfile(self.fileptr, numpy.uint16, 1)[0])
        self.fileptr.seek(1446)
        self.number_frames = int(numpy.fromfile(self.fileptr, numpy.uint32, 1)[0])

        self.fileptr.seek(108)
        image_mode = int(numpy.fromfile(self.fileptr, numpy.uint16, 1)[0])
        if (image_mode == 0):
            self.image_size = 4 * self.image_width * self.image_height
            self.image_mode = numpy.float32
        elif (image_mode == 1):
            self.image_size = 4 * self.image_width * self.image_height
            self.image_mode = numpy.uint32
        elif (image_mode == 2):
            self.image_size = 2 * self.image_width * self.image_height
            self.image_mode = numpy.int16
        elif (image_mode == 3):
            self.image_size = 2 * self.image_width * self.image_height
            self.image_mode = numpy.uint16
        else:
            print("unrecognized spe image format: ", image_mode)

def sorted_indexes(filename,data_type=np.uint32):

    if os.path.exists(filename):
        curfilename = filename
    elif os.path.exists(filename+".0"):
        curfilename = filename+".0"
    else:
        print "file not found:", filename
        sys.exit()

    f=open(curfilename,'rb')
    number_of_files=np.fromfile(f,dtype=np.uint32,count=1)[0]
    dims=np.fromfile(f,dtype=np.uint32,count=1)[0]
    dims3=dims**3
    total_size=np.fromfile(f,dtype=data_type,count=dims3)
    total_array=[]
    for j in range(dims3):
        total_array.append(np.empty(total_size[j],dtype=data_type))
    f.close()
    total_array=np.array(total_array)

    offset=np.zeros(dims3,dtype=data_type)
    for i in range(number_of_files):
        curfilename=filename+'.'+str(i)
        f=open(curfilename,'rb')
        f.seek(4*(2+dims3),os.SEEK_CUR)
        for j in range(dims3):
            size=np.fromfile(f,dtype=data_type,count=1)[0]
            array=np.fromfile(f,dtype=data_type,count=size)
            total_array[j][offset[j]:offset[j]+size]=array
            offset[j]+=size
        f.close()

    return total_array

def parseData(dataset='testing', path='.'):
	'''
	parseData - Parses a file into matrices
	Input - the name of file to be parsed
	Output - The data in matrix representation

	'''

	if dataset == 'training':
		image_file = os.path.join(path, 'train-images-idx3-ubyte')
		label_file = os.path.join(path, 'train-labels-idx1-ubyte')
	elif dataset == 'testing':
		image_file = os.path.join(path, 't10k-images-idx3-ubyte')
		label_file = os.path.join(path, 't10k-labels-idx1-ubyte')
	else:
		raise(ValueError, "'dataset' must be in testing or 'training'")

	# get the matrix for image data
	f_img = open(image_file, 'rb')
	magic_nr, size = struct.unpack(">II", f_img.read(8))  # parse the magic number, & size of dataset
	dim_x, dim_y = struct.unpack(">II", f_img.read(8))  # get the dimensions of each handwritten num
	X = np.fromfile(f_img, dtype=np.dtype('B'))
	X = X.reshape(size, dim_x * dim_y)


	# get the matrix for label data
	f_lbl = open(label_file, 'rb')
	magic_nr, size = struct.unpack(">II", f_lbl.read(8)) # only magic # and size of dataset
	y = np.fromfile(f_lbl, dtype=np.dtype('B'))
	#X[X > 1] = 1
	return X, y

    def __init__(self, filename, xml):
        DataReader.__init__(self, filename, xml)
        
        # Open the file & read the header.
        self.header_size = 4100
        self.fileptr = open(filename, "rb")

        # FIXME: Should check that these match the XML file.        
        self.fileptr.seek(42)
        self.image_width = int(numpy.fromfile(self.fileptr, numpy.uint16, 1)[0])
        self.fileptr.seek(656)
        self.image_height = int(numpy.fromfile(self.fileptr, numpy.uint16, 1)[0])
        self.fileptr.seek(1446)
        self.number_frames = int(numpy.fromfile(self.fileptr, numpy.uint32, 1)[0])

        self.fileptr.seek(108)
        image_mode = int(numpy.fromfile(self.fileptr, numpy.uint16, 1)[0])
        if (image_mode == 0):
            self.image_size = 4 * self.image_width * self.image_height
            self.image_mode = numpy.float32
        elif (image_mode == 1):
            self.image_size = 4 * self.image_width * self.image_height
            self.image_mode = numpy.uint32
        elif (image_mode == 2):
            self.image_size = 2 * self.image_width * self.image_height
            self.image_mode = numpy.int16
        elif (image_mode == 3):
            self.image_size = 2 * self.image_width * self.image_height
            self.image_mode = numpy.uint16
        else:
            print "unrecognized spe image format: ", image_mode

def readFlowFile(file_name,flip=False):
    data2D=None
    with open(file_name,'rb') as f:
        magic = np.fromfile(f, np.float32, count=1)
        if 202021.25 != magic:
            print 'Magic number incorrect. Invalid .flo file'
        else:
            w = np.fromfile(f, np.int32, count=1)
            h = np.fromfile(f, np.int32, count=1)
            if w.size==0 or h.size==0:
                # print type(w),type(h),w,h
                data2D=None;
            else:               
                # print (w, h)
                data = np.fromfile(f, np.float32, count=2*w*h)
                # Reshape data into 3D array (columns, rows, bands)
                # if flip is True:
                #     data2D = np.resize(data, (w, h, 2))
                #     data2D = data2D.
                #     data2D = np.reshape(data, (h, w, 2))
                #     # ,order='F')
                # else:
                data2D = np.reshape(data, (h, w, 2))
                # print data2D.shape
    return data2D

def getPressureData(fullFileName):
    itemsize = os.path.getsize(fullFileName)/8/2
    fp = open(fullFileName, 'r')
    time = np.fromfile(fp, dtype='Float64', count=itemsize);
    pressure = np.fromfile(fp, dtype='Float64', count=itemsize);
    fp.close()
    return (time, pressure)

def readFLO(path):
    f = open(path, 'rb')
    
    # Read magic number ("PIEH" in ASCII = float 202021.25)
    magic = np.fromfile(f, np.float32, count=1)

    if magic != 202021.25:
        raise Exception('Invalid .flo file') 

    # Read width
    f.seek(4)
    w = np.fromfile(f, np.int32, count=1)
    
    # Read height
    f.seek(8)
    h = np.fromfile(f, np.int32, count=1)
    
    # Read (u,v) coordinates
    f.seek(12)
    data = np.fromfile(f, np.float32, count=w*h*2)

    # Close file (.flo)
    f.close()

    # Reshape data into 3D array (columns, rows, bands)
    dataM = np.resize(data, (h, w, 2))

    # Extract u and v coordinates
    u = dataM[:,:,0]
    v = dataM[:,:,1]
    
    return w,h,u,v

    def test_masked_gauss(self):
        data = numpy.ones((50, 10))
        data[:, 5:] = 2
        lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
        lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
        swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
        mask = numpy.ones((50, 10))
        mask[:, :5] = 0
        masked_data = numpy.ma.array(data, mask=mask)
        res = kd_tree.resample_gauss(swath_def, masked_data.ravel(),
                                     self.area_def, 50000, 25000, segments=1)
        expected_mask = numpy.fromfile(os.path.join(os.path.dirname(__file__),
                                                    'test_files',
                                                    'mask_test_mask.dat'),
                                       sep=' ').reshape((800, 800))
        expected_data = numpy.fromfile(os.path.join(os.path.dirname(__file__),
                                                    'test_files',
                                                    'mask_test_data.dat'),
                                       sep=' ').reshape((800, 800))
        expected = expected_data.sum()
        cross_sum = res.data.sum()

        self.assertTrue(numpy.array_equal(expected_mask, res.mask),
                        msg='Gauss resampling of swath mask failed')
        self.assertAlmostEqual(cross_sum, expected, places=3,
                               msg='Gauss resampling of swath masked data failed')

def load_ply(fn):
    f = open(fn, 'rb')
    prop_to_dtype = { 'float' : np.float32, 'int' : np.int32, 'uchar' : np.uint8 }

    header = []
    while True:
        s = f.readline().split()
        header.append(s)
        if s[0] == 'end_header':
            break

    it = iter(header)
    s = it.next()
    elements = {}
    while True:
        if s[0] == 'end_header':
            break
        if s[0] == 'element':
            el_name, el_len = s[1], int(s[2])
            el_props = []
            s = it.next()
            while s[0] == 'property':
                el_props.append( s )
                s = it.next()
            if el_name == 'face':
                el_type = np.dtype( [('count', np.uint8), ('idx', np.int32, 3)] )
                elements[el_name] = np.fromfile(f, el_type, el_len)['idx'].copy()
            else:
                el_type = np.dtype( [(name, np.dtype(prop_to_dtype[tp])) for _, tp, name in el_props] )
                elements[el_name] = np.fromfile(f, el_type, el_len)
            continue
        s = it.next()
    return elements

def load_matrix_csr(path, verbose=False):
    t_start = time.time()
    data = np.fromfile(
        open(
            os.path.join(
                path,
                "bigrams.data.bin")),
        dtype=np.float32)
    col_ind = np.fromfile(
        open(
            os.path.join(
                path,
                "bigrams.col_ind.bin")),
        dtype=np.int64)
    row_ptr = np.fromfile(
        open(
            os.path.join(
                path,
                "bigrams.row_ptr.bin")),
        dtype=np.int64)
    dim = row_ptr.shape[0] - 1
    cooccurrence = scipy.sparse.csr_matrix(
        (data, col_ind, row_ptr), shape=(
            dim, dim), dtype=np.float32)
    t_end = time.time()
    if verbose:
        print("Matrix loaded in {0:0.2f} sec".format(t_end - t_start))
        print_stats(cooccurrence)
    return cooccurrence

def read(dataset = "training", path = "."):
    """
    Python function for importing the MNIST data set.  It returns an iterator
    of 2-tuples with the first element being the label and the second element
    being a numpy.uint8 2D array of pixel data for the given image.
    """

    if dataset is "training":
        fname_img = os.path.join(path, 'train-images-idx3-ubyte')
        fname_lbl = os.path.join(path, 'train-labels-idx1-ubyte')
    elif dataset is "testing":
        fname_img = os.path.join(path, 't10k-images-idx3-ubyte')
        fname_lbl = os.path.join(path, 't10k-labels-idx1-ubyte')
    else:
        raise ValueError, "dataset must be 'testing' or 'training'"

    # Load everything in some numpy arrays
    with open(fname_lbl, 'rb') as flbl:
        magic, num = struct.unpack(">II", flbl.read(8))
        lbl = np.fromfile(flbl, dtype=np.int8)

    with open(fname_img, 'rb') as fimg:
        magic, num, rows, cols = struct.unpack(">IIII", fimg.read(16))
        img = np.fromfile(fimg, dtype=np.uint8).reshape(len(lbl), rows, cols)

    get_img = lambda idx: (lbl[idx], img[idx])

    # Create an iterator which returns each image in turn
    for i in xrange(len(lbl)):
        yield get_img(i)

    def __parse_data(self):
        print 'Reading ', self._filepath
        f = open(self._filepath, 'rb')

        samp_per_segment = 64
        bytes_per_sample = 2
        channels = 2
        tcd_dtype= 'int16'
        f_size = os.path.getsize(self._filepath)
        segments = f_size / ( samp_per_segment * bytes_per_sample * channels )
        self._progress_bar.setMinimum(0)
        self._progress_bar.setMaximum(segments)
        self._value = 0
        self._progress_bar.setValue(self._value)
        chan1 = numpy.array([], dtype=tcd_dtype)
        chan2 = numpy.array([], dtype=tcd_dtype)
        data  = numpy.zeros((samp_per_segment), dtype=tcd_dtype)
        for seg in xrange(segments):
            self._value = self._value + 1
            self._progress_bar.setValue(self._value)
            data = numpy.fromfile(f, dtype=tcd_dtype, count=samp_per_segment)
            chan1 = numpy.concatenate((chan1, data.copy()) )
            data = numpy.fromfile(f, dtype=tcd_dtype, count=samp_per_segment)
            chan2 = numpy.concatenate((chan2, data.copy()) )

        f.close()

        chan1 = chan1.astype(float) / 2.0**11 * self._prf/2.0 *154000.0 / self._doppler_freq_1/10**3
        chan2 = chan2.astype(float) / 2.0**11 * self._prf/2.0 *154000.0 / self._doppler_freq_2/ 10**3

        self.chan1 = chan1
        self.chan2 = chan2

 def parseDataFile(self, uLong, precision, encoding, dataFile):
   assert uLong in ['uint32', 'uint64']
   assert precision in ['single', 'double']
   assert encoding in ['BigEndian', 'LittleEndian']
   
   require_numpy()
   
   fd = file(dataFile, 'rb')
   
   byteorder = {'LittleEndian': '<', 'BigEndian': '>'}[encoding]
   unsignedLongTypeString = {'uint32': 'u4', 'uint64': 'u8'}[uLong]
   realTypeString = {'single': 'f4', 'double': 'f8'}[precision]
   
   ulongDType = numpy.dtype(byteorder + unsignedLongTypeString)
   floatDType = numpy.dtype(byteorder + realTypeString)
   
   independentGeometry = []
   
   for independentVariable in self.independentVariables:
     size = numpy.fromfile(fd, dtype=ulongDType, count=1)
     independentGeometry.append(size)
     assert size == independentVariable['length']
     a = numpy.fromfile(fd, dtype=floatDType, count=size)
     independentVariable['array'] = a
   
   if len(independentGeometry) == 0:
     independentGeometry.append(1)
   
   for dependentVariable in self.dependentVariables:
     size = numpy.fromfile(fd, dtype=ulongDType, count=1)
     a = numpy.fromfile(fd, dtype=floatDType, count=size)
     assert a.size == size, "Data file %s has incorrect size. Variable '%s' wasn't written completely." % (dataFile, dependentVariable['name'])
     dependentVariable['array'] = a.reshape(*independentGeometry)