<jupyter_start><jupyter_text>Inkjet Figures<jupyter_code>from matplotlib.patches import ConnectionPatch
fig = plt.figure(figsize=(8, 7))

gs = fig.add_gridspec(3, 4*5)
gs.update(hspace=0.5,wspace=0.1)
dmin=-2.7
dlim=4
pal = ['#147edb','#c93b0c','#0b9e0e'] # EI, MPI, LCB
init_set=4 # grid spaces used for initialization round
optim_set=4 # grid spaces used for optimization round
inita=0
optima=0
rounds=['Initialization','Round 1','Round 2','Round 3', 'Round 4']
var=['Pressure','Frequency','Speed']
l=['p','f','s']

for u in range(5):
    if u > 0:
        inita=1
    if u > 1:
        optima+=1
    lower=init_set*inita+optima*optim_set
    upper=init_set+u*optim_set

    globals()['r'+str(u)]=bo.loc[bo['Round']==rounds[u],:]
    # position [0,u]
    if u == 0:
        globals()['p1'+str(u)]= fig.add_subplot(gs[0, lower:upper])
    else:
        globals()['p1'+str(u)]= fig.add_subplot(gs[0, lower:upper],sharey=p10)
    globals()['p1'+str(u)].axis('off')
    globals()['p2'+str(u)] = fig.add_subplot(gs[0, lower:upper-1],sharey=p10)
    globals()['p2'+str(u)].patch.set_alpha(0.0)
    globals()['p2'+str(u)].axis('off')
    s1=sns.kdeplot(data=globals()['r'+str(u)],y='Pressure',ax=globals()['p1'+str(u)],hue='Method',palette=pal,fill=True,alpha=0.8,clip=[0.,1.],multiple='stack',bw_adjust=0.5)
    s1.legend_.remove()
    globals()['p1'+str(u)].set_xlim([dmin,dlim])
    globals()['p2'+str(u)].set_xlim([-3.5,9])
    globals()['p1'+str(u)].plot((-0.1,-0.1),(0,1),c='k',linewidth=3)
    globals()['p2'+str(u)].text(-1.5,0,r'$\bf{0.0}$',fontsize=12)
    globals()['p2'+str(u)].text(-1.5,0.9,r'$\bf{1.0}$',fontsize=12)
    globals()['p2'+str(u)].text(-1.5,0.45,r'$\bf{0.5}$',fontsize=12)
    
    # position [1,u]
    if u == 0:
        globals()['f1'+str(u)]= fig.add_subplot(gs[1, lower:upper])
    else:
        globals()['f1'+str(u)]= fig.add_subplot(gs[1, lower:upper],sharey=f10)
    globals()['f1'+str(u)].axis('off')
    globals()['f2'+str(u)] = fig.add_subplot(gs[1, lower:upper-1],sharey=f10)
    globals()['f2'+str(u)].patch.set_alpha(0.0)
    globals()['f2'+str(u)].axis('off')
    s2=sns.kdeplot(data=globals()['r'+str(u)],y='Frequency',ax=globals()['f1'+str(u)],hue='Method',palette=pal,fill=True,alpha=0.8,clip=[0.,1.],multiple='stack',bw_adjust=0.5)
    if u > 0:
        s2.legend_.remove()
    else:
        s2.legend_.set_bbox_to_anchor((5.5, -0.8))
    globals()['f1'+str(u)].set_xlim([dmin,dlim])
    globals()['f2'+str(u)].set_xlim([-3.5,9])
    globals()['f1'+str(u)].plot((-0.1,-0.1),(0,1),c='k',linewidth=3)
    globals()['f2'+str(u)].text(-1.5,0,r'$\bf{0.0}$',fontsize=12)
    globals()['f2'+str(u)].text(-1.5,0.9,r'$\bf{1.0}$',fontsize=12)
    globals()['f2'+str(u)].text(-1.5,0.45,r'$\bf{0.5}$',fontsize=12)
    
    # position [2,u]
    if u == 0:
        globals()['s1'+str(u)]= fig.add_subplot(gs[2, lower:upper])
    else:
        globals()['s1'+str(u)]= fig.add_subplot(gs[2, lower:upper],sharey=s10)
    globals()['s1'+str(u)].axis('off')
    globals()['s2'+str(u)] = fig.add_subplot(gs[2, lower:upper-1],sharey=s10)
    globals()['s2'+str(u)].patch.set_alpha(0.0)
    globals()['s2'+str(u)].axis('off')
    s3=sns.kdeplot(data=globals()['r'+str(u)],y='Speed',ax=globals()['s1'+str(u)],hue='Method',palette=pal,fill=True,alpha=0.8,clip=[0.,1.],multiple='stack',bw_adjust=0.5)
    s3.legend_.remove()
    globals()['s1'+str(u)].set_xlim([dmin,dlim])
    globals()['s2'+str(u)].set_xlim([-3.5,9])
    globals()['s1'+str(u)].plot((-0.1,-0.1),(0,1),c='k',linewidth=3)
    globals()['s2'+str(u)].text(-1.5,0,r'$\bf{0.0}$',fontsize=12)
    globals()['s2'+str(u)].text(-1.5,0.9,r'$\bf{1.0}$',fontsize=12)
    globals()['s2'+str(u)].text(-1.5,0.45,r'$\bf{0.5}$',fontsize=12)
