def run(self, N=100):
        r'''
        '''
        im = self.image
        # Create a list of N random points to use as box centers
        pad = [0.1,0.1,0.45]  # Ensure points are near middle
        Cx = sp.random.randint(pad[0]*sp.shape(im)[0],(1-pad[0])*sp.shape(im)[0],N)
        Cy = sp.random.randint(pad[1]*sp.shape(im)[1],(1-pad[1])*sp.shape(im)[1],N)
        Cz = sp.random.randint(pad[2]*sp.shape(im)[2],(1-pad[2])*sp.shape(im)[2],N)
        C = sp.vstack((Cx,Cy,Cz)).T

        # Find maximum radius allowable for each point
        Rmax = sp.array(C>sp.array(sp.shape(im))/2)
        Rlim = sp.zeros(sp.shape(Rmax))
        Rlim[Rmax[:,0],0] = sp.shape(im)[0]
        Rlim[Rmax[:,1],1] = sp.shape(im)[1]
        Rlim[Rmax[:,2],2] = sp.shape(im)[2]
        R = sp.absolute(C-Rlim)
        R = R.astype(sp.int_)
        Rmin = sp.amin(R,axis=1)

        vol = []
        size = []
        porosity = []
        for i in range(0,N):
            for r in sp.arange(Rmin[i],1,-10):
                imtemp = im[C[i,0]-150:C[i,0]+150,C[i,1]-150:C[i,1]+150:,C[i,2]-r:C[i,2]+r]
                vol.append(sp.size(imtemp))
                size.append(2*r)
                porosity.append(sp.sum(imtemp==1)/(sp.size(imtemp)))

        vals = namedtuple('REV', ('porosity', 'size'))
        vals.porosity = porosity
        vals.size = size
        return vals

 def test_check_network_health_bidirectional_throat(self):
     net = OpenPNM.Network.Cubic(shape=[5, 5, 5])
     P12 = net['throat.conns'][0]
     net['throat.conns'][0] = [P12[1], P12[0]]
     a = net.check_network_health()
     assert sp.size(a['bidirectional_throats']) == 1
     assert sp.size(a['duplicate_throats']) == 0

    def __init__(self, network=None, phase=None, geometry=None,
                 pores=[], throats=[], **kwargs):
        super().__init__(**kwargs)
        logger.name = self.name

        # Associate with Network
        if network is None:
            self._net = GenericNetwork()
        else:
            self._net = network  # Attach network to self
            self._net._physics.append(self)  # Register self with network

        # Associate with Phase
        if phase is None:
            phase = GenericPhase(network=self._net)
        phase._physics.append(self)  # Register self with phase
        self._phases.append(phase)  # Register phase with self

        if geometry is not None:
            if (sp.size(pores) > 0) or (sp.size(throats) > 0):
                raise Exception('Cannot specify a Geometry AND pores or throats')
            pores = self._net.toindices(self._net['pore.' + geometry.name])
            throats = self._net.toindices(self._net['throat.' + geometry.name])

        # Initialize a label dictionary in the associated phase and network
        self._phases[0]['pore.'+self.name] = False
        self._phases[0]['throat.'+self.name] = False
        self._net['pore.'+self.name] = False
        self._net['throat.'+self.name] = False
        try:
            self.set_locations(pores=pores, throats=throats)
        except:
            self.controller.purge_object(self)
            raise Exception('Provided locations are in use, instantiation cancelled')

