Calculate enrichment based on hypergeometric distribution
179416
def calc_motif_enrichment(sample, background, mtc=None, len_sample=None, len_back=None):
    """"""
    
    INF = "Inf"


    if mtc not in [None, "Bonferroni", "Benjamini-Hochberg", "None"]:
        raise RuntimeError("Unknown correction: %s" % mtc)

    sig = {}
    p_value  = {}
    n_sample = {}
    n_back = {}
    
    if not(len_sample):
        len_sample = sample.seqn()
    if not(len_back):
        len_back = background.seqn()

    for motif in sample.motifs.keys():
        p = "NA"
        s = "NA"
        q = len(sample.motifs[motif])
        m = 0
        if(background.motifs.get(motif)):
            m = len(background.motifs[motif])
            n = len_back - m
            k = len_sample
            p = phyper(q - 1, m, n, k) 
            if p != 0:
                s = -(log(p)/log(10))
            else:
                s = INF
        else:
            s = INF
            p = 0.0

        sig[motif] = s
        p_value[motif] = p
        n_sample[motif] = q
        n_back[motif] = m
    
    if mtc == "Bonferroni":
        for motif in p_value.keys():
            if  p_value[motif] != "NA":
                p_value[motif] = p_value[motif] * len(p_value.keys())
                if p_value[motif] > 1:
                    p_value[motif] = 1
    elif mtc == "Benjamini-Hochberg":
        motifs = sorted(p_value.keys(), key=lambda x: -p_value[x])
        l = len(p_value)
        c = l
        for m in motifs:
            if  p_value[m] != "NA":
                p_value[m] = p_value[m] * l / c 
            c -= 1

    return (sig, p_value, n_sample, n_back)
539281
function bfrac(a, b, x, y) {
      var lambda;

      lambda = a - (a + b) * x;
      return brcomp(a, b, x, y) + Math.log(contFrac(
         function(i) { // betas from DiDonato
            var n;

            n = i - 1;
            return i === 0 ? 0
                           : n + n * (b - n) * x / (a + 2 * n - 1) +
                             (a + n) / (a + 2 * n + 1) *
                                (lambda + 1 + n * (1 + y));
         },
         function(i) { // alphas from DiDonato
            var n;

            n = i - 1;
            return i === 1 ? 1
                           : (a + n - 1) * (a + b + n - 1) /
                             (a + 2 * n - 1) / (a + 2 * n - 1) *
                             n * (b - n) * x * x;
         }
      ));
   }
141545
double[][] computeStatsFillModel(LloydsIterationTask task, KMeansModel model, final Vec[] vecs, final double[] means, final double[] mults, final int[] modes, int k) {
      // Fill in the model based on original destandardized centers
      if (model._parms._standardize) {
        model._output._centers_std_raw = task._cMeans;
      }
      model._output._centers_raw = destandardize(task._cMeans, _isCats, means, mults);
      model._output._size = task._size;
      model._output._withinss = task._cSqr;
      double ssq = 0;       // sum squared error
      for( int i=0; i<k; i++ )
        ssq += model._output._withinss[i]; // sum squared error all clusters
      model._output._tot_withinss = ssq;

      // Sum-of-square distance from grand mean
      if(k == 1) {
        model._output._totss = model._output._tot_withinss;
      }
      else {
        // If data already standardized, grand mean is just the origin
        TotSS totss = new TotSS(means,mults,modes, train().domains(), train().cardinality()).doAll(vecs);
        model._output._totss = totss._tss;
      }
      model._output._betweenss = model._output._totss - model._output._tot_withinss;  // MSE between-cluster
      model._output._iterations++;

      model._output._history_withinss = ArrayUtils.copyAndFillOf( model._output._history_withinss, model._output._history_withinss.length+1, model._output._tot_withinss);
      model._output._k = ArrayUtils.copyAndFillOf(model._output._k, model._output._k.length+1, k);
      model._output._training_time_ms = ArrayUtils.copyAndFillOf(model._output._training_time_ms, model._output._training_time_ms.length+1, System.currentTimeMillis());
      model._output._reassigned_count = ArrayUtils.copyAndFillOf(model._output._reassigned_count, model._output._reassigned_count.length+1, task._reassigned_count);

      // Two small TwoDimTables - cheap
      model._output._model_summary = createModelSummaryTable(model._output);
      model._output._scoring_history = createScoringHistoryTable(model._output);

      // Take the cluster stats from the model, and assemble them into a model metrics object
      model._output._training_metrics = makeTrainingMetrics(model);

      return task._cMeans;      // New centers
    }
775238
def compute_enriched_threshold(average_window_readcount):
    # type: (float) -> int
    

    current_threshold, survival_function = 0, 1
    for current_threshold in count(start=0, step=1):
        survival_function -= poisson.pmf(current_threshold,
                                         average_window_readcount)
        if survival_function <= WINDOW_P_VALUE:
            break

    island_enriched_threshold = current_threshold + 1

    return island_enriched_threshold
