def detectorgeometry(ddef, npdtype='float32'):
    """ calculates the geometric properties of the detector
        from its definitions
        parameters:
            ddef is array detector definitions

        returns:
            pixelpositions    The endpoint of all the vectors
            unitnormalvector  of the detector  """
    c0 = np.array([ddef[0],ddef[1],ddef[2]]) # corner c1-c0-c2
    c1 = np.array([ddef[3],ddef[4],ddef[5]]) # corner c0-c1-c3 or 1st axis endposition
    c2 = np.array([ddef[6],ddef[7],ddef[8]]) # corner c0-c2-c3 or 2nd axis endposition
    r1 = ddef[9]   # resolution in 1st dimension
    r2 = ddef[10]  # resolution in 2nd dimension

    dshape = (r2,r1)  # CONTROL of SEQUENCE
    # unit direction vectors of detector sides
    di = (c1 - c0) * (1.0/r1)
    dj = (c2 - c0) * (1.0/r2)
    def pcfun(j,i,k): return  c0[k] + (i+0.5) * di[k] + (j+0.5) * dj[k]
    def pcfunx(j,i):  return pcfun(j,i,0)
    def pcfuny(j,i):  return pcfun(j,i,1)
    def pcfunz(j,i):  return pcfun(j,i,2)
    pxs  = numpy.fromfunction(pcfunx, shape=dshape, dtype=npdtype )
    pys  = numpy.fromfunction(pcfuny, shape=dshape, dtype=npdtype )
    pzs  = numpy.fromfunction(pcfunz, shape=dshape, dtype=npdtype )
    pixelpositions = np.array(numpy.dstack((pxs,pys,pzs))) # shape = (r2,r1,3)

    pixelareavector = np.array(numpy.cross(di, dj))
    result = np.zeros(dshape)

    return  pixelpositions, pixelareavector, dshape, result

    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 test_masked_multi_from_sample(self):
     data = numpy.ones((50, 10))
     data[:, 5:] = 2
     mask1 = numpy.ones((50, 10))
     mask1[:, :5] = 0
     mask2 = numpy.ones((50, 10))
     mask2[:, 5:] = 0
     mask3 = numpy.ones((50, 10))
     mask3[:25, :] = 0
     data_multi = numpy.column_stack(
         (data.ravel(), data.ravel(), data.ravel()))
     mask_multi = numpy.column_stack(
         (mask1.ravel(), mask2.ravel(), mask3.ravel()))
     masked_data = numpy.ma.array(data_multi, mask=mask_multi)
     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)
     valid_input_index, valid_output_index, index_array, distance_array = \
         kd_tree.get_neighbour_info(swath_def,
                                    self.area_def,
                                    50000, neighbours=1, segments=1)
     res = kd_tree.get_sample_from_neighbour_info('nn', (800, 800),
                                                  masked_data,
                                                  valid_input_index,
                                                  valid_output_index, index_array,
                                                  fill_value=None)
     expected_fill_mask = numpy.fromfile(os.path.join(os.path.dirname(__file__),
                                                      'test_files',
                                                      'mask_test_full_fill_multi.dat'),
                                         sep=' ').reshape((800, 800, 3))
     fill_mask = res.mask
     self.assertTrue(numpy.array_equal(fill_mask, expected_fill_mask),
                     msg='Failed to create fill mask on masked data')

 def ListPlanes_Cubic(self,latticeConstant,k):
     """
     According to the apparatus geometry, make a list of avaible
     reflection planes.
     """
     
     #wave length
     lamb = 2*pi/k
     maxtheta = arctan(self.L/sqrt(2)/self.D)/2
     
     #Bragg's law
     d = latticeConstant
     nmax = int(2*d*sin(maxtheta)/lamb)
     assert nmax>0
     N = nmax*2+1
     nsize = [N, N, N]
     h = fromfunction(lambda a,b,c: (a-nmax), nsize).astype(int)
     k = fromfunction(lambda a,b,c: (b-nmax), nsize).astype(int)
     l = fromfunction(lambda a,b,c: (c-nmax), nsize).astype(int)
     
     #Interplanar spacing calculated
     ds = d/sqrt(h**2+k**2+l**2)
     hkl = {} 
     thetas = arcsin(lamb/2/ds)
     
     #Check for allowed hkls and put in a list
     for i in range(N):

    def test_custom(self):
        def wf(dist):
            return 1 - dist / 100000.0

        data = numpy.fromfunction(lambda y, x: (y + x) * 10 ** -5, (5000, 100))
        lons = numpy.fromfunction(
            lambda y, x: 3 + (10.0 / 100) * x, (5000, 100))
        lats = numpy.fromfunction(
            lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100))
        swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
        if (sys.version_info < (2, 6) or
                (sys.version_info >= (3, 0) and sys.version_info < (3, 4))):
            res = kd_tree.resample_custom(swath_def, data.ravel(),
                                          self.area_def, 50000, wf, segments=1)
        else:
            with warnings.catch_warnings(record=True) as w:
                res = kd_tree.resample_custom(swath_def, data.ravel(),
                                              self.area_def, 50000, wf, segments=1)
                self.assertFalse(
                    len(w) != 1, 'Failed to create neighbour radius warning')
                self.assertFalse(('Possible more' not in str(
                    w[0].message)), 'Failed to create correct neighbour radius warning')
        cross_sum = res.sum()
        expected = 4872.81050729
        self.assertAlmostEqual(cross_sum, expected,
                               msg='Swath custom resampling failed')

 def test_gauss_multi_uncert(self):
     data = numpy.fromfunction(lambda y, x: (y + x)*10**-6, (5000, 100))        
     lons = numpy.fromfunction(lambda y, x: 3 + (10.0/100)*x, (5000, 100))
     lats = numpy.fromfunction(lambda y, x: 75 - (50.0/5000)*y, (5000, 100))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     data_multi = numpy.column_stack((data.ravel(), data.ravel(),\
                                      data.ravel()))
     if sys.version_info < (2, 6):
         res, stddev, counts = kd_tree.resample_gauss(swath_def, data_multi,\
                                             self.area_def, 50000, [25000, 15000, 10000], 
                                             segments=1, with_uncert=True)
     else:
         with warnings.catch_warnings(record=True) as w:
             res, stddev, counts = kd_tree.resample_gauss(swath_def, data_multi,\
                                                 self.area_def, 50000, [25000, 15000, 10000], 
                                                 segments=1, with_uncert=True)
             self.failIf(len(w) != 1, 'Failed to create neighbour radius warning')
             self.failIf(('Possible more' not in str(w[0].message)), 'Failed to create correct neighbour radius warning') 
     cross_sum = res.sum()
     cross_sum_stddev = stddev.sum()
     cross_sum_counts = counts.sum()
     expected = 1461.84313918
     expected_stddev = 0.446204424799
     expected_counts = 4934802.0
     self.assertTrue(res.shape == stddev.shape and stddev.shape == counts.shape and counts.shape == (800, 800, 3))
     self.assertAlmostEqual(cross_sum, expected,
                             msg='Swath multi channel resampling gauss failed on data')
     self.assertAlmostEqual(cross_sum_stddev, expected_stddev,
                             msg='Swath multi channel resampling gauss failed on stddev')
     self.assertAlmostEqual(cross_sum_counts, expected_counts,
                             msg='Swath multi channel resampling gauss failed on counts')

 def __init__(self, angle, col, dangle=5, parent=None):
     super().__init__(parent)
     color = QColor(0, 0, 0) if col else QColor(128, 128, 128)
     angle_d = np.rad2deg(angle)
     angle_2 = 90 - angle_d - dangle
     angle_1 = 270 - angle_d + dangle
     dangle = np.deg2rad(dangle)
     arrow1 = pg.ArrowItem(
         parent=self, angle=angle_1, brush=color, pen=pg.mkPen(color)
     )
     arrow1.setPos(np.cos(angle - dangle), np.sin(angle - dangle))
     arrow2 = pg.ArrowItem(
         parent=self, angle=angle_2, brush=color, pen=pg.mkPen(color)
     )
     arrow2.setPos(np.cos(angle + dangle), np.sin(angle + dangle))
     arc_x = np.fromfunction(
         lambda i: np.cos((angle - dangle) + (2 * dangle) * i / 120.),
         (121,), dtype=int
     )
     arc_y = np.fromfunction(
         lambda i: np.sin((angle - dangle) + (2 * dangle) * i / 120.),
         (121,), dtype=int
     )
     pg.PlotCurveItem(
         parent=self, x=arc_x, y=arc_y, pen=pg.mkPen(color), antialias=False
     )

def vertical_flip_map(rows, cols):
    map_x = numpy.fromfunction(
        lambda r, c: c, (rows, cols), dtype=numpy.float32)
    map_y = numpy.fromfunction(
        lambda r, c: rows - 1 - r, (rows, cols), dtype=numpy.float32)

    return map_x, map_y

def horizontal_flip_map(rows, cols):
    map_x = numpy.fromfunction(
        lambda r, c: cols - 1 - c, (rows, cols), dtype=numpy.float32)
    map_y = numpy.fromfunction(
        lambda r, c: r, (rows, cols), dtype=numpy.float32)

    return map_x, map_y

    def test_custom_multi(self):
        def wf1(dist):
            return 1 - dist / 100000.0

        def wf2(dist):
            return 1

        def wf3(dist):
            return np.cos(dist) ** 2

        data = np.fromfunction(lambda y, x: (y + x) * 10 ** -6, (5000, 100))
        lons = np.fromfunction(
            lambda y, x: 3 + (10.0 / 100) * x, (5000, 100))
        lats = np.fromfunction(
            lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100))
        swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
        data_multi = np.column_stack((data.ravel(), data.ravel(),
                                      data.ravel()))
        with catch_warnings(UserWarning) as w:
            res = kd_tree.resample_custom(swath_def, data_multi,
                                          self.area_def, 50000, [wf1, wf2, wf3], segments=1)
            self.assertFalse(len(w) != 1)
            self.assertFalse('Possible more' not in str(w[0].message))
        cross_sum = res.sum()
        expected = 1461.8428378742638
        self.assertAlmostEqual(cross_sum, expected)

