Skip to content

update normalisation for average kersize #14

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nm19000 merged 4 commits into from
Jun 8, 2026
Merged

update normalisation for average kersize #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
82756a7
update normalisation for average kersize
nm19000 Jun 2, 2026
986012b
update comments
nm19000 Jun 4, 2026
50e9534
correct output of appertenance computation
nm19000 Jun 5, 2026
670a9e4
ran tests for pytorch version, used num_workers=0 for mac, test passed
nm19000 Jun 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 8 additions & 8 deletions examples/example_linalg.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

230 changes: 119 additions & 111 deletions examples/example_localizationmicroscopy.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

16 changes: 9 additions & 7 deletions src/accuracy_bounds/inverseproblems/kersize_compute.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,8 @@ def diams_feasibleset(feasible_set_y, p_1 ,p):
- p: Order of the average kernel size. Set to p=2 for the MSE lower bound computation and p=1 for MAE lower bound computation.

Returns:
- diameter_mean_y, num_feas, max_diam_Fy: diameter_mean_y of dim(0)= shape(input_data), the estimated mean diameter of the feasible set to the power p,
- diameter_mean_y, num_feas, max_diam_Fy: diameter_mean_y of dim(0)= shape(input_data), the estimated "mean"
(normalized by num_feas-1, i.e. lacking a factor num_feas) diameter of the feasible set to the power p,
consisting of all possible target data points, for one input point.
num_feas is the number of samples in the feasible set and will be used for statistics later on.
max_diam_Fy the maximum diameter of the feasible set,
Expand Down Expand Up @@ -44,10 +45,10 @@ def diams_feasibleset(feasible_set_y, p_1 ,p):

# get mean over diams, with factor 2 due to symmetry of the norm of the compute vectors in null space of F (norm(x-z)=norm(z-x))
# and divided by num_feas^2 ad we have that many terms
if num_feas > 0:
# compute 2 times sum over diams to the power p divided by num_feas^2
diameter_mean_y = 2*np.divide(np.sum(np.power(diam_y,p)), np.power(num_feas,2))
elif num_feas==0:
if num_feas > 1:

@JakobCode JakobCode Jun 4, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happens for case num_feas == 1?

# compute 2 times sum over diams to the power p divided by num_feas-1
diameter_mean_y = 2*np.divide(np.sum(np.power(diam_y,p)), num_feas-1)

@JakobCode JakobCode Jun 4, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please double check that this is correct and align with comment above

@JakobCode JakobCode Jun 4, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not (num_feas-1)^p ?

else:
diameter_mean_y = 0

return diameter_mean_y, num_feas, max_diam_Fy
Expand Down Expand Up @@ -173,16 +174,18 @@ def average_kernelsize(feasible_sets_list, p_1, p):

average_kersize = 0
num_samples = len(feasible_sets_list)
normalization_m=0

for feasible_set_y in feasible_sets_list:
# compute diameter of feasible set for one input data point (num_feas will be used for statistics later on)
diameter_mean_y, num_feas, max_diam_Fy = diams_feasibleset(feasible_set_y, p_1 ,p)
#add diameters means for obtaining average kersize to the power p
average_kersize = average_kersize + diameter_mean_y
normalization_m = normalization_m+num_feas

# get mean over input data
if average_kersize>0 and num_samples > 0:
average_kersize = np.divide(average_kersize, num_samples)
average_kersize = np.divide(average_kersize, normalization_m)
else:
average_kersize = 0
# take power 1/p to obtain average kersize
Expand Down Expand Up @@ -225,4 +228,3 @@ def average_kernelsize_sym(A, input_data, target_data, p_1, p_2, p, epsilon):
average_kersize_sym = np.power(average_kersize_sym, 1/p)

return average_kersize_sym

19 changes: 16 additions & 3 deletions src/accuracy_bounds/inverseproblems/kersize_compute_dataloader.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -262,15 +262,28 @@ def compute_mean_distance(y_idx, feasible_appertinance_matrix, dXX):
subdistXX = subdistXX.toarray() # Convert to dense matrix for max computation

size_feas = len(valid_idx)
if size_feas <= 1:
return 0.0, 0

subdistXX = dXX[valid_idx, :][:, valid_idx].toarray()

# compute average over sums of differences over twice the feasible sets
return np.divide(np.nansum(np.power(subdistXX,p_X)), np.power(size_feas,2))
value = np.nansum(np.power(subdistXX, p_X)) / (size_feas - 1)
return value, size_feas


n,p = feasible_appartenance.shape
results = Parallel(n_jobs=-1)(delayed(compute_mean_distance)(y_idx, feasible_appartenance, distsXX) for y_idx in tqdm(range(p)))
output = Parallel(n_jobs=-1)(delayed(compute_mean_distance)(y_idx, feasible_appartenance, distsXX) for y_idx in tqdm(range(p)))
results, size_feas_list = zip(*output)

