<jupyter_start><jupyter_text>Imports<jupyter_code>import tensorflow as tf
print(tf.version)
print(tf.test.is_built_with_cuda())
print(tf.test.is_gpu_available(cuda_only=False, min_cuda_compute_capability=None))
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# from plotting import newfig, savefig
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.animation import FuncAnimation, PillowWriter
import numpy as np
import scipy.io
from scipy.interpolate import griddata
import time
from pyDOE import lhs
import pickle as pkl
%matplotlib widget<jupyter_output><empty_output><jupyter_text>Equation<jupyter_code>k = 1<jupyter_output><empty_output><jupyter_text>$$\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2}$$  Load Data  Model<jupyter_code># Layers
u_layers = [2, 50, 50, 50, 50, 1]
pde_layers = [3, 100, 100, 1]

layers = [2, 50, 50, 50, 50, 1]
# tf placeholders for Identification
t_tf = tf.placeholder(tf.float32, shape=[None, 1])
x_tf = tf.placeholder(tf.float32, shape=[None, 1])
u_tf = tf.placeholder(tf.float32, shape=[None, 1])
t_tf, x_tf, u_tf
def initialize_NN(layers):
    weights = []
    biases = []
    num_layers = len(layers)
    for l in range(0, num_layers - 1):
        W = xavier_init(size=[layers[l], layers[l + 1]])
        b = tf.Variable(tf.zeros([1, layers[l + 1]], dtype=tf.float32),
                        dtype=tf.float32)
        weights.append(W)
        biases.append(b)
    return weights, biases

def xavier_init(size):
    in_dim = size[0]
    out_dim = size[1]
    xavier_stddev = np.sqrt(2 / (in_dim + out_dim))
    return tf.Variable(tf.truncated_normal([in_dim, out_dim],
                                           stddev=xavier_stddev,
                                           dtype=tf.float32),
                       dtype=tf.float32)
def neural_net(X, weights, biases):
    num_layers = len(weights) + 1
    H = X
    for l in range(0, num_layers - 2):
        W = weights[l]
        b = biases[l]
        H = tf.sin(tf.add(tf.matmul(H, W), b))
    W = weights[-1]
    b = biases[-1]
    Y = tf.add(tf.matmul(H, W), b)
    return Y
weights, biases = initialize_NN(layers)
# weights, biases
# load weights and biases

with open(__NAME + '/weights.pkl', 'rb') as db_file:
    W_pkl = pkl.load(db_file)

with open(__NAME + '/biases.pkl', 'rb') as db_file:
    B_pkl = pkl.load(db_file)

W = []
B = []
for w, b in zip(W_pkl, B_pkl):
    W.append(tf.Variable(w))
    B.append(tf.Variable(b))

weights = W
biases = B
lb_tf = tf.placeholder(tf.float32, shape=[2])
ub_tf = tf.placeholder(tf.float32, shape=[2])
# tf placeholders for Solution
t0_tf = tf.placeholder(tf.float32, shape=[None, 1])
x0_tf = tf.placeholder(tf.float32, shape=[None, 1])
u0_tf = tf.placeholder(tf.float32, shape=[None, 1])

t_lb_tf = tf.placeholder(tf.float32, shape=[None, 1])
x_lb_tf = tf.placeholder(tf.float32, shape=[None, 1])
t_ub_tf = tf.placeholder(tf.float32, shape=[None, 1])
x_ub_tf = tf.placeholder(tf.float32, shape=[None, 1])
u_ub_tf = tf.placeholder(tf.float32, shape=[None, 1])
u_lb_tf = tf.placeholder(tf.float32, shape=[None, 1])

t_f_tf = tf.placeholder(tf.float32, shape=[None, 1])
x_f_tf = tf.placeholder(tf.float32, shape=[None, 1])
def sol_net_u(t, x):
    X = tf.concat([t, x], 1)
    H = 2.0 * (X - lb_tf) / (ub_tf - lb_tf) - 1.0
    u = neural_net(H, weights, biases)
    u_x = tf.gradients(u, x)[0]
    return u, u_x

def sol_net_f(t, x):
    u, u_x = sol_net_u(t, x)

    u_t = tf.gradients(u, t)[0]

    u_xx = tf.gradients(u_x, x)[0]

    f = u_t - k * u_xx

    return f
