For 3D points, computes the median value and variance along each dimension.
433414
void computePointStatistics(Point[] points ) {
		final int length = points.length;
		double v[] = new double[length];

		for (int axis = 0; axis < 3; axis++) {
			double maxAbs = 0;
			for (int i = 0; i < length; i++) {
				v[i] = points[i].coordinate[axis];
				maxAbs = Math.max( maxAbs , Math.abs(v[i]));
			}

			double median = QuickSelect.select(v,length/2,length);
			switch( axis ) {
				case 0: medianPoint.x = median; break;
				case 1: medianPoint.y = median; break;
				case 2: medianPoint.z = median; break;
			}
		}

		for (int i = 0; i < length; i++) {
			v[i] = points[i].distanceSq(medianPoint);
		}
		medianDistancePoint = Math.sqrt(QuickSelect.select(v,length/2,length));

//		System.out.println("Median P ="+ medianPoint);
//		System.out.println("Median R ="+ medianDistancePoint);
//		System.out.println("Scale    ="+ (desiredDistancePoint / medianDistancePoint));
	}
1047550
def _evaluate(self, indices, norm_distances, out=None):
        """"""
        # slice for broadcasting over trailing dimensions in self.values
        vslice = (slice(None),) + (None,) * (self.values.ndim - len(indices))

        if out is None:
            out_shape = out_shape_from_meshgrid(norm_distances)
            out_dtype = self.values.dtype
            out = np.zeros(out_shape, dtype=out_dtype)
        else:
            out[:] = 0.0

        # Weights and indices (per axis)
        low_weights, high_weights, edge_indices = _create_weight_edge_lists(
            indices, norm_distances, self.schemes, self.nn_variants)

        # Iterate over all possible combinations of [i, i+1] for each
        # axis, resulting in a loop of length 2**ndim
        for lo_hi, edge in zip(product(*([['l', 'h']] * len(indices))),
                               product(*edge_indices)):
            weight = 1.0
            # TODO: determine best summation order from array strides
            for lh, w_lo, w_hi in zip(lo_hi, low_weights, high_weights):

                # We don't multiply in-place to exploit the cheap operations
                # in the beginning: sizes grow gradually as following:
                # (n, 1, 1, ...) -> (n, m, 1, ...) -> ...
                # Hence, it is faster to build up the weight array instead
                # of doing full-size operations from the beginning.
                if lh == 'l':
                    weight = weight * w_lo
                else:
                    weight = weight * w_hi
            out += np.asarray(self.values[edge]) * weight[vslice]
        return np.array(out, copy=False, ndmin=1)
705171
private static int aggregateDimValue(
      final int[] positions,
      final BufferAggregator[] theAggregators,
      final ByteBuffer resultsBuf,
      final int numBytesPerRecord,
      final int[] aggregatorOffsets,
      final int aggSize,
      final int aggExtra,
      final int dimIndex,
      int currentPosition
  )
  {
    if (SKIP_POSITION_VALUE == positions[dimIndex]) {
      return currentPosition;
    }
    if (INIT_POSITION_VALUE == positions[dimIndex]) {
      positions[dimIndex] = currentPosition * numBytesPerRecord;
      currentPosition++;
      final int pos = positions[dimIndex];
      for (int j = 0; j < aggSize; ++j) {
        theAggregators[j].init(resultsBuf, pos + aggregatorOffsets[j]);
      }
    }
    final int position = positions[dimIndex];

    switch (aggExtra) {
      case 7:
        theAggregators[6].aggregate(resultsBuf, position + aggregatorOffsets[6]);
        // fall through
      case 6:
        theAggregators[5].aggregate(resultsBuf, position + aggregatorOffsets[5]);
        // fall through
      case 5:
        theAggregators[4].aggregate(resultsBuf, position + aggregatorOffsets[4]);
        // fall through
      case 4:
        theAggregators[3].aggregate(resultsBuf, position + aggregatorOffsets[3]);
        // fall through
      case 3:
        theAggregators[2].aggregate(resultsBuf, position + aggregatorOffsets[2]);
        // fall through
      case 2:
        theAggregators[1].aggregate(resultsBuf, position + aggregatorOffsets[1]);
        // fall through
      case 1:
        theAggregators[0].aggregate(resultsBuf, position + aggregatorOffsets[0]);
    }
    for (int j = aggExtra; j < aggSize; j += AGG_UNROLL_COUNT) {
      theAggregators[j].aggregate(resultsBuf, position + aggregatorOffsets[j]);
      theAggregators[j + 1].aggregate(resultsBuf, position + aggregatorOffsets[j + 1]);
      theAggregators[j + 2].aggregate(resultsBuf, position + aggregatorOffsets[j + 2]);
      theAggregators[j + 3].aggregate(resultsBuf, position + aggregatorOffsets[j + 3]);
      theAggregators[j + 4].aggregate(resultsBuf, position + aggregatorOffsets[j + 4]);
      theAggregators[j + 5].aggregate(resultsBuf, position + aggregatorOffsets[j + 5]);
      theAggregators[j + 6].aggregate(resultsBuf, position + aggregatorOffsets[j + 6]);
      theAggregators[j + 7].aggregate(resultsBuf, position + aggregatorOffsets[j + 7]);
    }
    return currentPosition;
  }
