Compute the average luminosity an experiment was sensitive to

    Assumes that luminosity is uniformly distributed in space.
    Input is the sets of found and missed injections.
180900
def compute_volume_vs_mass(found, missed, mass_bins, bin_type, dbins=None):
    
    # mean and std estimate for luminosity
    volArray = bin_utils.BinnedArray(mass_bins)
    vol2Array = bin_utils.BinnedArray(mass_bins)

    # found/missed stats
    foundArray = bin_utils.BinnedArray(mass_bins)
    missedArray = bin_utils.BinnedArray(mass_bins)

    # compute the mean luminosity in each mass bin
    effvmass = []
    errvmass = []
    # 2D case first
    if bin_type == "Mass1_Mass2":
        for j, mc1 in enumerate(mass_bins.centres()[0]):
            for k, mc2 in enumerate(mass_bins.centres()[1]):
                newfound = filter_injections_by_mass(
                                              found, mass_bins, j, bin_type, k)
                newmissed = filter_injections_by_mass(
                                             missed, mass_bins, j, bin_type, k)

                foundArray[(mc1, mc2)] = len(newfound)
                missedArray[(mc1, mc2)] = len(newmissed)

                # compute the volume using this injection set
                meaneff, efferr, meanvol, volerr = mean_efficiency_volume(
                                                    newfound, newmissed, dbins)
                effvmass.append(meaneff)
                errvmass.append(efferr)
                volArray[(mc1, mc2)] = meanvol
                vol2Array[(mc1, mc2)] = volerr

        return volArray, vol2Array, foundArray, missedArray, effvmass, errvmass

    for j, mc in enumerate(mass_bins.centres()[0]):

        # filter out injections not in this mass bin
        newfound = filter_injections_by_mass(found, mass_bins, j, bin_type)
        newmissed = filter_injections_by_mass(missed, mass_bins, j, bin_type)

        foundArray[(mc, )] = len(newfound)
        missedArray[(mc, )] = len(newmissed)

        # compute the volume using this injection set
        meaneff, efferr, meanvol, volerr = mean_efficiency_volume(
                                                    newfound, newmissed, dbins)
        effvmass.append(meaneff)
        errvmass.append(efferr)
        volArray[(mc, )] = meanvol
        vol2Array[(mc, )] = volerr

    return volArray, vol2Array, foundArray, missedArray, effvmass, errvmass
93092
def mavg(data, xseq, **params):
    
    window = params['method_args']['window']

    # The first average comes after the full window size
    # has been swept over
    rolling = data['y'].rolling(**params['method_args'])
    y = rolling.mean()[window:]
    n = len(data)
    stderr = rolling.std()[window:]
    x = data['x'][window:]
    data = pd.DataFrame({'x': x, 'y': y})
    data.reset_index(inplace=True, drop=True)

    if params['se']:
        df = n - window  # Original - Used
        data['ymin'], data['ymax'] = tdist_ci(
            y, df, stderr, params['level'])
        data['se'] = stderr

    return data
1335550
def example_ylm(m=0, n=2, shape=128, limits=[-4, 4], draw=True, show=True, **kwargs):
    """"""
    import ipyvolume.pylab as p3

    __, __, __, r, theta, phi = xyz(shape=shape, limits=limits, spherical=True)
    radial = np.exp(-(r - 2) ** 2)
    data = np.abs(scipy.special.sph_harm(m, n, theta, phi) ** 2) * radial  # pylint: disable=no-member
    if draw:
        vol = p3.volshow(data=data, **kwargs)
        if show:
            p3.show()
        return vol
    else:
        return data
431402
def dist(self):
        
        if self.func is None:
            return random.uniform(self.min, self.max)
        if self._dist is None:
            x = numpy.arange(self.min, self.max, (self.max-self.min)/1000.0)
            p = self.func(x).cumsum()
            p -= p.min()
            p /= p.max()
            self._dist = interpolate.interp1d(p, x)
        return self._dist(random.random())
1214939
def sample(self, params, argmax_sampling=False):
        
        means = params[:, :, 0]
        log_std = params[:, :, 1]

        if argmax_sampling:
            return means
        else:
            return torch.randn_like(means) * torch.exp(log_std) + means
406880
def calcSMAfromT(self, epsilon=0.7):
        

        return eq.MeanPlanetTemp(self.albedo(), self.star.T, self.star.R, epsilon, self.T).a
