Calculate a PSD using Bartlett's method
484076
def bartlett(timeseries, segmentlength, **kwargs):
    
    kwargs.pop('noverlap', None)
    return welch(timeseries, segmentlength, noverlap=0, **kwargs)
291460
def _dD_dR(self, R):
        """"""
        HNn_R_sina = self._HNn / R / self._sin_alpha

        return HNn_R_sina * (0.3 * (HNn_R_sina + 0.3 * HNn_R_sina**2. + 1) / R / (0.3 * HNn_R_sina + 1)**2
                             - (1/R * (1 + 0.6 * HNn_R_sina) / (0.3 * HNn_R_sina + 1)))
1222534
def vel_disp(self, kwargs_profile, kwargs_aperture, kwargs_light, kwargs_anisotropy, num=1000):
        
        sigma_s2_sum = 0
        for i in range(0, num):
            sigma_s2_draw = self._vel_disp_one(kwargs_profile, kwargs_aperture, kwargs_light, kwargs_anisotropy)
            sigma_s2_sum += sigma_s2_draw
        sigma_s2_average = sigma_s2_sum/num
        return np.sqrt(sigma_s2_average)
757602
function getSqSegDist(p, p1, p2) {

    var x = p1.x,
        y = p1.y,
        dx = p2.x - x,
        dy = p2.y - y;

    if (dx !== 0 || dy !== 0) {

        var t = ((p.x - x) * dx + (p.y - y) * dy) / (dx * dx + dy * dy);

        if (t > 1) {
            x = p2.x;
            y = p2.y;

        } else if (t > 0) {
            x += dx * t;
            y += dy * t;
        }
    }

    dx = p.x - x;
    dy = p.y - y;

    return dx * dx + dy * dy;
}
743841
def recalculate_psd(self):
        

        seg_len = self.sample_rate * self.psd_segment_length
        e = len(self.strain)
        s = e - ((self.psd_samples + 1) * self.psd_segment_length / 2) * self.sample_rate
        psd = pycbc.psd.welch(self.strain[s:e], seg_len=seg_len, seg_stride=seg_len / 2)

        psd.dist = spa_distance(psd, 1.4, 1.4, self.low_frequency_cutoff) * pycbc.DYN_RANGE_FAC

        # If the new psd is similar to the old one, don't replace it
        if self.psd and self.psd_recalculate_difference:
            if abs(self.psd.dist - psd.dist) / self.psd.dist < self.psd_recalculate_difference:
                logging.info("Skipping recalculation of %s PSD, %s-%s",
                             self.detector, self.psd.dist, psd.dist)
                return True

        # If the new psd is *really* different than the old one, return an error
        if self.psd and self.psd_abort_difference:
            if abs(self.psd.dist - psd.dist) / self.psd.dist > self.psd_abort_difference:
                logging.info("%s PSD is CRAZY, aborting!!!!, %s-%s",
                             self.detector, self.psd.dist, psd.dist)
                self.psd = psd
                self.psds = {}
                return False

        # If the new estimate replaces the current one, invalide the ineterpolate PSDs
        self.psd = psd
        self.psds = {}
        logging.info("Recalculating %s PSD, %s", self.detector, psd.dist)
        return True
229475
def _Backward3_v_Ps(P, s):
    
    if s <= sc:
        return _Backward3a_v_Ps(P, s)
    else:
        return _Backward3b_v_Ps(P, s)
1290787
def computeORFMatrix(psr):
    

    # begin loop over all pulsar pairs and calculate ORF
    npsr = len(psr)
    ORF = N.zeros((npsr, npsr))
    phati = N.zeros(3)
    phatj = N.zeros(3)
    ptheta = [N.pi/2 - p['DECJ'].val for p in psr]
    pphi = [p['RAJ'].val for p in psr]
    for ll in range(0, npsr):
        phati[0] = N.cos(pphi[ll]) * N.sin(ptheta[ll])
        phati[1] = N.sin(pphi[ll]) * N.sin(ptheta[ll])
        phati[2] = N.cos(ptheta[ll])

        for kk in range(0, npsr):
            phatj[0] = N.cos(pphi[kk]) * N.sin(ptheta[kk])
            phatj[1] = N.sin(pphi[kk]) * N.sin(ptheta[kk])
            phatj[2] = N.cos(ptheta[kk])
           
            if ll != kk:
                xip = (1.-N.sum(phati*phatj)) / 2.
                ORF[ll, kk] = 3.*( 1./3. + xip * ( N.log(xip) -1./6.) )
            else:
                ORF[ll, kk] = 2.0

    return ORF