def MNEfit(stim,resp,order):
    # in order for dlogloss to work, we need to know -<g(yt(n),xt)>data
    # == calculate the constrained averages over the data set
    Nsamples = sp.size(stim,0)
    Ndim = sp.size(stim,1)
    psp = sp.mean(sp.mean(resp)) #spike probability (first constraint)
    avg = (1.0*stim.T*resp)/(Nsamples*1.0)
    avgs = sp.vstack((psp,avg))
    if(order > 1):
        avgsqrd = (stim.T*1.0)*(sp.array(sp.tile(resp,(1,Ndim)))*sp.array(stim))/(Nsamples*1.0)
        avgsqrd = sp.reshape(avgsqrd,(Ndim**2,1))
        avgs = sp.vstack((avgs,avgsqrd))
    
    #initialize params:
    pstart = sp.log(1/avgs[0,0] - 1)
    pstart = sp.hstack((pstart,(.001*(2*sp.random.rand(Ndim)-1))))
    if(order > 1):
        temp = .0005*(2*sp.random.rand(Ndim,Ndim)-1)
        pstart = sp.hstack((pstart,sp.reshape(temp+temp.T,(1,Ndim**2))[0]))
    
    #redefine functions with fixed vals:
    def logLoss(p):
        return LLF.log_loss(p, stim, resp, order)
    def dlogLoss(p):
        return LLF.d_log_loss(p, stim, avgs, order)
    #run the function:
    #pfinal = opt.fmin_tnc(logLoss,pstart,fprime=dlogLoss)
    # conjugate-gradient:
    pfinal = opt.fmin_cg(logLoss,pstart,fprime=dlogLoss)
    #pfinal = opt.fmin(logLoss,pstart,fprime=dlogLoss)
    return pfinal

def cdf2d(x, y, z):
    """
    2D Cumulative Density Function extracted from cdf2d.m
    Assumes the grid is uniformly defined, i.e. dx and dy are
    constant.
    """
    if size(x) < 2 or size(y) < 2:
        logger.critical("X or Y grids are not arrays")
        raise TypeError, "X or Y grids are not arrays"
    grid_area = abs(x[1] - x[0])*abs(y[1] - y[0])
    grid_volume = grid_area*z

    cz = zeros([x.size, y.size], 'd')
    cz[:, 0] = (grid_volume[:,0]).cumsum()
    cz[0, :] = (grid_volume[0,:]).cumsum()

    for i in xrange(1, len(x)):
        for j in xrange(1, len(y)):
            cz[i,j] = cz[i-1,j] + cz[i,j-1] - cz[i-1,j-1] + grid_volume[i,j]

    if cz[-1,-1] == 0:
        return cz
    #normalize cy
    cz = cz/cz[-1,-1]
    return cz