242325
function (Y) {
			var N = this.N;
			var dim = this.dim; // dim of output space
			var P = this.P;

			var pmul = this.iter < 100 ? 4 : 1; // trick that helps with local optima

			// compute current Q distribution, unnormalized first
			var Qu = zeros(N * N);
			var qsum = 0.0;
			for (var i = 0; i < N; i++) {
				for (var j = i + 1; j < N; j++) {
					var dsum = 0.0;
					for (var d = 0; d < dim; d++) {
						var dhere = Y[i][d] - Y[j][d];
						dsum += dhere * dhere;
					}
					var qu = 1.0 / (1.0 + dsum); // Student t-distribution
					Qu[i * N + j] = qu;
					Qu[j * N + i] = qu;
					qsum += 2 * qu;
				}
			}
			// normalize Q distribution to sum to 1
			var NN = N * N;
			var Q = zeros(NN);
			for (var q = 0; q < NN; q++) {
				Q[q] = Math.max(Qu[q] / qsum, 1e-100);
			}

			var cost = 0.0;
			var grad = [];
			for (var i = 0; i < N; i++) {
				var gsum = new Array(dim); // init grad for point i
				for (var d = 0; d < dim; d++) {
					gsum[d] = 0.0;
				}
				for (var j = 0; j < N; j++) {
					cost += -P[i * N + j] * Math.log(Q[i * N + j]); // accumulate cost (the non-constant portion at least...)
					var premult = 4 * (pmul * P[i * N + j] - Q[i * N + j]) * Qu[i * N + j];
					for (var d = 0; d < dim; d++) {
						gsum[d] += premult * (Y[i][d] - Y[j][d]);
					}
				}
				grad.push(gsum);
			}

			return {
				cost: cost,
				grad: grad
			};
		}
795457
def rdirichlet(theta, size=1):
    
    gammas = np.vstack([rgamma(theta, 1) for i in xrange(size)])
    if size > 1 and np.size(theta) > 1:
        return (gammas.T / gammas.sum(1))[:-1].T
    elif np.size(theta) > 1:
        return (gammas[0] / gammas[0].sum())[:-1]
    else:
        return 1.
371546
function(data, isHorizontal, precision) {
    var visibles = data.datasets.map(function(dataset) {
      if (!dataset._meta) return true;

      for (var i in dataset._meta) {
        return !dataset._meta[i].hidden;
      }
    });
    
    var datasetDataLength = 0;
    if (data && data.datasets && data.datasets[0] && data.datasets[0].data) {
      datasetDataLength = data.datasets[0].data.length;
    }
    var totals = Array.apply(null, new Array(datasetDataLength)).map(function(el, i) {
      return data.datasets.reduce(function(sum, dataset, j) {
        var key = dataset.stack;
        if (!sum[key]) sum[key] = 0;
        sum[key] += Math.abs(dataValue(dataset.data[i], isHorizontal)) * visibles[j];

        return sum;
      }, {});
    });

    data.calculatedData = data.datasets.map(function(dataset, i) {
      return dataset.data.map(function(val, i) {
        var total = totals[i][dataset.stack];
        var dv = dataValue(val, isHorizontal);
        return dv && total ? round(dv / total, precision) : 0;
      });
    });
  }
