function [hfig, hax] = grt_plotfit_st(fitps,fM,label,dbtype,xpt)

% function [hfig, hax] = grt_plotfit(fitps,fM,label,dbtype,xpt)
%
% fitps = output structure from grtfit
% xpt = string indicating the experiment: pxt, fxd, mbr, pbr
% dbtype = 'linear','piecewise', NOT 'maxpost'
% label = string to appear in top left of main subplot

% n = index to pick out which set of parameters from fitps to use
%       i.e., fitps(1).whatever will be plotted

bounds = 1; % if bounds, plot bounds in addition to contours
lw = 2; %linewidth
mrg = 1; %plot marginals
fit = 1; %plot pred-obs scatter

if length(fitps) > 1
    error('please give me a single set of parameters, please.')
end

m = [fitps(1).maa fitps(1).mab fitps(1).mba fitps(1).mbb];
c = zeros(2,2,4);
c(:,:,1) = fitps(1).caa;
c(:,:,2) = fitps(1).cab;
c(:,:,3) = fitps(1).cba;
c(:,:,4) = fitps(1).cbb;
if strcmp(dbtype,'maxpost')
    dbp(1,:) = [fitps.waa fitps.wab fitps.wba fitps.wbb];    
else
    dbp(1,:) = [0 0 fitps.cx fitps.dx 0];
    dbp(2,:) = [0 0 fitps.cy fitps.dy 0];
end
[nstim,nresp,nbr] = size(fM);

% bin widths
bin = .05;

[dist,x,y] = bivar_norm(m,c,bin,1);

lenx = length(x);
leny = length(y);
xlims = [x(1) x(lenx)];
ylims = [y(1) y(leny)];

if strcmp(dbtype,'maxpost')
    dec = dec_reg(dbtype,dist,x,y,dbp,nbr);
else
    [dbx,dby] = bound_spec(dbp,x,y,dbtype,nbr);
end

level = [.08 .08]; % determines the single contour to plot
minax = -2.5;
maxax = 4.5;

hfig = figure;
if mrg || fit
    subplot('Position',[.3 .3 .65 .65])
end
hax = gca;
pbaspect([1 1 1])
hold on
set(hax,'XTick',[],'YTick',[]);
for id = 1:nstim
    [crap,h] = contour(x,y,dist(:,:,id),level);
    set(h,'EdgeColor',[0 0 0],'LineWidth',lw)
end
plot(m(1,:),m(2,:),'k+','MarkerSize',10)
hax = gca;
set(hax,'XLim', xlims, 'YLim', ylims)

if ~strcmp(label,[])
    text(x(5),y(end-10),label,'FontName','Times','FontSize',16);
end

if strcmp(dbtype,'linear') || strcmp(dbtype,'piecewise')

    if bounds
        plot(x,dbx(2,:),'k','LineWidth',lw)%,'Color',[1 0 0])
        plot(dby(2,:),y,'k','LineWidth',lw)%,'Color',[0 0 1])
    end

elseif strcmp(dbtype,'maxpost')
   
    if bounds
        for is = 1:4
            [crap,h] = contour(x,y,dec(:,:,is,3),[.5 .5]);
            set(h,'EdgeColor',[0 0 0],'LineWidth',lw);
        end
    end
end

if strcmp(xpt,'pxt')
    xlabel('low - [f0] - high')
    ylabel('low - [spectral prominence] - high')
elseif strcmp(xpt,'fxd')
    xlabel('low - [freq.] - high')
    ylabel('short - [dur.] - long')
elseif strcmp(xpt,'mbr')
    text(xlims(1)+abs(.2*xlims(1)),ylims(1)+abs(.2*ylims(1)),'p','FontSize',14)
    text(xlims(1)+abs(.2*xlims(1)),ylims(2)-abs(.1*ylims(2)),'b','FontSize',14)
    text(xlims(2)-abs(.2*xlims(2)),ylims(1)+abs(.2*ylims(1)),'f','FontSize',14)
    text(xlims(2)-abs(.2*xlims(2)),ylims(2)-abs(.1*ylims(2)),'v','FontSize',14)
    xlabel('stop     fricative','FontSize',14)
    ylabel('voiceless    voiced','FontSize',14)
elseif strcmp(xpt,'pbr')
    text(xlims(1)+abs(.2*xlims(1)),ylims(1)+abs(.2*ylims(1)),'p','FontSize',14)
    text(xlims(1)+abs(.2*xlims(1)),ylims(2)-abs(.1*ylims(2)),'b','FontSize',14)
    text(xlims(2)-abs(.2*xlims(2)),ylims(1)+abs(.2*ylims(1)),'t','FontSize',14)
    text(xlims(2)-abs(.2*xlims(2)),ylims(2)-abs(.1*ylims(2)),'d','FontSize',14)
    %text(0,0,'p','FontSize',14)
    xlabel('labial      coronal','FontSize',14)
    ylabel('voiceless    voiced','FontSize',14)
end

