def add_boundary_pores(self, labels=['top', 'bottom', 'front', 'back',
                                         'left', 'right'], offset=None):
        r"""
        Add boundary pores to the specified faces of the network

        Pores are offset from the faces of the domain.

        Parameters
        ----------
        labels : string or list of strings
            The labels indicating the pores defining each face where boundary
            pores are to be added (e.g. 'left' or ['left', 'right'])

        offset : scalar or array_like
            The spacing of the network (e.g. [1, 1, 1]).  This must be given
            since it can be quite difficult to infer from the network,
            for instance if boundary pores have already added to other faces.

        """
        offset = sp.array(offset)
        if offset.size == 1:
            offset = sp.ones(3)*offset
        for item in labels:
            Ps = self.pores(item)
            coords = sp.absolute(self['pore.coords'][Ps])
            axis = sp.count_nonzero(sp.diff(coords, axis=0), axis=0) == 0
            ax_off = sp.array(axis, dtype=int)*offset
            if sp.amin(coords) == sp.amin(coords[:, sp.where(axis)[0]]):
                ax_off = -1*ax_off
            topotools.add_boundary_pores(network=self, pores=Ps, offset=ax_off,
                                         apply_label=item + '_boundary')

 def infer_diag_post(self,X_ii,D_i):
     
     X_i = dc(X_ii)
     ns = len(D_i)
     
     X_i.resize([ns,self.D])
     [m,V] = self.infer_diag(X_i,D_i)
     if sp.amin(V)<=-0.:
         class MJMError(Exception):
             pass
         print "negative/eq variance"
         print [m,V,X_i,D_i]
         print "_______________"
         #self.printc()
         raise(MJMError)
     if sp.amin(sp.var(m,axis=0))<-0.:
         class MJMError(Exception):
             pass
         print "negativevar of mean"
         print X_i.shape
         print [m,V,sp.var(m,axis=0),X_i,D_i]
         print "_______________"
         #self.printc()
         raise(MJMError)
     
     return [sp.mean(m,axis=0).reshape([1,ns]),(sp.mean(V,axis=0)+sp.var(m,axis=0)).reshape([1,ns])]

 def compute(i, j):
     if i == j:
         return 1.0
     elif trains[i].size <= 0 or trains[j].size <= 0:
         return 0.0
     else:
         diff_matrix = sp.absolute(trains[i] - sp.atleast_2d(trains[j]).T)
         return 0.5 * (
             sp.sum(kernel(sp.amin(diff_matrix, axis=0))) / trains[i].size +
             sp.sum(kernel(sp.amin(diff_matrix, axis=1))) / trains[j].size)

 def __init__(self, evaluator, hof1, hof2, **args):
     if "symmetric" in args:
         M = CiaoPlot.generateData(evaluator, hof1, hof2, symmetric=args["symmetric"])
         del args["symmetric"]
     else:
         M = CiaoPlot.generateData(evaluator, hof1, hof2)
     M *= 1 / (amin(M) - amax(M))
     M -= amin(M)
     self.relData = M
     ColorMap.__init__(self, M, minvalue=0, maxvalue=1, **args)

 def test_largest_sphere(self):
     net = OpenPNM.Network.Cubic(shape=[5, 5, 5], spacing=[0.1, 0.2, 0.3])
     geo = OpenPNM.Geometry.GenericGeometry(network=net, pores=net.Ps,
                                            throats=net.Ts)
     geo.models.add(propname='pore.diameter',
                    model=mods.largest_sphere,
                    iters=1)
     dmin = sp.amin(geo['pore.diameter'])
     assert dmin <= 0.1
     geo.models['pore.diameter']['iters'] = 5
     geo.regenerate()
     assert dmin < sp.amin(geo['pore.diameter'])

 def test_neighbor_max(self):
     catch = self.geo.pop('pore.seed', None)
     catch = self.geo.models.pop('pore.seed', None)
     catch = self.geo.models.pop('throat.seed', None)
     mod = gm.pore_misc.neighbor
     self.geo['throat.seed'] = sp.rand(self.net.Nt,)
     self.geo.models.add(model=mod,
                         propname='pore.seed',
                         throat_prop='throat.seed',
                         mode='max')
     assert sp.all(sp.in1d(self.geo['pore.seed'], self.geo['throat.seed']))
     pmin = sp.amin(self.geo['pore.seed'])
     tmin = sp.amin(self.geo['throat.seed'])
     assert pmin >= tmin

