pmf values = 0 cause cut-out regions in plots

[Stefan-Heimersheim]

Describe the bug
I have a sample from my prior (all weights = 1) of a interpolation function (FlexKnot, free-form parameterization etc.) and plot them using fgivenx. (Lists, as I have multiple, evidence-weighted, samples.)

cbar = fgivenx.plot_contours(function_list, zplot, chain_list, ax=ax, weights=weight_list, logZ=logZ_list, cache=cachestr, contour_line_levels=[1,2], fineness=0.2, smooth=0)

and I get plots like these (different sample sizes, 3rd different fineness/levels):

Did you see something like this before? Is it just too small sample size or is there some issue computing the pmfs?

The plots break due to infinite values (erfinv(1)) of z=0(pmf=0) in plot.py, plot(). Setting inf to 1e10 allows the lines to be plotted but does not fix the issue of why the pmf is computed as 0.

# Convert to sigmas
infs = numpy.where(z==0)
z = numpy.sqrt(2) * scipy.special.erfinv(1 - z)
z[infs] = 1e10

Adding e.g. smooth=2 does not improve the situation (here applied to the 2nd plot, with the z[infs] = 1e10 "fix", it makes the cut-out larger as it seems to propagate the problem over the smoothing area):

Edit: I just noticed the smoothing code, of course it cannot work as the smoothing is applied to the contours (where I have put an 1e10) and not the PDF

I'm not sure if the problem is due to my sample size (1st plot 1e5 samples, 2nd plot 1e4) or some issue, plotting for 1e5 samples however already takes an hour so it is not really feasible to increase the sample size much.

To Reproduce
If you think this is an important bug I can work out an MWE; right now I was mainly wondering if this is a common problem you see often or something that shouldn't happen.

Expected behavior
a) Smooth contours (physics dependent)
b) Continuous contour lines (by definition)

[williamjameshandley]

Hi @Stefan-Heimersheim, this isn't behaviour I've seen before, but if you can make a MWE I'll have a go at debugging

[Stefan-Heimersheim]

Okay, I got an MWE below and I played around a bit -- the problem seems to be dependent on the x-axis samples. With a straightforward

zplot = np.linspace(5,30,10)

Only the first and last bin are empty:

But if I try to get around this by adding smaller bins at the start and end I get weird empty patches

zplot = np.concatenate(([5.1], np.linspace(5.2,29.8,10), [29.9]))

Here's the full example code, running in ~5s, without any files needed,

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as sip

from fgivenx import plot_contours
from anesthetic import MCMCSamples

Nprior = 1000

# Generate prior samples to be able to correct for the priors
# Generate ordered points as in Handley et al. 2015 A2 (doi:10.1093/mnras/stv1911)
def generatePriorFlexKnotSamples(size=1000):
    # Note: Not the original FlexKnot, leaving start and end point to vary in z
    tmp = {}
    tmp['z0'] = 5.001
    tmp['x1'] = 1
    for i in range(3):
        j = str(i+1)
        tmp['v'+j] = np.random.uniform(low=0, high=1, size=size)
        tmp['z'+j] = eval("tmp['z"+str(i)+"']+(30-tmp['z"+str(i)+"'])*(1-tmp['v"+j+"']**(1/("+str(2+1)+"-"+j+"+1)))")
        tmp.pop('v'+j)
    for i in range(1):
        j = str(i+1)
        tmp['u'+j] = np.random.uniform(low=0, high=1, size=size)
        tmp['x'+str(i+2)] = eval("tmp['x"+str(i+1)+"']*tmp['u"+j+"']**(1/("+str(2-1)+"-"+j+"+1))")
        tmp.pop('u'+j)
    tmp.pop('x1')
    tmp.pop('z0')
    return MCMCSamples(tmp)

chain = MCMCSamples(generatePriorFlexKnotSamples(size=Nprior))
chain.weights = np.ones(len(chain))


def xifunc(z, p, keys=chain.keys(), usePCHIP=False):
    # Interpolate FlexKnot
    xfull = [1,1]
    zfull = [0, *p[np.where(keys=='z1')]]
    zfull.append(*p[np.where(keys=="z2")])
    xfull.append(*p[np.where(keys=="x2")])
    zfull.append(*p[np.where(keys=="z3")])
    xfull.append(0)
    zfull.append(50)
    xfull.append(0)
    zfull = np.array(zfull)
    xfull = np.array(xfull)
    assert np.all(np.diff(xfull)<=0), xfull
    assert np.all(np.diff(zfull)>=0), zfull
    if usePCHIP:
        interpolation = sip.pchip(zfull, xfull)
    else:
        interpolation = sip.interp1d(zfull, xfull)
    return interpolation(z)


fig, ax_fgivenx = plt.subplots(figsize=(12, 4))

#zplot = np.linspace(5,30,10)
zplot = np.concatenate(([5.1], np.linspace(5.2,29.8,10), [29.9]))
#zplot = np.concatenate(([5, 5.01], np.linspace(5.02,29.98,10), [29.99, 30]))
#zplot = np.concatenate(([5.1], np.linspace(5.2,29.8,100), [29.9]))

plot_contours(xifunc, zplot, chain,
    weights=chain.weights, ny=100, contour_line_levels=[1,2], rasterize_contours=True,
    fineness=0.25, ax=ax_fgivenx)

#plt.savefig("fgivenx_plot.png", dpi=600)
plt.show()

The plots in the previous posts above use up to 10^6 prior samples, and I also tried the MWE with those -- the empty spaces stayed even with massively increased sample size: