Check convergence based on mean squared error

        Returns
        -------

        converged : boolean
            Whether the parameter estimation converged.

        mse : float
            Mean squared error between prior and posterior.
179612
def _mse_converged(self):
        

        mse = mean_squared_error(self.local_prior, self.local_posterior_,
                                 multioutput='uniform_average')
        if mse > self.threshold:
            return False, mse
        else:
            return True, mse
952066
def _meanprecision(D, tol=1e-7, maxiter=None):
    ''''''
    N, K = D.shape
    logp = log(D).mean(axis=0)
    a0 = _init_a(D)
    s0 = a0.sum()
    if s0 < 0:
        a0 = a0/s0
        s0 = 1
    elif s0 == 0:
        a0 = ones(a.shape) / len(a)
        s0 = 1
    m0 = a0/s0

    # Start updating
    if maxiter is None:
        maxiter = MAXINT
    for i in xrange(maxiter):
        a1 = _fit_s(D, a0, logp, tol=tol)
        s1 = sum(a1)
        a1 = _fit_m(D, a1, logp, tol=tol)
        m = a1/s1
        # if norm(a1-a0) < tol:
        if abs(loglikelihood(D, a1)-loglikelihood(D, a0)) < tol: # much faster
            return a1
        a0 = a1
    raise Exception('Failed to converge after {} iterations, values are {}.'
                    .format(maxiter, a1))
56599
def _is_converged(self):
        ''''''
        if self.is_relaxed():
            # relaxation run case
            return self._get_line(['End of', 'Geometry Optimization'], self.outputf, return_string=False)
        else:
            # static run case
            return self._get_line('convergence has been achieved', self.outputf, return_string=False)
762343
def estimate(self, maxiter=250, convergence=1e-7):

        

        self.loglik = np.zeros(maxiter)

        iter = 0

        while iter < maxiter:

            self.loglik[iter] = self.E_step()
            if np.isnan(self.loglik[iter]):
                    print("undefined log-likelihood")
                    break
            self.M_step()

            if self.loglik[iter] - self.loglik[iter - 1] < 0 and iter > 0:
                    print("log-likelihood decreased by %f at iteration %d"
                          % (self.loglik[iter] - self.loglik[iter - 1],
                             iter))
            elif self.loglik[iter] - self.loglik[iter - 1] < convergence \
                    and iter > 0:
                    print("convergence at iteration %d, loglik = %f" %
                           (iter, self.loglik[iter]))
                    self.loglik = self.loglik[self.loglik < 0]
                    break

            iter += 1
955951
def map_parameter_estimate(self):
        
        res = self.pst.res
        assert res is not None
        # build the prior expectation parameter vector
        prior_expt = self.pst.parameter_data.loc[:,["parval1"]].copy()
        islog = self.pst.parameter_data.partrans == "log"
        prior_expt.loc[islog] = prior_expt.loc[islog].apply(np.log10)
        prior_expt = Matrix.from_dataframe(prior_expt)
        prior_expt.col_names = ["prior_expt"]
        # build the residual vector
        res_vec = Matrix.from_dataframe(res.loc[:,["residual"]])

        # form the terms of Schur's complement
        b = self.parcov * self.jco.T
        c = ((self.jco * self.parcov * self.jco.T) + self.obscov).inv
        bc = Matrix((b * c).x, row_names=b.row_names, col_names=c.col_names)

        # calc posterior expectation
        upgrade = bc * res_vec
        upgrade.col_names = ["prior_expt"]
        post_expt = prior_expt + upgrade

        # post processing - back log transform
        post_expt = pd.DataFrame(data=post_expt.x,index=post_expt.row_names,
                                 columns=["post_expt"])
        post_expt.loc[:,"prior_expt"] = prior_expt.x.flatten()
        post_expt.loc[islog,:] = 10.0**post_expt.loc[islog,:]
        # unc_sum = self.get_parameter_summary()
        # post_expt.loc[:,"standard_deviation"] = unc_sum.post_var.apply(np.sqrt)
        post_expt.sort_index(inplace=True)
        return post_expt
843724
def predict(self, X):
        
        X = np.atleast_2d(X)
        m = np.empty(shape=(0,1))
        s = np.empty(shape=(0,1))

        for k in range(X.shape[0]):
            preds = []
            for pred in self.model.estimators_:
                preds.append(pred.predict(X[k,:])[0])
            m = np.vstack((m ,np.array(preds).mean()))
            s = np.vstack((s ,np.array(preds).std()))
        return m, s
