<jupyter_start><jupyter_text>On two ways to use determinantal point processes for Monte Carlo integration See the [documentation](https://dppy.readthedocs.io/en/latest/continuous_dpps/multivariate_jacobi_ope.html) on ReadTheDocs or the [paper](https://negative-dependence-in-ml-workshop.lids.mit.edu/wp-content/uploads/sites/29/2019/06/icml_camera_ready.pdf) at ICML'19 workshop on Negative Dependence in ML You can play with the different parameters :) Table of Contents1&nbsp;&nbsp;Sampling1.1&nbsp;&nbsp;Minimal working example1.2&nbsp;&nbsp;Plot a sample in 1D or 2D1.3&nbsp;&nbsp;Timing2&nbsp;&nbsp;Numerical integration2.1&nbsp;&nbsp;Estimators2.2&nbsp;&nbsp;Integrands2.3&nbsp;&nbsp;Estimation2.4&nbsp;&nbsp;Variance decay2.4.1&nbsp;&nbsp;Of an integrand $f$2.4.2&nbsp;&nbsp;$f(x) = \sum_{k=0}^{M-1} \frac{1}{k+1} P_k(x)$  Imports If you have not installed `DPPy` yet, you can install it with the following<jupyter_code># !pip install dppy<jupyter_output><empty_output><jupyter_text>💣 **Note: to make sure you have the lastest version of the package please uncomment and run the following cell**<jupyter_code># !rm -r DPPy
# !git clone https://github.com/guilgautier/DPPy.git
# !pip install scipy --upgrade
# !pip install DPPy/.<jupyter_output><empty_output><jupyter_text>💣 If you have chosen to clone the repo and now wish to interact with the source code while running this notebook.You can uncomment the following cell.<jupyter_code># %load_ext autoreload
# %autoreload 2
# import os
# import sys
# sys.path.insert(0, os.path.abspath('..'))
import numpy as np
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'

from scipy import stats
from scipy.integrate import quad

import multiprocessing as mp

from dppy.multivariate_jacobi_ope import MultivariateJacobiOPE<jupyter_output><empty_output><jupyter_text>Sampling  Minimal working example<jupyter_code>d, N = 2, 30  # dimension / number of points
jac_params = 0.5 - np.random.rand(d, 2)  # Jacobi ensemble parameters

dpp = MultivariateJacobiOPE(N, jac_params)
dpp.sample()<jupyter_output><empty_output><jupyter_text>Plot a sample in 1D or 2D<jupyter_code>N, d = 300, 1
    
jac_params = 0.5 - np.random.rand(d, 2)

dpp = MultivariateJacobiOPE(N, jac_params)

sampl = dpp.sample()

print('\n'.join(['Display {} points in {}D'.format(dpp.N, dpp.dim),
                 'with parameters i.i.d. uniformly on [-1/2, 1/2]^{}'.format(dpp.dim),
                 '{}'.format(jac_params)]))

dpp.plot(weighted=False)
dpp.plot(weighted=True)
N, d = 300, 2
    
jac_params = 0.5 - np.random.rand(d, 2)

dpp = MultivariateJacobiOPE(N, jac_params)

sampl = dpp.sample()

print('\n'.join(['Display {} points in {}D'.format(dpp.N, dpp.dim),
                 'with parameters i.i.d. uniformly on [-1/2, 1/2]^{}'.format(dpp.dim),
                 '{}'.format(jac_params)]))

dpp.plot(weighted=False)
dpp.plot(weighted=True)<jupyter_output><empty_output><jupyter_text>Timing To get a quick idea of the time to get a sample you can run the following cell<jupyter_code>d, N = 2, 100
jac_params = 0.5 - np.random.rand(d, 2)
# jac_params = -0.5 * np.ones((d, 2))

dpp = MultivariateJacobiOPE(N, jac_params)

%timeit dpp.sample()<jupyter_output><empty_output><jupyter_text>Numerical integration  Estimators<jupyter_code>def BH_estimator(integrand, dpp, sample=None):
    
    if sample is not None:
        return np.sum(integrand(sample).ravel() / dpp.K(sample))
    
    else:
        sample = dpp.sample()
        return BH_estimator(integrand, dpp, sample)

def EZ_estimator(integrand, dpp, sample=None):

    if sample is not None:
        
        phi_x = dpp.eval_poly_multiD(sample, normalize='norm')
        integrand_x = integrand(sample).ravel()

        EZ_estimator = np.linalg.solve(phi_x, integrand_x)[0]
        EZ_estimator *= np.sqrt(dpp.mass_of_mu)
        
        return EZ_estimator#, np.linalg.cond(phi_sample)

    else:
        
        sample = dpp.sample()
        return EZ_estimator(integrand, dpp, sample)

def both_estimators(integrand, dpp, sample=None):
    
    if sample is not None:
        return BH_estimator(integrand, dpp, sample), EZ_estimator(integrand, dpp, sample)
    
    else:
        np.random.seed(None)
        sample = dpp.sample()
        return both_estimators(integrand, dpp, sample)<jupyter_output><empty_output><jupyter_text>Integrands<jupyter_code>def bump_eps(X, eps=0.05):
    """ https://en.wikipedia.org/wiki/Bump_function
    """

    if type(X) is float:
        f_x = np.exp(-1 / (1.0 - eps - X**2)) if abs(X) < 1 - eps else 0.
        
    elif X.ndim == 1:
        in_I = np.abs(X) < 1 - eps
        f_x = np.zeros(in_I.size)
        f_x[in_I] = np.exp(-1 / (1.0 - eps - X[in_I]**2))

    else:
        in_I = np.all(np.abs(X) < 1 - eps, axis=1)
        f_x = np.zeros(in_I.size)
        f_x[in_I] = np.exp(-np.sum(1.0 / (1.0 - eps - X[in_I]**2), axis=1))
    
    return f_x

def sine(X):

    if type(X) is float: 
        f_x = np.sin(np.pi*X)
    elif X.ndim == 1:
        f_x = np.sin(np.pi*X)
    else:
        f_x = np.prod(np.sin(np.pi*X), axis=-1)
    
    return f_x

def cosine(X):

    if type(X) is float: 
        f_x = np.cos(np.pi*X)
    elif X.ndim == 1:
        f_x = np.cos(np.pi*X)
    else:
        f_x = np.prod(np.cos(np.pi*X), axis=-1)
    
    return f_x

def absolute(X):
    if type(X) is float: 
        f_x = np.abs(X)
    elif X.ndim == 1:
        f_x = np.abs(X)
    else:
        f_x = np.prod(np.abs(X), axis=-1)
    
    return f_x

def heaviside(X, shift=0):
    
    if type(X) is float: 
        f_x = np.heaviside(X - shift, 0)
    elif X.ndim == 1:
        f_x = np.heaviside(X - shift, 0)
    else:
        f_x = np.prod(np.heaviside(X - shift, 0), axis=-1)
    
    return f_x

def mix(X):
    
    return 0.5 * (heaviside(X) - 0.5) * (cosine(X) + cosine(2*X) + sine(5*X))<jupyter_output><empty_output><jupyter_text>Estimation<jupyter_code>def integrand(x):
    return bump_eps(x, eps=0.05)
#     return cosine(x)
#     return 2 * (heaviside(x) - 0.5)
#     return absolute(x)
#     return mix(x)
d, N = 1, 30

jac_params = -0.5 + np.random.rand(d, 2)
# jac_params = -0.5*np.ones((d, 2))

dpp = MultivariateJacobiOPE(N, jac_params)
sampl = dpp.sample()

print('Estimation of the integral\n')
for lab, est in zip(['BH', 'EZ'], both_estimators(integrand, dpp, sampl)):
    print(lab)
    print(est)
    
if d == 1:
    
    print('scipy quad')
    print(quad(lambda x: dpp.eval_w(x)*integrand(x), 
               -1, 1)[0])

    tol = 1e-4
    X_ = np.linspace(-1 + tol, 1 - tol, 300)[:, None]

    print('numpy trapz')
    print(np.trapz(dpp.eval_w(X_)*integrand(X_),
                   X_.ravel()))
    
    
    tols = np.zeros_like(dpp.jacobi_params)
    tols[tols < 0] = 5e-2
    X_ = np.linspace(-1 + tols[0, 1], 1 - tols[0, 0], 300)[:, None]
    fig, ax = plt.subplots()
    ax.plot(X_, integrand(X_), label='f')
    ax.scatter(sampl, np.zeros_like(sampl), label='sample')
    ax.scatter(sampl, integrand(sampl), label='f(sample)')
    
    plt.figure()
    plt.title(r'Base measure $\propto (1-x)^a (1+x)^b$')
    plt.plot(X_, 0.5*stats.beta(*(1+jac_params[0])).pdf(0.5*(1-X_)),
             label='\n'.join([r'$a \approx {:1.3f}$'.format(jac_params[0, 0]),
                              r'$b \approx {:1.3f}$'.format(jac_params[0, 1])]))
    plt.ylim(bottom=0)
    plt.legend(fontsize=15, loc='best')<jupyter_output><empty_output><jupyter_text>Variance decay  To repeat the estimations, we use the package `multiprocessing`  In this notebook, to estimate the variance of both BH and EZ estimators, we draw $20$ samples with up to $N=100$ points for $d=1,2$ (by default) You can change the parameters, but sampling may take some time  Of an integrand $f$<jupyter_code>def integrand(x):
    return bump_eps(x, eps=0.05)
#     return cosine(x)
#     return 2 * (heaviside(x) - 0.5)
#     return absolute(x)
#     return mix(x)
dim_max = 2
nb_repeats = 20

var_results = dict()

for d in range(1, dim_max+1):
    print('dimension =', d)
    
    jac_params = -0.5 + np.random.rand(d, 2)
    jac_params[0, :] = -0.5
    if d == 1:
        N_min, N_max, N_step = 20, 100, 20
    else:
        N_min, N_max, N_step = 20, 100, 20
    
    var_results[(d,)] = jac_params
    
    for N in range(N_min, N_max+1, N_step):
        print('#points =', N)

        dpp = MultivariateJacobiOPE(N, jac_params)

        pool = mp.Pool(mp.cpu_count())

        results = pool.starmap(both_estimators, [(integrand, dpp) for _ in range(nb_repeats)])
        results = np.array(results)
        var_results[(d, N)] = np.var(results, axis=0)

        pool.close()
j_par = {d: var_results.get((d, )) for d in range(1, dim_max + 1)}
var_N = {d: np.array([key[1] for key in var_results.keys() if len(key)==2 and key[0]==d])
            for d in range(1, dim_max + 1)}
var_res = {d: np.array([value for key, value in var_results.items() if len(key)==2 and key[0]==d]).T
            for d in range(1, dim_max + 1)}

cols = ['blue', 'green']
CB_cols = ['#377eb8', '#ff7f00', '#4daf4a',
                  '#f781bf', '#a65628', '#984ea3',
                  '#999999', '#e41a1c', '#dede00']
markers = ['o', '^']
labels = ['BH', 'EZ']
for d in range(1, dim_max + 1):
    
    fig, ax = plt.subplots()
#     plt.title(r'Dimension $d={}$'.format(d), fontsize=20)
    ax.set_xlabel(r'$N$', fontsize=22)
    ax.xaxis.set_label_coords(0.98, -0.025)
    ax.set_ylabel(r'$\mathrm{\mathbb{V}}$ar', fontsize=22, rotation='horizontal')
    ax.yaxis.set_label_coords(-0.06, 0.94)
    ax.tick_params(axis = 'both', which = 'major', labelsize = 17.5)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)


    for c, m, lab, var_estim in zip(CB_cols[:2], markers, labels, var_res[d]):

        ax.loglog(var_N[d], var_estim, m, c=c, markersize=8)

        x_plot = np.array([np.min(var_N[d]), np.max(var_N[d])])
        
        slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(var_N[d]), np.log(var_estim))
        lab += r' {:.1f}, {:.2f}'.format(slope, r_value**2)
        
        ax.loglog(x_plot, np.exp(intercept)*x_plot**slope, c=c, label=lab)
    
    leg = ax.legend(fontsize=20, frameon=False, handlelength=0.6, loc='lower left')

    for line in leg.get_lines():
        line.set_linewidth(4.0)

    plt.show()<jupyter_output><empty_output><jupyter_text>$f(x) = \sum_{k=0}^{M-1} \frac{1}{k+1} P_k(x)$ $M=70$EZ provides perfect estimation when $N\geq M$, see the drop in the variance plot<jupyter_code>dim_max = 2
