Skip to content
Snippets Groups Projects
Commit e8c92c71 authored by Murat Ambarkutuk's avatar Murat Ambarkutuk :robot:
Browse files

lets talk about sammich scripts

parent a2920820
Branches
No related merge requests found
Hide whitespace changes
Inline Side-by-side
code/sammich_katsidimas.m 0 → 100644
+ 19
− 0
View file @ e8c92c71
clc; close all; clear all; prepare_generic;
experiment = load('data/katsidimas_etal_organized.mat');
d_true = pdist2(experiment.gt, experiment.sensor_locations);
% nx = 30; ny = 30;
% space = generate_space(0, 0, 30, 30, nx, ny);
% xy = [space.xx(:), space.yy(:)];
z = experiment.z;
size(z)
z = z - mean(z, 3);
e = energy(z, 3);
% figure(1);
% for i=[1,3]
% ax(i) = subplot(2,2,i);
% scatter(d_true(:, i), e(:, i));
% end
% linkaxes(ax, 'xy');
scatter(d_true(:, [1,3]), e(:, [1,3]))
code/sammich_multi.m 0 → 100644
+ 106
− 0
View file @ e8c92c71
clc; clearvars; close all; prepare_generic;
load("../data/goodwin_step_data.mat");
experiment = experiment1;
sensors = [3, 4];
nx = 275; ny = 125;
d_true = pdist2(experiment.sensor_locations(:, sensors)', experiment.groundtruth');
space = generate_space(13, 35, 27, 33, nx, ny);
xy = [space.xx(:), space.yy(:)];
%%
for j=1:experiment.nsteps
e(:, j) = experiment.steps(j).e;
end
e = e(sensors, :);
ncases = experiment.nsteps*numel(sensors);
ntraining = round(0.6*ncases);
% training_idx = randsample(ncases, ntraining);
training_idx = round(linspace(1, ncases, ntraining));
ee = e(:);
dd = d_true(:);
etraining = ee(training_idx);
dtraining = dd(training_idx);
[model, gof] = fit(etraining, dtraining, 'power1', 'Robust', 'Bisquare', 'Lower', [0, -inf], 'Upper', [inf, 0]);
beta = [model.a; model.b];
ncases_testing = 1;
nk = 50;
kmax = 3;
k_space1 = linspace(0, kmax, nk);
k_space2 = linspace(0, kmax, nk);
[kk1, kk2] = meshgrid(k_space1, k_space2);
k12 = [kk1(:), kk2(:)];
%%
fxy = zeros(ny, nx, 2, nk*nk);
max_prob = zeros(nk*nk, 1);
step_idx = 25;
for k_idx=1:length(k12)
start_time = time();
k = k12(k_idx, :);
for sensor_idx = 1:2
et = experiment.steps(step_idx).e(sensors(sensor_idx));
nsamples = length(experiment.steps(step_idx).e(sensors(sensor_idx)));
xi = experiment.sensor_locations(:, sensors(sensor_idx))';
ei = experiment.steps(step_idx).e(sensors(sensor_idx));
prob = pdf_fx(xy, xi, ei*k(sensor_idx), nsamples, experiment.sigma_zeta*1500, beta);
fxy(:, :, sensor_idx, k_idx) = reshape(prob, [ny, nx]);
end
joint = squeeze(prod(fxy(:, :, :, k_idx), 3));
max_prob(k_idx) = max(joint(:));
duration = seconds(time()-start_time);
if mod(k_idx, 100) == 0
fprintf("Step ID: %d k ID: %d is completed in %.2f [secs]\n", step_idx, k_idx, duration);
end
end
max_prob_over_k_space = reshape(max_prob, nk, nk);
figure(1);
ax1 = subplot(121);
contourf(kk1, kk2, max_prob_over_k_space, 300, "edgecolor", "none")
shading interp; axis equal; hold on;
plot(k_space1, k_space2, 'r-.');
xlabel("$\kappa_1$"); ylabel("$\kappa_2$")
ax2 = subplot(122);
plot(k_space1, diag(max_prob_over_k_space), 'r-.', "Linewidth", 3);
xlabel("$\kappa \approx \kappa_1 \approx \kappa_2$")
ylabel("Maximum Probability $f(x)$");
grid on; grid minor;
linkaxes([ax1, ax2], 'x')
% parfor k_idx=1:nk
% start_time = time();
% margi = squeeze(sum(fxy(:, :, :, k_idx, step_idx), 3));
% joint = squeeze(prod(fxy(:, :, :, k_idx, step_idx), 3));
% % division = trapz(space.x, trapz(space.y, joint));
% % joint_pdf = joint / division;
% f = figure('visible', 'off');
% ax_1 = subplot(221);
% surf(space.xx, space.yy, margi); hold on;
% for sensor_idx = 1:experiment.nsensors
% scatter(experiment.sensor_locations(1, sensor_idx), experiment.sensor_locations(2, sensor_idx), 'rx');
% end
% % colorbar;
% scatter(experiment.groundtruth(1, step_idx), experiment1.groundtruth(2, step_idx), 'k.')
% axis equal; shading interp; view([0, 90]); drawnow;
% ax_2 = subplot(222);
% surf(space.xx, space.yy, joint); hold on;
% for sensor_idx = 1:experiment.nsensors
% scatter3(experiment.sensor_locations(1, sensor_idx), experiment.sensor_locations(2, sensor_idx), 100, 'rx');
% end
% % colorbar;
% scatter3(experiment.groundtruth(1, step_idx), experiment1.groundtruth(2, step_idx), 100, 'k.')
% axis equal; shading interp; view([0, 90]); drawnow;
% subplot(2,2, [3:4]);
% plot(k_space, [max_prob(step_idx, 1:k_idx), zeros(1, nk-k_idx)]); grid on; grid minor; xlabel('SNR $\kappa$'); ylabel('Max Joint-Probability');
% ylim([0 max(max_prob(step_idx, :))]);
% linkaxes([ax_1, ax_2], 'xy');
% fname = sprintf("../data/k_study/step_%d/result_%d.png", step_idx, k_idx);
% exportgraphics(gcf, fname, 'BackgroundColor', 'white', 'Resolution', 600);
% close(f);
% %
% duration = seconds(time()-start_time);
% fprintf("Step ID: %d k ID: %d is completed in %.2f [secs]\n", step_idx, k_idx, duration);
% end
% end
% save("../data/k_study/results.mat", "fxy", "space", "k_space", "beta");
code/sammich_single.m
+ 72
− 78
View file @ e8c92c71
clc; clearvars; close all; prepare_generic;
load("../data/goodwin_step_data.mat");
experiment = experiment1;
nx = 44; ny = 16;
nx = 275; ny = 125;
d_true = pdist2(experiment.sensor_locations', experiment.groundtruth');
space = generate_space(13, 35, 27, 35, nx, ny);
% dspace = linspace(0, 30, 100);
% theta_space = linspace(0, 2*pi, 100);
% p_theta = ones(size(theta_space))/(2*pi);
space = generate_space(13, 31, 27, 33, nx, ny);
xy = [space.xx(:), space.yy(:)];
%%
for j=1:experiment.nsteps
......@@ -20,102 +17,99 @@ ee = e(:);
dd = d_true(:);
etraining = ee(training_idx);
dtraining = dd(training_idx);
[model, gof] = fit(etraining, dtraining, 'power1', 'Robust', 'Bisquare', 'Lower', [0, -inf], 'Upper', [inf, 0]);
figure(666); hold on;
plot(linspace(0, 0.8, 10), model.a*linspace(0, 0.8, 10).^model.b);
plot(etraining, dtraining, 'rx');
grid on; grid minor;
% exportgraphics(gcf, '../figures/propagation_curve.png', 'BackgroundColor', 'none', 'Resolution', 60);
close(666);
beta = [model.a; model.b];
ncases_testing = 1;
%%
best_error = inf;
best_model = NaN;
iter_max = 1000;
n_batch = 10;
error_threshold = 2;
point_threshold = 50;
qq = linspace(0, .6, 100);
for iter=1:iter_max
confirmed_inliers = [];
maybe_inliers = sort(randsample(training_idx, n_batch));
[maybe_model, gof] = fit(ee(maybe_inliers), dd(maybe_inliers), 'power1', 'Robust', 'Bisquare', 'Lower', [0, -inf], 'Upper', [inf, 0]);
evals = feval(maybe_model, etraining);
error = evals - dtraining;
idx = find(error <= error_threshold);
confirmed_inliers = [confirmed_inliers; idx];
if length(confirmed_inliers) > point_threshold;
[better_model, better_gof] = fit(ee(confirmed_inliers), dd(confirmed_inliers), 'power1', 'Robust', 'Bisquare', 'Lower', [0, -inf], 'Upper', [inf, 0]);
if better_gof.sse < best_error
best_model = better_model;
best_error = better_gof.sse
end
end
mean_error = best_error/ncases;
if mean_error < 1
fprintf("Best Error: %3.3f ", mean_error);
break;
end
end
beta = [best_model.a; best_model.b];
nk = 100;
k_space = linspace(0, 1, nk);
sensors = [4, 5, 6, 7, 8, 9, 10];
nsensors = length(sensors);
% step_indices = sort(randsample(experiment.nsteps, ncases_testing));
k_space = linspace(0, 10, nk);
%%
fxy = zeros(ny, nx, experiment.nsensors, nk, experiment.nsteps);
max_prob = zeros(experiment.nsteps, nk);
for step_idx = 1:experiment.nsteps
fxy = zeros(ny, nx, nsensors, nk);
max_prob = zeros(1, nk);
folder_name = sprintf("../data/k_study/step_%d", step_idx);
mkdir(folder_name)
for qq=1:length(k_space)
k = k_space(qq);
mkdir(folder_name);
for k_idx=1:nk
start_time = time();
for pp = 1:nsensors
sensor_idx = sensors(pp);
k = k_space(k_idx);
parfor sensor_idx = 1:experiment.nsensors
et = experiment.steps(step_idx).e(sensor_idx);
nsamples = length(experiment.steps(step_idx).e(sensor_idx));
for xy_idx = 1:length(xy)
[I,J] = ind2sub([ny, nx], xy_idx);
xy_interest = xy(xy_idx, :);
xi = experiment.sensor_locations(:, sensor_idx)';
ei = experiment.steps(step_idx).e(sensor_idx);
di = vecnorm(xy_interest - xi);
fxy(I, J, pp, qq) = pdf_fx(xy_interest, xi, ei*k, nsamples, experiment.sigma_zeta*1500, beta);
end
xi = experiment.sensor_locations(:, sensor_idx)';
ei = experiment.steps(step_idx).e(sensor_idx);
prob = pdf_fx(xy, xi, ei*k, nsamples, experiment.sigma_zeta*1500, beta);
qwert = reshape(prob, [ny, nx]);
division = trapz(space.x, trapz(space.y, qwert));
qwert_pdf = qwert / division;
fxy(:, :, sensor_idx, k_idx, step_idx) = qwert_pdf;
end
margi = squeeze(sum(fxy(:, :, :, qq), 3));
joint = squeeze(prod(fxy(:, :, :, qq), 3));
division = trapz(space.x, trapz(space.y, joint));
joint_pdf = joint / division;
max_prob(qq) = max(joint(:));
joint = squeeze(prod(fxy(:, :, :, k_idx, step_idx), 3));
max_prob(step_idx, k_idx) = max(joint(:));
duration = seconds(time()-start_time);
fprintf("Step ID: %d k ID: %d is completed in %.2f [secs]\n", step_idx, k_idx, duration);
end
parfor k_idx=nk:nk
start_time = time();
margi = squeeze(sum(fxy(:, :, :, k_idx, step_idx), 3));
joint = squeeze(prod(fxy(:, :, :, k_idx, step_idx), 3));
% division = trapz(space.x, trapz(space.y, joint));
% joint_pdf = joint / division;
f = figure('visible', 'off');
ax(1) = subplot(221);
ax_1 = subplot(221);
surf(space.xx, space.yy, margi); hold on;
for pp = 1:nsensors
sensor_idx = sensors(pp);
for sensor_idx = 1:experiment.nsensors
scatter(experiment.sensor_locations(1, sensor_idx), experiment.sensor_locations(2, sensor_idx), 'rx');
end
% colorbar;
scatter(experiment.groundtruth(1, step_idx), experiment1.groundtruth(2, step_idx), 'k.')
axis equal; shading interp; view([0, 90]); drawnow;
ax(2) = subplot(222);
ax_2 = subplot(222);
surf(space.xx, space.yy, joint); hold on;
for pp = 1:nsensors
sensor_idx = sensors(pp);
for sensor_idx = 1:experiment.nsensors
scatter3(experiment.sensor_locations(1, sensor_idx), experiment.sensor_locations(2, sensor_idx), 100, 'rx');
end
% colorbar;
scatter3(experiment.groundtruth(1, step_idx), experiment1.groundtruth(2, step_idx), 100, 'k.')
axis equal; shading interp; view([0, 90]); drawnow;
subplot(2,2, [3:4]);
plot(k_space, max_prob); grid on; grid minor; xlabel('SNR $\kappa$'); ylabel('Max Joint-Probability');
linkaxes(ax, 'xy');
fname = sprintf("../data/k_study/step_%d/result_%d.png", step_idx, qq);
exportgraphics(gcf, fname, 'BackgroundColor', 'none', 'Resolution', 600);
plot(k_space, [max_prob(step_idx, 1:k_idx), zeros(1, nk-k_idx)]); grid on; grid minor; xlabel('SNR $\kappa$'); ylabel('Max Joint-Probability');
% ylim([0 max(max_prob(step_idx, :))]);
linkaxes([ax_1, ax_2], 'xy');
fname = sprintf("../data/k_study/a_%d_result_%d.png", step_idx, k_idx);
exportgraphics(gcf, fname, 'BackgroundColor', 'white', 'Resolution', 600);
close(f);
%
duration = seconds(time()-start_time);
fprintf("k ID: %d is completed in %.2f [secs]\n", qq, duration);
end
f = figure('visible', 'off');
plot(k_space, max_prob); grid on; grid minor; xlabel('SNR $\kappa$'); ylabel('Max Joint-Probability'); "experiment"
fname = sprintf('../data/k_study/step_%d/prob.png', step_idx);
exportgraphics(f, fname, 'BackgroundColor', 'none', 'Resolution', 600);
close(f);
fname = sprintf("../data/k_study/step_%d/results.mat", step_idx);
save(fname, "fxy", "joint", "joint_pdf", "margi", "space");
fname_out = sprintf('../data/k_study/step_%d/anim.gif', step_idx);
fname = sprintf("../data/k_study/step_%d/result_*.png", step_idx);
filepath = dir(fname);
filepath = natsortfiles(filepath);
for k = 1:numel(filepath)
fname = fullfile(filepath(k).folder, filepath(k).name);
img = imread(fname);
[A,map] = rgb2ind(img,256);
if k ==1
imwrite(A, map, fname_out, 'gif', 'LoopCount', Inf, 'DelayTime', 1);
else
imwrite(A, map, fname_out, 'gif', 'WriteMode', 'append', 'DelayTime', 1);
end
fprintf("Step ID: %d k ID: %d is completed in %.2f [secs]\n", step_idx, k_idx, duration);
end
end
% save("../data/k_study/results.mat", "fxy", "space", "k_space", "beta");
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