952962
def _calculate_delta_pos(adjacency_arr, pos, t, optimal):
    """"""
    # XXX eventually this should be refactored for the sparse case to only
    # do the necessary pairwise distances
    delta = pos[:, np.newaxis, :] - pos

    # Distance between points
    distance2 = (delta*delta).sum(axis=-1)
    # Enforce minimum distance of 0.01
    distance2 = np.where(distance2 < 0.0001, 0.0001, distance2)
    distance = np.sqrt(distance2)
    # Displacement "force"
    displacement = np.zeros((len(delta), 2))
    for ii in range(2):
        displacement[:, ii] = (
            delta[:, :, ii] *
            ((optimal * optimal) / (distance*distance) -
             (adjacency_arr * distance) / optimal)).sum(axis=1)

    length = np.sqrt((displacement**2).sum(axis=1))
    length = np.where(length < 0.01, 0.1, length)
    delta_pos = displacement * t / length[:, np.newaxis]
    return delta_pos
1102730
def _D(self, R):
        """"""
        return (0.3 * self._HNn**2 / self._sin_alpha / R
                + self._HNn + R * self._sin_alpha) / (0.3 * self._HNn + R * self._sin_alpha)
470981
def in_lamp_reach(p):
    '''  '''
    v1 = XYPoint(Lime.x - Red.x, Lime.y - Red.y)
    v2 = XYPoint(Blue.x - Red.x, Blue.y - Red.y)

    q = XYPoint(p.x - Red.x, p.y - Red.y)
    s = cross_product(q, v2) / cross_product(v1, v2)
    t = cross_product(v1, q) / cross_product(v1, v2)

    return (s >= 0.0) and (t >= 0.0) and (s + t <= 1.0)
848895
def check_psd(matrix, tolerance=1e-6):
    

    hermitian = (matrix + matrix.T.conjugate()) / 2
    eigenvalues = np.linalg.eigh(hermitian)[0]
    return (eigenvalues > -tolerance).all()
292880
function calculateLegacyCurvePoints(datapoints, curvedLinesOptions, yPos) {

			var points = datapoints.points;
			var ps = datapoints.pointsize;
			var num = Number(curvedLinesOptions.curvePointFactor) * (points.length / ps);

			var xdata = new Array;
			var ydata = new Array;

			var curX = -1;
			var curY = -1;
			var j = 0;

			if (curvedLinesOptions.fit) {
				//insert a point before and after the "real" data point to force the line
				//to have a max,min at the data point.

				var fpDist;
				if ( typeof curvedLinesOptions.fitPointDist == 'undefined') {
					//estimate it
					var minX = points[0];
					var maxX = points[points.length - ps];
					fpDist = (maxX - minX) / (500 * 100);
					//x range / (estimated pixel length of placeholder * factor)
				} else {
					//use user defined value
					fpDist = Number(curvedLinesOptions.fitPointDist);
				}

				for (var i = 0; i < points.length; i += ps) {

					var frontX;
					var backX;
					curX = i;
					curY = i + yPos;

					//add point X s
					frontX = points[curX] - fpDist;
					backX = points[curX] + fpDist;

					var factor = 2;
					while (frontX == points[curX] || backX == points[curX]) {
						//inside the ulp
						frontX = points[curX] - (fpDist * factor);
						backX = points[curX] + (fpDist * factor);
						factor++;
					}

					//add curve points
					xdata[j] = frontX;
					ydata[j] = points[curY];
					j++;

					xdata[j] = points[curX];
					ydata[j] = points[curY];
					j++;

					xdata[j] = backX;
					ydata[j] = points[curY];
					j++;
				}
			} else {
				//just use the datapoints
				for (var i = 0; i < points.length; i += ps) {
					curX = i;
					curY = i + yPos;

					xdata[j] = points[curX];
					ydata[j] = points[curY];
					j++;
				}
			}

			var n = xdata.length;

			var y2 = new Array();
			var delta = new Array();
			y2[0] = 0;
			y2[n - 1] = 0;
			delta[0] = 0;

			for (var i = 1; i < n - 1; ++i) {
				var d = (xdata[i + 1] - xdata[i - 1]);
				if (d == 0) {
					//point before current point and after current point need some space in between
					return [];
				}

				var s = (xdata[i] - xdata[i - 1]) / d;
				var p = s * y2[i - 1] + 2;
				y2[i] = (s - 1) / p;
				delta[i] = (ydata[i + 1] - ydata[i]) / (xdata[i + 1] - xdata[i]) - (ydata[i] - ydata[i - 1]) / (xdata[i] - xdata[i - 1]);
				delta[i] = (6 * delta[i] / (xdata[i + 1] - xdata[i - 1]) - s * delta[i - 1]) / p;
			}

			for (var j = n - 2; j >= 0; --j) {
				y2[j] = y2[j] * y2[j + 1] + delta[j];
			}

			//   xmax  - xmin  / #points
			var step = (xdata[n - 1] - xdata[0]) / (num - 1);

			var xnew = new Array;
			var ynew = new Array;
			var result = new Array;

			xnew[0] = xdata[0];
			ynew[0] = ydata[0];

			result.push(xnew[0]);
			result.push(ynew[0]);

			for ( j = 1; j < num; ++j) {
				//new x point (sampling point for the created curve)
				xnew[j] = xnew[0] + j * step;

				var max = n - 1;
				var min = 0;

				while (max - min > 1) {
					var k = Math.round((max + min) / 2);
					if (xdata[k] > xnew[j]) {
						max = k;
					} else {
						min = k;
					}
				}

				//found point one to the left and one to the right of generated new point
				var h = (xdata[max] - xdata[min]);

				if (h == 0) {
					//similar to above two points from original x data need some space between them
					return [];
				}

				var a = (xdata[max] - xnew[j]) / h;
				var b = (xnew[j] - xdata[min]) / h;

				ynew[j] = a * ydata[min] + b * ydata[max] + ((a * a * a - a) * y2[min] + (b * b * b - b) * y2[max]) * (h * h) / 6;

				result.push(xnew[j]);
				result.push(ynew[j]);
			}

			return result;
		}