# plot text
p20.text(-6,0.20,r'$\bf{Pressure}$',fontsize=12,rotation=90)
f20.text(-6,0.15,r'$\bf{Frequency}$',fontsize=12,rotation=90)
s20.text(-6,0.32,r'$\bf{Speed}$',fontsize=12,rotation=90)
s20.text(-6,-0.35,r'$\bf{Initialization}$',fontsize=12)
s20.text(16,-0.35,r'$\bf{Batch\ 1}$',fontsize=12)
s20.text(33,-0.35,r'$\bf{Batch\ 2}$',fontsize=12)
s20.text(50,-0.35,r'$\bf{Batch\ 3}$',fontsize=12)
s20.text(67,-0.35,r'$\bf{Batch\ 4}$',fontsize=12)
# plt.savefig(impath+'kde_flow_R3.png',dpi=300)
cv_x_tem=bo.loc[:,['Pressure','Frequency','Speed','Method']]
cv_y_tem=bo.loc[:,['Min Function Val']]
meth='EI'

cv_y_tem=cv_y_tem[cv_x_tem['Method']==meth]
cv_x_tem=cv_x_tem[cv_x_tem['Method']==meth]
cv_x_tem.reset_index(drop=True,inplace=True)
cv_y_tem.reset_index(drop=True,inplace=True)
if meth=='LCB':
    drop_list=[10,15]# list of rows to drop
    cv_x_tem=cv_x_tem.drop(drop_list)
    cv_y_tem=cv_y_tem.drop(drop_list)
elif meth=='EI':
    drop_list=[15]# list of rows to drop
    cv_x_tem=cv_x_tem.drop(drop_list)
    cv_y_tem=cv_y_tem.drop(drop_list)
elif meth=='MPI':
    drop_list=[15]# list of rows to drop
    cv_x_tem=cv_x_tem.drop(drop_list)
    cv_y_tem=cv_y_tem.drop(drop_list)

res=100
x1=np.linspace(0,1,res)
x2=np.linspace(0,1,res)
X1, X2 = np.meshgrid(x1,x2)

feas=0.75 # score threshold for feasibility
for n in range(3):
    cv_x_temp0=np.copy(cv_x_tem.iloc[:,n])
    cv_x_temp1=np.copy(cv_x_tem.iloc[:,(n+1)%3])
    cv_y_temp=np.copy(cv_y_tem.iloc[:,0])
    for i in range(cv_x_tem.shape[0]):
        for j in range(cv_x_tem.shape[0]):
            if cv_x_temp0[i]==cv_x_temp0[j] and cv_x_temp1[i]==cv_x_temp1[j] and i!=j:
                if cv_y_temp[i] < cv_y_temp[j]: # want to remove failed prints over good prints
                    cv_x_temp0[j]=np.nan
                    cv_x_temp1[j]=np.nan
                    cv_y_temp[j]=np.nan
                else:
                    cv_x_temp0[i]=np.nan
                    cv_x_temp1[i]=np.nan
                    cv_y_temp[i]=np.nan
    cv_x_temp0=cv_x_temp0[~np.isnan(cv_x_temp0)]
    cv_x_temp1=cv_x_temp1[~np.isnan(cv_x_temp1)]
    cv_y_temp=cv_y_temp[~np.isnan(cv_y_temp)]
    df_rbf=scipy.interpolate.Rbf(cv_x_temp0,cv_x_temp1,cv_y_temp,function='linear')
    globals()['surf_f'+str(n)]=df_rbf(X1,X2)
    globals()['surf_inf'+str(n)]=np.copy(globals()['surf_f'+str(n)])
    globals()['surf_inf'+str(n)][globals()['surf_inf'+str(n)]<feas]=-1 # feasible
    globals()['surf_inf'+str(n)][globals()['surf_inf'+str(n)]>=feas]=1 # infeasible
    globals()['surf_f'+str(n)][globals()['surf_f'+str(n)]>=feas]=np.nan # remove infeasible points from plot