# tf graphs for Solution
u0_pred, _ = sol_net_u(t0_tf, x0_tf)
u_lb_pred, u_x_lb_pred = sol_net_u(t_lb_tf, x_lb_tf)
u_ub_pred, u_x_ub_pred = sol_net_u(t_ub_tf, x_ub_tf)
sol_f_pred = sol_net_f(t_f_tf, x_f_tf)
# loss for Solution
sol_loss = tf.reduce_sum(tf.square(u0_tf - u0_pred)) + \
           tf.reduce_sum(tf.square(u_ub_tf - u_ub_pred)) + \
           tf.reduce_sum(tf.square(u_lb_tf - u_lb_pred)) + \
           tf.reduce_sum(tf.square(sol_f_pred))
# Optimizer for Solution
sol_optimizer = tf.contrib.opt.ScipyOptimizerInterface(
    sol_loss,
    var_list = weights + biases,
    method='L-BFGS-B',
    options={
        'maxiter': 50000,
        'maxfun': 50000,
        'maxcor': 50,
        'maxls': 50,
        'ftol': 1.0 * np.finfo(float).eps
    })
adam_optimizer = tf.train.AdamOptimizer()
sol_train_op_Adam = adam_optimizer.minimize(
            sol_loss,
            var_list= weights + biases)
# tf session
sess = tf.Session(config=tf.ConfigProto(
    allow_soft_placement=True, log_device_placement=True))
init = tf.global_variables_initializer()
sess.run(init)<jupyter_output>Device mapping:
/job:localhost/replica:0/task:0/device:GPU:0 -> device: 0, name: NVIDIA GeForce GTX 1650 Ti, pci bus id: 0000:01:00.0, compute capability: 7.5<jupyter_text>Training  Prepare data<jupyter_code>lb = np.array([0.0, 0.0])
ub = np.array([10.0, 5.0])
N = 10
fig = plt.figure()
ax = fig.gca()
ax.set_xlim(lb[0], ub[0])
ax.set_ylim(lb[1], ub[1])
ax.set_xticks(np.arange(lb[0],ub[0],(ub[0] - lb[0])/N))
ax.set_yticks(np.arange(lb[1],ub[1],(ub[1] - lb[1])/N))
plt.grid()
l = lb + (ub - lb) * lhs(2, N) 
plt.scatter(l[:, 0], l[:, 1], color="r", label="lhs")
plt.title("Latin Hypercube Sampling\nN=10")
ax.set_xlabel('$t$')
ax.set_ylabel('$x$')
fig.set_figheight(3.8)
fig.set_figwidth(6)
plt.tight_layout()

N0 = 200
N_b = 400
N_f = 20000
(N0, N_b, N_f)
t_data = np.linspace(lb[0], ub[0], N_b)[:, None]
x_data = np.linspace(lb[1], ub[1], N0)[:, None]
T_sol, X_sol = np.meshgrid(t_data, x_data)
# U_sol = u_data
X_sol_star = np.hstack(
    (T_sol.flatten()[:, None],
    X_sol.flatten()[:, None])
    )
# U_sol_star = U_sol.flatten()[:, None]

print(X_sol_star.shape, X_sol_star[0:10], sep='\n')
L = 5
u_data = (6 * np.sin(np.pi * x_data / L))
X0 = np.hstack((T_sol[:, 0:1], X_sol[:, 0:1])) # left boundary
X_lb = np.hstack((T_sol[0:1, :].T, X_sol[0:1, :].T)) # lower boundary
X_ub = np.hstack((T_sol[0:1, :].T, np.repeat(ub[1], t_data.shape[0])[:, None])) # upper boundary
# shuffled initial boundary data (left boundary)
idx_x = np.random.choice(x_data.shape[0], N0, replace=False)
X0_train = X0[idx_x, :]
u0_train = u_data[idx_x, 0:1]
# shuffle time data
idx_t = np.random.choice(t_data.shape[0], N_b, replace=False)
tb_train = t_data[idx_t, :]
X_f_train = lb + (ub - lb) * lhs(2, N_f) 
fig, ax = plt.subplots()
ax.set_xlim(lb[0] -0.1, ub[0])
ax.set_ylim(lb[1] - 0.4, ub[1] + 0.4)
fig.set_figheight(3.2)
fig.set_figwidth(6)

ax.scatter(X0_train[:, 0], X0_train[:, 1], s=4, marker='.')
ax.scatter(tb_train[:, 0], np.repeat(lb[1], N_b), s=4, marker='.')
ax.scatter(tb_train[:, 0], np.repeat(ub[1], N_b), s=4, marker='.')
ax.scatter(X_f_train[:, 0], X_f_train[:, 1], s=4, marker='.', edgecolors='none')
# ax.imshow(u0_train, extent=(t_data[0, 0], t_data[1, 0], x0_train.max(), x0_train.min()), aspect='auto')
plt.tight_layout()
fig, ax = plt.subplots()
ax.set_xlabel('$x$')
ax.set_ylabel('$u$')
ax.plot(x_data[:, 0], (np.sin(np.pi * x_data / L))[:, 0:1])
fig.set_figheight(3.2)
fig.set_figwidth(6)
plt.tight_layout()
def callback(loss):
    print('Loss: %e' % (loss))
