def calc_eigs(self, sort=True, symmetric=False, mandel=False):
        if symmetric:
            eigfun=np.linalg.eigvalsh
        else:
            eigfun=np.linalg.eigvals

        if self.order==2:
            eigs=np.zeros(self.shape[-1]+self.N)
            for ind in np.ndindex(self.N):
                mat=self.val[:, :][ind]
                eigs.append(eigfun(mat))
        elif self.order==4:
            if mandel:
                matrixfun=lambda x: ElasticTensor.create_mandel(x)
                d=self.shape[2]
                eigdim=d*(d+1)/2
            else:
                matshape=(self.shape[0]*self.shape[1], self.shape[2]*self.shape[3])
                matrixfun=lambda x: np.reshape(val[ind], matshape)
                eigdim=self.shape[2]*self.shape[3]

            eigs=np.zeros(self.N+(eigdim,))
            val=np.copy(self.val)
            for ii in range(self.dim):
                val=np.rollaxis(val, self.val.ndim-self.dim+ii, ii)

            for ind in np.ndindex(*self.N):
                mat=matrixfun(val[ind])
                eigs[ind]=eigfun(mat)

        eigs=np.array(eigs)
        if sort:
            eigs=np.sort(eigs, axis=0)
        return eigs

def garnet_gen_s( Ns, Na, Nb, sparsity, neighbor):
    ### generating the Kernel
    kernel = np.zeros((Ns, Ns, Na))  # p(s'|s,a)
    for i, j in np.ndindex((Ns, Na)):
        echantillon = rd.sample(list(set(range(Ns)).intersection(range(i-neighbor,i+neighbor))), Nb)
        cumulative = np.concatenate(([0], sort([rd.random() for k in range(Nb - 1)]), [1]), axis=0)
        for k in range(Nb):
            kernel[echantillon[k], i, j] = cumulative[k + 1] - cumulative[k]
    ### generating rewards at random
    reward0 = np.zeros((Ns, Na))
    reward1 = np.zeros((Ns, Na))

    biais0 = np.random.randn(Ns)
    biais1 = np.random.randn(Ns)
    for i, j in np.ndindex((Ns, Na)):
        reward0[i,j] = biais0[i]
        reward1[i,j] = biais1[i]

    masque_reward = np.zeros((Ns, Na))
    N_sparsity = int(Ns * sparsity)
    i = 0
    while i < N_sparsity:
        i_ = rd.randint(0, Ns - 1)
        if masque_reward[i_, 0] == 0:
            masque_reward[i_, :] = 1
            i += 1
    reward0 = reward0 * masque_reward
    reward1 = reward1 * masque_reward
    control = np.random.randint(2, size=Ns)

    return Ns, Na, kernel, reward0, reward1, control

 def __call__(self, distinguish=False):
     """Get all positions as a list. If *distinguish* is *True*, nest the
     position list inside a 1-tuple (as the positions cannot be
     distinguished any further in this case)."""
     # 2012-05-03 - 2012-09-03
     return (list(numpy.ndindex(self.shape)),) if distinguish \
         else list(numpy.ndindex(self.shape))

def ndmeshgrid(grids, hnode=None):
    """
    Converts a list of (start, stop, n) tuples to an 'n-dimensional meshgrid'.
    In two dimensions, this would be:
    
        x = linspace(*grids[0])
        y = linspace(*grids[1])
        x,y = meshgrid(x,y)
        z = concatenate(x,y,axis=-1)
    
    or something like that. Also returns the number of locations in each direction
    as a list.
    """
    ndim = len(grids)
    grids = np.asarray(grids)
    ns = grids[:,2]
    axes = [np.linspace(*grid) for grid in grids]
    if hnode is None:
        x = np.empty(list(ns)+[ndim])
        for index in np.ndindex(*ns):
            x[index+(None,)] = [axes[i][index[i]] for i in xrange(ndim)]
        return np.atleast_2d(x.squeeze()), ns            
    else:
        for index in np.ndindex(*ns):
            hnode[index] = [axes[i][index[i]] for i in xrange(ndim)]
        return ns

    def __init__(self):

        n = 3 # colors
        m = len(kernels2)
        (h, w) = kernels2[0].shape
        X = numpy.zeros(( (2**m)*(2**n), w*h*3),dtype='float32')
        idx = 0
        for i in numpy.ndindex(*([2]*m)):
            for j in numpy.ndindex(*([2]*n)):

                example = numpy.zeros((n,h,w), dtype='float32')
                for k in xrange(m):
                    if i[k]:
                        example[0,:,:] += j[0] * kernels2[k]
                        example[1,:,:] += j[1] * kernels2[k]
                        example[2,:,:] += j[2] * kernels2[k]

                X[idx,:] = example.reshape(h * w * n)
                idx += 1

        view_converter = dense_design_matrix.DefaultViewConverter((h,w,1))

        super(SuperImposedShapes,self).__init__(X = X, view_converter = view_converter)

        assert not numpy.any(numpy.isnan(self.X))