nb_repeats = 20

M = 70
N_min, N_max, N_step = 10, 100, 10

var_results = dict()

for d in range(1, dim_max+1):
    print('dimension =', d)
    
    jac_params = -0.5 + np.random.rand(d, 2)
    jac_params[0, :] = -0.5
    
    dpp_gp = MultivariateJacobiOPE(M, jac_params)
    coefs = 1.0 / np.arange(1, dpp_gp.N + 1)
    
    def f_gp(X):
        return np.sum(coefs*dpp_gp.eval_poly_multiD(X, normalize='norm'), axis=-1)

    var_results[(d, )] = jac_params
    
    for N in range(N_min, N_max+1, N_step):
        print('#points =', N)

        dpp = MultivariateJacobiOPE(N, jac_params)

        pool = mp.Pool(mp.cpu_count())

        results = pool.starmap(both_estimators, [(f_gp, dpp) for _ in range(nb_repeats)])
        results = np.array(results)
        var_results[(d, N)] = np.var(results, axis=0)

        pool.close()
j_par = {d: var_results.get((d, )) for d in range(1, dim_max + 1)}
var_N = {d: np.array([key[1] for key in var_results.keys() if len(key)==2 and key[0]==d])
            for d in range(1, dim_max + 1)}
var_res = {d: np.array([value for key, value in var_results.items() if len(key)==2 and key[0]==d]).T
            for d in range(1, dim_max + 1)}


