def filter_introns(introns, genes, options):
    
    ### build interval trees of all genes starts and ends
    chrms = sp.array([_.strand for _ in genes])
    strands = sp.array([_.chr for _ in genes])
    gene_trees = dict()
    for c in sp.unique(chrms):
        for s in sp.unique(strands):
            gene_trees[(c, s)] = it.IntervalTree()
            c_idx = sp.where((chrms == c) & (strands == s))[0]
            for i in c_idx:
                gene_trees[(c, s)][genes[i].start:genes[i].stop] = i

    ### match all introns agains trees and remove elements overlapping
    ### more than one gene on the same chr/strand
    cnt_tot = 0
    cnt_rem = 0
    strand_list = ['+', '-']
    offset = options.intron_edges['append_new_terminal_exons_len']
    for si, s in enumerate(strand_list):
        for i in range(introns.shape[0]):
            if introns[i, si].shape[0] == 0:
                continue
            k_idx = []
            cnt_tot += introns[i, si].shape[0]
            for j in range(introns[i, si].shape[0]):
                if len(gene_trees[(s, genes[i].chr)].overlap(introns[i, si][j, 0] - offset, introns[i, si][j, 1] + offset)) == 1:
                    k_idx.append(j)
            if len(k_idx) < introns[i, si].shape[0]:
                cnt_rem += (introns[i, si].shape[0] - len(k_idx))
                introns[i, si] = introns[i, si][k_idx, :]
    print('removed %i of %i (%.2f percent) introns overlapping to no or multiple genes' % (cnt_rem, cnt_tot, cnt_rem / float(max(cnt_tot, 1)) * 100))

    return introns

def main():
    args = getArguments(getParser())

    # prepare logger
    logger = Logger.getInstance()
    if args.debug:
        logger.setLevel(logging.DEBUG)
    elif args.verbose:
        logger.setLevel(logging.INFO)

    # check if output image exists
    if not args.force:
        if os.path.exists(args.output):
            logger.warning("The output image {} already exists. Exiting.".format(args.output))
            exit(-1)

    # load input image
    input_data, input_header = load(args.input)

    logger.debug("Old number of regions={}.".format(len(scipy.unique(input_data))))

    # cut and relabel along the required dimension
    logger.info("Cutting and relabeling...")
    dimensions = range(input_data.ndim)
    del dimensions[args.dimension]
    __split_along(input_data, dimensions)

    logger.debug("New number of regions={}.".format(len(scipy.unique(input_data))))

    # save result contour volume
    save(input_data, args.output, input_header, args.force)

    logger.info("Successfully terminated.")

 def _read_sky_logfile(self):
     #TODO : expand to read errors, msgs etc
     # read in the whole sky log file, shouldn't be big
     f = open(self.skylogfile)
     lines = f.readlines()
     f.close()
     dust = [line.split()[1:] for line in lines if line.startswith('dtau_dust')]
     line = [line.split()[1:] for line in lines if line.startswith('dtau_line')]
     dust = _sp.array(dust, dtype='float')
     line = _sp.array(line, dtype='float')
     transitions = _sp.unique(dust[:,0])
     shells = _sp.unique(dust[:,1])
     dtau_dust = dict()
     dtau_line = dict()
     dtau_tot = dict()
     for t in transitions:
         d = []
         l = []
         for s in shells:
             d.append( _sp.mean([i[2] for i in dust if ((i[0]==t) * (i[1]==s))]) )
             l.append( _sp.mean([i[2] for i in line if ((i[0]==t) * (i[1]==s))]) )
         dtau_dust[t] = _sp.copy(d)
         dtau_line[t] = _sp.copy(l)
         dtau_tot[t] = _sp.array(d) + _sp.array(l)
     # create object to store in main class
     class Tau(object):pass
     Tau.dtau_dust = dtau_dust
     Tau.dtau_line = dtau_line
     Tau.dtau_tot = dtau_tot
     Tau.transitions = transitions
     Tau.shells = shells
     self.Tau = Tau

