def _generate_pores(self):
        r"""
        Generate the pores (coordinates, numbering and types)
        """
        self._logger.info("generate_pores: Create specified number of pores")

        #Find non-zero elements in image
        template = self._template
        Np = np.sum(template > 0)
        #Add pores to data and ifo
        pind = np.arange(0, Np)
        self.set_pore_info(label='all', locations=pind)
        self.set_pore_data(prop='numbering', data=pind)  # Remove eventually

        
        img_ind = np.ravel_multi_index(sp.nonzero(template), dims=sp.shape(template), order='F')
        self.set_pore_data(prop='voxel_index', data=img_ind)

        #This voxel_to_pore map is messy but works
        temp = sp.prod(sp.shape(template))*sp.ones(np.prod(sp.shape(template),),dtype=sp.int32)
        temp[img_ind] = pind
        self._voxel_to_pore_map = temp

        coords = self._Lc*(0.5 + np.transpose(np.nonzero(template)))
        self.set_pore_data(prop='coords', data=coords)
        self._logger.debug("generate_pores: End of method")

def check_if_click_is_on_an_existing_point(mouse_x_coord,mouse_y_coord):
    # First, figure out how many points we have.
    # Each point is one row in the coords_array,
    # so we count the number of rows, which is dimension-0 for Python
    number_of_points = scipy.shape(coords_array)[0]    
    this_coord = scipy.array([[ mouse_x_coord, mouse_y_coord ]]) 
            # The double square brackets above give the this_coord array 
            # an explicit structure of having rows and also columns
    if number_of_points > 0:  
        # If there are some points, we want to calculate the distance
        # of the new mouse-click location from every existing point.
        # One way to do this is to make an array which is the same size
        # as coords_array, and which contains the mouse x,y-coords on every row.
        # Then we can subtract that xy_coord_matchng_matrix from coords_array
        ones_vec = scipy.ones((number_of_points,1))
        xy_coord_matching_matrix = scipy.dot(ones_vec,this_coord)
        distances_from_existing_points = (coords_array - xy_coord_matching_matrix)
        squared_distances_from_existing_points = distances_from_existing_points**2
        sum_sq_dists = scipy.sum(squared_distances_from_existing_points,axis=1) 
                   # The axis=1 means "sum over dimension 1", which is columns for Python          
        euclidean_dists = scipy.sqrt(sum_sq_dists)
        distance_threshold = 0.5
        within_threshold_points = scipy.nonzero(euclidean_dists < distance_threshold )
        num_within_threshold_points = scipy.shape(within_threshold_points)[1]
        if num_within_threshold_points > 0:
            # We only want one matching point.
            # It's possible that more than one might be within threshold.
            # So, we take the unique smallest distance
            point_to_be_deleted = scipy.argmin(euclidean_dists)
            return point_to_be_deleted
        else: # If there are zero points, then we are not deleting any 
            point_to_be_deleted = -1
            return point_to_be_deleted

def full_obs(sys,poles):
    """Full order observer of the system sys

    Call:
    obs=full_obs(sys,poles)

    Parameters
    ----------
    sys : System in State Space form
    poles: desired observer poles

    Returns
    -------
    obs: ss
    Observer

    """
    if isinstance(sys, TransferFunction):
        "System must be in state space form"
        return
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    poles=mat(poles)
    L=place(a.T,c.T,poles)
    L=mat(L).T
    Ao=a-L*c
    Bo=hstack((b-L*d,L))
    n=shape(Ao)
    m=shape(Bo)
    Co=eye(n[0],n[1])
    Do=zeros((n[0],m[1]))
    obs=ss(Ao,Bo,Co,Do,sys.Tsamp)
    return obs

    def	results(self):
	""" This method computes the results as a dictionary """

	result = {}

	# here you can either use a perturbed fcc lattice or completely random positions        
        #positions=self.makeRandomPositions()
        positions = self.makePerturbedLattice()

        gradE = self.gradE(positions)
        force = self.getForce(positions)

        result['gradE'] = gradE
        result['force'] = force

        # what should the error tolerance be on this?
        error = abs((gradE+force)/force)

	max_error = error.max()
	
	# this is to identify where the max error is
	index1 = error.argmax()/scipy.shape(error)[1]  
	index2 = error.argmax()-scipy.shape(error)[1]*(index1)

        if (error > 10**(-8)).any(): # expected error according to numerical recipes
                result['Equal'] = False
        else:
                result['Equal'] = True

        result['errors'] = error
	result['max error'] = max_error
	result['max error location']=scipy.array([index1,index2])

        return result

