def getGenoIndex(self,pos_start=None,pos_end=None,chrom=None,windowsize=0):
        """computes 0-based genotype index from position of cumulative position.
        Positions can be given in one out of two ways:
        - position (pos_start-pos_end on chrom)
        - cumulative position (pos_cum_start-pos_cum_end)
        If all these are None (default), then all genotypes are returned

        Args:
            pos_start:       position based selection (start position)
            pos_end:       position based selection (end position)
            chrom:      position based selection (chromosome)

        Returns:
            idx_start:         genotype index based selection (start index)
            idx_end:         genotype index based selection (end index)
        """
        if (pos_start is not None) & (pos_end is not None):
            assert pos_start[0] == pos_end[0], "chromosomes don't match between start and end position"
            I = self.position["chrom"]==pos_start[0]
            I = I & (self.position["pos"]>=(pos_start[1]-windowsize)) & (self.position["pos"]<(pos_end[1]+windowsize))
            I =  sp.nonzero(I)[0]
            idx_start = I.min()
            idx_end = I.max()
        elif chrom is not None:
            I = self.position["chrom"]==chrom
            I =  sp.nonzero(I)[0]
            if I.size==0:
                return None
            idx_start = I.min()
            idx_end = I.max()
        else:
            idx_start=None
            idx_end=None
        return idx_start,idx_end

def graph_connected_components(graph_mat):
    """takes graph as matrix and return list of connected components
          
       arguments:
       graph_mat - matrix of graph
              
       output:
       list_of_comp - list of components, list_of_comp[i] - list of node numbers in (i-1)-th component"""
    component_of_graph = scipy.zeros(graph_mat.shape[0],dtype = scipy.int8) # component_of_graph[i] is the number of component of i-th node.
    cur_comp = 1 #the number of current component (see below)
    try:
        tmp_nodes_to_process = [scipy.nonzero(component_of_graph==0)[0][0]] #node indexes to process
    except IndexError:
            #exceptional situation, when graph_mat is empty
            return []

    #kind of breadth first search
    while(len(tmp_nodes_to_process)>0): #while there is nodes to process
        cur_node = tmp_nodes_to_process.pop() #take one from array
        lnodes_numbers = scipy.nonzero(graph_mat[cur_node,:])[0] #take indexes of all of it linked nodes
        #and choose that corresponds to the non processed nodes of them, the node is non processed if its component is zero
        lnodes_numbers = scipy.extract(component_of_graph[lnodes_numbers] == 0,lnodes_numbers)
        tmp_nodes_to_process +=lnodes_numbers.tolist()
        component_of_graph[lnodes_numbers] = cur_comp
        # if there is no linked nodes, start processing of new connected component, and next unprocessed node
        if (len(tmp_nodes_to_process) == 0):
            cur_comp+=1
            tmp_arr = scipy.nonzero(component_of_graph==0)[0]
            if (len(tmp_arr)>0):tmp_nodes_to_process = [tmp_arr[0]] 

    list_of_comp = [] #collect list
    for i in range(cur_comp+1):
        tmp_arr=scipy.nonzero(component_of_graph==(i+1))[0].tolist()
        if (len(tmp_arr)>0):list_of_comp+=[tmp_arr]
    return list_of_comp

 def cut(self):
   average = sp.sum(sp.absolute(self.data))/sp.size(self.data)
   head = sp.nonzero(sp.absolute(self.data)>average)[0][5]
   bottom = sp.nonzero(sp.absolute(self.data)>average)[0][-1]
   self.data = self.data[head:bottom]
   self.duration_list = self.duration_list[head:bottom]
   self.duration = self.duration_list[-1] - self.duration_list[0]

 def setup_Q(self):
     self.Q = scipy.zeros((self.nperiods, self.ndists, self.ndists))
     #D = self.D * self.Dmask
     for p in range(self.nperiods):
         for i, d1 in enumerate(self.dists):
             s1 = sum(d1)
             if s1 > 0:
                 for j, d2 in enumerate(self.dists):
                     s2 = sum(d2)
                     xor = scipy.logical_xor(d1, d2)
                     # only consider transitions between dists that are
                     # 1 step (dispersal or extinction) apart
                     if sum(xor) == 1:
                         dest = scipy.nonzero(xor)[0]
                         #prior = self.dist_priors[i]
                         if s1 < s2: # dispersal
                             rate = 0.0
                             for src in scipy.nonzero(d1)[0]:
                                 rate += self.D[p,src,dest] * \
                                         self.Dmask[p,src,dest]
                             # for each area in d1, add rate of
                             # dispersal to dest
                         else: # extinction
                             rate = self.E[p,dest]
                         #self.Q[i][j] = (prior * rate)
                         self.Q[p,i,j] = rate
         self.set_Qdiag(p)

    def getGenoIndex(self,pos0=None,pos1=None,chrom=None,pos_cum0=None,pos_cum1=None):
        """computes 0-based genotype index from position of cumulative position. 
        Positions can be given in one out of two ways: 
        - position (pos0-pos1 on chrom)
        - cumulative position (pos_cum0-pos_cum1)
        If all these are None (default), then all genotypes are returned

        Args:
            pos0:       position based selection (start position)
            pos1:       position based selection (stop position)
            chrom:      position based selection (chromosome)
            pos_cum0:   cumulative position based selection (start position)
            pos_cum1:   cumulative position based selection (stop position)
        
        Returns:
            i0:         genotype index based selection (start index)
            i1:         genotype index based selection (stop index)
        """
        if (pos0 is not None) & (pos1 is not None) & (chrom is not None):
            I = self.gneoChrom==chrom
            I = I & (self.genoPos>=p0) & (self.genoPos<p1)
            I = SP.nonzero(I)[0]
            i0 = I.min()
            i1 = I.max()
        elif (pos_cum0 is not None) & (pos_cum1 is not None):
            I = (self.genoPos_cum>=pos_cum0) & (self.genoPos_cum<pos_cum1)
            I = SP.nonzero(I)[0]
            if I.size==0:
                return None
            i0 = I.min()
            i1 = I.max()
        else:
            i0=None
            i1=None
        return i0,i1

