Provides a stochastic estimate of the variational lower bound,
       for any variational family and model density.
422194
def variational_lower_bound(params, t, logprob, sampler, log_density,
                            num_samples, rs):
    """"""
    samples = sampler(params, num_samples, rs)
    log_qs = log_density(params, samples)
    log_ps = logprob(samples, t)
    log_ps = np.reshape(log_ps, (num_samples, -1))
    log_qs = np.reshape(log_qs, (num_samples, -1))
    return np.mean(log_ps - log_qs)
495777
def _update_estimate_and_sampler(self, ell, ell_hat, weight, extra_info,
                                     **kwargs):
        """"""
        self._TP += ell_hat * ell * weight
        self._FP += ell_hat * (1 - ell) * weight
        self._FN += (1 - ell_hat) * ell * weight
        self._TN += (1 - ell_hat) * (1 - ell) * weight

        self._estimate[self.t_] = \
                self._F_measure(self.alpha, self._TP, self._FP, self._FN)
982964
def estimate_parameters(self, max_dist_kb, size_bin_kb, display_graph):
        
        logger.info("estimation of the parameters of the model")
        self.bins = np.arange(
            size_bin_kb, max_dist_kb + size_bin_kb, size_bin_kb
        )

        self.mean_contacts = np.zeros_like(self.bins, dtype=np.float32)
        self.dict_collect = dict()
        self.gpu_vect_frags.copy_from_gpu()
        epsi = self.mean_value_trans
        for k in self.bins:
            self.dict_collect[k] = []
        for i in range(0, 2000):
            # print "frag i = ", i
            start = self.sparse_matrix.indptr[i]
            end = self.sparse_matrix.indptr[i + 1]
            id_j = self.sparse_matrix.indices[start:end]
            data = self.sparse_matrix.data[start:end]

            info_i = self.np_sub_frags_2_frags[i]
            init_id_fi = info_i[0]
            # pos_i = self.S_o_A_frags["pos"][init_id_fi]
            id_c_i = self.S_o_A_frags["id_c"][init_id_fi]
            s_i = (
                self.S_o_A_frags["start_bp"][init_id_fi] / 1000.0
                + self.np_sub_frags_2_frags[i][1]
            )
            len_kb_c_i = self.S_o_A_frags["l_cont_bp"][init_id_fi] / 1000
            local_bins = np.arange(
                size_bin_kb,
                min(len_kb_c_i, max_dist_kb) + size_bin_kb,
                size_bin_kb,
            )
            local_storage = np.zeros_like(local_bins, dtype=np.int32)

            for fj, dj in zip(id_j, data):
                info_j = self.np_sub_frags_2_frags[fj]
                init_id_fj = info_j[0]
                id_c_j = self.S_o_A_frags["id_c"][init_id_fj]
                if id_c_i == id_c_j:
                    # pos_j = self.S_o_A_frags["pos"][init_id_fj]
                    s_j = (
                        self.S_o_A_frags["start_bp"][init_id_fj] / 1000.0
                        + self.np_sub_frags_2_frags[fj][1]
                    )
                    d = np.abs(s_i - s_j)
                    if d < max_dist_kb:
                        id_bin = d / size_bin_kb
                        local_storage[id_bin] += dj
                        # self.dict_collect[self.bins[id_bin]].append(dj)
            # we have to add also the zeros
            for my_bin, val in zip(local_bins, local_storage):
                # print "bin = ", bin
                self.dict_collect[my_bin].append(val)

        for id_bin in range(0, len(self.bins)):
            k = self.bins[id_bin]
            self.mean_contacts[id_bin] = np.mean(self.dict_collect[k])

        for id_bin in range(0, len(self.bins)):
            k = self.bins[id_bin]
            tmp = np.mean(self.dict_collect[k])
            if np.isnan(tmp) or tmp == 0:
                # if np.isnan(tmp):
                # if np.isnan(tmp):
                #     print "removing nan"
                self.mean_contacts[id_bin] = np.nan
            else:
                self.mean_contacts[id_bin] = tmp + epsi

        self.mean_contacts_upd = []
        self.bins_upd = []

        for count, ele in zip(self.mean_contacts, self.bins):
            if not np.isnan(count):
                self.bins_upd.append(ele)
                self.mean_contacts_upd.append(count)

        self.bins_upd = np.array(self.bins_upd)
        self.mean_contacts_upd = np.array(self.mean_contacts_upd)

        # self.mean_contacts_upd =
        # ndi.filters.gaussian_filter1d(self.mean_contacts_upd,
        # sigma=len(self.bins_upd) / 5.)

        p, self.y_estim = nuis.estimate_param_hic(
            self.mean_contacts_upd, self.bins_upd
        )
        ##########################################
        fit_param = p.x
        ##########################################
        logger.info("mean value trans = {}".format(self.mean_value_trans))
        ##########################################
        estim_max_dist = nuis.estimate_max_dist_intra(
            fit_param, self.mean_value_trans
        )
        logger.info("max distance cis/trans = {}".format(estim_max_dist))
        ##########################################
        self.param_simu = self.setup_model_parameters(
            fit_param, estim_max_dist
        )
        self.gpu_param_simu = cuda.mem_alloc(self.param_simu.nbytes)
        self.gpu_param_simu_test = cuda.mem_alloc(self.param_simu.nbytes)

        cuda.memcpy_htod(self.gpu_param_simu, self.param_simu)

        if display_graph:
            plt.loglog(self.bins_upd, self.mean_contacts_upd, "-*b")
            plt.loglog(self.bins_upd, self.y_estim, "-*r")
            plt.xlabel("genomic distance (kb)")
            plt.ylabel("frequency of contact")
            plt.legend(["obs", "fit"])
            plt.show()
