def test_domain_dimensions(self):
     face1 = self.net.pores('bottom_boundary')
     face2 = self.net.pores('top_boundary')
     assert sp.around(self.net.domain_length(face1, face2) / self.scale,
                      3) == 1.0
     assert sp.around(self.net.domain_area(face1) / self.scale**2,
                      3) == 1.0

def test_distance_center():
    shape = sp.array([7, 5, 9])
    spacing = sp.array([2, 1, 0.5])
    pn = OpenPNM.Network.Cubic(shape=shape, spacing=spacing)
    sx, sy, sz = spacing
    center_coord = sp.around(topology.find_centroid(pn['pore.coords']), 7)
    cx, cy, cz = center_coord
    coords = pn['pore.coords']
    x, y, z = coords.T
    coords = sp.concatenate((coords, center_coord.reshape((1, 3))))
    pn['pore.center'] = False
    mask1 = (x <= (cx + sx/2)) * (y <= (cy + sy/2)) * (z <= (cz + sz/2))
    mask2 = (x >= (cx - sx/2)) * (y >= (cy - sy/2)) * (z >= (cz - sz/2))
    center_pores_mask = pn.Ps[mask1 * mask2]
    pn['pore.center'][center_pores_mask] = True
    center = pn.Ps[pn['pore.center']]
    L1 = sp.amax(topology.find_pores_distance(network=pn,
                                              pores1=center,
                                              pores2=pn.Ps))
    L2 = sp.amax(topology.find_pores_distance(network=pn,
                                              pores1=pn.Ps,
                                              pores2=pn.Ps))
    l1 = ((shape[0] - 1) * sx) ** 2
    l2 = ((shape[1] - 1) * sy) ** 2
    l3 = ((shape[2] - 1) * sz) ** 2
    L3 = sp.sqrt(l1 + l2 + l3)
    assert sp.around(L1 * 2, 7) == sp.around(L2, 7)
    assert sp.around(L2, 7) == sp.around(L3, 7)

def xover(chrom,N,p):
    """Single point crossover with probability N,precision p
    """
    N = round(chrom.shape[0]*N)
    index1 = scipy.arange(chrom.shape[0])
    index2 = scipy.unique(scipy.around(scipy.rand(chrom.shape[0],)*chrom.shape[0]))[0:chrom.shape[0]/2]
    sel1,sel2 = [],[]
    for i in range(len(index1)):
        if index1[i] not in index2:
            sel1.append(index1[i])
        else:
            sel2.append(index1[i])
    select1 = sel1[0:min([int(round(len(sel1)*N)),int(round(len(sel2)*N))])]
    select2 = sel2[0:min([int(round(len(sel1)*N)),int(round(len(sel2)*N))])]       
    
    # set xover points
    xoverpnt = scipy.around(scipy.rand(len(select1),)*(chrom.shape[1]-1))
    
    # perform xover
    nchrom = copy.deepcopy(chrom)
    for i in range(len(select1)):
        try:
            slice1 = chrom[select1[i],0:int(xoverpnt[i])]
            slice2 = chrom[select2[i],0:int(xoverpnt[i])]
            nchrom[select2[i],0:int(xoverpnt[i])] = slice1
            nchrom[select1[i],0:int(xoverpnt[i])] = slice2
        except:
            nchrom = nchrom
    
    return nchrom

 def generate(self, **params):
     r'''
     radius: int
         The physical radius of the netowrk.  The number of pores determined
         from this and the lattice_spacing parameter.
     lattice_spacing : float
         The lattice constant for the network, used to scale distance between pores.
     '''
     self._logger.info(sys._getframe().f_code.co_name+": Start of network topology generation")
     #create spherical template, then call normal template procedure
     self._radius = params['radius']
     self._length = params['length']
     Nr = sp.around(self._radius/params['lattice_spacing'])
     Lz = sp.around(self._length/params['lattice_spacing'])
     temp = sp.ones((2*Nr, 2*Nr, Lz))
     temp[Nr, Nr, :] = 0
     temp = spim.distance_transform_edt(temp)
     temp = temp < Nr
     params['template'] = temp
     #Call standard generation protocol
     self._generate_setup(**params)
     self._generate_pores()
     self._generate_throats()
     self._add_boundaries()
     self._add_labels()
     self._logger.debug(sys._getframe().f_code.co_name+": Network generation complete")        
     return self

    def _get_fibre_slice(self, plane=None, index=None):
        r"""
        Plot an image of a slice through the fibre image
        plane contains percentage values of the length of the image in each axis

        Parameters
        ----------
        plane : array_like
        List of 3 values, [x,y,z], 2 must be zero and the other must be between
        zero and one representing the fraction of the domain to slice along
        the non-zero axis

        index : array_like
        similar to plane but instead of the fraction an index of the image is used
        """
        if hasattr(self, '_fibre_image') is False:
            logger.warning('This method only works when a fibre image exists, ' +
                           'please run make_fibre_image')
            return None
        if plane is None and index is None:
            logger.warning('Please provide either a plane array or index array')
            return None
        if self._fibre_image is None:
            self.make_fibre_image()

        if plane is not None:
            if 'array' not in plane.__class__.__name__:
                plane = sp.asarray(plane)
            if sp.sum(plane == 0) != 2:
                logger.warning('Plane argument must have two zero valued ' +
                               'elements to produce a planar slice')
                return None
            l = sp.asarray(sp.shape(self._fibre_image))
            s = sp.around(plane*l).astype(int)
        elif index is not None:
            if 'array' not in index.__class__.__name__:
                index = sp.asarray(index)
            if sp.sum(index == 0) != 2:
                logger.warning('Index argument must have two zero valued ' +
                               'elements to produce a planar slice')
                return None
            if 'int' not in str(index.dtype):
                index = sp.around(index).astype(int)
            s = index

        if s[0] != 0:
            slice_image = self._fibre_image[s[0], :, :]
        elif s[1] != 0:
            slice_image = self._fibre_image[:, s[1], :]
        else:
            slice_image = self._fibre_image[:, :, s[2]]

        return slice_image

    def write_extraction_mask(self):
        """ write both the .roi file and the tif pages
        ROI file format specification: each row a ROI,
        first col: kw for roi type
        if circle: label, layer, pos x, pos y, diameter
        if poly: label, layer, pos x_1, pos_y1, ... pos x_n, pos y_n
        """
        outpath = self.SaveFileDialog(title='saving ROIs',default_dir = self.Main.Options.general['cwd'], extension='*.roi')
        outpath = self.append_extension(outpath, '.roi')

        fh = open(outpath,'w')

        # iterate over ROIs
        for i,ROI in enumerate(self.Main.ROIs.ROI_list):
            label = str(ROI.label)

            if type(ROI) == myCircleROI:
#                pos = sp.array([ROI.pos().x(),ROI.pos().y()])
#                pos = self.get_ROI_position(ROI)
                x = str(sp.around(ROI.center[0],decimals=2))
                y = str(sp.around(ROI.center[1],decimals=2))
                d = str(sp.around(ROI.diameter,decimals=2))

                fh.write('\t'.join(['circle',label,x,y,d,'\n']))

            if type(ROI) == myPolyLineROI:
                 fh.write('\t'.join(['polygon',label]))
                 fh.write('\t')

                 handle_pos = [tup[1] for tup in ROI.getSceneHandlePositions()]
                 pos_mapped = [ROI.ViewBox.mapToView(pos) for pos in handle_pos]

                 for pos in (pos_mapped):

                     x = pos.x()
                     y = pos.y()

                     fh.write('\t'.join([str(sp.around(x,decimals=2)),str(sp.around(y,decimals=2))]))
                     fh.write('\t')
                 fh.write('\n')

        fh.close()

        print("saved ROIs in .roi format to", outpath)