normalization_m= np.sum(np.asarray(list(size_feas_list)))
normalization_m = np.sum(size_feas_list)

if normalization_m == 0:
return 0.0

# get average over K input data samples (y)
return np.power(np.nanmean(np.asarray(list(results))),np.divide(1,p_X))
return np.power(np.divide(np.sum(np.asarray(list(results))),normalization_m),np.divide(1,p_X))



Expand Down
21 changes: 12 additions & 9 deletions test/test_toy_example_torch.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
from accuracy_bounds.inverseproblems.utils import apply_forwardmodel
from accuracy_bounds.inverseproblems.kersize_compute_dataloader import worstcase_kernelsize_appartenance, average_kernelsize_appartenance, target_distances_cuda_V2
from accuracy_bounds.inverseproblems.feasible_sets_dataloader import feasible_appartenance_additive_noise_dataloader_cuda, feasible_appartenance_additive_noise_cuda
from accuracy_bounds.inverseproblems.utils import torch_sparse_to_scipy_csr, torch_csr_to_scipy
from accuracy_bounds.inverseproblems.utils_torch import torch_sparse_to_scipy_csr, torch_csr_to_scipy


@pytest.mark.skipif(not has_torch, reason="requires torch")
Expand All @@ -40,12 +40,13 @@ def test_worstcase_cuda():

max_k = 3000
batch_size = 100
num_workers = 0

input_data_k = input_data[0:max_k,:]
target_data_k = target_data[0:max_k,:]
input_loader1 = DataLoader(input_data_k, batch_size=batch_size, num_workers=batch_size, drop_last=False)
input_loader2 = DataLoader(input_data_k, batch_size=batch_size, num_workers=batch_size, drop_last=False)
target_loader1 = DataLoader(target_data_k, batch_size=batch_size, num_workers=batch_size, drop_last=False)
input_loader1 = DataLoader(input_data_k, batch_size=batch_size, num_workers=num_workers, drop_last=False)
input_loader2 = DataLoader(input_data_k, batch_size=batch_size, num_workers=num_workers, drop_last=False)
target_loader1 = DataLoader(target_data_k, batch_size=batch_size, num_workers=num_workers, drop_last=False)

feasible_appartenance = feasible_appartenance_additive_noise_dataloader_cuda(input_loader1, input_loader2, p_Y=p_2, epsilon= epsilon)
feasible_appartenance = feasible_appartenance.to(dtype=torch.float32).to_sparse_coo()
Expand Down Expand Up @@ -83,12 +84,13 @@ def test_average_kersize_cuda():

max_k = 3000
batch_size = 100
num_workers = 0

input_data_k = input_data[0:max_k,:]
target_data_k = target_data[0:max_k,:]
input_loader1 = DataLoader(input_data_k, batch_size=batch_size, num_workers=batch_size, drop_last=False)
input_loader2 = DataLoader(input_data_k, batch_size=batch_size, num_workers=batch_size, drop_last=False)
target_loader1 = DataLoader(target_data_k, batch_size=batch_size, num_workers=batch_size, drop_last=False)
input_loader1 = DataLoader(input_data_k, batch_size=batch_size, num_workers=num_workers, drop_last=False)
input_loader2 = DataLoader(input_data_k, batch_size=batch_size, num_workers=num_workers, drop_last=False)
target_loader1 = DataLoader(target_data_k, batch_size=batch_size, num_workers=num_workers, drop_last=False)

feasible_appartenance = feasible_appartenance_additive_noise_dataloader_cuda(input_loader1, input_loader2, p_Y=p_2, epsilon= epsilon)
feasible_appartenance = feasible_appartenance.to(dtype=torch.float32).to_sparse_coo()
Expand Down Expand Up @@ -122,9 +124,10 @@ def test_feas_appartance_w_and_wo_dataloader():
## Test if feasible appartenance matrix computation versions produce the same results

batch_size = 100
num_workers = 0

input_loader1 = DataLoader(input_data, batch_size=batch_size, num_workers=batch_size, drop_last=False)
input_loader2 = DataLoader(input_data, batch_size=batch_size, num_workers=batch_size, drop_last=False)
input_loader1 = DataLoader(input_data, batch_size=batch_size, num_workers=num_workers, drop_last=False)
input_loader2 = DataLoader(input_data, batch_size=batch_size, num_workers=num_workers, drop_last=False)

feas_app_1 = feasible_appartenance_additive_noise_dataloader_cuda(input_data=input_loader1, forwarded_target= input_loader2, p_Y=2, epsilon=epsilon)
feas_app_2 = feasible_appartenance_additive_noise_cuda(input_loader1, input_loader2, p_Y=2, epsilon= epsilon, batchsize=50)
Expand Down
Loading