429156
def free_param_bounds(self):
        
        fpb = CombinedBounds(self.k.free_param_bounds, self.noise_k.free_param_bounds)
        if self.mu is not None:
            fpb = CombinedBounds(fpb, self.mu.free_param_bounds)
        return fpb
623265
def inference(self, kern_r, kern_c, Xr, Xc, Zr, Zc, likelihood, Y, qU_mean ,qU_var_r, qU_var_c, indexD, output_dim):
        

        N, D, Mr, Mc, Qr, Qc = Y.shape[0], output_dim,Zr.shape[0], Zc.shape[0], Zr.shape[1], Zc.shape[1]

        uncertain_inputs_r = isinstance(Xr, VariationalPosterior)
        uncertain_inputs_c = isinstance(Xc, VariationalPosterior)
        uncertain_outputs = isinstance(Y, VariationalPosterior)

        grad_dict = self._init_grad_dict(N,D,Mr,Mc)

        beta = 1./likelihood.variance
        if len(beta)==1:
            beta = np.zeros(D)+beta

        psi0_r, psi1_r, psi2_r = self.gatherPsiStat(kern_r, Xr, Zr, uncertain_inputs_r)
        psi0_c, psi1_c, psi2_c = self.gatherPsiStat(kern_c, Xc, Zc, uncertain_inputs_c)

        #======================================================================
        # Compute Common Components
        #======================================================================

        Kuu_r = kern_r.K(Zr).copy()
        diag.add(Kuu_r, self.const_jitter)
        Lr = jitchol(Kuu_r)

        Kuu_c = kern_c.K(Zc).copy()
        diag.add(Kuu_c, self.const_jitter)
        Lc = jitchol(Kuu_c)

        mu, Sr, Sc = qU_mean, qU_var_r, qU_var_c
        LSr = jitchol(Sr)
        LSc = jitchol(Sc)

        LcInvMLrInvT = dtrtrs(Lc,dtrtrs(Lr,mu.T)[0].T)[0]
        LcInvLSc = dtrtrs(Lc, LSc)[0]
        LrInvLSr = dtrtrs(Lr, LSr)[0]
        LcInvScLcInvT = tdot(LcInvLSc)
        LrInvSrLrInvT = tdot(LrInvLSr)
        tr_LrInvSrLrInvT = np.square(LrInvLSr).sum()
        tr_LcInvScLcInvT = np.square(LcInvLSc).sum()

        mid_res = {
            'psi0_r': psi0_r,
            'psi1_r': psi1_r,
            'psi2_r': psi2_r,
            'psi0_c': psi0_c,
            'psi1_c': psi1_c,
            'psi2_c': psi2_c,
            'Lr':Lr,
            'Lc':Lc,
            'LcInvMLrInvT': LcInvMLrInvT,
            'LcInvScLcInvT': LcInvScLcInvT,
            'LrInvSrLrInvT': LrInvSrLrInvT,
        }

        #======================================================================
        # Compute log-likelihood
        #======================================================================

        logL = 0.
        for d in range(D):
            logL += self.inference_d(d, beta, Y, indexD, grad_dict, mid_res, uncertain_inputs_r, uncertain_inputs_c, Mr, Mc)

        logL += -Mc * (np.log(np.diag(Lr)).sum()-np.log(np.diag(LSr)).sum())  -Mr * (np.log(np.diag(Lc)).sum()-np.log(np.diag(LSc)).sum()) \
               - np.square(LcInvMLrInvT).sum()/2. - tr_LrInvSrLrInvT * tr_LcInvScLcInvT/2. + Mr*Mc/2.

        #======================================================================
        # Compute dL_dKuu
        #======================================================================

        tmp =  tdot(LcInvMLrInvT)/2. + tr_LrInvSrLrInvT/2. * LcInvScLcInvT - Mr/2.*np.eye(Mc)

        dL_dKuu_c = backsub_both_sides(Lc, tmp, 'left')
        dL_dKuu_c += dL_dKuu_c.T
        dL_dKuu_c *= 0.5

        tmp =  tdot(LcInvMLrInvT.T)/2. + tr_LcInvScLcInvT/2. * LrInvSrLrInvT - Mc/2.*np.eye(Mr)

        dL_dKuu_r = backsub_both_sides(Lr, tmp, 'left')
        dL_dKuu_r += dL_dKuu_r.T
        dL_dKuu_r *= 0.5

        #======================================================================
        # Compute dL_dqU
        #======================================================================

        tmp = - LcInvMLrInvT
        dL_dqU_mean = dtrtrs(Lc, dtrtrs(Lr, tmp.T, trans=1)[0].T, trans=1)[0]

        LScInv = dtrtri(LSc)
        tmp = -tr_LrInvSrLrInvT/2.*np.eye(Mc)
        dL_dqU_var_c = backsub_both_sides(Lc, tmp, 'left') + tdot(LScInv.T) * Mr/2.

        LSrInv = dtrtri(LSr)
        tmp =  -tr_LcInvScLcInvT/2.*np.eye(Mr)
        dL_dqU_var_r = backsub_both_sides(Lr, tmp, 'left') + tdot(LSrInv.T) * Mc/2.

        #======================================================================
        # Compute the Posterior distribution of inducing points p(u|Y)
        #======================================================================

        post = PosteriorMultioutput(LcInvMLrInvT=LcInvMLrInvT, LcInvScLcInvT=LcInvScLcInvT,
                LrInvSrLrInvT=LrInvSrLrInvT, Lr=Lr, Lc=Lc, kern_r=kern_r, Xr=Xr, Zr=Zr)

        #======================================================================
        # Compute dL_dpsi
        #======================================================================

        grad_dict['dL_dqU_mean'] += dL_dqU_mean
        grad_dict['dL_dqU_var_c'] += dL_dqU_var_c
        grad_dict['dL_dqU_var_r'] += dL_dqU_var_r
        grad_dict['dL_dKuu_c'] += dL_dKuu_c
        grad_dict['dL_dKuu_r'] += dL_dKuu_r

        if not uncertain_inputs_c:
            grad_dict['dL_dKdiag_c'] = grad_dict['dL_dpsi0_c']
            grad_dict['dL_dKfu_c'] = grad_dict['dL_dpsi1_c']

        if not uncertain_inputs_r:
            grad_dict['dL_dKdiag_r'] = grad_dict['dL_dpsi0_r']
            grad_dict['dL_dKfu_r'] = grad_dict['dL_dpsi1_r']

        return post, logL, grad_dict
1222940
def noise_op(latents, hparams):
  
  if hparams.latent_noise == 0 or hparams.mode != tf.estimator.ModeKeys.TRAIN:
    return latents
  latent_shape = common_layers.shape_list(latents)
  return latents + tf.random_normal(latent_shape, stddev=hparams.latent_noise)
634509
def compute_lower_limit(mu_in, post, alpha=0.9):
    
    if 0 < alpha < 1:
        dp = integral_element(mu_in, post)
        low_idx = bisect.bisect_right(dp.cumsum() / dp.sum(), 1 - alpha)
        # if alpha is in [0,1) and post is non-negative, bisect_right
        # will always return an index in the range of mu since
        # post.cumsum()/post.sum() will always begin at 0 and end at 1
        mu_low = mu_in[low_idx]
    elif alpha == 1:
        mu_low = numpy.min(mu_in[post > 0])
    else:
        raise ValueError("Confidence level must be in (0,1].")

    return mu_low
676385
def estimateHeritabilities(self, K, verbose=False):
        
        
        # Fit single trait model
        varg  = SP.zeros(self.P)
        varn  = SP.zeros(self.P)
        fixed = SP.zeros((1,self.P))
        
        for p in range(self.P):
            y = self.Y[:,p:p+1]
            lmm = limix.CLMM()
            lmm.setK(K)
            lmm.setSNPs(SP.ones((K.shape[0],1)))
            lmm.setPheno(y)
            lmm.setCovs(SP.zeros((K.shape[0],1)))
            lmm.setVarcompApprox0(-20, 20, 1000)
            lmm.process()
            delta = SP.exp(lmm.getLdelta0()[0,0])
            Vtot  = SP.exp(lmm.getLSigma()[0,0])
    
            varg[p] = Vtot
            varn[p] = delta*Vtot
            fixed[:,p] = lmm.getBetaSNP()
    
            if verbose: print(p)
    
        sth = {}
        sth['varg']  = varg
        sth['varn']  = varn
        sth['fixed'] = fixed

        return sth
167350
function inferSettings(graph) {
  var order = graph.order;

  return {
    barnesHutOptimize: order > 2000,
    strongGravityMode: true,
    gravity: 0.05,
    scalingRatio: 10,
    slowDown: 1 + Math.log(order)
  };
}
1151256
function estimateRootSize(seriesModel, targetInfo, viewRoot, containerWidth, containerHeight) {
    // If targetInfo.node exists, we zoom to the node,
    // so estimate whold width and heigth by target node.
    var currNode = (targetInfo || {}).node;
    var defaultSize = [containerWidth, containerHeight];

    if (!currNode || currNode === viewRoot) {
        return defaultSize;
    }

    var parent;
    var viewArea = containerWidth * containerHeight;
    var area = viewArea * seriesModel.option.zoomToNodeRatio;

    while (parent = currNode.parentNode) { // jshint ignore:line
        var sum = 0;
        var siblings = parent.children;

        for (var i = 0, len = siblings.length; i < len; i++) {
            sum += siblings[i].getValue();
        }
        var currNodeValue = currNode.getValue();
        if (currNodeValue === 0) {
            return defaultSize;
        }
        area *= sum / currNodeValue;

        // Considering border, suppose aspect ratio is 1.
        var parentModel = parent.getModel();
        var borderWidth = parentModel.get(PATH_BORDER_WIDTH);
        var upperHeight = Math.max(borderWidth, getUpperLabelHeight(parentModel, borderWidth));
        area += 4 * borderWidth * borderWidth
            + (3 * borderWidth + upperHeight) * Math.pow(area, 0.5);

        area > MAX_SAFE_INTEGER && (area = MAX_SAFE_INTEGER);

        currNode = parent;
    }

    area < viewArea && (area = viewArea);
    var scale = Math.pow(area / viewArea, 0.5);

    return [containerWidth * scale, containerHeight * scale];
}
1076524
def _hessfcn(self, x, lmbda):
        
        Pgen = x[self._Pg.i1:self._Pg.iN + 1] # Active generation in p.u.
        Qgen = x[self._Qg.i1:self._Qg.iN + 1] # Reactive generation in p.u.

        for i, g in enumerate(self._gn):
            g.p = Pgen[i] * self._base_mva # active generation in MW
            g.q = Qgen[i] * self._base_mva # reactive generation in MVAr

        Vang = x[self._Va.i1:self._Va.iN + 1]
        Vmag = x[self._Vm.i1:self._Vm.iN + 1]
        V = Vmag * exp(1j * Vang)
        nxtra = self._nxyz - 2 * self._nb

        #------------------------------------------------------------------
        #  Evaluate d2f.
        #------------------------------------------------------------------

        d2f = self._d2f(x) * self.opt["cost_mult"]
        # TODO: Generalised cost model.

        #------------------------------------------------------------------
        #  Evaluate Hessian of power balance constraints.
        #------------------------------------------------------------------

        nlam = len(lmbda["eqnonlin"]) / 2
        lamP = lmbda["eqnonlin"][:nlam]
        lamQ = lmbda["eqnonlin"][nlam:nlam + nlam]
        Gpaa, Gpav, Gpva, Gpvv = self.om.case.d2Sbus_dV2(self._Ybus, V, lamP)
        Gqaa, Gqav, Gqva, Gqvv = self.om.case.d2Sbus_dV2(self._Ybus, V, lamQ)

        d2G = vstack([
            hstack([
                vstack([hstack([Gpaa, Gpav]),
                        hstack([Gpva, Gpvv])]).real +
                vstack([hstack([Gqaa, Gqav]),
                        hstack([Gqva, Gqvv])]).imag,
                csr_matrix((2 * self._nb, nxtra))]),
            hstack([
                csr_matrix((nxtra, 2 * self._nb)),
                csr_matrix((nxtra, nxtra))
            ])
        ], "csr")

        #------------------------------------------------------------------
        #  Evaluate Hessian of flow constraints.
        #------------------------------------------------------------------

        nmu = len(lmbda["ineqnonlin"]) / 2
        muF = lmbda["ineqnonlin"][:nmu]
        muT = lmbda["ineqnonlin"][nmu:nmu + nmu]
        if self.flow_lim == "I":
            dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It = \
                self.om.case.dIbr_dV(self._Yf, self._Yt, V)
            Hfaa, Hfav, Hfva, Hfvv = \
                self.om.case.d2AIbr_dV2(dIf_dVa, dIf_dVm, If, self._Yf, V, muF)
            Htaa, Htav, Htva, Htvv = \
                self.om.case.d2AIbr_dV2(dIt_dVa, dIt_dVm, It, self._Yt, V, muT)
        else:
            f = [e.from_bus._i for e in self._ln]
            t = [e.to_bus._i for e in self._ln]
            # Line-bus connection matrices.
            Cf = csr_matrix((ones(self._nl), (range(self._nl), f)), (self._nl, self._nb))
            Ct = csr_matrix((ones(self._nl), (range(self._nl), t)), (self._nl, self._nb))
            dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St = \
                self.om.case.dSbr_dV(self._Yf, self._Yt, V)
            if self.flow_lim == PFLOW:
                Hfaa, Hfav, Hfva, Hfvv = \
                    self.om.case.d2ASbr_dV2(dSf_dVa.real(), dSf_dVm.real(),
                                            Sf.real(), Cf, self._Yf, V, muF)
                Htaa, Htav, Htva, Htvv = \
                    self.om.case.d2ASbr_dV2(dSt_dVa.real(), dSt_dVm.real(),
                                            St.real(), Ct, self._Yt, V, muT)
            elif self.flow_lim == SFLOW:
                Hfaa, Hfav, Hfva, Hfvv = \
                    self.om.case.d2ASbr_dV2(
                        dSf_dVa, dSf_dVm, Sf, Cf, self._Yf, V, muF)
                Htaa, Htav, Htva, Htvv = \
                    self.om.case.d2ASbr_dV2(
                        dSt_dVa, dSt_dVm, St, Ct, self._Yt, V, muT)
            else:
                raise ValueError

        d2H = vstack([
            hstack([
                vstack([hstack([Hfaa, Hfav]),
                        hstack([Hfva, Hfvv])]) +
                vstack([hstack([Htaa, Htav]),
                        hstack([Htva, Htvv])]),
                csr_matrix((2 * self._nb, nxtra))
            ]),
            hstack([
                csr_matrix((nxtra, 2 * self._nb)),
                csr_matrix((nxtra, nxtra))
            ])
        ], "csr")

        return d2f + d2G + d2H