fs=16
fs1=15
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(6,6))
ax.set_xticks([])
ax.set_yticks([])
ax.set_frame_on(False)
# add axes
ax = fig.add_subplot(111,projection='3d')
ax.plot([1,1],[0,1],[1,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)
ax.plot([1,0],[0,0],[1,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)
ax.plot([1,1],[0,0],[0,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)

lev=[-1,-0.9,1]
colr='#a3fff9'
ax.contourf(X1, X2,surf_inf0, lev,alpha=1,offset=1,zdir='z',colors=(colr,'k'),zorder=0)
ax.contourf(surf_inf1, X1, X2, lev,alpha=1,offset=1,zdir='x',colors=(colr,'k'),zorder=0)
ax.contourf(X2, surf_inf2, X1, lev,alpha=1,offset=0,zdir='y',colors=(colr,'k'),zorder=0)
ax.grid(False)
ax.w_xaxis.pane.fill = False
ax.w_yaxis.pane.fill = False
ax.w_zaxis.pane.fill = False
ax.view_init(azim=45)
ax.set_xlabel(r'$\bf{Pressure}$'+'\n[Normalized]',fontsize=fs)
ax.set_ylabel(r'$\bf{Frequency}$'+'\n[Normalized]',fontsize=fs)
ax.zaxis.set_rotate_label(False)
ax.set_zlabel(r'$\bf{Translation}$'+'\n[Normalized]',rotation=90,fontsize=fs)
ax.set_xticks([0,1])
ax.set_xlim([0,1])
ax.set_yticks([0,1])
ax.set_ylim([1,0])
ax.set_zticks([0,1])
ax.set_zlim([0,1])
ax.tick_params(axis='z',labelsize=fs1,pad=0)
ax.tick_params(axis='x',labelsize=fs1,pad=-1)
ax.tick_params(axis='y',labelsize=fs1,pad=-1)
ax.zaxis.labelpad = -1
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(False)
ax.dist=12

legend_elements = [Patch(facecolor='#a3fff9',label='Feasible'),
                  Patch(facecolor='k',label='Infeasible')]
ax.legend(handles=legend_elements, bbox_to_anchor=(0.85, 0.),prop={'size': fs},ncol=2)
plt.tight_layout()
plt.show()

print('Feasible % space for '+meth+' = '+str(100*(np.sum(surf_inf0==-1)+np.sum(surf_inf1==-1)+np.sum(surf_inf2==-1))/(np.sum(np.abs(surf_inf0))+np.sum(np.abs(surf_inf1))+np.sum(np.abs(surf_inf2)))))
fs=16
fs1=15
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(6,5))
ax.set_xticks([])
ax.set_yticks([])
ax.set_frame_on(False)
# add axes
ax = fig.add_subplot(111,projection='3d')
ax.plot([1,1],[0,1],[1,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)
ax.plot([1,0],[0,0],[1,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)
ax.plot([1,1],[0,0],[0,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)

vmin=0.4
vmax=feas
lev=np.linspace(vmin,vmax,100)
clev=np.linspace(vmin,vmax,3)
pl=ax.contourf(X1, X2,surf_f0, lev,alpha=1,offset=1,zdir='z',zorder=0,cmap='magma_r',vmin=vmin,vmax=vmax, antialiased=False)
ax.contour(X1, X2,surf_f0, clev,offset=1,zdir='z',colors='w',linewidths=1)
ax.contourf(surf_f1, X1, X2, lev,alpha=1,offset=1,zdir='x',zorder=0,cmap='magma_r',vmin=vmin,vmax=vmax, antialiased=False)
ax.contour(surf_f1, X1,X2,clev, offset=1,zdir='x',colors='w',linewidths=1)
ax.contourf(X2, surf_f2, X1, lev,alpha=1,offset=0,zdir='y',zorder=0,cmap='magma_r',vmin=vmin,vmax=vmax, antialiased=False)
ax.contour(X2, surf_f2,X1,clev, offset=0,zdir='y',colors='w',linewidths=1)
ax.grid(False)
ax.w_xaxis.pane.fill = False
ax.w_yaxis.pane.fill = False
ax.w_zaxis.pane.fill = False
ax.view_init(azim=45)
ax.set_xlabel(r'$\bf{Pressure}$'+'\n[Normalized]',fontsize=fs)
ax.set_ylabel(r'$\bf{Frequency}$'+'\n[Normalized]',fontsize=fs)
ax.zaxis.set_rotate_label(False)
ax.set_zlabel(r'$\bf{Translation}$'+'\n[Normalized]',rotation=90,fontsize=fs)
ax.set_xticks([0,1])
ax.set_xlim([0,1])
ax.set_yticks([0,1])
ax.set_ylim([1,0])
ax.set_zticks([0,1])
ax.set_zlim([0,1])
ax.tick_params(axis='z',labelsize=fs1,pad=0)
ax.tick_params(axis='x',labelsize=fs1,pad=-1)
ax.tick_params(axis='y',labelsize=fs1,pad=-1)
ax.zaxis.labelpad = -1
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(False)
ax.dist=12
plt.tight_layout()
plt.show()
fs=16
fs1=15
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(6,5))
ax.set_xticks([])
ax.set_yticks([])
ax.set_frame_on(False)
# add axes
ax = fig.add_subplot(111,projection='3d')
ax.plot([1,1],[0,1],[1,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)
ax.plot([1,0],[0,0],[1,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)
ax.plot([1,1],[0,0],[0,1],color='grey',zorder=1000,linestyle='--',alpha=0.8,linewidth=2)

vmin=0.35
vmax=feas
lev=np.linspace(vmin,vmax,100)
clev=np.linspace(vmin,vmax,3)
ax.scatter(cv_x_tem['Pressure'],cv_x_tem['Frequency'],cv_x_tem['Speed'],s=50)
ax.grid(False)
ax.w_xaxis.pane.fill = False
ax.w_yaxis.pane.fill = False
ax.w_zaxis.pane.fill = False
ax.view_init(azim=45)
ax.set_xlabel(r'$\bf{Pressure}$'+'\n[Normalized]',fontsize=fs)
ax.set_ylabel(r'$\bf{Frequency}$'+'\n[Normalized]',fontsize=fs)
ax.zaxis.set_rotate_label(False)
ax.set_zlabel(r'$\bf{Translation}$'+'\n[Normalized]',rotation=90,fontsize=fs)
ax.set_xticks([0,1])
ax.set_xlim([0,1])
ax.set_yticks([0,1])
ax.set_ylim([1,0])
ax.set_zticks([0,1])
ax.set_zlim([0,1])
ax.tick_params(axis='z',labelsize=fs1,pad=0)
ax.tick_params(axis='x',labelsize=fs1,pad=-1)
ax.tick_params(axis='y',labelsize=fs1,pad=-1)
ax.zaxis.labelpad = -1
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(False)
ax.dist=12

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Microfluidics Figures<jupyter_code>fig = plt.figure(figsize=(9.5, 5.5))

gs = fig.add_gridspec(2, 4*5)
gs.update(hspace=0.5,wspace=0.1)
dmin=-2.7
dlim=4
pal = ['#147edb','#c93b0c','#0b9e0e'] # EI, MPI, LCB
init_set=4 # grid spaces used for initialization round
optim_set=4 # grid spaces used for optimization round
inita=0
optima=0
rounds=['Initialization','Round 1','Round 2','Round 3', 'Round 4']
var=['Water Pressure','Oil Pressure']
l=['p','f']

for u in range(5):
    if u > 0:
        inita=1
    if u > 1:
        optima+=1
    lower=init_set*inita+optima*optim_set
    upper=init_set+u*optim_set

    globals()['r'+str(u)]=bom.loc[bo['Round']==rounds[u],:]
    # position [0,u]
    if u == 0:
        globals()['p1'+str(u)]= fig.add_subplot(gs[0, lower:upper])
    else:
        globals()['p1'+str(u)]= fig.add_subplot(gs[0, lower:upper],sharey=p10)
    globals()['p1'+str(u)].axis('off')
    globals()['p2'+str(u)] = fig.add_subplot(gs[0, lower:upper-1],sharey=p10)
    globals()['p2'+str(u)].patch.set_alpha(0.0)
    globals()['p2'+str(u)].axis('off')
    s1=sns.kdeplot(data=globals()['r'+str(u)],y='WaterPressure',ax=globals()['p1'+str(u)],hue='Method',palette=pal,fill=True,alpha=0.8,clip=[0.,1.],multiple='stack',bw_adjust=0.5)
    s1.legend_.remove()
    globals()['p1'+str(u)].set_xlim([dmin,dlim])
    globals()['p2'+str(u)].set_xlim([-3.5,9])
    globals()['p1'+str(u)].plot((-0.1,-0.1),(0,1),c='k',linewidth=3)
    globals()['p2'+str(u)].text(-1.5,0,r'$\bf{0.0}$',fontsize=12)
    globals()['p2'+str(u)].text(-1.5,0.9,r'$\bf{1.0}$',fontsize=12)
    globals()['p2'+str(u)].text(-1.5,0.45,r'$\bf{0.5}$',fontsize=12)
    
    # position [1,u]
    if u == 0:
        globals()['f1'+str(u)]= fig.add_subplot(gs[1, lower:upper])
    else:
        globals()['f1'+str(u)]= fig.add_subplot(gs[1, lower:upper],sharey=f10)
    globals()['f1'+str(u)].axis('off')
    globals()['f2'+str(u)] = fig.add_subplot(gs[1, lower:upper-1],sharey=f10)
    globals()['f2'+str(u)].patch.set_alpha(0.0)
    globals()['f2'+str(u)].axis('off')
    s2=sns.kdeplot(data=globals()['r'+str(u)],y='OilPressure',ax=globals()['f1'+str(u)],hue='Method',palette=pal,fill=True,alpha=0.8,clip=[0.,1.],multiple='stack',bw_adjust=0.5)
    if u > 0:
        s2.legend_.remove()
    else:
        s2.legend_.set_bbox_to_anchor((5.5, -0.8))
    globals()['f1'+str(u)].set_xlim([dmin,dlim])
    globals()['f2'+str(u)].set_xlim([-3.5,9])
    globals()['f1'+str(u)].plot((-0.1,-0.1),(0,1),c='k',linewidth=3)
    globals()['f2'+str(u)].text(-1.5,0,r'$\bf{0.0}$',fontsize=12)
    globals()['f2'+str(u)].text(-1.5,0.9,r'$\bf{1.0}$',fontsize=12)
    globals()['f2'+str(u)].text(-1.5,0.45,r'$\bf{0.5}$',fontsize=12)

# plot text
p20.text(-6,0.20,r'$\bf{Water\ Pressure}$',fontsize=12,rotation=90)
f20.text(-6,0.15,r'$\bf{Oil\ Pressure}$',fontsize=12,rotation=90)
s20.text(-6,-0.35,r'$\bf{Initialization}$',fontsize=12)
s20.text(16,-0.35,r'$\bf{Batch\ 1}$',fontsize=12)
s20.text(33,-0.35,r'$\bf{Batch\ 2}$',fontsize=12)
s20.text(50,-0.35,r'$\bf{Batch\ 3}$',fontsize=12)
s20.text(67,-0.35,r'$\bf{Batch\ 4}$',fontsize=12)
cv_x_tem=bom.loc[:,['WaterPressure','OilPressure','Method']]
cv_y_tem=bom.loc[:,['Min Function Val']]
meth='EI'
cv_y_tem=cv_y_tem[cv_x_tem['Method']==meth]
cv_x_tem=cv_x_tem[cv_x_tem['Method']==meth]

res=100
x1=np.linspace(0,1,res)
x2=np.linspace(0,1,res)
X1, X2 = np.meshgrid(x1,x2)

feas=0.75 # score threshold for feasibility
for n in range(1):
    cv_x_temp0=np.copy(cv_x_tem.iloc[:,n])
    cv_x_temp1=np.copy(cv_x_tem.iloc[:,(n+1)%2])
    cv_y_temp=np.copy(cv_y_tem.iloc[:,0])
    for i in range(cv_x_tem.shape[0]):
        for j in range(cv_x_tem.shape[0]):
            if cv_x_temp0[i]==cv_x_temp0[j] and cv_x_temp1[i]==cv_x_temp1[j] and i!=j:
                if cv_x_temp0[i] < cv_x_temp0[j]: # want to remove failed prints over good prints
                    cv_x_temp0[j]=np.nan
                    cv_x_temp1[j]=np.nan
                    cv_y_temp[j]=np.nan
                else:
                    cv_x_temp0[i]=np.nan
                    cv_x_temp1[i]=np.nan
                    cv_y_temp[i]=np.nan
    cv_x_temp0=cv_x_temp0[~np.isnan(cv_x_temp0)]
    cv_x_temp1=cv_x_temp1[~np.isnan(cv_x_temp1)]
    cv_y_temp=cv_y_temp[~np.isnan(cv_y_temp)]
    df_rbf=scipy.interpolate.Rbf(cv_x_temp0,cv_x_temp1,cv_y_temp,function='linear')
    globals()['surf_f'+str(n)]=df_rbf(X1,X2)
    globals()['surf_inf'+str(n)]=np.copy(globals()['surf_f'+str(n)])
    globals()['surf_inf'+str(n)][globals()['surf_inf'+str(n)]<feas]=-1 # feasible
    globals()['surf_inf'+str(n)][globals()['surf_inf'+str(n)]>=feas]=1 # infeasible
    globals()['surf_f'+str(n)][globals()['surf_f'+str(n)]>=feas]=np.nan # remove infeasible points from plot
fs=21
fs1=22
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(6,6))
ax.set_xticks([])
ax.set_yticks([])
ax.set_frame_on(False)
lev=[-1,-0.9,1]
colr='#a3fff9'
ax.contourf(X1, X2,surf_inf0, lev,alpha=1,colors=(colr,'k'),zorder=0)
ax.grid(False)
ax.set_xlabel(r'$\bf{Water\ Pressure}$'+'\n[Normalized]',fontsize=fs)
ax.set_ylabel(r'$\bf{Oil\ Pressure}$'+'\n[Normalized]',fontsize=fs)

ax.set_xticks([0,1])
ax.set_xlim([0,1])
ax.set_yticks([0,1])
ax.set_ylim([0,1])
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1,pad=5)
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(False)
ax.dist=12
plt.tight_layout()

# plt.savefig(impath+'feasibility_microf_LCB_R1.png',dpi=300,bbox_inches='tight',transparent=True)
plt.show()

print('Feasible % space for '+meth+' = '+str(100*(np.sum(surf_inf0==-1))/(np.sum(np.abs(surf_inf0)))))
fs=21
fs1=22
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(6,6))
ax.set_xticks([])
ax.set_yticks([])
ax.set_frame_on(False)
vmin=0.1
vmax=feas
lev=np.linspace(vmin,vmax,100)
clev=np.linspace(vmin,0.5,3)
pl=ax.contourf(X1, X2,surf_f0, lev,alpha=1,zorder=0,cmap='magma_r',vmin=vmin,vmax=vmax, antialiased=False)
ax.contour(X1, X2,surf_f0, clev,colors='w',linewidths=1)

ax.grid(False)
ax.set_xlabel(r'$\bf{Water\ Pressure}$'+'\n[Normalized]',fontsize=fs)
ax.set_ylabel(r'$\bf{Oil\ Pressure}$'+'\n[Normalized]',fontsize=fs)

ax.set_xticks([0,1])
ax.set_xlim([0,1])
ax.set_yticks([0,1])
ax.set_ylim([0,1])
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1,pad=5)
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(True)
ax.dist=12
plt.tight_layout()
# plt.savefig(impath+'learned_microf_EI_R1.png',dpi=300,bbox_inches='tight',transparent=True)
plt.show()
fs=21
fs1=22
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(6,6))
ax.set_xticks([])
ax.set_yticks([])
ax.set_frame_on(False)
ax.scatter(cv_x_tem['WaterPressure'],cv_x_tem['OilPressure'],s=60,alpha=0.7)

ax.grid(False)
ax.set_xlabel(r'$\bf{Water\ Pressure}$'+'\n[Normalized]',fontsize=fs)
ax.set_ylabel(r'$\bf{Oil\ Pressure}$'+'\n[Normalized]',fontsize=fs)

ax.set_xticks([0,1])
ax.set_xlim([0,1])
ax.set_yticks([0,1])
ax.set_ylim([0,1])
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1,pad=5)
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(True)
ax.dist=12
plt.tight_layout()
# plt.savefig(impath+'sample_microf_LCB_R1.png',dpi=300,bbox_inches='tight',transparent=True)
plt.show()
fs=21
fs1=22
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(8,6))
ax.set_frame_on(False)
pl=plt.scatter(cawe_geom['We*10^-3'],cawe_geom['Ca*10^-2'],c=cawe_geom['avg_diam'],cmap=parula_map,s=200,edgecolor='k',alpha=0.8) # swap axes because data is labelled incorrectly

