Inner loop of MM algorithm for ranking data.
346292
def _mm_rankings(n_items, data, params):
    """"""
    weights = exp_transform(params)
    wins = np.zeros(n_items, dtype=float)
    denoms = np.zeros(n_items, dtype=float)
    for ranking in data:
        sum_ = weights.take(ranking).sum()
        for i, winner in enumerate(ranking[:-1]):
            wins[winner] += 1
            val = 1.0 / sum_
            for item in ranking[i:]:
                denoms[item] += val
            sum_ -= weights[winner]
    return wins, denoms
161818
public void learn(RNNSearch<double[], double[]> range, double[][] samples, double radius, int Nz, double sigma) {
        TrainSample[] training = new TrainSample[samples.length];
        for (int i = 0; i < samples.length; i++) {
            training[i] = new TrainSample();
            training[i].query = samples[i];
            training[i].neighbors = new ArrayList<>();
            ArrayList<Neighbor<double[], double[]>> neighbors = new ArrayList<>();
            range.range(training[i].query, radius, neighbors);
            for (Neighbor<double[], double[]> n : neighbors) {
                training[i].neighbors.add(keys.get(n.index));
            }
        }

        model = new ArrayList<>(hash.size());
        for (int i = 0; i < hash.size(); i++) {
            model.add(new PosterioriModel(hash.get(i), training, Nz, sigma));
        }
    }
128060
def vq_nearest_neighbor(x, means,
                        soft_em=False, num_samples=10, temperature=None):
  """"""
  bottleneck_size = common_layers.shape_list(means)[0]
  x_norm_sq = tf.reduce_sum(tf.square(x), axis=-1, keepdims=True)
  means_norm_sq = tf.reduce_sum(tf.square(means), axis=-1, keepdims=True)
  scalar_prod = tf.matmul(x, means, transpose_b=True)
  dist = x_norm_sq + tf.transpose(means_norm_sq) - 2 * scalar_prod
  if soft_em:
    x_means_idx = tf.multinomial(-dist, num_samples=num_samples)
    x_means_hot = tf.one_hot(
        x_means_idx, depth=common_layers.shape_list(means)[0])
    x_means_hot = tf.reduce_mean(x_means_hot, axis=1)
  else:
    if temperature is None:
      x_means_idx = tf.argmax(-dist, axis=-1)
    else:
      x_means_idx = tf.multinomial(- dist / temperature, 1)
      x_means_idx = tf.squeeze(x_means_idx, axis=-1)
    if (common_layers.should_generate_summaries() and
        not common_layers.is_xla_compiled()):
      tf.summary.histogram("means_idx", tf.reshape(x_means_idx, [-1]))
    x_means_hot = tf.one_hot(x_means_idx, bottleneck_size)
  x_means_hot_flat = tf.reshape(x_means_hot, [-1, bottleneck_size])
  x_means = tf.matmul(x_means_hot_flat, means)
  e_loss = tf.reduce_mean(tf.squared_difference(x, tf.stop_gradient(x_means)))
  return x_means_hot, e_loss, dist