839475
def evaluate(self, epoch, iteration, summary):
        
        batch_time = AverageMeter(False)
        tot_tok_per_sec = AverageMeter(False)
        iterations = AverageMeter(False)
        enc_seq_len = AverageMeter(False)
        dec_seq_len = AverageMeter(False)
        stats = {}

        output = []

        for i, (src, indices) in enumerate(self.loader):
            translate_timer = time.time()
            src, src_length = src

            batch_size = self.loader.batch_size
            global_batch_size = batch_size * get_world_size()
            beam_size = self.beam_size

            bos = [self.insert_target_start] * (batch_size * beam_size)
            bos = torch.LongTensor(bos)
            if self.batch_first:
                bos = bos.view(-1, 1)
            else:
                bos = bos.view(1, -1)

            src_length = torch.LongTensor(src_length)
            stats['total_enc_len'] = int(src_length.sum())

            if self.cuda:
                src = src.cuda()
                src_length = src_length.cuda()
                bos = bos.cuda()

            with torch.no_grad():
                context = self.model.encode(src, src_length)
                context = [context, src_length, None]

                if beam_size == 1:
                    generator = self.generator.greedy_search
                else:
                    generator = self.generator.beam_search
                preds, lengths, counter = generator(batch_size, bos, context)

            stats['total_dec_len'] = lengths.sum().item()
            stats['iters'] = counter

            indices = torch.tensor(indices).to(preds)
            preds = preds.scatter(0, indices.unsqueeze(1).expand_as(preds), preds)

            preds = gather_predictions(preds).cpu()

            for pred in preds:
                pred = pred.tolist()
                detok = self.tokenizer.detokenize(pred)
                output.append(detok + '\n')

            elapsed = time.time() - translate_timer
            batch_time.update(elapsed, batch_size)

            total_tokens = stats['total_dec_len'] + stats['total_enc_len']
            ttps = total_tokens / elapsed
            tot_tok_per_sec.update(ttps, batch_size)

            iterations.update(stats['iters'])
            enc_seq_len.update(stats['total_enc_len'] / batch_size, batch_size)
            dec_seq_len.update(stats['total_dec_len'] / batch_size, batch_size)

            if i % self.print_freq == 0:
                log = []
                log += f'TEST '
                if epoch is not None:
                    log += f'[{epoch}]'
                if iteration is not None:
                    log += f'[{iteration}]'
                log += f'[{i}/{len(self.loader)}]\t'
                log += f'Time {batch_time.val:.3f} ({batch_time.avg:.3f})\t'
                log += f'Decoder iters {iterations.val:.1f} ({iterations.avg:.1f})\t'
                log += f'Tok/s {tot_tok_per_sec.val:.0f} ({tot_tok_per_sec.avg:.0f})'
                log = ''.join(log)
                logging.info(log)

        tot_tok_per_sec.reduce('sum')
        enc_seq_len.reduce('mean')
        dec_seq_len.reduce('mean')
        batch_time.reduce('mean')
        iterations.reduce('sum')

        if summary and get_rank() == 0:
            time_per_sentence = (batch_time.avg / global_batch_size)
            log = []
            log += f'TEST SUMMARY:\n'
            log += f'Lines translated: {len(self.loader.dataset)}\t'
            log += f'Avg total tokens/s: {tot_tok_per_sec.avg:.0f}\n'
            log += f'Avg time per batch: {batch_time.avg:.3f} s\t'
            log += f'Avg time per sentence: {1000*time_per_sentence:.3f} ms\n'
            log += f'Avg encoder seq len: {enc_seq_len.avg:.2f}\t'
            log += f'Avg decoder seq len: {dec_seq_len.avg:.2f}\t'
            log += f'Total decoder iterations: {int(iterations.sum)}'
            log = ''.join(log)
            logging.info(log)

        return output
168418
def validateSatellite(config, isochrone, kernel, stellar_mass, distance_modulus, trials=1, debug=False, seed=0):
    
    logger.info('=== Validate Satellite ===')

    config.params['kernel']['params'] = [kernel.r_h] # TODO: Need better solution to update size??
    logger.debug('Using Plummer profile spatial model with half-light radius %.2f deg'%(config.params['kernel']['params'][0]))
    roi = ugali.observation.roi.ROI(config, kernel.lon, kernel.lat)
    simulator = ugali.simulation.simulator.Simulator(config, roi=roi)
    catalog_base = ugali.observation.catalog.Catalog(config,roi=roi)
    mask = ugali.observation.mask.Mask(config, roi)

    coords = (kernel.lon, kernel.lat)
    
    results = {'mc_lon': [],
               'mc_lat': [],
               'mc_distance_modulus': [],
               'mc_stellar_mass': [],
               'mc_radius': [],
               'richness': [],
               'log_likelihood': [],
               'richness_lower': [],
               'richness_upper': [],
               'richness_limit': [],
               'f': [],
               'stellar_mass': []}

    numpy.random.seed(seed)

    for ii in range(0, trials):
        logger.info('=== Running Satellite %i ==='%ii)

        # Simulate
        catalog_satellite = simulator.satellite(isochrone, kernel, stellar_mass, distance_modulus, mc_source_id=1)
        #catalog_bootstrap = catalog_base.bootstrap()
        #catalog_merge = ugali.observation.catalog.mergeCatalogs([catalog_bootstrap, catalog_satellite])
        catalog_background = simulator.satellite(mc_source_id=1)
        catalog_merge = ugali.observation.catalog.mergeCatalogs([catalog_background, catalog_satellite])

        # Analyze
        likelihood = ugali.analysis.likelihood.Likelihood(config, roi, mask, catalog_merge, isochrone, kernel)
                                                               
        likelihood.precomputeGridSearch([distance_modulus])
        richness, log_likelihood, richness_lower, richness_upper, richness_upper_limit, richness_raw, log_likelihood_raw, p, f = likelihood.gridSearch(coords=coords, distance_modulus_index=0)

        results['mc_lon'].append(kernel.lon)
        results['mc_lat'].append(kernel.lat)
        results['mc_distance_modulus'].append(distance_modulus)
        results['mc_stellar_mass'].append(stellar_mass)
        results['mc_radius'].append(kernel.r_h)
        results['richness'].append(richness)
        results['log_likelihood'].append(log_likelihood)
        results['richness_lower'].append(richness_lower)
        results['richness_upper'].append(richness_upper)
        results['richness_limit'].append(richness_upper_limit)
        results['f'].append(f)
        results['stellar_mass'].append(richness * isochrone.stellarMass())

        logger.info('MC Stellar Mass = %.2f, Measured Stellar Mass = %.2f'%(stellar_mass,richness * isochrone.stellarMass()))
        if debug:
            return likelihood, richness, log_likelihood, richness_lower, richness_upper, richness_upper_limit, richness_raw, log_likelihood_raw, p, f

    return results