% [crap,h] = contour(x,y,dist(:,:,3),level);
% set(h,'EdgeColor',[1 0 0],'LineWidth',lw)
% [crap,h] = contour(x,y,dist(:,:,4),level);
% set(h,'EdgeColor',[0 0 1],'LineWidth',lw)

xtm = m(1,:);
xtv = ones(1,4);
ytm = m(2,:);
ytv = ones(1,4);
for i = 1:4
    xdif = xtm(i) == xtm;
    for j = 1:4
        if xdif(j)==1 & i<j
            xtv(i)=0;
        end
    end
end
for i = 1:4
    ydif = ytm(i) == ytm;
    for j = 1:4
        if ydif(j)==1 & i<j
            ytv(i)=0;
        end
    end
end
xtv = logical(xtv);
ytv = logical(ytv);

if mrg
    aamgx = (1/sqrt(2*pi*c(1,1,1)))*exp(-.5*((x-m(1,1))/sqrt(c(1,1,1))).^2);
    abmgx = (1/sqrt(2*pi*c(1,1,2)))*exp(-.5*((x-m(1,2))/sqrt(c(1,1,2))).^2);
    bamgx = (1/sqrt(2*pi*c(1,1,3)))*exp(-.5*((x-m(1,3))/sqrt(c(1,1,3))).^2);
    bbmgx = (1/sqrt(2*pi*c(1,1,4)))*exp(-.5*((x-m(1,4))/sqrt(c(1,1,4))).^2);
    subplot('Position',[.3 .04 .65 .2])
    hold on
    plot(x,aamgx,'k','LineWidth',lw)%,'Color',[1 0 0])
    plot(x,abmgx,'k--','LineWidth',lw)%,'Color',[0 0 1])
    plot(x,bamgx,'k','LineWidth',lw)%,'Color',[1 0 0])
    plot(x,bbmgx,'k--','LineWidth',lw)%,'Color',[0 0 1])
    hax = gca;
    set(hax,'XLim',xlims,'YLim',[0 .7])
    %set(hax,'XTick', sort(xtm(xtv)), 'YTick', [])
    set(hax,'XTick', [], 'YTick', [])
    pbaspect([(.65/.2) 1 1])
    %daspect([1 1 1])

    aamgy = (1/sqrt(2*pi*c(2,2,1)))*exp(-.5*((y-m(2,1))/sqrt(c(2,2,1))).^2);
    abmgy = (1/sqrt(2*pi*c(2,2,2)))*exp(-.5*((y-m(2,2))/sqrt(c(2,2,2))).^2);
    bamgy = (1/sqrt(2*pi*c(2,2,3)))*exp(-.5*((y-m(2,3))/sqrt(c(2,2,3))).^2);
    bbmgy = (1/sqrt(2*pi*c(2,2,4)))*exp(-.5*((y-m(2,4))/sqrt(c(2,2,4))).^2);
    subplot('Position',[.1 .3 .2 .65])
    hold on
    plot(aamgy,y,'k','LineWidth',lw)%,'Color',[1 0 0])
    plot(abmgy,y,'k','LineWidth',lw)%,'Color',[1 0 0])
    plot(bamgy,y,'k--','LineWidth',lw)%,'Color',[0 0 1])
    plot(bbmgy,y,'k--','LineWidth',lw)%,'Color',[0 0 1])
    hax = gca;
    set(hax,'XLim',[0 .7],'YLim',ylims)
    %set(hax,'XTick',[], 'YTick', sort(ytm(ytv)))
    set(hax,'XTick', [], 'YTick', [])
    pbaspect([1 (.65/.2) 1])
end

if fit
    if ~sum(size(fM)==0)
        subplot('Position',[.1 .04 .2 .2])
        predobs(fitps,fM,dbtype);
    end
end

surfplots = 0;
if surfplots
    figure
    hold on
    for id = 1:4
        surf(x,y,dist(:,:,id))
        shading interp
    end
    hax = gca;
    set(hax,'Zlim', [0 .2])
    set(hax,'XTick',[m(1,1) m(1,3) m(1,4)],'YTick', [m(2,1) m(2,3) m(2,2)])
 %   clim_up = max(max(max(dist)));
 %   clim_dn = 0;
%    set(hax,'CLim',[clim_dn clim_up])
    view([-45 45])
    
    figure
    hold on
    dist2 = level(1)*ones(leny,lenx,4);
    for id = 1:4
        for iy = 1:leny
            for ix = 1:lenx
                if dist(iy,ix,id) <= level(1)
                    dist2(iy,ix,id) = dist(iy,ix,id);
                end
            end
        end
    end
    for id = 1:4
        surf(x,y,dist2(:,:,id))
        shading interp
    end
    hax = gca;
    set(hax,'ZLim', [0 .2])
    set(hax,'XTick',[m(1,1) m(1,3) m(1,4)],'YTick', [m(2,1) m(2,3) m(2,2)])
 %   set(hax,'CLim',[clim_dn clim_up])
    view([-45 45])
end