Improved compressed sensing

examples/fixed_source/cooper1/input.py

@@ -51,12 +51,20 @@
 
 mcdc.tally.mesh_tally(
     scores=["flux"],
-    x=np.linspace(0.0, 4.0, 40),
-    y=np.linspace(0.0, 4.0, 40),
+    x=np.linspace(0.0, 4.0, 41),
+    y=np.linspace(0.0, 4.0, 41),
+)
+
+mcdc.tally.cs_tally(
+    N_cs_bins=[150],
+    cs_bin_size=[5.0, 5.0],
+    x=np.linspace(0.0, 4.0, 41),
+    y=np.linspace(0.0, 4.0, 41),
+    scores=["flux"],
 )
 
 # Setting
-mcdc.setting(N_particle=50)
+mcdc.setting(N_particle=1e3)
 mcdc.implicit_capture()
 
 # Run

examples/fixed_source/cooper2/input.py

@@ -43,8 +43,8 @@
 
 mcdc.tally.mesh_tally(
     scores=["flux"],
-    x=np.linspace(0.0, 4.0, 40),
-    y=np.linspace(0.0, 4.0, 40),
+    x=np.linspace(0.0, 4.0, 41),
+    y=np.linspace(0.0, 4.0, 41),
 )
 
 # Setting

examples/fixed_source/cooper_combo/cs_processing/cs_process.py

@@ -0,0 +1,102 @@
+import h5py
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.fft as spfft
+import time
+from mcdc.main import cs_reconstruct, construct_cs_sampling_matrix_S
+
+try:
+    import cvxpy as cp
+except ImportError:
+    print("CVXPY has not been installed. Please install it with 'pip install cvxpy'")
+
+
+def rel_l2_error(recon, real):
+    recon = recon.flatten()
+    real = real.flatten()
+
+    return np.linalg.norm(real - recon, ord=2) / np.linalg.norm(real, ord=2)
+
+
+Nx = 40
+Ny = 40
+with h5py.File("output.h5", "r") as f:
+    S, N_dims, bin_size_pixels = construct_cs_sampling_matrix_S(f, Nx, Ny)
+    N_cs_bins = f["tallies"]["cs_tally_0"]["N_cs_bins"][()]
+    cs_results = f["tallies"]["cs_tally_0"]["flux"]["mean"][:]
+    mesh_results = f["tallies"]["mesh_tally_0"]["flux"]["mean"][:]
+    mesh_sdev = f["tallies"]["mesh_tally_0"]["flux"]["sdev"][:]
+
+lambda_array = [
+    "mesh",
+    0,
+    5e-7,
+    1e-6,
+    5e-6,
+    1e-5,
+    5e-5,
+    1e-4,
+    5e-4,
+    1e-3,
+    5e-3,
+    1e-2,
+    5e-2,
+    1e-1,
+    2e-1,
+    3e-1,
+]
+recon_array = [
+    mesh_results if lam == "mesh" else cs_reconstruct(S, cs_results, lam, Nx, Ny)
+    for lam in lambda_array
+]
+rel_errors = [rel_l2_error(recon, mesh_results) for recon in recon_array]
+
+
+# Plotting the reconstructions for different lambda values
+fig, axes = plt.subplots(4, 4, figsize=(12, 12))
+for i, ax in enumerate(axes.flat):
+    l = lambda_array[i]
+    reconstruction = recon_array[i]
+    im = ax.imshow(reconstruction.T, origin="lower", extent=[0, 4, 0, 4])
+
+    ax.set_title(f"MC/DC Solution" if l == "mesh" else f"$\lambda$ = {l:.5g}")
+    ax.set_ylabel("y [cm]" if i % 4 == 0 else "")
+    ax.set_xlabel("x [cm]" if i >= 12 else "")
+
+    # Add colorbar
+    cbar = fig.colorbar(im, ax=ax, orientation="vertical", shrink=1)
+    cbar.formatter.set_powerlimits((0, 0))
+
+plt.suptitle(
+    "Basis Pursuit Denoising Reconstructions with Different Values of $\lambda$",
+    fontsize=16,
+)
+plt.tight_layout()
+plt.savefig(
+    f"{N_dims}D_cooper_recons_{N_cs_bins}_{bin_size_pixels}{bin_size_pixels}.png"
+)
+plt.show()
+
+
+# Plotting the relative errors
+plt.plot(lambda_array[1:], rel_errors[1:], label="Reconstruction Errors")
+plt.hlines(
+    np.linalg.norm(mesh_sdev.flatten(), ord=2)
+    / np.linalg.norm(mesh_results.flatten(), ord=2),
+    plt.xlim()[0],
+    plt.xlim()[1],
+    color="black",
+    linestyle="--",
+    label="Reference Std Dev",
+)
+plt.xscale("log")
+plt.yscale("log")
+plt.xlabel("$\lambda$")
+plt.ylabel("Relative L$^2$ Error")
+plt.title("Relative Error vs $\lambda$ - 2D Modified Cooper-Larsen Reconstructions")
+plt.legend()
+plt.tight_layout()
+plt.savefig(
+    f"{N_dims}D_cooper_errors_{N_cs_bins}_{bin_size_pixels}{bin_size_pixels}.png"
+)
+plt.show()