def __interpolateBetweenBinaryObjects(obj1, obj2, slices):
    """
    Takes two binary objects and puts slices slices in-between them, each of which
    contains a smooth binary transition between the objects.
    @note private inner function
    """
    if not obj1.shape == obj2.shape:
        raise AttributeError(
            "The two supplied objects have to be of the same shape, not {} and {}.".format(obj1.shape, obj2.shape)
        )

    # constant
    offset = 0.5  # must be a value smaller than the minimal distance possible
    temporal_dimension = 3

    # get all voxel position
    obj1_voxel = scipy.nonzero(obj1)
    obj2_voxel = scipy.nonzero(obj2)

    # get smallest pairwise distances between all object voxels
    distances = cdist(scipy.transpose(obj1_voxel), scipy.transpose(obj2_voxel))

    # keep for each True voxel of obj1 only the smallest distance to a True voxel in obj2
    min_distances = distances.min(1)

    # test if all seems to work
    if len(min_distances) != len(obj1_voxel[0]):
        raise Exception("Invalid number of minimal distances received.")

    # replace True voxels in obj1 with their respective distances to the True voxels in obj2
    thr_obj = obj1.copy()
    thr_obj = thr_obj.astype(scipy.float_)
    thr_obj[obj1_voxel] = min_distances
    thr_obj[obj1_voxel] += offset  # previous steps distances include zeros, therefore this is required

    # compute the step size for each slice that is added
    maximum = min_distances.max()
    step = maximum / float(slices + 1)
    threshold = maximum

    # control step: see if thr_obj really corresponds to obj1
    if not scipy.all(thr_obj.astype(scipy.bool_) == obj1.astype(scipy.bool_)):
        raise Exception("First created object does not correspond to obj1.")

    # assemble return volume
    return_volume = [thr_obj.astype(scipy.bool_)]  # corresponds to obj1
    for _ in range(slices):
        threshold -= step
        # remove all value higher than the threshold
        thr_obj[thr_obj > threshold] = 0
        # add binary volume to list /makes a copy)
        return_volume.append(thr_obj.astype(scipy.bool_))

    # add last slice (corresponds to es obj2 slice)
    thr_obj[thr_obj > offset] = 0
    return_volume.append(thr_obj.astype(scipy.bool_))

    # return binary scipy array
    return scipy.rollaxis(scipy.asarray(return_volume, dtype=scipy.bool_), 0, temporal_dimension + 1)