def jplot2d(*args):
	arg=' -2D -l SOUTH '
	i=0
	while (i<len(args)):
		i_arg=0
		q=False
		arg_txt=True
		while ((len(args)>(i_arg+i+1)) & arg_txt):
			if (type(args[i+i_arg+1]) is str):
				i_arg+=1
				arg=arg+' '+args[i+i_arg]
			else:
				arg_txt=False
				
		name='.'+str(scipy.stats.random_integers(10000))+'.plot.tmp'
		
		if ((len(args[i])==scipy.size(args[i]))):
			a=scipy.dstack((scipy.arange(1,scipy.size(args[i])+1),args[i]))[0]
			scipy.io.write_array(name,a)
		else:
			scipy.io.write_array(name,args[i])
			
		i+=i_arg+1
		arg=arg+' '+name
	
	os.system('java -cp jmathplot.jar org.math.plot.PlotPanel '+arg)

    def spheres(shape, radius, porosity):
        r"""
        Generate a packing of overlapping mono-disperse spheres

        Parameters
        ----------
        shape : list
            The size of the image to generate in [Nx, Ny, Nz] where N is the
            number of voxels.

        radius : scalar
            The radius of spheres in the packing

        porosity : scalar
            The porosity of the final image.  This number is approximated by
            the method so the returned result may not have exactly the
            specified value.

        Notes
        -----
        This method can also be used to generate a dispersion of hollows by
        treating ``porosity`` as solid volume fraction and inverting the
        returned image.
        """
        if sp.size(shape) == 1:
            shape = sp.full((3, ), int(shape))
        im = sp.zeros(shape, dtype=bool)
        while sp.sum(im)/sp.size(im) < (1 - porosity):
            temp = sp.rand(shape[0], shape[1], shape[2]) < 0.9995
            im = im + (spim.distance_transform_edt(temp) < radius)
        return ~im

    def blobs(shape, porosity, blobiness=8):
        """
        Generates an image containing amorphous blobs

        Parameters
        ----------
        shape : list
            The size of the image to generate in [Nx, Ny, Nz] where N is the
            number of voxels

        blobiness : scalar
            Controls the morphology of the image.  A higher number results in
            a larger number of smaller blobs.

        porosity : scalar
            The porosity of the final image.  This number is approximated by
            the method so the returned result may not have exactly the
            specified value.

        """
        if sp.size(shape) == 1:
            shape = sp.full((3, ), int(shape))
        [Nx, Ny, Nz] = shape
        sigma = sp.mean(shape)/(4*blobiness)
        mask = sp.rand(Nx, Ny, Nz)
        mask = spim.gaussian_filter(mask, sigma=sigma)
        hist = sp.histogram(mask, bins=1000)
        cdf = sp.cumsum(hist[0])/sp.size(mask)
        xN = sp.where(cdf >= porosity)[0][0]
        im = mask <= hist[1][xN]
        return im

 def test_far_apart_clusters_estimate_all(self):
     cluster1 = sp.randn(40,1000)
     cluster2 = sp.randn(40,1000) * 2
     cluster2[0,:] += 10
     clusterList1 = [cluster1[:,i]
                     for i in xrange(sp.size(cluster1,1))]
     clusterList2 = [cluster2[:,i]
                     for i in xrange(sp.size(cluster2,1))]
     total, pair = qa.overlap_fp_fn(
         {1: clusterList1, 2: clusterList2})
     self.assertLess(total[1][0], 1e-4)
     self.assertLess(total[1][1], 1e-4)
     self.assertLess(total[2][0], 1e-4)
     self.assertLess(total[2][1], 1e-4)
     self.assertLess(pair[1][2][0], 1e-4)
     self.assertLess(pair[1][2][1], 1e-4)
     self.assertLess(pair[2][1][0], 1e-4)
     self.assertLess(pair[2][1][1], 1e-4)
     self.assertGreater(total[1][0], 0.0)
     self.assertGreater(total[1][1], 0.0)
     self.assertGreater(total[2][0], 0.0)
     self.assertGreater(total[2][1], 0.0)
     self.assertGreater(pair[1][2][0], 0.0)
     self.assertGreater(pair[1][2][1], 0.0)
     self.assertGreater(pair[2][1][0], 0.0)
     self.assertGreater(pair[2][1][1], 0.0)

    def test_set_boundary_conditions_bctypes(self):
        self.alg.setup(invading_phase=self.water,
                       defending_phase=self.air,
                       trapping=True)
        Ps = sp.random.randint(0, self.net.Np, 10)

        self.alg.set_boundary_conditions(pores=Ps, bc_type='inlets')
        assert sp.sum(self.alg['pore.inlets']) == sp.size(sp.unique(Ps))
        self.alg['pore.inlets'] = False

        self.alg.set_boundary_conditions(pores=Ps, bc_type='outlets')
        assert sp.sum(self.alg['pore.outlets']) == sp.size(sp.unique(Ps))
        self.alg['pore.outlets'] = False

        self.alg.set_boundary_conditions(pores=Ps, bc_type='residual')
        assert sp.sum(self.alg['pore.residual']) == sp.size(sp.unique(Ps))
        self.alg['pore.residual'] = False

        flag = False
        try:
            self.alg.set_boundary_conditions(pores=Ps, bc_type='bad_type')
        except:
            flag = True
        assert flag

        flag = False
        try:
            self.alg.set_boundary_conditions(bc_type=None, mode='bad_type')
        except:
            flag = True
        assert flag

