multi_nested_integrator.py 文件源码

def integrate_remainder(sampler, logwidth, logVolremaining, logZ, H, globalLmax):
    # logwidth remains the same now for each sample
    remainder = list(sampler.remainder())
    logV = logwidth
    L0 = remainder[-1][2]
    L0 = globalLmax
    logLs = [Li - L0 for ui, xi, Li in remainder]
    Ls = numpy.exp(logLs)
    LsMax = Ls.copy()
    LsMax[-1] = numpy.exp(globalLmax - L0)
    Lmax = LsMax[1:].sum(axis=0) + LsMax[-1]
    #Lmax = Ls[1:].sum(axis=0) + Ls[-1]
    Lmin = Ls[:-1].sum(axis=0) + Ls[0]
    logLmid = log(Ls.sum(axis=0)) + L0
    logZmid = logaddexp(logZ, logV + logLmid)
    logZup  = logaddexp(logZ, logV + log(Lmax) + L0)
    logZlo  = logaddexp(logZ, logV + log(Lmin) + L0)
    logZerr = logZup - logZlo
    assert numpy.isfinite(H).all()
    assert numpy.isfinite(logZerr).all(), logZerr

    for i in range(len(remainder)):
        ui, xi, Li = remainder[i]
        wi = logwidth + Li
        logZnew = logaddexp(logZ, wi)
        #Hprev = H
        H = exp(wi - logZnew) * Li + exp(logZ - logZnew) * (H + logZ) - logZnew
        H[H < 0] = 0
        #assert (H>0).all(), (H, Hprev, wi, Li, logZ, logZnew)
        logZ = logZnew

    #assert numpy.isfinite(logZerr + (H / sampler.nlive_points)**0.5), (H, sampler.nlive_points, logZerr)

    return logV + logLmid, logZerr, logZmid, logZerr + (H / sampler.nlive_points)**0.5, logZerr + (H / sampler.nlive_points)**0.5