def kalman_filter(b,
                  V,
                  Phi,
                  y,
                  X,
                 sigma,
                  Sigma,
                  switch = 0,
                  D = None,
                  d = None,
                  G = None,
                  a = None,
                  c = None):
    r"""
    
    .. math::
       :nowrap:

       \begin{eqnarray*}
       \beta_{t|t-1} = \Phi \: \beta_{t-1|t-1}\\
       V_{t|t-1} = \Phi  V_{t-1|t-1} \Phi ^T + \Sigma \\
       e_t = y_t -  X_t \beta_{t|t-1}\\
       K_t =  V_{t|t-1} X_t^T (\sigma + X_t V_{t|t-1} X_t )^{-1}\\
       \beta_{t|t} = \beta_{t|t-1} + K_t e_t\\
       V_{t|t} = (I - K_t X_t^T) V_{t|t-1}\\
       \end{eqnarray*}

    """

    n = scipy.shape(X)[1]
    beta = scipy.empty(scipy.shape(X))
    n = len(b)
    if D is None:
        D = scipy.ones((1, n))
    if d is None:
        d = scipy.matrix(1.)
    if G is None:
        G = scipy.identity(n)
    if a is None:
        a = scipy.zeros((n, 1))
    if c is None:
        c = scipy.ones((n, 1))
#        import code; code.interact(local=locals())
    (b, V) = kalman_predict(b, V, Phi, Sigma)
    for i in xrange(len(X)):
        beta[i] = scipy.array(b).T
        (b, V, e, K) = kalman_upd(b,
                                V,
                                y[i],
                                X[i],
                                sigma,
                                Sigma,
                                switch,
                                D,
                                d,
                                G,
                                a,
                                c)
        (b, V) = kalman_predict(b, V, Phi, Sigma)
    return beta

    def __init__(self, params={}, **kwargs):
        """Initialize the AVSwThetaFB instance.  params should have
        the following keys:

        'ks'   - spring constant (required)
        'c'    - damper coeffiecent (defaults to 0)
        'Ka'   - actuator gain (defaults to 1)
        'Gc'   - proportional feedback gain (defaults to 1)
                 (Gc can be a transfer function modelled as a
                 ratio of polynomials when passed to FORTRAN.)
        'axis' - axis about which the actuator rotates (defaults to 1)
        'tau'  - first order pole of actuator."""
        if not params.has_key("Ka"):
            params["Ka"] = 1
        if not params.has_key("axis"):
            params["axis"] = 1
        if not params.has_key("c"):
            params["c"] = 0
        if not params.has_key("Gc"):
            params["Gc"] = 1
        if shape(params["Ka"]):
            params["Ka"] = params["Ka"][0]
        if params.has_key("tau"):
            if shape(params["tau"]):
                params["tau"] = params["tau"][0]
        TMMElementIHT.__init__(self, "avsthfb", params, **kwargs)

  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 autolim(self, myattr, freqvect, margin=0.1,db=0):
     if self.freqlim:
         ind1=thresh(freqvect,self.freqlim[0])
         ind2=thresh(freqvect,self.freqlim[1])
     else:
         ind1=0
         ind2=-1
     mymatrix=getattr(self,myattr)
     if len(shape(mymatrix))==1:
         submat=mymatrix[ind1:ind2]
     else:
         mymatrix=colwise(mymatrix)
         submat=mymatrix[ind1:ind2,:]
     if db:
         submat=20*log10(submat)
     # max and min need to be done columnwise
     # (maybe a Krauss matrix max)
     if len(shape(submat))==2:
         mymax=[]
         mymin=[]
         for q in range(shape(submat)[1]):
             mymax.append(max(submat[:,q]))
             mymin.append(min(submat[:,q]))
     else:
         mymax=max(submat)
         mymin=min(submat)
     if len(shape(mymax))>0:
         mymax=max(mymax)
         mymin=min(mymin)
     myspan=mymax-mymin
     mymargin=margin*myspan
     limout=[mymin-mymargin, mymax+mymargin]
     setattr(self,myattr+"lim",limout)
     return limout

def minreal(sys):
    """Minimal representation for state space systems

    Usage
    =====
    [sysmin]=minreal[sys]

    Inputs
    ------

    sys: system in ss or tf form

    Outputs
    -------
    sysfin: system in state space form
    """
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    nx=shape(a)[0]
    ni=shape(b)[1]
    no=shape(c)[0]

    out=tb03ad(nx,no,ni,a,b,c,d,'R')

    nr=out[3]
    A=out[0][:nr,:nr]
    B=out[1][:nr,:ni]
    C=out[2][:no,:nr]
    sysf=ss(A,B,C,sys.D,sys.Tsamp)
    return sysf