13506
def collective_dr_squared( self ):
        
        return sum( np.square( sum( [ atom.dr for atom in self.atoms ] ) ) )
1008457
def _get_split_rhat(values, round_to=2):
    """"""
    shape = values.shape
    if len(shape) != 2:
        raise TypeError("Effective sample size calculation requires 2 dimensional arrays.")
    _, num_samples = shape
    num_split = num_samples // 2
    # Calculate split chain mean
    split_chain_mean1 = np.mean(values[:, :num_split], axis=1)
    split_chain_mean2 = np.mean(values[:, num_split:], axis=1)
    split_chain_mean = np.concatenate((split_chain_mean1, split_chain_mean2))
    # Calculate split chain variance
    split_chain_var1 = np.var(values[:, :num_split], axis=1, ddof=1)
    split_chain_var2 = np.var(values[:, num_split:], axis=1, ddof=1)
    split_chain_var = np.concatenate((split_chain_var1, split_chain_var2))
    # Calculate between-chain variance
    between_chain_variance = num_samples / 2 * np.var(split_chain_mean, ddof=1)
    # Calculate within-chain variance
    within_chain_variance = np.mean(split_chain_var)
    # Estimate of marginal posterior variance
    split_rhat = np.sqrt(
        (between_chain_variance / within_chain_variance + num_samples / 2 - 1) / (num_samples / 2)
    )

    return round(split_rhat, round_to)
325526
public List<IntegerSolution> doCrossover(
          double probability, IntegerSolution parent1, IntegerSolution parent2) {
    List<IntegerSolution> offspring = new ArrayList<IntegerSolution>(2);

    offspring.add((IntegerSolution) parent1.copy()) ;
    offspring.add((IntegerSolution) parent2.copy()) ;

    int i;
    double rand;
    double y1, y2, yL, yu;
    double c1, c2;
    double alpha, beta, betaq;
    int valueX1, valueX2;

    if (randomGenerator.getRandomValue() <= probability) {
      for (i = 0; i < parent1.getNumberOfVariables(); i++) {
        valueX1 = parent1.getVariableValue(i);
        valueX2 = parent2.getVariableValue(i);
        if (randomGenerator.getRandomValue() <= 0.5) {
          if (Math.abs(valueX1 - valueX2) > EPS) {

            if (valueX1 < valueX2) {
              y1 = valueX1;
              y2 = valueX2;
            } else {
              y1 = valueX2;
              y2 = valueX1;
            }

            yL = parent1.getLowerBound(i);
            yu = parent1.getUpperBound(i);
            rand = randomGenerator.getRandomValue();
            beta = 1.0 + (2.0 * (y1 - yL) / (y2 - y1));
            alpha = 2.0 - Math.pow(beta, -(distributionIndex + 1.0));

            if (rand <= (1.0 / alpha)) {
              betaq = Math.pow((rand * alpha), (1.0 / (distributionIndex + 1.0)));
            } else {
              betaq = Math
                .pow(1.0 / (2.0 - rand * alpha), 1.0 / (distributionIndex + 1.0));
            }

            c1 = 0.5 * ((y1 + y2) - betaq * (y2 - y1));
            beta = 1.0 + (2.0 * (yu - y2) / (y2 - y1));
            alpha = 2.0 - Math.pow(beta, -(distributionIndex + 1.0));

            if (rand <= (1.0 / alpha)) {
              betaq = Math.pow((rand * alpha), (1.0 / (distributionIndex + 1.0)));
            } else {
              betaq = Math
                .pow(1.0 / (2.0 - rand * alpha), 1.0 / (distributionIndex + 1.0));
            }

            c2 = 0.5 * (y1 + y2 + betaq * (y2 - y1));

            if (c1 < yL) {
              c1 = yL;
            }

            if (c2 < yL) {
              c2 = yL;
            }

            if (c1 > yu) {
              c1 = yu;
            }

            if (c2 > yu) {
              c2 = yu;
            }

            if (randomGenerator.getRandomValue() <= 0.5) {
              offspring.get(0).setVariableValue(i, (int)c2);
              offspring.get(1).setVariableValue(i, (int)c1);
            } else {
              offspring.get(0).setVariableValue(i, (int)c1);
              offspring.get(1).setVariableValue(i, (int)c2);
            }
          } else {
            offspring.get(0).setVariableValue(i, valueX1);
            offspring.get(1).setVariableValue(i, valueX2);
          }
        } else {
          offspring.get(0).setVariableValue(i, valueX2);
          offspring.get(1).setVariableValue(i, valueX1);
        }
      }
    }

    return offspring;
  }
804646
def convenience_calc_fisher_approx(self, params):
        
        shapes, intercepts, betas = self.convenience_split_params(params)

        args = [betas,
                self.design,
                self.alt_id_vector,
                self.rows_to_obs,
                self.rows_to_alts,
                self.choice_vector,
                self.utility_transform,
                self.calc_dh_d_shape,
                self.calc_dh_dv,
                self.calc_dh_d_alpha,
                intercepts,
                shapes,
                self.ridge,
                self.weights]

        return cc.calc_fisher_info_matrix(*args)
