function bayes_mri_plot(data,mri,a_mri,xpt)

% function bayes_mri_plot(data,mri,a_mri,xpt)
%
% nsubj and axis limits defined internally

br = 3; %3

[nothing nsubj] = size(data);
if strcmp(xpt,'pbr')
    xlab_a = '\Delta_{labial}';
    ylab_a = '\Delta_{alveolar}';
else
    xlab_a = '\Delta_{stop}';
    ylab_a = '\Delta_{fricative}';
end
xlab_b = '\Delta_{voiceless}';
ylab_b = '\Delta_{voiced}';
msize_a = 12;
msize_b = 8;
lwid = 1;

cloud = .5*ones(1,3);

lenny = length(mri);
rvec = randperm(lenny);
n_sub = floor(0.1*lenny);
mri_sub = mri(rvec(1:n_sub),:,:,:,:);
a_mri_sub = a_mri(rvec(1:n_sub),:,:,:);

mri_x = [];
mri_y = [];
for s = 1:nsubj
    mri_x = [mri_x; mri(:,1,1,br,s)];
    mri_y = [mri_y; mri(:,2,1,br,s)];
end
x_min = min(mri_x);
x_max = max(mri_x);
y_min = min(mri_y);
y_max = max(mri_y);
x_dif = abs(x_max - x_min);
y_dif = abs(y_max - y_min);
x_min = x_min - .1*x_dif;
x_max = x_max + .1*x_dif;
y_min = y_min - .1*y_dif;
y_max = y_max + .1*y_dif;
x_dif = x_max - x_min;
y_dif = y_max - y_min;
if x_dif > y_dif
    adj = (x_dif - y_dif)/2;
    y_min = y_min-adj;
    y_max = y_max+adj;
else
    adj = (y_dif - x_dif)/2;
    x_min = x_min-adj;
    x_max = x_max+adj;
end
ax = [x_min x_max y_min y_max];
% a_1, a_2
figure
for s = 1:nsubj
    subplot(3,3,s)
    hold on
    axis square
    box on
    axis(ax)
    plot(mri_sub(:,1,1,br,s),mri_sub(:,2,1,br,s),'.',...
        'MarkerSize',msize_a,'MarkerEdgeColor',cloud)
    plot([-1 1],[0 0],'k:','LineWidth',lwid)
    plot([0 0],[-1 1],'k:','LineWidth',lwid)
    hpd_x = cilim(mri(:,1,1,br,s));
    hpd_y = cilim(mri(:,2,1,br,s));
    plot([hpd_x(1) hpd_x(1)],ax(3:4),'k')
    plot([hpd_x(2) hpd_x(2)],ax(3:4),'k')
    plot(ax(1:2),[hpd_y(1) hpd_y(1)],'k')
    plot(ax(1:2),[hpd_y(2) hpd_y(2)],'k')
    plot( ( ( data(s).pM(1,1,br) + data(s).pM(1,2,br) ) - ...
        ( data(s).pM(2,1,br) + data(s).pM(2,2,br) ) ), ...
        ( ( data(s).pM(3,3,br) + data(s).pM(3,4,br) ) - ...
        ( data(s).pM(4,3,br) + data(s).pM(4,4,br) ) ),'k+',...
        'MarkerSize',msize_b)
    plot( ( ( data(s).pM(1,1,br) + data(s).pM(1,2,br) ) - ...
        ( data(s).pM(2,1,br) + data(s).pM(2,2,br) ) ), ...
        ( ( data(s).pM(3,3,br) + data(s).pM(3,4,br) ) - ...
        ( data(s).pM(4,3,br) + data(s).pM(4,4,br) ) ),'ko',...
        'MarkerSize',msize_b)
    text(x_min+.1*x_dif,y_min+.1*y_dif,['S' num2str(s)])
    %text(.58,-.75,axlab(1,:))
    %text(-.9,.75,axlab(2,:))
    if s == 1
        set(gca,'XTick',[],'YTick',[floor(y_max*100)/100])
        ylabel(ylab_a)
    elseif s == 4
        set(gca,'XTick',[],'YTick',[0])
        ylabel(ylab_a)
    elseif s == 8
        set(gca,'XTick',[0],'YTick',[])
        xlabel(xlab_a)
    elseif s == 7
        set(gca,'XTick',[ceil(x_min*100)/100],'YTick',[ceil(y_min*100)/100])
        ylabel(ylab_a)
        xlabel(xlab_a)
    elseif s == 2 || s == 3 || s == 5 || s == 6
        set(gca,'XTick',[],'YTick',[])
    end