tf_dict = {
    lb_tf: lb,
    ub_tf: ub,
    t0_tf: X0_train[:, 0:1],
    x0_tf: X0_train[:, 1:2],
    u0_tf: u0_train,
    t_lb_tf: X_lb[:, 0:1],
    x_lb_tf: X_lb[:, 1:2],
    t_ub_tf: X_ub[:, 0:1],
    x_ub_tf: X_ub[:, 1:2],
    u_lb_tf: np.repeat(0, N_b)[:, None],
    u_ub_tf: np.repeat(0, N_b)[:, None],
    t_f_tf: X_f_train[:, 0:1],
    x_f_tf: X_f_train[:, 1:2]
}
start_time = time.time()
it = 0
end = False
while not(end):

    sess.run(sol_train_op_Adam, tf_dict)

    # Print
    if it % 10 == 0:
        elapsed = time.time() - start_time
        loss_value = sess.run(sol_loss, tf_dict)
        print('It: %d, Loss: %.3e, Time: %.2f' %
                (it, loss_value, elapsed))
        start_time = time.time()

        if loss_value < 5 * 10**(-3):
            end = True

    it = it + 1
sol_optimizer.minimize(sess,
                       feed_dict=tf_dict,
                       fetches=[sol_loss],
                       loss_callback=callback)
sess.run(sol_loss, feed_dict=tf_dict)
with open(__NAME + '/weights.pkl', 'wb') as db_file:
    pkl.dump(obj=sess.run(weights), file=db_file)

with open(__NAME + '/biases.pkl', 'wb') as db_file:
    pkl.dump(obj=sess.run(biases), file=db_file)
X_sol_star[:, 0:1]
u_pred = sess.run(u0_pred, {
        lb_tf: lb,
        ub_tf: ub,
        t0_tf: X_sol_star[:, 0:1],
        x0_tf: X_sol_star[:, 1:2]
    })
fig = plt.figure(figsize=(4*1.75,4), dpi=200)
ax = fig.gca()

ax.set_xlim(lb[0], ub[0])
ax.set_ylim(lb[1], ub[1])

# plt.subplots_adjust(bottom=0.17)
# plt.subplots_adjust(left=0.17)

plt.title('T')
ax.set_xlabel('$t$')
ax.set_ylabel('$x$')

plt.pcolormesh(np.reshape(X_sol_star[:, 0], (N0, -1)), 
               np.reshape(X_sol_star[:, 1], (N0, -1)), 
               np.reshape(u_pred[:, 0], (N0, -1)), 
               shading='gouraud',  cmap='jet')
plt.colorbar()

plt.tight_layout()

# plt.legend()

fig.savefig('Figures\\Heat 1.png')
t = np.reshape(X_sol_star[:, 0], (N0, -1))
x = np.reshape(X_sol_star[:, 1], (N0, -1))
u = np.reshape(u_pred[:, 0], (N0, -1))

x_init = x[:, 0]
u_init = u[:, 0]

fig = plt.figure(figsize=(4*1.75,4), dpi=200)
ax = fig.gca()

ax.set_xlim(lb[1], ub[1])

ax.yaxis.grid(color='gainsboro', linestyle='dotted', linewidth=1.5)
ax.xaxis.grid(color='gainsboro', linestyle='dotted', linewidth=0.8)
ax.axhline(0,linestyle='dotted', color='grey')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')

plt.title('t = 0')

ax.set_xlabel('$x$')
ax.set_ylabel('$T$')

plt.tight_layout()

ln, = ax.plot(x_init, u_init)

def update(frame):
    plt.title('t = {time:.2f}'.format(time = t[0, frame]))
    ln.set_data(x[:, frame], u[:, frame])

ani = FuncAnimation(fig, update, list(range(0, N_b)))
ani.event_source.stop()
writer = PillowWriter(fps=25) 
ani.save("Figures\\Heat 1.gif", writer=writer)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17825918010491548, 'fraction_numerical': 0.036428987759860115, 'mean_word_length': 3.293160967472894, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '132330', 'license': 'MIT', 'n_tokens_mistral': 4674, 'n_tokens_neox': 4458, 'n_words': 821, 'path': 'Heat eqn/Heat eqn.ipynb', 'repo_name': 'Ragav-KS/PINNs-Review', 'size': '10294'}