def dsimul(sys,u):
    """Simulate the discrete system sys
    Only for discrete systems!!!

    Call:
    y=dsimul(sys,u)

    Parameters
    ----------
    sys : Discrete System in State Space form
    u   : input vector
    Returns
    -------
    y: ndarray
    Simulation results

    """
    a=mat(sys.A)
    b=mat(sys.B)
    c=mat(sys.C)
    d=mat(sys.D)
    nx=shape(a)[0]
    ns=shape(u)[1]
    xk=zeros((nx,1))
    for i in arange(0,ns):
        uk=u[:,i]
        xk_1=a*xk+b*uk
        yk=c*xk+d*uk
        xk=xk_1
        if i==0:
            y=yk
        else:
            y=hstack((y,yk))
    y=array(y).T
    return y

 def GetMat(self,s,sym=False):
     """Return the element transfer matrix for the
     TorsionalSpringDamper element.  If sym=True, 's' must be a
     symbolic string and a matrix of strings will be returned.
     Otherwise, 's' is a numeric value (probably complex) and the
     matrix returned will be complex."""
     N=self.maxsize
     if sym:
         myparams=self.symparams
     else:
         myparams=self.params
     k=myparams['k']
     c=myparams['c']
     springterm=1/(k[0]+c[0]*s)
     if sym:
         maxlen=len(springterm)+10
         matout=eye(N,dtype='f')
         matout=matout.astype('S%d'%maxlen)
     else:
         matout=eye(N,dtype='D')
     matout[1,2]=springterm
     if max(shape(k))>1 and self.maxsize>=8:
         matout[5,6]=1/(k[1]+c[1]*s)
     if max(shape(k))>2 and self.maxsize>=12:
         matout[9,10]=1/(k[2]+c[2]*s)
     return matout