end

subplot(3,3,9)
hold on
axis square
box on
axis(ax)
plot(a_mri_sub(:,1,1,br),a_mri_sub(:,2,1,br),'.',...
    'MarkerSize',msize_a,'MarkerEdgeColor',cloud)
plot([-1 1],[0 0],'k:','LineWidth',lwid)
plot([0 0],[-1 1],'k:','LineWidth',lwid)
hpd_a_x = cilim(a_mri(:,1,1,br));
hpd_a_y = cilim(a_mri(:,2,1,br));
plot([hpd_a_x(1) hpd_a_x(1)],ax(3:4),'k')
plot([hpd_a_x(2) hpd_a_x(2)],ax(3:4),'k')
plot(ax(1:2),[hpd_a_y(1) hpd_a_y(1)],'k')
plot(ax(1:2),[hpd_a_y(2) hpd_a_y(2)],'k')
text(x_min+.1*x_dif,y_min+.1*y_dif,'\alpha')
xlabel(xlab_a)

axsize = .26;
axposx = [.15 .375 .6];
axposy = [.105 .405 .705];
for i = 1:9
   subplot(3,3,i)
   hax = gca;
   if i==1
       set(hax,'Position',[axposx(1) axposy(3) axsize axsize])
   elseif i==2
       set(hax,'Position',[axposx(2) axposy(3) axsize axsize])
   elseif i==3
       set(hax,'Position',[axposx(3) axposy(3) axsize axsize])
   elseif i==4
       set(hax,'Position',[axposx(1) axposy(2) axsize axsize])
   elseif i==5
       set(hax,'Position',[axposx(2) axposy(2) axsize axsize])
   elseif i==6
       set(hax,'Position',[axposx(3) axposy(2) axsize axsize])
   elseif i==7
       set(hax,'Position',[axposx(1) axposy(1) axsize axsize])
   elseif i==8
       set(hax,'Position',[axposx(2) axposy(1) axsize axsize])
   elseif i==9
       set(hax,'Position',[axposx(3) axposy(1) axsize axsize])
   end
end

%text(.58,-.75,axlab(1,:))
%text(-.9,.75,axlab(2,:))
%text(.65,.8,'\alpha')
set(gca,'XTick',[floor(x_max*100)/100],'YTick',[])

mri_x = [];
mri_y = [];
for s = 1:nsubj
    mri_x = [mri_x; mri(:,3,1,br,s)];
    mri_y = [mri_y; mri(:,4,1,br,s)];
end
x_min = min(mri_x);
x_max = max(mri_x);
y_min = min(mri_y);
y_max = max(mri_y);
x_dif = abs(x_max - x_min);
y_dif = abs(y_max - y_min);
x_min = x_min - .1*x_dif;
x_max = x_max + .1*x_dif;
y_min = y_min - .1*y_dif;
y_max = y_max + .1*y_dif;
x_dif = x_max - x_min;
y_dif = y_max - y_min;
if x_dif > y_dif
    adj = (x_dif - y_dif)/2;
    y_min = y_min-adj;
    y_max = y_max+adj;
else
    adj = (y_dif - x_dif)/2;
    x_min = x_min-adj;
    x_max = x_max+adj;
