aare/python/tests/test_interpolation.ipynb at dev/pipe

mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-04-20 21:14:36 +02:00

Files

T

mazzol_a 267ca87ab0 Dev/rosenblatttransform (#241 )

- added rosenblatttransform 
- added 3x3 eta methods 
- interpolation can be used with various eta functions
- added documentation for interpolation, eta calculation 
- exposed full eta struct in python 
- disable ClusterFinder for 2x2 clusters 
- factory function for ClusterVector

---------

Co-authored-by: Dhanya Thattil <dhanya.thattil@psi.ch>
Co-authored-by: Erik Fröjdh <erik.frojdh@psi.ch>

2025-11-21 14:48:46 +01:00

919 KiB

Raw Permalink Blame History

In [1]:

import numpy as np 
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt
from aare import ClusterFinder, calculate_eta2, Interpolator
import boost_histogram as bh
from aare._aare import Cluster2x2d, ClusterVector_Cluster2x2d
import pickle
import os

from test_Interpolation import create_photon_hit_with_gaussian_distribution, create_2x2cluster_from_frame, create_3x3cluster_from_frame, calculate_eta_distribution, photon_hit_in_euclidean_space

In [2]:

### plotting functions

def plot_eta_distribution(eta_distribution):
    plt.imshow(
    eta_distribution[:, :, 0].view().transpose(),   
    origin='lower',                  
    extent=[eta_distribution.axes[0].edges[0], eta_distribution.axes[0].edges[-1], eta_distribution.axes[1].edges[0], eta_distribution.axes[1].edges[-1]],
    aspect='auto')
    plt.colorbar()
    plt.xlabel('eta_x')
    plt.ylabel('eta_y')

def plot_uniform_eta_distribution(uniform_etax, uniform_etay, eta_extent): 
    fig, axes = plt.subplots(1, 2, figsize=(10, 5))  # (rows, cols)

    im1 = axes[0].imshow(
    uniform_etax[:, :, 0].view().transpose(),   
    origin='lower',                  
    extent=[eta_extent[0], eta_extent[1], eta_extent[0], eta_extent[1]],
    aspect='auto'
    )
    axes[0].set_title("uniform distribution of etax")
    fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04)
    axes[0].set_xlabel("eta_x")
    axes[0].set_ylabel("eta_y")

    im2 = axes[1].imshow(
    uniform_etay[:, :, 0].view().transpose(),   
    origin='lower',                  
    extent=[eta_extent[0], eta_extent[1], eta_extent[0], eta_extent[1]],
    aspect='auto'
    )
    axes[1].set_title("uniform distribution of etay")
    fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04)
    axes[1].set_xlabel("eta_x")
    axes[1].set_ylabel("eta_y")

    plt.tight_layout()
    plt.show()

def plot_photon_hit_distribution(eta_distribution, marg_CDF_etax, cond_CDF_etay): 
    normalized_eta_distribution = eta_distribution.values()[:,:, 0] / eta_distribution.values()[: , :, 0].sum() #normnalize 

    n_samples = 1000
    rng = np.random.default_rng(42)

    flat_indices = rng.choice(
        normalized_eta_distribution.size,
        size=n_samples,
        p=normalized_eta_distribution.ravel())

    x_indices, y_indices = np.unravel_index(flat_indices, normalized_eta_distribution.shape)

    photon_position_x = marg_CDF_etax[x_indices, y_indices]
    photon_position_y = cond_CDF_etay[x_indices, y_indices]

    fig, axes = plt.subplots(1, 2, figsize=(10, 5))  # (rows, cols)

    axes[0].scatter(eta_distribution.axes[0].edges[:-1][x_indices], eta_distribution.axes[1].edges[:-1][y_indices], s=1, alpha=0.5)
    axes[0].set_title("Eta distribution (eta_x, eta_y)")
    #plt.gca().set_aspect('equal')
    axes[0].set_xlabel("eta_x")
    axes[0].set_ylabel("eta_y")
    #axes[0].set_xlim(0,1)
    #axes[0].set_ylim(0,1)

    axes[1].scatter(photon_position_x, photon_position_y, s=1, alpha=0.5)
    axes[1].set_title("uniform Photon positions")
    #plt.gca().set_aspect('equal')
    axes[1].set_xlabel("x")
    axes[1].set_ylabel("y")
    #axes[1].set_xlim(0,1)
    #axes[1].set_ylim(0,1)

