<jupyter_start><jupyter_text>Comparision of  MC estimate of statistics using Fipy(ground truth) vs statistics by Deep Neural Network with 'n samples' each<jupyter_code>#load libraries 
from __future__ import division
import argparse

import matplotlib
matplotlib.use('PS')
%matplotlib inline
import tensorflow as tf
import random 
import numpy as np
import os
os.environ['KERAS_BACKEND'] = 'tensorflow'
from keras.models import Model
from keras.layers import Dense, Activation, Input, concatenate, Lambda, Add
from keras.utils import plot_model
from keras import backend as K
from keras.utils.generic_utils import get_custom_objects

import matplotlib.pyplot as plt 
import GPy
from fipy import *
# import matplotlib as mpl
# mpl.rcParams['figure.dpi']= 300
import seaborn as sns 
sns.set_context('paper')
sns.set_style('white')
from pdb import set_trace as keyboard
import sys
import time
# plt.rcParams["patch.force_edgecolor"] = True
# maindir='/Users/sharmila/Desktop/research/work/03_NN_PDE/Level3/part8_image/2_2dim_c=0_f=0'
# Taking needed inputs and loading files

# Load files
# load cellcenters data
cellcenters = np.load('cellcenters_nx1=32_nx2=32.npy')
# load MC input samples [conductivity] # these inputs are already bounded.
MC_inputfield = np.load('MC_samples_inputfield_exp_nx1=32_nx2=32_lx1=0.06_lx2=0.15_v=0.75_num_samples=100000.npy')
# load MC samples fipy # the below fipy file is for bounded inputs itself.
MC_u_fipy = np.load('MC_samples_u_fipy_exp_nx1=32_nx2=32_lx1=0.06_lx2=0.15_v=0.75_num_samples=100000.npy')

# DNN surrogate inputs
############################
DNN_type = 'Resnet'
act_func = 'swish'
n = 500 # Number of neurons in each block
num_block = 3 # Number of blocks
############################
DNN_weights_path='my_model_weights.h5'
print cellcenters.shape
print MC_inputfield.shape
print MC_u_fipy.shape
print np.shape(cellcenters)
# print cellcenters

nx1=int(np.sqrt(cellcenters.shape[0]))
nx2=int(np.sqrt(cellcenters.shape[0]))
print nx1
print nx2

x1_f = cellcenters[:,0][:,None] # x1-coordinates of cell centers (nx1*nx2,1) matrix
x2_f = cellcenters[:,1][:,None] # x2-coordinates of cell centers (nx1*nx2,1) matrix
print np.shape(x1_f)
print np.shape(x2_f)
# print x1_f
# print x2_f
cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", ["navy","blue","deepskyblue","limegreen","yellow","darkorange","red","maroon"])<jupyter_output><empty_output><jupyter_text>A.ground_using_fipy<jupyter_code>print MC_u_fipy.shape
print MC_u_fipy

ground_mean_fipy=np.mean(MC_u_fipy,axis=0)
ground_var_fipy=np.var(MC_u_fipy,axis=0)
print ground_mean_fipy.shape
print ground_var_fipy.shape
fig, ax = plt.subplots(figsize=(10,6))
c = ax.contourf(x1_f.reshape((nx1,nx2)), x2_f.reshape((nx1,nx2)), ground_mean_fipy.reshape((nx1,nx2)),100,cmap=cmap)   # set levels as previous levels
# This is the fix for the white lines between contour levels (https://stackoverflow.com/questions/8263769/hide-contour-linestroke-on-pyplot-contourf-to-get-only-fills)
for j in c.collections:
    j.set_edgecolor("face")
plt.colorbar(c)
# plt.title('mean: fipy')
plt.xlabel('$x_1$', fontsize = 14)
plt.ylabel('$x_2$', fontsize = 14)
plt.legend(loc='best')
plt.tight_layout()
plt.savefig(os.path.join(os.getcwd(),'mean_fipy.pdf'))
plt.show()

fig, ax = plt.subplots(figsize=(10,6))
c = ax.contourf(x1_f.reshape((nx1,nx2)), x2_f.reshape((nx1,nx2)), ground_var_fipy.reshape((nx1,nx2)),100,cmap=cmap)   # set levels as previous levels
# This is the fix for the white lines between contour levels (https://stackoverflow.com/questions/8263769/hide-contour-linestroke-on-pyplot-contourf-to-get-only-fills)
for j in c.collections:
    j.set_edgecolor("face")
plt.colorbar(c)
# plt.title('variance: fipy')
plt.xlabel('$x_1$', fontsize = 14)
plt.ylabel('$x_2$', fontsize = 14)
plt.legend(loc='best')
plt.tight_layout()
plt.savefig(os.path.join(os.getcwd(),'variance_fipy.pdf'))
plt.show()
## mid point x1=0.515,x2=0.515 statistics #location A
print MC_u_fipy[:,528].shape
print x1_f[528]
print x2_f[528]

fig, ax = plt.subplots()
ax.hist(MC_u_fipy[:,528], normed=True, alpha=0.5)
ax.set_xlabel('$U@x_1=0.515,x_2=0.515$')
ax.set_ylabel('ground_fipy_PDF')
plt.show()
## x1=0.359,x2=0.734 statistics #location B
print MC_u_fipy[:,747].shape
print x1_f[747]
print x2_f[747]

fig, ax = plt.subplots()
ax.hist(MC_u_fipy[:,747], normed=True, alpha=0.5)
ax.set_xlabel('$U@x_1=0.359,x_2=0.734$')
ax.set_ylabel('ground_fipy_PDF')
plt.show()<jupyter_output>(100000,)
[0.359375]
[0.734375]<jupyter_text>B.Neural_network<jupyter_code># ======================
##giving_DNN_architecture
# ------------------------------------------------------------

# ------------------------------------------------------------
# needs to be defined as activation class otherwise error
# AttributeError: 'Activation' object has no attribute '__name__'    
class Swish(Activation):
    
    def __init__(self, activation, **kwargs):
        super(Swish, self).__init__(activation, **kwargs)
        self.__name__ = 'swish'

def swish(x):
    return (K.sigmoid(x) * x)

get_custom_objects().update({'swish': Swish(swish)})
# ------------------------------------------------------------
# BUILDING DNN APPROXIMATOR
# ======================
# DNN network i/p:x1,x2,A o/p:prediction

x1 = Input(shape=(1,))
x2 = Input(shape=(1,))
A = Input(shape=(nx1*nx2,)) # input field image: conductivity image
a_val = Input(shape=(1,)) # input field value: conductivity value at the corresponding 'input(x1,x2)' location

if DNN_type == 'Resnet':
    x1_x2_A = concatenate([x1,x2,A])
    o = Dense(n)(x1_x2_A)
    for i in range(num_block):
        z = Dense(n, activation = act_func)(o)
        z = Dense(n, activation = act_func)(z)
        o = Add()([z, o])
    prediction = Dense(1)(o)
    print DNN_type

elif DNN_type == 'FC':
    num_neurons_per_layer = map(int, d.strip('[]').split(','))
    x1_x2_A = concatenate([x1,x2,A])
    z = Dense(num_neurons_per_layer[0], activation=act_func)(x1_x2_A)
    for n in num_neurons_per_layer[1:]:
        z = Dense(n, activation=act_func)(z)
    prediction = Dense(1)(z)
    print DNN_type

def myFunc(t):
    return ((1-t[0])+(t[0]*(1-t[0])*t[2])) # dirichlet on left and right sides of plate

u = Lambda(myFunc, output_shape=(1,))([x1,x2,prediction]) # field of interest : temperature
model = Model(inputs=[x1,x2,A], outputs=u)

init = tf.global_variables_initializer()   
sess = tf.Session()
K.set_session(sess)
sess.run(init) 

print('-------')
model.load_weights(DNN_weights_path)##loading saved weights of the model
# for layer in model.layers:
#     g=layer.get_config()
#     h=layer.get_weights()
#     print (g)
#     print (h)
# start = time.time()

# nsamples = MC_inputfield.shape[0]

# MC_samples_u_DNN = np.zeros((nsamples, nx1*nx2))  ##evaluating at cell centers

# for i in xrange(nsamples):
#     #display
#     if (i+1)%10000 == 0:
#         print "Predicting MC sample "+str(i+1)
        
#     a = MC_inputfield[ i , : ].reshape(-1,1)

#     u_DNN = sess.run(u, feed_dict={x1:x1_f, x2:x2_f,  A:np.repeat(a.T, np.shape(x1_f)[0], axis=0)})
#     MC_samples_u_DNN[i] = u_DNN.flatten()
    
# finish = time.time() - start
# print "Time (sec) to generate "+str(MC_inputfield.shape[0])+" MC samples using DNN: " +str(finish)

# np.save('MC_samples_u_DNN_exp_nx1=32_nx2=32_lx1=0.06_lx2=0.15_v=0.75_num_samples=100000.npy',MC_samples_u_DNN)
MC_u_DNN = np.load('MC_samples_u_DNN_exp_nx1=32_nx2=32_lx1=0.06_lx2=0.15_v=0.75_num_samples=100000.npy')
print MC_u_DNN.shape
print MC_u_DNN