528739
def estimate(self, points):
        """"""
        points = list(points)
        densities = callMLlibFunc(
            "estimateKernelDensity", self._sample, self._bandwidth, points)
        return np.asarray(densities)
1308294
def posterior(self, x, sigma=1.):
        """"""
        pr0 = 1. / self.scale**2  # prior precision
        prd = x.size / sigma**2  # data precision
        varp = 1. / (pr0 + prd)  # posterior variance
        mu = varp * (pr0 * self.loc + prd * x.mean())
        return Normal(loc=mu, scale=np.sqrt(varp))
745665
def _penalized_acquisition(self, x,  model, X_batch, r_x0, s_x0):
        '''
        
        '''
        fval = -self.acq.acquisition_function(x)[:,0]

        if self.transform=='softplus':
            fval_org = fval.copy()
            fval[fval_org>=40.] = np.log(fval_org[fval_org>=40.])
            fval[fval_org<40.] = np.log(np.log1p(np.exp(fval_org[fval_org<40.])))
        elif self.transform=='none':
            fval = np.log(fval+1e-50)

        fval = -fval
        if X_batch is not None:
            h_vals = self._hammer_function(x, X_batch, r_x0, s_x0)
            fval += -h_vals.sum(axis=-1)
        return fval
907830
def posterior(self, x):
        """"""
        return Gamma(a=self.a + 0.5 * x.size,
                     b=self.b + 0.5 * np.sum(x**2))
