Commit cb1802f5 authored by Leon Merten Lohse's avatar Leon Merten Lohse

tut nun alles

parent fd1e8306
......@@ -194,7 +194,7 @@ ArrayXX compute_image(const std::string & input_filename,
throw std::runtime_error( "inconsistent number of energy bins");
}
if (! (psi_x.cols() == hitmap_resolution)) {
if (! (psi_x.cols() == hitmap_resolution * hitmap_resolution)) {
throw std::runtime_error( "inconsistent resolution of psi_x");
}
......
import etamap
import numpy as np
from scipy import ndimage
import _moench
def compute_hitmaps(filename, hitmap_resolution, min_energy, max_energy, n_bins, gainmap):
gainmap_raw = gainmap.flatten().astype(np.float32)
hitmaps_raw = _moench.compute_hitmaps(filename,
hitmap_resolution,
min_energy,
max_energy,
n_bins,
gainmap_raw)
hitmaps = hitmaps_raw.reshape([n_bins, hitmap_resolution, hitmap_resolution])
return hitmaps
def compute_images(filename, hitmap_resolution, min_energy, max_energy, n_bins, gainmap, n_int, psi_x, psi_y):
gainmap_raw = gainmap.flatten().astype(np.float32)
psi_x_raw = psi_x.reshape([n_bins, -1]).astype(np.float32)
psi_y_raw = psi_y.reshape([n_bins, -1]).astype(np.float32)
images_raw = _moench.compute_image( filename,
hitmap_resolution,
min_energy,
max_energy,
n_bins,
gainmap_raw,
n_int,
psi_x_raw,
psi_y_raw )
images = images_raw.reshape([n_bins, _moench.dy * n_int, _moench.dx * n_int])
return images
def compute_spectrum(filename, min_energy, max_energy, n_bins, gainmap):
spectrum = compute_hitmaps(filename, 1, min_energy, max_energy, n_bins, gainmap)
energy_edges = np.linspace(min_energy, max_energy, n_bins+1)
energy_centers = (energy_edges[1:] + energy_edges[:-1])/2
return energy_centers, spectrum.reshape([n_bins,])
def compute_psi(eta):
# TODO: ueberarbeiten!!
# normalize
mu_0 = eta * eta.size / np.sum(eta.ravel())
u_0 = etamap.calc_a_and_b(mu_0)
iterations = 1000
step = 1e-2
print(f"Running gradient descent with max_iterations={iterations} and step_size={step}")
u, M, i_last = etamap.calc_u_gradient_descent(u_0, mu_0, max_iterations=iterations, gamma=step)
print(f"\nGradient descent stopped after {i_last} iterations")
# compare initial and final curl
curl0_f = etamap.curl(u_0)
curl1_f = etamap.curl(u)
curl0 = np.sqrt((curl0_f**2).sum())
curl1 = np.sqrt((curl1_f**2).sum())
print("Initial solution:")
print(f" Monge Kantorovich: {M[0]:.2e}, total L2-curl: {curl0:.2e}")
print("Final solution:")
print(f" Monge Kantorovich: {M[i_last]:.2e}, total L2-curl: {curl1:.2e}")
# remove values exceeding 1
for i in range(2):
u[i][u[i] > mu_0.shape[i]] = mu_0.shape[i]
u[i][u[i] < 0] = 0
u_norm_0 = (u[0] / mu_0.shape[0]).astype(np.float32)
u_norm_1 = (u[1] / mu_0.shape[1]).astype(np.float32)
return u_norm_0.T.copy(), u_norm_1.T.copy()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment