set sample and trait designs
        F:      NxK sample design
        A:      LxP sample design
        fast_computations:   False deactivates the fast computations for any and common effects (for debugging)
168815
def addFixedEffect(self,F=None,A=None,index=None):
        
        if F is None:
            F = np.ones((self.N,1))
        else:
            assert F.shape[0]==self.N, "F dimension mismatch"

        if ((A is None) or ( (A.shape == (self.P,self.P)) and (A==np.eye(self.P)).all() )):
            #case any effect
            self.F_any = np.hstack((self.F_any,F))
        elif (index is not None) and  ((A==self.A[index]).all()):
            #case common effect
            self.F[index] = np.hstack((self.F_index,F))
        else:
            #case general A
            assert A.shape[1]==self.P, "A dimension mismatch"
            self.F.append(F)
            self.A.append(A)

        self.clear_cache()
1246413
def quantitate_expression_noparallel(samples, run_parallel):
    
    data = samples[0][0]
    if "express" in dd.get_expression_caller(data):
        samples = run_parallel("run_express", samples)
    if "dexseq" in dd.get_expression_caller(data):
        samples = run_parallel("run_dexseq", samples)
    return samples
325784
def sampleset(self):
        """"""
        x1, x4 = self.bounding_box
        dw = self.width * 0.5
        x2 = self.x_0 - dw
        x3 = self.x_0 + dw

        if self._n_models == 1:
            w = [x1, x2, x3, x4]
        else:
            w = list(zip(x1, x2, x3, x4))

        return np.asarray(w)
706452
def _cache_key_select_sample_type(method, self, allow_blank=True, multiselect=False, style=None):
    
    key = update_timer(), allow_blank, multiselect, style
    return key
580441
def run_checks(number_samples, k_choices):
        
        assert isinstance(k_choices, int), \
            "Number of optimal trajectories should be an integer"

        if k_choices < 2:
            raise ValueError(
                "The number of optimal trajectories must be set to 2 or more.")
        if k_choices >= number_samples:
            msg = "The number of optimal trajectories should be less than the \
                    number of samples"
            raise ValueError(msg)
586554
def _xval_choose_funcs(self, fit_spectra, reject_outliers, fit_lb, fit_ub,
                           fitters=[ana.fit_gaussian,ana.fit_two_gaussian],
                           funcs = [ut.gaussian, ut.two_gaussian]):
        

        set1 = fit_spectra[::2]
        set2 = fit_spectra[1::2]

        errs = []
        signal_select = [] 
        # We can loop over functions and try each one out, checking the
        # error in each:
        for fitter in fitters:
            models = []
            signals = []
            for this_set in [set1, set2]:
                choose_transients, model, signal, params, this_idx =\
                    self._fit_helper(this_set, reject_outliers,
                                    fit_lb, fit_ub, fitter)
                models.append(np.nanmean(model[choose_transients], 0))
                signals.append(np.nanmean(signal[choose_transients], 0))

                signal_select.append(signal[choose_transients])
                
            #Cross-validate!
            errs.append(np.mean([ut.rmse(models[0], signals[1]),
                                 ut.rmse(models[1], signals[0])]))
        # We really only need to look at the first two:
        signal_err = ut.rmse(np.nanmean(signal_select[0], 0),
                             np.nanmean(signal_select[1], 0))
        # Based on the errors, choose a function. Also report errors:
        return (fitters[np.argmin(errs)], funcs[np.argmin(errs)], np.min(errs),
                signal_err)
1278048
def _read_or_calc_samples(
    sampler,
    modelidx=0,
    n_samples=100,
    last_step=False,
    e_range=None,
    e_npoints=100,
    threads=None,
):
    

    if not e_range:
        # return the results saved in blobs
        modelx, model = _process_blob(sampler, modelidx, last_step=last_step)
    else:
        # prepare bogus data for calculation
        e_range = validate_array(
            "e_range", u.Quantity(e_range), physical_type="energy"
        )
        e_unit = e_range.unit
        energy = (
            np.logspace(
                np.log10(e_range[0].value),
                np.log10(e_range[1].value),
                e_npoints,
            )
            * e_unit
        )
        data = {
            "energy": energy,
            "flux": np.zeros(energy.shape) * sampler.data["flux"].unit,
        }
        # init pool and select parameters
        chain = sampler.chain[-1] if last_step else sampler.flatchain
        pars = chain[np.random.randint(len(chain), size=n_samples)]
        blobs = []

        p = Pool(threads)
        modelouts = p.map(partial(sampler.modelfn, data=data), pars)
        p.close()
        p.terminate()

        for modelout in modelouts:
            if isinstance(modelout, np.ndarray):
                blobs.append([modelout])
            else:
                blobs.append(modelout)

        modelx, model = _process_blob(
            blobs, modelidx=modelidx, energy=data["energy"]
        )

    return modelx, model
413933
def set_rf_samples(n):
    
    forest._generate_sample_indices = (lambda rs, n_samples:
        forest.check_random_state(rs).randint(0, n_samples, n))
109895
public void cancelStackTraceSample(int sampleId, Throwable cause) {
		synchronized (lock) {
			if (isShutDown) {
				return;
			}

			PendingStackTraceSample sample = pendingSamples.remove(sampleId);
			if (sample != null) {
				if (cause != null) {
					LOG.info("Cancelling sample " + sampleId, cause);
				} else {
					LOG.info("Cancelling sample {}", sampleId);
				}

				sample.discard(cause);
				rememberRecentSampleId(sampleId);
			}
		}
	}
470664
def sample_given_context(self, c, c_dims):
        '''
        
        '''
        index = self.discrete_progress.sample_given_context(c, c_dims, self.space)
        return self.space.rand_value(index).flatten()[list(set(range(len(self.space.cardinalities))) - set(c_dims))]
37609
def xstep(self):
        
        """

        # This test reflects empirical evidence that two slightly
        # different implementations are faster for single or
        # multi-channel data. This kludge is intended to be temporary.
        if self.cri.Cd > 1:
            for i in range(self.Nb):
                self.xistep(i)
        else:
            self.YU[:] = self.Y[..., np.newaxis] - self.U
            b = np.swapaxes(self.ZSf[..., np.newaxis], self.cri.axisK, -1) \
                + self.rho*sl.rfftn(self.YU, None, self.cri.axisN)
            for i in range(self.Nb):
                self.Xf[..., i] = sl.solvedbi_sm(
                    self.Zf[..., [i], :], self.rho, b[..., i],
                    axis=self.cri.axisM)
            self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)


        if self.opt['LinSolveCheck']:
            ZSfs = np.sum(self.ZSf, axis=self.cri.axisK, keepdims=True)
            YU = np.sum(self.Y[..., np.newaxis] - self.U, axis=-1)
            b = ZSfs + self.rho*sl.rfftn(YU, None, self.cri.axisN)
            Xf = self.swapaxes(self.Xf)
            Zop = lambda x: sl.inner(self.Zf, x, axis=self.cri.axisM)
            ZHop = lambda x: np.conj(self.Zf) * x
            ax = np.sum(ZHop(Zop(Xf)) + self.rho*Xf, axis=self.cri.axisK,
                        keepdims=True)
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
1350884
def initGP(self,fast=False):
        
        if fast:
            assert self.n_terms==2, 'CVarianceDecomposition: for fast inference number of terms must be == 2'
            assert self.P>1,        'CVarianceDecomposition: for fast inference number of traits must be > 1'
            self.vd.initGPkronSum()
        else:
            self.vd.initGP()
        self.gp=self.vd.getGP()
        self.init=True
        self.fast=fast
451195
def setSampleTypes(self, value, **kw):
        
        bsc = getToolByName(self, 'bika_setup_catalog')
        # convert value to objects
        if value and type(value) == str:
            value = [bsc(UID=value)[0].getObject(), ]
        elif value and type(value) in (list, tuple) and type(value[0]) == str:
            value = [bsc(UID=uid)[0].getObject() for uid in value if uid]
        if not type(value) in (list, tuple):
            value = [value, ]
        # Find all SampleTypes that were removed
        existing = self.Schema()['SampleTypes'].get(self)
        removed = existing and [s for s in existing if s not in value] or []
        added = value and [s for s in value if s not in existing] or []
        ret = self.Schema()['SampleTypes'].set(self, value)

        # finally be sure that we aren't trying to set None values here.
        removed = [x for x in removed if x]
        added = [x for x in added if x]

        for st in removed:
            samplepoints = st.getSamplePoints()
            if self in samplepoints:
                samplepoints.remove(self)
                st.setSamplePoints(samplepoints)

        for st in added:
            st.setSamplePoints(list(st.getSamplePoints()) + [self, ])

        return ret
1220426
public static void sampleWithoutReplacement(int[] array, int numArgClasses, Random rand) {
    int[] temp = new int[numArgClasses];
    for (int i = 0; i < temp.length; i++) {
      temp[i] = i;
    }
    shuffle(temp, rand);
    System.arraycopy(temp, 0, array, 0, array.length);
  }
1015282
def Downsampled(cls, stats, interval=None):
    
    interval = interval or cls.DEFAULT_SAMPLING_INTERVAL

    result = cls(stats)
    result.cpu_samples = cls._Downsample(
        kind=CpuSample, samples=stats.cpu_samples, interval=interval)
    result.io_samples = cls._Downsample(
        kind=IOSample, samples=stats.io_samples, interval=interval)
    return result
645723
def infer(self, data, initial_proposal=None, full_output=False,**kwargs):

        

        # Apply data masks now so we don't have to do it on the fly.
        data, pixels_affected = self._apply_data_mask(data)

        # Any channels / parameters to ignore?
        matched_channels, missing_channels, ignore_parameters \
            = self._match_channels_to_data(data)
        parameters = [p for p in self.parameters if p not in ignore_parameters]
        #parameters = list(set(self.parameters).difference(ignore_parameters))

        logger.debug("Inferring {0} parameters: {1}".format(len(parameters),
            ", ".join(parameters)))

        # What sampling behaviour will we have?
        # - Auto-convergence:
        #       + Sample for `minimum_sample` (default 2000, 200 walkers)
        #       + Calculate the maximum exponential autocorrelation time for
        #         all parameters
        #       + For the rest of the chain, calculate the autocorrelation time
        #       + Ensure that the number of samples we have is more than 
        #         `effectively_independent_samples` (default 100) times.
        # - Specified convergence:
        #       + Burn for `burn` (default 2000) steps
        #       + Sample for `sample` (default 2000) steps

        kwd = {
            "auto_convergence": False, # TODO CHANGE ME
            "walkers": 100,
            "burn": 2000,
            "sample": 2000,
            # The minimum_sample, n_tau_exp_as_burn_in, minimum_eis are only
            # used if auto_convergence is turned on.
            "minimum_sample": 2000,
            "maximum_sample": 100000,
            "n_tau_exp_as_burn_in": 3,
            "minimum_effective_independent_samples": 100,
            "check_convergence_frequency": 1000,
            "a": 2.0,
            "threads": 1
        }

        # Update from the model, then update from any keyword arguments given.
        kwd.update(self._configuration.get("infer", {}).copy())
        kwd.update(**kwargs)

        # Make some checks.
        if kwd["walkers"] % 2 > 0 or kwd["walkers"] < 2 * len(parameters):
            raise ValueError("the number of walkers must be an even number and "
                "be at least twice the number of model parameters")

        check_keywords = ["threads", "a"]
        if kwd["auto_convergence"]:
            logger.info("Convergence will be estimated automatically.")
            check_keywords += ["minimum_sample", "check_convergence_frequency",
                "minimum_effective_independent_samples", "n_tau_exp_as_burn_in",
                "maximum_sample"]
        
        else:
            check_keywords += ["burn", "sample"]
            logger.warn("No convergence checks will be done!")
            logger.info("Burning for {0} steps and sampling for {1} with {2} "\
                "walkers".format(kwd["burn"], kwd["sample"], kwd["walkers"]))

        for keyword in check_keywords:
            if kwd[keyword] < 1:
                raise ValueError("keyword {} must be a positive value".format(
                    keyword))

        # Check for non-standard proposal scales.
        if kwd["a"] != 2.0:
            logger.warn("Using proposal scale of {0:.2f}".format(kwd["a"]))

        # If no initial proposal given, estimate the model parameters.
        if initial_proposal is None:
            initial_proposal = self.estimate(data)

        # Initial proposal could be:
        #   - an array (N_walkers, N_dimensions)
        #   - a dictionary containing key/value pairs for the dimensions
        if isinstance(initial_proposal, dict):

            wavelengths_required = []
            for channel, spectrum in zip(matched_channels, data):
                if channel is None: continue
                z = initial_proposal.get("z",
                    initial_proposal.get("z_{}".format(channel), 0))
                wavelengths_required.append(
                    [spectrum.disp[0] * (1 - z), spectrum.disp[-1] * (1 - z)])

            closest_point = [initial_proposal[p] \
                for p in self.grid_points.dtype.names]
            subset_bounds = self._initialise_approximator(
                closest_point=closest_point, 
                wavelengths_required=wavelengths_required, force=True, **kwargs)

            initial_proposal = self._initial_proposal_distribution(
                parameters, initial_proposal, kwd["walkers"])

        elif isinstance(initial_proposal, np.ndarray):
            initial_proposal = np.atleast_2d(initial_proposal)
            if initial_proposal.shape != (kwd["walkers"], len(parameters)):
                raise ValueError("initial proposal must be an array of shape "\
                    "(N_parameters, N_walkers) ({0}, {1})".format(kwd["walkers"],
                        len(parameters)))

        # Prepare the convolution functions.
        self._create_convolution_functions(matched_channels, data, parameters)

        # Create the sampler.
        logger.info("Creating sampler with {0} walkers and {1} threads".format(
            kwd["walkers"], kwd["threads"]))
        debug = kwargs.get("debug", False)
        sampler = emcee.EnsembleSampler(kwd["walkers"], len(parameters),
            inference.ln_probability, a=kwd["a"], threads=kwd["threads"],
            args=(parameters, self, data, debug),
            kwargs={"matched_channels": matched_channels})

        # Regardless of whether we automatically check for convergence or not,
        # we will still need to burn in for some minimum amount of time.
        if kwd["auto_convergence"]:
            # Sample for `minimum_sample` period.
            descr, iterations = "", kwd["minimum_sample"]
        else:
            # Sample for `burn` period
            descr, iterations = "burn-in", kwd["burn"]

        # Start sampling.
        t_init = time()
        acceptance_fractions = []
        progress_bar = kwargs.get("__show_progress_bar", True)
        sampler, init_acceptance_fractions, pos, lnprob, rstate, init_elapsed \
            = self._sample(sampler, initial_proposal, iterations, descr=descr,
                parameters=parameters, __show_progress_bar=progress_bar)
        acceptance_fractions.append(init_acceptance_fractions)

        # If we don't have to check for convergence, it's easy:
        if not kwd["auto_convergence"]:

            # Save the chain and log probabilities before we reset the chain.
            burn, sample = kwd["burn"], kwd["sample"]
            converged = None # we don't know!
            burn_chains = sampler.chain
            burn_ln_probabilities = sampler.lnprobability

            # Reset the chain.
            logger.debug("Resetting chain...")
            sampler.reset()

            # Sample the posterior.
            sampler, prod_acceptance_fractions, pos, lnprob, rstate, t_elapsed \
                = self._sample(sampler, pos, kwd["sample"], lnprob0=lnprob,
                    rstate0=rstate, descr="production", parameters=parameters,
                    __show_progress_bar=progress_bar)

            production_chains = sampler.chain
            production_ln_probabilities = sampler.lnprobability
            acceptance_fractions.append(prod_acceptance_fractions)

        else:

            # Start checking for convergence at a frequency
            # of check_convergence_frequency
            last_state = [pos, lnprob, rstate]
            converged, total_steps = False, 0 + iterations
            min_eis_required = kwd["minimum_effective_independent_samples"]
            while not converged and kwd["maximum_sample"] > total_steps:
                
                # Check for convergence.
                # Estimate the exponential autocorrelation time.
                try:
                    tau_exp, rho, rho_max_fit \
                        = utils.estimate_tau_exp(sampler.chain)

                except:
                    logger.exception("Exception occurred when trying to "
                        "estimate the exponential autocorrelation time:")

                    logger.info("To recover, we are temporarily setting tau_exp"
                        " to {0}".format(total_steps))
                    tau_exp = total_steps

                logger.info("Estimated tau_exp at {0} is {1:.0f}".format(
                    total_steps, tau_exp))

                # Grab everything n_tau_exp_as_burn_in times that.
                burn = int(np.ceil(tau_exp)) * kwd["n_tau_exp_as_burn_in"]
                sample = sampler.chain.shape[1] - burn

                if 1 > sample:
                    logger.info("Sampler has not converged because {0}x the "
                        "estimated exponential autocorrelation time of {1:.0f}"
                        " is step {2}, and we are only at step {3}".format(
                            kwd["n_tau_exp_as_burn_in"], tau_exp, burn,
                            total_steps))
        
                else:

                    # Calculate the integrated autocorrelation time in the 
                    # remaining sample, for every parameter.
                    tau_int = utils.estimate_tau_int(sampler.chain[:, burn:])

                    # Calculate the effective number of independent samples in 
                    # each parameter.
                    num_effective = (kwd["walkers"] * sample)/(2*tau_int)
                    
                    logger.info("Effective number of independent samples in "
                        "each parameter:")
                    for parameter, n_eis in zip(parameters, num_effective):
                        logger.info("\t{0}: {1:.0f}".format(parameter, n_eis))

                    if num_effective.min() > min_eis_required:
                        # Converged.
                        converged = True
                        logger.info("Convergence achieved ({0:.0f} > {1:.0f})"\
                            .format(num_effective.min() > min_eis_required))

                        # Separate the samples into burn and production..
                        burn_chains = sampler.chain[:, :burn, :]
                        burn_ln_probabilities = sampler.lnprobability[:burn]

                        production_chains = sampler.chain[:, burn:, :]
                        production_ln_probabilities = sampler.lnprobability[burn:]
                        break

                    else:
                        # Nope.
                        logger.info("Sampler has not converged because it did "
                            "not meet the minimum number of effective "
                            "independent samples ({0:.0f})".format(kwd["n"]))
                
                # Keep sampling.
                iterations = kwd["check_convergence_frequency"]
                logger.info("Trying for another {0} steps".format(iterations))

                pos, lnprob, rstate = last_state
                sampler, af, pos, lnprob, rstate, t_elapsed = self._sample(
                    sampler, pos, iterations, lnprob0=lnprob, rstate0=rstate,
                    descr="", parameters=parameters,
                    __show_progress_bar=progress_bar)

                total_steps += iterations
                acceptance_fractions.append(af)
                last_state.extend(pos, lnprob, rstate)
                del last_state[:3]

            if not converged:
                logger.warn("Maximum number of samples ({:.0f}) reached without"
                    "convergence!".format(kwd["maximum_sample"]))

        logger.info("Total time elapsed: {0} seconds".format(time() - t_init))
        
        if sampler.pool:
            sampler.pool.close()
            sampler.pool.join()

        # Stack burn and production information together.
        chains = np.hstack([burn_chains, production_chains])
        lnprobability = np.hstack([
            burn_ln_probabilities, production_ln_probabilities])
        acceptance_fractions = np.hstack(acceptance_fractions)

        chi_sq, dof, model_fluxes = self._chi_sq(dict(zip(parameters, 
            [np.percentile(chains[:, burn:, i], 50) 
                for i in range(len(parameters))])), data)

        # Convert velocity scales.
        symbol, scale, units = self._preferred_redshift_scale
        labels = [] + parameters
        scales = np.ones(len(parameters))
        if symbol != "z":
            for i, parameter in enumerate(parameters):
                if parameter == "z" or parameter.startswith("z_"):
                    chains[:, :, i] *= scale
                    scales[i] = scale
                    if "_" in parameter:
                        labels[i] = "_".join([symbol, parameter.split("_")[1:]])
                    else:
                        labels[i] = symbol
                    logger.debug("Scaled {0} (now {1}) to units of {2}".format(
                        parameter, labels[i], units))

        # Calculate MAP values and associated uncertainties.
        theta = OrderedDict()
        for i, label in enumerate(labels):
            l, c, u = np.percentile(chains[:, burn:, i], [16, 50, 84])
            theta[label] = (c, u-c, l-c)

        # Re-arrange the chains to be in the same order as the model parameters.
        indices = np.array([parameters.index(p) \
            for p in self.parameters if p in parameters])
        chains = chains[:, :, indices]

        # Remove the convolution functions.
        if not kwargs.get("__keep_convolution_functions", False):
            self._destroy_convolution_functions()

        if full_output:
            metadata = {
                "burn": burn,
                "walkers": kwd["walkers"],
                "sample": sample,
                "parameters": labels,
                "scales": scales,
                "chi_sq": chi_sq,
                "dof": dof
            }
            return (theta, chains, lnprobability, acceptance_fractions, sampler,
                metadata)
        return theta
843652
private static boolean sample(Set<Long> rootIds) {
        if (Double.compare(sampleRate, 1.0d) >= 0)
            return true;
        int threshold = (int) (sampleRate * PRECISION);
        for (Long id : rootIds) {
            int mod = (int) (Math.abs(id) % PRECISION);
            if (mod < threshold) {
                return true;
            }
        }
        return false;
    }