The first phase of inference:
        Compute: log-likelihood, dL_dKmm

        Cached intermediate results: Kmm, KmmInv,
518376
def inference_likelihood(self, kern, X, Z, likelihood, Y):
        

        num_data, output_dim = Y.shape
        input_dim = Z.shape[0]
        if self.mpi_comm != None:
            from mpi4py import MPI
            num_data_all = np.array(num_data,dtype=np.int32)
            self.mpi_comm.Allreduce([np.int32(num_data), MPI.INT], [num_data_all, MPI.INT])
            num_data = num_data_all

        if isinstance(X, VariationalPosterior):
            uncertain_inputs = True
        else:
            uncertain_inputs = False

        #see whether we've got a different noise variance for each datum
        beta = 1./np.fmax(likelihood.variance, 1e-6)
        het_noise = beta.size > 1
        if het_noise:
            self.batchsize = 1

        psi0_full, psi1Y_full, psi2_full, YRY_full = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs)

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

        Kmm = kern.K(Z).copy()
        diag.add(Kmm, self.const_jitter)
        if not np.isfinite(Kmm).all():
            print(Kmm)
        Lm = jitchol(Kmm)
        LmInv = dtrtri(Lm)

        LmInvPsi2LmInvT = LmInv.dot(psi2_full.dot(LmInv.T))
        Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
        LL = jitchol(Lambda)
        LLInv = dtrtri(LL)
        logdet_L = 2.*np.sum(np.log(np.diag(LL)))
        LmLLInv = LLInv.dot(LmInv)
        
        b  = psi1Y_full.dot(LmLLInv.T)
        bbt = np.square(b).sum()
        v = b.dot(LmLLInv).T
        LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)
        
        tmp = -LLInv.T.dot(LLinvPsi1TYYTPsi1LLinvT+output_dim*np.eye(input_dim)).dot(LLInv)
        dL_dpsi2R = LmInv.T.dot(tmp+output_dim*np.eye(input_dim)).dot(LmInv)/2.
        
        # Cache intermediate results
        self.midRes['dL_dpsi2R'] = dL_dpsi2R
        self.midRes['v'] = v

        #======================================================================
        # Compute log-likelihood
        #======================================================================
        if het_noise:
            logL_R = -np.log(beta).sum()
        else:
            logL_R = -num_data*np.log(beta)
        logL = -(output_dim*(num_data*log_2_pi+logL_R+psi0_full-np.trace(LmInvPsi2LmInvT))+YRY_full-bbt)/2.-output_dim*logdet_L/2.

        #======================================================================
        # Compute dL_dKmm
        #======================================================================

        dL_dKmm =  dL_dpsi2R - output_dim*LmInv.T.dot(LmInvPsi2LmInvT).dot(LmInv)/2.

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

        if not self.Y_speedup or het_noise:
            wd_inv = backsub_both_sides(Lm, np.eye(input_dim)- backsub_both_sides(LL, np.identity(input_dim), transpose='left'), transpose='left')
            post = Posterior(woodbury_inv=wd_inv, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=Lm)
        else:
            post = None

        #======================================================================
        # Compute dL_dthetaL for uncertian input and non-heter noise
        #======================================================================

        if not het_noise:
            dL_dthetaL = (YRY_full*beta + beta*output_dim*psi0_full - num_data*output_dim*beta)/2. - beta*(dL_dpsi2R*psi2_full).sum() - beta*(v.T*psi1Y_full).sum()
            self.midRes['dL_dthetaL'] = dL_dthetaL

        return logL, dL_dKmm, post