#        outpath = os.path.splitext(outpath)[0] + '_mask.tif'
#        outpath = self.MainWindow.SaveFileDialog(title='saving ROIs',defaultdir=self.path,extension='.tif')

        # extraction mask
        self.Main.Processing.calc_extraction_mask()
        io.save_tstack(self.Main.Data.extraction_mask.astype('uint16'),os.path.splitext(outpath)[0] + '_mask.tif')
        print("saved ROIs in .tif format to", outpath)

        pass

def get_cb_ticks(values):
    min_tick = sp.nanmin(values)
    max_tick = sp.nanmax(values)
    med_tick = min_tick + (max_tick - min_tick) / 2.0
    if max_tick > 1.0:
        min_tick = sp.ceil(min_tick)
        max_tick = sp.floor(max_tick)
        med_tick = sp.around(med_tick)
    else:
        min_tick = sp.ceil(min_tick * 100.0) / 100.0
        max_tick = sp.floor(max_tick * 100.0) / 100.0
        med_tick = sp.around(med_tick, 2)
    return [min_tick, med_tick, max_tick]

    def on_shifted_dwp_curves(self, t):
        a = P4Rm()
        if a.AllDataDict['model'] == 0:
            temp_1 = arange(2, len(a.ParamDict['dwp'])+1)
            temp_2 = temp_1 * t / (len(a.ParamDict['dwp']))
            P4Rm.ParamDict['x_dwp'] = t - temp_2
            shifted_dwp = a.ParamDict['dwp'][:-1:]
            temp_3 = in1d(around(a.ParamDict['depth'], decimals=3),
                          around(a.ParamDict['x_dwp'], decimals=3))
            temp_4 = a.ParamDict['DW_i'][temp_3]
            P4Rm.ParamDict['scale_dw'] = shifted_dwp / temp_4
            P4Rm.ParamDict['scale_dw'][a.ParamDict['scale_dw'] == 0] = 1.

            P4Rm.ParamDict['DW_shifted'] = shifted_dwp/a.ParamDict['scale_dw']
            P4Rm.ParamDict['dw_out'] = a.ParamDict['dwp'][-1]

        elif a.AllDataDict['model'] == 1:
            temp_1 = arange(0, len(a.ParamDict['dwp'])+1-3)
            temp_2 = temp_1 * t / (len(a.ParamDict['dwp'])-3)
            P4Rm.ParamDict['x_dwp'] = t - temp_2
            shifted_dwp = a.ParamDict['dwp'][1:-1:]
            temp_3 = in1d(around(a.ParamDict['depth'], decimals=3),
                          around(a.ParamDict['x_dwp'], decimals=3))
            temp_4 = a.ParamDict['DW_i'][temp_3]
            P4Rm.ParamDict['scale_dw'] = shifted_dwp / temp_4
            P4Rm.ParamDict['scale_dw'][a.ParamDict['scale_dw'] == 0] = 1.

            P4Rm.ParamDict['DW_shifted'] = shifted_dwp/a.ParamDict['scale_dw']
            temp_5 = array([a.ParamDict['dwp'][0], a.ParamDict['dwp'][-1]])
            P4Rm.ParamDict['dw_out'] = temp_5

        elif a.AllDataDict['model'] == 2:
            x_dw_temp = []
            x_dw_temp.append(t*(1-a.ParamDict['dwp'][1]))
            x_dw_temp.append(t*(1-a.ParamDict['dwp'][1] +
                             a.ParamDict['dwp'][2]/2))
            x_dw_temp.append(t*(1-a.ParamDict['dwp'][1] -
                             a.ParamDict['dwp'][3]/2))
            x_dw_temp.append(t*0.05)
            P4Rm.ParamDict['x_dwp'] = x_dw_temp

            y_dw_temp = []
            y_dw_temp.append(a.ParamDict['dwp'][0])
            y_dw_temp.append(1. - (1-a.ParamDict['dwp'][0])/2)
            y_dw_temp.append(1. - (1-a.ParamDict['dwp'][0])/2 -
                             (1-a.ParamDict['dwp'][6])/2)
            y_dw_temp.append(a.ParamDict['dwp'][6])
            P4Rm.ParamDict['DW_shifted'] = y_dw_temp