540252
private void doMutation(double probability, DoubleSolution solution) {
    double rnd, deltaq, tempDelta;
    double y, yl, yu;

    for (int i = 0; i < solution.getNumberOfVariables(); i++) {
      if (randomGenerator.nextDouble() <= probability) {
        y = solution.getVariableValue(i);
        yl = solution.getLowerBound(i) ;
        yu = solution.getUpperBound(i) ;
        rnd = randomGenerator.nextDouble();
          
        tempDelta = Math.pow(rnd, -delta);
        deltaq = 0.5 * (rnd - 0.5) * (1 - tempDelta);
          
        y = y + deltaq * (yu - yl);
        y = solutionRepair.repairSolutionVariableValue(y, yl, yu);
        solution.setVariableValue(i, y);
      }
    }
  }
59872
private double computeExplainedVariance(double[] eigenValues, int filteredEigenPairs) {
    double strongsum = 0., weaksum = 0.;
    for(int i = 0; i < filteredEigenPairs; i++) {
      strongsum += eigenValues[i];
    }
    for(int i = filteredEigenPairs; i < eigenValues.length; i++) {
      weaksum += eigenValues[i];
    }
    return strongsum / (strongsum + weaksum);
  }
128016
def generateAmplitudes(k):
    """"""
    res = np.array([np.random.random() for i in range(k)])
    return res/res.sum()
547818
def rmultivariate_hypergeometric(n, m, size=None):
    

    N = len(m)
    urn = np.repeat(np.arange(N), m)

    if size:
        draw = np.array([[urn[i] for i in np.random.permutation(len(urn))[:n]]
                         for j in range(size)])

        r = [[np.sum(draw[j] == i) for i in range(len(m))]
             for j in range(size)]
    else:
        draw = np.array([urn[i] for i in np.random.permutation(len(urn))[:n]])

        r = [np.sum(draw == i) for i in range(len(m))]
    return np.asarray(r)
1133423
private int sampleFullConditional(int m, int n) {

        // remove z_i from the count variables
        int topic = z[m][n];
        word_topic_matrix[documents[m][n]][topic]--;
        nd[m][topic]--;
        nwsum[topic]--;
        ndsum[m]--;

        // do multinomial sampling via cumulative method:
        float[] p = new float[K];
        for (int k = 0; k < K; k++) {
            p[k] = (word_topic_matrix[documents[m][n]][k] + beta) / (nwsum[k] + V * beta)
                * (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
        }
        
        topic = drawFromProbability(p);

        // add newly estimated z_i to count variables
        word_topic_matrix[documents[m][n]][topic]++;
        nd[m][topic]++;
        nwsum[topic]++;
        ndsum[m]++;

        return topic;
    }
714730
function collectAndCombineLikeTerms(node) {
  if (node.op === '+') {
    const status = collectAndCombineOperation(node);
    if (status.hasChanged()) {
      return status;
    }
    // we might also be able to just combine if they're all the same term
    // e.g. 2x + 4x + x (doesn't need collecting)
    return addLikeTerms(node, true);
  }
  else if (node.op === '*') {
    // collect and combine involves there being coefficients pulled the front
    // e.g. 2x * x^2 * 5x => (2*5) * (x * x^2 * x) => ... => 10 x^4
    const status = collectAndCombineOperation(node);
    if (status.hasChanged()) {
      // make sure there's no * between the coefficient and the symbol part
      status.newNode.implicit = true;
      return status;
    }
    // we might also be able to just combine polynomial terms
    // e.g. x * x^2 * x => ... => x^4
    return multiplyLikeTerms(node, true);
  }
  else {
    return Node.Status.noChange(node);
  }
}
550243
private double delta(double y, double bMutationParameter) {
    double rand = randomGenenerator.getRandomValue();
    int it, maxIt;
    it = currentIteration;
    maxIt = maxIterations;

    return (y * (1.0 -
        Math.pow(rand,
            Math.pow((1.0 - it / (double) maxIt), bMutationParameter)
        )));
  }
318583
def calcChebyshev(coeffs, validDomain, freqs):
    

    logger = logging.getLogger(__name__)

    domain = (validDomain[1] - validDomain[0])[0]
    bins = -1 + 2* n.array([ (freqs[i]-validDomain[0,i])/domain for i in range(len(freqs))])
    ncoeffs = len(coeffs[0])/2
    rr = n.array([n.polynomial.chebyshev.chebval(bins[i], coeffs[i,:ncoeffs]) for i in range(len(coeffs))])
    ll = n.array([n.polynomial.chebyshev.chebval(bins[i], coeffs[i,ncoeffs:]) for i in range(len(coeffs))])

    return rr,ll
1096874
def bernstein(n, t):
    """"""
    t1 = 1-t
    return [n_choose_k(n, k) * t1**(n-k) * t**k for k in range(n+1)]