From 7ed47ed4ad87a68ba4955e4eca2d68fd236f5df3 Mon Sep 17 00:00:00 2001
From: JanHab <jan.habscheid@rwth-aachen.de>
Date: Wed, 26 Mar 2025 10:35:32 +0100
Subject: [PATCH] =?UTF-8?q?add=20given=20code=C2=A7?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
.../mcmc_task2_u0X1_final_student_mars21.py | 422 ++++++++++++++++++
Project3/src/GivenCode/simulations_Feb20.py | 75 ++++
Project3/src/GivenCode/solver_Feb20.py | 112 +++++
3 files changed, 609 insertions(+)
create mode 100644 Project3/src/GivenCode/mcmc_task2_u0X1_final_student_mars21.py
create mode 100644 Project3/src/GivenCode/simulations_Feb20.py
create mode 100644 Project3/src/GivenCode/solver_Feb20.py
diff --git a/Project3/src/GivenCode/mcmc_task2_u0X1_final_student_mars21.py b/Project3/src/GivenCode/mcmc_task2_u0X1_final_student_mars21.py
new file mode 100644
index 0000000..c46cea7
--- /dev/null
+++ b/Project3/src/GivenCode/mcmc_task2_u0X1_final_student_mars21.py
@@ -0,0 +1,422 @@
+import matplotlib.pyplot as plt
+import numpy as np
+#import scipy.io as sio
+import scipy.sparse as sparse
+
+try:
+ from simulations_Feb20 import simulate_g
+except ImportError:
+ from python_scripts.simulations_Feb20 import simulate_g
+
+
+##############################################
+# read observation data provided by teacher #
+##############################################
+
+# use only the first data set
+observation = np.load("data/True_Observation_data_1_upd.npz")
+#print( list(observation.keys()) )
+
+
+######################################################################################
+# In the following different parameters required for solving conslaw are extracted #
+######################################################################################
+
+
+# number of observation data set
+Ant_Init = 1
+
+#y_k_true = observation["y_k_true"]
+#y_k_true = np.zeros_like(y_k_true)
+
+# provide information about grid,
+x = observation["x"]
+xh = observation["xh"]
+
+# grid cells for cons law
+N = observation["M"].item()
+
+# Final time
+T = observation["T"].item()
+
+# describe time points for which we have observations
+time_obs = observation["t_obspoint"]
+JX = time_obs <= T
+
+Ldom = observation["Ldom"].item()
+dx = Ldom/N
+
+# positions of point values for k(x)
+x_k = observation["x_k"]
+
+# index vector for numbering positions in space and time
+Nt = np.sum(JX)
+Ix = np.arange(0, N)
+It = np.arange(0, Nt)
+
+# number of unknown points to be determined in k(x) which is assumed piecewise linear
+Mr = observation["Mr"].item()
+
+###############################
+# Observation data h(theta) #
+###############################
+u_true_X1 = np.zeros((Nt))
+u_true_X2 = np.zeros((Nt))
+u_true_X3 = np.zeros((Nt))
+u_true_X4 = np.zeros((Nt))
+
+# correct when Ant_Init = 1
+u_true_X1 = observation["u_true_X1"]
+u_true_X2 = observation["u_true_X2"]
+u_true_X3 = observation["u_true_X3"]
+u_true_X4 = observation["u_true_X4"]
+
+
+# provide initial data set: each initialdata (of total=Ant_Init) is a column
+initial_true = np.zeros((N,2))
+initial_true[:,0] = observation["u0"]
+
+
+# provide position of observation points X_1, X_2, X_3, X_4
+X_position = observation["X_position"]
+X_1 = X_position[0]
+X_2 = X_position[1]
+X_3 = X_position[2]
+X_4 = X_position[3]
+
+
+####################################
+# define parameters used in MCMC #
+####################################
+add_unc = 0.03 # noise in observation data
+error_val = 100 # lower limit in Loss function
+kk_tolerance = 1000 # number of steps in MCMC algorithm
+teller = 1
+tau = 20 # extract "theta^k" with step size "tau" to represent a posteriori "theta"
+
+N_ens = int(kk_tolerance/tau)
+
+# function to create one long column vector to represent the matrix of
+# 4 column vectors representing observation data at X1, X2, X3, X4
+def gen_measurements(u, weights):
+ measurement = []
+ for i, u in enumerate(u):
+ m = u + add_unc * np.random.normal(size=Nt * Ant_Init)
+ m[m < 0] = 0
+ m[m > 1] = 1
+ measurement.append(weights[i] * m)
+ return np.concatenate(measurement)
+
+
+# weight factor used to adjust the importance of the 4 different time series
+weights = [1, 1, 1, 1]
+
+# compute one single column vector to represent all observation data
+measurement = gen_measurements([u_true_X1, u_true_X2, u_true_X3, u_true_X4], weights)
+
+
+# parametervector "theta" which is updated by the mcmc algorithm
+# theta is composed of both k(x) and f(u) - however, only k(x) unknown
+initial_ensemble_kx = np.zeros((Mr, 1))
+
+# we generate an ensemble parameter vector by extracting every 20th member from "initial_ensemble_F"
+updated_ensemble_kx = np.zeros((Mr, N_ens))
+
+# mcmc parameter
+#beta = 0.05 * 4 #
+beta = 0.05 #
+#beta = 0.05 * 0.25 #
+
+# we generate initial k(x) ~ {k_i}_i=1^9 as random normal distribution N(mu,sigma)
+mu = 1
+sigma = 0.25
+
+
+
+#########################################################
+# Choose what should be the initial parameter vector #
+#########################################################
+
+
+
+mcmc_start = True
+#mcmc_start = False
+
+if mcmc_start:
+ #####################################
+ # choice 1: generate theta^{(0)} #
+ #####################################
+
+ print("Initial guess of k(x)")
+ kk_x = np.random.normal(loc=mu, scale=sigma, size=Mr)
+ kk_x[kk_x < 0.95] = 0.95
+ kk_x[kk_x > 1.3] = 1.3
+ kkx_u = kk_x
+else:
+ #######################################
+ # choice 2: start with theta^{(k)} #
+ #######################################
+
+ with np.load('data/mcmc_u0X1_sample100_it1_student_Mars21.npz') as data_kx:
+ #with np.load('data/mcmc_u0X1_sample100_it2_student_Mars21.npz') as data_kx:
+ #with np.load('data/mcmc_u0X1_sample100_it3_student_Mars21.npz') as data_kx:
+ #with np.load('data/mcmc_u0X1_sample100_it4_student_Mars21.npz') as data_kx:
+ kkx_u = data_kx["theta_1"]
+
+
+# Create the y_k vector with padding 1s at the beginning and the end
+y_k = np.concatenate(([1], kkx_u, [1]))
+
+
+# number of parameter vectors (in case of more than one ensemble)
+K = 0
+
+
+# Create the plot
+plt.figure(1)
+plt.plot(x_k, y_k)
+plt.axis([-0.1, Ldom + 0.1, 0.1, 1.9])
+plt.grid(True)
+plt.xlabel("X")
+plt.ylabel("k(x)")
+plt.title("Resistance k(x)")
+
+# plt.ioff()
+
+# Show the plot and pause execution to view it
+plt.show(block=False)
+plt.pause(2) # Pause for a short period of time; adjust as needed
+
+# Close the plot after pausing
+plt.close(1)
+
+# "initial_ensemble_kx" is composed of vector k(x) (M:M+Mr)
+# which is hidden and we seek to recover
+initial_ensemble_kx[0:Mr, K] = kkx_u
+
+#######################################################
+# Run simulation to generate h(kkx_u)
+#######################################################
+
+###############################
+# Step S3: Compute h(theta) #
+###############################
+
+# h ~ simulate forward model
+# theta ~ initial_ensemble_kx
+
+# u_X =
+# ensembleOfPredictedObservations =
+
+ensembleOfPredictedObservations = simulate_g(
+ N, x, xh, Ldom, T,
+ initial_true,
+ Ant_Init,
+ Nt,
+ time_obs,
+ It,
+ JX,
+ K,
+ weights,
+ initial_ensemble_kx,
+ X_1, X_2, X_3, X_4,
+ x_k, kkx_u,
+ u_true_X1, u_true_X2, u_true_X3, u_true_X4,
+ #to_plot=False,
+ #to_plot=True,
+)
+
+
+# do a numbering of the column vector of predicted observation = 1 long column vector
+# in case we want to throw away some observation data
+measIndex = np.arange(0, ensembleOfPredictedObservations.shape[0])
+weight_obs = sparse.diags(1 / (add_unc * np.ones(measIndex.shape[0])) ** 2)
+
+
+# Compute loss function
+meas_diff = ensembleOfPredictedObservations[measIndex, K] - measurement[measIndex]
+Phi_kkx = meas_diff.T @ weight_obs @ meas_diff
+
+
+# Initialization
+theta_1 = kkx_u.copy()
+Phi_theta_1 = Phi_kkx.copy()
+
+Loss = -1 * np.ones(kk_tolerance + 1)
+Loss[0] = Phi_theta_1.copy()
+
+input("Now entering main loop: \n Press enter to continue!")
+
+# counter of steps in MCMC
+kk = 0
+
+while Loss[kk] > error_val and kk < kk_tolerance:
+
+
+ #####################################################################
+ # Step S2: compute "phi=gu" as a perturbation of "theta=fu_kp1" #
+ #####################################################################
+
+ phi_1 = np.sqrt((1 - beta**2)) * theta_1 + beta * sigma * np.random.normal(size = Mr)
+
+ phi_1[phi_1 < 0.25] = 0.25
+ phi_1[phi_1 > 1.75] = 1.75
+
+ # prepare for solving the conservation law to obtain h(phi)
+ initial_ensemble_kx[np.arange(0, Mr), K] = phi_1
+
+ # Run simulation to generate G(gu)
+ # u_X is the generated time series to be compared with data
+
+
+ ###############################
+ # Step S3: Compute h(phi) #
+ ###############################
+ ensembleOfPredictedObservations = simulate_g(
+ N, x, xh, Ldom, T,
+ initial_true,
+ Ant_Init,
+ Nt,
+ time_obs,
+ It,
+ JX,
+ K,
+ weights,
+ initial_ensemble_kx,
+ X_1, X_2, X_3, X_4,
+ x_k, theta_1,
+ u_true_X1, u_true_X2, u_true_X3, u_true_X4
+ )
+
+ # Compute Loss
+ meas_diff = ensembleOfPredictedObservations[measIndex, K] - measurement[measIndex]
+ Phi_phi_1 = meas_diff.T @ weight_obs @ meas_diff
+
+ ########################################
+ # Step S3: Compute a(theta_1, phi_1) #
+ ########################################
+
+ # MCMC a(theta,phi)
+ a_theta_1_phi_1 = np.min([1, np.exp(Phi_theta_1 - Phi_phi_1)] )
+
+ ##########################################################
+ # Step S4 and S5: compare "choice" and "a(theta, phi)" #
+ #########################################################
+
+ choice = np.random.uniform(0, 1)
+
+ if choice <= a_theta_1_phi_1:
+ theta_1 = phi_1.copy()
+ Phi_theta_1 = Phi_phi_1.copy()
+
+ else:
+ theta_1 = theta_1.copy()
+ Phi_theta_1 = Phi_theta_1.copy()
+
+ #####################################################
+ # check if we should extract this parameter vector #
+ #####################################################
+
+ check = kk / tau
+
+ if check == teller:
+ updated_ensemble_kx[0 : Mr, teller] = theta_1
+ teller += 1
+
+ fig, axs = plt.subplots(1, 2, figsize=(15, 5))
+
+ #axs[0].plt.figure()
+ axs[0].to_plot = Loss[Loss >= 0]
+ axs[0].plot(Loss)
+ axs[0].set_title("Loss function")
+ axs[0].set_ylabel("Loss")
+ axs[0].set_xlabel("Iteration number")
+ axs[0].grid(True)
+ #axs[0].show()
+
+ y_k_ident = np.concatenate([[1], theta_1, [1]])
+
+ #axs[1].plt.figure()
+ #axs[1].plot(x_k, y_k_true, '-r')
+ axs[1].plot(x_k, y_k_ident,'-b')
+ axs[1].set_title("k(x) function")
+ axs[1].set_ylabel("k(x)")
+ axs[1].set_xlabel("x")
+ #axs[1].legend(['True', 'Identified'])
+ axs[1].legend(['Identified'])
+ plt.axis([0, Ldom, 0.25, 1.75])
+ axs[1].grid(True)
+ #axs[1].show()
+
+ # Update the plot
+ #plt.draw()
+
+ # Show the plots
+ plt.show()
+ #plt.pause(2)
+ #input("Extracting this parameter vector: \n Press enter to continue!")
+
+ #############################################
+ # Step S6: move to next iteration in MCMC #
+ #############################################
+
+ kk += 1
+ Loss[kk] = Phi_theta_1.copy()
+
+ y_k_ident = np.concatenate([[1], theta_1, [1]])
+
+ print("Iteration:", kk, " Loss:", Loss[kk])
+
+
+######################################################
+# Save to file variables from this MCMC compuation #
+######################################################
+
+np.savez(
+ "data/mcmc_u0X1_sample100_it1_student_Mars21.npz",
+ #"data/mcmc_u0X1_sample100_it2_student_Mars21.npz",
+ #"data/mcmc_u0X1_sample100_it3_student_Mars21.npz",
+ # "data/mcmc_u0X1_sample100_it4_student_Mars21.npz",
+ # "data/mcmc_u0X1_sample100_it5_student_Mars21.npz",
+ #"data/mcmc_u0X1_sample100_it6_student_Mars21.npz",
+ Loss=Loss,
+ #fu_kp1=fu_kp1,
+ #Phi_fu_kp1_y=Phi_fu_kp1_y,
+ theta_1=theta_1,
+ Phi_theta_1=Phi_theta_1,
+ updatedEnsemble_kx = updated_ensemble_kx,
+ kk=kk,
+ teller=teller,
+ x_k=x_k,
+ #y_k_true=y_k_true,
+ y_k_ident=y_k_ident,
+ beta=beta,
+)
+print("Sampling of 100 members done")
+print("Final loss was:", Loss)
+print("And the lower error val was set to:", error_val)
+
+
+# Plot the loss function
+plt.figure()
+to_plot = Loss[Loss >= 0]
+plt.plot(Loss)
+plt.title("Loss function")
+plt.ylabel("Loss")
+plt.xlabel("Iteration number")
+#plt.axis([0, Ldom, 0.5, 2.0])
+plt.grid(True)
+plt.show()
+
+plt.figure()
+#plt.plot(x_k, y_k_true, '-r')
+plt.plot(x_k, y_k_ident, '-b')
+plt.title("k(x) function")
+plt.ylabel("k(x)")
+plt.xlabel("x")
+plt.axis([0, Ldom, 0.25, 1.75])
+#plt.legend(['True', 'Identified'])
+plt.legend(['Identified'])
+plt.grid(True)
+plt.show()
+
diff --git a/Project3/src/GivenCode/simulations_Feb20.py b/Project3/src/GivenCode/simulations_Feb20.py
new file mode 100644
index 0000000..571b913
--- /dev/null
+++ b/Project3/src/GivenCode/simulations_Feb20.py
@@ -0,0 +1,75 @@
+#import matplotlib.pyplot as plt
+import numpy as np
+#import scipy.sparse as sparse
+
+try:
+ from python_scripts.solver_Feb20 import my_solver
+except ImportError:
+ from solver_Feb20 import my_solver
+
+
+def simulate_g(
+ N, x, xh, Ldom, Tend,
+ initial_true,
+ Ant_Init,
+ Nt,
+ time_obs,
+ It,
+ JX,
+ K,
+ weights,
+ initial_ensemble_kx,
+ X_1, X_2, X_3, X_4,
+ x_k, kx_old,
+ u_true_X1, u_true_X2, u_true_X3, u_true_X4,
+ #to_plot=False,
+ #to_plot = True,
+):
+ numEns = 1
+
+ # vector that can store u(x1,t) for "Ant_Init" different initial data in one column vector
+ ensembleOfPredictedObs_X1_0toT3 = np.zeros((Nt * Ant_Init, numEns))
+ ensembleOfPredictedObs_X2_0toT3 = np.zeros((Nt * Ant_Init, numEns))
+ ensembleOfPredictedObs_X3_0toT3 = np.zeros((Nt * Ant_Init, numEns))
+ ensembleOfPredictedObs_X4_0toT3 = np.zeros((Nt * Ant_Init, numEns))
+
+ # list that can store u(x1,t), u(x2,t), u(x3,t), u(x4,t) in column vector for each initial_data #K0
+ u_X = [np.zeros((len(time_obs), Ant_Init)) for _ in range(4)]
+
+ #####################################
+ # Loop to simulate Ant_Init cases #
+ #####################################
+
+ for K_0 in range(0, Ant_Init):
+
+ u_X[0][:, K_0], u_X[1][:, K_0], u_X[2][:, K_0], u_X[3][:, K_0], y_k = my_solver(
+ N, x, xh, Ldom, Tend,
+ initial_ensemble_kx[:, K],
+ initial_true[:, K_0],
+ time_obs,
+ [X_1, X_2, X_3, X_4],
+ x_k,
+ #to_plot=to_plot,
+ #to_plot = False,
+ #to_plot = True,
+ )
+
+ # Store the obtained 4 time series "u_X" in the "ensembleOfPredictedObs_X1_0toT3"
+
+ ensembleOfPredictedObs_X1_0toT3[It + (K_0) * Nt, K] = u_X[0][JX, K_0]
+ ensembleOfPredictedObs_X2_0toT3[It + (K_0) * Nt, K] = u_X[1][JX, K_0]
+ ensembleOfPredictedObs_X3_0toT3[It + (K_0) * Nt, K] = u_X[2][JX, K_0]
+ ensembleOfPredictedObs_X4_0toT3[It + (K_0) * Nt, K] = u_X[3][JX, K_0]
+
+
+ # generate "ensembleOfPredictedObservations" which possess same structure as "measurement"
+ ensembleOfPredictedObservations = np.concatenate(
+ [
+ weights[0] * ensembleOfPredictedObs_X1_0toT3,
+ weights[1] * ensembleOfPredictedObs_X2_0toT3,
+ weights[2] * ensembleOfPredictedObs_X3_0toT3,
+ weights[3] * ensembleOfPredictedObs_X4_0toT3,
+ ]
+ )
+ # return u_X, ensembleOfPredictedObservations
+ return ensembleOfPredictedObservations
diff --git a/Project3/src/GivenCode/solver_Feb20.py b/Project3/src/GivenCode/solver_Feb20.py
new file mode 100644
index 0000000..5721f0d
--- /dev/null
+++ b/Project3/src/GivenCode/solver_Feb20.py
@@ -0,0 +1,112 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def my_solver(
+ #M,
+ N,
+ x,
+ xh,
+ Ldom,
+ Tend,
+ Fluxkx_ensem,
+ initial_true,
+ time_obs,
+ X,
+ x_k,
+ #to_plot=False,
+ #to_plot=True,
+):
+ """
+ Solution of u_t + (k(x)f(u))_x = 0
+ """
+
+ # Final time
+ T = Tend
+
+ # Delta x
+ dx = Ldom / N
+
+ # observation points in spatial domain
+ s = [ (np.sum(x <= X[i]) - 1) for i in range(4)]
+
+ # output predicted observations
+ u_X_out = [np.zeros(len(time_obs)) for _ in range(4)]
+
+ # terrain function
+ k_val = Fluxkx_ensem[:]
+ #y_k = np.concatenate((np.array([1]), k_val, np.array([1])))
+ y_k = np.concatenate(([1.0], k_val, [1.0] ))
+
+ # Define the flux function f(u)
+ def fun_flux(u):
+ f = np.zeros_like(u)
+ f = u * (1 - u)
+ return f
+
+
+ # estimate the time step dt based om CFL condition
+ dt_estimat = dx * (3/7) # based on max |f'| = 1, max k(x) < 4/3
+ NTime = round(T / dt_estimat)
+
+ dt = T / NTime
+
+ #MA = max(maxfu_prime, 0.1)
+
+ MA = 1 # Rusanov flux is used
+ lmd = dt / dx
+
+ t_obs = np.linspace(0, T, NTime + 1)
+
+ J1 = np.arange(1, N - 1)
+ J2 = np.arange(0, N - 1)
+
+ u0 = initial_true[0:N]
+
+ u = np.zeros(N)
+
+ # teller = 1
+
+ u_old = u0.copy()
+
+ # Compute predicted time serie state at t=0
+ u_X = [np.zeros(NTime + 1) for i in range(4)]
+ for ind in range(4):
+ u_X[ind][0] = u0[s[ind]]
+
+
+
+ for j in range(NTime):
+
+ # fluxes at grid interface
+ F_half = np.zeros(N - 1)
+ Flux = np.zeros(N) # compute flux for all cells
+
+ Flux = fun_flux(u_old)
+ F_half[J2] = 0.5 * (Flux[J2] + Flux[J2 + 1]) - 0.5 * MA * (u_old[J2 + 1] - u_old[J2] )
+
+ k_half = np.zeros(N - 1)
+ k_half[J2] = np.interp(xh[J2 + 1], x_k, y_k)
+
+ # The numerical scheme
+
+ u[J1] = u_old[J1] - lmd * (k_half[J1] * F_half[J1] - k_half[J1 - 1] * F_half[J1 - 1])
+
+ u[u > 1] = 1.0
+ u[u < 0] = 0.0
+
+ u[0] = u[1]
+ u[-1] = u[-2]
+
+ # Extract predicted observation data
+ for ind in range(4):
+ u_X[ind][j+1] = u[s[ind]]
+
+
+ u_old = u.copy()
+
+ # transform predicted observation over to same time grid as used for true observation
+ for ind in range(4):
+ u_X_out[ind][:] = np.interp(time_obs[:], t_obs[:], u_X[ind][:])
+
+ return *u_X_out, y_k
\ No newline at end of file
--
GitLab