def est_condprob2(data, val, given):
    """Calculate the probability of P(X|Y,Z)

    est_condprob2(data, 'A', ['M', 'LC'])"""

    if not isinstance(given, list):
        raise IndexError("Given must be a list or tuple of givens")
    elif len(given) != 2:
        raise IndexError("I need multiple givens!  Give me more...give me more!")

    gcols = []
    for g in given:
        if g in ['M', 'F']:
            gcols.append(1)
        elif g in ['LC', 'SC', 'T']:
            gcols.append(2)
        elif g in ['A', 'B', 'C']:
            gcols.append(0)

    if val in ['M', 'F']:
        vcol = 1
    elif val in ['LC', 'SC', 'T']:
        vcol = 2
    elif val in ['A', 'B', 'C']:
        vcol = 0

    datsize = data.shape[0]
    needed = [val, given[0], given[1]]
    t = sp.where([sp.all(data[i]==needed) for i in range(datsize)])[0]

    t2 = sp.where([sp.all(data[i,1:]==given) for i in range(datsize)])[0]
    
    return float(t.size)/t2.size

    def cross_validation(self, x, y, v=5, sig_r=2.0 ** sp.arange(-8, 0), mu_r=10.0 ** sp.arange(-15, 0)):
        # Get parameters
        n = x.shape[0]
        ns = sig_r.size
        nm = mu_r.size
        err = sp.zeros((ns, nm))

        # Initialization of the indices for the cross validation
        cv = CV()
        cv.split_data_class(y, v=v)

        for i in range(ns):
            for j in range(nm):
                for k in range(v):
                    model_temp = KDA()
                    model_temp.train(x[cv.it[k], :], y[cv.it[k]], sig=sig_r[i], mu=mu_r[j])
                    yp = model_temp.predict(x[cv.iT[k], :], x[cv.it[k], :], y[cv.it[k]])
                    yp.shape = y[cv.iT[k]].shape
                    t = sp.where(yp != y[cv.iT[k]])[0]
                    err[i, j] += float(t.size) / yp.size
                    del model_temp
        err /= v
        t = sp.where(err == err.min())
        self.sig = sig_r[t[0][0]]
        self.mu = mu_r[t[1][0]]
        return sig_r[t[0][0]], mu_r[t[1][0]], err