def planeRegression(imgarray):
	"""
	performs a two-dimensional linear regression and subtracts it from an image
	essentially a fast high pass filter
	"""
	size = (imgarray.shape)[0]
	count = float((imgarray.shape)[0]*(imgarray.shape)[1])
	def retx(y,x):
		return x
	def rety(y,x):
		return y
	xarray = numpy.fromfunction(retx, imgarray.shape)
	yarray = numpy.fromfunction(rety, imgarray.shape)
	xsum = float(xarray.sum())
	xsumsq = float((xarray*xarray).sum())
	ysum = xsum
	ysumsq = xsumsq
	xysum = float((xarray*yarray).sum())
	xzsum = float((xarray*imgarray).sum())
	yzsum = float((yarray*imgarray).sum())
	zsum = imgarray.sum()
	zsumsq = (imgarray*imgarray).sum()
	xarray = xarray.astype(numpy.float32)
	yarray = yarray.astype(numpy.float32)
	leftmat = numpy.array( [[xsumsq, xysum, xsum], [xysum, ysumsq, ysum], [xsum, ysum, count]] )
	rightmat = numpy.array( [xzsum, yzsum, zsum] )
	resvec = linalg.solve(leftmat,rightmat)
	#print " ... plane_regress: x-slope:",round(resvec[0]*size,5),\
	#	", y-slope:",round(resvec[1]*size,5),", xy-intercept:",round(resvec[2],5)
	newarray = imgarray - xarray*resvec[0] - yarray*resvec[1] - resvec[2]
	#del imgarray,xarray,yarray,resvec
	return newarray

    def test_gauss_multi_uncert(self):
        data = np.fromfunction(lambda y, x: (y + x) * 10 ** -6, (5000, 100))
        lons = np.fromfunction(
            lambda y, x: 3 + (10.0 / 100) * x, (5000, 100))
        lats = np.fromfunction(
            lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100))
        swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
        data_multi = np.column_stack((data.ravel(), data.ravel(),
                                      data.ravel()))
        with catch_warnings(UserWarning) as w:
            # The assertion below checks if there is only one warning raised
            # and whether it contains a specific message from pyresample
            # On python 2.7.9+ the resample_gauss method raises multiple deprecation warnings
            # that cause to fail, so we ignore the unrelated warnings.
            res, stddev, counts = kd_tree.resample_gauss(swath_def, data_multi,
                                                         self.area_def, 50000, [
                                                             25000, 15000, 10000],
                                                         segments=1, with_uncert=True)
            self.assertTrue(len(w) >= 1)
            self.assertTrue(
                any(['Possible more' in str(x.message) for x in w]))
        cross_sum = res.sum()
        cross_sum_counts = counts.sum()
        expected = 1461.8429990248171
        expected_stddev = [0.44621800779801657, 0.44363137712896705,
                           0.43861019464274459]
        expected_counts = 4934802.0
        self.assertTrue(res.shape == stddev.shape and stddev.shape ==
                        counts.shape and counts.shape == (800, 800, 3))
        self.assertAlmostEqual(cross_sum, expected)

        for i, e_stddev in enumerate(expected_stddev):
            cross_sum_stddev = stddev[:, :, i].sum()
            self.assertAlmostEqual(cross_sum_stddev, e_stddev)
        self.assertAlmostEqual(cross_sum_counts, expected_counts)