1203803
private static double[] computePerDimensionVariances(Relation<? extends NumberVector> relation, double[] center, DBIDs neighborhood) {
    final int dim = center.length;
    double[] variances = new double[dim];
    for(DBIDIter iter = neighborhood.iter(); iter.valid(); iter.advance()) {
      NumberVector databaseObject = relation.get(iter);
      for(int d = 0; d < dim; d++) {
        final double deviation = databaseObject.doubleValue(d) - center[d];
        variances[d] += deviation * deviation;
      }
    }
    return VMath.timesEquals(variances, 1. / neighborhood.size());
  }
727358
def averages(self):
        
        if not self.mesh._compute_at_vertices:
            return None

        return np.mean(self.vertices_per_triangle, axis=1)
806427
def pld3(pos, line_vertex, line_dir):
    """"""
    pos = np.atleast_2d(pos)
    line_vertex = np.atleast_1d(line_vertex)
    line_dir = np.atleast_1d(line_dir)
    c = np.cross(line_dir, line_vertex - pos)
    n1 = np.linalg.norm(c, axis=1)
    n2 = np.linalg.norm(line_dir)
    out = n1 / n2
    if out.ndim == 1 and len(out) == 1:
        return out[0]
    return out
799510
def _compute(self):
        """"""
        delta = 2 * pi / self._len if self._len else 0
        self._x_pos = [.5 * pi + i * delta for i in range(self._len + 1)]
        for serie in self.all_series:
            serie.points = [(v, self._x_pos[i])
                            for i, v in enumerate(serie.values)]
            if self.interpolate:
                extended_x_pos = ([.5 * pi - delta] + self._x_pos)
                extended_vals = (serie.values[-1:] + serie.values)
                serie.interpolated = list(
                    map(
                        tuple,
                        map(
                            reversed,
                            self._interpolate(extended_x_pos, extended_vals)
                        )
                    )
                )

        # x labels space
        self._box.margin *= 2
        self._rmin = self.zero
        self._rmax = self._max or 1
        self._box.set_polar_box(self._rmin, self._rmax)
        self._self_close = True
636783
void compute_Rv() {
		double t = Math.sqrt(v1*v1 + v2*v2);
		double s = Math.sqrt(t*t + 1);
		double cosT = 1.0/s;
		double sinT = Math.sqrt(1-1.0/(s*s));

		K_x.data[0] = 0;   K_x.data[1] = 0;   K_x.data[2] = v1;
		K_x.data[3] = 0;   K_x.data[4] = 0;   K_x.data[5] = v2;
		K_x.data[6] = -v1; K_x.data[7] = -v2; K_x.data[8] = 0;
		CommonOps_DDRM.divide(K_x,t);
		CommonOps_DDRM.setIdentity(R_v);
		CommonOps_DDRM.addEquals(R_v,sinT,K_x);
		CommonOps_DDRM.multAdd(1.0-cosT,K_x,K_x,R_v);
	}
