%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Author: Ivan Chang
%%   Date: March 2001
%%    v1.1
%%
%% This is the voting function for the evolved neural network.
%% The function NNorgVote takes an evolved neuronal organization and test data
%% as input, and return the classification results.
%%
%% Further work for Spring quarter
%% -------------------------------
%% 1. Add graphical layout routines.
%% 2. Generate evolved code in Matlab or other languages like C.
%% 3. Enhance efficicency, may need to rewrite code in a faster language.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ sonno_class_probabilities, votes, classification, confusion_m, correct, declare ] = ...
   NNorgVote( organization, NLAYER, NNEURON, topology, testdata, Class_index, ...
   B, nclass)

global ...
   unary_kernels ...
   binary_kernels ...
   operand_size ...
   neurons;
   
operand_size = B;
neurons = organization;

%% 2D matrix whose rows represent test sample indices and columns represent
%% SONNOs.  We use the following convention to order the SONNOs:
%%    Odd-numbered columns denote TCOT of neuron ID# (column+1)/2+offset
%%    Even-numbered columns denote TCOTI of neuron ID# (column/2+offset)
%% where offset = number of raw input neurons.
votes = zeros(size(testdata,2),2*NLAYER*NNEURON);

%% 2D matrix representing the class probabilities of SONNOS
sonno_class_probabilities = zeros(2*NLAYER*NNEURON,nclass);

%% Names of kernel functions that are treated as low-level, canonical neural
%% functions.  Ideally they should satisfy closure and completeness (though
%% no need to be orthogonal, i.e. they can depend on each another in Hilbert
%% Space in some non-trivial manner).  Completeness guarantees that the set is
%% sufficient to implement any function, closure guarantees there's no outlier
%% function that cannot be approximated by an appropriate series expansion
%% based on the kernel set.  Current research in neuroscience suggests that
%% many operations have origins in specific molecular dynamics mechanisms,
%% e.g. division can be implemented by so-called shunt inhibition (Mead, Koch).
%% As suggested by Koch, understanding how multiplication is implemented at
%% the molecular level can give us a ground example of how computations
%% implemented by neuronal hardware.
unary_kernels = { 'ident', 'log_2', 'exp_2' };
binary_kernels = { 'f1', 'f2', 'f3', 'f4', 'f5', 'f6', 'f7' };

%% Initialize input neurons with test data
input_neuron_IDs = 1:topology(1,1)-1;
for i=input_neuron_IDs,
   %	added ,2
   neurons{i}.out = num2cell(testdata(i,:),2);
end

%% Use the evolved feedforward network to build sonno_class_probabilities
for i_layer=1:NLAYER,
   for neuron=1:NNEURON,
      neuronID = topology(neuron,i_layer);
      nsc = neurons{neuronID}.nsc;
      if length(nsc) > 0
         nt = sum(sum(nsc,1),2);
         if nt > 0
            sonno_n = 2*(neuronID-max(input_neuron_IDs))-1;
            sonno_class_probabilities(sonno_n:sonno_n+1,:) = nsc/nt;              
         end
      end
   end
end

%% Starting from the input neurons, use the evolved feedforward network to
%% compute the set of voting neurons (those neurons whose SONNOs are
%% activated as a result of processing the test data, SONNOs are threshold
%% control outputs of neurons whose 'sot' or 'soti' are set to 1).
for i_layer=1:NLAYER,
   for neuron=1:NNEURON,
      neuronID = topology(neuron,i_layer);
      
      %% Dummy neurons are not eligible to vote, forget about them.
      if ~strcmp(neurons{neuronID}.func,'')
         %% For this neuron, compute the system output based on selected input
         %% features.
         %	added ,2
         neurons{neuronID}.out = num2cell(computeNeuronFunc(neuronID),2);
      
         %% Determine which subspace TCOT or TCOTI the selected input features
         %% belong: 1 indicates TCOT, 0 indicates TCOTI
         subspace = cat(2,neurons{neuronID}.out{:}) >= neurons{neuronID}.threshold;
      
         %% Is this neuron's TCOT a SONNO?  Yes if 'sot' is 1, no otherwise.
         %% If yes, activate this SONNO by setting the corresponding sample entries
         %% above threshold to 1 in 'votes'.
         if neurons{neuronID}.sot == 1 & confidence_criteria(neuronID)
            votes(:,2*(neuronID-max(input_neuron_IDs))-1) = subspace';
         end
      
         %% Is this neuron's TCOTI a SONNO?  Yes if 'soti' is 1, no otherwise.
         %% If yes, activate this SONNO by setting the corresponding sample entries
         %% above threshold to 1 in 'votes'.
         if neurons{neuronID}.soti == 1 & confidence_criteria(neuronID)
            votes(:,2*(neuronID-max(input_neuron_IDs))) = ~subspace';
         end
      end
   end
end

