Vectorized QMC¶
In [1]:
Copied!
%%capture
# @title Execute this cell to install dependancies
try:
import google.colab
import os
!pip install -q qmcpy >> /dev/null
!apt-get update && apt-get install -y --no-install-recommends texlive-latex-base texlive-fonts-recommended texlive-latex-extra cm-super dvipng
except:
pass
import matplotlib.pyplot as plt
plt.rcParams.update({
"text.usetex": True,
"font.family": "serif",
"text.latex.preamble": r"\usepackage{amsmath}\usepackage{amssymb}\newcommand{\bx}{\boldsymbol{x}}"
})
%%capture
# @title Execute this cell to install dependancies
try:
import google.colab
import os
!pip install -q qmcpy >> /dev/null
!apt-get update && apt-get install -y --no-install-recommends texlive-latex-base texlive-fonts-recommended texlive-latex-extra cm-super dvipng
except:
pass
import matplotlib.pyplot as plt
plt.rcParams.update({
"text.usetex": True,
"font.family": "serif",
"text.latex.preamble": r"\usepackage{amsmath}\usepackage{amssymb}\newcommand{\bx}{\boldsymbol{x}}"
})
In [2]:
Copied!
import qmcpy as qp
import numpy as np
import qmcpy as qp
import numpy as np
In [3]:
Copied!
from matplotlib import pyplot
%matplotlib inline
root = None#'./'
from matplotlib import pyplot
%matplotlib inline
root = None#'./'
LD Sequence¶
In [4]:
Copied!
n = 2**6
s = 10
fig,ax = pyplot.subplots(figsize=(8,4),nrows=1,ncols=3)
for i,(dd,name) in enumerate(zip(
[qp.IIDStdUniform(2,seed=7),qp.DigitalNetB2(2,seed=7),qp.Lattice(2,seed=7)],
['IID','Randomized Digital Net (LD)','Randomized Lattice (LD)'])):
pts = dd.gen_samples(n)
ax[i].scatter(pts[0:n//4,0],pts[0:n//4,1],color='k',marker='s',s=s)
ax[i].scatter(pts[n//4:n//2,0],pts[n//4:n//2,1],color='k',marker='o',s=s)
ax[i].scatter(pts[n//2:n,0],pts[n//2:n,1],color='k',marker='^',s=s)
ax[i].set_aspect(1)
ax[i].set_xlabel(r'$x_{1}$')
ax[i].set_ylabel(r'$x_{2}$')
ax[i].set_xlim([0,1])
ax[i].set_ylim([0,1])
ax[i].set_xticks([0,.25,.5,.75,1])
ax[i].set_xticklabels([r'$0$',r'$1/4$',r'$1/2$',r'$3/4$',r'$1$'])
ax[i].set_yticks([0,.25,.5,.75,1])
ax[i].set_yticklabels([r'$0$',r'$1/4$',r'$1/2$',r'$3/4$',r'$1$'])
ax[i].set_title(name)
if root: fig.savefig(root+'ld_seqs.pdf',transparent=True)
n = 2**6
s = 10
fig,ax = pyplot.subplots(figsize=(8,4),nrows=1,ncols=3)
for i,(dd,name) in enumerate(zip(
[qp.IIDStdUniform(2,seed=7),qp.DigitalNetB2(2,seed=7),qp.Lattice(2,seed=7)],
['IID','Randomized Digital Net (LD)','Randomized Lattice (LD)'])):
pts = dd.gen_samples(n)
ax[i].scatter(pts[0:n//4,0],pts[0:n//4,1],color='k',marker='s',s=s)
ax[i].scatter(pts[n//4:n//2,0],pts[n//4:n//2,1],color='k',marker='o',s=s)
ax[i].scatter(pts[n//2:n,0],pts[n//2:n,1],color='k',marker='^',s=s)
ax[i].set_aspect(1)
ax[i].set_xlabel(r'$x_{1}$')
ax[i].set_ylabel(r'$x_{2}$')
ax[i].set_xlim([0,1])
ax[i].set_ylim([0,1])
ax[i].set_xticks([0,.25,.5,.75,1])
ax[i].set_xticklabels([r'$0$',r'$1/4$',r'$1/2$',r'$3/4$',r'$1$'])
ax[i].set_yticks([0,.25,.5,.75,1])
ax[i].set_yticklabels([r'$0$',r'$1/4$',r'$1/2$',r'$3/4$',r'$1$'])
ax[i].set_title(name)
if root: fig.savefig(root+'ld_seqs.pdf',transparent=True)
Simple Example¶
In [5]:
Copied!
def cantilever_beam_function(T,compute_flags): # T is (n x 3)
Y = np.zeros((2,len(T)),dtype=float) # (n x 2)
l,w,t = 100,4,2
T1,T2,T3 = T[:,0],T[:,1],T[:,2] # Python is indexed from 0
if compute_flags[0]: # compute D. x^2 is "x**2" in Python
Y[0] = 4*l**3/(T1*w*t)*np.sqrt(T2**2/t**4+T3**2/w**4)
if compute_flags[1]: # compute S
Y[1] = 600*(T2/(w*t**2)+T3/(w**2*t))
return Y
true_measure = qp.Gaussian(
sampler = qp.DigitalNetB2(dimension=3,seed=7),
mean = [2.9e7,500,1000],
covariance = np.diag([(1.45e6)**2,(100)**2,(100)**2]))
integrand = qp.CustomFun(true_measure,
g = cantilever_beam_function,
dimension_indv = 2)
qmc_stop_crit = qp.CubQMCNetG(integrand,
abs_tol = 1e-3,
rel_tol = 1e-6)
solution,data = qmc_stop_crit.integrate()
print(solution)
# [2.42575885e+00 3.74999973e+04]
print(data)
def cantilever_beam_function(T,compute_flags): # T is (n x 3)
Y = np.zeros((2,len(T)),dtype=float) # (n x 2)
l,w,t = 100,4,2
T1,T2,T3 = T[:,0],T[:,1],T[:,2] # Python is indexed from 0
if compute_flags[0]: # compute D. x^2 is "x**2" in Python
Y[0] = 4*l**3/(T1*w*t)*np.sqrt(T2**2/t**4+T3**2/w**4)
if compute_flags[1]: # compute S
Y[1] = 600*(T2/(w*t**2)+T3/(w**2*t))
return Y
true_measure = qp.Gaussian(
sampler = qp.DigitalNetB2(dimension=3,seed=7),
mean = [2.9e7,500,1000],
covariance = np.diag([(1.45e6)**2,(100)**2,(100)**2]))
integrand = qp.CustomFun(true_measure,
g = cantilever_beam_function,
dimension_indv = 2)
qmc_stop_crit = qp.CubQMCNetG(integrand,
abs_tol = 1e-3,
rel_tol = 1e-6)
solution,data = qmc_stop_crit.integrate()
print(solution)
# [2.42575885e+00 3.74999973e+04]
print(data)
[2.42570423e+00 3.75000056e+04]
Data (Data)
solution [2.426e+00 3.750e+04]
comb_bound_low [2.425e+00 3.750e+04]
comb_bound_high [2.426e+00 3.750e+04]
comb_bound_diff [0.001 0.073]
comb_flags [ True True]
n_total 2^(17)
n [ 2048 131072]
time_integrate 0.215
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.001
rel_tol 1.00e-06
n_init 2^(10)
n_limit 2^(35)
CustomFun (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean [2.9e+07 5.0e+02 1.0e+03]
covariance [[2.102e+12 0.000e+00 0.000e+00]
[0.000e+00 1.000e+04 0.000e+00]
[0.000e+00 0.000e+00 1.000e+04]]
decomp_type PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 3
replications 1
randomize LMS_DS
gen_mats_source joe_kuo.6.21201.txt
order NATURAL
t 63
alpha 1
n_limit 2^(32)
entropy 7
BO QEI¶
See the QEI Demo in QMCPy or the BoTorch Acquisition documentation for details on Bayesian Optimization using q-Expected Improvement.
In [6]:
Copied!
import scipy
from sklearn.gaussian_process import GaussianProcessRegressor,kernels
f = lambda x: np.cos(10*x)*np.exp(.2*x)+np.exp(-5*(x-.4)**2)
xplt = np.linspace(0,1,100)
yplt = f(xplt)
x = np.array([.1, .2, .4, .7, .9])
y = f(x)
ymax = y.max()
gp = GaussianProcessRegressor(kernel=kernels.RBF(length_scale=1.0,length_scale_bounds=(1e-2, 1e2)),
n_restarts_optimizer = 16).fit(x[:,None],y)
yhatplt,stdhatplt = gp.predict(xplt[:,None],return_std=True)
tpax = 32
x0mesh,x1mesh = np.meshgrid(np.linspace(0,1,tpax),np.linspace(0,1,tpax))
post_mus = np.zeros((tpax,tpax,2),dtype=float)
post_sqrtcovs = np.zeros((tpax,tpax,2,2),dtype=float)
for j0 in range(tpax):
for j1 in range(tpax):
candidate = np.array([[x0mesh[j0,j1]],[x1mesh[j0,j1]]])
post_mus[j0,j1],post_cov = gp.predict(candidate,return_cov=True)
evals,evecs = scipy.linalg.eig(post_cov)
post_sqrtcovs[j0,j1] = np.sqrt(np.maximum(evals.real,0))*evecs
def qei_acq_vec(x,compute_flags):
xgauss = scipy.stats.norm.ppf(x)
n = len(x)
qei_vals = np.zeros((tpax,tpax,n),dtype=float)
for j0 in range(tpax):
for j1 in range(tpax):
if compute_flags[j0,j1]==False: continue
sqrt_cov = post_sqrtcovs[j0,j1]
mu_post = post_mus[j0,j1]
for i in range(len(x)):
yij = sqrt_cov@xgauss[i]+mu_post
qei_vals[j0,j1,i] = max((yij-ymax).max(),0)
return qei_vals
qei_acq_vec_qmcpy = qp.CustomFun(
true_measure = qp.Uniform(qp.DigitalNetB2(2,seed=7)),
g = qei_acq_vec,
dimension_indv = (tpax,tpax),
parallel=False)
qei_vals,qei_data = qp.CubQMCNetG(qei_acq_vec_qmcpy,abs_tol=.025,rel_tol=0).integrate() # .0005
print(qei_data)
a = np.unravel_index(np.argmax(qei_vals,axis=None),qei_vals.shape)
xnext = np.array([x0mesh[a[0],a[1]],x1mesh[a[0],a[1]]])
fnext = f(xnext)
import scipy
from sklearn.gaussian_process import GaussianProcessRegressor,kernels
f = lambda x: np.cos(10*x)*np.exp(.2*x)+np.exp(-5*(x-.4)**2)
xplt = np.linspace(0,1,100)
yplt = f(xplt)
x = np.array([.1, .2, .4, .7, .9])
y = f(x)
ymax = y.max()
gp = GaussianProcessRegressor(kernel=kernels.RBF(length_scale=1.0,length_scale_bounds=(1e-2, 1e2)),
n_restarts_optimizer = 16).fit(x[:,None],y)
yhatplt,stdhatplt = gp.predict(xplt[:,None],return_std=True)
tpax = 32
x0mesh,x1mesh = np.meshgrid(np.linspace(0,1,tpax),np.linspace(0,1,tpax))
post_mus = np.zeros((tpax,tpax,2),dtype=float)
post_sqrtcovs = np.zeros((tpax,tpax,2,2),dtype=float)
for j0 in range(tpax):
for j1 in range(tpax):
candidate = np.array([[x0mesh[j0,j1]],[x1mesh[j0,j1]]])
post_mus[j0,j1],post_cov = gp.predict(candidate,return_cov=True)
evals,evecs = scipy.linalg.eig(post_cov)
post_sqrtcovs[j0,j1] = np.sqrt(np.maximum(evals.real,0))*evecs
def qei_acq_vec(x,compute_flags):
xgauss = scipy.stats.norm.ppf(x)
n = len(x)
qei_vals = np.zeros((tpax,tpax,n),dtype=float)
for j0 in range(tpax):
for j1 in range(tpax):
if compute_flags[j0,j1]==False: continue
sqrt_cov = post_sqrtcovs[j0,j1]
mu_post = post_mus[j0,j1]
for i in range(len(x)):
yij = sqrt_cov@xgauss[i]+mu_post
qei_vals[j0,j1,i] = max((yij-ymax).max(),0)
return qei_vals
qei_acq_vec_qmcpy = qp.CustomFun(
true_measure = qp.Uniform(qp.DigitalNetB2(2,seed=7)),
g = qei_acq_vec,
dimension_indv = (tpax,tpax),
parallel=False)
qei_vals,qei_data = qp.CubQMCNetG(qei_acq_vec_qmcpy,abs_tol=.025,rel_tol=0).integrate() # .0005
print(qei_data)
a = np.unravel_index(np.argmax(qei_vals,axis=None),qei_vals.shape)
xnext = np.array([x0mesh[a[0],a[1]],x1mesh[a[0],a[1]]])
fnext = f(xnext)
Data (Data)
solution [[0.06 0.079 0.072 ... 0.06 0.061 0.067]
[0.08 0.064 0.067 ... 0.064 0.065 0.071]
[0.073 0.067 0.032 ... 0.032 0.033 0.039]
...
[0.06 0.064 0.032 ... 0. 0.001 0.007]
[0.061 0.065 0.033 ... 0.001 0.001 0.007]
[0.067 0.071 0.039 ... 0.007 0.007 0.007]]
comb_bound_low [[0.059 0.078 0.071 ... 0.059 0.06 0.066]
[0.078 0.063 0.066 ... 0.064 0.064 0.07 ]
[0.071 0.066 0.032 ... 0.032 0.033 0.038]
...
[0.059 0.063 0.032 ... 0. 0.001 0.006]
[0.06 0.064 0.033 ... 0.001 0.001 0.006]
[0.065 0.07 0.038 ... 0.006 0.006 0.006]]
comb_bound_high [[0.061 0.081 0.074 ... 0.061 0.062 0.068]
[0.081 0.065 0.068 ... 0.065 0.066 0.072]
[0.074 0.068 0.033 ... 0.033 0.034 0.04 ]
...
[0.061 0.065 0.033 ... 0. 0.001 0.008]
[0.062 0.066 0.034 ... 0.001 0.001 0.008]
[0.068 0.072 0.04 ... 0.008 0.008 0.008]]
comb_bound_diff [[0.002 0.002 0.003 ... 0.002 0.002 0.003]
[0.002 0.002 0.002 ... 0.002 0.002 0.002]
[0.002 0.002 0.001 ... 0.001 0.001 0.002]
...
[0.002 0.002 0.001 ... 0. 0. 0.002]
[0.002 0.002 0.001 ... 0. 0. 0.002]
[0.003 0.002 0.002 ... 0.002 0.002 0.002]]
comb_flags [[ True True True ... True True True]
[ True True True ... True True True]
[ True True True ... True True True]
...
[ True True True ... True True True]
[ True True True ... True True True]
[ True True True ... True True True]]
n_total 2^(10)
n [[1024 1024 1024 ... 1024 1024 1024]
[1024 1024 1024 ... 1024 1024 1024]
[1024 1024 1024 ... 1024 1024 1024]
...
[1024 1024 1024 ... 1024 1024 1024]
[1024 1024 1024 ... 1024 1024 1024]
[1024 1024 1024 ... 1024 1024 1024]]
time_integrate 7.262
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.025
rel_tol 0
n_init 2^(10)
n_limit 2^(35)
CustomFun (AbstractIntegrand)
Uniform (AbstractTrueMeasure)
lower_bound 0
upper_bound 1
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 2^(1)
replications 1
randomize LMS_DS
gen_mats_source joe_kuo.6.21201.txt
order NATURAL
t 63
alpha 1
n_limit 2^(32)
entropy 7
In [7]:
Copied!
from matplotlib import cm
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(8,3.25))
ax[0].scatter(x,y,color='k',label='Query Points')
ax[0].plot(xplt,yplt,color='k',linestyle='--',label='True function',linewidth=1)
ax[0].plot(xplt,yhatplt,color='k',label='GP Mean',linewidth=1)
ax[0].fill_between(xplt,yhatplt-1.96*stdhatplt,yhatplt+1.96*stdhatplt,color='k',alpha=.25,label='95% CI')
ax[0].scatter(xnext,fnext,color='k',marker='*',s=200,zorder=10)
ax[0].set_xlim([0,1])
ax[0].set_xticks([0,1])
ax[0].set_xlabel(r'$x$')
ax[0].set_ylabel(r'$y$')
fig.legend(labels=['data','true function','posterior mean',r'95\% CI','next points by qEI'],loc='lower center',bbox_to_anchor=(.5,-.05),ncol=5)
contour = ax[1].contourf(x0mesh,x1mesh,qei_vals,cmap=cm.Greys_r)
ax[1].scatter([xnext[0]],[xnext[1]],color='k',marker='*',s=200)
fig.colorbar(contour,ax=None,shrink=1,aspect=5)
ax[1].scatter(x0mesh.flatten(),x1mesh.flatten(),color='w',s=1)
ax[1].set_xlim([0,1])
ax[1].set_xticks([0,1])
ax[1].set_ylim([0,1])
ax[1].set_yticks([0,1])
ax[1].set_xlabel(r'$x_1$')
ax[1].set_ylabel(r'$x_2$')
if root: fig.savefig(root+'gp.pdf',transparent=True)
from matplotlib import cm
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(8,3.25))
ax[0].scatter(x,y,color='k',label='Query Points')
ax[0].plot(xplt,yplt,color='k',linestyle='--',label='True function',linewidth=1)
ax[0].plot(xplt,yhatplt,color='k',label='GP Mean',linewidth=1)
ax[0].fill_between(xplt,yhatplt-1.96*stdhatplt,yhatplt+1.96*stdhatplt,color='k',alpha=.25,label='95% CI')
ax[0].scatter(xnext,fnext,color='k',marker='*',s=200,zorder=10)
ax[0].set_xlim([0,1])
ax[0].set_xticks([0,1])
ax[0].set_xlabel(r'$x$')
ax[0].set_ylabel(r'$y$')
fig.legend(labels=['data','true function','posterior mean',r'95\% CI','next points by qEI'],loc='lower center',bbox_to_anchor=(.5,-.05),ncol=5)
contour = ax[1].contourf(x0mesh,x1mesh,qei_vals,cmap=cm.Greys_r)
ax[1].scatter([xnext[0]],[xnext[1]],color='k',marker='*',s=200)
fig.colorbar(contour,ax=None,shrink=1,aspect=5)
ax[1].scatter(x0mesh.flatten(),x1mesh.flatten(),color='w',s=1)
ax[1].set_xlim([0,1])
ax[1].set_xticks([0,1])
ax[1].set_ylim([0,1])
ax[1].set_yticks([0,1])
ax[1].set_xlabel(r'$x_1$')
ax[1].set_ylabel(r'$x_2$')
if root: fig.savefig(root+'gp.pdf',transparent=True)
Bayesian Logistic Regression¶
In [8]:
Copied!
import pandas as pd
from sklearn.model_selection import train_test_split
df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/haberman/haberman.data',header=None)
df.columns = ['Age','1900 Year','Axillary Nodes','Survival Status']
df.loc[df['Survival Status']==2,'Survival Status'] = 0
x,y = df[['Age','1900 Year','Axillary Nodes']],df['Survival Status']
xt,xv,yt,yv = train_test_split(x,y,test_size=.33,random_state=7)
import pandas as pd
from sklearn.model_selection import train_test_split
df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/haberman/haberman.data',header=None)
df.columns = ['Age','1900 Year','Axillary Nodes','Survival Status']
df.loc[df['Survival Status']==2,'Survival Status'] = 0
x,y = df[['Age','1900 Year','Axillary Nodes']],df['Survival Status']
xt,xv,yt,yv = train_test_split(x,y,test_size=.33,random_state=7)
In [9]:
Copied!
print(df.head(),'\n')
print(df[['Age','1900 Year','Axillary Nodes']].describe(),'\n')
print(df['Survival Status'].astype(str).describe())
print('\ntrain samples: %d test samples: %d\n'%(len(xt),len(xv)))
print('train positives %d train negatives: %d'%(np.sum(yt==1),np.sum(yt==0)))
print(' test positives %d test negatives: %d'%(np.sum(yv==1),np.sum(yv==0)))
xt.head()
print(df.head(),'\n')
print(df[['Age','1900 Year','Axillary Nodes']].describe(),'\n')
print(df['Survival Status'].astype(str).describe())
print('\ntrain samples: %d test samples: %d\n'%(len(xt),len(xv)))
print('train positives %d train negatives: %d'%(np.sum(yt==1),np.sum(yt==0)))
print(' test positives %d test negatives: %d'%(np.sum(yv==1),np.sum(yv==0)))
xt.head()
Age 1900 Year Axillary Nodes Survival Status
0 30 64 1 1
1 30 62 3 1
2 30 65 0 1
3 31 59 2 1
4 31 65 4 1
Age 1900 Year Axillary Nodes
count 306.000000 306.000000 306.000000
mean 52.457516 62.852941 4.026144
std 10.803452 3.249405 7.189654
min 30.000000 58.000000 0.000000
25% 44.000000 60.000000 0.000000
50% 52.000000 63.000000 1.000000
75% 60.750000 65.750000 4.000000
max 83.000000 69.000000 52.000000
count 306
unique 2
top 1
freq 225
Name: Survival Status, dtype: object
train samples: 205 test samples: 101
train positives 151 train negatives: 54
test positives 74 test negatives: 27
Out[9]:
| Age | 1900 Year | Axillary Nodes | |
|---|---|---|---|
| 46 | 41 | 58 | 0 |
| 199 | 57 | 64 | 1 |
| 115 | 49 | 64 | 10 |
| 128 | 50 | 61 | 0 |
| 249 | 63 | 63 | 0 |
In [10]:
Copied!
blr = qp.BayesianLRCoeffs(
sampler = qp.DigitalNetB2(4,seed=1),
feature_array = xt, # np.ndarray of shape (n,d-1)
response_vector = yt, # np.ndarray of shape (n,)
prior_mean = 0, # normal prior mean = (0,0,...,0)
prior_covariance = 5) # normal prior covariance = 5I
qmc_sc = qp.CubQMCNetG(blr,
abs_tol = .1,
rel_tol = .5,
error_fun = "BOTH",
n_limit=2**18)
blr_coefs,blr_data = qmc_sc.integrate()
print(blr_data)
# LDTransformData (AccumulateData Object)
# solution [-0.004 0.13 -0.157 0.008]
# comb_bound_low [-0.006 0.092 -0.205 0.007]
# comb_bound_high [-0.003 0.172 -0.109 0.012]
# comb_flags [ True True True True]
# n_total 2^(18)
# n [[ 1024. 1024. 262144. 2048.]
# [ 1024. 1024. 262144. 2048.]]
# time_integrate 2.229
blr = qp.BayesianLRCoeffs(
sampler = qp.DigitalNetB2(4,seed=1),
feature_array = xt, # np.ndarray of shape (n,d-1)
response_vector = yt, # np.ndarray of shape (n,)
prior_mean = 0, # normal prior mean = (0,0,...,0)
prior_covariance = 5) # normal prior covariance = 5I
qmc_sc = qp.CubQMCNetG(blr,
abs_tol = .1,
rel_tol = .5,
error_fun = "BOTH",
n_limit=2**18)
blr_coefs,blr_data = qmc_sc.integrate()
print(blr_data)
# LDTransformData (AccumulateData Object)
# solution [-0.004 0.13 -0.157 0.008]
# comb_bound_low [-0.006 0.092 -0.205 0.007]
# comb_bound_high [-0.003 0.172 -0.109 0.012]
# comb_flags [ True True True True]
# n_total 2^(18)
# n [[ 1024. 1024. 262144. 2048.]
# [ 1024. 1024. 262144. 2048.]]
# time_integrate 2.229
Data (Data)
solution [ 0.262 -0.043 -0.226 -1.203]
comb_bound_low [ 0.185 -0.067 -0.306 -1.656]
comb_bound_high [ 0.347 -0.035 -0.163 -0.75 ]
comb_bound_diff [0.162 0.031 0.143 0.906]
comb_flags [ True True True False]
n_total 2^(18)
n [[ 1024 1024 1024 262144]
[ 1024 1024 1024 262144]]
time_integrate 2.856
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.100
rel_tol 2^(-1)
n_init 2^(10)
n_limit 2^(18)
BayesianLRCoeffs (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean 0
covariance 5
decomp_type PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 2^(2)
replications 1
randomize LMS_DS
gen_mats_source joe_kuo.6.21201.txt
order NATURAL
t 63
alpha 1
n_limit 2^(32)
entropy 1
/Users/alegresor/Desktop/QMCSoftware/qmcpy/stopping_criterion/abstract_cub_qmc_ld_g.py:215: MaxSamplesWarning:
Already generated 262144 samples.
Trying to generate 262144 new samples would exceeds n_limit = 262144.
No more samples will be generated.
Note that error tolerances may not be satisfied.
warnings.warn(warning_s, MaxSamplesWarning)
In [11]:
Copied!
from sklearn.linear_model import LogisticRegression
def metrics(y,yhat):
y,yhat = np.array(y),np.array(yhat)
tp = np.sum((y==1)*(yhat==1))
tn = np.sum((y==0)*(yhat==0))
fp = np.sum((y==0)*(yhat==1))
fn = np.sum((y==1)*(yhat==0))
accuracy = (tp+tn)/(len(y))
precision = tp/(tp+fp)
recall = tp/(tp+fn)
return [accuracy,precision,recall]
results = pd.DataFrame({name:[] for name in ['method','Age','1900 Year','Axillary Nodes','Intercept','Accuracy','Precision','Recall']})
for i,l1_ratio in enumerate([0,.5,1]):
lr = LogisticRegression(random_state=7,penalty="elasticnet",solver='saga',l1_ratio=l1_ratio).fit(xt,yt)
results.loc[i] = [r'Elastic-Net \lambda=%.1f'%l1_ratio]+lr.coef_.squeeze().tolist()+[lr.intercept_.item()]+metrics(yv,lr.predict(xv))
blr_predict = lambda x: 1/(1+np.exp(-np.array(x)@blr_coefs[:-1]-blr_coefs[-1]))>=.5
blr_train_accuracy = np.mean(blr_predict(xt)==yt)
blr_test_accuracy = np.mean(blr_predict(xv)==yv)
results.loc[len(results)] = ['Bayesian']+blr_coefs.squeeze().tolist()+metrics(yv,blr_predict(xv))
import warnings
warnings.simplefilter('ignore',FutureWarning)
results.set_index('method',inplace=True)
print(results.head())
#root: results.to_latex(root+'lr_table.tex',formatters={'%s'%tt:lambda v:'%.1f'%(100*v) for tt in ['accuracy','precision','recall']},float_format="%.2e")
from sklearn.linear_model import LogisticRegression
def metrics(y,yhat):
y,yhat = np.array(y),np.array(yhat)
tp = np.sum((y==1)*(yhat==1))
tn = np.sum((y==0)*(yhat==0))
fp = np.sum((y==0)*(yhat==1))
fn = np.sum((y==1)*(yhat==0))
accuracy = (tp+tn)/(len(y))
precision = tp/(tp+fp)
recall = tp/(tp+fn)
return [accuracy,precision,recall]
results = pd.DataFrame({name:[] for name in ['method','Age','1900 Year','Axillary Nodes','Intercept','Accuracy','Precision','Recall']})
for i,l1_ratio in enumerate([0,.5,1]):
lr = LogisticRegression(random_state=7,penalty="elasticnet",solver='saga',l1_ratio=l1_ratio).fit(xt,yt)
results.loc[i] = [r'Elastic-Net \lambda=%.1f'%l1_ratio]+lr.coef_.squeeze().tolist()+[lr.intercept_.item()]+metrics(yv,lr.predict(xv))
blr_predict = lambda x: 1/(1+np.exp(-np.array(x)@blr_coefs[:-1]-blr_coefs[-1]))>=.5
blr_train_accuracy = np.mean(blr_predict(xt)==yt)
blr_test_accuracy = np.mean(blr_predict(xv)==yv)
results.loc[len(results)] = ['Bayesian']+blr_coefs.squeeze().tolist()+metrics(yv,blr_predict(xv))
import warnings
warnings.simplefilter('ignore',FutureWarning)
results.set_index('method',inplace=True)
print(results.head())
#root: results.to_latex(root+'lr_table.tex',formatters={'%s'%tt:lambda v:'%.1f'%(100*v) for tt in ['accuracy','precision','recall']},float_format="%.2e")
Age 1900 Year Axillary Nodes Intercept \
method
Elastic-Net \lambda=0.0 -0.012279 0.034401 -0.115153 0.001990
Elastic-Net \lambda=0.5 -0.012041 0.034170 -0.114770 0.002025
Elastic-Net \lambda=1.0 -0.011803 0.033940 -0.114387 0.002061
Bayesian 0.262086 -0.043253 -0.225631 -1.202777
Accuracy Precision Recall
method
Elastic-Net \lambda=0.0 0.742574 0.766667 0.932432
Elastic-Net \lambda=0.5 0.742574 0.766667 0.932432
Elastic-Net \lambda=1.0 0.742574 0.766667 0.932432
Bayesian 0.732673 0.737374 0.986486
Sensitivity Indices¶
Ishigami Function¶
In [12]:
Copied!
a,b = 7,0.1
dnb2 = qp.DigitalNetB2(3,seed=7)
ishigami = qp.Ishigami(dnb2,a,b)
idxs = np.array([
[True,False,False],
[False,True,False],
[False,False,True],
[True,True,False],
[True,False,True],
[False,True,True]],dtype=bool)
ishigami_si = qp.SensitivityIndices(ishigami,idxs)
qmc_algo = qp.CubQMCNetG(ishigami_si,abs_tol=.05)
solution,data = qmc_algo.integrate()
print(data)
si_closed = solution[0].squeeze()
si_total = solution[1].squeeze()
ci_comb_low_closed = data.comb_bound_low[0].squeeze()
ci_comb_high_closed = data.comb_bound_high[0].squeeze()
ci_comb_low_total = data.comb_bound_low[1].squeeze()
ci_comb_high_total = data.comb_bound_high[1].squeeze()
print("\nApprox took %.1f sec and n = 2^(%d)"%
(data.time_integrate,np.log2(data.n_total)))
print('\t si_closed:',si_closed)
print('\t si_total:',si_total)
print('\t ci_comb_low_closed:',ci_comb_low_closed)
print('\t ci_comb_high_closed:',ci_comb_high_closed)
print('\t ci_comb_low_total:',ci_comb_low_total)
print('\t ci_comb_high_total:',ci_comb_high_total)
true_indices = qp.Ishigami._exact_sensitivity_indices(idxs,a,b)
si_closed_true = true_indices[0]
si_total_true = true_indices[1]
a,b = 7,0.1
dnb2 = qp.DigitalNetB2(3,seed=7)
ishigami = qp.Ishigami(dnb2,a,b)
idxs = np.array([
[True,False,False],
[False,True,False],
[False,False,True],
[True,True,False],
[True,False,True],
[False,True,True]],dtype=bool)
ishigami_si = qp.SensitivityIndices(ishigami,idxs)
qmc_algo = qp.CubQMCNetG(ishigami_si,abs_tol=.05)
solution,data = qmc_algo.integrate()
print(data)
si_closed = solution[0].squeeze()
si_total = solution[1].squeeze()
ci_comb_low_closed = data.comb_bound_low[0].squeeze()
ci_comb_high_closed = data.comb_bound_high[0].squeeze()
ci_comb_low_total = data.comb_bound_low[1].squeeze()
ci_comb_high_total = data.comb_bound_high[1].squeeze()
print("\nApprox took %.1f sec and n = 2^(%d)"%
(data.time_integrate,np.log2(data.n_total)))
print('\t si_closed:',si_closed)
print('\t si_total:',si_total)
print('\t ci_comb_low_closed:',ci_comb_low_closed)
print('\t ci_comb_high_closed:',ci_comb_high_closed)
print('\t ci_comb_low_total:',ci_comb_low_total)
print('\t ci_comb_high_total:',ci_comb_high_total)
true_indices = qp.Ishigami._exact_sensitivity_indices(idxs,a,b)
si_closed_true = true_indices[0]
si_total_true = true_indices[1]
Data (Data)
solution [[0.314 0.443 0.004 0.756 0.567 0.44 ]
[0.582 0.443 0.245 0.984 0.567 0.694]]
comb_bound_low [[0.293 0.426 0. 0.722 0.539 0.418]
[0.552 0.432 0.236 0.969 0.544 0.67 ]]
comb_bound_high [[0.334 0.459 0.008 0.79 0.595 0.463]
[0.612 0.453 0.253 1. 0.59 0.718]]
comb_bound_diff [[0.041 0.033 0.008 0.069 0.056 0.046]
[0.06 0.021 0.017 0.031 0.045 0.048]]
comb_flags [[ True True True True True True]
[ True True True True True True]]
n_total 2^(10)
n [[[1024 1024 1024 1024 1024 1024]
[1024 1024 1024 1024 1024 1024]
[1024 1024 1024 1024 1024 1024]]
[[1024 1024 1024 1024 1024 1024]
[1024 1024 1024 1024 1024 1024]
[1024 1024 1024 1024 1024 1024]]]
time_integrate 0.021
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.050
rel_tol 0
n_init 2^(10)
n_limit 2^(35)
SensitivityIndices (AbstractIntegrand)
indices [[ True False False]
[False True False]
[False False True]
[ True True False]
[ True False True]
[False True True]]
Uniform (AbstractTrueMeasure)
lower_bound -3.142
upper_bound 3.142
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 3
replications 1
randomize LMS_DS
gen_mats_source joe_kuo.6.21201.txt
order NATURAL
t 63
alpha 1
n_limit 2^(32)
entropy 7
Approx took 0.0 sec and n = 2^(10)
si_closed: [0.3135849 0.44255327 0.00393187 0.75607985 0.56701027 0.44044942]
si_total: [0.58195085 0.44252811 0.24471329 0.98433863 0.56700193 0.69402006]
ci_comb_low_closed: [0.29324641 0.42598379 0. 0.72174769 0.53904217 0.4175601 ]
ci_comb_high_closed: [0.33392338 0.45912275 0.00786373 0.79041202 0.59497838 0.46333875]
ci_comb_low_total: [0.55205604 0.43181773 0.23639082 0.96867726 0.54432148 0.67010028]
ci_comb_high_total: [0.61184567 0.45323849 0.25303576 1. 0.58968238 0.71793983]
In [13]:
Copied!
fig,ax = pyplot.subplots(figsize=(6,3))
ax.grid(False)
for spine in ['top','left','right','bottom']: ax.spines[spine].set_visible(False)
width = .75
ax.errorbar(fmt='none',color='k',
x = 1-si_total_true,
y = np.arange(len(si_closed)),
xerr = 0,
yerr = width/2,
alpha = 1)
bar_closed = ax.barh(np.arange(len(si_closed)),np.flip(si_closed),width,label='Closed SI',color='w',edgecolor='k',alpha=.75,linestyle='--')
ax.errorbar(fmt='none',color='k',
x = si_closed,
y = np.flip(np.arange(len(si_closed)))+width/4,
xerr = np.vstack((si_closed-ci_comb_low_closed,ci_comb_high_closed-si_closed)),
yerr = 0,
#elinewidth = 5,
alpha = .75)
bar_total = ax.barh(np.arange(len(si_closed)),si_total,width,label='Total SI',color='w',alpha=.25,edgecolor='k',left=1-si_total,zorder=10,linestyle='-.')
ax.errorbar(fmt='none',color='k',
x = 1-si_total,
y = np.arange(len(si_closed))-width/4,
xerr = np.vstack((si_total-ci_comb_low_total,ci_comb_high_total-si_total)),
yerr = 0,
#elinewidth = 5,
alpha = .25)
closed_labels = [r'$\underline{s}_{\{%s\}} = %.2f$'%(','.join([str(i+1) for i in idx]),c) for idx,c in zip(idxs[::-1],np.flip(si_closed))]
closed_labels[3] = ''
total_labels = [r'$\overline{s}_{\{%s\}} = %.2f$'%(','.join([str(i+1) for i in idx]),t) for idx,t in zip(idxs,si_total)]
ax.bar_label(bar_closed,label_type='center',labels=closed_labels)
ax.bar_label(bar_total,label_type='center',labels=total_labels)
ax.set_xlim([-.001,1.001])
ax.axvline(x=0,ymin=0,ymax=len(si_closed),color='k',alpha=.25)
ax.axvline(x=1,ymin=0,ymax=len(si_closed),color='k',alpha=.25)
ax.set_yticklabels([])
if root: fig.savefig(root+'ishigami.pdf',transparent=True)
fig,ax = pyplot.subplots(figsize=(6,3))
ax.grid(False)
for spine in ['top','left','right','bottom']: ax.spines[spine].set_visible(False)
width = .75
ax.errorbar(fmt='none',color='k',
x = 1-si_total_true,
y = np.arange(len(si_closed)),
xerr = 0,
yerr = width/2,
alpha = 1)
bar_closed = ax.barh(np.arange(len(si_closed)),np.flip(si_closed),width,label='Closed SI',color='w',edgecolor='k',alpha=.75,linestyle='--')
ax.errorbar(fmt='none',color='k',
x = si_closed,
y = np.flip(np.arange(len(si_closed)))+width/4,
xerr = np.vstack((si_closed-ci_comb_low_closed,ci_comb_high_closed-si_closed)),
yerr = 0,
#elinewidth = 5,
alpha = .75)
bar_total = ax.barh(np.arange(len(si_closed)),si_total,width,label='Total SI',color='w',alpha=.25,edgecolor='k',left=1-si_total,zorder=10,linestyle='-.')
ax.errorbar(fmt='none',color='k',
x = 1-si_total,
y = np.arange(len(si_closed))-width/4,
xerr = np.vstack((si_total-ci_comb_low_total,ci_comb_high_total-si_total)),
yerr = 0,
#elinewidth = 5,
alpha = .25)
closed_labels = [r'$\underline{s}_{\{%s\}} = %.2f$'%(','.join([str(i+1) for i in idx]),c) for idx,c in zip(idxs[::-1],np.flip(si_closed))]
closed_labels[3] = ''
total_labels = [r'$\overline{s}_{\{%s\}} = %.2f$'%(','.join([str(i+1) for i in idx]),t) for idx,t in zip(idxs,si_total)]
ax.bar_label(bar_closed,label_type='center',labels=closed_labels)
ax.bar_label(bar_total,label_type='center',labels=total_labels)
ax.set_xlim([-.001,1.001])
ax.axvline(x=0,ymin=0,ymax=len(si_closed),color='k',alpha=.25)
ax.axvline(x=1,ymin=0,ymax=len(si_closed),color='k',alpha=.25)
ax.set_yticklabels([])
if root: fig.savefig(root+'ishigami.pdf',transparent=True)
In [14]:
Copied!
fig,ax = pyplot.subplots(nrows=1,ncols=3,figsize=(8,3))
x_1d = np.linspace(0,1,num=128)
x_1d_mat = np.tile(x_1d,(3,1)).T
y_1d = qp.Ishigami._exact_fu_functions(x_1d_mat,indices=[[0],[1],[2]],a=a,b=b)
for i in range(2):
ax[i].plot(x_1d,y_1d[:,i],color='k')
ax[i].set_xlim([0,1])
ax[i].set_xticks([0,1])
ax[i].set_xlabel(r'$x_{%d}$'%(i+1))
ax[i].set_title(r'$f_{\{%d\}} \in [%.1f,%.1f]$'%(i+1,y_1d[:,i].min(),y_1d[:,i].max()))
x_mesh,y_mesh = np.meshgrid(x_1d,x_1d)
xquery = np.zeros((x_mesh.size,3))
for i,idx in enumerate([[1,2]]): # [[0,1],[0,2],[1,2]]
xquery[:,idx[0]] = x_mesh.flatten()
xquery[:,idx[1]] = y_mesh.flatten()
zquery = qp.Ishigami._exact_fu_functions(xquery,indices=[idx],a=a,b=b)
z_mesh = zquery.reshape(x_mesh.shape)
ax[2+i].contourf(x_mesh,y_mesh,z_mesh,cmap=cm.Greys_r)
ax[2+i].set_xlabel(r'$x_{%d}$'%(idx[0]+1))
ax[2+i].set_ylabel(r'$x_{%d}$'%(idx[1]+1))
ax[2+i].set_title(r'$f_{\{%d,%d\}} \in [%.1f,%.1f]$'%(tuple([i+1 for i in idx])+(z_mesh.min(),z_mesh.max())))
ax[2+i].set_xlim([0,1])
ax[2+i].set_ylim([0,1])
ax[2+i].set_xticks([0,1])
ax[2+i].set_yticks([0,1])
if root: fig.savefig(root+'ishigami_fu.pdf')
fig,ax = pyplot.subplots(nrows=1,ncols=3,figsize=(8,3))
x_1d = np.linspace(0,1,num=128)
x_1d_mat = np.tile(x_1d,(3,1)).T
y_1d = qp.Ishigami._exact_fu_functions(x_1d_mat,indices=[[0],[1],[2]],a=a,b=b)
for i in range(2):
ax[i].plot(x_1d,y_1d[:,i],color='k')
ax[i].set_xlim([0,1])
ax[i].set_xticks([0,1])
ax[i].set_xlabel(r'$x_{%d}$'%(i+1))
ax[i].set_title(r'$f_{\{%d\}} \in [%.1f,%.1f]$'%(i+1,y_1d[:,i].min(),y_1d[:,i].max()))
x_mesh,y_mesh = np.meshgrid(x_1d,x_1d)
xquery = np.zeros((x_mesh.size,3))
for i,idx in enumerate([[1,2]]): # [[0,1],[0,2],[1,2]]
xquery[:,idx[0]] = x_mesh.flatten()
xquery[:,idx[1]] = y_mesh.flatten()
zquery = qp.Ishigami._exact_fu_functions(xquery,indices=[idx],a=a,b=b)
z_mesh = zquery.reshape(x_mesh.shape)
ax[2+i].contourf(x_mesh,y_mesh,z_mesh,cmap=cm.Greys_r)
ax[2+i].set_xlabel(r'$x_{%d}$'%(idx[0]+1))
ax[2+i].set_ylabel(r'$x_{%d}$'%(idx[1]+1))
ax[2+i].set_title(r'$f_{\{%d,%d\}} \in [%.1f,%.1f]$'%(tuple([i+1 for i in idx])+(z_mesh.min(),z_mesh.max())))
ax[2+i].set_xlim([0,1])
ax[2+i].set_ylim([0,1])
ax[2+i].set_xticks([0,1])
ax[2+i].set_yticks([0,1])
if root: fig.savefig(root+'ishigami_fu.pdf')
Neural Network¶
In [15]:
Copied!
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.neural_network import MLPClassifier
data = load_iris()
feature_names = data["feature_names"]
feature_names = [fn.replace('sepal ','S')\
.replace('length ','L')\
.replace('petal ','P')\
.replace('width ','W')\
.replace('(cm)','') for fn in feature_names]
target_names = data["target_names"]
xt,xv,yt,yv = train_test_split(data["data"],data["target"],
test_size = 1/3,
random_state = 7)
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.neural_network import MLPClassifier
data = load_iris()
feature_names = data["feature_names"]
feature_names = [fn.replace('sepal ','S')\
.replace('length ','L')\
.replace('petal ','P')\
.replace('width ','W')\
.replace('(cm)','') for fn in feature_names]
target_names = data["target_names"]
xt,xv,yt,yv = train_test_split(data["data"],data["target"],
test_size = 1/3,
random_state = 7)
In [16]:
Copied!
mlpc = MLPClassifier(random_state=7,max_iter=1024).fit(xt,yt)
yhat = mlpc.predict(xv)
print("accuracy: %.1f%%"%(100*(yv==yhat).mean()))
# accuracy: 98.0%
sampler = qp.DigitalNetB2(4,seed=7)
true_measure = qp.Uniform(sampler,
lower_bound = xt.min(0),
upper_bound = xt.max(0))
fun = qp.CustomFun(
true_measure = true_measure,
g = lambda x: mlpc.predict_proba(x).T,
dimension_indv = 3)
si_fun = qp.SensitivityIndices(fun,indices="all")
qmc_algo = qp.CubQMCNetG(si_fun,abs_tol=.005)
nn_sis,nn_sis_data = qmc_algo.integrate()
mlpc = MLPClassifier(random_state=7,max_iter=1024).fit(xt,yt)
yhat = mlpc.predict(xv)
print("accuracy: %.1f%%"%(100*(yv==yhat).mean()))
# accuracy: 98.0%
sampler = qp.DigitalNetB2(4,seed=7)
true_measure = qp.Uniform(sampler,
lower_bound = xt.min(0),
upper_bound = xt.max(0))
fun = qp.CustomFun(
true_measure = true_measure,
g = lambda x: mlpc.predict_proba(x).T,
dimension_indv = 3)
si_fun = qp.SensitivityIndices(fun,indices="all")
qmc_algo = qp.CubQMCNetG(si_fun,abs_tol=.005)
nn_sis,nn_sis_data = qmc_algo.integrate()
accuracy: 98.0%
In [17]:
Copied!
#print(nn_sis_data.flags_indv.shape)
#print(nn_sis_data.flags_comb.shape)
print('samples: 2^(%d)'%np.log2(nn_sis_data.n_total))
print('time: %.1e'%nn_sis_data.time_integrate)
print('indices:\n%s'%nn_sis_data.integrand.indices)
import pandas as pd
df_closed = pd.DataFrame(nn_sis[0],columns=target_names,index=[str(np.where(idx)[0]) for idx in nn_sis_data.integrand.indices])
print('\nClosed Indices')
print(df_closed)
df_total = pd.DataFrame(nn_sis[1],columns=target_names,index=[str(np.where(idx)[0]) for idx in nn_sis_data.integrand.indices])
print('\nTotal Indices')
df_closed_singletons = df_closed.loc[['[%d]'%i for i in range(4)]].T
df_closed_singletons['sum singletons'] = df_closed_singletons[['[%d]'%i for i in range(4)]].sum(1)
df_closed_singletons.columns = data['feature_names']+['sum']
df_closed_singletons = df_closed_singletons*100
df_closed_singletons
#print(nn_sis_data.flags_indv.shape)
#print(nn_sis_data.flags_comb.shape)
print('samples: 2^(%d)'%np.log2(nn_sis_data.n_total))
print('time: %.1e'%nn_sis_data.time_integrate)
print('indices:\n%s'%nn_sis_data.integrand.indices)
import pandas as pd
df_closed = pd.DataFrame(nn_sis[0],columns=target_names,index=[str(np.where(idx)[0]) for idx in nn_sis_data.integrand.indices])
print('\nClosed Indices')
print(df_closed)
df_total = pd.DataFrame(nn_sis[1],columns=target_names,index=[str(np.where(idx)[0]) for idx in nn_sis_data.integrand.indices])
print('\nTotal Indices')
df_closed_singletons = df_closed.loc[['[%d]'%i for i in range(4)]].T
df_closed_singletons['sum singletons'] = df_closed_singletons[['[%d]'%i for i in range(4)]].sum(1)
df_closed_singletons.columns = data['feature_names']+['sum']
df_closed_singletons = df_closed_singletons*100
df_closed_singletons
samples: 2^(15)
time: 1.5e+00
indices:
[[ True False False False]
[False True False False]
[False False True False]
[False False False True]
[ True True False False]
[ True False True False]
[ True False False True]
[False True True False]
[False True False True]
[False False True True]
[ True True True False]
[ True True False True]
[ True False True True]
[False True True True]]
Closed Indices
setosa versicolor virginica
[0] 0.000000 0.067568 0.099035
[1] 0.051979 0.002589 0.000000
[2] 0.712711 0.326664 0.498664
[3] 0.051251 0.018706 0.120796
[0 1] 0.052284 0.080725 0.112412
[0 2] 0.714204 0.462147 0.641450
[0 3] 0.050927 0.086925 0.207602
[1 2] 0.840864 0.433329 0.514042
[1 3] 0.102538 0.004948 0.126151
[2 3] 0.822405 0.583743 0.705192
[0 1 2] 0.843055 0.571913 0.661004
[0 1 3] 0.103351 0.103411 0.217246
[0 2 3] 0.824347 0.816010 0.946992
[1 2 3] 0.995490 0.739267 0.726782
Total Indices
Out[17]:
| sepal length (cm) | sepal width (cm) | petal length (cm) | petal width (cm) | sum | |
|---|---|---|---|---|---|
| setosa | 0.000000 | 5.197885 | 71.271100 | 5.125145 | 81.594130 |
| versicolor | 6.756844 | 0.258851 | 32.666361 | 1.870561 | 41.552616 |
| virginica | 9.903452 | 0.000000 | 49.866351 | 12.079616 | 71.849419 |
In [18]:
Copied!
nindices = len(nn_sis_data.integrand.indices)
fig,ax = pyplot.subplots(figsize=(9,5))
ticks = np.arange(nindices)
width = .25
for i,(alpha,species) in enumerate(zip([.25,.5,.75],data['target_names'])):
cvals = df_closed[species].to_numpy()
tvals = df_total[species].to_numpy()
ticks_i = ticks+i*width
ax.bar(ticks_i,cvals,width=width,align='edge',color='k',alpha=alpha,label=species)
#ax.bar(ticks_i,np.flip(tvals),width=width,align='edge',bottom=1-np.flip(tvals),color=color,alpha=.1)
ax.set_xlim([0,13+3*width])
ax.set_xticks(ticks+1.5*width)
# closed_labels = [r'$\underline{s}_{\{%s\}}$'%(','.join([r'\text{%s}'%feature_names[i] for i in idx])) for idx in nn_sis_data.integrand.indices]
closed_labels = ['\n'.join([feature_names[i] for i in np.where(idx)[0]]) for idx in nn_sis_data.integrand.indices]
ax.set_xticklabels(closed_labels,rotation=0)
ax.set_ylim([0,1]); ax.set_yticks([0,1])
ax.grid(False)
for spine in ['top','right','bottom']: ax.spines[spine].set_visible(False)
ax.legend(frameon=False,loc='lower center',bbox_to_anchor=(.5,-.2),ncol=3);
fig.tight_layout()
if root: fig.savefig(root+'nn_si.pdf')
nindices = len(nn_sis_data.integrand.indices)
fig,ax = pyplot.subplots(figsize=(9,5))
ticks = np.arange(nindices)
width = .25
for i,(alpha,species) in enumerate(zip([.25,.5,.75],data['target_names'])):
cvals = df_closed[species].to_numpy()
tvals = df_total[species].to_numpy()
ticks_i = ticks+i*width
ax.bar(ticks_i,cvals,width=width,align='edge',color='k',alpha=alpha,label=species)
#ax.bar(ticks_i,np.flip(tvals),width=width,align='edge',bottom=1-np.flip(tvals),color=color,alpha=.1)
ax.set_xlim([0,13+3*width])
ax.set_xticks(ticks+1.5*width)
# closed_labels = [r'$\underline{s}_{\{%s\}}$'%(','.join([r'\text{%s}'%feature_names[i] for i in idx])) for idx in nn_sis_data.integrand.indices]
closed_labels = ['\n'.join([feature_names[i] for i in np.where(idx)[0]]) for idx in nn_sis_data.integrand.indices]
ax.set_xticklabels(closed_labels,rotation=0)
ax.set_ylim([0,1]); ax.set_yticks([0,1])
ax.grid(False)
for spine in ['top','right','bottom']: ax.spines[spine].set_visible(False)
ax.legend(frameon=False,loc='lower center',bbox_to_anchor=(.5,-.2),ncol=3);
fig.tight_layout()
if root: fig.savefig(root+'nn_si.pdf')