1344583
def likelihood(x, m=None, Cinv=None, sigma=1, detC=None):
        """"""
        # testing: MC integrate must be one: mean(p(x_i)) * volume(where x_i are uniformely sampled)
        # for i in xrange(3): print mean([cma.likelihood(20*r-10, dim * [0], None, 3) for r in rand(10000,dim)]) * 20**dim
        if m is None:
            dx = x
        else:
            dx = x - m  # array(x) - array(m)
        n = len(x)
        s2pi = (2 * np.pi)**(n / 2.)
        if Cinv is None:
            return exp(-sum(dx**2) / sigma**2 / 2) / s2pi / sigma**n
        if detC is None:
            detC = 1. / np.linalg.linalg.det(Cinv)
        return  exp(-np.dot(dx, np.dot(Cinv, dx)) / sigma**2 / 2) / s2pi / abs(detC)**0.5 / sigma**n
461514
def partition_inference(tctc_mat, comp, tau, sigma, kappa):
     
    """
    communityinfo = {}
    communityinfo['community'] = []
    communityinfo['start'] = np.empty(0)
    communityinfo['end'] = np.empty(0)
    communityinfo['size'] = np.empty(0)
    for i, tcomp in enumerate(comp):
        # This can go in parallel loop
        if len(tcomp) > 0:
            for traj in tcomp:
                # Check it does not already exist.
                ignore = 0
                preexisting = 0
                if i != 0:
                    cutoff = i-1-kappa
                    if cutoff < 0:
                        cutoff = 0
                    
                    if np.any(np.sum(np.sum(tctc_mat[traj, :, cutoff:i][:, traj], axis=0), axis=0) == np.power(len(traj), 2)):
                        # Make sure that a small trajectory could exist 
                        for checknode in np.where(communityinfo['end']>=cutoff)[0]:
                            if traj == communityinfo['community'][checknode]:
                                ignore = 1
                        if ignore == 0: 
                            for checknode in np.where(communityinfo['end']>=cutoff)[0]:
                                if set(communityinfo['community'][checknode]).issuperset(traj):
                                    preexisting = 1                              
                if ignore == 0:
                    # Check how long it continues
                    # For efficiency, increase in blocks
                    approxmaxlength = tau*2
                    a = np.sum(
                        np.sum(tctc_mat[traj, :, i:i+approxmaxlength][:, traj], axis=0), axis=0)
                    if len(traj)*len(traj)*approxmaxlength == a.sum():
                        ok = 0
                        ii = 1
                        while ok == 0:
                            b = np.sum(np.sum(
                                tctc_mat[traj, :, i+(approxmaxlength*ii):i+(approxmaxlength*(ii+1))][:, traj], axis=0), axis=0)
                            a = np.append(a, b)
                            if len(traj)*len(traj)*approxmaxlength != b.sum():
                                ok = 1
                            if i+(approxmaxlength*(ii+1)) > tctc_mat.shape[-1]:
                                ok = 1
                            ii += 1
                    a = np.where(a == np.power(len(traj), 2))[0]
                    # Add an additional value that is false in case end of time series
                    if len(a) == 1:
                        stopind = i + 1
                    else:
                        a = np.append(a, a.max()+kappa+2)
                        # Find the stop index (if stopind = 4 and start = 0, then tctc_mat[:,:,start:stopind]==1)
                        stopind = i + np.split(a, np.where(
                                np.diff(a) > kappa+1)[0]+1)[0][-1] + 1
                    # Add trajectory to dictionary
                    if ((stopind - i) >= tau or preexisting == 1) and len(traj) >= sigma:
                        communityinfo['community'].append(sorted(traj))
                        communityinfo['start'] = np.append(communityinfo['start'], int(i))
                        communityinfo['end'] = np.append(
                            communityinfo['end'], int(stopind))
                        communityinfo['size'] = np.append(communityinfo['size'], len(traj))

    communityinfo = pd.DataFrame(communityinfo)

    communityinfo['start'] = communityinfo['start'].astype(int)
    communityinfo['end'] = communityinfo['end'].astype(int)
    # First check that there is not already a trajectory that is ongoing
    badrows = []
    for v in communityinfo.iterrows():
        skipselrule = (communityinfo['end'] == v[1]['end'])
        for u in communityinfo[skipselrule].iterrows():
            a = 1
            if u[1]['start'] > v[1]['start'] and sorted(u[1]['community']) == sorted(v[1]['community']):
                badrows.append(u[0])
    communityinfo = communityinfo.drop(badrows)

    # Then see if any subset trajectory can be placed earlier in time.
    for v in communityinfo.iterrows():
        skipselrule = (communityinfo['end'] <= v[1]['start']) & (
            communityinfo['end']+kappa >= v[1]['start'])
        for u in communityinfo[skipselrule].iterrows():
            a = 1
            if set(u[1]['community']).issuperset(v[1]['community']):
                communityinfo.loc[v[0], 'start'] = u[1]['start']

    # It is possible to make the condition below effective_length
    communityinfo['length'] = np.array(communityinfo['end']) - np.array(communityinfo['start'])
    communityinfo = communityinfo[communityinfo['length'] >= tau]
    communityinfo = communityinfo[communityinfo['size'] >= sigma]

    # Make sure that the traj is not completely enguled by another
    badrows = []
    if kappa > 0:
        for v in communityinfo.iterrows():
            skipselrule = (communityinfo['end'] == v[1]['end']) & (
                communityinfo['start'] < v[1]['start'])
            for u in communityinfo[skipselrule].iterrows():
                if set(v[1]['community']).issubset(u[1]['community']):
                    badrows.append(v[0])
        communityinfo = communityinfo.drop(badrows)

    return communityinfo
1368338
def NegLnLike(x, time, flux, errors, kernel):
    '''
    

    '''

    gp = GP(kernel, x, white=True)
    gp.compute(time, errors)
    if OLDGEORGE:
        nll = -gp.lnlikelihood(flux)
        # NOTE: There was a bug on this next line! Used to be
        #
        #    ngr = -gp.grad_lnlikelihood(flux) / gp.kernel.pars
        #
        # But I think we want
        #
        # dlogL/dx =     dlogL/dlogx^2       * dlogx^2/dx^2 * dx^2/dx
        #          = gp.grad_lnlikelihood()  *     1/x^2    *   2x
        #          = 2 * gp.grad_lnlikelihood() / x
        #          = 2 * gp.grad_lnlikelihood() / np.sqrt(x^2)
        #          = 2 * gp.grad_lnlikelihood() / np.sqrt(gp.kernel.pars)
        #
        # (with a negative sign out front for the negative gradient).
        # So we probably weren't optimizing the GP correctly! This affects
        # all campaigns through C13. It's not a *huge* deal, since the sign
        # of the gradient was correct and the model isn't that sensitive to
        # the value of the hyperparameters, but it may have contributed to
        # the poor performance on super variable stars. In most cases it means
        # the solver takes longer to converge and isn't as good at finding
        # the minimum.
        ngr = -2 * gp.grad_lnlikelihood(flux) / np.sqrt(gp.kernel.pars)
    else:
        nll = -gp.log_likelihood(flux)
        ngr = -2 * gp.grad_log_likelihood(flux) / \
            np.sqrt(np.exp(gp.get_parameter_vector()))

    return nll, ngr
190907
def expected_log_likelihood(self, mus, sigmas, y):
        

        # Flatten the covariance
        T = mus.shape[0]
        D = self.D_in
        sigs_vec = sigmas.reshape((T, D ** 2))

        # Compute the log likelihood of each column
        ll = np.zeros((T, self.D_out))
        for n in range(self.D_out):

            an = self.A[n]

            E_loglmbda = np.dot(mus, an)
            ll[:,n] += y[:,n] * E_loglmbda

            # Vectorized log likelihood calculation
            aa_vec = np.outer(an, an).reshape((D ** 2,))
            ll[:,n] = -np.exp(E_loglmbda + 0.5 * np.dot(sigs_vec, aa_vec))

        return ll
352592
def dM(self, k, t, param, Mkt, tips=None, gaps=None):
        """"""
        assert 0 <= k < self.ncats
        assert ((param in self.freeparams) or (param == 't') or (
                param == self.distributedparam))
        assert param not in self.distributionparams
        return self._models[k].dM(t, param, Mkt, tips=tips, gaps=gaps)
610643
def log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # prior log p(x) -- quadratic terms
        J_diag, J_upper_diag = self.sparse_J_prior
        lp = -0.5 * np.sum(x * symm_block_tridiag_matmul(J_diag, J_upper_diag, x))

        # prior log p(x) -- linear terms
        _, h_init, log_Z_init = self.info_init_params
        _, _, _, h1, h2, log_Z_dyn = self.info_dynamics_params
        lp += x[0].dot(h_init)
        lp += np.sum(x[:-1] * h1)
        lp += np.sum(x[1:] * h2)

        # prior log p(x) -- normalization constants
        lp += log_Z_init
        lp += np.sum(log_Z_dyn)

        # likelihood log p(y | x)
        lp += self.log_conditional_likelihood(x)

        return lp
1266834
def init_discrete_hmm_ml(C_full, nstates, reversible=True, stationary=True, active_set=None, P=None,
                         eps_A=None, eps_B=None, separate=None):
    """"""
    raise NotImplementedError('ML-initialization not yet implemented')
1336322
def Dlmk(l,m,k,phi1,phi2,theta1,theta2):
    
    
    return exp(complex(0.,-m*phi1)) * dlmk(l,m,k,theta1) * \
      exp(complex(0.,-k*gamma(phi1,phi2,theta1,theta2)))
668800
def estimate_L(model,bounds,storehistory=True):
    
    def df(x,model,x0):
        x = np.atleast_2d(x)
        dmdx,_ = model.predictive_gradients(x)
        res = np.sqrt((dmdx*dmdx).sum(1)) # simply take the norm of the expectation of the gradient
        return -res

    samples = samples_multidimensional_uniform(bounds,500)
    samples = np.vstack([samples,model.X])
    pred_samples = df(samples,model,0)
    x0 = samples[np.argmin(pred_samples)]
    res = scipy.optimize.minimize(df,x0, method='L-BFGS-B',bounds=bounds, args = (model,x0), options = {'maxiter': 200})
    minusL = res.fun[0][0]
    L = -minusL
    if L<1e-7: L=10  ## to avoid problems in cases in which the model is flat.
    return L
167945
def sample_from_likelihood(self, n_timesteps=10):
        

        self.latent_state_sequences = lmap(
            lambda A: ltake(
                n_timesteps,
                iterate(
                    lambda s: pd.Series(A @ s.values, index=s.index), self.s0
                ),
            ),
            self.transition_matrix_collection,
        )

        self.observed_state_sequences = [
            [self.sample_observed_state(s) for s in latent_state_sequence]
            for latent_state_sequence in self.latent_state_sequences
        ]
913300
def log_likelihood(C, T):
    
    C = C.tocsr()
    T = T.tocsr()
    ind = scipy.nonzero(C)
    relT = np.array(T[ind])[0, :]
    relT = np.log(relT)
    relC = np.array(C[ind])[0, :]

    return relT.dot(relC)
851071
def data_log_likelihood(self, successes, trials, beta):
        ''''''
        return binom.logpmf(successes, trials, 1.0 / (1 + np.exp(-beta))).sum()
328890
def infer_dims(self, x, y, dims_x, dims_y, dims_out):
        
        OptimizedInverseModel.infer_x(self, y)
        assert len(x) == len(dims_x)
        assert len(y) == len(dims_y)
        if len(self.fmodel.dataset) == 0:
            return [[0.0]*self.dim_out]
        else:
            _, index = self.fmodel.dataset.nn_dims(x, y, dims_x, dims_y, k=1)
            guesses = [self.fmodel.dataset.get_dims(index[0], dims_out)]

            result = []
            for g in guesses:
                res = cma.fmin(lambda q:self._error_dims(q, dims_x, dims_y, dims_out), g, self.cmaes_sigma,
                               options={'bounds':[self.lower, self.upper],
                               'verb_log':0,
                               'verb_disp':False,
                               'maxfevals':self.maxfevals,
                               'seed': self.seed})
                result.append((res[1], res[0]))

            return sorted(result)[0][1]
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
1002897
def infer_cm(tpm):
    
    network_size = tpm.shape[-1]
    all_contexts = tuple(all_states(network_size - 1))
    cm = np.empty((network_size, network_size), dtype=int)
    for a, b in np.ndindex(cm.shape):
        cm[a][b] = infer_edge(tpm, a, b, all_contexts)
    return cm
419981
def inference_minibatch(self, kern, X, Z, likelihood, Y):
        

        num_data, output_dim = Y.shape

        if isinstance(X, VariationalPosterior):
            uncertain_inputs = True
        else:
            uncertain_inputs = False

        #see whether we've got a different noise variance for each datum
        beta = 1./np.fmax(likelihood.variance, 1e-6)
        het_noise = beta.size > 1
        # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
        #self.YYTfactor = beta*self.get_YYTfactor(Y)
        if self.Y_speedup and not het_noise:
            YYT_factor = self.get_YYTfactor(Y)
        else:
            YYT_factor = Y

        n_start = self.batch_pos
        batchsize = num_data if self.batchsize is None else self.batchsize
        n_end = min(batchsize+n_start, num_data)
        if n_end==num_data:
            isEnd = True
            self.batch_pos = 0
        else:
            isEnd = False
            self.batch_pos = n_end

        if batchsize==num_data:
            Y_slice = YYT_factor
            X_slice =X
        else:
            Y_slice = YYT_factor[n_start:n_end]
            X_slice = X[n_start:n_end]

        if not uncertain_inputs:
            psi0 = kern.Kdiag(X_slice)
            psi1 = kern.K(X_slice, Z)
            psi2 = None
            betapsi1 = np.einsum('n,nm->nm',beta,psi1)
        elif het_noise:
            psi0 = kern.psi0(Z, X_slice)
            psi1 = kern.psi1(Z, X_slice)
            psi2 = kern.psi2(Z, X_slice)
            betapsi1 = np.einsum('n,nm->nm',beta,psi1)

        if het_noise:
            beta = beta[n_start] # assuming batchsize==1

        betaY = beta*Y_slice

        #======================================================================
        # Load Intermediate Results
        #======================================================================

        dL_dpsi2R = self.midRes['dL_dpsi2R']
        v = self.midRes['v']

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

        dL_dpsi0 = -output_dim * (beta * np.ones((n_end-n_start,)))/2.

        dL_dpsi1 = np.dot(betaY,v.T)

        if uncertain_inputs:
            dL_dpsi2 = beta* dL_dpsi2R
        else:
            dL_dpsi1 += np.dot(betapsi1,dL_dpsi2R)*2.
            dL_dpsi2 = None

        #======================================================================
        # Compute dL_dthetaL
        #======================================================================

        if het_noise:
            if uncertain_inputs:
                psiR = np.einsum('mo,mo->',dL_dpsi2R,psi2)
            else:
                psiR = np.einsum('nm,no,mo->',psi1,psi1,dL_dpsi2R)

            dL_dthetaL = ((np.square(betaY)).sum(axis=-1) + np.square(beta)*(output_dim*psi0)-output_dim*beta)/2. - np.square(beta)*psiR- (betaY*np.dot(betapsi1,v)).sum(axis=-1)
        else:
            if isEnd:
                dL_dthetaL = self.midRes['dL_dthetaL']
            else:
                dL_dthetaL = 0.

        if uncertain_inputs:
            grad_dict = {'dL_dpsi0':dL_dpsi0,
                         'dL_dpsi1':dL_dpsi1,
                         'dL_dpsi2':dL_dpsi2,
                         'dL_dthetaL':dL_dthetaL}
        else:
            grad_dict = {'dL_dKdiag':dL_dpsi0,
                         'dL_dKnm':dL_dpsi1,
                         'dL_dthetaL':dL_dthetaL}

        return isEnd, (n_start,n_end), grad_dict