% deblur1_gcv_template.m
% GCV for (nonlinear!) edge-preserving image deblurring

if ~isvar('y') % read data
	load deblur1_data % loads image variable 'y' and PSF variable 'psf'
	[nx ny] = size(y);
	nd = nx * ny; % # of data points
	A = Gblur(true(nx,ny), 'psf', psf); % system matrix
	C = Cfd2(nx,ny); % 2D finite differences
end

im plc 2 3
clim = [-10 130]; % display all images on same gray scale
im(1, y, clim)
drawnow


delta = 5; % small compared to max(y) - min(y).  todo: vary if time permits
dpot = @(t) t ./ sqrt(1 + (t / delta).^2); % derivative of hyperbola potential
niter = 200; % max # of iter
tol = 1e-5; % stop tolerance for change in x

%% try several values of the regularization parameter

sig = 1; % we don't know the noise level, so just use sigma=1

reglist = 2 .^ [-5:-2]; % todo: adjust to find a reasonable range
nr = numel(reglist);
gcv = nan(nr,1);
rss = nan(nr,1);
gcv_best = inf;
x_best = [];
ir_best = [];

% w = 2 * (rand(nx,ny) > 0.5) - 1; % bernoulli +/- 1
% weps = 1; % todo: find a reasonable setting

for ir=1:nr
	reg = reglist(ir);
	step = 1 / (1 + reg * 8); % todo: is this appropriate step size?

	xinit = y;
	x = gradient_descent1(A, y, xinit, reg, C, dpot, step, tol, niter);

	rss(ir) = norm(col(A*x) - y(:))^2 / sig^2; % RSS
	gcv(ir) = rss(ir); % todo: modify this to be the correct GCV formula!

	if gcv(ir) < gcv_best
		ir_best = ir;
		gcv_best = gcv(ir);
		x_best = x;
	end


	im(2, x, clim)
	xlabelf('log_2(\beta) = %4.1f', log2(reg))

	im subplot 3
	plot(log2(reglist), rss, '-o')
	titlef 'RSS', xlabelf 'log_2(\beta)'

	im subplot 6
	plot(log2(reglist), gcv, '-o')
	titlef 'GCV', xlabelf 'log_2(\beta)'

	drawnow
end

	im(4, x_best, clim)
	xlabelf('log_2(\beta) = %4.1f', log2(reglist(ir_best)))