examples/fixed_source/cooper_combo/cs_processing/single_recon_cs_process.py

@@ -0,0 +1,28 @@
+import h5py
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.fft as spfft
+import time
+from mcdc.main import cs_reconstruct, construct_cs_sampling_matrix_S
+
+try:
+    import cvxpy as cp
+except ImportError:
+    print("CVXPY has not been installed. Please install it with 'pip install cvxpy'")
+
+Nx = 40
+Ny = 40
+lambda_ = 1e-4
+with h5py.File("output.h5", "r") as f:
+    S, N_dims, bin_size_pixels = construct_cs_sampling_matrix_S(f, Nx, Ny)
+    N_cs_bins = f["tallies"]["cs_tally_0"]["N_cs_bins"][()]
+    cs_results = f["tallies"]["cs_tally_0"]["flux"]["mean"][:]
+    mesh_results = f["tallies"]["mesh_tally_0"]["flux"]["mean"][:]
+    mesh_sdev = f["tallies"]["mesh_tally_0"]["flux"]["sdev"][:]
+
+start_time = time.time()
+recon = cs_reconstruct(S, cs_results, lambda_, Nx, Ny)
+recon_time = time.time() - start_time
+
+plt.imshow(recon)
+plt.show()

examples/fixed_source/cooper_combo/input.py

@@ -60,8 +60,16 @@
     y=np.linspace(0.0, 4.0, 41),
 )
 
+mcdc.tally.cs_tally(
+    N_cs_bins=[300],
+    cs_bin_size=[3.0, 3.0],
+    scores=["flux"],
+    x=np.linspace(0.0, 4.0, 41),
+    y=np.linspace(0.0, 4.0, 41),
+)
+
 # Setting
-mcdc.setting(N_particle=1e2)
+mcdc.setting(N_particle=1e3)
 mcdc.implicit_capture()
 
 # Run

examples/fixed_source/cooper_combo/process.py

@@ -2,7 +2,6 @@
 import h5py
 import numpy as np
 
-
 # Load result
 with h5py.File("output.h5", "r") as f:
     tally = f["tallies/mesh_tally_0"]

examples/fixed_source/kobayashi3-TD/cs_processing/cs_process.py

@@ -0,0 +1,105 @@
+import h5py
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.fft as spfft
+import time
+from mcdc.main import cs_reconstruct, construct_cs_sampling_matrix_S
+
+try:
+    import cvxpy as cp
+except ImportError:
+    print("CVXPY has not been installed. Please install it with 'pip install cvxpy'")
+
+
+def rel_l2_error(recon, real):
+    recon = recon.flatten()
+    real = real.flatten()
+
+    return np.linalg.norm(real - recon, ord=2) / np.linalg.norm(real, ord=2)
+
+
+Nx = 30
+Ny = 50
+with h5py.File("output.h5", "r") as f:
+    S, N_dims, bin_size_pixels = construct_cs_sampling_matrix_S(
+        f, Nx, Ny, problem_lims=[60, 100, 60]
+    )
+    print(N_dims, bin_size_pixels)
+    N_cs_bins = f["tallies"]["cs_tally_0"]["N_cs_bins"][()]
+    cs_results = f["tallies"]["cs_tally_0"]["flux"]["mean"][:]
+    mesh_results = f["tallies"]["mesh_tally_0"]["flux"]["mean"][:]
+    mesh_sdev = f["tallies"]["mesh_tally_0"]["flux"]["sdev"][:]
+
+lambda_array = [
+    "mesh",
+    0,
+    5e-7,
+    1e-6,
+    5e-6,
+    1e-5,
+    5e-5,
+    1e-4,
+    5e-4,
+    1e-3,
+    5e-3,
+    1e-2,
+    5e-2,
+    1e-1,
+    2e-1,
+    3e-1,
+]
+recon_array = [
+    mesh_results if lam == "mesh" else cs_reconstruct(S, cs_results, lam, Nx, Ny)
+    for lam in lambda_array
+]
+rel_errors = [rel_l2_error(recon, mesh_results) for recon in recon_array]
+
+
+# Plotting the reconstructions for different lambda values
+fig, axes = plt.subplots(4, 4, figsize=(12, 12))
+for i, ax in enumerate(axes.flat):
+    l = lambda_array[i]
+    reconstruction = recon_array[i]
+    im = ax.imshow(reconstruction.T, origin="lower", extent=[0, 60, 0, 100])
+
+    ax.set_title(f"MC/DC Solution" if l == "mesh" else f"$\lambda$ = {l:.5g}")
+    ax.set_ylabel("y [cm]" if i % 4 == 0 else "")
+    ax.set_xlabel("x [cm]" if i >= 12 else "")
+
+    # Add colorbar
+    cbar = fig.colorbar(im, ax=ax, orientation="vertical", shrink=1)
+    cbar.formatter.set_powerlimits((0, 0))
+
+plt.suptitle(
+    "Basis Pursuit Denoising Reconstructions with Different Values of $\lambda$",
+    fontsize=16,
+)
+plt.tight_layout()
+plt.savefig(
+    f"{N_dims}D_kobayashi_recons_{N_cs_bins}_{bin_size_pixels}{bin_size_pixels}.png"
+)
+plt.show()
+
+
+# Plotting the relative errors
+plt.plot(lambda_array[1:], rel_errors[1:], label="Reconstruction Errors")
+plt.hlines(
+    np.linalg.norm(mesh_sdev.flatten(), ord=2)
+    / np.linalg.norm(mesh_results.flatten(), ord=2),
+    plt.xlim()[0],
+    plt.xlim()[1],
+    color="black",
+    linestyle="--",
+    label="Reference Std Dev",
+)
+plt.xscale("log")
+plt.yscale("log")
+plt.xlabel("$\lambda$")
+plt.ylabel("Relative L$^2$ Error")
+plt.title(f"Relative Error vs $\lambda$ - {N_dims}D Kobayashi Reconstructions")
+plt.legend()
+plt.tight_layout()
+plt.savefig(
+    f"{N_dims}D_kobayashi_errors_{N_cs_bins}_{bin_size_pixels}{bin_size_pixels}.png"
+)
+plt.show()