def plot_median_errors(RefinementLevels):
        for i in RefinementLevels[0].cases:
            x =[];
            y =[];
            print "Analyzing median error on: ", i ;
            for r in RefinementLevels:                
                x.append(r.LUT.D_dim*r.LUT.P_dim)
                r.get_REL_ERR_SU2(i)
                y.append(r.SU2[i].median_ERR*100)
            
            x = sp.array(x)
            y = sp.array(y)            
            y = y[sp.argsort(x)]
            x = x[sp.argsort(x)]
                                    
            LHM = sp.ones((len(x),2))
            RHS = sp.ones((len(x),1))            
            LHM[:,1] = sp.log10(x)
            RHS[:,0] = sp.log10(y)

            sols = sp.linalg.lstsq(LHM,RHS)
            b = -sols[0][1]
            plt.loglog(x,y, label='%s, %s'%(i,r'$O(\frac{1}{N})^{%s}$'%str(sp.around(b,2))), basex=10, basey=10, \
                       subsy=sp.linspace(10**(-5), 10**(-2),20),\
                       subsx=sp.linspace(10**(2), 10**(5),50))
            
            #for r in RefinementLevels:                
               # x.append(r.LUT.D_dim*r.LUT.P_dim)
              #  r.get_REL_ERR_SciPy(i)
             #   y.append(r.SciPy[i].median_ERR*100)
            #plt.plot(x,y, label='SciPy: %s'%i)
        plt.grid(which='both')
        plt.xlabel('Grid Nodes (N)')
        plt.ylabel('Median relative error [%]')
        return;

 def __init__(self, shape, num_points=None, **kwargs):
     # Generate Delaunay tessellation from super class, then trim
     super().__init__(shape=shape, num_points=num_points, **kwargs)
     points = self['pore.coords']
     conns = self['throat.conns']
     # Find centroid of each pair of nodes
     c = points[conns]
     m = (c[:, 0, :] + c[:, 1, :])/2
     # Find radius of circle connecting each pair of nodes
     r = sp.sqrt(sp.sum((c[:, 0, :] - c[:, 1, :])**2, axis=1))/2
     # Use KD-Tree to find distance to nearest neighbors
     tree = sptl.cKDTree(points)
     n = tree.query(x=m, k=1)[0]
     # Identify throats whose centroid is not near an unconnected node
     g = sp.around(n, decimals=5) == sp.around(r, decimals=5)
     trim(self, throats=~g)

 def test_compress_geom(self):
     b1 = self.net.pores('bottom_boundary')
     b2 = self.net.pores('top_boundary')
     height1 = self.net.domain_length(b1, b2)
     self.geo_vox.compress_geometry(factor=[1, 1, 0.5])
     height2 = self.net.domain_length(b1, b2)
     assert sp.around(height1/height2, 5) == 2.0

    def rekernel(self):
        dcontrol = self.control[ord('d')]
        econtrol = self.control[ord('e')]
        rcontrol = self.control[ord('r')]
        radius = rcontrol.val
        dvalue = dcontrol.val
        evalue = econtrol.val
        rmax   = rcontrol.limit[2]
        if self.rlast != radius:
            inner, outer = float(radius-1), float(radius)
            shape      = (self.edge, self.edge)
            self.radii = list(product(arange(-rmax,rmax+1,1.0), repeat=2))
            self.radii = array([sqrt(x*x+y*y) for x,y in self.radii]).reshape(shape)

            if True:
                self.negative = -exp(-dvalue*(self.radii-outer)**2)
                self.positive = +exp(-dvalue*(self.radii-inner)**2)
            else:
                self.radii = around(self.radii)
                self.negative = zeros((self.edge,self.edge),dtype=float)
                self.negative[self.radii == outer] = -1.0
                self.positive = zeros(shape,dtype=float)
                self.positive[self.radii == inner] = +1.0

            self.negative /= fabs(self.negative.sum())
            self.positive /= fabs(self.positive.sum())

            self.kernel = self.negative + self.positive
            self.rlast = radius
        if self.elast != evalue:
            self.gauss = exp(-evalue * self.radii**2)
            self.gauss /= self.gauss.sum()
            self.elast = evalue