def _update_network(network, net):
    # Infer Np and Nt from length of given prop arrays in file
    for element in ['pore', 'throat']:
        N = [_sp.shape(net[i])[0] for i in net.keys() if i.startswith(element)]
        if N:
            N = _sp.array(N)
            if _sp.all(N == N[0]):
                if (network._count(element) == N[0]) \
                        or (network._count(element) == 0):
                    network.update({element+'.all': _sp.ones((N[0],),
                                                             dtype=bool)})
                    net.pop(element+'.all', None)
                else:
                    raise Exception('Length of '+element+' data in file ' +
                                    'does not match network')
            else:
                raise Exception(element+' data in file have inconsistent ' +
                                'lengths')

    # Add data on dummy net to actual network
    for item in net.keys():
        # Try to infer array types and change if necessary
        # Chcek for booleans disguised and 1's and 0's
        num0s = _sp.sum(net[item] == 0)
        num1s = _sp.sum(net[item] == 1)
        if (num1s + num0s) == _sp.shape(net[item])[0]:
            net[item] = net[item].astype(bool)
        # Write data to network object
        if item not in network:
            network.update({item: net[item]})
        else:
            logger.warning('\''+item+'\' already present')
    return network

    def num_neighbors(self,pnums,labels=['all']):
        r"""
        Returns an ndarray containing the number of neigbhor pores for each 
        element in Pnums

        Parameters
        ----------
        pnums : array_like
            Pores whose neighbors are to be counted
        labels : list of string, optional
            The pore labels that should be included in the count

        Returns
        -------
        num_neighbors : 1D array with number of neighbors in each element, 
        useful for finding the number of neighbors of a certain type

        Examples
        --------
        >>> pn = OpenPNM.Network.TestNet()
        >>> Pnum = [0,1]
        >>> pn.num_neighbors(Pnum,flatten=False)
        array([3, 4], dtype=int8)
        >>> pn.num_neighbors(Pnum)
        7
        """
        #Convert string to list, if necessary
        if type(labels) == str: labels = [labels]

        #Count number of neighbors
        neighborPs = self.find_neighbor_pores(pnums,labels=labels,flatten=False)
        num = sp.zeros(sp.shape(neighborPs),dtype=sp.int8)
        for i in range(0,sp.shape(num)[0]):
            num[i] = sp.size(neighborPs[i])
        return num

 def getHDF5Description(self):
     u_shape = scipy.shape(self.u)
     lambd_shape = scipy.shape(self.lambd)
     class SystemSave(tables.IsDescription):
         u = tables.FloatCol(shape = u_shape)
         lambd = tables.FloatCol(shape = lambd_shape)
     return SystemSave

    def __init__(self, U, Y, statedim, reg=None):
        if size(shape(U)) == 1:
            U = reshape(U, (-1,1))
        if size(shape(Y)) == 1:
            Y = reshape(Y, (-1,1))
        if reg is None:
            reg = 0

        yDim = size(Y,1)
        uDim = size(U,1)

        self.output_size = size(Y,1) # placeholder

        # number of samples of past/future we'll mash together into a 'state'
        width = 1
        # total number of past/future pairings we get as a result
        K = size(U,0) - 2 * width + 1

        # build hankel matrices containing pasts and futures
        U_p = array([ravel(U[t : t + width]) for t in range(K)]).T
        U_f = array([ravel(U[t + width : t + 2 * width]) for t in range(K)]).T
        Y_p = array([ravel(Y[t : t + width]) for t in range(K)]).T
        Y_f = array([ravel(Y[t + width : t + 2 * width]) for t in range(K)]).T

        # solve the eigenvalue problem
        YfUfT = dot(Y_f, U_f.T)
        YfUpT = dot(Y_f, U_p.T)
        YfYpT = dot(Y_f, Y_p.T)
        UfUpT = dot(U_f, U_p.T)
        UfYpT = dot(U_f, Y_p.T)
        UpYpT = dot(U_p, Y_p.T)
        F = bmat([[None, YfUfT, YfUpT, YfYpT],
                  [YfUfT.T, None, UfUpT, UfYpT],
                  [YfUpT.T, UfUpT.T, None, UpYpT],
                  [YfYpT.T, UfYpT.T, UpYpT.T, None]])
        Ginv = bmat([[pinv(dot(Y_f,Y_f.T)), None, None, None],
                     [None, pinv(dot(U_f,U_f.T)), None, None],
                     [None, None, pinv(dot(U_p,U_p.T)), None],
                     [None, None, None, pinv(dot(Y_p,Y_p.T))]])
        F = F - eye(size(F, 0)) * reg

        # Take smallest eigenvalues
        _, W = eigs(Ginv.dot(F), k=statedim, which='SR')

        # State sequence is a weighted combination of the past
        W_U_p = W[ width * (yDim + uDim) : width * (yDim + uDim + uDim), :]
        W_Y_p = W[ width * (yDim + uDim + uDim):, :]
        X_hist = dot(W_U_p.T, U_p) + dot(W_Y_p.T, Y_p)

        # Regress; trim inputs to match the states we retrieved
        R = concatenate((X_hist[:, :-1], U[width:-width].T), 0)
        L = concatenate((X_hist[:, 1: ], Y[width:-width].T), 0)
        RRi = pinv(dot(R, R.T))
        RL  = dot(R, L.T)
        Sys = dot(RRi, RL).T
        self.A = Sys[:statedim, :statedim]
        self.B = Sys[:statedim, statedim:]
        self.C = Sys[statedim:, :statedim]
        self.D = Sys[statedim:, statedim:]