def test_2():
    Temp1, Temp2 = Temp, Temp/temp_ratio
    double_temp = np.sqrt(Temp1*Temp2)
    ss_temp = np.sqrt(Temp1) + np.sqrt(Temp2)
    delta_temp = np.abs(Temp2 - Temp1)

    Rho1, Rho2 = 2*Rho*np.sqrt(Temp2)/ss_temp, 2*Rho*np.sqrt(Temp1)/ss_temp
    f = Rho1/(np.pi*Temp1)**1.5 * np.fromfunction(lambda i,j,k: np.exp(-sqr_xi(i,j,k,Speed)/(Temp1))*dxi(i,j,k), dimension)
    negative = np.fromfunction(lambda i,j,k: j<N_y, dimension)
    f[negative] = Rho2/(np.pi*Temp2)**1.5 * np.fromfunction(lambda i,j,k: np.exp(-sqr_xi(i,j,k,-Speed)/(Temp2))*dxi(i,j,k), dimension)[negative]

    rho, temp, speed, qflow, tau = calc_macro(f)

    Rho_ = Rho
    Temp_ = double_temp * (1 + 8./3*np.dot(Speed,Speed)/ss_temp**2)
    Speed_ = -Speed * delta_temp / ss_temp**2
    Qflow_ = [ 0, 2*Rho*double_temp*delta_temp/ss_temp/np.sqrt(np.pi) * (1 + 2*np.dot(Speed,Speed)/ss_temp**2), 0 ]
    Tau_ = [ 0, 0, 4*Rho*Speed[0]*double_temp/ss_temp/np.sqrt(np.pi) ]

    #splot(f)
    #splot(f*c[0]*c[1])

    print "\n-- Test #2: free molecular flows - sum of 2 half-Maxwellians"
    print "rho =", err(Rho_, rho)
    print "temp =", err(Temp_, temp)
    print "speed =", err(Speed_, speed)
    print "qflow =", err(Qflow_, qflow/rho)
    print "tau =", err(Tau_, tau/rho)

def __set_pm__(L, window, taps, fwidth):
    global pm
    if pm['L'] == L and pm['taps'] == taps and pm['fwidth'] == fwidth:
        if type(window) == str and pm['window_name'] == window: return
        elif window is pm['window']: return
    else:
        pm['L'] = L
        pm['taps'] = taps
        pm['fwidth'] = fwidth
        def sinx_x(x):
            t = n.pi * taps * fwidth * (x/float(L) - .5)
            v = n.where(t != 0, t, 1)
            return n.where(t != 0, n.sin(v) / v, 1)
        pm['sinx_x'] = n.fromfunction(sinx_x, (L,))
    if type(window) == str:
        wf = {}
        wf['hamming'] = lambda x: .54 + .46 * cos(2*n.pi*x/L - n.pi)
        wf['hanning'] = lambda x: .5 + .5 * n.cos(2*n.pi*x/(L+1) - n.pi)
        wf['none'] = lambda x: 1
        pm['window'] = n.fromfunction(wf[window], (L,))
        pm['window_name'] = window
    else:
        pm['window'] = window
        pm['window_name'] = None
    pm['window_sinx_x'] = pm['window'] * pm['sinx_x']