980034
@Reference(authors = "D. Pelleg, A. Moore", //
      title = "X-means: Extending K-means with Efficient Estimation on the Number of Clusters", //
      booktitle = "Proc. 17th Int. Conf. on Machine Learning (ICML 2000)", //
      url = "http://www.pelleg.org/shared/hp/download/xmeans.ps", //
      bibkey = "DBLP:conf/icml/PellegM00")
  public static <V extends NumberVector> double logLikelihood(Relation<V> relation, Clustering<? extends MeanModel> clustering, NumberVectorDistanceFunction<? super V> distanceFunction) {
    List<? extends Cluster<? extends MeanModel>> clusters = clustering.getAllClusters();
    // number of clusters
    final int m = clusters.size();

    // number of objects in the clustering
    int n = 0;
    // cluster sizes
    int[] n_i = new int[m];
    // total variance
    double d = 0.;
    // variances
    double[] d_i = new double[m];

    // Iterate over clusters:
    Iterator<? extends Cluster<? extends MeanModel>> it = clusters.iterator();
    for(int i = 0; it.hasNext(); ++i) {
      Cluster<? extends MeanModel> cluster = it.next();
      n += n_i[i] = cluster.size();
      d += d_i[i] = varianceOfCluster(cluster, distanceFunction, relation);
    }

    // No remaining variance, if every point is on its own:
    if(n <= m) {
      return Double.NEGATIVE_INFINITY;
    }

    // Total variance (corrected for bias)
    final double logv = FastMath.log(d / (n - m));

    final int dim = RelationUtil.dimensionality(relation);
    // log likelihood of this clustering
    double logLikelihood = 0.;

    // Aggregate
    for(int i = 0; i < m; i++) {
      logLikelihood += n_i[i] * FastMath.log(n_i[i]) // Post. entropy Rn log Rn
          - n_i[i] * .5 * MathUtil.LOGTWOPI // Rn/2 log2pi
          - n_i[i] * dim * .5 * logv // Rn M/2 log sigma^2
          - (d_i[i] - m) * .5; // (Rn-K)/2
    }
    logLikelihood -= n * FastMath.log(n); // Prior entropy, sum_i Rn log R
    return logLikelihood;
  }
1096142
def compute_ratings(data=None):
    
    if data is None:
        with sqlite3.connect("ratings.db") as db:
            data = db.execute("select model_winner, model_loser from wins").fetchall()
    model_ids = set([d[0] for d in data]).union(set([d[1] for d in data]))

    # Map model_ids to a contiguous range.
    ordered = sorted(model_ids)
    new_id = {}
    for i, m in enumerate(ordered):
        new_id[m] = i

    # A function to rewrite the model_ids in our pairs
    def ilsr_data(d):
        p1, p2 = d
        p1 = new_id[p1]
        p2 = new_id[p2]
        return (p1, p2)

    pairs = list(map(ilsr_data, data))
    ilsr_param = choix.ilsr_pairwise(
        len(ordered),
        pairs,
        alpha=0.0001,
        max_iter=800)

    hessian = choix.opt.PairwiseFcts(pairs, penalty=.1).hessian(ilsr_param)
    std_err = np.sqrt(np.diagonal(np.linalg.inv(hessian)))

    # Elo conversion
    elo_mult = 400 / math.log(10)

    min_rating = min(ilsr_param)
    ratings = {}

    for model_id, param, err in zip(ordered, ilsr_param, std_err):
        ratings[model_id] = (elo_mult * (param - min_rating), elo_mult * err)

    return ratings
154684
def _rank_optimize(self):
        
        assert self.dag
        for l in reversed(self.layers):
            for v in l:
                gv = self.grx[v]
                for x in v.N(-1):
                   if all((self.grx[y].rank>=gv.rank for y in x.N(+1))):
                        gx = self.grx[x]
                        self.layers[gx.rank].remove(x)
                        gx.rank = gv.rank-1
                        self.layers[gv.rank-1].append(x)
95213
public Object[] mocom(int n, int s, double minn[], double maxn[], int m, int MaxIter) {
        // Start with generating the initial population
        double D[][] = lhsu(minn, maxn, s);
        // Compute the objective function value for each point
        double F[][] = compf(D, m);
        // Now save some important variables
        int nobj = F[0].length;
        // Now do Pareto ranking
        Object ret[] = parrank(F, nobj);
        int R[] = (int[]) ret[0];
        int Rmax = ((Integer) ret[1]).intValue();
        // Now start optimization loop
        int loopcounter = 1;
        while (Rmax > 1) {
            // Assign selection probability P(i) to each of the members
            double P[] = asswght(R, Rmax, s);

            // Construct A to be the points having largest ranks
            ret = worst(D, n, R, Rmax);
            double A[][] = (double[][]) ret[0]; //A <- worst Points
            int L[] = (int[]) ret[1];
            int nA = ((Integer) ret[2]).intValue();

            // Select n points from D to generate Simplex
            ret = choose(P, D, F, R, n);
            double S1[][] = (double[][]) ret[0];
            double F1[][] = (double[][]) ret[1];
            int R1[] = (int[]) ret[2];

            double FA[][] = new double[nA][F[0].length];

            for (int j = 0; j < nA; j++) {
                //addsim
                if (S1.length != F1.length || F1.length != R1.length) {
                    throw new RuntimeException("Component " + this + ": internal error.");
                }
                //build S,SF,SR
                double S[][] = new double[S1.length + 1][S1[0].length];
                double SF[][] = new double[F1.length + 1][F1[0].length];
                int SR[] = new int[R1.length + 1];

                for (int c1 = 0; c1 < S1.length; c1++) {
                    for (int c2 = 0; c2 < S1[0].length; c2++) {
                        S[c1][c2] = S1[c1][c2];
                    }
                    for (int c2 = 0; c2 < F1[0].length; c2++) {
                        SF[c1][c2] = F1[c1][c2];
                    }
                    SR[c1] = R1[c1];
                }
                for (int c1 = 0; c1 < S1[0].length; c1++) {
                    S[S1.length][c1] = A[j][c1];
                }

                for (int c1 = 0; c1 < F1[0].length; c1++) {
                    SF[F1.length][c1] = F[L[j]][c1];
                }
                SR[R1.length] = R[L[j]];

                //mosim
                Object res[] = mosim(S, SF, SR, n, minn, maxn);

                S = (double[][]) res[0];
                SF = (double[][]) res[1];

                for (int c1 = 0; c1 < n; c1++) {
                    A[j][c1] = S[S.length - 1][c1];
                }
                for (int c1 = 0; c1 < FA[j].length; c1++) {
                    FA[j][c1] = SF[SF.length - 1][c1];
                }
            }
            ret = update(D, L, A, n, nA, F, FA);
            D = (double[][]) ret[0];
            F = (double[][]) ret[1];

            // Compute paretorank for each of the parameter sets according to Goldberg, 1989
            ret = parrank(F, nobj);
            R = (int[]) ret[0];
            Rmax = ((Integer) ret[1]).intValue();
            //System.out.println("Evolution Loop:" + loopcounter + " - Trial - " + icall);

            double currentResult[][] = new double[s][N + M + 1];
            for (int i = 0; i < s; i++) {
                for (int j = 0; j < N; j++) {
                    currentResult[i][j] = D[i][j];
                }
                for (int j = 0; j < M; j++) {
                    currentResult[i][N + j] = F[i][j];
                }
                currentResult[i][N + M] = R[i];
            }
            sort(currentResult, N + M);

            int c = 0;
            if (continousOutput) {
                while (c < s && loopcounter % 10 == 0) {
                    //writer.write("Parameter: ");
                    writer.write("" + currentResult[c][N + M] + "\t");
                    for (int i = 0; i < N; i++) {
                        writer.write("" + currentResult[c][i] + "\t");
                    }
                    //System.out.print("Function - Values:");
                    for (int i = 0; i < M; i++) {
                        writer.write("" + currentResult[c][N + i] + "\t");
                    }
                    //System.out.print("Rank:");
                    //System.out.println("");
                    c++;
                    writer.println("");
                    writer.flush();
                }
            }
            loopcounter++;
            if (MaxIter < loopcounter) {
                System.out.println("********************************************");
                System.out.println("---------->OPTIMIZATION STOP<---------------");
                System.out.println("--->MAXIMUM NUMBER OF LOOPS HAS REACHED<----");
                System.out.println("********************************************");
                c = 0;
                System.out.print("Rank\t");
                for (int i = 0; i < N; i++) {
                    System.out.print(this.parameterNames[i] + "\t");
                }
                for (int i = 0; i < M; i++) {
                    System.out.print(this.effNames[i] + "\t");
                }
                System.out.println("");
                while (c < s) {
                    System.out.print("" + currentResult[c][N + M] + "\t");
                    for (int i = 0; i < N; i++) {
                        System.out.print("" + currentResult[c][i] + "\t");
                    }
                    for (int i = 0; i < M; i++) {
                        System.out.print("" + currentResult[c][N + i] + "\t");
                    }
                    System.out.println("");
                    c++;
                }
                c = 0;
                writer.println("Number of model runs: " + icall);
                writer.write("Rank\t");
                for (int i = 0; i < N; i++) {
                    writer.write(this.parameterNames[i] + "\t");
                }
                for (int i = 0; i < M; i++) {
                    writer.write(this.effNames[i] + "\t");
                }
                writer.write("\n");
                while (c < s) {
                    writer.write("" + currentResult[c][N + M] + "\t");
                    for (int i = 0; i < N; i++) {
                        writer.write("" + currentResult[c][i] + "\t");
                    }
                    for (int i = 0; i < M; i++) {
                        writer.write("" + currentResult[c][N + i] + "\t");
                    }
                    c++;
                    writer.println("");
                    writer.flush();
                }
                writer.flush();
                break;
            }
            if (Rmax <= 1) {
                System.out.println("********************************************");
                System.out.println("---------->OPTIMIZATION STOP<---------------");
                System.out.println("-------------->SUCCESSFUL<------------------");
                System.out.println("********************************************");
                c = 0;
                System.out.print("Rank\t");
                for (int i = 0; i < N; i++) {
                    System.out.print(this.parameterNames[i] + "\t");
                }
                for (int i = 0; i < M; i++) {
                    System.out.print(this.effNames[i] + "\t");
                }
                System.out.println("");

                while (c < s) {
                    System.out.print("" + currentResult[c][N + M] + "\t");
                    for (int i = 0; i < N; i++) {
                        System.out.print("" + currentResult[c][i] + "\t");
                    }
                    for (int i = 0; i < M; i++) {
                        System.out.print("" + currentResult[c][N + i] + "\t");
                    }
                    System.out.println("");
                    c++;
                }
                c = 0;
                writer.println("Number of model runs: " + this.icall);
                writer.write("Rank\t");
                for (int i = 0; i < N; i++) {
                    writer.write(this.parameterNames[i] + "\t");
                }
                for (int i = 0; i < M; i++) {
                    writer.write(this.effNames[i] + "\t");
                }
                writer.println();
                while (c < s) {
                    writer.write("" + currentResult[c][N + M] + "\t");
                    for (int i = 0; i < N; i++) {
                        writer.write("" + currentResult[c][i] + "\t");
                    }
                    for (int i = 0; i < M; i++) {
                        writer.write("" + currentResult[c][N + i] + "\t");
                    }
                    c++;
                    writer.println();
                    writer.flush();
                }
                writer.flush();
            }
        }
        return new Object[]{D, F, R, new Integer(Rmax)};
    }