def process(destinations_list):
    # out
    out = [idx for idx in np.ndindex(destinations_list[0].arr.shape)]
    # generate slices list of dictionaries
    if len(destinations_list[-1].arr.shape) > 1:
        # multidimensional scan
        slices = []
        for i, idx in enumerate(np.ndindex(destinations_list[-1].arr.shape[:-1])):
            s = {}
            s['index'] = destinations_list[-1].arr.shape[-1]*i
            s['name'] = destinations_list[-1].hardware.name
            s['units'] = destinations_list[-1].units
            s['points'] = destinations_list[-1].arr[idx]
            if destinations_list[-1].method == 'set_position':
                s['use actual'] = True
            else:
                s['use actual'] = False
            slices.append(s)
    else:
        # 1D scan
        s = {}
        s['index'] = 0
        s['name'] = destinations_list[0].hardware.name
        s['units'] = destinations_list[0].units
        s['points'] = destinations_list[0].arr
        if destinations_list[0].method == 'set_position':
            s['use actual'] = True
        else:
            s['use actual'] = False
        slices = [s]
    return out, slices

def make_reg_masks(regfile,shape):
    r = pyregion.open(regfile)

    if len(r) != 2:
        raise Exception('Exactly two box regions required')

    paths = []
    for reg in r:
        paths.append(get_region_box(reg.coord_list))

    #Always have A be the top half
    if paths[0].get_extents().ymax > paths[1].get_extents().ymax:
        pathA = paths[0]
        pathB = paths[1]
    else:
        pathA = paths[1]
        pathB = paths[2]


    print 'Building skymasks'
    maskA = np.array([True if pathA.contains_point([x,y]) else False for x,y in np.ndindex(shape)])
    maskA = maskA.reshape(shape).T
    
    maskB = np.array([True if pathB.contains_point([x,y]) else False for x,y in np.ndindex(shape)])
    maskB = maskB.reshape(shape).T

    return (~maskA, ~maskB)

    def backprop(self, targets, mu):
        if len(targets) != self.nout:
            raise ValueError('wrong number of targets')

        # output deltas first
        # partial E wrt v_k = sum w_jk z_j where a_k = sigmoid(v_k)
        # odelta = (self.aout - targets) * dsigmoid(self.aout)
        # hdelta = np.dot(odelta, self.wout)[:-1] * dsigmoid(self.ain[:-1])
        
        # matrix ops not working for some reason :(, I have a bug
        # time to be more straightforward

        odelta = np.zeros(self.nout)
        for k in range(self.nout):
            odelta[k] = (self.aout[k] - targets[k]) * dsigmoid(self.aout[k])

        for k, j in np.ndindex(self.wout.shape):
            if self.nhid:
                self.wout[k, j] -= mu * odelta[k] * self.ahid[j]
            else:
                self.wout[k, j] -= mu * odelta[k] * self.ain[j]


        if self.nhid:
            hdelta = np.zeros(self.nhid)
            for j in range(self.nhid):
                hdelta[j] = dsigmoid(self.ahid[j]) * np.dot(self.wout[:, j], odelta)

            for j, i in np.ndindex(self.win.shape):
                self.win[j, i] -= mu * hdelta[j] * self.ain[i]

        # self.wout -= mu * np.outer(odelta, self.aout)
        # self.win -= mu * np.outer(hdelta, self.ain)

        return 0.5 * np.linalg.norm(targets - self.aout)

