def _sample(pts, shape, norm):
        (x, y, z), floor = np.modf(pts.T)
        floor = floor.astype(int)
        ceil = floor + 1
        x[x < 0] = 0
        y[y < 0] = 0
        z[z < 0] = 0

        i000 = np.ravel_multi_index((floor[2], floor[1], floor[0]), shape, mode='clip')
        i100 = np.ravel_multi_index((floor[2], floor[1],  ceil[0]), shape, mode='clip')
        i010 = np.ravel_multi_index((floor[2],  ceil[1], floor[0]), shape, mode='clip')
        i001 = np.ravel_multi_index(( ceil[2], floor[1], floor[0]), shape, mode='clip')
        i101 = np.ravel_multi_index(( ceil[2], floor[1],  ceil[0]), shape, mode='clip')
        i011 = np.ravel_multi_index(( ceil[2],  ceil[1], floor[0]), shape, mode='clip')
        i110 = np.ravel_multi_index((floor[2],  ceil[1],  ceil[0]), shape, mode='clip')
        i111 = np.ravel_multi_index(( ceil[2],  ceil[1],  ceil[0]), shape, mode='clip')

        v000 = (1-x)*(1-y)*(1-z)
        v100 = x*(1-y)*(1-z)
        v010 = (1-x)*y*(1-z)
        v110 = x*y*(1-z)
        v001 = (1-x)*(1-y)*z
        v101 = x*(1-y)*z
        v011 = (1-x)*y*z
        v111 = x*y*z

        allj = np.vstack([i000, i100, i010, i001, i101, i011, i110, i111]).T.ravel()
        data = np.vstack([v000, v100, v010, v001, v101, v011, v110, v111]).T.ravel()

        uniquej = np.unique(allj)
        uniquejdata = np.array([data[allj==j].sum() for j in uniquej])
        
        return uniquej, uniquejdata / float(norm)

    def _compute_multipliers(self, X, y):
        n_samples, n_features = X.shape

        K = self._gram_matrix(X)
        # Solves
        # min 1/2 x^T P x + q^T x
        # s.t.
        #  Gx \coneleq h
        #  Ax = b

        P = cvxopt.matrix(np.outer(y, y) * K)
        q = cvxopt.matrix(-1 * np.ones(n_samples))

        # -a_i \leq 0
        # TODO(tulloch) - modify G, h so that we have a soft-margin classifier
        G_std = cvxopt.matrix(np.diag(np.ones(n_samples) * -1))
        h_std = cvxopt.matrix(np.zeros(n_samples))

        # a_i \leq c
        G_slack = cvxopt.matrix(np.diag(np.ones(n_samples)))
        h_slack = cvxopt.matrix(np.ones(n_samples) * self._c)

        G = cvxopt.matrix(np.vstack((G_std, G_slack)))
        h = cvxopt.matrix(np.vstack((h_std, h_slack)))

        A = cvxopt.matrix(y, (1, n_samples))
        b = cvxopt.matrix(0.0)

        solution = cvxopt.solvers.qp(P, q, G, h, A, b)

        # Lagrange multipliers
        return np.ravel(solution['x'])

    def insert(self,A,B,news=None):
        '''
        Insert two blocks into the center of the cylinder.

        Parameters
        ----------
        A,B : any hashable object
            The scopes of the insert block points.
        news : list of any hashable object, optional
            The new scopes for the points of the cylinder before the insertion.
            If None, the old scopes remain unchanged.
        '''
        aspids,bspids,asrcoords,bsrcoords=[],[],[],[]
        for i,rcoord in enumerate(self.block):
            aspids.append(PID(scope=A,site=i))
            bspids.append(PID(scope=B,site=i))
            asrcoords.append(rcoord-self.translation/2)
            bsrcoords.append(rcoord+self.translation/2)
        if len(self)==0:
            self.pids=aspids+bspids
            self.rcoords=np.vstack([asrcoords,bsrcoords])
        else:
            if news is not None:
                assert len(news)*len(self.block)==len(self)
                self.pids=[PID(scope=scope,site=i) for scope in news for i in range(len(self.block))]
            apids,bpids=self.pids[:len(self)//2],self.pids[len(self)//2:]
            arcoords,brcoords=self.rcoords[:len(self)//2]-self.translation,self.rcoords[len(self)//2:]+self.translation
            self.pids=apids+aspids+bspids+bpids
            self.rcoords=np.vstack([arcoords,asrcoords,bsrcoords,brcoords])
        self.icoords=np.zeros(self.rcoords.shape)
        if np.any(np.asarray(list(self.neighbours.values()))==np.inf):
            self.neighbours={i:length for i,length in enumerate(minimumlengths(self.rcoords,self.vectors,self.nneighbour,Lattice.ZMAX))}

def standardize_polygons_str(data_str):
    """Given a POLYGON string, standardize the coordinates to a 1x1 grid.

    Input : data_str (taken from above)
    Output: tuple of polygon objects
    """
    # find all of the polygons in the letter (for instance an A
    # needs to be constructed from 2 polygons)
    path_strs = re.findall("\(\(([^\)]+?)\)\)", data_str.strip())

    # convert the data into a numpy array
    polygons_data = []
    for path_str in path_strs:
        data = np.array([
            tuple(map(float, x.split())) for x in path_str.strip().split(",")])
        polygons_data.append(data)

    # standardize the coordinates
    min_coords = np.vstack(data.min(0) for data in polygons_data).min(0)
    max_coords = np.vstack(data.max(0) for data in polygons_data).max(0)
    for data in polygons_data:
        data[:, ] -= min_coords
        data[:, ] /= (max_coords - min_coords)

    polygons = []
    for data in polygons_data:
        polygons.append(load_wkt(
            "POLYGON((%s))" % ",".join(" ".join(map(str, x)) for x in data)))

    return tuple(polygons)

def display_layer(X, filename="../images/layer.png"):
    """
    Produces an image, composed of the given N images, patches or neural network weights,
    stored in the array X. Saves it with the given filename.
    :param X: numpy array of size (NxD) — N images, patches or neural network weights
    :param filename: a string, the name of the produced file
    :return: None
    """
    if not isinstance(X, np.ndarray):
        raise TypeError("'X' must be a numpy array")
    N, D = X.shape
    d = get_reshaped_image_size(D)

    if N == 1:
        return X.reshape(d, d, 3)
    divizors = [n for n in range(1, N) if N % n == 0]
    im_sizes = divizors[int(len(divizors) / 2)], int(N / divizors[int(len(divizors) / 2)])
    for i in range(im_sizes[0]):
        # img_row = np.hstack((img_row, np.zeros((d, 1, 3))))
        img_row = np.hstack((np.zeros((d, 1, 3)), np.array(X[i * im_sizes[0], :].reshape(d, d, 3))))
        img_row = np.hstack((img_row, np.zeros((d, 1, 3))))
        for j in range(1, im_sizes[1]):
            img_row = np.hstack((img_row, X[i * im_sizes[1] + j, :].reshape(d, d, 3)))
            img_row = np.hstack((img_row, np.zeros((d, 1, 3))))
        if i == 0:
            img = img_row
        else:
            img = np.vstack((img, img_row))
        img = np.vstack((img, np.zeros((1, img.shape[1], 3))))
    img = np.vstack((np.zeros((1, img.shape[1], 3)), img))
    imsave(filename, img)
    return img

    def train(self, training_set, test_set=None, pool=None, N=200):
        M = map if pool is None else pool.map

        user_items = [[] for i in range(self.nusers)]
        item_users = [[] for i in range(self.nitems)]
        [(user_items[u].append(a), item_users[a].append(u))
         for u, a in training_set]

        if test_set is not None:
            test_user_items = defaultdict(list)
            [test_user_items[u].append(a) for u, a in test_set]
            test_args = [(u, t, user_items[u])
                         for u, t in test_user_items.items()]

        for i in range(10):
            print("Updating users")
            vtv = np.dot(self.V.T, self.V)
            self.U = np.vstack(M(_function_wrapper(self,
                                                   "compute_user_update",
                                                   vtv), user_items))

            print("Updating items")
            utu = np.dot(self.U.T, self.U)
            self.V = np.vstack(M(_function_wrapper(self,
                                                   "compute_item_update",
                                                   utu), item_users))

            # Compute the held out recall.
            if test_set is not None:
                print("Computing held out recall")
                yield np.mean(M(_function_wrapper(self, "compute_recall", N=N),
                                test_args))
            else:
                yield 0.0

    def add(self, key, p):
        """Add a new semantic pointer to the vocabulary.

        The pointer value can be a `.SemanticPointer` or a vector.
        """
        if self.readonly:
            raise ReadonlyError(attr='Vocabulary',
                                msg="Cannot add semantic pointer '%s' to "
                                    "read-only vocabulary." % key)

        if not key[0].isupper():
            raise SpaParseError(
                "Semantic pointers must begin with a capital letter.")
        if not isinstance(p, pointer.SemanticPointer):
            p = pointer.SemanticPointer(p)

        if key in self.pointers:
            raise ValidationError("The semantic pointer %r already exists"
                                  % key, attr='pointers', obj=self)

        self.pointers[key] = p
        self.keys.append(key)
        self.vectors = np.vstack([self.vectors, p.v])

        # Generate vector pairs
        if self.include_pairs and len(self.keys) > 1:
            for k in self.keys[:-1]:
                self.key_pairs.append('%s*%s' % (k, key))
                v = (self.pointers[k] * p).v
                self.vector_pairs = np.vstack([self.vector_pairs, v])

def image_border(rgb, left=0, right=0, top=0, bottom=0, color=[1, 1, 1]):
    orig_shape = rgb.shape

    if left > 0:
        # note: do this every time because it changes throughout
        height, width = rgb.shape[0:2]
        pad = rgb_pad(height, left, color)
        rgb = np.hstack((pad, rgb))

    if right > 0:
        height, width = rgb.shape[0:2]
        pad = rgb_pad(height, right, color)
        rgb = np.hstack((rgb, pad))

    if top > 0:
        height, width = rgb.shape[0:2]
        pad = rgb_pad(top, width, color)
        rgb = np.vstack((pad, rgb))
        assert rgb.shape[0] == height + top

    if bottom > 0:
        height, width = rgb.shape[0:2]
        pad = rgb_pad(bottom, width, color)
        rgb = np.vstack((rgb, pad))
        assert rgb.shape[0] == height + bottom

    assert rgb.shape[0] == orig_shape[0] + top + bottom
    assert rgb.shape[1] == orig_shape[1] + left + right

    return rgb

def generate_delta(params):
    print params
    D = params["D"]
    beta = params["beta"]
    npts = params["npts"]
    print "DOS half-bandwidth : ", D
    Gamma = params["gamma"]
    dos_prec = 0.1 

    omega_grid = np.arange(-2*D,2*D,dos_prec)
    dos_vals = dos_function(omega_grid) 
    data = np.vstack([omega_grid,dos_vals]).transpose()
    np.savetxt("dos.dat",data)

    fermi = lambda w : 1. / (1.+np.exp(beta*w))
    delta_wt = lambda tau, w : -fermi(w) * np.exp(tau*w) * dos_function(w) * Gamma
    delta_f = lambda tau : integrate.quad(lambda w: -fermi(w) * np.exp(tau*w) * dos_function(w) * Gamma, -2*D, 2*D) 

    tau_grid = np.linspace(0,beta,npts+1)
    delta_vals = np.array([delta_f(x)[0] for x in tau_grid])
    data_out = np.vstack([range(npts+1), delta_vals, delta_vals]) 
    fname = "delta_tau.dat"
    np.savetxt("delta_tau.dat", data_out.transpose())
    print "Saved", fname

    kramers_kronig_imag = lambda z : integrate.quad(lambda w: np.imag(dos_function(w) / (1j*z - w)), -2*D, 2*D)
    kramers_kronig_real = lambda z : integrate.quad(lambda w: np.real(dos_function(w) / (1j*z - w)), -2*D, 2*D)
    matsubara_grid = (2*np.arange(0, npts, 1, dtype=np.float) + 1)*np.pi/beta
    delta_iw = np.array([[x, kramers_kronig_real(x)[0], kramers_kronig_imag(x)[0]] for x in matsubara_grid])
    fname = "delta_iw.dat"
    np.savetxt(fname, delta_iw)
    print "Saved", fname

def eventPhase(iyd,n=[1,0],eps=1e-9):
  '''
  Construct a analytic signal like version of the system
  '''
  # Work out the angle of the normal vector
  a = arctan2(n[1],n[0])
  z = (iyd[0,:]+1j*iyd[1,:])*exp(1j*a)
  # Find zero crossings of first co-ordinate
  idx0 = where(logical_and(z.real[:-1]<0,z.real[1:]>0))[0]
  # Produce a two dimensional array which we can be interpolated to 
  # estimate phase by assuming a linear trend between events
  idx2P = idx0-eps
  tme = arange(z.size)
  fv = hstack([vstack([idx0,zeros_like(idx0)]),vstack([idx2P,2*pi*ones_like(idx2P)])])
  fv = array(sorted(fv.T,key=lambda x: x[0])).T
  # Fix edge condition
  if fv[1,0] > pi:
    fv = fv[:,1:]
  # Interpolate phase
  nph = interp1d(*fv,bounds_error=False)(tme)
  # We now have nans on the ends, do a linear interpolation to get average 
  # phase velocity and add on to the ends a linear trend like this
  gdIdx = logical_not(isnan(nph))
  nph[gdIdx]=unwrap(nph[gdIdx])
  m,c = polyfit(tme[gdIdx],nph[gdIdx],1)
  if fv[0,-1]+1!= nph.size:
    nph[fv[0,-1]+1:] = nph[fv[0,-1]]+m*(arange(nph[fv[0,-1]+1:].size)+1)
  if fv[0,0]!=0:
    nph[:fv[0,0]] = nph[fv[0,0]]-m*(nph[:fv[0,0]].size-arange(nph[:fv[0,0]].size))
  return(nph)

    def process_current_split(self):
        # Training of node on training data
        for data, label in self.input_node.request_data_for_training(False):
            self.train(data, label)

        # Compute performance metrics SSNR_AS and SSNR_vs on the training
        # data
        performance = {"ssnr_as" : self.ssnr.ssnr_as(),
                       "ssnr_vs" : self.ssnr.ssnr_vs()}

        # Collect test data (if any)
        X_test = None
        D_test = None
        for data, label in self.input_node.request_data_for_testing():
            if label == self.erp_class_label:
                D = numpy.diag(numpy.ones(data.shape[0]))
            else:
                D = numpy.zeros((data.shape[0], data.shape[0]))

            if X_test is None:
                X_test = deepcopy(data)
                D_test = D
            else:
                X_test = numpy.vstack((X_test, data))
                D_test = numpy.vstack((D_test, D))

        # If there was separate test data:
        # compute metrics that require test data
        if X_test is not None:
            performance["ssnr_vs_test"] = self.ssnr.ssnr_vs_test(X_test, D_test)

        # Add SSNR-based metrics computed in this split to result collection
        self.ssnr_collection.add_split(performance, train=False,
                                       split=self.current_split,
                                       run=self.run_number)

    def r_log_spiral(self,phi):
	"""
	return distance from center for angle phi of logarithmic spiral

	Parameters
	----------
	phi: scalar or np.array with polar angle values

	Returns
	-------
	r(phi) = rx * exp(b * phi) as np.array

	Notes
	-----
	see http://en.wikipedia.org/wiki/Logarithmic_spiral
	"""
	if np.isscalar(phi):
	    phi = np.array([phi])
	ones = np.ones(phi.shape[0])

	# self.rx.shape = 8
	# phi.shape = p
	# then result is given as (8,p)-dim array, each row stands for one rx

	result = np.tensordot(self.rx , np.exp((phi - 3.*pi*ones) / np.tan(pi/2. - self.idisk)),axes = 0)
	result = np.vstack((result, np.tensordot(self.rx , np.exp((phi - pi*ones) / np.tan(pi/2. - self.idisk)),axes = 0) ))
	result = np.vstack((result, np.tensordot(self.rx , np.exp((phi + pi*ones) / np.tan(pi/2. - self.idisk)),axes = 0) ))
	return np.vstack((result, np.tensordot(self.rx , np.exp((phi + 3.*pi*ones) / np.tan(pi/2. - self.idisk)),axes = 0) ))

	def _normals(self, hits, directs):
		"""
		Finds the normal to the parabola in a bunch of intersection points, by
		taking the derivative and rotating it. Used internally by quadric.
		
		Arguments:
		hits - the coordinates of intersections, as an n by 3 array.
		directs - directions of the corresponding rays, n by 3 array.
		"""
		hit = N.dot(N.linalg.inv(self._working_frame), N.vstack((hits.T, N.ones(hits.shape[0]))))
		dir_loc = N.dot(self._working_frame[:3,:3].T, directs.T)

		partial_x = 2.*hit[0]*self.a+self.c*hit[1]+self.d
		partial_y = 2.*hit[1]*self.b+self.c*hit[0]+self.e
		partial_z = -1*N.ones(N.shape(hits)[0])
		
		local_normal = N.vstack((partial_x, partial_y, partial_z))
		local_unit = local_normal/N.sqrt(N.sum(local_normal**2, axis=0))

		down = N.sum(dir_loc * local_unit, axis=0) > 0.
		local_unit[:,down] *= -1

		normals = N.dot(self._working_frame[:3,:3], local_unit)

		return normals

def generate_seq_to_one_hot(X, Y, vocab_size, batch_size):
    n_samples = len(X)
    seq_len = len(X[0])
    start = 0
    while 1:
        stop = start + batch_size
        chunk = X[start: stop]
        slices = []
        for i, seq_indexes in enumerate(chunk):
            x_slice = np.zeros([seq_len, vocab_size])
            x_slice[np.arange(seq_len), seq_indexes] = 1
            slices.append(x_slice)
        x_out = np.stack(slices, axis=0)
        y_out = Y[start: stop]
        start += batch_size
        if (start + batch_size) > n_samples:
            print 'reshuffling, %s + %s > %s' % (start, batch_size, n_samples)
            remaining_X = X[start: start + batch_size]
            remaining_Y = Y[start: start + batch_size]
            random_index = np.random.permutation(n_samples)
            X = np.vstack((remaining_X, X[random_index, :]))
            Y = np.vstack((remaining_Y, Y[random_index, :]))
            start = 0
            n_samples = len(X)
        yield x_out, y_out

        def get_new_cell(self):
            """Returns new basis vectors"""
            a = np.sqrt(self.a)
            b = np.sqrt(self.b)
            c = np.sqrt(self.c)

            ad = self.atoms.cell[0] / np.linalg.norm(self.atoms.cell[0])

            Z = np.cross(self.atoms.cell[0], self.atoms.cell[1])
            Z /= np.linalg.norm(Z)
            X = ad - np.dot(ad, Z) * Z
            X /= np.linalg.norm(X)
            Y = np.cross(Z, X)

            alpha = np.arccos(self.x / (2 * b * c))
            beta = np.arccos(self.y / (2 * a * c))
            gamma = np.arccos(self.z / (2 * a * b))

            va = a * np.array([1, 0, 0])
            vb = b * np.array([np.cos(gamma), np.sin(gamma), 0])
            cx = np.cos(beta)
            cy = (np.cos(alpha) - np.cos(beta) * np.cos(gamma)) \
                / np.sin(gamma)
            cz = np.sqrt(1. - cx * cx - cy * cy)
            vc = c * np.array([cx, cy, cz])

            abc = np.vstack((va, vb, vc))
            T = np.vstack((X, Y, Z))
            return np.dot(abc, T)

def study_redmapper_lrg_3d(hemi='north'):
    # create 3d grid object
    grid = grid3d(hemi=hemi)
    
    # load SDSS data
    sdss = load_sdss_data_both_catalogs(hemi)
    
    # load redmapper catalog
    rm = load_redmapper(hemi=hemi)
    
    # get XYZ positions (Mpc) of both datasets
    x_sdss, y_sdss, z_sdss = grid.xyz_from_radecz(sdss['ra'], sdss['dec'], sdss['z'], applyzcut=False)
    x_rm, y_rm, z_rm = grid.xyz_from_radecz(rm['ra'], rm['dec'], rm['z_spec'], applyzcut=False)
    pos_sdss = np.vstack([x_sdss, y_sdss, z_sdss]).T
    pos_rm = np.vstack([x_rm, y_rm, z_rm]).T

    # build a couple of KDTree's, one for SDSS, one for RM.
    from sklearn.neighbors import KDTree
    tree_sdss = KDTree(pos_sdss, leaf_size=30)
    tree_rm = KDTree(pos_rm, leaf_size=30)

    lrg_counts = tree_sdss.query_radius(pos_rm, 100., count_only=True)
    pl.clf()
    pl.hist(lrg_counts, bins=50)
    
    
    ipdb.set_trace()

    def frame_based(self, annotated_ground_truth, system_output, resolution=0.01):
        # Convert event list into frame-based representation
        system_event_roll = self.list_to_roll(data=system_output, time_resolution=resolution)
        annotated_event_roll = self.list_to_roll(data=annotated_ground_truth, time_resolution=resolution)

        # Fix durations of both event_rolls to be equal
        if annotated_event_roll.shape[0] > system_event_roll.shape[0]:
            padding = numpy.zeros((annotated_event_roll.shape[0] - system_event_roll.shape[0], len(self.class_list)))
            system_event_roll = numpy.vstack((system_event_roll, padding))

        if system_event_roll.shape[0] > annotated_event_roll.shape[0]:
            padding = numpy.zeros((system_event_roll.shape[0] - annotated_event_roll.shape[0], len(self.class_list)))
            annotated_event_roll = numpy.vstack((annotated_event_roll, padding))

        # Compute frame-based metrics
        Nref = sum(sum(annotated_event_roll))
        Ntot = sum(sum(system_event_roll))
        Ntp = sum(sum(system_event_roll + annotated_event_roll > 1))
        Nfp = sum(sum(system_event_roll - annotated_event_roll > 0))
        Nfn = sum(sum(annotated_event_roll - system_event_roll > 0))
        Nsubs = min(Nfp, Nfn)

        eps = numpy.spacing(1)

        results = dict()
        results['Rec'] = Ntp / (Nref + eps)
        results['Pre'] = Ntp / (Ntot + eps)
        results['F'] = 2 * ((results['Pre'] * results['Rec']) / (results['Pre'] + results['Rec'] + eps))
        results['AEER'] = (Nfn + Nfp + Nsubs) / (Nref + eps)

        return results

 def offsetPlane(plane, x, y):
     """
     Takes a numpy 2D array and returns the same plane offset by x and y,
     adding rows and columns of 0 values
     """
     height, width = plane.shape
     dataType = plane.dtype
     # shift x by cropping, creating a new array of columns and stacking
     # horizontally
     if abs(x) > 0:
         newCols = zeros((height, abs(x)), dataType)
         x1 = max(0, 0 - x)
         x2 = min(width, width - x)
         crop = plane[0:height, x1:x2]
         if x > 0:
             plane = hstack((newCols, crop))
         else:
             plane = hstack((crop, newCols))
     # shift y by cropping, creating a new array of rows and stacking
     # vertically
     if abs(y) > 0:
         newRows = zeros((abs(y), width), dataType)
         y1 = max(0, 0 - y)
         y2 = min(height, height - y)
         crop = plane[y1:y2, 0:width]
         if y > 0:
             plane = vstack((newRows, crop))
         else:
             plane = vstack((crop, newRows))
     return plane

def upsample_swc(swc):

    tswc = swc.copy()

    id_idx = {}
    # Build a nodeid->idx hash table
    for nodeidx in range(tswc.shape[0]):
        id_idx[tswc[nodeidx, 0]] = nodeidx

    newid = tswc[:,0].max() + 1
    newnodes = []
    for nodeidx in range(tswc.shape[0]):
        pid = tswc[nodeidx, -1] # parent id

        if pid not in id_idx:
            # raise Exception('Parent with id %d not found' % pid)
            continue

        nodepos = tswc[nodeidx, 2:5]
        parentpos = tswc[id_idx[pid], 2:5]

        if np.linalg.norm(nodepos - parentpos) > 1.: # Add a node in the middle if too far
            mid_pos = nodepos + 0.5 * (parentpos - nodepos)
            newnodes.append( np.asarray([newid, 2, mid_pos[0], mid_pos[1], mid_pos[2], 1, pid]) )
            newid += 1
            tswc[nodeidx, -1] = newid

    # Stack the new nodes to the end of the swc file
    newnodes = np.vstack(newnodes)
    tswc = np.vstack((tswc, newnodes))
    return tswc

def dobetterstuff(inpath):
    data_files = [f for f in os.listdir(inpath) if f.endswith('.mel.npy')]
    random.shuffle(data_files)
    artists = set([f[:18] for f in data_files])
    artist_string_to_id = dict([(s,i) for i, s in enumerate(artists)])

    def get_split(datafiles___, splitpercent):
        # gen = filtered_stratified_split(datafiles___,
        #                                 sklearn.cross_validation.StratifiedShuffleSplit,
        #                                 [1] * len(datafiles___), n_iterations=1, test_size=splitpercent)
        gen = sklearn.cross_validation.ShuffleSplit(len(datafiles___), 1, splitpercent)
        for i_trs, i_tes in gen:
            return [datafiles___[i] for i in i_trs],  [datafiles___[i] for i in i_tes]

    training_files, test_files =  get_split(data_files, .2)
    training_files, validation_files = get_split(training_files, .2)

    print training_files
    print test_files
    print validation_files

    train_set_y = np.hstack([[artist_string_to_id[f[:18]]] * 129 for f in training_files])
    train_set_x = np.vstack([np.load(os.path.join(inpath, f)) for f in training_files])
    test_set_y = np.hstack([[artist_string_to_id[f[:18]]] * 129 for f in test_files])
    test_set_x = np.vstack([np.load(os.path.join(inpath, f)) for f in test_files])
    validation_set_y = np.hstack([[artist_string_to_id[f[:18]]] * 129 for f in validation_files])
    validation_set_x = np.vstack([np.load(os.path.join(inpath, f)) for f in validation_files])

    datasets = [(train_set_x, train_set_y), (validation_set_x, validation_set_y), (test_set_x, test_set_y)]
    return datasets