221167
def next_hyperparameter_lowest_mu(fun_prediction,
                                  fun_prediction_args,
                                  x_bounds, x_types,
                                  minimize_starting_points,
                                  minimize_constraints_fun=None):
    '''
    
    '''
    best_x = None
    best_acquisition_value = None
    x_bounds_minmax = [[i[0], i[-1]] for i in x_bounds]
    x_bounds_minmax = numpy.array(x_bounds_minmax)

    for starting_point in numpy.array(minimize_starting_points):
        res = minimize(fun=_lowest_mu,
                       x0=starting_point.reshape(1, -1),
                       bounds=x_bounds_minmax,
                       method="L-BFGS-B",
                       args=(fun_prediction, fun_prediction_args, \
                             x_bounds, x_types, minimize_constraints_fun))

        if (best_acquisition_value is None) or (res.fun < best_acquisition_value):
            res.x = numpy.ndarray.tolist(res.x)
            res.x = lib_data.match_val_type(res.x, x_bounds, x_types)
            if (minimize_constraints_fun is None) or (minimize_constraints_fun(res.x) is True):
                best_acquisition_value = res.fun
                best_x = res.x

    outputs = None
    if best_x is not None:
        mu, sigma = fun_prediction(best_x, *fun_prediction_args)
        outputs = {'hyperparameter': best_x, 'expected_mu': mu,
                   'expected_sigma': sigma, 'acquisition_func': "lm"}
    return outputs
953392
public double repairSolutionVariableValue(double value, double lowerBound, double upperBound) {
    if (lowerBound > upperBound) {
      throw new JMetalException("The lower bound (" + lowerBound + ") is greater than the "
          + "upper bound (" + upperBound+")") ;
    }
    double result = value ;
    if (value < lowerBound) {
      result = randomGenerator.getRandomValue(lowerBound, upperBound) ;
    }
    if (value > upperBound) {
      result = randomGenerator.getRandomValue(lowerBound, upperBound) ;
    }

    return result ;
  }