937924
def many_psds(k=2,fs=1.0, b0=1.0, N=1024):
    
    psd=[]
    for j in range(k):
        print j
        x = noise.white(N=2*4096,b0=b0,fs=fs)
        f, tmp = noise.numpy_psd(x,fs)
        if j==0:
            psd = tmp
        else:
            psd = psd + tmp
    return f, psd/k
4858
function getShapeTwo(x,y,options,modules,mod_matrix){

  var EOL = '\r\n';
  var width = options.width;
  var height = options.height;
  var length = modules.length;
  var style = options.style;
  var xsize = width / (length + 2 * options.padding);
  var ysize = height / (length + 2 * options.padding);
  var px = (x * xsize + options.padding * xsize);
  var py = (y * ysize + options.padding * ysize);

  //If in the middle of a straight line, return a 1b3b
  if( (mod_matrix.top && mod_matrix.bottom) || (mod_matrix.left && mod_matrix.right) ){
    var orientation = (mod_matrix.top && mod_matrix.bottom) ? 0 : 90;
    switch(options.style){
      case 'square':
        return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')"  x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
      case 'circle':
        return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')" rx="'+xsize+'" x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
      case 'sharp-edge':
        return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')" x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
    }
  }

  //If a diagonal is in between the two surrounding (module is on a corner), return a 2a1b1a
  if(mod_matrix.topLeft || mod_matrix.topRight || mod_matrix.bottomLeft || mod_matrix.bottomRight){

    var orientation = 0;
    if(mod_matrix.top && mod_matrix.topRight && mod_matrix.right){orientation=90;}
    else if(mod_matrix.right && mod_matrix.bottomRight && mod_matrix.bottom){orientation=180;}
    else if(mod_matrix.left && mod_matrix.bottomLeft && mod_matrix.bottom){orientation=270;}

    switch(options.style){
      case 'square':
        return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')"  x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
      case 'circle':
        return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')" rx="'+xsize+'" x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
      case 'sharp-edge':
        return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')" x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
    }
  }

  //Else in the center of a right angle
  var orientation = 0;
  if(mod_matrix.top && mod_matrix.right){orientation=90;}
  else if(mod_matrix.right && mod_matrix.bottom){orientation=180;}
  else if(mod_matrix.bottom && mod_matrix.left){orientation=270;}

  switch(options.style){
    case 'square':
      return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')"  x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
    case 'circle':
      return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')"  rx="'+xsize+'" x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
    case 'sharp-edge':
      return '<rect transform="rotate('+orientation+' '+(px+xsize/2)+' '+(py+ysize/2)+')" x="' + px + '" y="' + py + '" width="' + xsize + '" height="' + ysize + '" style="fill:' + options.color + ';shape-rendering:crispEdges;"/>' + EOL;
  }
}
1287038
def vecdiff(htilde, hinterp, sample_points, psd=None):
    
    vecdiffs = numpy.zeros(sample_points.size-1, dtype=float)
    for kk,thisf in enumerate(sample_points[:-1]):
        nextf = sample_points[kk+1]
        vecdiffs[kk] = abs(_vecdiff(htilde, hinterp, thisf, nextf, psd=psd))
    return vecdiffs
852266
function powerSpectrum(amplitudes) {
  var N = amplitudes.length;

  return amplitudes.map(function (a) {
    return (a * a) / N;
  });
}
893125
private void computeOctaveScales() {
		octaveImages[0] = tempImage0;
		for (int i = 1; i < numScales+3; i++) {
			octaveImages[i].reshape(tempImage0.width, tempImage0.height);
			applyGaussian(octaveImages[i - 1], octaveImages[i], kernelSigmaToK[i-1]);
		}

		for (int i = 1; i < numScales+3; i++) {
			differenceOfGaussian[i-1].reshape(tempImage0.width, tempImage0.height);
			PixelMath.subtract(octaveImages[i],octaveImages[i - 1],differenceOfGaussian[i-1]);
		}
	}