def procedure(ticks):
    n = 500
    b = .000006662 
    D = 1
    alpha = 2

    n_types = ticks**D
    #print 'Number of types: {}'.format(n_types)
    M = np.zeros([ticks**D, ticks**D])
    registry = {}
    
    next_id = 0

    for index in np.ndindex(tuple([ticks] * D)):
        i = index[:D]
        registry[i] = next_id 
        next_id += 1

    for index in np.ndindex(tuple([ticks]* D * 2)):
        i = index[:D]
        j = index[D:]

        if i != j:
            pos_i = [float(_i) / (ticks - 1) for _i in i]
            pos_j = [float(_j) / (ticks - 1) for _j in j]

            M[registry[i], registry[j]] = .5 * n**2 / n_types**2 *\
                b / (b + model.distance(None, pos_i, pos_j)**alpha) 

    eigvals = scipy.linalg.eigvals(M) 
    return max(eigvals)

def set_permutation_symmetry_fc3_deprecated(fc3):
    fc3_sym = np.zeros(fc3.shape, dtype='double')
    for (i, j, k) in list(np.ndindex(fc3.shape[:3])):
        fc3_sym[i, j, k] = set_permutation_symmetry_fc3_elem(fc3, i, j, k)

    for (i, j, k) in list(np.ndindex(fc3.shape[:3])):
        fc3[i, j, k] = fc3_sym[i, j, k]

    def run(self, triplet, is_sym_fc3_q=False):
        num_patom = self._primitive.get_number_of_atoms()
        if is_sym_fc3_q:
            index_exchage = np.array([[0,1,2],[1,2,0],[2,0,1],[0,2,1],[1,0,2],[2,1,0]])
            fc3_reciprocal = np.zeros(
                (num_patom, num_patom, num_patom, 3, 3, 3), dtype='complex128')
            for e, index in enumerate(index_exchage):
                self._triplet = triplet[index]
                self._fc3_reciprocal = np.zeros(
                    (num_patom, num_patom, num_patom, 3, 3, 3), dtype='complex128')
                self._real_to_reciprocal()
                for patoms in np.ndindex((num_patom, num_patom, num_patom)):
                    i,j,k = np.array(patoms)
                    ii, ji, ki = np.array(patoms)[index]
                    for cart in np.ndindex((3,3,3)):
                        l, m, n = np.array(cart)
                        li, mi, ni = np.array(cart)[index]
                        fc3_reciprocal[i,j,k,l,m,n] += self._fc3_reciprocal[ii, ji, ki, li, mi, ni] / 6
            self._fc3_reciprocal[:] = fc3_reciprocal

        else:
            self._triplet = triplet
            self._fc3_reciprocal = np.zeros(
                (num_patom, num_patom, num_patom, 3, 3, 3), dtype='complex128')
            self._real_to_reciprocal()