703482
def log_likelihood_rankings(data, params):
    """"""
    loglik = 0
    params = np.asarray(params)
    for ranking in data:
        for i, winner in enumerate(ranking[:-1]):
            loglik -= logsumexp(params.take(ranking[i:]) - params[winner])
    return loglik
347096
function(x, y, z, w) {
      var perm = this.perm;
      var grad4 = this.grad4;

      var n0, n1, n2, n3, n4; // Noise contributions from the five corners
      // Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in
      var s = (x + y + z + w) * F4; // Factor for 4D skewing
      var i = Math.floor(x + s);
      var j = Math.floor(y + s);
      var k = Math.floor(z + s);
      var l = Math.floor(w + s);
      var t = (i + j + k + l) * G4; // Factor for 4D unskewing
      var X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space
      var Y0 = j - t;
      var Z0 = k - t;
      var W0 = l - t;
      var x0 = x - X0; // The x,y,z,w distances from the cell origin
      var y0 = y - Y0;
      var z0 = z - Z0;
      var w0 = w - W0;
      // For the 4D case, the simplex is a 4D shape I won't even try to describe.
      // To find out which of the 24 possible simplices we're in, we need to
      // determine the magnitude ordering of x0, y0, z0 and w0.
      // Six pair-wise comparisons are performed between each possible pair
      // of the four coordinates, and the results are used to rank the numbers.
      var rankx = 0;
      var ranky = 0;
      var rankz = 0;
      var rankw = 0;
      if (x0 > y0) rankx++;
      else ranky++;
      if (x0 > z0) rankx++;
      else rankz++;
      if (x0 > w0) rankx++;
      else rankw++;
      if (y0 > z0) ranky++;
      else rankz++;
      if (y0 > w0) ranky++;
      else rankw++;
      if (z0 > w0) rankz++;
      else rankw++;
      var i1, j1, k1, l1; // The integer offsets for the second simplex corner
      var i2, j2, k2, l2; // The integer offsets for the third simplex corner
      var i3, j3, k3, l3; // The integer offsets for the fourth simplex corner
      // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
      // Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w
      // impossible. Only the 24 indices which have non-zero entries make any sense.
      // We use a thresholding to set the coordinates in turn from the largest magnitude.
      // Rank 3 denotes the largest coordinate.
      i1 = rankx >= 3 ? 1 : 0;
      j1 = ranky >= 3 ? 1 : 0;
      k1 = rankz >= 3 ? 1 : 0;
      l1 = rankw >= 3 ? 1 : 0;
      // Rank 2 denotes the second largest coordinate.
      i2 = rankx >= 2 ? 1 : 0;
      j2 = ranky >= 2 ? 1 : 0;
      k2 = rankz >= 2 ? 1 : 0;
      l2 = rankw >= 2 ? 1 : 0;
      // Rank 1 denotes the second smallest coordinate.
      i3 = rankx >= 1 ? 1 : 0;
      j3 = ranky >= 1 ? 1 : 0;
      k3 = rankz >= 1 ? 1 : 0;
      l3 = rankw >= 1 ? 1 : 0;
      // The fifth corner has all coordinate offsets = 1, so no need to compute that.
      var x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords
      var y1 = y0 - j1 + G4;
      var z1 = z0 - k1 + G4;
      var w1 = w0 - l1 + G4;
      var x2 = x0 - i2 + 2.0 * G4; // Offsets for third corner in (x,y,z,w) coords
      var y2 = y0 - j2 + 2.0 * G4;
      var z2 = z0 - k2 + 2.0 * G4;
      var w2 = w0 - l2 + 2.0 * G4;
      var x3 = x0 - i3 + 3.0 * G4; // Offsets for fourth corner in (x,y,z,w) coords
      var y3 = y0 - j3 + 3.0 * G4;
      var z3 = z0 - k3 + 3.0 * G4;
      var w3 = w0 - l3 + 3.0 * G4;
      var x4 = x0 - 1.0 + 4.0 * G4; // Offsets for last corner in (x,y,z,w) coords
      var y4 = y0 - 1.0 + 4.0 * G4;
      var z4 = z0 - 1.0 + 4.0 * G4;
      var w4 = w0 - 1.0 + 4.0 * G4;
      // Work out the hashed gradient indices of the five simplex corners
      var ii = i & 255;
      var jj = j & 255;
      var kk = k & 255;
      var ll = l & 255;
      // Calculate the contribution from the five corners
      var t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
      if (t0 < 0) n0 = 0.0;
      else {
        var gi0 = (perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32) * 4;
        t0 *= t0;
        n0 = t0 * t0 * (grad4[gi0] * x0 + grad4[gi0 + 1] * y0 + grad4[gi0 + 2] * z0 + grad4[gi0 + 3] * w0);
      }
      var t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
      if (t1 < 0) n1 = 0.0;
      else {
        var gi1 = (perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32) * 4;
        t1 *= t1;
        n1 = t1 * t1 * (grad4[gi1] * x1 + grad4[gi1 + 1] * y1 + grad4[gi1 + 2] * z1 + grad4[gi1 + 3] * w1);
      }
      var t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
      if (t2 < 0) n2 = 0.0;
      else {
        var gi2 = (perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32) * 4;
        t2 *= t2;
        n2 = t2 * t2 * (grad4[gi2] * x2 + grad4[gi2 + 1] * y2 + grad4[gi2 + 2] * z2 + grad4[gi2 + 3] * w2);
      }
      var t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
      if (t3 < 0) n3 = 0.0;
      else {
        var gi3 = (perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32) * 4;
        t3 *= t3;
        n3 = t3 * t3 * (grad4[gi3] * x3 + grad4[gi3 + 1] * y3 + grad4[gi3 + 2] * z3 + grad4[gi3 + 3] * w3);
      }
      var t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
      if (t4 < 0) n4 = 0.0;
      else {
        var gi4 = (perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32) * 4;
        t4 *= t4;
        n4 = t4 * t4 * (grad4[gi4] * x4 + grad4[gi4 + 1] * y4 + grad4[gi4 + 2] * z4 + grad4[gi4 + 3] * w4);
      }
      // Sum up and scale the result to cover the range [-1,1]
      return 27.0 * (n0 + n1 + n2 + n3 + n4);
    }
