function [marg,pfah,dotr,nrml,npfh,ndot,cy_idx,cx_idx] = condp(conts);

% function [marg,pfah,dotr,nrml,npfh,ndot] = condp(conts);

m = [-.75 .75; 0 0]; % mean vector
c_b = [1 0; 0 1]; % baseline covariance matrix
c_y = [-.5 0 .5]; % horizontal bound
c_x = [-.5 0 .5]; % vertical bound
n_c_y = length(c_y);
n_c_x = length(c_x);

bin = .05;
x = -5:bin:5;
y = -4:bin:4;
lenx = length(x);
leny = length(y);

%r = -.95:bin:.95;
%r = [-.7 -.35 0 .35 .7];
%r = [-.5 .5];
r = -.9:.45:.9;
lenr = length(r);

c = zeros(2,2,2);

n_reg = 2; %number of regions (above/below)
n_dist = 2; %number of marginals/dist
n_prob = 2; %number of probabilities (fa, h)

dist = zeros(leny,lenx,n_dist);
marg = zeros(lenr*lenr,lenx,n_dist,n_reg,n_c_y); % 2,2,3 -> 2 marg, above/below, 3 c_y 
pfah = zeros(lenr*lenr,n_prob,n_reg,n_c_x,n_c_y); % 2,2,3,3 -> fa/h, above/below, 3 c_x, 3 c_y
dotr = zeros(lenr,lenr,n_c_x+1,n_reg,n_c_y); % 4,2,3 -> d' obs/3 c_x, above/below, 3 c_y
nrml = zeros(lenr*lenr,lenx,n_dist,n_reg,n_c_y);
npfh = zeros(lenr*lenr,n_prob,n_reg,n_c_x,n_c_y);
ndot = zeros(lenr,lenr,n_c_x+1,n_reg,n_c_y);


r_idx = 0;
for i = 1:n_c_y
    [buh,cy_idx(i)] = min(abs(c_y(i)-y));
end
for i = 1:n_c_x
    [feh,cx_idx(i)] = min(abs(c_x(i)-x));
end

zd = (1/sqrt(2*pi))*exp(-.5*(x.^2));
zc = zeros(1,lenx);
for i = 1:lenx
    zc(i) = bin*sum(zd(1:i));
end
if conts
    figure
end
for ra_idx = 1:lenr
    for rb_idx = 1:lenr
        r_idx = r_idx+1;
        disp([num2str(r_idx) ' of ' num2str(lenr*lenr)])
        c(:,:,1) = c_b + [0 r(ra_idx); r(ra_idx) 0];
        c(:,:,2) = c_b + [0 r(rb_idx); r(rb_idx) 0];
        dist = zeros(leny,lenx,2);
        for id = 1:2
            rdcov = sqrt(det(c(:,:,id))); % for normalizing
            incov = inv(c(:,:,id));
            moo = m(:,id);
            for iy = leny:-1:1
                for ix = 1:lenx
                    xy = [x(ix); y(iy)];
                    zz = xy-moo;
                    dist(iy,ix,id) = (1/(2*pi*rdcov)) * exp(-0.5*zz'*incov*zz);
                end
            end
        end
        
        if conts
            level = [.005 .005];
            subplot(lenr,lenr,r_idx)
            hold on
            axis square
            hax = gca;
            set(hax,'XLim',[-5 5],'YLim',[-5 5],'XTick',[],'YTick',[]);%'%axis([-5 5 -4 4])
            contour(x,y,dist(:,:,1),level,'EdgeColor',[0 0 0]);
            contour(x,y,dist(:,:,2),level,'EdgeColor',[0 0 0]);
        end

        for ic = 1:n_c_y
            for im = 1:2
                for ix = 1:lenx
                    % above the horizontal decision bound
                    marg(r_idx,ix,im,1,ic) = sum(dist(cy_idx(ic):leny,ix,im));
                    % below the horizontal decision bound
                    marg(r_idx,ix,im,2,ic) = sum(dist(1:cy_idx(ic)-1,ix,im));
                end
            end
        end
        % above the bound
        for ic = 1:n_c_y
            for ir = 1:2
                for im = 1:2
                    marg(r_idx,:,im,ir,ic) = marg(r_idx,:,im,ir,ic)/sum(marg(r_idx,:,im,ir,ic));
                    mm(im) = sum(x.*marg(r_idx,:,im,ir,ic));
                    vm(im) = sum((x.^2).*marg(r_idx,:,im,ir,ic))-mm(im)^2;
                    nrml(r_idx,:,im,ir,ic) = (1/sqrt(2*pi*vm(im)))*exp(-(1/(2*vm(im)))*((x-mm(im)).^2));
                    nrml(r_idx,:,im,ir,ic) = nrml(r_idx,:,im,ir,ic)/sum(nrml(r_idx,:,im,ir,ic));
                end
                dotr(ra_idx,rb_idx,n_c_x+1,ir,ic) = (mm(2)-mm(1))/(.5*(sqrt(vm(1))+sqrt(vm(2))));
            end
        end
        for icy = 1:n_c_y
            for ir = 1:2
                for icx  = 1:n_c_x
                    for ifh = 1:2
                        pfah(r_idx,ifh,ir,icx,icy) = sum(marg(r_idx,cx_idx(icx):lenx,ifh,ir,icy));
                        npfh(r_idx,ifh,ir,icx,icy) = sum(nrml(r_idx,cx_idx(icx):lenx,ifh,ir,icy));
                    end
                    [fm,fi] = min(abs(pfah(r_idx,1,ir,icx,icy)-zc));
                    [hm,hi] = min(abs(pfah(r_idx,2,ir,icx,icy)-zc));
                    dotr(ra_idx,rb_idx,icx,ir,icy) = x(hi)-x(fi);
                    [fm,fi] = min(abs(npfh(r_idx,1,ir,icx,icy)-zc));
                    [hm,hi] = min(abs(npfh(r_idx,2,ir,icx,icy)-zc));
                    ndot(ra_idx,rb_idx,icx,ir,icy) = x(hi)-x(fi);
                end
            end
        end
    end
end