465040
def BatchNorm(x, params, axis=(0, 1, 2), epsilon=1e-5,
              center=True, scale=True, **unused_kwargs):
  """"""
  mean = np.mean(x, axis, keepdims=True)
  # Fast but less numerically-stable variance calculation than np.var.
  m1 = np.mean(x**2, axis, keepdims=True)
  var = m1 - mean**2
  z = (x - mean) / np.sqrt(var + epsilon)

  # Expand the parameters to have the right axes.
  beta, gamma = params
  # TODO(phawkins): np.expand_dims should accept an axis tuple.
  # (https://github.com/numpy/numpy/issues/12290)
  ed = tuple(None if i in axis else slice(None) for i in range(np.ndim(x)))
  beta = beta[ed]
  gamma = gamma[ed]

  # Return the z rescaled by the parameters if requested.
  if center and scale:
    return gamma * z + beta
  if center:
    return z + beta
  if scale:
    return gamma * z
  return z
1242564
def numpy_weighted_median(data, weights=None):
    """"""
    import numpy as np
    if weights is None:
        return np.median(np.array(data).flatten())
    data, weights = np.array(data).flatten(), np.array(weights).flatten()
    if any(weights > 0):
        sorted_data, sorted_weights = map(np.array, zip(*sorted(zip(data, weights))))
        midpoint = 0.5 * sum(sorted_weights)
        if any(weights > midpoint):
            return (data[weights == np.max(weights)])[0]
        cumulative_weight = np.cumsum(sorted_weights)
        below_midpoint_index = np.where(cumulative_weight <= midpoint)[0][-1]
        if cumulative_weight[below_midpoint_index] - midpoint < sys.float_info.epsilon:
            return np.mean(sorted_data[below_midpoint_index:below_midpoint_index+2])
        return sorted_data[below_midpoint_index+1]
185454
def centroid(data, method=median):
    ""
    out = []
    for i in range(len(data[0])):
        out.append(method([x[i] for x in data]))
    return tuple(out)
207781
def median(values):
    
    values.sort()
    n = int(len(values) / 2)
    return values[n]
571831
def median(self):
        """"""
        mu = self.mean()
        ret_val = math.exp(mu)
        if math.isnan(ret_val):
            ret_val = float("inf")
        return ret_val
620931
def computeNormals(self):
        
        poly = self.polydata(False)
        pnormals = poly.GetPointData().GetNormals()
        cnormals = poly.GetCellData().GetNormals()
        if pnormals and cnormals:
            return self

        pdnorm = vtk.vtkPolyDataNormals()
        pdnorm.SetInputData(poly)
        pdnorm.ComputePointNormalsOn()
        pdnorm.ComputeCellNormalsOn()
        pdnorm.FlipNormalsOff()
        pdnorm.ConsistencyOn()
        pdnorm.Update()
        return self.updateMesh(pdnorm.GetOutput())
666434
def var(self):
        
        mn = self.mean
        vr = np.mean((self._mcpts - mn) ** 2)
        return vr
498994
def _amp2d_to_3d(self, amp, sigma_x, sigma_y):
        
        return amp / (np.sqrt(np.pi) * np.sqrt(sigma_x * sigma_y * 2))
703964
protected void computeStatistics(T input, int innerWidth, int innerHeight) {
		original.init(blockWidth,blockHeight,thresholdFromLocalBlocks);

		int statPixelStride = stats.getImageType().getNumBands();

		int indexStats = 0;
		for (int y = 0; y < innerHeight; y += blockHeight) {
			for (int x = 0; x < innerWidth; x += blockWidth, indexStats += statPixelStride) {
				original.computeBlockStatistics(x,y,blockWidth,blockHeight,indexStats,input,stats);
			}
			// handle the case where the image's width isn't evenly divisible by the block's width
			if( innerWidth != input.width ) {
				original.computeBlockStatistics(innerWidth,y,input.width-innerWidth,blockHeight,indexStats,input,stats);
				indexStats += statPixelStride;
			}
		}
		// handle the case where the image's height isn't evenly divisible by the block's height
		if( innerHeight != input.height ) {
			int y = innerHeight;
			int blockHeight = input.height-innerHeight;
			for (int x = 0; x < innerWidth; x += blockWidth, indexStats += statPixelStride) {
				original.computeBlockStatistics(x,y,blockWidth,blockHeight,indexStats,input,stats);
			}
			if( innerWidth != input.width ) {
				original.computeBlockStatistics(innerWidth,y,input.width-innerWidth,blockHeight,indexStats,input,stats);
			}
		}
	}