class_docstr: Simple Nested sampler for reference.
class_docstr_tok: ["Simple","Nested","sampler","for","reference","."]
docstr: Explore parameter space.
  
  Parameters
  ----------
  update_interval_iter:
      Update region after this many iterations.
  update_interval_ncall:
      Update region after update_interval_ncall likelihood calls.
  log_interval:
      Update stdout status line every log_interval iterations
  dlogz:
      Target evidence uncertainty.
  max_iters:
      maximum number of integration iterations.
docstr_tok: ["Explore","parameter","space",".","Parameters","--","--","--","--","--","update_interval_iter",":","Update","region","after","this","many","iterations",".","update_interval_ncall",":","Update","region","after","update_interval_ncall","likelihood","calls",".","log_interval",":","Update","stdout","status","line","every","log_interval","iterations","dlogz",":","Target","evidence","uncertainty",".","max_iters",":","maximum","number","of","integration","iterations","."]
code: def run(
      self,
      update_interval_iter=None,
      update_interval_ncall=None,
      log_interval=None,
      dlogz=0.001,
      max_iters=None,
  ):
      """Explore parameter space.
  
      Parameters
      ----------
      update_interval_iter:
          Update region after this many iterations.
      update_interval_ncall:
          Update region after update_interval_ncall likelihood calls.
      log_interval:
          Update stdout status line every log_interval iterations
      dlogz:
          Target evidence uncertainty.
      max_iters:
          maximum number of integration iterations.
  
      """
      if update_interval_ncall is None:
          update_interval_ncall = max(1, round(self.num_live_points))
  
      if update_interval_iter is None:
          if update_interval_ncall == 0:
              update_interval_iter = max(1, round(self.num_live_points))
          else:
              update_interval_iter = max(
                  1, round(0.2 * self.num_live_points)
              )
  
      if log_interval is None:
          log_interval = max(1, round(0.2 * self.num_live_points))
      else:
          log_interval = round(log_interval)
          if log_interval < 1:
              raise ValueError("log_interval must be >= 1")
  
      viz_callback = get_default_viz_callback()
  
      prev_u = []
      prev_v = []
      prev_logl = []
      if self.log:
          # try to resume:
          self.logger.info("Resuming...")
          for i in range(self.num_live_points):
              _, row = self.pointstore.pop(-np.inf)
              if row is not None:
                  prev_logl.append(row[1])
                  prev_u.append(row[3 : 3 + self.x_dim])
                  prev_v.append(
                      row[
                          3
                          + self.x_dim : 3
                          + self.x_dim
                          + self.num_params
                      ]
                  )
              else:
                  break
  
          prev_u = np.array(prev_u)
          prev_v = np.array(prev_v)
          prev_logl = np.array(prev_logl)
  
          num_live_points_missing = self.num_live_points - len(
              prev_logl
          )
      else:
          num_live_points_missing = -1
  
      if self.use_mpi:
          num_live_points_missing = self.comm.bcast(
              num_live_points_missing, root=0
          )
          prev_u = self.comm.bcast(prev_u, root=0)
          prev_v = self.comm.bcast(prev_v, root=0)
          prev_logl = self.comm.bcast(prev_logl, root=0)
  
      use_point_stack = True
  
      assert num_live_points_missing >= 0
      if num_live_points_missing > 0:
          if self.use_mpi:
              # self.logger.info('Using MPI with rank [%d]', self.mpi_rank)
              if self.mpi_rank == 0:
                  active_u = np.random.uniform(
                      size=(num_live_points_missing, self.x_dim)
                  )
              else:
                  active_u = np.empty(
                      (num_live_points_missing, self.x_dim),
                      dtype=np.float64,
                  )
              active_u = self.comm.bcast(active_u, root=0)
          else:
              active_u = np.random.uniform(
                  size=(num_live_points_missing, self.x_dim)
              )
          active_v = self.transform(active_u)
  
          if self.use_mpi:
              if self.mpi_rank == 0:
                  chunks = [[] for _ in range(self.mpi_size)]
                  for i, chunk in enumerate(active_v):
                      chunks[i % self.mpi_size].append(chunk)
              else:
                  chunks = None
              data = self.comm.scatter(chunks, root=0)
              active_logl = self.loglike(data)
              recv_active_logl = self.comm.gather(active_logl, root=0)
              recv_active_logl = self.comm.bcast(
                  recv_active_logl, root=0
              )
              active_logl = np.concatenate(recv_active_logl, axis=0)
          else:
              active_logl = self.loglike(active_v)
  
          if self.log_to_disk:
              for i in range(num_live_points_missing):
                  self.pointstore.add(
                      _listify(
                          [-np.inf, active_logl[i], 0.0],
                          active_u[i, :],
                          active_v[i, :],
                      ),
                      num_live_points_missing,
                  )
  
          if len(prev_u) > 0:
              active_u = np.concatenate((prev_u, active_u))
              active_v = np.concatenate((prev_v, active_v))
              active_logl = np.concatenate((prev_logl, active_logl))
          assert active_u.shape == (self.num_live_points, self.x_dim)
          assert active_v.shape == (
              self.num_live_points,
              self.num_params,
          )
          assert active_logl.shape == (self.num_live_points,)
      else:
          active_u = prev_u
          active_v = prev_v
          active_logl = prev_logl
  
      saved_u = []
      saved_v = []  # Stored points for posterior results
      saved_logl = []
      saved_logwt = []
      h = 0.0  # Information, initially 0.
      logz = -1e300  # ln(Evidence Z), initially Z=0
      logvol = log(1.0 - exp(-1.0 / self.num_live_points))
      logz_remain = np.max(active_logl)
      fraction_remain = 1.0
      ncall = num_live_points_missing  # number of calls we already made
      first_time = True
      if self.x_dim > 1:
          transformLayer = AffineLayer(wrapped_dims=self.wrapped_axes)
      else:
          transformLayer = ScalingLayer(wrapped_dims=self.wrapped_axes)
      transformLayer.optimize(active_u, active_u)
      region = MLFriends(active_u, transformLayer)
  
      if self.log:
          self.logger.info("Starting sampling ...")
      ib = 0
      samples = []
      ndraw = 100
      it = 0
      next_update_interval_ncall = -1
      next_update_interval_iter = -1
  
      while max_iters is None or it < max_iters:
          # Worst object in collection and its weight (= volume * likelihood)
          worst = np.argmin(active_logl)
          logwt = logvol + active_logl[worst]
  
          # Update evidence Z and information h.
          logz_new = np.logaddexp(logz, logwt)
          h = (
              exp(logwt - logz_new) * active_logl[worst]
              + exp(logz - logz_new) * (h + logz)
              - logz_new
          )
          logz = logz_new
  
          # Add worst object to samples.
          saved_u.append(np.array(active_u[worst]))
          saved_v.append(np.array(active_v[worst]))
          saved_logwt.append(logwt)
          saved_logl.append(active_logl[worst])
  
          # expected_vol = np.exp(-it / self.num_live_points)
  
          # The new likelihood constraint is that of the worst object.
          loglstar = active_logl[worst]
  
          if (
              ncall > next_update_interval_ncall
              and it > next_update_interval_iter
          ):
              if first_time:
                  nextregion = region
              else:
                  # rebuild space
                  # print()
                  # print("rebuilding space...", active_u.shape, active_u)
                  nextTransformLayer = transformLayer.create_new(
                      active_u, region.maxradiussq
                  )
                  nextregion = MLFriends(active_u, nextTransformLayer)
  
              # print("computing maxradius...")
              r, f = _update_region_bootstrap(
                  nextregion,
                  30,
                  0.0,
                  self.comm if self.use_mpi else None,
                  self.mpi_size,
              )
  
              nextregion.maxradiussq = r
              nextregion.enlarge = f
              # force shrinkage of volume
              # this is to avoid re-connection of dying out nodes
              if (
                  nextregion.estimate_volume()
                  < region.estimate_volume()
              ):
                  region = nextregion
                  transformLayer = region.transformLayer
              region.create_ellipsoid(
                  minvol=exp(-it / self.num_live_points)
                  * self.volfactor
              )
  
              if self.log:
                  viz_callback(
                      points=dict(
                          u=active_u, p=active_v, logl=active_logl
                      ),
                      info=dict(
                          it=it,
                          ncall=ncall,
                          logz=logz,
                          logz_remain=logz_remain,
                          paramnames=self.paramnames
                          + self.derivedparamnames,
                          logvol=logvol,
                      ),
                      region=region,
                      transformLayer=transformLayer,
                  )
                  self.pointstore.flush()
  
              next_update_interval_ncall = ncall + update_interval_ncall
              next_update_interval_iter = it + update_interval_iter
              first_time = False
  
          while True:
              if ib >= len(samples) and use_point_stack:
                  # root checks the point store
                  next_point = np.zeros(
                      (1, 3 + self.x_dim + self.num_params)
                  )
  
                  if self.log_to_disk:
                      _, stored_point = self.pointstore.pop(loglstar)
                      if stored_point is not None:
                          next_point[0, :] = stored_point
                      else:
                          next_point[0, :] = -np.inf
                      use_point_stack = not self.pointstore.stack_empty
  
                  if self.use_mpi:  # and informs everyone
                      use_point_stack = self.comm.bcast(
                          use_point_stack, root=0
                      )
                      next_point = self.comm.bcast(next_point, root=0)
  
                  # assert not use_point_stack
  
                  # unpack
                  likes = next_point[:, 1]
                  samples = next_point[:, 3 : 3 + self.x_dim]
                  samplesv = next_point[
                      :,
                      3 + self.x_dim : 3 + self.x_dim + self.num_params,
                  ]
                  # skip if we already know it is not useful
                  ib = 0 if np.isfinite(likes[0]) else 1
  
              while ib >= len(samples):
                  # get new samples
                  ib = 0
  
                  nc = 0
                  u = region.sample(nsamples=ndraw)
                  nu = u.shape[0]
                  if nu == 0:
                      v = np.empty((0, self.x_dim))
                      logl = np.empty((0,))
                  else:
                      v = self.transform(u)
                      logl = self.loglike(v)
                      nc += nu
                      accepted = logl > loglstar
                      u = u[accepted, :]
                      v = v[accepted, :]
                      logl = logl[accepted]
                      # father = father[accepted]
  
                  # collect results from all MPI members
                  if self.use_mpi:
                      recv_samples = self.comm.gather(u, root=0)
                      recv_samplesv = self.comm.gather(v, root=0)
                      recv_likes = self.comm.gather(logl, root=0)
                      recv_nc = self.comm.gather(nc, root=0)
                      recv_samples = self.comm.bcast(
                          recv_samples, root=0
                      )
                      recv_samplesv = self.comm.bcast(
                          recv_samplesv, root=0
                      )
                      recv_likes = self.comm.bcast(recv_likes, root=0)
                      recv_nc = self.comm.bcast(recv_nc, root=0)
                      samples = np.concatenate(recv_samples, axis=0)
                      samplesv = np.concatenate(recv_samplesv, axis=0)
                      likes = np.concatenate(recv_likes, axis=0)
                      ncall += sum(recv_nc)
                  else:
                      samples = np.array(u)
                      samplesv = np.array(v)
                      likes = np.array(logl)
                      ncall += nc
  
                  if self.log:
                      for ui, vi, logli in zip(
                          samples, samplesv, likes
                      ):
                          self.pointstore.add(
                              _listify([loglstar, logli, 0.0], ui, vi),
                              ncall,
                          )
  
              if likes[ib] > loglstar:
                  active_u[worst] = samples[ib, :]
                  active_v[worst] = samplesv[ib, :]
                  active_logl[worst] = likes[ib]
  
                  # if we keep the region informed about the new live points
                  # then the region follows the live points even if maxradius is not updated
                  region.u[worst, :] = active_u[worst]
                  region.unormed[
                      worst, :
                  ] = region.transformLayer.transform(
                      region.u[worst, :]
                  )
  
                  # if we track the cluster assignment, then in the next round
                  # the ids with the same members are likely to have the same id
                  # this is imperfect
                  # transformLayer.clusterids[worst] = transformLayer.clusterids[father[ib]]
                  # so we just mark the replaced ones as "unassigned"
                  transformLayer.clusterids[worst] = 0
                  ib = ib + 1
                  break
              else:
                  ib = ib + 1
  
          # Shrink interval
          logvol -= 1.0 / self.num_live_points
          logz_remain = np.max(active_logl) - it / self.num_live_points
          fraction_remain = np.logaddexp(logz, logz_remain) - logz
  
          if it % log_interval == 0 and self.log:
              # nicelogger(self.paramnames, active_u, active_v, active_logl, it, ncall, logz, logz_remain, region=region)
              sys.stdout.write(
                  "Z=%.1g+%.1g | Like=%.1g..%.1g | it/evals=%d/%d eff=%.4f%%  \r"
                  % (
                      logz,
                      logz_remain,
                      loglstar,
                      np.max(active_logl),
                      it,
                      ncall,
                      np.inf if ncall == 0 else it * 100 / ncall,
                  )
              )
              sys.stdout.flush()
  
              # if efficiency becomes low, bulk-process larger arrays
              ndraw = max(
                  128,
                  min(
                      16384,
                      round((ncall + 1) / (it + 1) / self.mpi_size),
                  ),
              )
  
          # Stopping criterion
          if fraction_remain < dlogz:
              break
          it = it + 1
  
      logvol = -len(saved_v) / self.num_live_points - log(
          self.num_live_points
      )
      for i in range(self.num_live_points):
          logwt = logvol + active_logl[i]
          logz_new = np.logaddexp(logz, logwt)
          h = (
              exp(logwt - logz_new) * active_logl[i]
              + exp(logz - logz_new) * (h + logz)
              - logz_new
          )
          logz = logz_new
          saved_u.append(np.array(active_u[i]))
          saved_v.append(np.array(active_v[i]))
          saved_logwt.append(logwt)
          saved_logl.append(active_logl[i])
  
      saved_u = np.array(saved_u)
      saved_v = np.array(saved_v)
      saved_wt = exp(np.array(saved_logwt) - logz)
      saved_logl = np.array(saved_logl)
      logzerr = np.sqrt(h / self.num_live_points)
  
      if self.log_to_disk:
          with open(
              os.path.join(self.logs["results"], "final.csv"), "w"
          ) as f:
              writer = csv.writer(f)
              writer.writerow(
                  ["niter", "ncall", "logz", "logzerr", "h"]
              )
              writer.writerow([it + 1, ncall, logz, logzerr, h])
          self.pointstore.close()
  
      if not self.use_mpi or self.mpi_rank == 0:
          print()
          print(
              "niter: {:d}\n ncall: {:d}\n nsamples: {:d}\n logz: {:6.3f} +/- {:6.3f}\n h: {:6.3f}".format(
                  it + 1, ncall, len(saved_v), logz, logzerr, h
              )
          )
  
      self.results = dict(
          samples=resample_equal(saved_v, saved_wt / saved_wt.sum()),
          ncall=ncall,
          niter=it,
          logz=logz,
          logzerr=logzerr,
          weighted_samples=dict(
              upoints=saved_u,
              points=saved_v,
              weights=saved_wt,
              logweights=saved_logwt,
              logl=saved_logl,
          ),
      )
  
      return self.results
  
imports: from __future__ import print_function, division
  import os
  import sys
  import csv
  import json
  import operator
  import time
  import warnings
  from numpy import log, exp, logaddexp
  import numpy
  from .utils import create_logger, make_run_dir, resample_equal, vol_prefactor, vectorize, listify
  from .utils import is_affine_transform, normalised_kendall_tau_distance, distributed_work_chunk_size
  from ultranest.mlfriends import MLFriends, AffineLayer, ScalingLayer, find_nearby, WrappingEllipsoid, RobustEllipsoidRegion
  from .store import HDF5PointStore, TextPointStore, NullPointStore
  from .viz import get_default_viz_callback
  from .ordertest import UniformOrderAccumulator
  from .netiter import PointPile, SingleCounter, MultiCounter, BreadthFirstIterator, TreeNode, count_tree_between, find_nodes_before, logz_sequence
  from .netiter import dump_tree, combine_results
  from .hotstart import get_auxiliary_contbox_parameterization
label_desc: Use image node_id 2 for calling the NestedSampler obj's underlying member method code with example usage: obj.run(update_interval_iter, update_interval_ncall, log_interval, dlogz, max_iters) and returns: self