def move_nodes(d, lfun, nodes):#, dely):
    " Move nodes one timestep (projecting lost points to the boundary) "

    # get delauney 
    dely = get_dely(d, nodes)

    # force constant 
    dt = .1520015
    deps = h0 * sqrt(finfo(float64).eps)
    restlength_factor = 1.400025
    
    # bars and midpoints
    bars, barmids = dely[-2:]
    barvecs = nodes[:,bars[:,1]] - nodes[:,bars[:,0]]
    barls = sqrt((barvecs**2).sum(0))
    u = barvecs / barls       # unit vectors
    
    # force from each bar
    restlen = restlength_factor * lfun(*barmids)

    # print('restlen: {0} \nbarls : {1}'.format(restlen.shape, barls.shape))
    logic = restlen > barls
    f = where(logic, restlen-barls, 0.0)
    # f = where(f<h0/2.0, f, h0/2.0)
    # 
    ns = nodes.shape
    spmat = sparse.csc_matrix
    # print(ns)
    # print(u[0].shape)
    # print(f.shape)
    # print(bars[:,0].shape)

    dp =  (-spmat((u[0]*f, (0*bars[:,0], bars[:,0])), shape=ns).todense() +
            spmat((u[0]*f, (0*bars[:,1], bars[:,1])), shape=ns).todense())
    dp += (-spmat((u[1]*f, (0*bars[:,0]+1, bars[:,0])), shape=ns).todense() +
            spmat((u[1]*f, (0*bars[:,1]+1, bars[:,1])), shape=ns).todense())

    nodes = array(nodes + dt*dp)

    # project boundary points back into the domain
    d_ = d(*nodes)
    ix = nonzero(d_>0)
    some_out = True
    count = 0
    while some_out:
        gradx = 1.0/deps * (d(nodes[0,ix]+deps, nodes[1,ix]) - d_[ix])
        grady = 1.0/deps * (d(nodes[0,ix], nodes[1,ix] + deps) - d_[ix])
        norm = sqrt(gradx**2 + grady**2)
        nodes[0,ix] -= d_[ix]*gradx / norm
        nodes[1,ix] -= d_[ix]*grady / norm
        d_ = d(*nodes)
        ix = nonzero(d_>geps)
        some_out = ix[0].size
        count+=1
        if count>5: #raise ValueError("counted "+str(ix[0].size)+" nodes oob")
            print("counted ",str(ix[0].size)," nodes oob")
            break
    return nodes

 def __call__(self, dx):
     """
     hat function evaluation
     """
     x = convarg(dx) - self.c
     y = sc.zeros_like(x)
     i = sc.nonzero((x>self.a)&(x<=.0))
     y[i] = 1. - x[i]/self.a
     i = sc.nonzero((x>.0)&(x<self.b))
     y[i] = 1. - x[i]/self.b
     return(y)

def gen_single_trial(interval_lengths, rates):
    """ Generates a single spike train with intervals of length
    `interval_lengths` and the firing rates given in `rates`.
    """
    boundaries = sp.ones(len(interval_lengths) + 1) * pq.s
    boundaries[1:] = [l.rescale(boundaries.units) for l in interval_lengths]
    rates = rates[sp.nonzero(boundaries[1:])]
    boundaries = boundaries[sp.nonzero(boundaries)]
    boundaries[0] = 0.0 * pq.s
    boundaries = sp.cumsum(boundaries)
    return stools.st_concatenate([stg.gen_homogeneous_poisson(
        rate, t_start=boundaries[i], t_stop=boundaries[i + 1])
        for i, rate in enumerate(rates)])