epsilon = 1e-6;
%% Now the 2D matrix 'votes' have all the voting information for the test samples.
%% Let's compute the final weightings of the votes.  
for i=1:size(testdata,2),
   %% It may be possible that no SONNO is ever activated, in this case, either
   %% they are all not confident enough to vote, or that the presented sample
   %% does not belong to the known class, all unknown classes will be given the
   %% identifier 0 by default.
   if any(votes(i,:))
      nvotes = sum(votes(i,:));
      %% pcc(i,c) = probability of ith SONNO's vote for class c
      pcc = sonno_class_probabilities.*repmat(votes(i,:)',1,nclass);
      %% pccmax(1,c) = maximum probability of all SONNOs voting for class c
      pccmax = max(pcc,[],1);
      %% numerator of fusion function
      numerator = pccmax + (nvotes - sum(pcc,1)).*(1 - pccmax);
      %% denominator of fusion function
      denominator = sum((1./(1 - pcc + epsilon)).*repmat(votes(i,:)',1,nclass),1);
      Wc = 1 - numerator./denominator;
      [ Wcmax identifier ] = max(Wc);
      classification(1,i) = identifier;
   else
      classification(1,i) = 0;
   end
end

%%	test classification results 
	[confusion_m correct declare]=confusion_matrix(classification, Class_index, nclass);
%%	end NNorgVote



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% The ident(x) function simply returns its argument as output.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = ident(x)
z = x;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% The log_2(x) function works like the conventional log to base 2, but also
%% implicitly contains the additonal rectification nonlinearity - the function
%% is turned on only when the signal x (acting both as an input and a control
%% signal) is larger than the threshold value 1, and turned off otherwise.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = log_2(x)
z = log2(max(1,x));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% The exp_2(x) function works like the function 2^x, however it is designed to
%% work within the range of input data, and its output is also confined to
%% within the same range.  If output exceeds the range the most significant
%% bits will be dropped.  The exp_2 function also implicitly assumes
%% rectification nonlinearity so that when input signal falls below certain
%% reference level, inhibition occurs.  Signal above the reference level,
%% on the other hand, excites the function.  The Matlab code below emulates
%% the behavior we just described.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = exp_2(x)
global operand_size;
negative_x_mask = x>=0;
z = 2.^x;
z = z.*negative_x_mask;
frac_z = z - floor(z);
z = bitshift(floor(z),0,operand_size);
z = z + frac_z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% The add(x,y) function works like the conventional function '+', which is
%% linear, except that it's scaled by a factor 1/2 so that the sum of the two
%% input signals is within the dynamic range of B bits (provided inputs are).
%% This function is behaviorally equivalent to linear bidirectional devices of
%% finite gain.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = add(x,y)
z = 0.5*(x+y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% The sub(x,y) function works like the conventional function '-' which is
%% linearly, but implicitly contains the rectification nonlinearity so it
%% behaviorally resembles physical devices like diodes or MOS transistors in
%% diode-connected mode in the linear region.  Conceptually, the sign of (x-y)
%% is a binary control signal which turns off the function when the sign is -ve,
%% on otherwise.  Note that 'sub' is distinct from 'add' we defined above,
%% since they model different devices, the former being nonlinear and the latter
%% linear.  This is clearly not the case for the conventional '-' and '+'
%% functions - one can obtain '-' from '+', or vice versa.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = sub(x,y)
z = max(0,x-y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Various functions that use 'add' and 'sub'
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = f1(x,y)
z = add(x,y);

function z = f2(x,y)
z = sub(x,y);

function z = f3(x,y)
z = sub(y,x);

function z = f4(x,y)
z = add(x,y/2);

function z = f5(x,y)
z = add(y,x/2);

function z = f6(x,y)
z = sub(x,y/2);

function z = f7(x,y)
z = sub(y,x/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Given neuronID, this Matlab function computes the system output from the
%% selected function and feature set.  It'll be used by neurons for voting.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = computeNeuronFunc(neuronID)
global neurons unary_kernels binary_kernels;

func = neurons{neuronID}.func;
features = neurons{neuronID}.features;
for i=1:length(features),
   arg{i} = strcat(strcat('cat(2,neurons{',num2str(features{i})),'}.out{:})');
   func = strrep(func,strcat('@',num2str(i)),arg{i});
end
z = eval(func);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Any additional criteria for activating the SONNOs of neurons are lumped here.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pass = confidence_criteria(neuronID)
global neurons;

pass = 1;
if sum(sum(neurons{neuronID}.nsc,1),2) < 5
   pass = 0;
end



function [confusion_m, correct, declare]=confusion_matrix(classification, Class_index, ncl)
%%	this function verifies classification result by computing the confusion matrix
%%	the confusion matrix has first column corresponding to samples which were not classified
%%	remaining columns correspond to the class index+1

%	find confusion matrix
%	and probability of correct classification
for i=1:ncl
   result=classification(find(Class_index==i));
   total_from_class_i=size(result,2);
   for j=0:ncl
      confusion_m(i,j+1)=size(find(result==j),2)/total_from_class_i;
   end; 
end;
total=size(Class_index,2);
correct=size(find(classification==Class_index),2)/total;
declare=1-size(find(classification==0),2)/total;