def __set_pm__(L, window, taps, fwidth):
    global pm
    if pm['L'] == L and pm['taps'] == taps and \
       pm['fwidth'] == fwidth:
        if type(window) == str and pm['window_name'] == window: return
        elif window is pm['window']: return
    else:
        pm['L'] = L
        pm['taps'] = taps
        pm['fwidth'] = fwidth
        def sinx_x(x):
            t = n.pi * taps * fwidth * (x/float(L) - .5)
            v = n.where(t != 0, t, 1)
            return n.where(t != 0, n.sin(v) / v, 1)
        pm['sinx_x'] = n.fromfunction(sinx_x, (L,))
    if type(window) == str: 
        wf = {}
        wf['blackman'] = lambda x: .42-.5*n.cos(2*n.pi*x/(L-1))+.08*n.cos(4*n.pi*x/(L-1))
        wf['blackman-harris'] = lambda x: .35875 - .48829*n.cos(2*n.pi*x/(L-1)) + .14128*n.cos(4*n.pi*x/(L-1)) - .01168*n.cos(6*n.pi*x/(L-1))
        wf['gaussian0.4'] = lambda x: n.exp(-0.5 * ((x - (L-1)/2)/(0.4 * (L-1)/2))**2)
        wf['kaiser2'] = lambda x: i0(n.pi * 2 * n.sqrt(1-(2*x/(L-1) - 1)**2)) / i0(n.pi * 2)
        wf['kaiser3'] = lambda x: i0(n.pi * 3 * n.sqrt(1-(2*x/(L-1) - 1)**2)) / i0(n.pi * 3)
        wf['hamming'] = lambda x: .54 - .46 * n.cos(2*n.pi*x/(L-1))
        wf['hanning'] = lambda x: .5 - .5 * n.cos(2*n.pi*x/(L-1))
        wf['parzen'] = lambda x: 1 - n.abs(L/2. - x) / (L/2.)
        wf['none'] = lambda x: 1
        pm['window'] = n.fromfunction(wf[window], (L,))
        pm['window_name'] = window
    else: 
        pm['window'] = window
        pm['window_name'] = None
    pm['window_sinx_x'] = pm['window'] * pm['sinx_x']

	def kerph2im(self,theta=0):
		"""
		adds up the sine transform for uv phases & then multiplies Kmat,
		the transfer matrix from uv phases to kernel phases.
	
		Returns image to kernel phase transfer matrix.
		"""
		# To make the image pixel centered:
		self.off = np.array([self.fov/2.-0.5, self.fov/2.-0.5])
		# empty sine transform matrix:
		self.ph2im = np.zeros((len(self.uv), self.fov,self.fov))
		self.sym2im = np.zeros((len(self.uv), self.fov,self.fov))
		rotMatrix = np.array([[np.cos(theta), -np.sin(theta)], 
						 [np.sin(theta),  np.cos(theta)]])
		for q,one_uv in enumerate(self.uv):
			self.rcoord = np.dot(rotMatrix,one_uv)
			self.ph2im[q,:,:] = self.red[q]*np.fromfunction(self.ffs, (self.fov, self.fov))
			self.sym2im[q,:,:] = self.red[q]*np.fromfunction(self.ffc, (self.fov,self.fov))
		# flatten for matrix multiplication
		#self.ph2im = self.ph2im.reshape(len(self.uv), self.fov*self.fov)
		# Matrix multiply Kmat & flattened sin transform matrix
		self.kerim = np.dot(self.Kmat, self.ph2im.reshape(len(self.uv), self.fov*self.fov))
		# Reshape back to image dimensions
		self.kerim = self.kerim.reshape(len(self.Kmat), self.fov,self.fov)
		return self.kerim, self.sym2im

def multi_2_func():
    """
    使用函数创建2维数组
    :return: 
    """
    print np.fromfunction(lambda x, y: (x + 1) * y, (10, 5))
    split_line()

