%
%
% Robots - Copyright 2011, Profs. W. R. Ashurst and T. D. Placek.
% All rights reserved.
% Auburn University
% Material for CHEN 3600 - SU11
%
%

clear all; clc; clf; close all force;
format compact
format long

tic

N_a = 500;
N_b = 500;
Sig_a = 1.0;
Sig_b = 1.0;
Mu_a = 10;
Mu_b = 11;
n_bins = 50;
Value_a = 20;
Value_b = 1;

crit = 11.0;

Sample_A = normrnd(Mu_a, Sig_a, 1, N_a);
Sample_B = normrnd(Mu_b, Sig_b, 1, N_b);
N_aa = 0;
N_bb = 0;
Sample_AA = zeros(N_a+N_b,2);
Sample_BB = zeros(N_b+N_a,2);
Actual_total_value = N_a*Value_a + N_b*Value_b;

for i = 1:N_a
if Sample_A(i) < crit
N_aa = N_aa + 1;
Sample_AA(N_aa,1) = Sample_A(i); % correct judgement
else
N_bb = N_bb + 1;
Sample_BB(N_bb,2) = Sample_A(i); % type I error
end
end

for i = 1:N_b
if Sample_B(i) < crit
N_aa = N_aa + 1;
Sample_AA(N_aa,2) = Sample_B(i); % type II error
else
N_bb = N_bb + 1;
Sample_BB(N_bb,1) = Sample_B(i); % correct judgement
end
end

figure(1);
init_pool = [Sample_A', Sample_B'];
hist(init_pool, n_bins);
xlim([Mu_a - 5*Sig_a, Mu_b + 5*Sig_b]);
xlabel('Mass Bin');
ylabel('Number');
title('Initial Pooled Samples')

figure(2)

subplot(2,2,1)
hist(Sample_A, 50)
xlim([Mu_a - 5*Sig_a, Mu_a + 5*Sig_a])
xlabel('Mass Bin');
ylabel('Number');
title('Initial Barrel A Samples')

subplot(2,2,2)
hist(Sample_B, 50)
xlim([Mu_b - 5*Sig_b, Mu_b + 5*Sig_b])
xlabel('Mass Bin');
ylabel('Number');
title('Initial Barrel B Samples')

subplot(2,2,3)
hist(Sample_AA, 50)
xlim([min(Mu_a - 5*Sig_a, Mu_b - 5*Sig_b),max(Mu_a + 5*Sig_a, Mu_b + 5*Sig_b )])
xlabel('Mass Bin');
ylabel('Number');
title('Contents of Barrel AA')

subplot(2,2,4)
hist(Sample_BB, 50)
xlim([min(Mu_a - 5*Sig_a, Mu_b - 5*Sig_b),max(Mu_a + 5*Sig_a, Mu_b + 5*Sig_b )])
xlabel('Mass Bin');
ylabel('Number');
title('Contents of Barrel BB')

%
% Theoretical Treatment
%

x_min = min(Mu_a, Mu_b) - 4*max(Sig_a, Sig_b);
x_max = max(Mu_a, Mu_b) + 4*max(Sig_a, Sig_b);
x = linspace(x_min, x_max);
pdf_a = normpdf(x, Mu_a, Sig_a);
pdf_b = normpdf(x, Mu_b, Sig_b);

figure(3)
plot(x, pdf_a, x, pdf_b)
title('Theoretical Analysis')
t2x = x(x<=crit);
t1x = x(x>=crit);
t1y = normpdf(t1x, Mu_a, Sig_a);
patch([crit, crit, t1x], [0, normpdf(crit, Mu_a, Sig_a), t1y], 'r','FaceAlpha',0.15);
t2y = normpdf(t2x, Mu_b, Sig_b);
patch([t2x, crit, crit], [t2y, normpdf(crit, Mu_b, Sig_b), 0],'b','FaceAlpha',0.15);

%
% Analysis of Sorting Results:
% Use data structure to allow another calling script to
% run multiple times to track trends in metrics.
%

sorted_sample(1).num_in_a_pile = N_aa;
sorted_sample(1).num_in_b_pile = N_bb;
sorted_sample(1).num_a_correctly_placed = sum(Sample_AA(:,1)>0);
sorted_sample(1).num_b_correctly_placed = sum(Sample_BB(:,1)>0);
sorted_sample(1).num_of_b_in_the_a_pile = sum(Sample_AA(:,2)>0);
sorted_sample(1).num_of_a_in_the_b_pile = sum(Sample_BB(:,2)>0);

sorted_sample(1).recovery_of_a = sum(Sample_AA(:,1)>0)/N_a;
sorted_sample(1).recovery_of_b = sum(Sample_BB(:,1)>0)/N_b;
sorted_sample(1).type_I_error = sum(Sample_BB(:,2)>0)/N_a;
sorted_sample(1).type_II_error = sum(Sample_AA(:,2)>0)/N_b;
sorted_sample(1).valuation_of_a = N_aa*Value_a;
sorted_sample(1).valuation_of_b = N_bb*Value_b;
sorted_sample(1).profit_by_deception = -Actual_total_value + N_aa*Value_a + N_bb*Value_b;
sorted_sample(1).purity_a = sorted_sample(1).num_a_correctly_placed /sorted_sample(1).num_in_a_pile;
sorted_sample(1).purity_a_theoretical = normcdf(crit, Mu_a, Sig_a)/(normcdf(crit, Mu_a, Sig_a)+normcdf(crit, Mu_b, Sig_b));
sorted_sample(1).purity_b = sorted_sample(1).num_b_correctly_placed /sorted_sample(1).num_in_b_pile;
sorted_sample(1).purity_b_theoretical = (1-normcdf(crit, Mu_b, Sig_b))/((1-normcdf(crit, Mu_a, Sig_a))+(1-normcdf(crit, Mu_b, Sig_b)))

toc