def confidence_scores(raw_counts, perm_counts, num_features):
    """Return confidence scores.
    
    """
    logging.debug(("Getting confidence scores for shape {shape} with "
                   "{num_features} features").format(
            shape=np.shape(raw_counts),
            num_features=num_features))
    if np.shape(raw_counts) != np.shape(perm_counts):
        raise Exception((
                "raw_counts and perm_counts must have same shape. "
                "raw_counts is {raw} and perm_counts is {perm}").format(
                raw=np.shape(raw_counts), perm=np.shape(perm_counts)))
    
    shape = np.shape(raw_counts)
    adjusted = np.zeros(shape)
    for idx in np.ndindex(shape[:-1]):
        adjusted[idx] = adjust_num_diff(perm_counts[idx], raw_counts[idx], num_features)

    # (unpermuted counts - mean permuted counts) / unpermuted counts
    res = (raw_counts - adjusted) / raw_counts

    for idx in np.ndindex(res.shape[:-1]):
        res[idx] = ensure_scores_increase(res[idx])

    return res

 def __init__(self, filename):
     print('Loading \"%s\" ..' % filename)
     with open(filename) as f:
         if f.readline().strip() != 'OFF': raise Exception("Invalid format")
         self.nverts, self.nfaces, _ = map(int, f.readline().split())
         self.vertices, self.faces = np.zeros((self.nverts, 3)), np.zeros((self.nfaces, 3), np.uint32)
         for i in range(self.nverts):
             self.vertices[i, :] = np.fromstring(f.readline(), sep=' ')
         for i in range(self.nfaces):
             self.faces[i, :] = np.fromstring(f.readline(), sep=' ', dtype=np.uint32)[1:]
     print('Computing face and vertex normals ..')
     v = [self.vertices[self.faces[:, i], :] for i in range(3)]
     face_normals = np.cross(v[2] - v[0], v[1] - v[0])
     face_normals /= np.linalg.norm(face_normals, axis=1)[:, None]
     self.normals = np.zeros((self.nverts, 3))
     for i, j in np.ndindex(self.faces.shape):
         self.normals[self.faces[i, j], :] += face_normals[i, :]
     self.normals /= np.linalg.norm(self.normals, axis=1)[:, None]
     print('Building adjacency matrix ..')
     self.adjacency = [set() for _ in range(self.nfaces)]
     for i, j in np.ndindex(self.faces.shape):
         e0, e1 = self.faces[i, j], self.faces[i, (j+1)%3]
         self.adjacency[e0].add(e1)
         self.adjacency[e1].add(e0)
     print('Randomly initializing fields ..')
     self.o_field = np.zeros((self.nverts, 3))
     self.p_field = np.zeros((self.nverts, 3))
     min_pos, max_pos = self.vertices.min(axis=0), self.vertices.max(axis=0)
     np.random.seed(0)
     for i in range(self.nverts):
         d, p = np.random.standard_normal(3), np.random.random(3)
         d -= np.dot(d, self.normals[i]) * self.normals[i]
         self.o_field[i] = d / np.linalg.norm(d)
         self.p_field[i] = (1-p) * min_pos + p * max_pos

    def numpy_max_pool_nd(input, ds, ignore_border=False, mode='max'):
        '''Helper function, implementing pool_nd in pure numpy'''
        if len(input.shape) < len(ds):
            raise NotImplementedError('input should have at least %s dim,'
                                      ' shape is %s'
                                      % (str(ds), str(input.shape)))
        nd = len(ds)
        si = [0] * nd
        if not ignore_border:
            for i in range(nd):
                if input.shape[-nd + i] % ds[i]:
                    si[i] += 1
        out_shp = list(input.shape[:-nd])
        for i in range(nd):
            out_shp.append(input.shape[-nd + i] // ds[i] + si[i])
        output_val = numpy.zeros(out_shp)
        func = numpy.max
        if mode == 'sum':
            func = numpy.sum
        elif mode != 'max':
            func = numpy.average

        for l in numpy.ndindex(*input.shape[:-nd]):
            for r in numpy.ndindex(*output_val.shape[-nd:]):
                patch = input[l][tuple(slice(r[i] * ds[i], (r[i] + 1) * ds[i])
                                       for i in range(nd))]
                output_val[l][r] = func(patch)
        return output_val

    def _make_default_data(geom, shape_np, dtype):
        # Check whether corners of each image plane are valid
        coords = []
        if not geom.is_regular:
            for idx in np.ndindex(geom.shape):
                pix = (np.array([0.0, float(geom.npix[0][idx] - 1)]),
                       np.array([0.0, float(geom.npix[1][idx] - 1)]))
                pix += tuple([np.array(2 * [t]) for t in idx])
                coords += geom.pix_to_coord(pix)

        else:
            pix = (np.array([0.0, float(geom.npix[0] - 1)]),
                   np.array([0.0, float(geom.npix[1] - 1)]))
            pix += tuple([np.array(2 * [0.0]) for i in range(geom.ndim - 2)])
            coords += geom.pix_to_coord(pix)

        if np.all(np.isfinite(np.vstack(coords))):
            if geom.is_regular:
                data = np.zeros(shape_np, dtype=dtype)
            else:
                data = np.full(shape_np, np.nan, dtype=dtype)
                for idx in np.ndindex(geom.shape):
                    data[idx,
                         slice(geom.npix[0][idx]),
                         slice(geom.npix[1][idx])] = 0.0
        else:
            data = np.full(shape_np, np.nan, dtype=dtype)
            idx = geom.get_idx()
            m = np.all(np.stack([t != -1 for t in idx]), axis=0)
            data[m] = 0.0

        return data

    def _nppma_ndindex(self, tup, N=None, ecut=None):
        """Generates vibrational signatures in NP per mode Aproximation (NPPMA) 
        
        """
        if N is None:
            raise Exception("Number of states to be used must be defined")

        if ecut is not None:
            vec = self.convert_energy_2_internal_u(ecut)

        two = numpy.zeros(len(tup), dtype=numpy.int)
        two[:] = N + 1

        for i in range(len(two)):
            if two[i] > tup[i] + 1:
                two[i] = tup[i] + 1
                
        shp = tuple(two)
        
        if ecut is not None:
            #with energy_units("int"):    
            for sig in numpy.ndindex(shp):
                en = self.convert_energy_2_internal_u(self.vibenergy(vsig=sig))
                if en <= vec:
                    yield sig
                        
        else:
            for sig in numpy.ndindex(shp):
                yield sig        

 def _iter_indexes(self, array):
     if self.select is None:
         for indexes in np.ndindex(array.shape[1:]):
             yield indexes
     else:
         for i0 in self.select:
             for irest in np.ndindex(array.shape[2:]):
                 yield (i0,) + irest

def python_local_maxima(data, wsize, mode=wrap):
  result = np.ones(shape=data.shape,dtype='bool')
  for pos in np.ndindex(data.shape):
    myval = data[pos]  
    for offset in np.ndindex(wsize):
      neighbor_idx = tuple(mode(p, o-w/2, w) for (p, o, w) in zip(pos, offset, wsize))
      result[pos] &= (data[neighbor_idx] <= myval)
  return result 

 def _geo_ts_to_vec(self, data, pts):
     res = {}
     for name, ts in data.items():
         tpe = self.source_type_map[name] 
         ids = [idx for idx in np.ndindex(pts.shape[:-1])]
         res[name] = tpe.vector_t([tpe(api.GeoPoint(*pts[idx]),
                                   ts[idx]) for idx in np.ndindex(pts.shape[:-1])])
     return res

def read_css():

    """I may never have written a more painful function in my life.

    If you want data, and are thinking of using numpy or pandas --
    read it in by hand.
    """
    # Can't delete array elements.
    #ix = 1 #so as to skip the header row.
    n = 21
    chunk = 2*n + 3 # Both matrices, and 3 extra rows
    advice = np.zeros((n, n, n))
    friendship = np.zeros((n, n, n))

    pdf = pd.read_csv("/Users/alexloewi/Documents/Data/cognitive social structures/rearranged_cogsocstr.txt", sep="\t")

    matrix_columns = pdf[pdf.columns[0:21]]
    #print 'matrix columns!!!!!!,', matrix_columns

    for perceiver in range(n):
        # This gets all the data for one person
        x = (chunk*perceiver)
        y = (chunk*(perceiver+1))-1

        a = np.array(matrix_columns.ix[x:x+20])
        np.fill_diagonal(a, 0)
        f = np.array(matrix_columns.ix[x+21:x+41])
        np.fill_diagonal(f, 0)
        advice[perceiver,:,:]     = a #np.array(matrix_columns.ix[0:20])
        friendship[perceiver,:,:] = f #np.array(matrix_columns.ix[21:41])

    # Consensus matrices (AND rule)
    ca = np.zeros((n,n))
    cf = np.zeros((n,n))

    for i,j in np.ndindex(ca.shape):
        if advice[i,i,j] + advice[j,i,j] == 2:
            ca[i,j] = 1

    for i,j in np.ndindex(cf.shape):
        if friendship[i,i,j] + friendship[j,i,j] == 2:
            cf[i,j] = 1

    # Self-proclaimed relationships (OR rule)
    sa = np.zeros((n,n))
    sf = np.zeros((n,n))

    for i,j in np.ndindex(sa.shape):
        if advice[i,i,j] + advice[j,i,j] >= 1:
            sa[i,j] = 1

    for i,j in np.ndindex(sf.shape):
        if friendship[i,i,j] + friendship[j,i,j] >= 1:
            sf[i,j] = 1


    return advice, friendship, ca, cf, sa, sf