end
ax = [x_min x_max y_min y_max];
% b_1,b_2
figure
for s = 1:nsubj
    subplot(3,3,s)
    hold on
    axis square
    box on
    axis(ax)
    plot(mri_sub(:,3,1,br,s),mri_sub(:,4,1,br,s),'.',...
        'MarkerSize',msize_a,'MarkerEdgeColor',cloud)
    plot([-1 1],[0 0],'k:','LineWidth',lwid)
    plot([0 0],[-1 1],'k:','LineWidth',lwid)
    hpd_x = cilim(mri(:,3,1,br,s));
    hpd_y = cilim(mri(:,4,1,br,s));
    plot([hpd_x(1) hpd_x(1)],ax(3:4),'k')
    plot([hpd_x(2) hpd_x(2)],ax(3:4),'k')
    plot(ax(1:2),[hpd_y(1) hpd_y(1)],'k')
    plot(ax(1:2),[hpd_y(2) hpd_y(2)],'k')
    plot( ( ( data(s).pM(1,1,br) + data(s).pM(1,3,br) ) - ...
        ( data(s).pM(3,1,br) + data(s).pM(3,3,br) ) ), ...
        ( ( data(s).pM(2,2,br) + data(s).pM(2,4,br) ) - ...
        ( data(s).pM(4,2,br) + data(s).pM(4,4,br) ) ),'k+',...
        'MarkerSize',msize_b)
    plot( ( ( data(s).pM(1,1,br) + data(s).pM(1,3,br) ) - ...
        ( data(s).pM(3,1,br) + data(s).pM(3,3,br) ) ), ...
        ( ( data(s).pM(2,2,br) + data(s).pM(2,4,br) ) - ...
        ( data(s).pM(4,2,br) + data(s).pM(4,4,br) ) ),'ko',...
        'MarkerSize',msize_b)
    text(x_min+.1*x_dif,y_min+.1*y_dif,['S' num2str(s)])
    %text(.58,-.75,axlab(3,:))
    %text(-.9,.75,axlab(4,:))
    if s == 1
        set(gca,'XTick',[],'YTick',[floor(y_max*100)/100])
        ylabel(ylab_b)
    elseif s == 4
        set(gca,'XTick',[],'YTick',[0])
        ylabel(ylab_b)
    elseif s == 8
        set(gca,'XTick',[0],'YTick',[])
        xlabel(xlab_b)
    elseif s == 7
        set(gca,'XTick',[ceil(x_min*100)/100],'YTick',[ceil(y_min*100)/100])
        ylabel(ylab_b)
        xlabel(xlab_b)
    elseif s == 2 || s == 3 || s == 5 || s == 6
        set(gca,'XTick',[],'YTick',[])
    end
end

subplot(3,3,9)
hold on
axis square
box on
axis(ax)
plot(a_mri_sub(:,3,1,br),a_mri_sub(:,4,1,br),'.',...
    'MarkerSize',msize_a,'MarkerEdgeColor',cloud)
plot([-1 1],[0 0],'k:','LineWidth',lwid)
plot([0 0],[-1 1],'k:','LineWidth',lwid)
hpd_a_x = cilim(a_mri(:,3,1,br));
hpd_a_y = cilim(a_mri(:,4,1,br));
plot([hpd_a_x(1) hpd_a_x(1)],ax(3:4),'k')
plot([hpd_a_x(2) hpd_a_x(2)],ax(3:4),'k')
plot(ax(1:2),[hpd_a_y(1) hpd_a_y(1)],'k')
plot(ax(1:2),[hpd_a_y(2) hpd_a_y(2)],'k')
text(x_min+.1*x_dif,y_min+.1*y_dif,'\alpha')
set(gca,'XTick',[floor(x_max*100)/100],'YTick',[])
xlabel(xlab_b)

axsize = .26;
axposx = [.15 .375 .6];
axposy = [.105 .405 .705];
for i = 1:9
   subplot(3,3,i)
   hax = gca;
   if i==1
       set(hax,'Position',[axposx(1) axposy(3) axsize axsize])
   elseif i==2
       set(hax,'Position',[axposx(2) axposy(3) axsize axsize])
   elseif i==3
       set(hax,'Position',[axposx(3) axposy(3) axsize axsize])
   elseif i==4
       set(hax,'Position',[axposx(1) axposy(2) axsize axsize])
   elseif i==5
       set(hax,'Position',[axposx(2) axposy(2) axsize axsize])
   elseif i==6
       set(hax,'Position',[axposx(3) axposy(2) axsize axsize])
   elseif i==7
       set(hax,'Position',[axposx(1) axposy(1) axsize axsize])
   elseif i==8
       set(hax,'Position',[axposx(2) axposy(1) axsize axsize])
   elseif i==9
       set(hax,'Position',[axposx(3) axposy(1) axsize axsize])
   end
end