mean_DNN = np.mean(MC_u_DNN,axis=0)
var_DNN = np.var(MC_u_DNN,axis=0)
print mean_DNN.shape
print var_DNN.shape
fig, ax = plt.subplots(figsize=(10,6))
c = ax.contourf(x1_f.reshape((nx1,nx2)), x2_f.reshape((nx1,nx2)), mean_DNN.reshape((nx1,nx2)),100,cmap=cmap)   # set levels as previous levels
# This is the fix for the white lines between contour levels (https://stackoverflow.com/questions/8263769/hide-contour-linestroke-on-pyplot-contourf-to-get-only-fills)
for j in c.collections:
    j.set_edgecolor("face")
plt.colorbar(c)
# plt.title('mean: fipy')
plt.xlabel('$x_1$', fontsize = 14)
plt.ylabel('$x_2$', fontsize = 14)
plt.legend(loc='best')
plt.tight_layout()
plt.savefig(os.path.join(os.getcwd(),'mean_DNN.pdf'))
plt.show()

fig, ax = plt.subplots(figsize=(10,6))
c = ax.contourf(x1_f.reshape((nx1,nx2)), x2_f.reshape((nx1,nx2)), var_DNN.reshape((nx1,nx2)),100,cmap=cmap)   # set levels as previous levels
# This is the fix for the white lines between contour levels (https://stackoverflow.com/questions/8263769/hide-contour-linestroke-on-pyplot-contourf-to-get-only-fills)
for j in c.collections:
    j.set_edgecolor("face")
plt.colorbar(c)
# plt.title('variance: fipy')
plt.xlabel('$x_1$', fontsize = 14)
plt.ylabel('$x_2$', fontsize = 14)
plt.legend(loc='best')
plt.tight_layout()
plt.savefig(os.path.join(os.getcwd(),'variance_DNN.pdf'))
plt.show()
## mid point x1=0.515,x2=0.515 statistics #location A
print MC_u_DNN[:,528].shape
print x1_f[528]
print x2_f[528]

fig, ax = plt.subplots()
ax.hist(MC_u_DNN[:,528], normed=True, alpha=0.5)
ax.set_xlabel('$U@x_1=0.515,x_2=0.515$')
ax.set_ylabel('DNN_PDF')
plt.show()
## x1=0.359,x2=0.734 statistics #location B
print MC_u_DNN[:,747].shape
print x1_f[747]
print x2_f[747]

fig, ax = plt.subplots()
ax.hist(MC_u_DNN[:,747], normed=True, alpha=0.5)
ax.set_xlabel('$U@x_1=0.359,x_2=0.734$')
ax.set_ylabel('DNN_PDF')
plt.show()<jupyter_output>(100000,)
[0.359375]
[0.734375]<jupyter_text>C.Comparing Fipy vs NN<jupyter_code>plt.rc('legend', fontsize=16) 
fig, ax = plt.subplots(figsize=(10,6))
sns.distplot(MC_u_fipy[:,528], bins =100, norm_hist = True,label='FVM', color="green", kde_kws=dict(linewidth=2),hist=False, hist_kws=dict(fill=False))
sns.distplot(MC_u_DNN[:,528], bins =100, norm_hist = True,label='DNN',color="red", kde_kws=dict(linewidth=2),hist=False, hist_kws=dict(fill=False))
ax.set_xlabel(r'$U(x=[0.515,0.515])$', fontsize = 16)
ax.set_ylabel('$PDF$', fontsize = 16)
plt.legend(loc="best")
plt.tight_layout()
plt.savefig(os.path.join(os.getcwd(),'PDF@x=locA_fipy_vs_DNN.pdf'))
plt.rc('legend', fontsize=16) 
fig, ax = plt.subplots(figsize=(10,6))
sns.distplot(MC_u_fipy[:,747], bins =100, norm_hist = True,label='FVM', color="green", kde_kws=dict(linewidth=2),hist=False, hist_kws=dict(fill=False))
sns.distplot(MC_u_DNN[:,747], bins =100, norm_hist = True,label='DNN',color="red", kde_kws=dict(linewidth=2),hist=False, hist_kws=dict(fill=False))
ax.set_xlabel(r'$U(x=[0.359,0.734])$', fontsize = 16)
ax.set_ylabel('$PDF$', fontsize = 16)
plt.legend(loc="best")
plt.tight_layout()
plt.savefig(os.path.join(os.getcwd(),'PDF@x=locB_fipy_vs_DNN.pdf'))<jupyter_output><empty_output><jupyter_text>Relative L2 ErrorLet $\mathbf{y}_{\mbox{true}}$ be the *ground truth vector* from fipy and $\hat{\mathbf{y}}$ be the *estimated vector* from neural network.The relattive L2 error is defined to be:$$L_2[\mathbf{y}_{\mbox{true}},\hat{\mathbf{y}}] = \frac{\parallel\mathbf{y}_{\mbox{true}}-\hat{\mathbf{y}} \parallel_2}{\parallel \mathbf{y}_{\mbox{true}}\parallel_2},$$where $\parallel\cdot\parallel_2$ is the standard Euclidean norm.<jupyter_code>## Relative L2 error in the mean and variance
y_true_mean=ground_mean_fipy
y_true_var=ground_var_fipy