In [3]:

pixel_width = 1e-4

values = np.arange(0.5*pixel_width, 0.1, pixel_width)

num_pixels = values.size

print(f"num pixels: ({values.size} x {values.size})")

X, Y = np.meshgrid(values, values)

data_points = np.stack([X.ravel(), Y.ravel()], axis=1)

num pixels: (1000 x 1000)

In [4]:

# mean is in bottom right quadrant 
variance = 10*pixel_width

covariance_matrix = np.array([[variance, 0], [0, variance]])

mean = 650*pixel_width

mean = np.array([mean, mean])

base_frame = create_photon_hit_with_gaussian_distribution(mean, variance, data_points)

pixels_per_superpixel = int(num_pixels*0.5)

plt.imshow(base_frame) 

plt.colorbar()

plt.axvline(x=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axhline(y=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)

plt.xlabel('x_0')
plt.ylabel('y_0')

Out[4]:

Text(0, 0.5, 'y_0')

No description has been provided for this image

In [5]:

### create eta distribution
num_frames = 5000
random_number_generator = np.random.default_rng(42)
eta_distribution = calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, bh.axis.Regular(100, -0.1, 0.6), bh.axis.Regular(100, -0.1, 0.6))

test_data_path = os.getenv("AARE_TEST_DATA") + "/eta_distributions"
filename = test_data_path + "/eta_distribution_2x2cluster_gaussian.pkl"
with open(filename, "wb") as f:
    pickle.dump(eta_distribution, f)

In [7]:

filename = os.getenv("AARE_TEST_DATA") + "/eta_distributions/eta_distribution_2x2cluster_gaussian.pkl"
with open(filename, "rb") as f:
        eta_distribution = pickle.load(f)

In [8]:

plot_eta_distribution(eta_distribution)

In [9]:

interpolator = Interpolator(eta_distribution, eta_distribution.axes[0].edges, eta_distribution.axes[1].edges, eta_distribution.axes[2].edges[:-1])
#interpolator = Interpolator(eta_distribution.axes[0].edges, eta_distribution.axes[1].edges, eta_distribution.axes[2].edges[:-1])

#interpolator.rosenblatttransform(eta_distribution.values())

In [18]:

marg_CDF_etax = interpolator.get_ietax()
cond_CDF_etay = interpolator.get_ietay()

plot_photon_hit_distribution(eta_distribution, marg_CDF_etax, cond_CDF_etay)

In [10]:

marg_CDF_etax = interpolator.get_ietax()
cond_CDF_etay = interpolator.get_ietay()

#plot_uniform_eta_distribution(uniform_etax, uniform_etay, [eta_distribution.axes[0].edges[0], eta_distribution.axes[0].edges[-1]])

In [11]:

cluster = create_2x2cluster_from_frame(base_frame, pixels_per_superpixel)

print("cluster.x:", cluster.x)
print("cluster.y:", cluster.y)
clustervec = ClusterVector_Cluster2x2d()
clustervec.push_back(cluster)

eta = calculate_eta2(cluster)
print("eta: ", eta)
bin_size = (eta_distribution.axes[0].edges[-1] - eta_distribution.axes[0].edges[0])/eta_distribution.axes[0].edges.shape[0]

bin_index_x = int((eta[0] - eta_distribution.axes[0].edges[0])/bin_size)
bin_index_y = int((eta[1] - eta_distribution.axes[1].edges[0])/bin_size)

print("distance x:", marg_CDF_etax[bin_index_x, bin_index_y, 0])
print("distance y:", cond_CDF_etay[bin_index_x, bin_index_y, 0])

photon_hit = interpolator.interpolate(clustervec)

print(photon_hit)

cluster.x: 1
cluster.y: 1
eta:  (0.3519491547025806, 0.3519491547025806, 71553511.0)
distance x: 0.7083333333333334
distance y: 0.7564102564102564
[(0.78166653, 0.81284991, 71553511.)]

In [12]:

# scale to cluster pixel width 
#cluster_center = num_pixels*0.75*pixel_width
cluster_center = (pixels_per_superpixel*pixel_width*(cluster.x + 0.5), pixels_per_superpixel*pixel_width*(cluster.y + 0.5))
#scaled_photon_hit = photon_hit_in_euclidean_space(cluster_center, pixels_per_superpixel, photon_hit)

print(photon_hit[0])