CB_cols = ['#377eb8', '#ff7f00', '#4daf4a',
                  '#f781bf', '#a65628', '#984ea3',
                  '#999999', '#e41a1c', '#dede00']
markers = ['o', '^']
labels = ['BH', 'EZ']
for d in range(1, dim_max + 1):
    
    fig, ax = plt.subplots()
    plt.title(r'Dimension $d={}$, $M$ = {}'.format(d, M), fontsize=20)
    ax.set_xlabel(r'$N$', fontsize=22)
    ax.xaxis.set_label_coords(1.03, -0.0)
    ax.set_ylabel(r'$\mathrm{\mathbb{V}}$ar', fontsize=22, rotation='horizontal')
    ax.yaxis.set_label_coords(-0.06, 0.95)
    ax.tick_params(axis = 'both', which = 'major', labelsize = 17.5)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)

    for c, m, lab, var_estim in zip(CB_cols[:2], markers, labels, var_res[d]):

        ax.loglog(var_N[d], var_estim, m, c=c, markersize=8)

        x_plot = np.array([np.min(var_N[d]), np.max(var_N[d])])
        
        if lab == 'BH':
            slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(var_N[d]), np.log(var_estim))
            lab += r' {:.1f}, {:.2f}'.format(slope, r_value**2)
        
            ax.loglog(x_plot, np.exp(intercept)*x_plot**slope, c=c, label=lab)
        
    leg = ax.legend(fontsize=20, frameon=False, handlelength=0.6, loc='lower left')

    for line in leg.get_lines():
        line.set_linewidth(4.0)
        
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16436158019408573, 'fraction_numerical': 0.031175976159547644, 'mean_word_length': 3.10282131661442, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 45, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '142466', 'license': 'MIT', 'n_tokens_mistral': 5569, 'n_tokens_neox': 5228, 'n_words': 1100, 'path': 'notebooks/multivariate_jacobi_ensemble.ipynb', 'repo_name': 'derezin/DPPy', 'size': '13087'}