y_cap_mean = mean_DNN
y_cap_var = var_DNN

L2_mean=(np.linalg.norm(y_true_mean-y_cap_mean) / np.linalg.norm(y_true_mean))
L2_var=(np.linalg.norm(y_true_var-y_cap_var) / np.linalg.norm(y_true_var))

print L2_mean
print L2_var 
from sklearn.metrics import r2_score
print r2_score(y_true_mean,y_cap_mean) 
print r2_score(y_true_var,y_cap_var)
### checking
# checking
print('*******')
h=[]
for i in xrange(MC_u_DNN.shape[0]):
    j=np.amax(MC_u_DNN[i,:])
    if j>=1.05:
        print j
        h.append(i)
print('*******')
print h
print len(h)
    
from matplotlib.backends.backend_pdf import PdfPages
pdf_pages = PdfPages('badplots.pdf')
for p in xrange(len(h)):
    # Initialize the plot
    fig = plt.figure(figsize=(18,6))

    try:
        ax1.lines.remove(c1[0])
        ax2.lines.remove(c2[0])
        ax3.lines.remove(c3[0])

    except:
        pass

    ax1 = fig.add_subplot(1, 3, 1)
    c1 = ax1.contourf( x1_f.reshape((nx1,nx2)), x2_f.reshape((nx1,nx2)), np.log(MC_inputfield)[h[p],:].reshape((nx1,nx2)), 100, cmap=cmap)   # set levels automatically
    # This is the fix for the white lines between contour levels (https://stackoverflow.com/questions/8263769/hide-contour-linestroke-on-pyplot-contourf-to-get-only-fills)
    for j in c1.collections:
        j.set_edgecolor("face")
    plt.colorbar(c1)
    plt.title('$log(Input \ field)$', fontsize=14)
    plt.xlabel('$x1$', fontsize=12)
    plt.ylabel('$x2$', fontsize=12)


    ax2 = fig.add_subplot(1, 3, 2)
    c2 = ax2.contourf( x1_f.reshape((nx1,nx2)), x2_f.reshape((nx1,nx2)), MC_u_fipy[h[p],:].reshape((nx1,nx2)), 100, cmap=cmap)  # set levels as previous levels
    # This is the fix for the white lines between contour levels (https://stackoverflow.com/questions/8263769/hide-contour-linestroke-on-pyplot-contourf-to-get-only-fills)
    for j in c2.collections:
        j.set_edgecolor("face")
    plt.colorbar(c2)
    plt.title('$FVM \ solution$',fontsize=14)
    plt.xlabel('$x1$', fontsize=12)
    plt.ylabel('$x2$', fontsize=12)


    ax3 = fig.add_subplot(1, 3, 3)
    c3 = ax3.contourf( x1_f.reshape((nx1,nx2)), x2_f.reshape((nx1,nx2)), MC_u_DNN[h[p],:].reshape((nx1,nx2)), 100, cmap=cmap)   # set levels automatically
    # This is the fix for the white lines between contour levels (https://stackoverflow.com/questions/8263769/hide-contour-linestroke-on-pyplot-contourf-to-get-only-fills)
    for j in c3.collections:
        j.set_edgecolor("face")
    plt.colorbar(c3)
    plt.title('$DNN \ solution$', fontsize=14)
    plt.xlabel('$x1$', fontsize=12)
    plt.ylabel('$x2$', fontsize=12)

    plt.tight_layout()
    plt.suptitle('test_case='+str(h[p]), fontsize=12) #indexed from 0
#     plt.savefig(os.path.join(resultdir,'test_case='+str(i+1)+'_DNN_type='+str(DNN_type)+'_nx1*nx2='+str(nx1*nx2)+'_seed='+str(seed)+'_'+str(variation)+'_nnpred-fipy.pdf')) 
    pdf_pages.savefig(fig)
    plt.show()
# # Write the PDF document to the disk
pdf_pages.close()
#check<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1827214431682052, 'fraction_numerical': 0.04883376788105137, 'mean_word_length': 5.15704989154013, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 7, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '247238', 'license': 'MIT', 'n_tokens_mistral': 6218, 'n_tokens_neox': 5861, 'n_words': 996, 'path': 'scripts/2d_spde_prob/uq/uqcase_2/2dim_MC_code.ipynb', 'repo_name': 'PredictiveScienceLab/variational-elliptic-SPDE', 'size': '14191'}