def crtpop(ni,nv,prec):
    """Create a random population array of size
    ni by nv in the range 0:preci-1.  Use prec = 2
    to create binary string
    """
    pop = scipy.around(scipy.rand(ni,nv)*(prec-1))
    
    return pop

def toRanks(A):
    """
    converts the columns of A to ranks
    """
    AA=sp.zeros_like(A)
    for i in range(A.shape[1]):
        AA[:,i] = st.rankdata(A[:,i])
    AA=sp.array(sp.around(AA),dtype="int")-1
    return AA

def mutate(chrom,N,p):
    """Mutation with probability N and precision p
    """
    index = []
    for x in range(int(scipy.around(chrom.shape[0]*chrom.shape[1]*N))):
        index.append((int(scipy.around(scipy.rand(1,)[0]*(chrom.shape[0]-1))),
        int(scipy.around(scipy.rand(1,)[0]*(chrom.shape[1]-1)))))

    for x in index:
        if p == 1:
            if chrom[x] == 1:
                chrom[x] = 0
            else:
                chrom[x] = 1
        else:
            chrom[x] = int(scipy.around(scipy.rand(1,)[0]*(p-1)))
    
    return chrom

 def test_general_toroidal(self):
     phys = self.phys
     r_tor = 1e-6
     phys.add_model(propname='throat.purcell_pressure',
                    model=pm.capillary_pressure.purcell,
                    r_toroid=r_tor)
     phys['throat.scale_a'] = r_tor
     phys['throat.scale_b'] = r_tor
     phys.add_model(propname='throat.general_pressure',
                    model=pm.meniscus.general_toroidal,
                    mode='max',
                    num_points=1000)
     a = sp.around(phys['throat.purcell_pressure'], 10)
     b = sp.around(phys['throat.general_pressure'], 10)
     assert sp.allclose(a, b)
     h = phys.check_data_health()
     for check in h.values():
         if len(check) > 0:
             assert 1 == 2

 def NormLog(self,X,Y,parameterValues, independentValues):
     """
     This kind of normalization is correct
     if the data are uniform in log scale,
     as prepared by our code toBinDistributions.py
     """
     lgX = scipy.log10(X)
     D = scipy.around(lgX[1] - lgX[0],2)
     bins = 10**(lgX+D/2.) - 10**(lgX-D/2.)
     return Y/sum(Y*bins)

def _assert_files_equal_testing(e, a):

    da = sp.loadtxt(a, dtype='str', delimiter='\t')
    de = sp.loadtxt(e, dtype='str', delimiter='\t')

    ### check header
    assert sp.all(da[0, :] == de[0, :])
    da = da[1:, ]
    de = de[1:, ]

    ### check text cols
    assert sp.all(da[:, [0, 1]] == de[:, [0, 1]])
    da = da[:, 2:]
    de = de[:, 2:]

    ### check p-values (up to certain precision)
    da = sp.around(da.astype('float'), decimals=6)
    de = sp.around(de.astype('float'), decimals=6)
    assert sp.all(da == de)

 def is_near_constant(self, min_num_diff=10):
     vals = sp.array(self.values)
     if sp.std(vals) > 0:
         vals = 50 * (vals - sp.mean(vals)) / sp.std(vals)
         vals = vals - vals.min() + 0.1
         b_counts = sp.bincount(sp.array(sp.around(vals), dtype='int'))
         b = b_counts.max() > len(vals) - min_num_diff
         return b
     else:
         return True

 def is_near_constant(self, pid, min_num_diff=10):
     vals = sp.array(self.phen_dict[pid]["values"])
     if sp.std(vals) > 0:
         vals = 50 * (vals - sp.mean(vals)) / sp.std(vals)
         vals = vals - vals.min() + 0.1
         b_counts = sp.bincount(sp.array(sp.around(vals), dtype="int"))
         b = b_counts.max() > len(vals) - min_num_diff
         return b
     else:
         return True