alpiges/LinConGauss

test_subset fails in the wrong case

doweichert opened this issue · 2 comments

doweichert commented

Hi,

to understand the behaviour of the subset simulation, I have slightly modified the code of test_subset.py to generate some plot for an 1D example.

import numpy as np
import matplotlib.pyplot as plt

from LinConGauss import LinearConstraints
from LinConGauss.multilevel_splitting import SubsetSimulation

def plot_nestings():
    """ Plot nesting history with 1D constraints. """
    x = np.linspace(-2, 2, 100).reshape(1, -1)
    constraint_tracker = subset_simulator.tracker.nestings
    num_nestings = len(constraint_tracker)
    if num_nestings > 1:
        fig, axs = plt.subplots(1, num_nestings, sharex = 'all', sharey = 'all')
        for j, nesting in enumerate(constraint_tracker):

            y = nesting.shifted_lincon.evaluate(x)
            axs[j].plot(x.T, y.T)

            y = nesting.lincon.evaluate(x)
            axs[j].plot(x.T, y.T, color = 'k')
            for c in range(nesting.lincon.N_constraints):
                if nesting.lincon.A[c] > 0:
                    axs[j].fill_between(x.flatten(), 0, y[c], where=(y[c] > 0), color='k', alpha=0.2)
                else:
                    axs[j].fill_between(x.flatten(), 0, y[c], where=(y[c] > 0), color='k', alpha=0.2)

            axs[j].scatter(nesting.x_in, np.zeros_like(nesting.x_in))
            axs[j].hlines(0, -2, 2, color='k')
        plt.title(str(nesting.lincon.integration_domain(subset_simulator.tracker.x_inits()[:,-1]) == 1.))
        plt.savefig(fname=str(seed)+'subset1D.png', format='png')
        plt.close()

if __name__ == '__main__':
    # define some linear constraints
    n_lc = 2
    n_dim = 1

    for seed in range(0, 50):
        np.random.seed(seed)
        lincon = LinearConstraints(2 * np.random.randn(n_lc, n_dim), np.random.randn(n_lc, 1))

        subset_simulator = SubsetSimulation(lincon, 16, 0.5)
        try:
            subset_simulator.run(verbose=False)
            plot_nestings()
        except: print('failed to run')

    print('finished')

Now I'm struggling with the interpretation of its' output for seeds like 7, 10, 13, please find the one for seed 7 attached.
7subset1D

What I don't get is, why the initilization of the last nesting with gamma = 0 (the blue point on the rightmost plot) is not fullfilling both conditions, as the title indicates. Obviously, this does not make sense, as both constraints are larger 0.

Do you have any suggestions? What have I missed?

Best,
Dorina

alpiges commented

Hi Dorina,

It seems to me that this is a nasty python broadcasting error that the code doesn't check for.
The function lincon.integration_domain takes values x with shape (D, N) and returns a 1 if x is in the domain, else 0. That is an array of length 1 no matter how many constraints.
In your case, the shape of x is (D,), which causes lincon.evaluate to broadcast to a (2, 2) array and erroneously suggests that lincon.integration_domainchecks for every constraint individually if it is fulfilled, which is not the case.

So you should write
plt.title(str(nesting.lincon.integration_domain(subset_simulator.tracker.x_inits()[:,-1]).reshape(1, -1) == 1.))
instead of
plt.title(str(nesting.lincon.integration_domain(subset_simulator.tracker.x_inits()[:,-1]) == 1.))

On a vaguely related note, your point reminded me that test_subset.py is incomplete and should be fixed. Furthermore it might be good to check for the shape issue that you encountered in the routines. I'll add this as issues and address them soon.

Let me know if that resolves your issue.

doweichert commented