add given code§

Project3/src/GivenCode/mcmc_task2_u0X1_final_student_mars21.py

+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()
+

Project3/src/GivenCode/simulations_Feb20.py

+#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

Project3/src/GivenCode/solver_Feb20.py

+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