def test_late_pore_and_throat_filling():
    phys.models.add(propname='pore.fractional_filling',
                    model=OpenPNM.Physics.models.multiphase.late_pore_filling,
                    Pc=0,
                    Swp_star=0.2,
                    eta=1)
    mod = OpenPNM.Physics.models.multiphase.late_throat_filling
    phys.models.add(propname='throat.fractional_filling',
                    model=mod,
                    Pc=0,
                    Swp_star=0.2,
                    eta=1)
    phys.regenerate()
    drainage.setup(invading_phase=water, defending_phase=air,
                   pore_filling='pore.fractional_filling',
                   throat_filling='throat.fractional_filling')
    drainage.set_inlets(pores=pn.pores('boundary_top'))
    drainage.run()
    data = drainage.get_drainage_data()
    assert sp.amin(data['invading_phase_saturation']) == 0.0
    assert sp.amax(data['invading_phase_saturation']) < 1.0

    drainage.return_results(Pc=5000)
    assert 'pore.occupancy' in water.keys()
    assert 'throat.occupancy' in water.keys()
    assert 'pore.partial_occupancy' in water.keys()
    assert 'throat.partial_occupancy' in water.keys()

def overf_power_spectrum(amp, index, f0, dt, n, cut_off=0):
    """Calculates the theoretical f**`index` power spectrum.
    """
    
    if cut_off < 0:
        raise ValueError("Low frequency cut off must not be negative.")
    # Sometimes the fitting routines do something weird that causes
    # an overflow from a ridiculous index.  Limit the index.
    index = max(index, -20)
    # Get the frequencies represented in the FFT.
    df = 1.0/dt/n
    freq = sp.arange(n, dtype=float)
    freq[n//2+1:] -= freq[-1] + 1
    freq = abs(freq)*df
    # 0th (mean) mode is meaningless has IR divergence.  Deal with it later (in
    # the cut off.
    freq[0] = 1
    # Make the power spectrum.
    power = (freq/f0)**index
    power *= amp
    # Restore frequency of mean mode.
    freq[0] = 0
    # Find the power just above the cut off frequency.
    p_cut = power[sp.amin(sp.where(freq > cut_off)[0])]
    # Flatten off the power spectrum.
    power[freq <= cut_off] = p_cut
    return power

    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 neighbor(geometry, pore_prop='pore.seed', mode='min', **kwargs):
    r"""
    Adopt a value based on the values in neighboring pores

    Parameters
    ----------
    geometry : OpenPNM Geometry Object
        The object containing the ``pore_prop`` to be used.

    pore_prop : string
        The dictionary key to the array containing the pore property to be
        used in the calculation.  Default is 'pore.seed'.

    mode : string
        Controls how the throat property is calculated.  Options are 'min',
        'max' and 'mean'.

    """
    network = geometry._net
    throats = network.throats(geometry.name)
    P12 = network.find_connected_pores(throats)
    pvalues = network[pore_prop][P12]
    if mode == 'min':
        value = _sp.amin(pvalues, axis=1)
    if mode == 'max':
        value = _sp.amax(pvalues, axis=1)
    if mode == 'mean':
        value = _sp.mean(pvalues, axis=1)
    return value

def neighbor(geometry, network, pore_prop='pore.seed', mode='min', **kwargs):
    r"""
    Adopt a value based on the values in the neighboring pores

    Parameters
    ----------
    mode : string
        Indicates how to select the values from the neighboring pores.  The
        options are:

        - min : (Default) Uses the minimum of the value found in the neighbors
        - max : Uses the maximum of the values found in the neighbors
        - mean : Uses an average of the neighbor values

    pore_prop : string
        The dictionary key containing the pore property to be used.
    """
    throats = network.throats(geometry.name)
    P12 = network.find_connected_pores(throats)
    pvalues = network[pore_prop][P12]
    if mode == 'min':
        value = _sp.amin(pvalues, axis=1)
    if mode == 'max':
        value = _sp.amax(pvalues, axis=1)
    if mode == 'mean':
        value = _sp.mean(pvalues, axis=1)
    return value

def scale(x, M=None, m=None, REVERSE=None):
    """ Function that standardize the data
        Input:
            x: the data
            M: the Max vector
            m: the Min vector
        Output:
            x: the standardize data
            M: the Max vector
            m: the Min vector
    """
    if not sp.issubdtype(x.dtype, float):
        do_convert = 1
    else:
        do_convert = 0
    if REVERSE is None:
        if M is None:
            M = sp.amax(x, axis=0)
            m = sp.amin(x, axis=0)
            if do_convert:
                xs = 2 * (x.astype("float") - m) / (M - m) - 1
            else:
                xs = 2 * (x - m) / (M - m) - 1
            return xs, M, m
        else:
            if do_convert:
                xs = 2 * (x.astype("float") - m) / (M - m) - 1
            else:
                xs = 2 * (x - m) / (M - m) - 1
            return xs
    else:
        return (1 + x) / 2 * (M - m) + m

 def _do_one_outer_iteration(self, **kwargs):
     r"""
     One iteration of an outer iteration loop for an algorithm
     (e.g. time or parametric study)
     """
     # Checking for the necessary values in Picard algorithm
     nan_tol = sp.isnan(self['pore.source_tol'])
     nan_max = sp.isnan(self['pore.source_maxiter'])
     self._tol_for_all = sp.amin(self['pore.source_tol'][~nan_tol])
     self._maxiter_for_all = sp.amax(self['pore.source_maxiter'][~nan_max])
     if self._guess is None:
         self._guess = sp.zeros(self._coeff_dimension)
     t = 1
     step = 0
     # The main Picard loop
     while t > self._tol_for_all and step <= self._maxiter_for_all:
         X, t, A, b = self._do_inner_iteration_stage(guess=self._guess,
                                                     **kwargs)
         logger.info('tol for Picard source_algorithm in step ' +
                     str(step) + ' : ' + str(t))
         self._guess = X
         step += 1
     # Check for divergence
     self._steps = step
     if t >= self._tol_for_all and step > self._maxiter_for_all:
         raise Exception('Iterative algorithm for the source term reached '
                         'to the maxiter: ' + str(self._maxiter_for_all) +
                         ' without achieving tol: ' +
                         str(self._tol_for_all))
     logger.info('Picard algorithm for source term converged!')
     self.A = A
     self.b = b
     self._tol_reached = t
     return X

 def test_calculates_frequencies(self) :
     self.Data.calc_freq()
     self.assertTrue(hasattr(self.Data, 'freq'))
     self.assertEqual(len(self.Data.freq), self.nfreq)
     self.assertAlmostEqual(self.Data.field['BANDWID'], 
                            sp.amax(self.Data.freq) -
                            sp.amin(self.Data.freq), -5)

def plot_pairwise_velocities_r(case,color,all_radial_distances,all_radial_velocities):
    dr = 0.3 # Mpc/h
    rmin, rmax = sp.amin(all_radial_distances), sp.amax(all_radial_distances) 
    rrange = rmax-rmin 
    N = int(sp.ceil(rrange/dr))
    rs = sp.linspace(rmin,rmax,N)
    v12_of_r = [[] for index in range(N)]
    
    for r,v12 in zip(all_radial_distances,all_pairwise_velocities):
    
        index = int(sp.floor((r-rmin)/dr))
        v12_of_r[index].append(v12)
            
    
    sigma_12s = sp.zeros(N)
    v12_means = sp.zeros(N)
    for index in range(len(sigma_12s)):
        v12_of_r_index = sp.array(v12_of_r[index])
        print "number of counts in the", index,"th bin:", len(v12_of_r_index)
        sigma_12 = sp.sqrt(sp.mean(v12_of_r_index**2))
        v12_mean = -sp.mean(v12_of_r_index)
        sigma_12s[index] = sigma_12
        v12_means[index] = v12_mean
    
    
    plt.plot(rs,sigma_12s,color=color,label='$\sigma_{12}$')
    plt.plot(rs,v12_means,color=color,label='$|v_{12}|$')
    plt.xlabel('r [Mpc/h]')
    plt.ylabel('[km/s]')
    plt.xscale('log')
    plt.axis([0.5,100,0,600])

    def __init__(self, mat, cmap=None, pixelspervalue=20, minvalue=None, maxvalue=None):
        """ Make a colormap image of a matrix or sequence of Matrix/Connection objects

        :key mat: the matrix to be used for the colormap.
        """
        if isinstance(mat, (ParameterContainer, Connection)):
            mat = reshape(mat.params, (mat.outdim, mat.indim))

        if not isinstance(mat, ndarray):
            raise ValueError("Don't know how to display a ColorMap for a matrix of type {}".format(type(mat)))
        if minvalue == None:
            minvalue = amin(mat)
        if maxvalue == None:
            maxvalue = amax(mat)
        if isinstance(cmap, basestring) and cmap.strip():
            cmap = getattr(cm, cmap.lower().strip())
        if not cmap:
            cmap = cm.hot

        figsize = (array(mat.shape) / 100. * pixelspervalue)[::-1]
        self.fig = figure(figsize=figsize)
        axes([0, 0, 1, 1]) # Make the plot occupy the whole canvas
        axis('off')
        self.fig.set_size_inches(figsize)
        imshow(mat, cmap=cmap, clim=(minvalue, maxvalue), interpolation='nearest')

def averageLogDistribution(values, log_step=0.2, first_point=None, last_point=None):
    """
    calculates the <values> vs. xVariable in log scale

    Parameters:
    ---------------
    values : dict
    A dictionary where the keys are the xValues
    and each element contains an array-like sequence of data
    Example:
    [1: array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]),
     2: array([1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2]),
     
    values : ndarray
    Two columns array of xValues and yValues,
    to be rearranged as above
     
    Returns:
    center point of the bin, average value within the bin

    """
    # Check the list of Values
    if not checkIfVoid(values):
        print("Error")
    if isinstance(values, dict):
        xValues = np.asarray(values.keys())
        yValues = np.asarray(values.values())
    elif isinstance(values, np.ndarray):
        xValues = np.unique(values[:, 0])
        yValues = []
        for xVal in xValues:
            index = values[:, 0] == xVal
            yValues.append(values[index, 1])
        yValues = scipy.array(yValues)
    else:
        print("Values shape not recognized")
        return
    if not first_point:
        first_point = scipy.amin(xValues) * 0.99
    if not last_point:
        last_point = scipy.amax(xValues) * 1.01

    xbins, bins = getLogBins(first_point, last_point, log_step)
    yAverage = []
    for i, j in zip(bins[:-1], bins[1:]):
        q1, q2 = np.greater_equal(xValues, i), np.less(xValues, j)
        q = np.logical_and(q1, q2)
        if sum(q) == 0:
            averageValue = np.NaN
        else:
            allElements = [val for val in itertools.chain(*yValues[q])]
            averageValue = sum(allElements) / float(len(allElements))
            # print averageValue, allElements
        yAverage.append(averageValue)
    yAverage = np.asanyarray(yAverage)
    # Check if there are NaN values
    iNan = np.isnan(yAverage)
    x = xbins[~iNan]
    y = yAverage[~iNan]
    return x, y

def late_pore_filling(physics,
                      phase,
                      network,
                      Pc,
                      Swp_star=0.2,
                      eta=3,
                      wetting_phase=False,
                      pore_occupancy='pore.occupancy',
                      throat_capillary_pressure='throat.capillary_pressure',
                      **kwargs):
    r'''
    Applies a late pore filling model to calculate fractional pore filling as 
    a function of applied capillary pressure.
    
    Parameters
    ----------
    Pc : float
        The capillary pressure in the non-wetting phase (Pc > 0)
    Swp_star : float
        The residual wetting phase in an invaded pore immediately after
        nonwetting phase invasion
    eta : float
        Exponent to control the rate at which wetting phase is displaced
    wetting_phase : boolean
        Indicates whether supplied phase is the wetting or non-wetting phase
    
    
    '''
    pores = phase.pores(physics.name)
    prop = phase[throat_capillary_pressure]
    neighborTs = network.find_neighbor_throats(pores,flatten=False)
    Pc_star = sp.array([sp.amin(prop[row]) for row in neighborTs])
    Swp = Swp_star*(Pc_star/Pc)**eta
    if wetting_phase:
        values = Swp*phase[pore_occupancy]*(Pc_star<Pc)
    else:
        values = (1-Swp)*(1-phase[pore_occupancy])*(Pc_star<Pc)
    return values
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    

def main(database):

    # Commits per committer limited to the 30 first with the highest accumulated activity
    query = "select count(*) from scmlog group by committer_id order by count(*) desc limit 40"

    # Connecting to the data base and retrieving data
    connector = connect(database)
    results = int(connector.execute(query))
    if results > 0:
        results_aux = connector.fetchall()
    else:
        print ("Error when retrieving data")
        return

    # Moving data to a list
    commits = []
    for commit in results_aux[5:]:
        #   for commits in results_aux:
        commits.append(int(commit[0]))

    # Calculating basic statistics
    print "max: " + str(sp.amax(commits))
    print "min: " + str(sp.amin(commits))
    print "mean: " + str(sp.mean(commits))
    print "median: " + str(sp.median(commits))
    print "std: " + str(sp.std(commits))
    print ".25 quartile: " + str(sp.percentile(commits, 25))
    print ".50 quartile: " + str(sp.percentile(commits, 50))
    print ".75 quartile: " + str(sp.percentile(commits, 75))

def test_residual_and_lpf():
    phys.models.add(propname='pore.fractional_filling',
                    model=OpenPNM.Physics.models.multiphase.late_pore_filling,
                    Pc=0,
                    Swp_star=0.2,
                    eta=1)
    phys.models.add(propname='throat.fractional_filling',
                    model=OpenPNM.Physics.models.multiphase.late_throat_filling,
                    Pc=0,
                    Swp_star=0.2,
                    eta=1)
    phys.regenerate()
    drainage.setup(invading_phase=water, defending_phase=air,
                   pore_filling='pore.fractional_filling',
                   throat_filling='throat.fractional_filling')
    drainage.set_inlets(pores=pn.pores('boundary_top'))
    resPs = pn.pores('internal')[sp.random.random(len(pn.pores('internal')))<0.1]
    resTs = pn.throats('internal')[sp.random.random(len(pn.throats('internal')))<0.1]
    drainage.set_residual(pores=resPs, throats=resTs)
    drainage.run()
    drainage.return_results(Pc=5000)
    data = drainage.get_drainage_data()
    assert sp.all(water["pore.partial_occupancy"][resPs] == 1.0)
    assert sp.all(water["throat.partial_occupancy"][resTs] == 1.0)
    assert sp.amin(data['invading_phase_saturation']) > 0.0
    assert sp.amax(data['invading_phase_saturation']) < 1.0
    assert sp.all(water["pore.occupancy"]+air["pore.occupancy"] == 1.0)
    total_pp = water["pore.partial_occupancy"]+air["pore.partial_occupancy"]
    assert sp.all(total_pp == 1.0)
    assert sp.all(water["throat.occupancy"]+air["throat.occupancy"] == 1.0)
    total_pt = water["throat.partial_occupancy"]+air["throat.partial_occupancy"] 
    assert sp.all(total_pt == 1.0)