def ideal_data(num, dimU, dimY, dimX, noise=1):
    """Linear system data"""
    # generate randomized linear system matrices
    A = randn(dimX, dimX)
    B = randn(dimX, dimU)
    C = randn(dimY, dimX)
    D = randn(dimY, dimU)

    # make sure state evolution is stable
    U, S, V = svd(A)
    A = dot(U, dot(diag(S / max(S)), V))
    U, S, V = svd(B)
    S2 = zeros((size(U,1), size(V,0)))
    S2[:,:size(U,1)] = diag(S / max(S))
    B = dot(U, dot(S2, V))

    # random input
    U = randn(num, dimU)

    # initial state
    X = reshape(randn(dimX), (1,-1))

    # initial output
    Y = reshape(dot(C, X[-1]) + dot(D, U[0]), (1,-1))

    # generate next state
    X = concatenate((X, reshape(dot(A, X[-1]) + dot(B, U[0]), (1,-1))))

    # and so forth
    for u in U[1:]:
        Y = concatenate((Y, reshape(dot(C, X[-1]) + dot(D, u), (1,-1))))
        X = concatenate((X, reshape(dot(A, X[-1]) + dot(B, u), (1,-1))))

    return U, Y + randn(num, dimY) * noise

  def convertToTimeAndFrequencyData(self, grain, target):
    d = self.sample_duration
    length = max(sp.shape(sp.arange(1, sp.size(target) - sp.size(self.gabor[1]), d)))
    scale = sp.zeros((88, length))
    datacapsule = sp.zeros((sp.shape(self.gabor)[1], grain))

    # 行列を束ねて処理
    #   個々にgabor*datadataを計算すると時間がかかる
    #   一気にdatadataを束ねるとメモリ消費量が半端ない
    #   よってdatadataに定義された数だけ束ねて処理
    m = 0
    datasize = sp.size(target) - sp.size(self.gabor[1])

    for k in sp.arange(1, datasize+1, d*grain):
      capsule_pointer = 0
      endl = k+d*grain

      if endl > datasize:
        endl = k + datasize%(d*grain)

      for l in sp.arange(k, endl, d):
        datadata = target[l:l+sp.size(self.gabor[1])]
        datacapsule[:, capsule_pointer] = datadata
        capsule_pointer += 1

      try:
        scale[:, m:m+grain] = sp.absolute(
            sp.dot(self.gabor,datacapsule[:, :capsule_pointer]))
      except ValueError:
        pass
      m += grain
    
    self.time_freq = scale

    def find_clusters(self, mask=[]):
        r"""
        Identify connected clusters of pores in the network.

        Parameters
        ----------
        mask : array_like, boolean
            A list of active nodes.  This method will automatically search
            for clusters based on site or bond connectivity depending on
            wheather the received mask is Np or Nt long.

        Returns
        -------
        clusters : array_like
            An Np long list of clusters numbers

        """
        if sp.size(mask) == self.num_throats():
            # Convert to boolean mask if not already
            temp = sp.zeros((self.num_throats(),), dtype=bool)
            temp[mask] = True
        elif sp.size(mask) == self.num_pores():
            conns = self.find_connected_pores(throats=self.throats())
            conns[:, 0] = mask[conns[:, 0]]
            conns[:, 1] = mask[conns[:, 1]]
            temp = sp.array(conns[:, 0]*conns[:, 1], dtype=bool)
        else:
            raise Exception('Mask received was neither Nt nor Np long')
        temp = self.create_adjacency_matrix(data=temp,
                                            sprsfmt='csr',
                                            dropzeros=True)
        clusters = sprs.csgraph.connected_components(csgraph=temp,
                                                     directed=False)[1]
        return clusters

    def __init__(self,filename='', path='', xtra_pore_data=None, xtra_throat_data=None,**kwargs):

        r"""
        """
        super(MatFile,self).__init__(**kwargs)
        if filename == '':
            return
        if path == '':
            path = os.path.abspath('.')
        self._path = path
        filepath = os.path.join(self._path,filename)
        self._xtra_pore_data=xtra_pore_data
        self._xtra_throat_data=xtra_throat_data
        self._dictionary=spio.loadmat(filepath)

        self._Np=sp.size(self._dictionary['pnumbering'])
        self._Nt=sp.size(self._dictionary['tnumbering'])

        #Run through generation steps
        self._add_pores()
        self._add_throats()
        self._remove_disconnected_clusters()
        self._add_xtra_pore_data()
        self._add_xtra_throat_data()
        self._add_geometry()

