function [hfig, hax] = grt_plotfit(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
sbounds = 1; % if sbounds, plot equal base rate bounds only
lw = 2; %linewidth
mrg = 1; %plot marginals
fit = 0; %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.c1aa fitps.c1ab fitps.c1ba fitps.c1bb];
    dbp(2,:) = [fitps.c2aa fitps.c2ab fitps.c2ba fitps.c2bb];
    dbp(3,:) = [fitps.c3aa fitps.c3ab fitps.c3ba fitps.c3bb];
    dbp(4,:) = [fitps.c4aa fitps.c4ab fitps.c4ba fitps.c4bb];
    dbp(5,:) = [fitps.c5aa fitps.c5ab fitps.c5ba fitps.c5bb];    
else
    dbp(1,:) = [fitps.ax fitps.bx fitps.cx1 fitps.dx1 fitps.sx];
    dbp(2,:) = [fitps.ay fitps.by fitps.cy2 fitps.dy2 fitps.sy];
    dbp(3,:) = [fitps.ax fitps.bx fitps.cx3 fitps.dx3 fitps.sx];
    dbp(4,:) = [fitps.ay fitps.by fitps.cy4 fitps.dy4 fitps.sy];
    dbp(5,:) = [fitps.ax fitps.bx fitps.cx5 fitps.dx5 fitps.sx];
    dbp(6,:) = [fitps.ay fitps.by fitps.cy3 fitps.dy3 fitps.sy];
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
        if sbounds
            plot(x,dbx(2,:),'k','LineWidth',lw)%,'Color',[1 0 0])
            plot(dby(2,:),y,'k','LineWidth',lw)%,'Color',[0 0 1])
        else
            plot(x,dbx(1,:),'k:','LineWidth',lw)
            plot(x,dbx(2,:),'k','LineWidth',lw)
            plot(x,dbx(3,:),'k:','LineWidth',lw)
            plot(dby(1,:),y,'k:','LineWidth',lw)
            plot(dby(2,:),y,'k','LineWidth',lw)
            plot(dby(3,:),y,'k:','LineWidth',lw)
        end
    end

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

if strcmp(xpt,'pxt')
    xlabel('low - [f0] - high')
    ylabel('dull - [Timbre] - sharp')
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