317026
def kde(data, npoints=_npoints):
    
    # Clipping of severe outliers to concentrate more KDE samples in the parameter range of interest
    mad = np.median(np.fabs(np.median(data) - data))
    cut = (data > np.median(data) - 5. * mad) & (data < np.median(data) + 5. * mad)
    x = data[cut]
    kde = scipy.stats.gaussian_kde(x)
    # No penalty for using a finer sampling for KDE evaluation except computation time
    values = np.linspace(np.min(x), np.max(x), npoints)
    kde_values = kde.evaluate(values)
    peak = values[np.argmax(kde_values)]
    return values[np.argmax(kde_values)], kde.evaluate(peak)
1184225
def sample(self, size=(), rule="R", antithetic=None, verbose=False, **kws):
        
        size_ = numpy.prod(size, dtype=int)
        dim = len(self)
        if dim > 1:
            if isinstance(size, (tuple,list,numpy.ndarray)):
                shape = (dim,) + tuple(size)
            else:
                shape = (dim, size)
        else:
            shape = size

        out = self.kernel.resample(size_)[0]
        try:
            out = out.reshape(shape)
        except:
            if len(self) == 1:
                out = out.flatten()
            else:
                out = out.reshape(dim, out.size/dim)

        return out
274594
def get_irradiance_value(self, month, day, hour):
        """"""
        dt = DateTime(month, day, hour, leap_year=self.is_leap_year)
        count = int(dt.hoy * self.timestep)
        return self.direct_normal_irradiance[count], \
            self.diffuse_horizontal_irradiance[count]
699002
function sample_covariance(x, y) {

        // The two datasets must have the same length which must be more than 1
        if (x.length <= 1 || x.length != y.length){
            return null;
        }

        // determine the mean of each dataset so that we can judge each
        // value of the dataset fairly as the difference from the mean. this
        // way, if one dataset is [1, 2, 3] and [2, 3, 4], their covariance
        // does not suffer because of the difference in absolute values
        var xmean = mean(x),
            ymean = mean(y),
            sum = 0;

        // for each pair of values, the covariance increases when their
        // difference from the mean is associated - if both are well above
        // or if both are well below
        // the mean, the covariance increases significantly.
        for (var i = 0; i < x.length; i++){
            sum += (x[i] - xmean) * (y[i] - ymean);
        }

        // the covariance is weighted by the length of the datasets.
        return sum / (x.length - 1);
    }
111765
def V(self, brightest=False):
        
        mags = self.get_photometry(brightest=brightest, convert=False)
        VT, dVT = mags['VT']
        BT, dBT = mags['BT']
        if (-0.25 < BT - VT < 2.0):
            (a, b, c, d) = (0.00097, 0.1334, 0.05486, 0.01998)
            V = (VT + a - b * (BT - VT) + c * (BT - VT)**2 -
                d * (BT - VT)**3)

            dVdVT = 1 + b - 2*c*(BT-VT) + 3*d*(BT-VT)**2
            dVdBT = -b + 2*c*(BT-VT) - 3*d*(BT-VT)**2
            dV = np.sqrt((dVdVT**2 * dVT**2) + (dVdBT**2*dBT**2))

        else:
            raise ValueError('BT-VT outside of range to convert')

        return V, dV
