From e88cec8fd6d80c64cfc0abb43bf856579ad00efd Mon Sep 17 00:00:00 2001 From: Kayla Clements Date: Tue, 14 Nov 2023 12:24:23 -0800 Subject: [PATCH] Merging variance deconvolution feature into main. UQ capability tested with numba and MPI. UQ available for fixed-source problems, materials/nuclides uncertainties. Added UQ regression test. Squashed commit of feature branch work from Sept 29 2023 to Nov 14 2023. --- mcdc/__init__.py | 1 + mcdc/card.py | 25 ++ mcdc/constant.py | 1 + mcdc/input_.py | 78 +++++ mcdc/kernel.py | 288 ++++++++++++++++++ mcdc/loop.py | 15 +- mcdc/main.py | 69 +++++ mcdc/type_.py | 126 ++++++++ .../fixed_source/variance_deconv/__init__.py | 0 .../fixed_source/variance_deconv/answer.h5 | Bin 0 -> 169648 bytes .../fixed_source/variance_deconv/test.py | 61 ++++ 11 files changed, 663 insertions(+), 1 deletion(-) create mode 100644 test/regression/fixed_source/variance_deconv/__init__.py create mode 100644 test/regression/fixed_source/variance_deconv/answer.h5 create mode 100644 test/regression/fixed_source/variance_deconv/test.py diff --git a/mcdc/__init__.py b/mcdc/__init__.py index d12edfcc..bc43b278 100644 --- a/mcdc/__init__.py +++ b/mcdc/__init__.py @@ -20,6 +20,7 @@ weight_roulette, IC_generator, dsm, + uq, print_card, reset_cards, ) diff --git a/mcdc/card.py b/mcdc/card.py index 88cd4056..98bd8a6c 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -31,6 +31,7 @@ def reset(self): "total": False, "current": False, "eddington": False, + "exit": False, "mesh": make_card_mesh(), } @@ -116,6 +117,14 @@ def reset(self): "IC_cycle_stretch": 1.0, "branchless_collision": False, "dsm_order": 1, + "uq": False, + } + + self.uq_deltas = { + "tag": "Uq", + "nuclides": [], + "materials": [], + "surfaces": [], } @@ -152,6 +161,7 @@ def make_card_nuclide(G, J): card["sensitivity"] = False card["sensitivity_ID"] = 0 card["dsm_Np"] = 1.0 + card["uq"] = False return card @@ -176,6 +186,7 @@ def make_card_material(N_nuclide, G, J): card["chi_s"] = np.zeros([G, G]) card["chi_p"] = np.zeros([G, G]) card["sensitivity"] = False + card["uq"] = False return card @@ -284,3 +295,17 @@ def make_card_mesh(): "azi": np.array([-PI, PI]), "g": np.array([-INF, INF]), } + + +def make_card_uq(): + return { + "tag": "t", + "ID": -1, + "key": "k", + "mean": 0.0, + "delta": 0.0, + "distribution": "d", + "rng_seed": 0, + "group": False, + "group_group": False, + } diff --git a/mcdc/constant.py b/mcdc/constant.py index f7d87cab..e85a13e2 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -48,3 +48,4 @@ SEED_SPLIT_SOURCE_PRECURSOR = nb.uint64(0x546F6464) SEED_SPLIT_BANK = nb.uint64(0x5279616E) SEED_SPLIT_PARTICLE = nb.uint64(0) +SEED_SPLIT_UQ = nb.uint(0x5368656261) diff --git a/mcdc/input_.py b/mcdc/input_.py index 250c1d93..1210c432 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -17,6 +17,7 @@ make_card_universe, make_card_lattice, make_card_source, + make_card_uq, ) from mcdc.constant import ( GYRATION_RADIUS_ALL, @@ -1279,6 +1280,83 @@ def dsm(order=1): card["dsm_order"] = order +def uq(**kw): + def append_card(delta_card, global_tag): + delta_card["distribution"] = dist + delta_card["flags"] = [] + for key in kw.keys(): + check_support(parameter["tag"] + " parameter", key, parameter_list, False) + delta_card["flags"].append(key) + delta_card[key] = kw[key] + mcdc.input_deck.uq_deltas[global_tag].append(delta_card) + + mcdc.input_deck.technique["uq"] = True + # Make sure N_batch > 1 + if mcdc.input_deck.setting["N_batch"] <= 1: + print_error("Must set N_batch>1 with mcdc.setting() prior to mcdc.uq() call.") + + # Check uq parameter + parameter_ = check_support( + "uq parameter", + list(kw)[0], + ["nuclide", "material", "surface", "source"], + False, + ) + parameter = kw[parameter_] + del kw[parameter_] + parameter["uq"] = True + + # Confirm supplied distribution + check_requirement("uq", kw, ["distribution"]) + dist = check_support("distribution", kw["distribution"], ["uniform"], False) + del kw["distribution"] + + # Only remaining keywords should be the parameter delta(s) + + if parameter["tag"] == "Material": + parameter_list = [ + "capture", + "scatter", + "fission", + "nu_s", + "nu_p", + "nu_d", + "chi_p", + "chi_d", + "speed", + "decay", + ] + global_tag = "materials" + if parameter["N_nuclide"] == 1: + nuc_card = make_card_nuclide(parameter["G"], parameter["J"]) + nuc_card["ID"] = parameter["nuclide_IDs"][0] + append_card(nuc_card, "nuclides") + delta_card = make_card_material( + parameter["N_nuclide"], parameter["G"], parameter["J"] + ) + for name in ["ID", "nuclide_IDs", "nuclide_densities"]: + delta_card[name] = parameter[name] + elif parameter["tag"] == "Nuclide": + parameter_list = [ + "capture", + "scatter", + "fission", + "nu_s", + "nu_p", + "nu_d", + "chi_p", + "chi_d", + "speed", + "decay", + ] + global_tag = "nuclides" + delta_card = make_card_nuclide(parameter["G"], parameter["J"]) + delta_card["ID"] = parameter["ID"] + append_card(delta_card, global_tag) + # elif parameter['tag'] is 'Surface': + # elif parameter['tag'] is 'Source': + + # ============================================================================== # Util # ============================================================================== diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 58aaf272..795a42c6 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -153,6 +153,16 @@ def rng_from_seed(seed): return rng_(seed) / RNG_MOD +@njit +def rng_array(seed, shape, size): + xi = np.zeros(size) + for i in range(size): + xi_seed = split_seed(i, seed) + xi[i] = rng_from_seed(xi_seed) + xi = xi.reshape(shape) + return xi + + # ============================================================================= # Particle source operations # ============================================================================= @@ -958,6 +968,14 @@ def surface_bc(P, surface, trans): surface_reflect(P, surface, trans) +@njit +def surface_exit_evaluate(P): + if P["surface_ID"] == 0: + return 0 + else: + return 1 + + @njit def surface_reflect(P, surface, trans): ux = P["ux"] @@ -1267,6 +1285,19 @@ def score_tracklength(P, distance, mcdc): score_eddington(s, g, t, x, y, z, flux, P, tally["score"]["eddington"]) +@njit +def score_exit(P, x, mcdc): + tally = mcdc["tally"] + material = mcdc["materials"][P["material_ID"]] + + s = P["sensitivity_ID"] + mu, azi = mesh_get_angular_index(P, tally["mesh"]) + g = mesh_get_energy_index(P, tally["mesh"]) + + flux = P["w"] / abs(P["ux"]) + score_flux(s, g, 0, x, 0, 0, mu, azi, flux, tally["score"]["exit"]) + + @njit def score_flux(s, g, t, x, y, z, mu, azi, flux, score): score["bin"][s, g, t, x, y, z, mu, azi] += flux @@ -1822,6 +1853,13 @@ def surface_crossing(P, mcdc): # Record old material for sensitivity quantification material_ID_old = P["material_ID"] + # Tally particle exit + if mcdc["tally"]["exit"] and not P["alive"]: + # Reflectance if P["surface_ID"] == 0, else transmittance + exit_idx = surface_exit_evaluate(P) + # Score on tally + score_exit(P, exit_idx, mcdc) + # Check new cell? if P["alive"] and not surface["reflective"]: cell = mcdc["cells"][P["cell_ID"]] @@ -3407,3 +3445,253 @@ def binary_search(val, grid): else: right = mid - 1 return int(right) + + +# ============================================================================= +# Variance Deconvolution +# ============================================================================= + + +@njit +def uq_resample(mean, delta, info): + # Currently only uniform distribution + shape = mean.shape + size = mean.size + xi = rng_array(info["rng_seed"], shape, size) + + return mean + (2 * xi - 1) * delta + + +@njit +def reset_material(mcdc, idm, material_uq): + # Assumes all nuclides have already been re-sampled + # Basic XS + material = mcdc["materials"][idm] + for tag in literal_unroll(("capture", "scatter", "fission", "total")): + if material_uq["flags"][tag]: + material[tag][:] = 0.0 + for n in range(material["N_nuclide"]): + nuc1 = mcdc["nuclides"][material["nuclide_IDs"][n]] + density = material["nuclide_densities"][n] + material[tag] += nuc1[tag] * density + + # Effective speed + if material_uq["flags"]["speed"]: + material["speed"][:] = 0.0 + for n in range(material["N_nuclide"]): + nuc2 = mcdc["nuclides"][material["nuclide_IDs"][n]] + density = material["nuclide_densities"][n] + material["speed"] += nuc2["speed"] * nuc2["total"] * density + if max(material["total"]) == 0.0: + material["speed"][:] = nuc2["speed"][:] + else: + material["speed"] /= material["total"] + + # Calculate effective spectra and multiplicities of scattering and prompt fission + G = material["G"] + if max(material["scatter"]) > 0.0: + shape = material["chi_s"].shape + nuSigmaS = np.zeros(shape) + for i in range(material["N_nuclide"]): + nuc3 = mcdc["nuclides"][material["nuclide_IDs"][i]] + density = material["nuclide_densities"][i] + SigmaS = np.diag(nuc3["scatter"]) * density + nu_s = np.diag(nuc3["nu_s"]) + chi_s = np.ascontiguousarray(nuc3["chi_s"].transpose()) + nuSigmaS += chi_s.dot(nu_s.dot(SigmaS)) + chi_nu_s = nuSigmaS.dot(np.diag(1.0 / material["scatter"])) + material["nu_s"] = np.sum(chi_nu_s, axis=0) + material["chi_s"] = np.ascontiguousarray( + chi_nu_s.dot(np.diag(1.0 / material["nu_s"])).transpose() + ) + if max(material["fission"]) > 0.0: + nuSigmaF = np.zeros((G, G), dtype=float) + for n in range(material["N_nuclide"]): + nuc4 = mcdc["nuclides"][material["nuclide_IDs"][n]] + density = material["nuclide_densities"][n] + SigmaF = np.diag(nuc4["fission"]) * density + nu_p = np.diag(nuc4["nu_p"]) + chi_p = np.ascontiguousarray(np.transpose(nuc4["chi_p"])) + nuSigmaF += chi_p.dot(nu_p.dot(SigmaF)) + chi_nu_p = nuSigmaF.dot(np.diag(1.0 / material["fission"])) + material["nu_p"] = np.sum(chi_nu_p, axis=0) + # Required because the below function otherwise returns an F-contiguous array + material["chi_p"] = np.ascontiguousarray( + np.transpose(chi_nu_p.dot(np.diag(1.0 / material["nu_p"]))) + ) + + # Calculate delayed and total fission multiplicities + if max(material["fission"]) > 0.0: + material["nu_f"][:] = material["nu_p"][:] + for j in range(material["J"]): + total = np.zeros(material["G"]) + for n in range(material["N_nuclide"]): + nuc5 = mcdc["nuclides"][material["nuclide_IDs"][n]] + density = material["nuclide_densities"][n] + total += nuc5["nu_d"][:, j] * nuc5["fission"] * density + material["nu_d"][:, j] = total / material["fission"] + material["nu_f"] += material["nu_d"][:, j] + + +@njit +def reset_nuclide(nuclide, nuclide_uq): + for name in literal_unroll( + ("speed", "decay", "capture", "fission", "nu_s", "nu_p") + ): + if nuclide_uq["flags"][name]: + nuclide[name] = uq_resample( + nuclide_uq["mean"][name], nuclide_uq["delta"][name], nuclide_uq["info"] + ) + + if nuclide_uq["flags"]["scatter"]: + scatter = uq_resample( + nuclide_uq["mean"]["scatter"], + nuclide_uq["delta"]["scatter"], + nuclide_uq["info"], + ) + nuclide["scatter"] = np.sum(scatter, 0) + nuclide["chi_s"][:, :] = np.swapaxes(scatter, 0, 1)[:, :] + for g in range(nuclide["G"]): + if nuclide["scatter"][g] > 0.0: + nuclide["chi_s"][g, :] /= nuclide["scatter"][g] + + if nuclide_uq["flags"]["total"]: + nuclide["total"][:] = ( + nuclide["capture"] + nuclide["scatter"] + nuclide["fission"] + ) + + if nuclide_uq["flags"]["nu_d"]: + nu_d = uq_resample( + nuclide_uq["mean"]["nu_d"], nuclide_uq["delta"]["nu_d"], nuclide_uq["info"] + ) + nuclide["nu_d"][:, :] = np.swapaxes(nu_d, 0, 1)[:, :] + + if nuclide_uq["flags"]["nu_f"]: # True if either nu_p or nu_d is true + nuclide["nu_f"] = nuclide["nu_p"] + for j in range(nuclide["J"]): + nuclide["nu_f"] += nuclide["nu_d"][:, j] + + # Prompt fission spectrum (If G == 1, all ones) + if nuclide_uq["flags"]["chi_p"]: + chi_p = uq_resample( + nuclide_uq["mean"]["chi_p"], + nuclide_uq["delta"]["chi_p"], + nuclide_uq["info"], + ) + nuclide["chi_p"][:, :] = np.swapaxes(chi_p, 0, 1)[:, :] + # Normalize + for g in range(nuclide["G"]): + if np.sum(nuclide["chi_p"][g, :]) > 0.0: + nuclide["chi_p"][g, :] /= np.sum(nuclide["chi_p"][g, :]) + + # Delayed fission spectrum (matrix of size JxG) + if nuclide_uq["flags"]["chi_d"]: + chi_d = uq_resample( + nuclide_uq["mean"]["chi_d"], + nuclide_uq["delta"]["chi_d"], + nuclide_uq["info"], + ) + # Transpose: [gout, dg] -> [dg, gout] + nuclide["chi_d"][:, :] = np.swapaxes(chi_d, 0, 1)[:, :] + # Normalize + for dg in range(nuclide["J"]): + if np.sum(nuclide["chi_d"][dg, :]) > 0.0: + nuclide["chi_d"][dg, :] /= np.sum(nuclide["chi_d"][dg, :]) + + +@njit +def uq_reset(mcdc, seed): + # Types of uq parameters: materials, nuclides + N = len(mcdc["technique"]["uq_"]["nuclides"]) + for i in range(N): + mcdc["technique"]["uq_"]["nuclides"][i]["info"]["rng_seed"] = split_seed( + i, seed + ) + idn = mcdc["technique"]["uq_"]["nuclides"][i]["info"]["ID"] + reset_nuclide(mcdc["nuclides"][idn], mcdc["technique"]["uq_"]["nuclides"][i]) + + M = len(mcdc["technique"]["uq_"]["materials"]) + for i in range(M): + mcdc["technique"]["uq_"]["materials"][i]["info"]["rng_seed"] = split_seed( + i, seed + ) + idm = mcdc["technique"]["uq_"]["materials"][i]["info"]["ID"] + reset_material(mcdc, idm, mcdc["technique"]["uq_"]["materials"][i]) + + +@njit +def uq_tally_closeout_history(mcdc): + tally = mcdc["tally"] + uq_tally = mcdc["technique"]["uq_tally"] + + for name in literal_unroll(score_list): + if uq_tally[name]: + uq_score_closeout_history(tally["score"][name], uq_tally["score"][name]) + + +@njit +def uq_score_closeout_history(score, uq_score): + # Assumes N_batch > 1 + # Accumulate square of history score, but continue to accumulate bin + history_bin = score["bin"] - uq_score["batch_bin"] + uq_score["batch_var"][:] += history_bin**2 + uq_score["batch_bin"] = score["bin"] + + +@njit +def uq_tally_closeout_batch(mcdc): + uq_tally = mcdc["technique"]["uq_tally"] + + for name in literal_unroll(score_list): + if uq_tally[name]: + # Reset bin + uq_tally["score"][name]["batch_bin"].fill(0.0) + uq_reduce_bin(uq_tally["score"][name]) + + +@njit +def uq_reduce_bin(score): + # MPI Reduce + buff = np.zeros_like(score["batch_var"]) + with objmode(): + MPI.COMM_WORLD.Reduce(np.array(score["batch_var"]), buff, MPI.SUM, 0) + score["batch_var"][:] = buff + + +@njit +def uq_tally_closeout(mcdc): + tally = mcdc["tally"] + uq_tally = mcdc["technique"]["uq_tally"] + + for name in literal_unroll(score_list): + # Uq_tally implies tally, but tally does not imply uq_tally + if uq_tally[name]: + uq_score_closeout(name, mcdc) + elif tally[name]: + score_closeout(tally["score"][name], mcdc) + + +@njit +def uq_score_closeout(name, mcdc): + score = mcdc["tally"]["score"][name] + uq_score = mcdc["technique"]["uq_tally"]["score"][name] + + N_history = mcdc["setting"]["N_particle"] + + # At this point, score["sdev"] is still just the sum of the squared mean from every batch + uq_score["batch_var"] = (uq_score["batch_var"] / N_history - score["sdev"]) / ( + N_history - 1 + ) + + # If we're here, N_batch > 1 + N_history = mcdc["setting"]["N_batch"] + + # Store results + score["mean"][:] = score["mean"] / N_history + uq_score["batch_var"] /= N_history + uq_score["batch_bin"] = (score["sdev"] - N_history * np.square(score["mean"])) / ( + N_history - 1 + ) + score["sdev"][:] = np.sqrt( + (score["sdev"] / N_history - np.square(score["mean"])) / (N_history - 1) + ) diff --git a/mcdc/loop.py b/mcdc/loop.py index b30528a3..8762098c 100644 --- a/mcdc/loop.py +++ b/mcdc/loop.py @@ -35,6 +35,9 @@ def loop_fixed_source(mcdc): if mcdc["setting"]["N_batch"] > 1: with objmode(): print_header_batch(mcdc) + if mcdc["technique"]["uq"]: + seed_uq = kernel.split_seed(seed_batch, SEED_SPLIT_UQ) + kernel.uq_reset(mcdc, seed_uq) # Loop over time censuses for idx_census in range(mcdc["setting"]["N_census"]): @@ -65,9 +68,15 @@ def loop_fixed_source(mcdc): # Tally history closeout kernel.tally_reduce_bin(mcdc) kernel.tally_closeout_history(mcdc) + # Uq closeout + if mcdc["technique"]["uq"]: + kernel.uq_tally_closeout_batch(mcdc) # Tally closeout - kernel.tally_closeout(mcdc) + if mcdc["technique"]["uq"]: + kernel.uq_tally_closeout(mcdc) + else: + kernel.tally_closeout(mcdc) # ========================================================================= @@ -180,6 +189,10 @@ def loop_source(seed, mcdc): if not mcdc["setting"]["mode_eigenvalue"] and mcdc["setting"]["N_batch"] == 1: kernel.tally_closeout_history(mcdc) + # Tally history closeout for multi-batch uq simulation + if mcdc["technique"]["uq"]: + kernel.uq_tally_closeout_history(mcdc) + # Progress printout percent = (idx_work + 1.0) / mcdc["mpi_work_size"] if mcdc["setting"]["progress_bar"] and int(percent * 100.0) > N_prog: diff --git a/mcdc/main.py b/mcdc/main.py index c8443312..00d3b9ee 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -168,6 +168,8 @@ def prepare(): type_.make_type_lattice(input_deck.lattices) type_.make_type_source(G) type_.make_type_tally(N_tally_scores, input_deck.tally) + type_.make_type_uq_tally(N_tally_scores, input_deck.tally) + type_.make_type_uq(input_deck.uq_deltas, G, J) type_.make_type_setting(input_deck) type_.make_type_technique(N_particle, G, input_deck.technique) type_.make_type_global(input_deck) @@ -306,6 +308,7 @@ def prepare(): "iQMC", "IC_generator", "branchless_collision", + "uq", ]: mcdc["technique"][name] = input_deck.technique[name] @@ -405,6 +408,61 @@ def prepare(): # Threshold mcdc["technique"]["dsm_order"] = input_deck.technique["dsm_order"] + # ========================================================================= + # Variance Deconvolution - UQ + # ========================================================================= + if mcdc["technique"]["uq"]: + # Assumes that all tallies will also be uq tallies + for name in type_.uq_tally.names: + if name != "score": + mcdc["technique"]["uq_tally"][name] = input_deck.tally[name] + + M = len(input_deck.uq_deltas["materials"]) + for i in range(M): + idm = input_deck.uq_deltas["materials"][i]["ID"] + mcdc["technique"]["uq_"]["materials"][i]["info"]["ID"] = idm + mcdc["technique"]["uq_"]["materials"][i]["info"][ + "distribution" + ] = input_deck.uq_deltas["materials"][i]["distribution"] + for name in input_deck.uq_deltas["materials"][i]["flags"]: + mcdc["technique"]["uq_"]["materials"][i]["flags"][name] = True + mcdc["technique"]["uq_"]["materials"][i]["delta"][ + name + ] = input_deck.uq_deltas["materials"][i][name] + flags = mcdc["technique"]["uq_"]["materials"][i]["flags"] + if flags["capture"] or flags["scatter"] or flags["fission"]: + flags["total"] = True + flags["speed"] = True + if flags["nu_p"] or flags["nu_d"]: + flags["nu_f"] = True + if mcdc["materials"][idm]["N_nuclide"] > 1: + for name in type_.uq_mat.names: + mcdc["technique"]["uq_"]["materials"][i]["mean"][ + name + ] = input_deck.materials[idm][name] + + N = len(input_deck.uq_deltas["nuclides"]) + for i in range(N): + mcdc["technique"]["uq_"]["nuclides"][i]["info"][ + "distribution" + ] = input_deck.uq_deltas["nuclides"][i]["distribution"] + idn = input_deck.uq_deltas["nuclides"][i]["ID"] + mcdc["technique"]["uq_"]["nuclides"][i]["info"]["ID"] = idn + for name in type_.uq_nuc.names: + mcdc["technique"]["uq_"]["nuclides"][i]["mean"][ + name + ] = input_deck.nuclides[idn][name] + for name in input_deck.uq_deltas["nuclides"][i]["flags"]: + mcdc["technique"]["uq_"]["nuclides"][i]["flags"][name] = True + mcdc["technique"]["uq_"]["nuclides"][i]["delta"][ + name + ] = input_deck.uq_deltas["nuclides"][i][name] + flags = mcdc["technique"]["uq_"]["nuclides"][i]["flags"] + if flags["capture"] or flags["scatter"] or flags["fission"]: + flags["total"] = True + if flags["nu_p"] or flags["nu_d"]: + flags["nu_f"] = True + # ========================================================================= # MPI # ========================================================================= @@ -578,6 +636,17 @@ def generate_hdf5(mcdc): "tally/" + name_h5 + "/sdev", data=np.squeeze(T["score"][name]["sdev"]), ) + if mcdc["technique"]["uq_tally"][name]: + mc_var = mcdc["technique"]["uq_tally"]["score"][name][ + "batch_var" + ] + tot_var = mcdc["technique"]["uq_tally"]["score"][name][ + "batch_bin" + ] + f.create_dataset( + "tally/" + name_h5 + "/uq_var", + data=np.squeeze(tot_var - mc_var), + ) # Eigenvalues if mcdc["setting"]["mode_eigenvalue"]: diff --git a/mcdc/type_.py b/mcdc/type_.py index 35a83d7b..2b6a851a 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -143,6 +143,7 @@ def make_type_nuclide(G, J): ("sensitivity", bool_), ("sensitivity_ID", int64), ("dsm_Np", float64), + ("uq", bool_), ] ) @@ -169,6 +170,7 @@ def make_type_material(G, J, Nmax_nuclide): ("chi_s", float64, (G, G)), ("chi_p", float64, (G, G)), ("sensitivity", bool_), + ("uq", bool_), ] ) @@ -326,6 +328,7 @@ def make_type_source(G): "total", "current", "eddington", + "exit", ) @@ -356,6 +359,7 @@ def make_type_score(shape): ["total", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], ["current", (Ns, Ng, Nt, Nx, Ny, Nz, 3)], ["eddington", (Ns, Ng, Nt, Nx, Ny, Nz, 6)], + ["exit", (Ns, Ng, Nt, 2, Ny, Nz, Nmu, N_azi)], ] # Add score flags to structure @@ -443,6 +447,7 @@ def make_type_technique(N_particle, G, card): ("iQMC", bool_), ("IC_generator", bool_), ("branchless_collision", bool_), + ("uq", bool_), ] # ========================================================================= @@ -569,10 +574,131 @@ def make_type_technique(N_particle, G, card): ("dsm_order", int64), ] + # ========================================================================= + # Variance Deconvolution + # ========================================================================= + struct += [("uq_tally", uq_tally), ("uq_", uq)] + # Finalize technique type technique = np.dtype(struct) +# UQ +def make_type_uq_tally(Ns, tally_card): + global uq_tally + + def make_type_uq_score(shape): + return np.dtype( + [ + ("batch_bin", float64, shape), + ("batch_var", float64, shape), + ] + ) + + # Tally estimator flags + struct = [] + + # Mesh, but doesn't need to be added + mesh, Nx, Ny, Nz, Nt, Nmu, N_azi, Ng = make_type_mesh(tally_card["mesh"]) + + # Scores and shapes + scores_shapes = [ + ["flux", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["density", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["fission", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["total", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["current", (Ns, Ng, Nt, Nx, Ny, Nz, 3)], + ["eddington", (Ns, Ng, Nt, Nx, Ny, Nz, 6)], + ["exit", (Ns, Ng, Nt, 2, Ny, Nz, Nmu, N_azi)], + ] + + # Add score flags to structure + for i in range(len(scores_shapes)): + name = scores_shapes[i][0] + struct += [(name, bool_)] + + # Add scores to structure + scores_struct = [] + for i in range(len(scores_shapes)): + name = scores_shapes[i][0] + shape = scores_shapes[i][1] + if not tally_card[name]: + shape = (0,) * len(shape) + scores_struct += [(name, make_type_uq_score(shape))] + scores = np.dtype(scores_struct) + struct += [("score", scores)] + + # Make tally structure + uq_tally = np.dtype(struct) + + +def make_type_uq(uq_deck, G, J): + global uq, uq_nuc, uq_mat + + # def make_type_parameter(shape): + # return np.dtype( + # [ + # ("tag", str_), # nuclides, materials, surfaces, sources + # ("ID", int64), + # ("key", str_), + # ("mean", float64, shape), + # ("delta", float64, shape), + # ("distribution", str_), + # ("rng_seed", uint64), + # ] + # ) + + def make_type_parameter(G, J, decay=False): + # Fields are things that can have deltas + struct = [ + ("speed", float64, (G,)), + ("capture", float64, (G,)), + ("scatter", float64, (G, G)), + ("fission", float64, (G,)), + ("nu_s", float64, (G,)), + ("nu_p", float64, (G,)), + ("nu_d", float64, (G, J)), + ("chi_p", float64, (G, G)), + ] + struct += [("decay", float64, (J,)), ("chi_d", float64, (J, G))] + return np.dtype(struct) + + uq_nuc = make_type_parameter(G, J, True) + uq_mat = make_type_parameter(G, J) + + flags = np.dtype( + [ + ("speed", bool_), + ("decay", bool_), + ("total", bool_), + ("capture", bool_), + ("scatter", bool_), + ("fission", bool_), + ("nu_s", bool_), + ("nu_f", bool_), + ("nu_p", bool_), + ("nu_d", bool_), + ("chi_s", bool_), + ("chi_p", bool_), + ("chi_d", bool_), + ] + ) + info = np.dtype([("distribution", str_), ("ID", int64), ("rng_seed", uint64)]) + + container = np.dtype( + [("mean", uq_nuc), ("delta", uq_mat), ("flags", flags), ("info", info)] + ) + + N_nuclide = len(uq_deck["nuclides"]) + N_material = len(uq_deck["materials"]) + uq = np.dtype( + [("nuclides", container, (N_nuclide,)), ("materials", container, (N_material,))] + ) + + +param_names = ["tag", "ID", "key", "mean", "delta", "distribution", "rng_seed"] + + # ============================================================================== # Global # ============================================================================== diff --git a/test/regression/fixed_source/variance_deconv/__init__.py b/test/regression/fixed_source/variance_deconv/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/test/regression/fixed_source/variance_deconv/answer.h5 b/test/regression/fixed_source/variance_deconv/answer.h5 new file mode 100644 index 0000000000000000000000000000000000000000..4a0b39085bf1c3f1aacdcc6896951d08714bb8ff GIT binary patch literal 169648 zcmeHw4V)uInRoB(2VoHih#EPMi>Mem#SsuADjSv$1IQAWuPA7e&15%`OeW-GcNXzw zL{0;Tn&U;y`9;6Hp3eDl8a?oGr=BBlI^Lm<9_py4?4fXTdg_RZquvWw-OuwR)0I^9 zq_>ml&UW^<^VG+`tE>L?)KgVm)zvpg&prRht;cUww4cwXY*PHXlm0x3zvg?#^dRz6 zd=ecSP}We+0$kak|5KCzh986J^SH})o_C(Yu#`SYyi7T;8C;YOg?sYQU?6(l*)gr8 zaV%~K{?#_cLdyTBzhEGfczaH;&Rwi1atQ5w6UlUV{FP z3H5jBaKIPohQEfKM9A;Sk!FDOw~3`B`{z?Z4@7#_!uZ8I<+$L-QJ?MtNn~ESs?Y1_V)W3B~Ji5OMXJpkWcNh5p62Lpa%C-*b()SONUSB1p=I+oh;}VuwM=irBu* z`myCB0Z;bJ?YwyT?~ehp1JsX&@$x$A52aS9rt`UpOoI9;^1lQD`8 za5n>F2a>4Hb6A|_c5sOFEXZa_+Hbnuel=wG{AzAOhbav*5a3n#@xVChEww6Mj{|sJ z1NEA4d_Jf@o8xezd|Uw?q) zXAcPA89_uZQs?>7nEn{1^L!~Df-=(+1!6!Yhohzad?Y^ZFV zgTEa~&uNbdZfm!n(~nS2e-a-+irZ09z*9}Xk5m41w-+NHd_Jx7!SC3jl-MoD}cLzlOG^#SeoU)zk~iee03?uzH}oS$8wJ>jO6 z2yJzNazN5}#qze(^xu75cQS?_Y)&+DLTerJocu_T6AHJI6W=OMXyrlUL6j3wxUqJ$ z6yRqJUqcQztj>|PMnHGuC8WQ+jl%JUC~T}fLyojn6~4B}zqTAvZFTYs`b6Pt;b^rc z)4Gk8UT-cscM`co_w0m!_VeHQ?nBpq|MM(`U_u?+PIp4t0Sa!hylr)!wK}F01Q(bF z3vW3`m)e_WZ9U6a({!J;OHrPWo1--B&9k-!N$aIn`4;6lU}K&&%R`o#)<-o)c)qC& zBsGl5Z*|hz5-VLwNq74x#3W%wUy-F;Ym1JQBriKPKryCCuQr0T(5qcUNl^lnS4vnj zPn0EEC?%||Qlu54n-E^J*ZlqnI^K!+MA z^{3>w-mi+D=a9vsI}BHX9X?En;BDh%7zNp=WB5$htlD^iir#v}7|#1pX_STM@-arg z54~ZN@mUOU#nX+~^S%}a_|N5Itoes$cV%G;{$B&uZ3yryb$M8&OiG#m!8MMF57vzCiDpLdet1ZD@s|vg8iwH z!x7S)kABy?zlh~##2J37!RJ(uIUn$zN&oAYUlH)qa4Hpsd{t6=bpHi79(py#UWl?> z%w#5*o1FxZ(%%11r;3$oDWjCrsS117K-rxum)YwEN@a?j5WwOK)kImllu6xGF{xdy z`FCqFoz3Y0O=%iBbGGO>eQP`y@I@2X286?SbY!ho7g`596=4v z?<_CB?{Ah2vLKsOy2b?p00sw9wl0$6!Wg}@LymKLthrwGI}Umn^~uK-UFH1F^xDtI zPdq+YRd=3;J~m=}1aw8Vf8XEl1w7Rp)o{Aoi;)jLv-Urqa2Fpm3MFZskFgg3w?2;I zLzR!M9v@5-I(GJ0SMWjW?{U;(8drsJenZMViRUlSK}Bspt}D0wr@|G~zJinhCq?V^ z!H)tkTIVi6FTXxNe|2vhh0dU(brps5meUwK`}*GT6grNM#fvGVw`^i?B-J~fLZj&L zy_Q0TBhC5^sy~v}?d0*UXa~dN<3x`Srs5qtd9q9S2g#Y^k$>`VbaG%~PQ!3TL+w|?}tJ3J5Mx6a3QSMou0Q~7znbv~SI zhdb-@kt?ahy=4)D8?WmfPoYV4j9)<^y=5GOgO~S?cS~q@JA2t`Jy9J|Ub9S?^I>&;{!#%i4YMpcBHmk{A51lv>ziHc^H;gT9+|Dr zFT*yG_PN(B+=mgx>Y2YD9MoS$_JO}OnCkOV9Wn#v?ZxW){MFb#{9K8gZJ+t|`RJRt z>~mf!T{y=@aNZ#L6R@FQg7^zvl22qyYsn`YG>3h@zQ!fb-#cg z8W~ur-GjfTTR(c+9i9jBTj%2eSMouOQ~7znbw1p&-aWpTO59u0FtcW&)=9iPSnZt- zg~rk0uTY3PWp_J+o`dwrzV+@0v0YKh#8-_Rj*#Yj^!uFmTLOEcI-U#GN z1iUoNvgC+(Z@qgk)nKk~cCB|m;tG3Y_I>dznok%V4W>@KRGrI!*?x$*TOjXQEJFA_ z;84-ux)U5>^e=n>9P{Y+gu5A__qvie>{`I$gzo@{=>4na{ScqU2N)t;{b9f@egYFl z0ne*9vgcj!3!#Cc(I5XEfdJy__`Vt8`1imt8O2|N3kb*O7n5_nFz*ko?q_)v`#(R& zk+a8NpjVrf80LBf^^+?^+e@YM-P#W7>4;;f@3x-N?LgAgz&qPPJsnvL4c*=|dffxf z^{U@_&f|jdyt?jXT}uY+=i`YH;{!Bbynmx#ec5s{)QgX+9+L*mp@-$|)(#% z67Yf@-mnU6av#ySVqW`}2>XfRXTJ~eqVdFm0`%xhSAZ1Z0yjXs?@Ag0K)Ry8-Sr20 zM^|rD`qkg33+hP~Wx-M5_VaPVi19I~_I;9oCugD3vo(5vC`ahAsiQ~J*SoAW(B41DbH z_%QHsg)8`=eW67(=riv{hd{I_{3GP(Et41=_}{(bDHKN0iqgE9LVC*?3|6Tl>kU)p zG2H)G0>FTRyX|mDztf?*&#jw}0Ptb^ca`7higre7$FKYDNEE8+d~PCx0Q9l5z+poj z`6uAq@C)rlnzRpYQ`6~8b&w&QjP;HQDCW3KCG@ss{&qXj^c&Zo5^lQ zF?TE6zO-_+v^$l~(DCI8`SJMjwjYl#Z~ITf^k_6oc3iYR-TX0Ju<&San!4)`>=dQn z|9bSSpq|v;8E{yr{d_b>jE_OJ@0NfkXSMHc{TLJ@2cu|3iBN|pfX`s~_}>Wt7(9nhxA)` z>(RjU24zs=zhzgnGtxJ2J(`%XSFW}F(R`hk*4Lv~1o}rKb9>9V*3{OI59LoGk7U^? z$Lq@nzaAa?d&*aDS;OE6b$GpD3RrnXY276NU~m}a{NKoNPLJ)@E2w8`Pxh@x`%sT4 zWnzYEhWzU*hsVd^uHeILJ^H|xp}w+MZL{A1hl>8tL*StH5?KWQJjQ#%-3*ww7pv>h zTd{r6lYZi(_pq zDJZu(tSUaJ`*iDvPfw3V8?mzc0uuC@Kq{FIm0*LQCb=wEweZf`-?nA-aB zahoIg;MaG@zauBIhS9-q%W+`k+4*??BFBX>I{tth=k(Zaz4GZnn%H#|K=&huQh|(f@$@qV}-9wddV^*bJDh|EvEAIf`M%0>1`F5&g~IgQJ1| zumbiWfPPQ7n*n<7kIuJW#NyP)pq@qV=Z*UyKK>Y7Ai^!b0C3ZYBerOK0KJ)_#PGdu zK(q6(xcC1^-`NXW20p&#Ni{mJ2<`ip)2)jaC|)K@xdBHryu>f$A^KB z$47(@cp}8_3m8K~NWC5Tt*=WyfqeNn6FbfMhwFuYX$$*06a1EbBELW7R{fnESHo!E z|Ci$yFgkb`G2I(q$Kb#dz2hmgfsW9V6fzvyZHEK;jp=;_(jSu6_2{SjY z`PWx=^WoDoM77-ie)JQt+>lPEbZW=xV3#=`c0C_D33(!>i4?c#ht>K2^96cBmVx+i zmIOQOQ&h~>41Nbfl*?I+Sl-Cjv+;@`dVN1LhV#BA1#Fh5H-|X?7CCMPaVvj7hd@wk zt}h7y7!1}_QC5B`$CmE2zZz{c5l!FH7A zpgb4lc_`0EN#EIsqW=Pv7ox17+>t2fbLotd&lNJMlA;u*A)>G!<#izW?`f>^vu}PV z%0ziGv7@NqDzl>0F#YusgVVQ{xl-TE`LL^>+cqNoT)*pN5rJN)>(vDMy05u@Shdel z0WZ$(_s#a-B#(f8Hij%OJ`McL1nR%yv%pWKTzNM5slP;+CtS+FTHDJr*?7^<%BcMZ z4SP=|8gO^(tgKOWnxU&BE{Tsn&5aYcmPGazQK#t2IF7iw{u86qU@p2r( z94iK=SGT!RuLj4*`#nBbE9ms2cY1sn`1tsU@Db8&-Qf3cd>Z-UpO;LptF7htZ_GYL z-ak-olYb<~(QI`71UYU2OMf6F$JG(H_$P8)193}1InHjwc?ar39kanxQHUn<98 zazhv9dev{gbR6mvjp(UJDmff2&H1ov-102siPwvoqxcyfAAdX|d{`Z~oFHg7&bGtm z`eBuilRQ3{CFra>o#OFf;NuLB4+9@(xsnflKg4*L?+=aLVG_NOm&91%W3x{(|Feu&-37ytZWdR>hy zzaJudioAcI+!kLc$I;ZYd(i(JW%+Y&p(u(!1?&K);3k`yuuU z^oSU@yJXhV>*m7;d|-0@?}vE1z}~3VGC1U`%=xhEdCR+yCtm;ErXN=4)9)7O30a0A zSNs}p$M1)Toy*ssF`-%X&YdsE1rb*pm*bjDwKy)HicJ?CR$}HX{98zw1u30=-ae z*6hI~XRaSs?Q>JWi?jQEv;7d`JD|UdVRg*L!CypwBmsUZ<@^NrslOadEqbXAi-EPa zm;U!dw6MMSxuTFJalsDp^7|oXcY?l)vKPjP)}?Y>1ab41$#GG{EnY6iQG0D%DaT=Q zLl@?H)o;HvFR;Jr(934!hUR?OHEvlzo_M`*+x*-vKE5#`d{`Z~d{fYFoNb5A+ly5` zzUT453|ME~>0yr#10O$eB_DKNMH4%x8I;PKp&t;%FW(IDqVJJpe+lumU&(Qc$Ky%b zPXLZ?FOq|6-;cSn{`mbH%deC3m?UDalH;<7^It8;(R?rX204yy>r^>7z52Z?^=fc@ zjGY2HKu_37BM#MZ#a53GrUIRD#jzv8M@U!f!SCNV9{J**mrPT`S{KRh-ZIstfkk@ z2lOLUEB)`^I8$J6RBIU=@>S-1*!4VQJMu(KQ$gLPA6Do8=LqzKECcc3EXnWR@bBj9 z&zRp?^v-1ExFF&hua)DPh^y!1xH-hBZ<6C?5LZjeabOh`C3vwM2iBB*epT<>(1WAW*Az|LI zkv*4#Ukx9W&?xyp0bEO0fb=7-$>9Pwpx)qq6S6mr$KV1=kNCg*yw)_$q~eNuVs zmxy_y%Gv!A=+)--cZ$^tX#nNlC`n(MRXHw>mo4v?aUT*QIbEemVN)XTH2jr zUjT+5`q(#p%X&a@s+`llLY&y0Pwinrlc`FkluPCHkNz@5K2@pc-wy^!!4#fIXA10s z&UbLe@jhn4etrJfi0QM@-fcueJJRa&kxwE|kbbA0oTbk&?%bqAF@p!D`S#bzx#kL^ zcVV9#7ed_JUOCQ(xM)d^LzrWQ*>i*FjZnV!`1w}wQ~4cRFfU`aFko(n!OrwQ-~BS$ z89$fL2IqH{m!I#}-zIPWWJT4Rhs;-jaYl(?T20Z zu;Im^sAQQ)ZENIklst6jyJJ{S{77AORQ-rDV`lF<^IS+_qhfWwyA^q&ks9U5S^CV+ zcNgBlw?E1^$S`{6-Y&<55ZC-OIgYd_{8l;6m<#9m%;vkzzku@5_&M}J@KgDnS0XP* zb}}G(KQjK8pa-J&BU^Vt{Ng9%xY#ElzW6DSl=p6X2O38xvTyy*vrk0(=I(j+%5B61 zF=muxy?wXKadEsX+K}Uti1WQ)jw3q}_<$T|%9yt!IiK3`=InJ01=lRUu zSBU*zC?A#Ac_r|2G)M+$o*l=VkX2~jVbuF-fo(o2l{1NP08(t^ATw8}rt`Up3~X{@ zY0Kap1{PgP54mTV=J&+ljWa)#}XD(q1znRK?0yS7UE ztf(Ah7evJ|Ye%zx<+URbcZWha-;S`EC@Xi8h5-xR<1jgO(oY%9pU+s`U-K3mZ~3|GG&sMryfprzTrYo=G^%$#sozj$ z<+vbT?*D5!j>aR6kH~SobK{&)bG_>KIoCT!Y@P!(O|{(r{+joU7$4A*=(%6_ysuw< z`H*0~S9eY$S(OtQXDPz>XIDWz=OM@dKhC0lB8lsMLN(pogSd-l*8Mp8QqhiBcxbsFA}F(p*36;^S)`AFReZ`%k_;B79ig2fiS%BMnEi z3sEM_+ly5`7Ck-;?c;ze`Jiz)jlZ*~NR1yuKS1#`zKCFam>Z8wW4N}asr{}4p;S;# zR>6_FE~h*8U^w+Q0rZ=RZheL;>A58Bjw4C?qK^q~!?mc0a{3X<=}+PV%*`xR%fr~% zn4Xo6y<5Wi^mAEqR(}Zdnazrdxm^5bVwB4^dZ~Vc_sVfa#E;)6$5FjU?w8}nF}+qx zj_aN4!G{h{uU1^4SHt7OhpOT1cGk^@PghFPX6tq)faQnuI-1_-I34I393MfC4_4ou zekkPeVc=ujO8g(VVB*)e9!q&VT7e?It-^+0k z#4Y}V92Z60(r4v3gxMwrXNM!M)T?ejjPnDk9g^1j?Wlks(p1~3_@ELp? z_Th6VBP`I6-eD|&@^K7vZ)N|rAIVi2(yrUCY+;df!LiKU1}tyWu9N@TV9145noep@ zrE4kq<0yd1iB)c`T;x5RISBtW=3aFNfV83H*c4$A1`aIv*(Gb|n1NUSSVmUY{?l2U zrU6unXE68CtNfZB!A-1f`Ocr+qyF|w&9hh(`%%eALa6qauG1<_`@MV|B>trf406!a z7;>;i56d7`tS0#^u1aN1Er#Lyuf3D6e)qk<34N}udFQ?UU(edV_hW7UJ@-BM$m56I z-1e`^l;`vq@EGtIU<~l@JM5a8)_KX5rz)kXVlIs=?aSsW8RmvqMd6=sv+ELCoF=t6 zO={1+t9U^bFR0=JaXmeD1<3shMGN19eAB=x`Uv>_gU_0G2L~{)R=d;x`A92*c`db62Nb?X7Rc^I zFfRM=@^&3UT=<7_ToiHB56f|J#MOQz$F)a<>^x&MTG2L_)N?+y-RpnDt+nqdDvSU=iHYB5AcjU0}UV0R8v z+Vja@$lG<>PAtoDAuL_^=W<*SaiO2dab(AuOL82s?fa=5N9{lGUvk`9D`|Lo6?BDO z4Udm;j}NA#o%wFW&)e3+dJ5q})|G3x*3ePgVYin65tW$U8P z0Dlzy3kvwB(eDX&GeGZ~%|W#)O&r!VkHdfg6_4hvNlZV=oeRN#zoR`atic5|`l>$+ z{>4v#%wl?VZhp(3V?p{GyAo2O^cIiDaz6nf5lnCHF>nyB_#trBx${vGJ~YuiUh?Vp zrhL-+fB)-nGXne87+Tq^+z9l_r$mu}nPcR797A0FSUIkVnTc+e^;Z1TJ*?IN1^8A_WRlny#)KQ;6KQX&ZP|iCcPITGNN5kWTRdi?m zp^i8oR^R)d74YJ0J8a%w>{`cfB0s!6!-0ZBjEOlPR?lzx6qcVXI~Ca3`a$E)8QiKL z!1Dpr4D??VAIIWENaDDRN;!=VexK;z&q3hFS$6637e_>&?c(EJSMWjO&JgP3-0ygM zL>db+i{2)u&qkD@%wvFmt9)D+#M>p0mgB||SMK`)sf}rGi5&?4|;qsk~{N~ z&{x4?@$JMd`k|fds<5wd!|xdVolDlQya4Qf5c?m|`LDrC0at$xNc;8e-(0P4v<{xc zZcTlvykAujH~&mIu86q!7CEklxK==p^CK>HlpME+dJ;TBj$1<9_*3LKgqaF+y|Sy{ z|1s(t64K6Dt?hDxV|aW#G9rBR+wVW>YP(~uAN}(2^AY0%dOX^nMXjAAt&8-3PT}JM zo(@&}@x+MO53A4Lp2YII>-n3x{jjUQJsQW))bglzbz6I}svmwV|ETFlzx_+l)$=IQ zT6!;FN1~W2RdVUPetLkS>`3fNRnl3dRM?X!XEN<$d$jWf61%3Vg^5&2dq&O98Q4?G zl*;!*?*&wKpS4d<-1Z)QnyY{s`h{Lg1f6|M<(nPwg3l zJ+V#^UwI7a&ta@+z&XU~!2!Wg}DwH!z7HJOs* z1lb;Vu(@8@HJ^*2o{<5fqB+aL43Cf3dVDYy=&YCR91%YH9k*PH<@Qm7K&f(mXL|Yl zR$wYV;5iG@wQQpoQ7F87nV@dzKoe+P zmb99rb-P+YJ47vp=y2A&g!T)KWAc4Q7=ql;b98$B%)e9KKP+N8wO7b-^j=$( zHiii9-@pzlet|r_2;%%_%W-(Tozl9*1SiOk)nIeG)308gG6H&an9{&%m`Tg)V+Ny( zyX3egM%O3gxFw7ZUL(gTxGprFm*W;tF{aOy<3fl_2IaW5>Y2G-S)I29zwX9)%366l zSlr-vKTFWAM`K*o|ND67yOIy;Hy5#AYFsRD@AHVOy+)3UBW~^tIWCC`9(k@D7f0O8 ziE>;7am8oLaR{?bn76Zj`^_B!eNt`9^B@zlpN}g%K3Iiy)_dOQ@nPU&x5tNpkHU!a zVb}Aqy~q=GSZ$+7Pa4|+>+0yXr%^pW>%Wh1pFq#3#~uuA&gJ8Vb>RcXahsJ{RO9T0 z^6^gyqwC!BcK)2DAg0$kOP*d9aR<(ow57jxt*iuZ1b;|veLR4{?#p zl6T{Q1J6xex!{g(V9v@7LJM-)>czhW6_`1i3fscn=$_Mlm zH0!l*AK@>M4@#N31ZVZr-S(8{qu=MX%UEvefvMaoIUFtR=VQgye84zvlhQ&pp1w#v z{_!Dh_EmD+65)23;80r{uax5)tyC;?=6Ypy-B)qty6@ol@OgYNed~8V1c zG?j-7ifQ;VPcfa?%}yx9)BY3%(ogKm6(*+kG3Y++k1ko*L?sJ8Ej%%qDQ5{^nyThA zl}bk0S4vc}>|n+0R9-)Ok@?fvR3S|n%DwjDbF?&a*G{Gr?DWJ8JN8h!x$fVcn`R}$ zBMZw@`Mnu9gplz1d5ikVi|t7LI_=29slo(1b#RI$m?-bdWQuq~q8@)uX@7odZ-Qm3 z!VZhni=c;Bhz-0y1L%=-Kju%cgCQ#{8#+E$q0bq${f79-)HFeK{Eoy#ZW2<|{i_b7 z1b@1en%t$I)TsMcmmf2thfk(h7ILY4A~!L;I(&D&YVhkd0swZd<0^Lb^B#933YjX? z09|ERyy2)v=Eqoe-lcYcBOWZm!rKQzu8!%(Cs^tA`*h(UrkyYS0L(eSV zi781oubs%5p33KQq}mhZ$;4D?qOAzkYvBSKUDJ6`l#VAo2JB?Oe7s`UI&290C(>*x zMU5Pekmh{!J8yh}pkG!UQC_o5nDb$me!RpL`eC+SHiPe<@lOfNc=qB)g8}}0Gg@D% zbLWfgDdqI>V}^~Ww33rmj98?>OfMI!@SGuNF@5hccQe52Wdw^8LVHD}qX__wyZQ6? zCaZuiN7UXBYPE_R(9M|wJ6k>iNP;66Ey+Rga1 z92dm=MfS^a2s0JK)2r>S(5vC`vD4#&X=mqo?kz|Tv9SoS$QpNFWJ%%t?^`sIntUM;M8ZDOzXWh}hf=$qd7B(3xQ*+Aho@Jcc7LZQ9_R@ z=O%IA1R&x2kP3VcaWYj%?a54L3YB*MO7`BYJyU;7()u~u^;>&uc#VqY2nqe>W`)Kd z>U$v9qRi3rs`~A69PvGSvm6)18?|nc-Y