ax.grid(False)
ax.set_xlabel('Capillary Number (Ca) $10^{-2}$',fontsize=fs)
ax.set_ylabel('Weber Number (We) $10^{-3}$',fontsize=fs)
ax.set_yticks([0.5,1.,1.5,2.,2.5])
ax.set_xticks([1.,2.,3.,4.])
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1,pad=5)
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(True)
ax.dist=12

cbar=plt.colorbar(pl,shrink=0.7,pad=0.1)
cbar.ax.tick_params(labelsize=fs) 
cbar.set_label('Droplet Diameter [$\mu$m]',size=fs)
plt.tight_layout()
plt.grid()
# plt.savefig(impath+'micro_diam_R1.png',dpi=300,bbox_inches='tight',transparent=True)
plt.show()
fs=21
fs1=22
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(8,6))
ax.set_frame_on(False)
pl=plt.scatter(cawe_geom['We*10^-3'][103-97:],cawe_geom['Ca*10^-2'][103-97:],c=drops_per_row['drop_avg'],cmap='cividis',s=200,edgecolor='k',alpha=0.8) # swap axes because data is labelled incorrectly

ax.grid(False)
ax.set_xlabel('Capillary Number (Ca) $10^{-2}$',fontsize=fs)
ax.set_ylabel('Weber Number (We) $10^{-3}$',fontsize=fs)
ax.set_yticks([0.5,1.,1.5,2.,2.5])
ax.set_xticks([1.,2.,3.,4.])
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1,pad=5)
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(True)
ax.dist=12

