
% Generate synthetic data from a slice to test vbglmar_slice

close all;
clear all;
clear slice;

% Set main parameters
update_F=0;
T=400;
maxits=32;
xdim=16;
slice.p=4;
priors.WA='Spatial - GMRF';
%priors.overA='Slice';
%priors.overA='Voxel';
priors.overA='None';

% Set up design matrix
X=[boxcars(T,1,5), ones(T,1)];
slice.X=X;
slice=spm_vb_init_volume(slice.X,slice.p);

show_Fthetas=0;
ydim=xdim;

% Set up Laplacian
xords = [1:xdim]'*ones(1,ydim); xords = xords(:)';  % plane X coordinates
yords = ones(xdim,1)*[1:ydim];  yords = yords(:)';  % plane Y coordinates
zords   = ones(xdim*ydim,1)';	%-plane Z coordinates
xyz   = [xords; yords; zords];		%-voxel coordinates
vxyz  = spm_vb_neighbors(xyz');


% Generate regression coefficients from Laplace prior
log_alpha=0;
%log_alpha=-3;
alpha=exp(log_alpha);

from_disk=0;
if from_disk
    load w_mean
else
    %r=randn(xdim*ydim,1);
    %z=r/(sqrt(alpha));
    %w=Linv*z;
    %w_mean=reshape(w,xdim,ydim);
    %disp(sprintf('alpha=%1.2f, mean(z)=%1.2f, std(z)=%1.2f, mean(w)=%1.2f, std(w)=%1.2f',alpha,mean(z),std(z),mean(w),std(w)));
    w_mean=0.5*ones(xdim,ydim);
end


% Generate fMRI time series
n=1;
lambda=1;
if ~from_disk
    for y=1:ydim,
        for x=1:xdim,
            %a1_mean(x,y)=(x+y-1)/(xdim+ydim);
            a1_mean(x,y)=0.5;
            zz=randn(T,1)*sqrt(1/lambda);
            ee(1)=zz(1);
            for t=2:T,
                ee(t,:)=a1_mean(x,y)*ee(t-1,:)+zz(t);
            end
            Y(:,n) = X*[w_mean(x,y);w_mean(x,y)]+ee;
            n=n+1;
        end
    end
end

true_a(:,:,1)=a1_mean;
for pp=2:slice.p,
    true_a(:,:,pp)=zeros(xdim,ydim);
end
% Info about true regression coefficients
cmax=max(w_mean(:));
cmin=min(w_mean(:));

cmax=1;
cmin=0;

% if cmax==cmin
%     cmax=2*cmin;
%     cmin=cmin/2;
% end
w_true_vec=[w_mean(:)',w_mean(:)'];
slice.w_true_vec=w_true_vec(:);

% Declare what parameters are to be optimised
slice.update_w=1;
slice=spm_vb_set_priors(slice,priors,vxyz);
slice.update_lambda=1;
slice.mean_lambda=lambda*ones(xdim*ydim,1); 
slice.update_F=update_F;
slice.maxits=maxits;

slice=spm_vb_glmar(Y,slice);

% display some results
figure
subplot(2,2,1);
imagesc(w_mean,[cmin cmax]);
colorbar
title('True w1_{mean}');
%colormap gray

subplot(2,2,2);
w_img=reshape(slice.wk_ols(1,:),xdim,ydim);
imagesc(w_img,[cmin cmax]);
colorbar
title('OLS w1_{mean}');
%colormap gray

subplot(2,2,3);
w_img=reshape(slice.wk_mean(1,:),xdim,ydim);
imagesc(w_img,[cmin cmax]);
colorbar
title('VB w1_{mean}');
%colormap gray
disp(' ');

e_ols=slice.wk_ols(:)-slice.w_true_vec;
e_vb=slice.wk_mean(:)-slice.w_true_vec;
ols_err=e_ols'*e_ols;
vb_err=e_vb'*e_vb;
p_err=100*(1-vb_err/ols_err);
disp(' ');
disp('Regression coefficients:');
disp(sprintf('Percentage reduction in error = %1.2f', p_err));
disp(' ');
disp('');

disp(sprintf('True alpha = %1.2f', alpha));
for kk=1:slice.k,
    disp(sprintf('Regressor %d',kk));
    disp(sprintf('Estimated alpha = %1.2f', slice.mean_alpha(kk)));
    disp('');
    %disp(sprintf('Alpha from true data %1.2f',slice.mean_alpha_true(kk)));
    disp(sprintf('Number of coefficients %d',slice.N));
    disp(sprintf('Number of data-dependent coeffs = %1.2f',slice.gamma_tot(kk)));
end

% Plot F as a function of VB iteration
figure;plot(slice.F_record);
title('F');

for pp=1:slice.p,
    
    true=squeeze(true_a(:,:,pp));
    true_vec=true(:);
    
    figure
    subplot(2,2,1);
    imagesc(true,[-1 1]);
    colorbar
    title('True');
    subplot(2,2,2);
    a1_img=reshape(slice.ap_ols(pp,:),xdim,ydim);
    imagesc(a1_img,[-1 1]);
    colorbar
    title('OLS');subplot(3,1,3);
    subplot(2,2,3);
    a1_img=reshape(slice.ap_mean(pp,:),xdim,ydim);
    imagesc(a1_img,[-1 1]);
    colorbar
    title('VB');
    
    disp(' ');
    disp(sprintf('Error estimating AR(%d) coefficient',pp));
    a1_ols_err=sum((slice.ap_ols(pp,:)'-true_vec).^2);
    disp(sprintf('OLS: %1.2f',a1_ols_err));
    a1_vb_err=sum((slice.ap_mean(pp,:)'-true_vec).^2);
    disp(sprintf('VB: %1.2f',a1_vb_err));
    disp(' ');
    
end