def compute_ndvi(im, r=0, ir=1, NODATA=-10000):
    """The function computes the NDVI of a multivalued image. It checks if
    there is NODATA value or division per zeros.
    
    Args:
    im: the image to process
    r: the number of the band that corresponds to the red band.
    ir: the number of the band that corresponds to the infra-red band.
    NODATA: the value of the NODATA
    
    Returns:
    ndvi =  the ndvi of the image
    """

    ## Get the size fo the image
    [nl, nc, nb] = im.shape

    ## Be sure that we can do 'floating operation'
    imf = im.astype(sp.float64)
    ndvi = sp.empty((nl, nc))

    if nb < 2:
        print "Two bands are needed to compute the NDVI"
        return None
    else:
        den = imf[:, :, ir - 1] + imf[:, :, r - 1]  # Pre compute the denominator
        t = sp.where((den > 0) & (imf[:, :, 1] != NODATA))
        ndvi[t] = (imf[t[0], t[1], ir - 1] - imf[t[0], t[1], r - 1]) / den[t]  # compute the ndvi on the safe samples

        if len(t[0]) < nl * nc:
            t = sp.where((den == 0) | (imf[:, :, 1] == NODATA))  # check for problematic pixels
            ndvi[t] = NODATA

        imf = []
        return ndvi

 def _generate_masked_mesh(self, cell_mask=None):
     r"""
     Generates the mesh based on the cell mask provided
     """
     #
     if cell_mask is None:
         cell_mask = sp.ones(self.data_map.shape, dtype=bool)
     #
     # initializing arrays
     self._edges = sp.ones(0, dtype=str)
     self._merge_patch_pairs = sp.ones(0, dtype=str)
     self._create_blocks(cell_mask)
     #
     # building face arrays
     mapper = sp.ravel(sp.array(cell_mask, dtype=int))
     mapper[mapper == 1] = sp.arange(sp.count_nonzero(mapper))
     mapper = sp.reshape(mapper, (self.nz, self.nx))
     mapper[~cell_mask] = -sp.iinfo(int).max
     #
     boundary_dict = {
         'bottom':
             {'bottom': mapper[0, :][cell_mask[0, :]]},
         'top':
             {'top': mapper[-1, :][cell_mask[-1, :]]},
         'left':
             {'left': mapper[:, 0][cell_mask[:, 0]]},
         'right':
             {'right': mapper[:, -1][cell_mask[:, -1]]},
         'front':
             {'front': mapper[cell_mask]},
         'back':
             {'back': mapper[cell_mask]},
         'internal':
             {'bottom': [], 'top': [], 'left': [], 'right': []}
     }
     #
     # determining cells linked to a masked cell
     cell_mask = sp.where(~sp.ravel(cell_mask))[0]
     inds = sp.in1d(self._field._cell_interfaces, cell_mask)
     inds = sp.reshape(inds, (len(self._field._cell_interfaces), 2))
     inds = inds[:, 0].astype(int) + inds[:, 1].astype(int)
     inds = (inds == 1)
     links = self._field._cell_interfaces[inds]
     #
     # adjusting order so masked cells are all on links[:, 1]
     swap = sp.in1d(links[:, 0], cell_mask)
     links[swap] = links[swap, ::-1]
     #
     # setting side based on index difference
     sides = sp.ndarray(len(links), dtype='<U6')
     sides[sp.where(links[:, 1] == links[:, 0]-self.nx)[0]] = 'bottom'
     sides[sp.where(links[:, 1] == links[:, 0]+self.nx)[0]] = 'top'
     sides[sp.where(links[:, 1] == links[:, 0]-1)[0]] = 'left'
     sides[sp.where(links[:, 1] == links[:, 0]+1)[0]] = 'right'
     #
     # adding each block to the internal face dictionary
     inds = sp.ravel(mapper)[links[:, 0]]
     for side, block_id in zip(sides, inds):
         boundary_dict['internal'][side].append(block_id)
     self.set_boundary_patches(boundary_dict, reset=True)

 def simulate(self,x0,lambd):
     Dt = self.param['Dt']
     # Dt needs to be a multiple of param['Dt']
     dt = self.param['dt']
     D = lambd[1]
     a = lambd[0]
     N = self.param['N']
     drift = self.param['drift']
     x = scipy.array(x0)
     
     tstart = 0
     tcur = tstart
     while (tcur < tstart + Dt + dt/2 ):
         tcur += dt
         # the random number
         dW=self.rand.normal(loc=0.,scale=scipy.sqrt(2*D*dt),size=N)
       #  if tcur == dt:    #only print random number for first time step
       #      print 'dW =',  dW            
         
         # the process
         drift_term = a * drift(x)
         x=x+drift_term*dt+dW
         # and reflecting boundary conditions
         scipy.where(x>self.domain[1],2*self.domain[1]-x,x)
         scipy.where(x<self.domain[0],2*self.domain[0]-x,x)
     return x

    def populate_out_of_dip_theta(self, n, dip):

        out_of_dip = asarray(self.populate_distribution(
            self.out_of_dip_theta_dist, n))

        (errorIndexes,) = where((out_of_dip > (175 - dip)) &
                                (out_of_dip < (185 - dip)))

        if len(errorIndexes) > 0:
            for i in errorIndexes:
                blnBadNum = True
                count = 0
                while blnBadNum:
                    newNum = self.populate_distribution(
                        self.out_of_dip_theta_dist, 1)
                    if ((newNum[0] <= (175 - dip)) | (newNum[0] >= (185 - dip))):
                        blnBadNum = False
                    count = count + 1
                    if count > 1000:
                        msg = "Bad out of dip theta range in fault \
                                     source file"
                        raise IOError(msg)
                out_of_dip[i] = newNum[0]
        (errorIndexes,) = where((out_of_dip > (175 - dip)) &
                                (out_of_dip < (185 - dip)))

        if len(errorIndexes) > 0:
            msg = "Bad out of dip theta range in fault \
                                     source file"
            raise IOError(msg)

        return out_of_dip

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 getEncodedData(filename,encoding="additive",phenotype_id=None,maf=0.0):
    f = h5py.File(filename,'r')
    phenotype_id = str(phenotype_id)
    if not phenotype_id==None:
        sample_ids = f['Genotype/sample_ids'][:]
        p_sample_ids = f['Phenotypes'][phenotype_id]['sample_ids'][:]
        y = f['Phenotypes'][phenotype_id]['y'][:]
        ind = sp.where(~sp.isnan(y))[0]
        y = y[ind]
        p_sample_ids = p_sample_ids[ind]
        ind = (sp.reshape(sample_ids,(sample_ids.shape[0],1))==p_sample_ids).nonzero()
        raw = f['Genotype/raw'][:]
        raw = raw[ind[0],:]
        [encoded,maf_v] = encodeHeterozygousData(raw)
        ind = sp.where(maf_v>=maf)[0]
        encoded = encoded[:,ind]
        identifiers = f['Genotype/identifiers'][:]
        identifiers = identifiers[ind]
        maf_v = maf_v[ind]
        f.close()
        return [encoded,maf_v,identifiers]
    if encoding=="additive":
        if 'encoded_additive' in f['Genotype'].keys():
            encoded = f['Genotype/encoded_additive'][:]
            maf_v = f['Genotype/global_maf'][:]
        else:
            [encoded,maf_v] = encodeHeterozygousData(f['Genotype/raw'][:])
    identifiers = f['Genotype/identifiers'][:]
    f.close()
    return [encoded,maf_v,identifiers]

    def newEpisode(self):
        if self.learning:
            params = ravel(self.explorationlayer.module.params)
            target = ravel(sum(self.history.getSequence(self.history.getNumSequences()-1)[2]) / 500)
        
            if target != 0.0:
                self.gp.addSample(params, target)
                if len(self.gp.trainx) > 20:
                    self.gp.trainx = self.gp.trainx[-20:, :]
                    self.gp.trainy = self.gp.trainy[-20:]
                    self.gp.noise = self.gp.noise[-20:]
                    
                self.gp._calculate()
                        
                # get new parameters where mean was highest
                max_cov = diag(self.gp.pred_cov).max()
                indices = where(diag(self.gp.pred_cov) == max_cov)[0]
                pick = indices[random.randint(len(indices))]
                new_param = self.gp.testx[pick]
            
                # check if that one exists already in gp training set
                if len(where(self.gp.trainx == new_param)[0]) > 0:
                    # add some normal noise to it
                    new_param += random.normal(0, 1, len(new_param))

                self.explorationlayer.module._setParameters(new_param)

            else:
                self.explorationlayer.drawRandomWeights()
        
        # don't call StateDependentAgent.newEpisode() because it randomizes the params
        LearningAgent.newEpisode(self)

 def cryptoInternal(self, data, base):
     addresses = scipy.array(range(base, base + (len(data) * 2), 2), scipy.uint32)
     for mask, xorVal in self.XOR_TABLE1:
         data = scipy.where((addresses & mask) == mask, data ^ xorVal, data)
     for mask, xorVal in self.XOR_TABLE2:
         data = scipy.where((addresses & mask) != 0,    data ^ xorVal, data)
     return data

 def smart_threshold(self):
     self.median = numpy.median(self.data)
     self.std = numpy.std(self.data)
     blank = scipy.where(self.data < self.median+0.25*self.std)
     signal = scipy.where(self.data > self.median+0.25*self.std)
     self.data[blank] = 0.0
     self.data[signal] = 1.0

 def maskLowStddVoxels(self, dds, nMeanDds, nStddDds):
     unique = np.unique(sp.where(nStddDds.subd.asarray() <= 1.0/3, dds.subd.asarray(), dds.mtype.maskValue()))
     unique = unique[sp.where(unique != dds.mtype.maskValue())]
     if (dds.mpi.comm != None):
         unique = dds.mpi.comm.allreduce(unique.tolist(), op=mpi.SUM)
         unique = np.unique(unique)
     rootLogger.info("Unique constant stdd values = %s" % (unique,))
     rootLogger.info("Creating mask from unique constant values...")
     mskDds = mango.zeros_like(dds, mtype="segmented")
     for uVal in unique:
         mskDds.asarray()[...] = sp.where(dds.asarray() == uVal, 1, mskDds.asarray())
     rootLogger.info("Done creating mask from unique constant values.")
 
     rootLogger.info("Labeling connected constant zero-stdd regions...")
     mskDds.updateHaloRegions()
     mskDds.mirrorOuterLayersToBorder(False)
     self.writeIntermediateDds("_000ZeroStddForLabeling", mskDds)
     lblDds = mango.image.label(mskDds, 1)
     rootLogger.info("Done labeling connected constant stdd regions.")
     self.writeIntermediateDds("_000ZeroStdd", lblDds)
     
     countThresh = 0.01 * sp.product(lblDds.shape)
     rootLogger.info("Eliminating large clusters...")
     lblDds = mango.image.eliminate_labels_by_size(lblDds, minsz=int(countThresh), val=lblDds.mtype.maskValue())
     self.writeIntermediateDds("_000ZeroStddLargeEliminated", lblDds)
 
     rootLogger.info("Assigning mask values...")
     mskDds.subd.asarray()[...] = \
         sp.where(lblDds.subd.asarray() == lblDds.mtype.maskValue(), True, False)
     self.writeIntermediateDds("_000ZeroStddMskd", mskDds)
     del lblDds
     for tmpDds in [dds, nMeanDds, nStddDds]:
         tmpDds.subd.asarray()[...] = \
             sp.where(mskDds.subd.asarray(), tmpDds.mtype.maskValue(), tmpDds.subd.asarray())

 def eliminatePercentileTails(self, mskDds, loPercentile=10.0, hiPercentile=90.0):
     """
     Trims lower and/or upper image histogram tails by replacing :samp:`mskDds`
     voxel values with :samp:`mskDds.mtype.maskValue()`. 
     """
     rootLogger.info("Eliminating percentile tails...")
     rootLogger.info("Calculating element frequencies...")
     elems, counts = elemfreq(mskDds)
     rootLogger.info("elems:\n%s" % (elems,))
     rootLogger.info("counts:\n%s" % (counts,))
     cumSumCounts = sp.cumsum(counts, dtype="float64")
     percentiles = 100.0*(cumSumCounts/float(cumSumCounts[-1]))
     percentileElems = elems[sp.where(sp.logical_and(percentiles > loPercentile, percentiles < hiPercentile))]
     loThresh = percentileElems[0]
     hiThresh = percentileElems[-1]
     rootLogger.info("Masking percentiles range (%s,%s) = (%s,%s)" % (loPercentile, hiPercentile, loThresh, hiThresh))
     mskDds.asarray()[...] = \
         sp.where(
             sp.logical_and(
                 sp.logical_and(mskDds.asarray() >= loThresh, mskDds.asarray() <= hiThresh),
                 mskDds.asarray() != mskDds.mtype.maskValue()
             ),
             mskDds.asarray(),
             mskDds.mtype.maskValue()
         )
     rootLogger.info("Done eliminating percentile tails.")

 def step(self, *args):
     """First update the step size, then actually take a step along the gradient."""
     g = self.model.grad(*args);
     
     # Update the weighted Exponential sq avg.
     self.sqExpAvgGrad *= self.exponentAvgM;
     self.sqExpAvgGrad += (1-self.exponentAvgM) * g**2;
     self.sqExpAvgGrad[:] = where(self.sqExpAvgGrad < EPSILON, EPSILON, self.sqExpAvgGrad);
     
     # Uodate the muVect
     possUpdate = 1 + self.qLearningRate * g * self.expAvgGrad / self.sqExpAvgGrad
     #log.debug('max(possUpdate): %.4f,  min(possUpdate): %.4f' % (max(possUpdate), min(possUpdate)))
     ## Keep this from going negative.
     possUpdate = where(possUpdate < 0.001, 0.001, possUpdate);
     self.muVect *= possUpdate
     
     # Do something to cap the update rate.  This is allowing the step rate to overpower the decay completely
     self.muVect = where(self.muVect > self.maxMuVect, self.maxMuVect, self.muVect);
     
     # Then update the exponential average
     self.expAvgGrad *= self.exponentAvgM;
     self.expAvgGrad += (1-self.exponentAvgM) * g;
     
     self.model.params -= self.muVect * g
     Trainer.step(self,*args)

 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 spectralSlope(wl, flux, dFlux, wlStart, wlStop, beta_guess, **kwargs):
    bm = scipy.where( (wl > wlStart) & (wl < wlStop) & numpy.isfinite(flux) )[0]

    if ( 'strongLines' in kwargs ):
       for line, width in zip(kwargs['strongLines'], kwargs['lineWidths']):
           new_bm = scipy.where( abs(wl[bm]-line) > width)
           bm = bm[new_bm[0]]

    x = wl[bm]
    y = flux[bm]
    dy = dFlux[bm]

    normalization = y[0]
    z = normalization*(x/wlStart)**beta_guess

    coeffs = [normalization, beta_guess]
    
    fitfunc = lambda p, x : p[0]*(x/wlStart)**(p[1])
    errfunc = lambda p, x, z, dz: numpy.abs((fitfunc(p, x) - z)/dz)
    pfit = scipy.optimize.leastsq(errfunc, coeffs, args=(numpy.asarray(x, dtype=numpy.float64),
    numpy.asarray(y,dtype=numpy.float64), numpy.asarray(dy,dtype=numpy.float64)), full_output = 1)

    if ( 'plt' in kwargs ):
        original = Gnuplot.Data(x, y, with_='lines')
        guess = Gnuplot.Data(x, z, with_='lines')
        new = Gnuplot.Data(x, pfit[0][0]*(x/wlStart)**(pfit[0][1]), with_='lines')
        kwargs['plt'].plot(original, guess, new)
        #raw_input()

    return pfit[0]

  def scanSound(self, source, minnotel):
    binarized = source
    scale = 60. / self.wavetempo * (binarized[0].size / self.duration)
    noise_length = scale*minnotel

    antinoised = sp.zeros_like(binarized)

    for i in range(sp.shape(binarized)[0]):
      new_line = binarized[i, :].copy()
      diffed = sp.diff(new_line)
      ones_keys = sp.where(diffed == 1)[0]
      minus_keys = sp.where(diffed == -1)[0]
      
      if(ones_keys.size != 0 and minus_keys.size != 0):
        if(ones_keys[0] > minus_keys[0]):
          new_line = self.cutNoise(
              (0, minus_keys[0]), noise_length, new_line)
          minus_keys = sp.delete(minus_keys, 0)

        if(ones_keys[-1] > minus_keys[-1]):
          new_line = self.cutNoise(
              (ones_keys[-1], new_line.size-1), noise_length, new_line)
          ones_keys = sp.delete(ones_keys, -1)

        for j in range(sp.size(ones_keys)):
          new_line = self.cutNoise(
              (ones_keys[j], minus_keys[j]), noise_length, new_line)

        antinoised[i, :] = new_line

    return antinoised