def Capillary_Pressure_Curve(net,
                             fluid,
                             capillary_pressure='capillary_pressure',
                             pore_volume='volume',
                             throat_volume='volume',
                             fig=None):
  r"""
  Plot drainage capillary pressure curve

  Parameters
  ----------
  net : OpenPNM Network Object
      The network for which the graphs are desired
  fig : Matplotlib figure object
      Canvas on which to draw plots

  """
  if type(fluid)==str: fluid = net.find_object_by_name(fluid)
  try:
    PcPoints = sp.unique(fluid.get_throat_data(prop=capillary_pressure))
  except KeyError:
    raise Exception('Capillary pressure simulation has not been run')

  PcPoints = sp.unique(fluid.get_throat_data(prop=capillary_pressure))
  Snwp = sp.zeros_like(PcPoints)
  Ps = sp.r_[0:net.num_pores('internal')]
  for i in range(1,sp.size(PcPoints)):
      Pc = PcPoints[i]
      Snwp[i] = sum((fluid.get_throat_data(prop=capillary_pressure)[Ps]<Pc)*(net.get_throat_data(prop=throat_volume)[Ps]))/sum(net.get_throat_data(prop=throat_volume)[Ps])
  
  if fig==None: fig = plt.figure()
  ax = fig.add_subplot(111)
  ax.plot(PcPoints,Snwp,'r.-')
  ax.set_xlabel('Capillary Pressure')
  ax.set_ylabel('Fluid Saturation')

def batch_metrics(unit_list, threshold, t_ref, t_cen):
    ''' This here function runs metrics on a batch of data.  Pass in units from the catalog.
    '''
    
    from scipy import unique
    
    samp_rate = 30000.
    n_samples = 30
    n_chans = 4
    
    # Find common Sessions
    sessions = unique([unit.session for unit in unit_list])
    
    for session in sessions:
        
        units = session.units
        tetrodes = unique([unit.tetrode for unit in units])
        
        for tetrode in tetrodes:
            data = load_spikes(session.path, tetrode,  samp_rate, n_samples, n_chans)
            f_p, f_n = metrics(data, threshold, t_ref, t_cen, session.duration)
            # Doing this because sometimes there is no cluster 0 sometimes
            f_p.setdefault(1)
            f_n.setdefault(1)
            units = [ unit for unit in session.units if unit.tetrode == tetrode] 
            for unit in units:
                unit.falsePositive = f_p[unit.cluster]
                unit.falseNegative = f_n[unit.cluster]

    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 _check_trapping(self, inv_val):
        r"""
        Determine which pores and throats are trapped by invading phase.  This
        method is called by ``run`` if 'trapping' is set to True.
        """
        # Generate a list containing boolean values for throat state
        Tinvaded = self['throat.inv_Pc'] < sp.inf
        # Add residual throats, if any, to list of invaded throats
        Tinvaded = Tinvaded + self['throat.residual']
        # Invert logic to find defending throats
        Tdefended = ~Tinvaded
        [pclusters, tclusters] = self._net.find_clusters2(mask=Tdefended,
                                                          t_labels=True)
        # See which outlet pores remain uninvaded
        outlets = self['pore.outlets']*(self['pore.inv_Pc'] == sp.inf)
        # Identify clusters connected to remaining outlet sites
        def_clusters = sp.unique(pclusters[outlets])
        temp = sp.in1d(sp.unique(pclusters), def_clusters, invert=True)
        trapped_clusters = sp.unique(pclusters)[temp]
        trapped_clusters = trapped_clusters[trapped_clusters >= 0]

        # Find defending clusters NOT connected to the outlet pores
        pmask = np.in1d(pclusters, trapped_clusters)
        # Store current applied pressure in newly trapped pores
        pinds = (self['pore.trapped'] == sp.inf) * (pmask)
        self['pore.trapped'][pinds] = inv_val

        # Find throats on the trapped defending clusters
        tinds = self._net.find_neighbor_throats(pores=pinds,
                                                mode='intersection')
        self['throat.trapped'][tinds] = inv_val
        self['throat.entry_pressure'][tinds] = 1000000