scaled_photon_hit = ((photon_hit[0][0])*pixels_per_superpixel*pixel_width, (photon_hit[0][1])*pixels_per_superpixel*pixel_width)

print(f"previous center: ({cluster_center[0]}, {cluster_center[1]})")
print(f"interpolated center: ({photon_hit[0][0]},{photon_hit[0][1]})")
print(f"scaled interpolated center: ({scaled_photon_hit[0]},{scaled_photon_hit[1]})")
print(f"actual center: ({mean},{mean})")

(0.7816665310999981, 0.8128499138725078, 71553511.0)
previous center: (0.07500000000000001, 0.07500000000000001)
interpolated center: (0.7816665310999981,0.8128499138725078)
scaled interpolated center: (0.03908332655499991,0.040642495693625394)
actual center: ([0.065 0.065],[0.065 0.065])

In [13]:

plt.imshow(base_frame) 

plt.colorbar()

plt.axvline(x=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axhline(y=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)

plt.scatter(mean/pixel_width, mean/pixel_width, color='red', s=40, marker='x', label='Actual photon hit')
plt.scatter(scaled_photon_hit[0]/pixel_width, scaled_photon_hit[1]/pixel_width, color = 'blue', s=40, marker='x', label='interpolated photon hit')

plt.xlabel('x_0')
plt.ylabel('y_0')

Out[13]:

Text(0, 0.5, 'y_0')

3x3 Cluster¶

In [14]:

from aare._aare import Cluster3x3d, ClusterVector_Cluster3x3d

In [15]:

# mean is in center quadrant 
variance = 10*pixel_width

covariance_matrix = np.array([[variance, 0], [0, variance]])

mean_x = (1 + 0.8)*(num_pixels/3)*pixel_width
mean_y = (1 + 0.2)*(num_pixels/3)*pixel_width

mean = np.array([mean_x, mean_y])

base_frame = create_photon_hit_with_gaussian_distribution(mean, variance, data_points)

pixels_per_superpixel = int(num_pixels/3)

plt.imshow(base_frame) 

plt.colorbar()

plt.axvline(x=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axhline(y=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axvline(x=2*pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axhline(y=2*pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)

plt.xlabel('x_0')
plt.ylabel('y_0')

Out[15]:

Text(0, 0.5, 'y_0')

In [ ]:

### calculate eta distribution

num_frames = 5000

random_number_generator = np.random.default_rng(42)
eta_distribution = calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, bh.axis.Regular(100, -0.1, 1.0), bh.axis.Regular(100, -0.1, 1.0), False)

filename = os.getenv("AARE_TEST_DATA") + "/eta_distributions/eta_distribution_3x3cluster_gaussian.pkl"
with open(filename, "wb") as f:
    pickle.dump(eta_distribution, f)

In [ ]:

filename = os.getenv("AARE_TEST_DATA") + "/eta_distributions/eta_distribution_3x3cluster_gaussian.pkl"
with open(filename, "rb") as f:
        eta_distribution = pickle.load(f)

In [7]:

plot_eta_distribution(eta_distribution)

In [ ]:

In [14]:

#interpolator = Interpolator(eta_distribution, eta_distribution.axes[0].edges, eta_distribution.axes[1].edges, eta_distribution.axes[2].edges[:-1])

interpolator = Interpolator(eta_distribution.axes[0].edges, eta_distribution.axes[1].edges, eta_distribution.axes[2].edges[:-1])
interpolator.rosenblatttransform(eta_distribution.values())
marg_CDF_etax = interpolator.get_ietax()
cond_CDF_etay = interpolator.get_ietay()

In [15]:

plot_photon_hit_distribution(eta_distribution, marg_CDF_etax, cond_CDF_etay)

In [ ]:

plot_uniform_eta_distribution(uniform_etax, uniform_etay, [eta_distribution.axes[0].edges[0], eta_distribution.axes[0].edges[-1]])

In [16]:

cluster = create_3x3cluster_from_frame(base_frame, pixels_per_superpixel)

eta = calculate_eta2(cluster)
print(eta)
print(eta_distribution.axes[0].edges)
clustervec = ClusterVector_Cluster3x3d()
#print(marg_CDF_etax)
print(marg_CDF_etax[93, 0,0])
clustervec.push_back(cluster)



photon_hit = interpolator.interpolate(clustervec)

(0.5499154381017892, 0.44933336606613455, 48591988.0)
[-0.1   -0.093 -0.086 -0.079 -0.072 -0.065 -0.058 -0.051 -0.044 -0.037
 -0.03  -0.023 -0.016 -0.009 -0.002  0.005  0.012  0.019  0.026  0.033
  0.04   0.047  0.054  0.061  0.068  0.075  0.082  0.089  0.096  0.103
  0.11   0.117  0.124  0.131  0.138  0.145  0.152  0.159  0.166  0.173
  0.18   0.187  0.194  0.201  0.208  0.215  0.222  0.229  0.236  0.243
  0.25   0.257  0.264  0.271  0.278  0.285  0.292  0.299  0.306  0.313
  0.32   0.327  0.334  0.341  0.348  0.355  0.362  0.369  0.376  0.383
  0.39   0.397  0.404  0.411  0.418  0.425  0.432  0.439  0.446  0.453
  0.46   0.467  0.474  0.481  0.488  0.495  0.502  0.509  0.516  0.523
  0.53   0.537  0.544  0.551  0.558  0.565  0.572  0.579  0.586  0.593
  0.6  ]
0.9999999999999998
dX: 1
dY: 0
u: 1
v: 0

In [17]:

## scale to cluster pixel width 
#cluster_center = (1+0.5)*pixels_per_superpixel*pixel_width

cluster_center = (pixels_per_superpixel*pixel_width*(cluster.x + 0.5), pixels_per_superpixel*pixel_width*(cluster.y + 0.5))

print(photon_hit[0])

scaled_photon_hit = (photon_hit[0][0]*pixels_per_superpixel*pixel_width, photon_hit[0][1]*pixels_per_superpixel*pixel_width)

print(f"previous center: ({cluster_center[0]}, {cluster_center[1]})")
print(f"interpolated center: ({photon_hit[0][0]},{photon_hit[0][1]})")
print(f"scaled interpolated center: ({scaled_photon_hit[0]},{scaled_photon_hit[1]})")
print(f"actual center: ({mean},{mean})")

(1.5000000000000002, 1.5, 48591988.0)
previous center: (0.04995000000000001, 0.04995000000000001)
interpolated center: (1.5000000000000002,1.5)
scaled interpolated center: (0.04995000000000001,0.04995)
actual center: ([0.06 0.04],[0.06 0.04])

In [18]:

plt.imshow(base_frame) 

plt.colorbar()

plt.axvline(x=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axhline(y=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axvline(x=2*pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axhline(y=2*pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)

plt.scatter(mean_x/pixel_width, mean_y/pixel_width, color='red', s=40, marker='x', label='Actual photon hit')
plt.scatter(scaled_photon_hit[0]/pixel_width, scaled_photon_hit[1]/pixel_width, color = 'blue', s=40, marker='x', label='interpolated photon hit')

plt.xlabel('x_0')
plt.ylabel('y_0')

Out[18]:

Text(0, 0.5, 'y_0')

Interpolation with Eta3¶

In [6]:

from aare._aare import Cluster3x3d, ClusterVector_Cluster3x3d, calculate_eta3, calculate_cross_eta3

In [4]:

# mean is in center quadrant 
variance = 10*pixel_width

covariance_matrix = np.array([[variance, 0], [0, variance]])

mean = (1 + 0.3)*(num_pixels/3)*pixel_width

mean = np.array([mean, mean])

base_frame = create_photon_hit_with_gaussian_distribution(mean, variance, data_points)

pixels_per_superpixel = int(num_pixels/3)

plt.imshow(base_frame) 

plt.colorbar()

plt.axvline(x=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axhline(y=pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axvline(x=2*pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)
plt.axhline(y=2*pixels_per_superpixel, color='black', linestyle='--', linewidth=0.5)

plt.xlabel('x_0')
plt.ylabel('y_0')

Out[4]:

Text(0, 0.5, 'y_0')

In [7]:

num_frames = 2000

random_number_generator = np.random.default_rng(42)
eta_distribution = calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, calculate_eta3, False)

In [8]:

plot_eta_distribution(eta_distribution)

In [9]:

num_frames = 2000

random_number_generator = np.random.default_rng(42)
eta_distribution = calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, calculate_cross_eta3, False)

In [10]:

plot_eta_distribution(eta_distribution)

In [ ]:

919 KiB Raw Permalink Blame History

3x3 Cluster¶

Interpolation with Eta3¶

919 KiB

Raw Permalink Blame History