>

Source code for minke.distribution

from scipy import random
import numpy

[docs]def sky_params(net, time, ra, dec, psi, incl=numpy.pi): if not (0 < ra < 2 * numpy.pi) or not (0 < psi < 2 * numpy.pi) or not(-numpy.pi/2 < dec < numpy.pi/2): return 0.0 fp, fx = 0, 0 for ifo in net: tp, tx, _, _ = response(time, ra, dec, incl, psi, 'radians', ifo) fp += tp**2 fx += fx**2 fp, fx = math.sqrt(fp), math.sqrt(fx) return (fp + fx) / len(net)
[docs]def uniform_dec(num): """ Declination distribution: uniform in sin(dec), which leads to a uniform distribution across all declinations. Parameters ---------- num : int The number of random declinations to produce. """ return (numpy.pi / 2.) - numpy.arccos(2 * random.random_sample(num) - 1)
[docs]def uniform_theta(num): """ Uniform in cos distribution. num controls the number of draws. """ return numpy.arccos(2 * random.random_sample(num) - 1)
[docs]def uniform_phi(num): """ Uniform in (0, 2pi) distribution. num controls the number of draws. """ return random.random_sample(num) * 2 * numpy.pi
[docs]def supernova_angle(num, divisions = 10): """ Draw from a discrete distribution of angles. """ theta = numpy.linspace(-1, 1, divisions) ctheta = numpy.arccos(theta) phi = numpy.linspace(0, numpy.pi*2, divisions) out_t = numpy.random.choice(ctheta, num, replace = True) out_p = numpy.random.choice(phi, num, replace = True) return zip(out_t, out_p)
[docs]def uniform_sky(number=1): """ Get a set of (RA, declination, polarization) drawn from an isotropic distribution over the whole sky. Parameters ---------- number : int The number of random sky locations and polarisations to be produced. Returns ------- ra : array of float Randomly drawn right ascensions. dec : array of float Randomly drawn declinations. pol : array of float Randomly drawn polarisations. """ expnum = number ra = uniform_phi(expnum) dec = uniform_dec(expnum) pol = uniform_phi(expnum) return ra, dec, pol
[docs]def favorable_sky(net, time): """ Wander through the skies, searching for a most favorable location. --- draw extrinsic parameters as if the network antenna pattern magnitude were the PDF. """ ndraws = len(time) ra_out, dec_out, psi_out = numpy.empty((3, len(time))) while ndraws > 0: ra = numpy.random.uniform(0, 2 * numpy.pi, 1000) dec = numpy.random.uniform(-numpy.pi / 2, numpy.pi / 2, 1000) psi = numpy.random.uniform(0, 2 * numpy.pi, 1000) rnd = numpy.random.uniform(0, 1, 1000) for r, d, p, n in zip(ra, dec, psi, rnd): if n < sky_params(net, time, r, d, p): ndraws -= 1 ra_out[ndraws] = r dec_out[ndraws] = d psi_out[ndraws] = p break return ra_out, dec_out, psi_out
[docs]def uniform_interval(interval, num): """ Return a number, or a list of numbers which are sampled from a uniform distribution. Parameters ---------- interval : tuple (lower, upper) The interval which the uniform distribution covers. num : int The number of samples which should be drawn from the distribution. Returns ------- sample : float or list A sample, or a list of samples. Notes ----- http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.random_sample.html#numpy.random.random_sample """ return (interval[1] - interval[0]) * random.random_sample(num) + interval[0]
[docs]def uniform_time(tstart, tstop, number): """ Get a set of randomized (integer) event times. """ return random.randint(tstart, tstop, number) + random.rand(number)
[docs]def log_uniform(lower, upper, number): """ Draw uniformly in the log of a predefined range. Parameters ---------- lower : float The lower hrss value (n.b. not the lower log(hrss) value) upper : float The upper hrss value (n.b. not the upper log(hrss) value) number : int The number of samples to be drawn from the distribution. Returns ------- sample : float An sample value drawn from the log uniform distribution. """ log10h = (numpy.log10(lower), numpy.log10(upper)) return 10**uniform_interval(log10h, number)
[docs]def even_time(start, stop, rate, jitter=0): """ Produce an evenly-distributed set of times which has some jitter. Parameters ---------- start : int The gps start time of the set. stop : int The gps end time of the set. rate : int The rate of events per year. jitter : float The number of seconds around which the time can "jitter". This has the effect of adding a small random time on to each time returned in the distribution, up to a maximum of half the jitter value, or taking away up to half the jitter value. Returns ------- times : ndarray A numpy array of gps times where injections should be made. """ expnum_exact = (stop-start) * rate / 365.0 / 24 / 3600 interval = (stop - start) / expnum_exact time = numpy.array([ start + i*interval + jitter *(random.rand()-0.5) for i in range(1,int(expnum_exact)+1) ]) return time
[docs]def burst_dist(minimum, maximum, size=1): """ Generate an hrss drawn from the distance distribution \[ r + 50/r \] used for Burst allsky mock data challenges. Parameters ---------- minimum : float The lowest hrss to be produced. maximum : float The largest hrss to be produced. size : int The number of draws to produce. Defaults to 1. """ #extent = numpy.log10(maximum) - numpy.log10(minimum) # Convert hrss to distance minimum, maximum = 1/minimum, 1/maximum #Set-up the distributions limits a = (minimum + 50.0/minimum)**2 b = (maximum + 50.0/maximum)**2 rnd = numpy.random.uniform(a,b, size) # Convert back to an hrss return 1.0/numpy.sqrt(rnd)