def subgraph_with_center(graph_mat, node_names, center_name):
    """takes number of center node, and returns new subgraph, that consists of all nodes connected with center.
        
        arguments:
        graph_mat,node_names - graph description as matrix (see matrix_from_tuple_list doc for details)
        center_num - name of the center node

        output:
        subgraph,sub_node_names - subgraph description as matrix (see matrix_from_tuple_list doc for details)"""
    center_num = scipy.nonzero(node_names==center_name)[0][0]
    center_friends_num = scipy.nonzero(graph_mat[center_num,:])[0] #indexes of nodes that linked with  central node including itself
    subgraph = graph_mat[center_friends_num,:][:,center_friends_num] # FIXME we consider part of graph which consists of nodes linked with center only 
    sub_node_names = node_names[center_friends_num]
    return (subgraph,sub_node_names)

def dup_fig(x, y):
    colors = [None, 'black','blue','brown','purple','orange','cyan','gray','yellow','black','red','green']
    k = len(sp.unique(y))
    for i in range(1, k+1):
        inds = sp.nonzero(y == i)
        plt.plot(x[inds, 0], x[inds, 1], 'o', color=colors[i])
    plt.show()

    def getPhenotypes(self,phenotype_IDs=None,phenotype_query=None,sample_idx=None,center=True,intersection=False):
        """load Phenotypes

        Args:
            phenotype_IDs:      names of phenotypes to load
            phenotype_query:    string hoding a pandas query (e.g. "(environment==1) & (phenotype_ID=='growth')"
                                selects all phenotypes that have a phenotype_ID equal to growth under environment 1.
            sample_idx:         Boolean sample index for subsetting individuals
            center:             Boolean: mean center (and mean-fill in missing values if intersection==False)? (default True)
            impute:             imputation of missing values (default: True)
            intersection:       restrict observation to those obseved in all phenotypes? (default: False)

        Returns:
            phenotypes:     [N x P] scipy.array of phenotype values for P phenotypes
            sample_idx_intersect:        index of individuals in phenotypes after filtering missing valuesS
        """
        if phenotype_IDs is not None:
            I = SP.array([SP.nonzero(self.phenotype_ID==n)[0][0] for n in phenotype_IDs])
        elif phenotype_query is not None:
            try:
                I = self.index_frame.query(phenotype_query).values[:,0]

                #if there are no results we won't actually get an exception, we just get an
                #empty response
                if len(I) == 0:
                    print "query '%s' yielded no results!" % (phenotype_query)
                    I = SP.zeros([0],dtype="int")

            except Exception, arg:

                print "query '%s' yielded no results: %s" % (phenotype_query, str(arg))

                I = SP.zeros([0],dtype="int")

	def __init__(self, grid, fArray, zDrawsSorted):
		assert(len(grid) == len(fArray))
		(self.grid, self.fArray) = (grid, fArray)
		self.zDraws = zDraws		
		self.slopes = scipy.zeros(len(grid) - 1)
		self.dx = grid[1] - grid[0]
		for i in range(len(grid) - 1):
			self.slopes[i] = (fArray[i+1] - fArray[i]) / self.dx
		# set up sums
		self.cellSums = scipy.zeros(len(grid) + 1)
		self.boundaryIndices = [len(zDraws)] * len(grid)
		for (i, x) in enumerate(grid):
			indices = scipy.nonzero(self.zDraws >= x)[0]
			if (len(indices) > 0):
				self.boundaryIndices[i] = indices[0]
		self.cellSums[0] = scipy.sum(self.zDraws[0:self.boundaryIndices[0]])
		for i in range(1, len(self.cellSums)-1):
			self.cellSums[i] = scipy.sum(self.zDraws[self.boundaryIndices[i-1] : self.boundaryIndices[i]])
		self.cellSums[-1] = scipy.sum(self.zDraws[self.boundaryIndices[-1] : ])
		
		diff = scipy.sum(self.zDraws) - scipy.sum(self.cellSums)
		print("diff: %f" % diff)
		for i in range(len(grid)):
			if (self.boundaryIndices[i] < len(self.zDraws)):
				print("grid point %f, boundary %f" % (self.grid[i], self.zDraws[self.boundaryIndices[i]]))
			else:
				print("grid point %f, no draws to right" % self.grid[i])

    def errorApproximation(self, ratio, dim=20):

        self.buildMatrix()

        sumNonzeros = (self.vxm !=0).sum()
        numTest = int(ratio*sumNonzeros)

        elementList = []

        nonZeroTuple = sp.nonzero(self.vxm)

        for x in range(int(numTest)):
            rInt = sp.random.randint(0,nonZeroTuple[0].size)
            randrow = nonZeroTuple[0][rInt]
            randcolumn = nonZeroTuple[1][rInt]

            valElementIndex = [randrow,randcolumn]
            elementList.append(valElementIndex)

        self.modvxm = sp.copy(self.vxm)

        for x in elementList:
            self.modvxm[x[0],x[1]] = 0

        self.modvmx = self.fillAverages(vxm = self.modvxm)
        self.newmodvxm = self.predict(dim,vxm=self.modvxm)

        sqDiff = 0
        for x in elementList:
            sqDiff += sp.square(self.newmodvxm[x[0],x[1]] - self.vxm[x[0],x[1]])
        self.rmse = sp.sqrt(sqDiff/len(elementList))

    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 _do_one_inner_iteration(self,inv_val):
     #Written by: Jeff Gostick ([email protected])
     r"""
     Determines which throats are invaded at a given applied capillary pressure
     
     This function uses the scipy.csgraph module for the graph theory cluster
     labeling algorithm (connected_components)
     
     Dependencies:
         - 
     Creates:
         - 
     """
     #Generate a tlist containing boolean values for throat state
     self._net.throat_properties['invaded'] = self._net.throat_properties['Pc_entry']<inv_val
     #Fill adjacency matrix with invasion state info
     self._net.create_adjacency_matrix('invaded',sprsfmt='csr',dropzeros=True)
     clusters = sprs.csgraph.connected_components(self._net._adjmatrix_csr)[1]
     #Clean up (not invaded = -2, invaded >=0)
     clusters = (clusters[0:]>=0)*(clusters[0:]+1)
     #Identify clusters connected to invasion sites
     if self._ALOP == 1:
         inj_clusters = self._inv_src*clusters
     elif self._OP == 1:
         temp1 = self._net.throat_properties['invaded']*((self._net.throat_properties['connections'][:,0]+1)-1)
         temp2 = self._net.throat_properties['invaded']*((self._net.throat_properties['connections'][:,1]+1)-1)
         inj_clusters = np.append(self._net.pore_properties['numbering'][temp1[temp1>=0]],self._net.pore_properties['numbering'][temp2[temp2>=0]])
     #Trim non-connected clusters
     temp = sp.unique(inj_clusters[sp.nonzero(inj_clusters)])
     inv_clusters = sp.zeros([np.size(clusters,0)],np.int32)
     for i in range(0,np.shape(temp)[0]):
         pores=sp.where(clusters==temp[i])[0]
         inv_clusters[pores] = temp[i]
     return(inv_clusters)