cbar=plt.colorbar(pl,shrink=0.7,pad=0.1)
cbar.ax.tick_params(labelsize=fs) 
cbar.set_label('# of Drops per Row',size=fs)
plt.tight_layout()
plt.grid()
# plt.savefig(impath+'micro_drops_per_row_R1.png',dpi=300,bbox_inches='tight',transparent=True)
plt.show()
fs=21
fs1=22
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(8,6))
ax.set_frame_on(False)
pl=plt.scatter(cawe_geom['We*10^-3'],cawe_geom['Ca*10^-2'],c=cawe_geom['objective'],cmap='magma_r',s=200,edgecolor='k',alpha=0.8) # swap axes because data is labelled incorrectly

ax.grid(False)
ax.set_xlabel('Capillary Number (Ca) $10^{-2}$',fontsize=fs)
ax.set_ylabel('Weber Number (We) $10^{-3}$',fontsize=fs)
ax.set_yticks([0.5,1.,1.5,2.,2.5])
ax.set_xticks([1.,2.,3.,4.])
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1,pad=5)
ax.xaxis.labelpad = 2
ax.yaxis.labelpad = 2
ax.set_frame_on(True)
ax.dist=12

cbar=plt.colorbar(pl,shrink=0.7,pad=0.1)
cbar.ax.tick_params(labelsize=fs) 
cbar.set_label('Objective Value',size=fs)
plt.tight_layout()
plt.grid()
# plt.savefig(impath+'CaWe_numbers_R2.png',dpi=300,bbox_inches='tight',transparent=True)
plt.show()<jupyter_output><empty_output><jupyter_text>Minimum Objective Value<jupyter_code>fs=21
fs1=19
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(7,6))