def plotConn():
    # Create plot
    figh = figure(figsize=(8,6))
    figh.subplots_adjust(left=0.02) # Less space on left
    figh.subplots_adjust(right=0.98) # Less space on right
    figh.subplots_adjust(top=0.96) # Less space on bottom
    figh.subplots_adjust(bottom=0.02) # Less space on bottom
    figh.subplots_adjust(wspace=0) # More space between
    figh.subplots_adjust(hspace=0) # More space between
    h = axes()
    totalconns = zeros(shape(f.connprobs))
    for c1 in range(size(f.connprobs,0)):
        for c2 in range(size(f.connprobs,1)):
            for w in range(f.nreceptors):
                totalconns[c1,c2] += f.connprobs[c1,c2]*f.connweights[c1,c2,w]*(-1 if w>=2 else 1)
    imshow(totalconns,interpolation='nearest',cmap=bicolormap(gap=0))

    # Plot grid lines
    hold(True)
    for pop in range(f.npops):
        plot(array([0,f.npops])-0.5,array([pop,pop])-0.5,'-',c=(0.7,0.7,0.7))
        plot(array([pop,pop])-0.5,array([0,f.npops])-0.5,'-',c=(0.7,0.7,0.7))

    # Make pretty
    h.set_xticks(range(f.npops))
    h.set_yticks(range(f.npops))
    h.set_xticklabels(f.popnames)
    h.set_yticklabels(f.popnames)
    h.xaxis.set_ticks_position('top')
    xlim(-0.5,f.npops-0.5)
    ylim(f.npops-0.5,-0.5)
    clim(-abs(totalconns).max(),abs(totalconns).max())
    colorbar()

 def __init__(self, start):
     """Initialize a Particle at the given start vector."""
     self.dim = scipy.size(start)
     self.position = start
     self.velocity = scipy.zeros(scipy.size(start))
     self.bestPosition = scipy.zeros(scipy.size(start))
     self._fitness = None
     self.bestFitness = -scipy.inf

def slepian_bw(windows,Fs):
    l = s.size(windows,axis = 0)
    N = s.size(windows,axis = 1)
    avg = s.zeros(2*N-1)
    for i in range(l):
        w_norm = windows[i,:]/s.sqrt(s.sum(windows[i,:]**2))
        avg += s.signal.correlate(w_norm,w_norm)**2        
    return Fs/s.sum(avg/l)

def maxCellNum(gridLimit, gridSpace):
    """
    Get maximum cell number based on grid and grid space
    """
    latCells = size(arange(gridLimit['yMax'], gridLimit['yMin'], -gridSpace['y']))
    lonCells = size(arange(gridLimit['xMin'], gridLimit['xMax'], gridSpace['x']))

    return latCells*lonCells - 1

def off(matrix):
    h_size = sp.size(matrix, 1) // 3
    v_size = sp.size(matrix, 0) // 4
    size = sp.size(matrix, 0) * sp.size(matrix, 1)
    T = sp.zeros(size)
    for i in range(size):
        T[i] = matrix[(i // 12) // h_size + (i % 12) // 3 * v_size, (i // 12) % h_size + i % 3 * h_size]
    return T

    def folderinfo(self, path):
        """
        accept folder path and return tuple:
            - target path
            - item count
            - list of items
        :param path:
        :return:
        """
        loc = self.pstruct[1]
        for p in path.split(":"): #drill down to the folder specified in the path...
            try:
                loc[p]
            except:
                loc=loc[p.encode()]
            else:
                loc=loc[p]
        items=[]
        for key in loc.keys():
            try:
                key.decode()
            except:
                items.append(key)
            else:
                items.append(key.decode())

        waves=[]
        variables=[]
        strings=[]
        folders=[]
        for item in items:
            try:
                loc[item]
            except:
                if type(loc[item.encode()])==igor.record.wave.WaveRecord:
                    waves.append(item)
                if type(loc[item.encode()])==int or type(loc[item.encode()])==float or type(loc[item.encode()])==numpy.float64:
                    variables.append(item)
                if type(loc[item.encode()])==str:
                    strings.append(item)
                if type(loc[item.encode()])==dict:
                    folders.append(item)
            else:
                if type(loc[item])==igor.record.wave.WaveRecord:
                    waves.append(item)
                if type(loc[item])==int or type(loc[item])==float:
                    variables.append(item)
                if type(loc[item])==str:
                    strings.append(item)
                if type(loc[item])==dict:
                    folders.append(item)

        # waves = [item for item in items if type(loc[item.encode()])==igor.record.wave.WaveRecord]
        # variables = [item for item in items if type(loc[item.encode()])==int or type(item)==float]
        # strings = [item for item in items if type(loc[item.encode()])==str]
        # folders = [item for item in items if type(loc[item.encode()])==dict]

        return dict(fullpath=path,waves=(scipy.size(waves), waves),variables=(scipy.size(variables),variables),strings=(scipy.size(strings),strings),folders=(scipy.size(folders),folders))