def nonna_lsq(target, aux, idx=(), names=(), order=2):
	"""
	This function returns the coefficients of the least square prediction of the target
	signal, using the auxiliary signals and their powers, as specified by the order argument.
	
	Input arguments:
	target = target signal
	aux    = matrix of auxiliary signals
	idx    = boolean vector to select a subset of the data for the LSQ fit
	order  = order of the polynomial of aux signals to be used in the fit, default is 2
	names  = list of the auxiliary signal names
	
	Output:
	p      = list of coefficients
	X      = matrix of the signals used in the reconstruction
	cnames = list of the corresponding signals
	
	Note that the mean will be removed from the auxiliary signals. 
	"""
	# number of auxiliary channels
	naux = scipy.shape(aux[1])
	
	if len(names) == 0:
		# since the user didn't provide signal names, let's build some
		names = map(lambda x: 'S'+str(x), scipy.arange(naux)+1)
		
	if len(idx) == 0:
		# no index means use all
		idx = numpy.array(target, dtype=bool)
		idx[:] = True
	
	##### PREPARE CHANNELS FOR LSQ PREDICTION 

	# prepare channels and their squared values
	X = scipy.zeros((scipy.shape(aux)[0], order*scipy.shape(aux)[1]+1))
	cnames = []
	for i in range(scipy.shape(aux)[1]):
		for j in range(order):
			# add the (j+1)th power of the signal after removing the mean
			X[:,order*i+j] = numpy.power((aux[:,i] - scipy.mean(aux[idx,i])), j+1)
			# then remove the mean of the result
			X[:,order*i+j] = X[:,order*i+j] - scipy.mean(X[idx,order*i+j])
			# save the name, including the power
			if j==0:
				cnames.append(names[i])
			else:
				cnames.append(names[i]+'^'+str(j+1))
				
	# add a constant at the end of the list
	X[:,-1] = 1
	cnames.append('1')
	# convert to matrix object for simpler manipulation
	X = scipy.mat(X)
	
	##### best estimate of coefficients to minimize the squared error
	p = scipy.linalg.inv(X[idx,:].T * X[idx,:]) * X[idx,:].T * scipy.mat(target[idx]).T

	# return all the results
	return p, X, cnames

    def __init__(self, x, y, axis=-1, makecopy=0, bounds_error=1,
                 fill_value=None):
        """Initialize a piecewise-constant interpolation class

        Description:
          x and y are arrays of values used to approximate some function f:
            y = f(x)
          This class returns a function whose call method uses piecewise-
          constant interpolation to find the value of new points.

        Inputs:
            x -- a 1d array of monotonically increasing real values.
                 x cannot include duplicate values. (otherwise f is
                 overspecified)
            y -- an nd array of real values.  y's length along the
                 interpolation axis must be equal to the length
                 of x.
            axis -- specifies the axis of y along which to 
                    interpolate. Interpolation defaults to the last
                    axis of y.  (default: -1)
            makecopy -- If 1, the class makes internal copies of x and y.
                    If 0, references to x and y are used. The default 
                    is to copy. (default: 0)
            bounds_error -- If 1, an error is thrown any time interpolation
                            is attempted on a value outside of the range
                            of x (where extrapolation is necessary).
                            If 0, out of bounds values are assigned the
                            NaN (#INF) value.  By default, an error is
                            raised, although this is prone to change.
                            (default: 1)
        """
        self.datapoints = (array(x, Float), array(y, Float))   # RHC -- for access from PyDSTool
        self.type = Float   # RHC -- for access from PyDSTool
        self.axis = axis
        self.makecopy = makecopy   # RHC -- renamed from copy to avoid nameclash
        self.bounds_error = bounds_error
        if fill_value is None:
            self.fill_value = NaN   # RHC -- was:   array(0.0) / array(0.0)
        else:
            self.fill_value = fill_value

        # Check that both x and y are at least 1 dimensional.
        if len(shape(x)) == 0 or len(shape(y)) == 0:
            raise ValueError, "x and y arrays must have at least one dimension."  
        # make a "view" of the y array that is rotated to the
        # interpolation axis.
        oriented_x = x
        oriented_y = swapaxes(y,self.interp_axis,axis)
        interp_axis = self.interp_axis
        len_x,len_y = shape(oriented_x)[interp_axis], \
                            shape(oriented_y)[interp_axis]
        if len_x != len_y:
            raise ValueError, "x and y arrays must be equal in length along "\
                              "interpolation axis."
        if len_x < 2 or len_y < 2:
            raise ValueError, "x and y arrays must have more than 1 entry"            
        self.x = array(oriented_x,copy=self.makecopy)
        self.y = array(oriented_y,copy=self.makecopy)

def augment_3_vector(v=array([0.0, 0.0, 0.0]), free=False):
    if free:
        freeval = 0.0
    else:
        freeval = 1.0
    assert shape(v) == (3,), "v argument must be 3-vector -- found " + str(shape(v))
    vaug = resize(v, (4,))
    vaug[3] = freeval
    return vaug

 def test_save_and_load_networkx_no_phases(self):
     G = op.io.NetworkX.to_networkx(network=self.net)
     project = op.io.NetworkX.from_networkx(G)
     assert len(project) == 1
     net = project.network
     assert net.Np == 8
     assert net.Nt == 12
     assert sp.shape(net['pore.coords']) == (8, 3)
     assert sp.shape(net['throat.conns']) == (12, 2)

 def test_save_and_load_csv_no_phases(self):
     fname = os.path.join(TEMP_DIR, 'test_save_csv_1')
     io.CSV.save(network=self.net, filename=fname)
     assert os.path.isfile(fname+'.csv')
     net = io.CSV.load(fname+'.csv')
     assert net.Np == 27
     assert net.Nt == 54
     assert sp.shape(net['pore.coords']) == (27, 3)
     assert sp.shape(net['throat.conns']) == (54, 2)