examples/fixed_source/kobayashi3-TD/cs_processing/single_recon_cs_process.py

@@ -0,0 +1,27 @@
+import h5py
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.fft as spfft
+import time
+from mcdc.main import cs_reconstruct, construct_cs_sampling_matrix_S
+
+try:
+    import cvxpy as cp
+except ImportError:
+    print("CVXPY has not been installed. Please install it with 'pip install cvxpy'")
+
+Nx = 30
+Ny = 50
+lambda_ = 1e-4
+with h5py.File("output.h5", "r") as f:
+    S, N_dims, bin_size_pixels = construct_cs_sampling_matrix_S(
+        f, Nx, Ny, problem_lims=[60, 100, 60]
+    )
+    N_cs_bins = f["tallies"]["cs_tally_0"]["N_cs_bins"][()]
+    cs_results = f["tallies"]["cs_tally_0"]["flux"]["mean"][:]
+    mesh_results = f["tallies"]["mesh_tally_0"]["flux"]["mean"][:]
+    mesh_sdev = f["tallies"]["mesh_tally_0"]["flux"]["sdev"][:]
+
+start_time = time.time()
+recon = cs_reconstruct(S, cs_results, lambda_, Nx, Ny)
+recon_time = time.time() - start_time

examples/fixed_source/kobayashi3-TD/input.py

@@ -60,26 +60,16 @@
     scores=["flux"],
     x=np.linspace(0.0, 60.0, 31),
     y=np.linspace(0.0, 100.0, 51),
-    # t=np.linspace(0.0, 200.0, 21),
-    # g=np.array([-0.5, 3.5, 6.5]) # fast (0, 1, 2, 3) and thermal (4, 5, 6) groups
 )
 
-mcdc.tally.cell_tally(source_cell, scores=["flux"])
-mcdc.tally.cell_tally(void_cell, scores=["flux"])
-mcdc.tally.cell_tally(shield_cell, scores=["flux"])
-
-
 mcdc.tally.cs_tally(
-    N_cs_bins=[150],
-    cs_bin_size=[8.0, 8.0],
-    x=np.linspace(0.0, 60.0, 31),
-    y=np.linspace(0.0, 100.0, 51),
+    N_cs_bins=[1000],
+    cs_bin_size=np.array([3.0, 3.0]),
     scores=["flux"],
 )
 
-
 # Setting
-mcdc.setting(N_particle=1e2)
+mcdc.setting(N_particle=1e3)
 
 # Run
 mcdc.run()

examples/fixed_source/kobayashi3-TD/process.py

@@ -4,17 +4,12 @@
 import h5py
 import matplotlib.animation as animation
 
-
 # =============================================================================
 # Plot results
 # =============================================================================
 
 # Results
 with h5py.File("output.h5", "r") as f:
-    cs_recon = f["tallies/cs_tally_0/flux/reconstruction"][:]
-    plt.imshow(cs_recon)
-    plt.show()
-
     tallies = f["tallies/mesh_tally_0"]
     flux = tallies["flux"]
     grid = tallies["grid"]
@@ -29,8 +24,7 @@
     phi = flux["mean"][:]
     phi_sd = flux["sdev"][:]
 
-    for i in range(len(f["input_deck"]["cell_tallies"])):
-        flux_score = f[f"tallies/cell_tally_{i}/flux"]
-        print(
-            f'cell {i+1} mean = {flux_score["mean"][()]}, sdev = {flux_score["sdev"][()]}'
-        )
+    print(len(x))
+
+    plt.imshow(phi[:, :, 0])
+    plt.show()