def binSyntheticSpectrum(spectrum, native_wl, new_wl):
    """
        This routine pixelates a synthetic spectrum, in effect simulating the 
        discrete nature of detector pixels.
    """
    retval = numpy.zeros(len(new_wl))
    for i in range(len(new_wl)-1):
        bm = scipy.where( (native_wl > new_wl[i]) & (
            native_wl <= new_wl[i+1]))[0]
        if (len(bm) > 1):
            num=scipy.integrate.simps(spectrum[bm], x=native_wl[bm])
            denom = max(native_wl[bm]) - min(native_wl[bm])
            retval[i] = num/denom
        elif (len(bm) == 1):
            retval[i] = 0.0#native_wl[bm]
        else:
            retval[i] = 0.0#retval[-1]

    bm = scipy.where(native_wl > new_wl[-1])[0]
    if len(bm) > 1:
        num = scipy.integrate.simps(spectrum[bm], x=native_wl[bm])
        denom = max(native_wl[bm]) - min(native_wl[bm])
        retval[-1] = num/denom
    else:
        if len(bm) == 1:
            retval[-1] = spectrum[bm]
        else:
            retval[-1] = spectrum[-1]

    return retval

def get_concentration_functions(composition_table_dict):

    meta = composition_table_dict['meta']
    composition_table = Table.from_dict(composition_table_dict['data'])
    elements = [col for col in composition_table.columns if col not in meta]
    x = composition_table["X"].values
    y = composition_table["Y"].values
    cats = composition_table["X"].unique()
    concentration, conc, d, y_c, functions = {}, {}, {}, {}, RecursiveDict()

    for el in elements:
        concentration[el] = to_numeric(composition_table[el].values)/100.
        conc[el], d[el], y_c[el] = {}, {}, {}

        if meta['X'] == 'category':
            for i in cats:
                k = '{:06.2f}'.format(float(i))
                y_c[el][k] = to_numeric(y[where(x==i)])
                conc[el][k] = to_numeric(concentration[el][where(x==i)])
                d[el][k] = interp1d(y_c[el][k], conc[el][k])

            functions[el] = lambda a, b, el=el: d[el][a](b)

        else:
            functions[el] = interp2d(float(x), float(y), concentration[el])

    return functions