def hzd_do_value(sa, r_nu, rtrn_rte):
    """
    parrams:
    sa       [vector (nx1)] response spectral accelerations
    r_nu     [vector (nx1)] event activity for the corresponding element
    in sa
    rtrn_rte  [vector (mx1)] return rates of interest.

    returns:
    hzd       [vector (1xm)] hazard value for each return rate
    """
    # Get rid of events with sa = 0, since they will effect the end of the
    # curve

    assert sa.shape == r_nu.shape, str(
        sa.shape) + 'should = ' + str(r_nu.shape)
    nonzero_ind = nonzero(sa)
    sa = sa[nonzero_ind]
    r_nu = r_nu[nonzero_ind]

    hzd, cumnu = _rte2cumrte(sa, r_nu)
    # annual exceedance rate = cumulative event activity
    # for exceedance rates larger than what we have data for, give 0.
    # for exceedance rates smaller than what we have data for, give hzd[0].
    if len(hzd) == 0:
        hzd_val = zeros(rtrn_rte.shape)
    else:
        hzd_val = interp(rtrn_rte, cumnu, hzd, left=hzd[0], right=0.0)
    return hzd_val

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 interpolateBetweenBinaryObjects(obj1, obj2, slices):
    """
    Takes two binary objects and puts slices slices in-between them, each of which
    contains a smooth binary transition between the objects.
    """
    # constants
    temporal_dimension = 3
    # flip second returned binary objects along temporal axis
    slicer = [slice(None) for _ in range(obj1.ndim + 1)]
    slicer[temporal_dimension] = slice(None, None, -1)
    # logical-and combination
    ret = (
        __interpolateBetweenBinaryObjects(obj1, obj2, slices)
        | __interpolateBetweenBinaryObjects(obj2, obj1, slices)[slicer]
    )
    # control step: see if last volume corresponds to obj2
    slicer[temporal_dimension] = slice(-1, None)
    if not scipy.all(scipy.squeeze(ret[slicer]) == obj2.astype(scipy.bool_)):
        raise Exception(
            "Last created object does not correspond to obj2. Difference of {} voxels.".format(
                len(scipy.nonzero(scipy.squeeze(ret[slicer]) & obj2.astype(scipy.bool_))[0])
            )
        )

    return ret

    def getRegion(self,size=3e4,min_nSNPs=1,chrom_i=None,pos_min=None,pos_max=None):
        """
        Sample a region from the piece of genotype X, chrom, pos
        minSNPnum:  minimum number of SNPs contained in the region
        Ichrom:  restrict X to chromosome Ichrom before taking the region
        cis:        bool vector that marks the sorted region
        region:  vector that contains chrom and init and final position of the region
        """
        if (self.chrom is None) or (self.pos is None):
            bim = plink_reader.readBIM(self.bfile,usecols=(0,1,2,3))
            chrom = SP.array(bim[:,0],dtype=int)
            pos   = SP.array(bim[:,3],dtype=int)
        else:
            chrom = self.chrom
            pos   = self.pos

        if chrom_i is None:
            n_chroms = chrom.max()
            chrom_i  = int(SP.ceil(SP.rand()*n_chroms))

        pos   = pos[chrom==chrom_i]
        chrom = chrom[chrom==chrom_i]

        ipos = SP.ones(len(pos),dtype=bool)
        if pos_min is not None:
            ipos = SP.logical_and(ipos,pos_min<pos)

        if pos_max is not None:
            ipos = SP.logical_and(ipos,pos<pos_max)

        pos = pos[ipos]
        chrom = chrom[ipos]

        if size==1:
            # select single SNP
            idx = int(SP.ceil(pos.shape[0]*SP.rand()))
            cis  = SP.arange(pos.shape[0])==idx
            region = SP.array([chrom_i,pos[idx],pos[idx]])
        else:
            while 1:
                idx = int(SP.floor(pos.shape[0]*SP.rand()))
                posT1 = pos[idx]
                posT2 = pos[idx]+size
                if posT2<=pos.max():
                    cis = chrom==chrom_i
                    cis*= (pos>posT1)*(pos<posT2)
                    if cis.sum()>min_nSNPs: break
            region = SP.array([chrom_i,posT1,posT2])

        start = SP.nonzero(cis)[0].min()
        nSNPs  = cis.sum()

        if self.X is None:
            rv = plink_reader.readBED(self.bfile,useMAFencoding=True,start = start, nSNPs = nSNPs,bim=bim)
            Xr = rv['snps']
        else:
            Xr = self.X[:,start:start+nSnps]

        return Xr, region