def test_linear_solvers():
    pn = OpenPNM.Network.Cubic([1, 40, 30], spacing=0.0001)
    geom = OpenPNM.Geometry.Toray090(network=pn,
                                     pores=pn.pores(),
                                     throats=pn.throats())
    air = OpenPNM.Phases.Air(network=pn)
    phys_air = OpenPNM.Physics.Standard(network=pn,
                                        phase=air,
                                        pores=pn.pores(),
                                        throats=pn.throats())

    BC1_pores = pn.pores(labels=['left'])
    BC2_pores = pn.pores(labels=['right'])

    alg_1 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
    alg_1.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=1,
                                  pores=BC1_pores)
    alg_1.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=0,
                                  pores=BC2_pores)
    alg_1.run(iterative_solver='gmres')

    alg_2 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
    alg_2.set_boundary_conditions(bctype='Neumann',
                                  bcvalue=-1e-11,
                                  pores=BC1_pores)
    alg_2.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=0,
                                  pores=BC2_pores)
    alg_2.run(iterative_solver='cg')

    alg_3 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
    alg_3.set_boundary_conditions(bctype='Neumann_group',
                                  bcvalue=-3e-10,
                                  pores=BC1_pores)
    alg_3.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=0,
                                  pores=BC2_pores)
    alg_3.run()

    alg_4 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
    alg_4.set_boundary_conditions(bctype='Neumann_group',
                                  bcvalue=-3e-10,
                                  pores=BC1_pores)
    alg_4.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=0,
                                  pores=BC2_pores)
    alg_4.setup()
    alg_4.solve()

    assert round(sp.absolute(alg_1.rate(BC1_pores))[0], 16) == round(sp.absolute(alg_1.rate(BC2_pores))[0], 16)
    assert round(sp.absolute(alg_2.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_2['pore.'+air.name+'_bcval_Neumann']))[0]*len(BC1_pores), 16)
    assert round(sp.absolute(alg_3.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_3['pore.'+air.name+'_bcval_Neumann_group']))[0], 16)
    assert round(sp.absolute(alg_4.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_4['pore.'+air.name+'_bcval_Neumann_group']))[0], 16)

    assert round(sp.absolute(sp.sum(alg_1.rate(BC1_pores,mode='single'))),16) == round(sp.absolute(alg_1.rate(BC1_pores))[0],16)
    assert round(sp.absolute(sp.sum(alg_2.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_2.rate(BC2_pores))[0],16)
    assert round(sp.absolute(sp.sum(alg_3.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_3.rate(BC2_pores))[0],16)
    assert round(sp.absolute(sp.sum(alg_4.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_4.rate(BC2_pores))[0],16)

 def _do_outer_iteration_stage(self):
     #Generate curve from points
     for inv_val in self._inv_points:
         #Apply one applied pressure and determine invaded pores
         logger.info('Applying capillary pressure: '+str(inv_val))
         self._do_one_inner_iteration(inv_val)
     #Store results using networks' get/set method
     self['pore.inv_Pc'] = self._p_inv
     self['throat.inv_Pc'] = self._t_inv
     #Find invasion sequence values (to correspond with IP algorithm)
     self._p_seq = sp.searchsorted(sp.unique(self._p_inv),self._p_inv)
     self._t_seq = sp.searchsorted(sp.unique(self._t_inv),self._t_inv)
     self['pore.inv_seq'] = self._p_seq
     self['throat.inv_seq'] = self._t_seq
     #Calculate Saturations
     v_total = sp.sum(self._net['pore.volume'])+sp.sum(self._net['throat.volume'])
     sat = 0.
     self['pore.inv_sat'] = 1.
     self['throat.inv_sat'] = 1.
     for i in range(self._npts):
         inv_pores = sp.where(self._p_seq==i)[0]
         inv_throats = sp.where(self._t_seq==i)[0]
         new_sat = (sum(self._net['pore.volume'][inv_pores])+sum(self._net['throat.volume'][inv_throats]))/v_total
         sat += new_sat
         self['pore.inv_sat'][inv_pores] = sat
         self['throat.inv_sat'][inv_throats] = sat

def __compute_affiliation(label_image, mask_image, bounding_boxes):
    """
    Computes which regions of the supplied label_image belong to the mask_image's foreground
    respectively background. When a region belongs to both, it is assigned to the foreground
    if more voxels belong to the foreground than in the background and vice-versa.
    In the case of equal affiliation, the region is assigned to the background.
    @return fg_ids, bg_ids
    """
    # simple extraction
    fg_ids = list(scipy.unique(label_image[mask_image]))
    bg_ids = list(scipy.unique(label_image[~mask_image]))
    # decide for overlapping regions whether they are 50 or more in fg or in bg
    for rid in set(fg_ids) & set(bg_ids):
        relevant_region_label_image = label_image[bounding_boxes[rid - 1]]
        relevant_region_mask_image = mask_image[bounding_boxes[rid - 1]]
        fg_part = 0
        bg_part = 0
        for affiliation, rid2 in zip(relevant_region_mask_image.ravel(), relevant_region_label_image.ravel()):
            if rid2 == rid:
                if affiliation: fg_part += 1
                else: bg_part += 1
        #fg_part = relevant_region_label_image[relevant_region_mask_image]
        #bg_part = relevant_region_label_image[~relevant_region_mask_image]
        if fg_part > bg_part: # if more voxels of region rid in fg than in bg
            bg_ids.remove(rid)
        else:
            fg_ids.remove(rid)
    # debug line, can be removed if the above code is final
    if 0 != len(set(fg_ids) & set(bg_ids)): raise Exception('Error making fg and bg ground truth distinct.') 
    return fg_ids, bg_ids

def fit_dispersion(counts, disp_raw, disp_conv, sf, CFG, dmatrix1):

    mean_count = sp.mean(counts / sf, axis=1)[:, sp.newaxis]
    index = sp.where(disp_conv)[0]

    lowerBound = sp.percentile(sp.unique(disp_raw[index]), 1)
    upperBound = sp.percentile(sp.unique(disp_raw[index]), 99)

    idx = sp.where((disp_raw > lowerBound) & (disp_raw < upperBound))[0]

    matrix = sp.ones((idx.shape[0], 2), dtype='float')
    matrix[:, 0] /= mean_count[idx].ravel()

    modGamma = sm.GLM(disp_raw[idx], matrix, family=sm.families.Gamma(sm.families.links.identity))
    res = modGamma.fit()
    Lambda = res.params

    disp_fitted = disp_raw.copy()
    ok_idx = sp.where(~sp.isnan(disp_fitted))[0]
    disp_fitted[ok_idx] = Lambda[0] / mean_count[ok_idx] + Lambda[1]

    if sp.sum(disp_fitted > 0) > 0:
        print "Found dispersion fit"

    if CFG['diagnose_plots']:
        plot.mean_variance_plot(counts=counts,
                                disp=disp_fitted,
                                matrix=dmatrix1,
                                figtitle='Fitted Dispersion Estimate',
                                filename=os.path.join(CFG['plot_dir'], 'dispersion_fitted.pdf'),
                                CFG=CFG)

    return (disp_fitted, Lambda, idx)

 def plot_setup(p):
     pylab.ylabel("Throughput gain [\%]")
     pylab.xscale('log', basex=2)
     pylab.xticks(
         list(scipy.unique(group['symbols'])),
         list(scipy.unique(group['symbols'])))
     plotter.set_markers(p)
     plotter.set_slave_info(slavename)

    def find_neighbor_throats(self,pnums,flatten=True,mode='union'):
        r"""
        Returns a list of throats neighboring the given pore(s)

        Parameters
        ----------
        pnums : array_like
            Indices of pores whose neighbors are sought
        flatten : boolean, optional
            If flatten is True (default) a 1D array of unique throat ID numbers
            is returned. If flatten is False the returned array contains arrays
            of neighboring throat ID numbers for each input pore, in the order
            they were sent.
        mode : string, optional
            Specifies which neighbors should be returned.  The options are: 
            
            * 'union' : All neighbors of the input pores

            * 'intersection' : Only neighbors shared by all input pores 
            
            * 'not_intersection' : Only neighbors not shared by any input pores

        Returns
        -------
        neighborTs : 1D array (if flatten is True) or ndarray of arrays (if
            flatten if False)

        Examples
        --------
        >>> pn = OpenPNM.Network.Cubic(name='doc_test').generate(divisions=[5,5,5],lattice_spacing=[1])
        >>> pn.find_neighbor_throats(pnums=[0,1])
        array([0, 1, 2, 3, 4, 5])
        >>> pn.find_neighbor_throats(pnums=[0,1],flatten=False)
        array([array([0, 1, 2]), array([0, 3, 4, 5])], dtype=object)
        """
        #Test for existance of incidence matrix
        try:
            neighborTs = self.incidence_matrix['lil']['connections'].rows[[pnums]]
        except:
            self._logger.info('Creating incidence matrix, please wait')
            self.create_incidence_matrix()
            neighborTs = self.incidence_matrix['lil']['connections'].rows[[pnums]]
        if flatten:
            #All the empty lists must be removed to maintain data type after hstack (numpy bug?)
            neighborTs = [sp.asarray(x) for x in neighborTs if x]
            neighborTs = sp.hstack(neighborTs)
            #Remove references to input pores and duplicates
            if mode == 'not_intersection':
                neighborTs = sp.unique(sp.where(sp.bincount(neighborTs)==1)[0])
            elif mode == 'union':
                neighborTs = sp.unique(neighborTs)
            elif mode == 'intersection':
                neighborTs = sp.unique(sp.where(sp.bincount(neighborTs)>1)[0])
        else:
            for i in range(0,sp.size(pnums)):
                neighborTs[i] = sp.array(neighborTs[i])
        return sp.array(neighborTs,ndmin=1)

 def plot_setup(p):
     pylab.ylabel("Throughput" + " [" + list(group['unit'])[0] + "]")
     pylab.yscale('log')
     pylab.xscale('log', basex=2)
     pylab.xticks(
         list(scipy.unique(group['symbols'])),
         list(scipy.unique(group['symbols'])))
     plotter.set_markers(p)
     plotter.set_legend_columns(3)

 def test_no_late_filling(self):
     mip = op.algorithms.Porosimetry(network=self.net)
     mip.setup(phase=self.hg)
     mip.set_inlets(pores=self.net.pores('left'))
     mip.run()
     assert len(sp.unique(mip['pore.invasion_pressure'])) > 1
     assert len(sp.unique(mip['pore.invasion_sequence'])) > 1
     assert len(sp.unique(mip['throat.invasion_pressure'])) > 1
     assert len(sp.unique(mip['throat.invasion_sequence'])) > 1

    def evaluate_trapping(self, p_outlets):
        r"""
        Finds trapped pores and throats after a full ordinary
        percolation simulation has been run.

        Parameters
        ----------
        p_outlets : array_like
            A list of pores that define the wetting phase outlets.
            Disconnection from these outlets results in trapping.

        Returns
        -------
        It creates arrays called ``pore.trapped`` and ``throat.trapped``, but
        also adjusts the ``pore.inv_Pc`` and ``throat.inv_Pc`` arrays to set
        trapped locations to have infinite invasion pressure.

        """
        self['pore.trapped'] = sp.zeros([self.Np, ], dtype=float)
        self['throat.trapped'] = sp.zeros([self.Nt, ], dtype=float)
        try:
            # Get points used in OP
            inv_points = sp.unique(self['pore.inv_Pc'])
        except:
            raise Exception('Orindary percolation has not been run!')
        tind = self._net.throats()
        conns = self._net.find_connected_pores(tind)
        for inv_val in inv_points[0:-1]:
            # Find clusters of defender pores
            Pinvaded = self['pore.inv_Pc'] <= inv_val
            Cstate = sp.sum(Pinvaded[conns], axis=1)
            Tinvaded = self['throat.inv_Pc'] <= inv_val
            # 0 = all open, 1=1 pore filled,
            # 2=2 pores filled 3=2 pores + 1 throat filled
            Cstate = Cstate + Tinvaded
            clusters = self._net.find_clusters(Cstate == 0)
            # Clean up clusters (invaded = -1, defended >=0)
            clusters = clusters * (~Pinvaded) - (Pinvaded)
            # Identify clusters connected to outlet sites
            out_clusters = sp.unique(clusters[p_outlets])
            trapped_pores = ~sp.in1d(clusters, out_clusters)
            trapped_pores[Pinvaded] = False
            if sum(trapped_pores) > 0:
                inds = (self['pore.trapped'] == 0) * trapped_pores
                self['pore.trapped'][inds] = inv_val
                trapped_throats = self._net.find_neighbor_throats(trapped_pores)
                trapped_throat_array = np.asarray([False] * len(Cstate))
                trapped_throat_array[trapped_throats] = True
                inds = (self['throat.trapped'] == 0) * trapped_throat_array
                self['throat.trapped'][inds] = inv_val
                inds = (self['throat.trapped'] == 0) * (Cstate == 2)
                self['throat.trapped'][inds] = inv_val
        self['pore.trapped'][self['pore.trapped'] > 0] = sp.inf
        self['throat.trapped'][self['throat.trapped'] > 0] = sp.inf
        self['pore.inv_Pc'][self['pore.trapped'] > 0] = sp.inf
        self['throat.inv_Pc'][self['throat.trapped'] > 0] = sp.inf

def centre_of_mass(geometry, vertices='throat.offset_vertices', **kwargs):
    r"""
    Calculate the centre of mass of the throat from the voronoi vertices.
    """
    Nt = geometry.num_throats()
    outer_verts = geometry['throat.vertices']
    offset_verts = geometry[vertices]
    normal = geometry['throat.normal']
    z_axis = [0, 0, 1]
    value = _sp.ndarray([Nt, 3])
    for i in range(Nt):
        if len(offset_verts[i]) > 2:
            verts = offset_verts[i]
        elif len(outer_verts[i]) > 2:
            verts = outer_verts[i]
        else:
            verts = []
        if len(verts) > 0:
            # For boundaries some facets will already be aligned with the axis -
            # if this is the case a rotation is unnecessary and could also cause
            # problems
            angle = tr.angle_between_vectors(normal[i], z_axis)
            if angle == 0.0 or angle == _sp.pi:
                "We are already aligned"
                rotate_input = False
                facet = verts
            else:
                rotate_input = True
                M = tr.rotation_matrix(tr.angle_between_vectors(normal[i], z_axis),
                                       tr.vector_product(normal[i], z_axis))
                facet = _sp.dot(verts, M[:3, :3].T)
            # Now we have a rotated facet aligned with the z axis - make 2D
            facet_2D = _sp.column_stack((facet[:, 0], facet[:, 1]))
            z = _sp.unique(_sp.around(facet[:, 2], 10))
            if len(z) == 1:
                # We need the vertices arranged in order so perform a convex hull
                hull = ConvexHull(facet_2D)
                ordered_facet_2D = facet_2D[hull.vertices]
                # Call the routine to calculate an area wighted centroid from the
                # 2D polygon
                COM_2D = vo.PolyWeightedCentroid2D(ordered_facet_2D)
                COM_3D = _sp.hstack((COM_2D, z))
                # If we performed a rotation we need to rotate back
                if (rotate_input):
                    MI = tr.inverse_matrix(M)
                    # Unrotate the offset coordinates using the inverse of the
                    # original rotation matrix
                    value[i] = _sp.dot(COM_3D, MI[:3, :3].T)
                else:
                    value[i] = COM_3D
            else:
                print('Rotation Failed: ' + str(_sp.unique(facet[:, 2])))

    return value

    def get_scan_IF_inds(self, scan_ind, IF_ind) :
        """Gets the record indices of the fits file that correspond to the
        given scan and IF.

        Note that the scans are numbered with 0 corresponding to the first scan
        in the file i.e., it is not the session scan number."""

        # TODO: Should check valid scan IF, and raise value errors as apropriate
        thescan = self.scan_set[scan_ind]
        theIF = self.IF_set[IF_ind]
        
        # Find all the records that correspond to this IF and this scan.
        # These indicies *should now be ordered in time, cal (on off)
        # and in polarization, once the IF is isolated.
        (inds_sif,) = sp.where(sp.logical_and(self._IFs_all==theIF, 
                                        self._scans_all==thescan))
        ncal = len(sp.unique(self.fitsdata.field('CAL')[inds_sif]))
        npol = len(sp.unique(self.fitsdata.field('CRVAL4')[inds_sif]))
        
        # Reform to organize by pol, cal, etc.
        ntimes = len(inds_sif)//npol//ncal
        inds_sif = sp.reshape(inds_sif, (ntimes, npol, ncal))

        if self.verify_ordering > 0:
            # We expect noise cal to be on for every second record.
            for thecal in range(ncal) :
                tmp = sp.unique(self.fitsdata.field('CAL')[inds_sif[:,:,thecal]])
                if len(tmp) > 1 :
                    raise ce.DataError("Calibration (ON/OFF) not in "
                                    "perfect order in file: "+self.fname)
            # Polarization should cycle through 4 modes (-5,-7,-8,-6)
            for thepol in range(npol) :
                tmp = sp.unique(self.fitsdata.field('CRVAL4')
                            [inds_sif[:,thepol,:]])
                if len(tmp) > 1 :
                    raise ce.DataError("Polarizations not in perfect order in "
                                    "file: "+self.fname)
            # We expect the entries to be sorted in time and for time to not
            # change across pol and cal.
            lastLST = 0
            for ii in range(ntimes) :
                # Sometimes won't have the LST.
                try :
                    thisLST = self.fitsdata.field('LST')[inds_sif[ii,0,0]]
                # If 'LST' is missing raises a KeyError in later versions of
                # pyfits, and a NameError in earlier ones.
                except (KeyError, NameError) :
                    break
                if not (sp.allclose(self.fitsdata.field('LST')
                        [inds_sif[ii,:,:]] - thisLST, 0)) :
                    raise ce.DataError("LST change across cal or pol in "
                                       "file: " + self.fname)

        return inds_sif

    def run(self, npts=25, inv_pressures=None):
        r"""
        Run the algorithm for specified number of points or at given capillary
        pressures.

        Parameters
        ----------
        npts : scalar
            The number of points to obtain on the curve.  The points are
            automatically selected to span the range of capillary pressures
            using a logarithmic spacing (more points are lower capillary
            pressure values).

        inv_pressures : array_like
            A list of capillary pressures to apply. List should contain
            increasing and unique values.
        """
        # If no invasion points are given then generate some
        if inv_pressures is None:
            logger.info('Generating list of invasion pressures')
            min_p = sp.amin(self['throat.entry_pressure']) * 0.98  # nudge down
            max_p = sp.amax(self['throat.entry_pressure']) * 1.02  # bump up
            inv_points = sp.logspace(sp.log10(min_p),
                                     sp.log10(max_p),
                                     npts)
        else:
            # Make sure the given invastion points are sensible
            inv_points = sp.unique(inv_pressures)
        self._inv_points = inv_points

        # Ensure inlets are set
        if sp.sum(self['pore.inlets']) == 0:
            raise Exception('Inlet pores have not been specified')

        # Ensure outlet pores are set if trapping is enabled
        if self._trapping:
            if sp.sum(self['pore.outlets']) == 0:
                raise Exception('Outlet pores have not been specified')

        # Generate curve from points
        for inv_val in self._inv_points:
            # Apply one applied pressure and determine invaded pores
            logger.info('Applying capillary pressure: ' + str(inv_val))
            self._apply_percolation(inv_val)
            if self._trapping:
                logger.info('Checking for trapping')
                self._check_trapping(inv_val)

        # Find invasion sequence values (to correspond with IP algorithm)
        Pinv = self['pore.inv_Pc']
        self['pore.inv_seq'] = sp.searchsorted(sp.unique(Pinv), Pinv)
        Tinv = self['throat.inv_Pc']
        self['throat.inv_seq'] = sp.searchsorted(sp.unique(Tinv), Tinv)

def make_plot_twoclass(X,Y,W,kernel):
    fig = plt.figure(figsize=(5,4))
    fig.clf()
    colors = "brymcwg"

    # Plot the decision boundary.
    h = .2 # stepsize in mesh
    x_min, x_max = X[0,:].min() - 1, X[0,:].max() + 1
    y_min, y_max = X[1,:].min() - 1, X[1,:].max() + 1
    xx, yy = sp.meshgrid(sp.arange(x_min, x_max, h),
            sp.arange(y_min, y_max, h))

    Z = predict_svm_kernel(sp.c_[sp.ones(xx.ravel().shape[-1]), xx.ravel(), yy.ravel()].T,sp.vstack((sp.ones((1,X.shape[-1])),X)),W,kernel).reshape(xx.shape)
    cs = plt.contourf(xx, yy, Z,alpha=.5)
    plt.axis('tight')
    plt.colorbar()
    plt.axis('equal')
    y = sp.maximum(0,-Y)+1
    # plot the data
    plt.hold(True)

    ypred = 	W.T.dot(kernel[0](X,X,kernel[1]).T)
    for ic in sp.unique(y):
        idx = (y == int(ic)).flatten()
        sv = (Y.flatten()[idx]*ypred[idx] < 1)
        plt.plot(X[0,idx.nonzero()[0][sv]], X[1,idx.nonzero()[0][sv]], colors[int(ic)]+'o',markersize=13)
        plt.plot(X[0,idx.nonzero()[0][~sv]], X[1,idx.nonzero()[0][~sv]], colors[int(ic)]+'o',markersize=7)
    plt.axis('tight')

    plt.xlabel('$X_1$')
    plt.ylabel('$X_2$')

    #plt.title('SVM, Accuracy=%0.2f'%(Y==sp.sign(ypred)).mean())

    #plt.show()
    plt.savefig('./svm_kernel.pdf')

    fig = plt.figure(figsize=(5,5))
    fig.clf()
    colors = "brymcwg"
    for ic in sp.unique(y):
        idx = (y == int(ic)).flatten()
        plt.plot(X[0,idx], X[1,idx], colors[int(ic)]+'o',markersize=8)
    plt.axis('tight')

    plt.xlabel('$X_1$')
    plt.ylabel('$X_2$')
    plt.xlim((x_min,x_max))
    plt.ylim((y_min,y_max))
    plt.grid()
    #plt.show()
    plt.savefig('./svm_kernel_xor_data.pdf')