function filename = naiveBayes
% set number of bins used
NUMBINS = 3;

% features: dist_cm, freq, ang, error

% Distance Vector (converted from in to cm)
dist_cm = [6 12 24 36]*2.54;
% Frequency Vector (in cycles/fov of captured image)
freq = [1 2 3 4 5 6 7 8 9 10]/0.5;
% Orientation Vector (in degrees where 0 is vertical, 90 is horizontal)
ang = 0:10:90;

[xe,xs,xh,y] = prepareData('NB');
%xh(:,4) = xh(:,4)/100;

% Discretize error feature
% method: find min and max values of error (cielab and scielab separately),
% divide range by 4, creating 4 range buckets

% cielab errors
xe1 = min(xe(:,4));
xe7 = max(xe(:,4));
xe_incr = round((xe7 - xe1)/6);    % Calculate increment for cielab errors
xe2 = xe1+xe_incr;
xe3 = xe2+xe_incr;
xe4 = xe3+xe_incr;
xe5 = xe4+xe_incr;
xe6 = xe5+xe_incr;

% scielab errors
xs1 = min(xs(:,4));
xs7 = max(xs(:,4));
xs_incr = round((xs7 - xs1)/6);     % Calculate increment for cielab errors
xs2 = xs1+xs_incr;
xs3 = xs2+xs_incr;
xs4 = xs3+xs_incr;
xs5 = xs4+xs_incr;
xs6 = xs5+xs_incr;

% hasler metric
xh1 = min(xh(:,4));
xh7 = max(xh(:,4));
xh_incr = round((xh7 - xh1)/6);     % Calculate increment for Hasler metric
xh2 = xh1+xh_incr;
xh3 = xh2+xh_incr;
xh4 = xh3+xh_incr;
xh5 = xh4+xh_incr;
xh6 = xh5+xh_incr;


% Probabilities

numPoints = length(y);

% Probability that a point is in each bin
for b = 1:NUMBINS,
    p_c(b) = length(find(y==(b-1)))/numPoints;
end

% Probabilities for each bin
for x = 1:NUMBINS,
    I = find(y==(x-1));     % indices of values in bin x
    numBin = length(I);     % num total points in bin x
    
    % Probability of specific distance in bin x
    for d = 1:4,
        eval([ 'p_dist_' num2str(x) '(d) = length(find( xe(I,1) == dist_cm(d)))/numBin;']);
        eval([ 'p_dist_' num2str(x) '_s(d) = length(find( xs(I,1) == dist_cm(d)))/numBin;']); % sCielab
        eval([ 'p_dist_' num2str(x) '_h(d) = length(find( xh(I,1) == dist_cm(d)))/numBin;']); % Hasler
    end

    % Probability of specific freq in bin x
    for f = 2:2:20,
        eval([ 'p_freq_' num2str(x) '(f/2) = length(find( xe(I,2) == f))/numBin;']);
        eval([ 'p_freq_' num2str(x) '_s(f/2) = length(find( xs(I,2) == f))/numBin;']); % sCielab
        eval([ 'p_freq_' num2str(x) '_h(f/2) = length(find( xh(I,2) == f))/numBin;']); % Hasler
    end
        
    % Probability of specific orientation in bin x
    for a = 0:10:90,
        eval([ 'p_ang_' num2str(x) '((a+10)/10) = length(find( xe(I,3) == a))/numBin;']);
        eval([ 'p_ang_' num2str(x) '_s((a+10)/10) = length(find( xs(I,3) == a))/numBin;']); % sCielab
        eval([ 'p_ang_' num2str(x) '_h((a+10)/10) = length(find( xh(I,3) == a))/numBin;']); % Hasler
    end
    
    % Probability of specific error range in bin 1
    eval([ 'p_err_' num2str(x) '(1) = length(find( xe(I,4) <= xe2))/numBin;']);
    eval([ 'p_err_' num2str(x) '(2) = length(find( (xe(I,4) > xe2) & (xe(I,4) <= xe3)))/numBin;']);
    eval([ 'p_err_' num2str(x) '(3) = length(find( (xe(I,4) > xe3) & (xe(I,4) <= xe4)))/numBin;']);
    eval([ 'p_err_' num2str(x) '(4) = length(find( (xe(I,4) > xe4) & (xe(I,4) <= xe5)))/numBin;']);
    eval([ 'p_err_' num2str(x) '(5) = length(find( (xe(I,4) > xe5) & (xe(I,4) <= xe6)))/numBin;']);
    eval([ 'p_err_' num2str(x) '(6) = length(find( xe(I,4) > xe6))/numBin;']);

    % sCielab error values
    eval([ 'p_err_' num2str(x) '_s(1) = length(find( xs(I,4) <= xs2))/numBin;']);
    eval([ 'p_err_' num2str(x) '_s(2) = length(find( (xs(I,4) > xs2) & (xs(I,4) <= xs3)))/numBin;']);
    eval([ 'p_err_' num2str(x) '_s(3) = length(find( (xs(I,4) > xs3) & (xs(I,4) <= xs4)))/numBin;']);
    eval([ 'p_err_' num2str(x) '_s(4) = length(find( (xs(I,4) > xs4) & (xs(I,4) <= xs5)))/numBin;']);
    eval([ 'p_err_' num2str(x) '_s(5) = length(find( (xs(I,4) > xs5) & (xs(I,4) <= xs6)))/numBin;']);
    eval([ 'p_err_' num2str(x) '_s(6) = length(find( xs(I,4) > xs6))/numBin;']);

    % Hasler metric values
    eval([ 'p_err_' num2str(x) '_h(1) = length(find( xh(I,4) <= xh2))/numBin;']);
    eval([ 'p_err_' num2str(x) '_h(2) = length(find( (xh(I,4) > xh2) & (xh(I,4) <= xh3)))/numBin;']);
    eval([ 'p_err_' num2str(x) '_h(3) = length(find( (xh(I,4) > xh3) & (xh(I,4) <= xh4)))/numBin;']);
    eval([ 'p_err_' num2str(x) '_h(4) = length(find( (xh(I,4) > xh4) & (xh(I,4) <= xh5)))/numBin;']);
    eval([ 'p_err_' num2str(x) '_h(5) = length(find( (xh(I,4) > xh5) & (xh(I,4) <= xh6)))/numBin;']);
    eval([ 'p_err_' num2str(x) '_h(6) = length(find( xh(I,4) > xh6))/numBin;']);

end

save NBmodels;

filename = 'NBmodels.mat';


% Make Pretty Plots of classification of training data
% modify the assignment to temp_x to display different vectors--either xe
% (CIELAB error), xs (sCIELAB error), or xh (Hasler metric)
data = {xe xs xh};
type = {'CIELAB' 'sCIELAB' 'hasler'};


for K = 1:3,
    bin = classifyNB(type{K},data{K});

    figure()
    x=data{K};
    index = find(bin==0);
    temp_x = x(index,:);

    % plot the data; change the index into temp_x to see comparisons among
    % different dimensions; temp_x(:,1) is constant so it's useless to plot
    % against it (just ignore it)
    plot(temp_x(:,4), temp_x(:,1), 'bo');
    hold on

    index = find(bin==1);
    temp_x = x(index,:);
    plot(temp_x(:,4), temp_x(:,1),'g^');
    hold on

    index = find(bin==2);
    temp_x = x(index,:);
    plot(temp_x(:,4), temp_x(:,1),'rx');
    hold on

    grid on

    hold off
    title([type{K} ' data for NB classifier']);
    xlabel('Metric');
    ylabel('Distance (cm)');
end