def make_introns_feasible(introns, genes, CFG):
# introns = make_introns_feasible(introns, genes, CFG)

    tmp1 = sp.array([x.shape[0] for x in introns[:, 0]])
    tmp2 = sp.array([x.shape[0] for x in introns[:, 1]])
    
    unfeas = sp.where((tmp1 > 200) | (tmp2 > 200))[0]
    print >> CFG['fd_log'], 'found %i unfeasible genes' % unfeas.shape[0]

    while unfeas.shape[0] > 0:
        ### make filter more stringent
        CFG['read_filter']['exon_len'] = min(36, CFG['read_filter']['exon_len'] + 4)
        CFG['read_filter']['mincount'] = 2 * CFG['read_filter']['mincount']
        CFG['read_filter']['mismatch'] = max(CFG['read_filter']['mismatch'] - 1, 0)

        ### get new intron counts
        tmp_introns = get_intron_list(genes[unfeas], CFG)
        introns[unfeas, :] = tmp_introns

        ### stil unfeasible?
        tmp1 = sp.array([x.shape[0] for x in introns[:, 0]])
        tmp2 = sp.array([x.shape[0] for x in introns[:, 1]])

        still_unfeas = sp.where((tmp1 > 200) | (tmp2 > 200))[0]
        idx = sp.where(~sp.in1d(unfeas, still_unfeas))[0]

        for i in unfeas[idx]:
            print >> CFG['fd_log'], '[feasibility] set criteria for gene %s to: min_ex %i, min_conf %i, max_mism %i' % (genes[i].name, CFG['read_filter']['exon_len'], CFG['read_filter']['mincount'], CFG['read_filter']['mismatch'])
        unfeas = still_unfeas;

    return introns