880894
def lranksums(x,y):
    
    n1 = len(x)
    n2 = len(y)
    alldata = x+y
    ranked = rankdata(alldata)
    x = ranked[:n1]
    y = ranked[n1:]
    s = sum(x)
    expected = n1*(n1+n2+1) / 2.0
    z = (s - expected) / math.sqrt(n1*n2*(n1+n2+1)/12.0)
    prob = 2*(1.0 -zprob(abs(z)))
    return z, prob
1026901
def rank(tensor: BKTensor) -> int:
    """"""
    if isinstance(tensor, np.ndarray):
        return len(tensor.shape)

    return len(tensor[0].size())
270164
def _build_mat(self, ranks, n_candidates, prefcounts):
        

        mat = zeros((n_candidates, n_candidates))
        for i, j in itertools.combinations(range(n_candidates), 2):
            preference = ranks[:, i] - ranks[:, j]
            h_ij = dot((preference < 0), prefcounts)  # prefers i to j
            h_ji = dot((preference > 0), prefcounts)  # prefers j to i
            mat[i, j] = h_ij - h_ji
            mat[j, i] = h_ji - h_ij
        return mat
111804
def EvalBinomialPmf(k, n, p):
    
    return scipy.stats.binom.pmf(k, n, p)
294565
def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None):
        

        # Select objects in annulus
        cut_annulus = self.roi.inAnnulus(catalog.lon,catalog.lat)
        mag_1 = catalog.mag_1[cut_annulus]
        mag_2 = catalog.mag_2[cut_annulus]

        # Units are (deg^2)
        solid_angle = ugali.utils.binning.take2D(self.solid_angle_mmd, mag_2, mag_1,
                                                 self.roi.bins_mag, self.roi.bins_mag)

        # Weight each object before binning
        # Divide by solid angle and bin size in magnitudes to get number density 
        # [objs / deg^2 / mag^2]
        if weights is None:
            number_density = (solid_angle*self.roi.delta_mag**2)**(-1)
        else:
            number_density = weights*(solid_angle*self.roi.delta_mag**2)**(-1)

        method = str(method).lower()
        if method == 'cloud-in-cells':
            # Apply cloud-in-cells algorithm
            mmd_background = ugali.utils.binning.cloudInCells(mag_2,mag_1,
                                                              [self.roi.bins_mag,self.roi.bins_mag],
                                                              weights=number_density)[0]
        elif method == 'bootstrap':
            # Not implemented
            raise ValueError("Bootstrap method not implemented")
            mag_1 + (mag_1_err * np.random.normal(0, 1., len(mag_1)))
            mag_2 + (mag_2_err * np.random.normal(0, 1., len(mag_2)))

        elif method == 'histogram':
            # Apply raw histogram
            mmd_background = np.histogram2d(mag_1,mag_2,bins=[self.roi.bins_mag,self.roi.bins_mag],
                                            weights=number_density)[0]

        elif method == 'kde':
            # Gridded kernel density estimator
            logger.warning("### KDE not implemented properly")
            mmd_background = ugali.utils.binning.kernelDensity(mag_2,mag_1,
                                                               [self.roi.bins_mag,self.roi.bins_mag],
                                                               weights=number_density)[0]
        elif method == 'uniform':
            logger.warning("### WARNING: Uniform CMD")
            hist = np.histogram2d(mag_1,mag_2,bins=[self.roi.bins_mag,self.roi.bins_mag], 
                                  weights=number_density)[0]
            mmd_background = np.mean(hist)*np.ones(hist.shape)
            observable = (self.solid_angle_mmd > self.minimum_solid_angle)
            mmd_background *= observable
            return mmd_background
        else:
            raise ValueError("Unrecognized method: %s"%method)
        ## Account for the objects that spill out of the observable space
        ## But what about the objects that spill out to red colors??
        #for index_color in range(0, len(self.roi.centers_color)):
        #    for index_mag in range(0, len(self.roi.centers_mag)):
        #        if self.solid_angle_cmd[index_mag][index_color] < self.minimum_solid_angle:
        #            cmd_background[index_mag - 1][index_color] += cmd_background[index_mag][index_color]
        #            cmd_background[index_mag][index_color] = 0.
        #            break

        mmd_area = self.solid_angle_mmd*self.roi.delta_mag**2 # [deg^2 * mag^2]

        # ADW: This accounts for leakage to faint magnitudes
        # But what about the objects that spill out to red colors??
        # Maximum obsevable magnitude index for each color (uses the fact that
        # np.argmin returns first minimum (zero) instance found.
        # NOTE: More complicated maps may have holes causing problems

        observable = (self.solid_angle_mmd > self.minimum_solid_angle)
        index_mag_1 = observable.argmin(axis=0) - 1
        index_mag_2 = np.arange(len(self.roi.centers_mag))
        # Add the cumulative leakage back into the last bin of the CMD
        leakage = (mmd_background * ~observable).sum(axis=0)
        ### mmd_background[[index_mag_1,index_mag_2]] += leakage
        # Zero out all non-observable bins
        ### mmd_background *= observable

        # Avoid dividing by zero by setting empty bins to the value of the 
        # minimum filled bin of the CMD. This choice is arbitrary and 
        # could be replaced by a static minimum, some fraction of the 
        # CMD maximum, some median clipped minimum, etc. However, should 
        # be robust against outliers with very small values.
        min_mmd_background = max(mmd_background[mmd_background > 0.].min(),
                                 1e-4*mmd_background.max())
        mmd_background[observable] = mmd_background[observable].clip(min_mmd_background)

        ### # ADW: This is a fudge factor introduced to renormalize the CMD
        ### # to the number of input stars in the annulus. While leakage
        ### # will still smooth the distribution, it shouldn't result in 
        ### fudge_factor = len(mag) / float((cmd_background*cmd_area).sum())
        ### cmd_background *= fudge_factor

        return mmd_background
1337793
function arrayMin(data) {
	var i = data.length,
		min = data[0];

	while (i--) {
		if (data[i] < min) {
			min = data[i];
		}
	}
	return min;
}
206597
function rank(A, tol) {
    var S = svd(A, false, false)[0];
    if (tol === undefined) {
        tol = common_1.EPSILON; // TODO : use numpy's formula
    }
    var rank = 0;
    for (var _i = 0, _a = S.toArray(); _i < _a.length; _i++) {
        var n = _a[_i];
        if (n > tol) {
            rank++;
        }
    }
    return rank;
}
303137
def rank(self, n, mu, sigma, crit=.5, upper=10000, xtol=1):
        

        return _make_rank(self, n, mu, sigma, crit=crit, upper=upper,
                                                                    xtol=xtol)