800198
def check_convergence(newX, oldX, e_rel):
    
    # Calculate the norm for columns and rows, which can be used for debugging
    # Otherwise skip, since it takes extra processing time
    new_old = newX*oldX
    old2 = oldX**2
    norms = [np.sum(new_old), np.sum(old2)]
    convergent = norms[0] >= (1-e_rel**2)*norms[1]
    return convergent, norms
288156
public static boolean isConvergenced(List<KmeansPoint> basePoints, List<KmeansPoint> newPoints,
            double convergenceThres)
    {
        boolean result = true;

        for (int pointIndex = 0; pointIndex < basePoints.size(); pointIndex++)
        {
            double distance = MathUtils.distance(basePoints.get(pointIndex).getDataPoint(),
                    newPoints.get(pointIndex).getDataPoint());
            if (distance > convergenceThres)
            {
                result = false;
                break;
            }
        }

        return result;
    }
85704
def _compute_mle(self):
        
        xmax = self._lnlfn.interp.xmax
        x0 = max(self._lnlfn.mle(), xmax * 1e-5)
        ret = opt.fmin(lambda x: np.where(
            xmax > x > 0, -self(x), np.inf), x0, disp=False)
        mle = float(ret[0])
        return mle
326373
def posterior(self, x, s=1.):
        """"""
        pr0 = 1. / self.sigma**2  # prior precision
        prd = x.size / s**2  # data precision
        varp = 1. / (pr0 + prd)  # posterior variance
        mu = varp * (pr0 * self.mu + prd * x.mean())
        return TruncNormal(mu=mu, sigma=np.sqrt(varp), a=self.a, b=self.b)
1144288
def _mean_absolute_error(y, y_pred, w):
    """"""
    return np.average(np.abs(y_pred - y), weights=w)
68527
def estimate_gaussian(X):
    
    mean = np.mean(X,0)
    variance = np.var(X,0)
    return Gaussian(mean,variance)
919926
def check_constraint_convergence(X, L, LX, Z, U, R, S, step_f, step_g, e_rel, e_abs):
    

    if isinstance(L, list):
        M = len(L)
        convergence = True
        errors = []
        # recursive call
        for i in range(M):
            c, e = check_constraint_convergence(X, L[i], LX[i], Z[i], U[i], R[i], S[i],
                                                step_f, step_g[i], e_rel, e_abs)
            convergence &= c
            errors.append(e)
        return convergence, errors
    else:
        # check convergence of prime residual R and dual residual S
        e_pri, e_dual = get_variable_errors(X, L, LX, Z, U, step_g, e_rel, e_abs)
        lR = l2(R)
        lS = l2(S)
        convergence = (lR <= e_pri) and (lS <= e_dual)
        return convergence, (e_pri, e_dual, lR, lS)
402570
def mean_squared_error(pred:Tensor, targ:Tensor)->Rank0Tensor:
    ""
    pred,targ = flatten_check(pred,targ)
    return F.mse_loss(pred, targ)
856522
def accuracy(batch, model_predictions):
  """"""
  _, targets = batch
  model_predictions, targets = _make_list(model_predictions, targets)
  correct = []
  for (prediction, target) in zip(model_predictions, targets):
    predicted_class = np.argmax(prediction, axis=-1)
    correct.append(np.equal(predicted_class, target))
  return masked_mean(correct, targets)
1008452
def _check_params(self):
        
        if self.n_estimators <= 0:
            raise ValueError("n_estimators must be greater than 0 but "
                             "was %r" % self.n_estimators)

        if not 0.0 < self.subsample <= 1.0:
            raise ValueError("subsample must be in ]0; 1] but "
                             "was %r" % self.subsample)

        if not 0.0 < self.learning_rate <= 1.0:
            raise ValueError("learning_rate must be within ]0; 1] but "
                             "was %r" % self.learning_rate)

        if not 0.0 <= self.dropout_rate < 1.0:
            raise ValueError("dropout_rate must be within [0; 1[, but "
                             "was %r" % self.dropout_rate)

        if self.loss not in LOSS_FUNCTIONS:
            raise ValueError("Loss '{0:s}' not supported. ".format(self.loss))
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