def make_gauss_init(pkinit, boxsize=100.0, ngrid=100, seed=314159, exactpk=True, smw=2.0):
    
    #mathematical artifact of Fourier transform, possibly incorrect? (ngrid+1)/2)
    thirdim = ngrid/2+1
    
    
    #inverse(?) of the cell size (in 2pi)
    kmin = 2*np.pi/np.float(boxsize)
    
    #shape of k grid
    sk = (ngrid,ngrid,thirdim)
    
    #No clue. it represents the grid, but the method / definition is non-trivial.
    a = np.fromfunction(lambda x,y,z:x, sk).astype(np.float)
    a[np.where(a > ngrid/2)] -= ngrid
    b=np.transpose(a, (1, 0, 2))
    c = np.fromfunction(lambda x,y,z:z, sk).astype(np.float)
    c[np.where(c > ngrid/2)] -= ngrid
    
    #From this, a/b/c = x/y/z but why are 3x10x10x6 arrays necessary?
    kgrid = kmin*np.sqrt(a**2+b**2+c**2).astype(np.float)
    a = 0
    b = 0
    c = 0
    #K-space is a bitch.
    
    
    #self-explanatory
    rs = np.random.RandomState(seed)
    
    #What are these complex numbers representing?
    if (exactpk):
        dk = np.exp(2j*np.pi*rs.rand(ngrid*ngrid*(thirdim))).reshape(sk).astype(np.complex64)
    else:
        dk = np.empty(sk,dtype=np.complex64)
        dk.real = rs.normal(size=ngrid*ngrid*(thirdim)).reshape(sk).astype(np.float32)
        dk.imag = rs.normal(size=ngrid*ngrid*(thirdim)).reshape(sk).astype(np.float32)
        dk /= np.sqrt(2.0)
    #This section provides a unit length, random(?) phase complex number associated with
    #each point on the grid.
    
    #Gaussian filter
    filt=np.exp(-kgrid**2*smw**2)
    
    #interpolate power spectrum to approximate at non-keyed values.
    pkinterp=ip.interp1d(pkinit[:,0], pkinit[:,1])
    
    #Pk undefined at 0, so this provides an explicit exception.
    if (kgrid[0,0,0]==0):
        dk[0,0,0]=0
        wn0=np.where(kgrid!=0)
        
        #The mathematics? 
        dk[wn0] *= np.sqrt(filt[wn0]*pkinterp(kgrid[wn0]))*ngrid**3/boxsize**1.5   
    #Why is this if statement necessary? We know the P(k) will be undefined.
    else:
        dk *= np.sqrt(filt*pkinterp(kgrid.flatten())).reshape(sk)*ngrid**3/boxsize**1.5
        
    dk=nyquist(dk)
    return dk

    def load_cifar(self,data_path):
        image_list = []
        label_list = []
        for i in xrange(5):
            fo = open(os.path.join(data_path,'data_batch_'+str(i+1)),'rb')
            dict = cPickle.load(fo)
            fo.close()
            image_list.append(dict['data'])
            label_list.append(dict['labels'])
        train_images = np.r_[image_list[0],image_list[1],image_list[2],image_list[3],image_list[4]]
        train_labels = np.r_[label_list[0],label_list[1],label_list[2],label_list[3],label_list[4]]
        fo = open(os.path.join(data_path,'test_batch'), 'rb')
        dict= cPickle.load(fo)
        fo.close()
        test_images = dict['data']
        test_labels = dict['labels']

        test_images = test_images / 255.0
        train_images = train_images / 255.0

        # resize images to 2D
        dsize = 32
        train_images = train_images.reshape(len(train_images),3,dsize,dsize)
        test_images = test_images.reshape(len(test_images),3,dsize,dsize)
        train_labels = np.array(train_labels)
        test_labels = np.array(test_labels)
        
        self.test_images = test_images
        self.test_labels = np.fromfunction(lambda i,j:j==test_labels[i],(len(test_labels),max(test_labels)+1),dtype=int)+0
        self.train_images = train_images
        self.train_labels = np.fromfunction(lambda i,j:j==train_labels[i],(len(train_labels),max(train_labels)+1),dtype=int)+0

def _Slices(Beta, legendre_orders, smoothing=0):
    """Convolve Beta with a Gaussian function of 1/e width smoothing.

    """

    pol = len(legendre_orders)
    NP = len(Beta[0])  # number of points in 3_d plot.
    index = range(NP)

    Beta_convol = np.zeros((pol, NP))
    Slice_3D = np.zeros((pol, 2*NP, 2*NP))

    # Convolve Beta's with smoothing function
    if smoothing > 0:
        # smoothing function
        Basis_s = np.fromfunction(lambda i: np.exp(-(i - (NP)/2)**2 /
                                  (2*smoothing**2))/(smoothing*2.5), (NP,))
        for i in range(pol):
            Beta_convol[i] = np.convolve(Basis_s, Beta[i], mode='same')
    else:
        Beta_convol = Beta

    # Calculate ordered slices:
    for i in range(pol):
        Slice_3D[i] = np.fromfunction(lambda k, l: _SL(i, (k-NP), (l-NP),
                       Beta_convol, index, legendre_orders), (2*NP, 2*NP))
    # Sum ordered slices up
    Slice = np.sum(Slice_3D, axis=0)

    return Slice, Beta_convol