1357769
def generate(self,rprs=None, mass=None, radius=None,
                n=2e4, fp_specific=0.01, u1=None, u2=None,
                 starmodel=None,
                Teff=None, logg=None, rbin_width=0.3,
                MAfn=None, lhoodcachefile=None):
        

        n = int(n)

        if starmodel is None:
            if type(mass) is type((1,)):
                mass = dists.Gaussian_Distribution(*mass)
            if isinstance(mass, dists.Distribution):
                mdist = mass
                mass = mdist.rvs(1e5)

            if type(radius) is type((1,)):
                radius = dists.Gaussian_Distribution(*radius)
            if isinstance(radius, dists.Distribution):
                rdist = radius
                radius = rdist.rvs(1e5)
        else:
            samples = starmodel.random_samples(1e5)
            mass = samples['mass_0_0'].values
            radius = samples['radius_0_0'].values
            Teff = samples['Teff_0_0'].mean()
            logg = samples['logg_0_0'].mean()

        logging.debug('star mass: {}'.format(mass))
        logging.debug('star radius: {}'.format(radius))
        logging.debug('Teff: {}'.format(Teff))
        logging.debug('logg: {}'.format(logg))

        if u1 is None or u2 is None:
            if Teff is None or logg is None:
                logging.warning('Teff, logg not provided; using solar limb darkening')
                u1 = 0.394; u2=0.296
            else:
                u1,u2 = ldcoeffs(Teff, logg)

        #use point estimate of rprs to construct planets in radius bin
        #rp = self.rprs*np.median(radius)
        #rbin_min = (1-rbin_width)*rp
        #rbin_max = (1+rbin_width)*rp

        rprs_bin_min = (1-rbin_width)*self.rprs
        rprs_bin_max = (1+rbin_width)*self.rprs

        radius_p = radius * (np.random.random(int(1e5))*(rprs_bin_max - rprs_bin_min) + rprs_bin_min)
        mass_p = (radius_p*RSUN/REARTH)**2.06 * MEARTH/MSUN #hokey, but doesn't matter

        logging.debug('planet radius: {}'.format(radius_p))

        stars = pd.DataFrame()
        #df_orbpop = pd.DataFrame() #for orbit population

        tot_prob = None; tot_dprob = None; prob_norm = None
        n_adapt = n
        while len(stars) < n:
            n_adapt = int(n_adapt)
            inds = np.random.randint(len(mass), size=n_adapt)

            #calculate eclipses.
            ecl_inds, df, (prob,dprob) = calculate_eclipses(mass[inds], mass_p[inds],
                                                        radius[inds], radius_p[inds],
                                                        15, np.inf, #arbitrary
                                                        u11s=u1, u21s=u2,
                                                        band=self.band,
                                                        period=self.period,
                                                        calc_mininc=True,
                                                        return_indices=True,
                                                        MAfn=MAfn)

            df['mass_A'] = mass[inds][ecl_inds]
            df['mass_B'] = mass_p[inds][ecl_inds]
            df['radius_A'] = radius[inds][ecl_inds]
            df['radius_B'] = radius_p[inds][ecl_inds]
            df['u1'] = u1 * np.ones_like(df['mass_A'])
            df['u2'] = u2 * np.ones_like(df['mass_A'])
            df['P'] = self.period * np.ones_like(df['mass_A'])

            ok = (df['dpri']>0) & (df['T14_pri'] > 0)

            stars = pd.concat((stars, df[ok]))

            logging.info('{} Transiting planet systems generated (target {})'.format(len(stars),n))
            logging.debug('{} nans in stars[dpri]'.format(np.isnan(stars['dpri']).sum()))

            if tot_prob is None:
                prob_norm = (1/dprob**2)
                tot_prob = prob
                tot_dprob = dprob
            else:
                prob_norm = (1/tot_dprob**2 + 1/dprob**2)
                tot_prob = (tot_prob/tot_dprob**2 + prob/dprob**2)/prob_norm
                tot_dprob = 1/np.sqrt(prob_norm)

            n_adapt = min(int(1.2*(n-len(stars)) * n_adapt//len(df)), 5e4)
            n_adapt = max(n_adapt, 100)

        stars = stars.reset_index()
        stars.drop('index', axis=1, inplace=True)
        stars = stars.iloc[:n]

        stars['mass_1'] = stars['mass_A']
        stars['radius_1'] = stars['radius_A']
        stars['mass_2'] = stars['mass_B']
        stars['radius_2'] = stars['radius_B']

        #make OrbitPopulation?

        #finish below.

        if fp_specific is None:
            rp = stars['radius_2'].mean() * RSUN/REARTH
            fp_specific = fp_fressin(rp)

        priorfactors = {'fp_specific':fp_specific}

        self._starmodel = starmodel

        EclipsePopulation.__init__(self, stars=stars,
                                   period=self.period, cadence=self.cadence,
                                   model=self.model,
                                   priorfactors=priorfactors, prob=tot_prob,
                                   lhoodcachefile=lhoodcachefile)