x_plot=np.arange(20,61,1)
plt.plot(x_plot,boEI_min[20:],label='GP-EI',c=pal[0],linewidth=2,alpha=0.8)
plt.plot(x_plot,boMPI_min[20:],label='GP-MPI',c=pal[1],linewidth=2,alpha=0.8)
plt.plot(x_plot,boLCB_min[20:],label='GP-LCB',c=pal[2],linewidth=2,alpha=0.8)
plt.legend(prop={'size':fs1})
plt.ylabel('Min Objective Value, $\ell^*$',fontsize=fs)
plt.xlabel('Sample',fontsize=fs)
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1-1,pad=5)
plt.ylim([0.,0.6])
# plt.savefig(impath+'inkjet_minval_R2.png',dpi=300,bbox_inches='tight')
plt.show()

fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(7,6))

x_plot=np.arange(20,61,1)
plt.plot(x_plot,bomEI_min[20:],label='GP-EI',c=pal[0],linewidth=2,alpha=0.8)
plt.plot(x_plot,bomMPI_min[20:],label='GP-MPI',c=pal[1],linewidth=2,alpha=0.8)
plt.plot(x_plot,bomLCB_min[20:],label='GP-LCB',c=pal[2],linewidth=2,alpha=0.8)
plt.legend(prop={'size':fs1})
plt.ylabel('Min Objective Value, $\ell^*$',fontsize=fs)
plt.xlabel('Sample',fontsize=fs)
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1,pad=5)
plt.ylim([0,0.3])
# plt.savefig(impath+'micro_minval_R2.png',dpi=300,bbox_inches='tight')
plt.show()

# compute interpolation field for all 3 runs
res=100
x1=np.linspace(0,1,res)
x2=np.linspace(0,1,res)
X1, X2 = np.meshgrid(x1,x2)

cv_x_tem=bo.loc[:,['Pressure','Frequency','Speed','Method']]
cv_y_tem=bo.loc[:,['Min Function Val']]
for n in range(3):
    cv_x_temp0=np.copy(cv_x_tem.iloc[:,n])
    cv_x_temp1=np.copy(cv_x_tem.iloc[:,(n+1)%3])
    cv_y_temp=np.copy(cv_y_tem.iloc[:,0])
    for i in range(cv_x_tem.shape[0]):
        for j in range(cv_x_tem.shape[0]):
            if cv_x_temp0[i]==cv_x_temp0[j] and cv_x_temp1[i]==cv_x_temp1[j] and i!=j:
                if cv_y_temp[i] < cv_y_temp[j]: # want to remove failed prints over good prints
                    cv_x_temp0[j]=np.nan
                    cv_x_temp1[j]=np.nan
                    cv_y_temp[j]=np.nan
                else:
                    cv_x_temp0[i]=np.nan
                    cv_x_temp1[i]=np.nan
                    cv_y_temp[i]=np.nan
    cv_x_temp0=cv_x_temp0[~np.isnan(cv_x_temp0)]
    cv_x_temp1=cv_x_temp1[~np.isnan(cv_x_temp1)]
    cv_y_temp=cv_y_temp[~np.isnan(cv_y_temp)]
    df_rbf=scipy.interpolate.Rbf(cv_x_temp0,cv_x_temp1,cv_y_temp,function='linear')
    globals()['surf_f'+str(n)]=df_rbf(X1,X2)
    
# run random sampling simulation
import random
min_ys_100 = []
for _ in range(3): # take 100 averages
    rand_xs = []
    rand_ys = []
    min_ys = [1]
    for r in range(42):
        rand_x = random.sample(range(0, 100), 3)
        rand_xs.append(rand_x)
        y0 = surf_f0[rand_x[0], rand_x[1]]
        y1 = surf_f1[rand_x[1], rand_x[2]]
        y2 = surf_f2[rand_x[2], rand_x[0]]
        Y = np.min([y0, y1, y2])
        rand_ys.append(Y)
        if r > 0:
            if Y < min_ys[r-1]:
                min_ys.append(Y)
            else:
                min_ys.append(min_ys[r-1])
        min_ys_100.append(min_ys)
        
# get average and bounds of random sampling
min_ys_100_avg = np.mean(np.array(min_ys_100), axis=0)[1:]
min_ys_lower = np.min(np.array(min_ys_100), axis=0)[1:]
min_ys_upper = np.max(np.array(min_ys_100), axis=0)[1:]

# get average and bounds of 3 experimental bo runs
min_bos = np.array([boEI_min[20:],boMPI_min[20:],boLCB_min[20:]])
min_bos_lower = np.min(min_bos, axis=0)
min_bos_upper = np.max(min_bos, axis=0)
min_bos_avg = np.mean(min_bos, axis=0)
fs=21
fs1=19
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(7,6))

x_plot=np.arange(20,61,1)
plt.fill_between(x_plot, min_ys_lower, min_ys_upper, color='k', alpha=0.1)
plt.plot(x_plot,min_ys_100_avg,label='3-trial avg. Random Search',c='k',linewidth=2,alpha=0.8, marker='|')

plt.fill_between(x_plot, min_bos_lower, min_bos_upper, color='#9467bd', alpha=0.1)
plt.plot(x_plot,min_bos_avg,label='3-trial avg. Bayesian Optimization',c='#9467bd',linewidth=2,alpha=0.8, marker='|')

plt.legend(prop={'size':fs1}, bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left',
           ncol=1, mode="expand", borderaxespad=0.)
plt.ylabel('Min Objective, $f_{min}$',fontsize=fs)
plt.xlabel('Sample',fontsize=fs)
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1-1,pad=5)
plt.ylim([0.3,0.85])
# plt.savefig(impath+'drops-bo-vs-rand1.png',dpi=300,bbox_inches='tight')
plt.show()

fs=21
fs1=19
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(7,6))

x_plot=np.arange(20,61,1)
plt.fill_between(x_plot, min_ys_lower, min_ys_upper, color='k', alpha=0.1)
plt.plot(x_plot,min_ys_100_avg,label='3-trial avg. Random Search',c='k',linewidth=2,alpha=0.8, marker='|')

plt.fill_between(x_plot, min_bos_lower, min_bos_upper, color='#9467bd', alpha=0.1)
plt.plot(x_plot,min_bos_avg,label='3-trial avg. Bayesian Optimization',c='#9467bd',linewidth=2,alpha=0.8, marker='|')

plt.legend(prop={'size':fs1}, bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left',
           ncol=1, mode="expand", borderaxespad=0.)
plt.ylabel('Min Objective, $f_{min}$',fontsize=fs)
plt.xlabel('Sample',fontsize=fs)
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1-1,pad=5)
plt.ylim([0.3,0.85])
# plt.savefig(impath+'drops-bo-vs-rand1.png',dpi=300,bbox_inches='tight')
plt.show()

print(min_ys_100_avg, '\n\n',min_bos_avg)
fs=21
fs1=19
fig,ax = plt.subplots(ncols=1,nrows=1,figsize=(7,6))

x_plot=np.arange(20,61,1)
plt.plot(x_plot,boEI_min[20:],label='EI',c=pal[0],linewidth=2,alpha=0.8)
plt.plot(x_plot,boMPI_min[20:],label='MPI',c=pal[1],linewidth=2,alpha=0.8)
plt.plot(x_plot,boLCB_min[20:],label='LCB',c=pal[2],linewidth=2,alpha=0.8)
plt.legend(prop={'size':fs1}, bbox_to_anchor=(0., 1.02, 0.75, .102), loc='lower left',
           ncol=3, mode="expand", borderaxespad=0.)
plt.ylabel('Minimization Objective, $\ell^*$',fontsize=fs)
plt.xlabel('Sample',fontsize=fs)
ax.tick_params(axis='x',labelsize=fs1,pad=5)
ax.tick_params(axis='y',labelsize=fs1-1,pad=5)
plt.ylim([0.35,0.85])
# plt.savefig(impath+'drops-zoom-legend.png',dpi=300,bbox_inches='tight')
plt.show()<jupyter_output><empty_output><jupyter_text>SHAP feature importance<jupyter_code>from sklearn.model_selection import train_test_split
from sklearn import preprocessing
from sklearn.ensemble import RandomForestRegressor
import shap
bo_notinit=bo[bo['Round']!='Initialization']
# The target variable is 'quality'.
Y = bo_notinit['Min Function Val']
X =  bo_notinit[['Pressure','Frequency','Speed']]
# Split the data into train and test data:
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2)
# Build the model with the random forest regression algorithm:
model = RandomForestRegressor(max_depth=6, random_state=0, n_estimators=10)
model.fit(X_train, Y_train)
shap_values = shap.TreeExplainer(model).shap_values(X_train)

fig = plt.figure()
shap.summary_plot(shap_values, X_train)
# fig.savefig(impath+'shap_vals.png',dpi=300,bbox_inches='tight')
bom_notinit=bom[bom['Round']!='Initialization']
# The target variable is 'quality'.
Y = bom_notinit['Min Function Val']
X =  bom_notinit[['WaterPressure','OilPressure']]
# Split the data into train and test data:
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2)
# Build the model with the random forest regression algorithm:
model = RandomForestRegressor(max_depth=6, random_state=0, n_estimators=10)
model.fit(X_train, Y_train)
shap_values = shap.TreeExplainer(model).shap_values(X_train)

fig = plt.figure()
shap.summary_plot(shap_values, X_train)
# fig.savefig(impath+'shap_micro_vals.png',dpi=300,bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22702943800178413, 'fraction_numerical': 0.0577094627050024, 'mean_word_length': 6.356638061585058, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 28, 'https://': 0, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '196994', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 14290, 'n_tokens_neox': 13595, 'n_words': 1298, 'path': 'Plot_Figures.ipynb', 'repo_name': 'PV-Lab/Automated-Fluid-Optimization', 'size': '29146'}