%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Author: Ivan Chang
%%   Date: March 2001
%%   v1.1
%%
%% Questions about program
%% -----------------------
%% Email me at ichang@bobcat.ent.ohiou.edu
%%
%% Introduction
%% ------------
%% Matlab program simulating an artificial evolution process that solves
%% the clustering problem approximately.  Since the covering problem [1]
%% is NPC (NP-complete) is a superset of the clustering problems, therefore
%% clustering is NPC as well.
%%
%% Overview
%% --------
%% This program demonstrates an approach in which clustering is an emergent
%% property of a distributed network of simple threshold computing elements.
%% The idea is based on the notion of a "fitness" functional, suitably chosen
%% to represent the clustering problem.  Mathematically the fitness functional
%% has the form f(S;C,T) where S is a covering of training set T and C is the
%% class distribution of T.  Since f must correlate to how well a function
%% can emerge from the network (the resulting clustering function), f depends
%% on the network organization and the goal is thus to optimize f for all possible
%% coverings S that can be represented by a given organization.  For any covering
%% S represented by the given organization (i.e., a solution to our clustering
%% problem), the function f must be chosen such that its value can be computed
%% with little overhead (say in polynomial time with respect to input size),
%% this way we can verify in a tractable manner whether this value is within an
%% acceptable range of values that correspond to "good" solutions, normalized
%% with respect to the given organization.  Since it is impracticable to
%% examine all possible coverings, the approach taken is that each computing
%% element (neuron) attempts to solve a low-dimension problem locally, the local
%% fitness function for each computing element is chosen in the statistical sense:
%% f(S';C',T') = H(C':S';T')/H(C';T') where H(C':S';T') is the mutual entropy of
%% the local distribution of classes C' among the chosen set T' (subset of T) of
%% samples of local features and the distribution of clusters S' among these
%% sample features; H(C';T') is the entropy of the distribution of classes among
%% the chosen samples.  The clustering is thus solved cooperatively by all the
%% resulting active neurons and f(S;C,T) reflects the cooperative nature of the
%% algorithm.  Refer to [6] for more information about the cooperative voting
%% procedure among active neurons.
%% 
%% T = nd*ms training data set (a total of ms nd-dimensional training vectors)
%% C = 1*ms class distribution corresponding to the ms training patterns
%%
%% Free Parameters: leaving these free allow us to experiment with different
%% ---------------  organizations for problem solving.
%%
%%       NLAYER = number of layers of neurons
%%      NNEURON = number of neurons per layer
%% NSUBSPACECTR = number of subspace control inputs per neuron
%%     NNEAREST = number of nearest neighbors the neuron can connect, 1 in 2D
%%                lattice
%%    NNNEAREST = number of next nearest neighbors the neuron can connect,
%%                at most 3 in finite 2D lattice
%%      NREMOTE = number of neighbors whose Manhattan distance > 2 that the
%%                neuron can connect
%%          pss = subspace selection probability
%%          idt = information deficiency threshold value
%%            B = this limits the dynamic range of input and output values in bits
%%
%% References
%% ----------
%% [1] M. R. Garey and C. S. Johnson, Computers and Intractability.
%% [2] C. Mead, Analog VLSI and Neural Systems.
%% [3] C. Koch, Biophysics of Computation: Information Processing in Single Neurons.
%% [4] M. P. Johnson, Evolving Visual Routines, MIT Media Lab M.Sc. Thesis.
%% [5] D. D. Lee and H. S. Seung, Learning the parts of objects by non-negative.
%%     matrix factorization, Nature, 401, 788-791 (1999).
%% [6] Matlab project document.
%% 
%% Example Usage [ Note: typical training time for 2 layers, 3 neurons/layer ~ 3 hrs
%% -------------                 training time for 20 layers, 3 neurons/layer ~ 18 hrs
%%                       on HP 333 MHz PC with virtual memory set to 512 MB or more.
%%
%%                 In general, you need to increase the number of layers and/or number
%%                 of neurons per layer in order to obtain sufficient number of active
%%                 neurons for the voting process for complex data-class distributions.
%%                 The parameter idt controls misclassification and indecision rates,
%%                 the smaller the idt, the smaller these rates are, idt ~ 0.01 should
%%                 be sufficient provided we use large number of layers, recommend
%%                 trying 30 or 40 layers with 5 neurons/layer.
%%
%%                 To see a classification example, follow these URLs:
%%                    http://www.ent.ohiou.edu/~ichang/EECS715/ez_train_1layer.jpg
%%                    http://www.ent.ohiou.edu/~ichang/EECS715/ez_train.htm
%%
%%                 I am currently investigating the structure of the function and
%%                 threshold selection process to see if an active selection process
%%                 could be implemented to speed up the learning, rather than the
%%                 exhaustive search as done currently.  This program will be ported
%%                 to C and large scale Monte Carlo simulation conducted on the design
%%                 parameters to search for the best design.  This will form part of
%%                 my thesis work. ]
%%
%% For training data set 0 (training0.mat), do this
%%		load training0.mat;
%%    td0 = Training_data;
%%    c0 = Class_index_training;
%%    [ neurons_0 topology_0 cpu_time_0 ] = selforgNN(td0,c0,3,5,1,2,2,2,0.3,0.01,8);
%%
%% For training data set 1 (training1.mat), do this
%%		load training1.mat;
%%    td1 = Training_data;
%%    c1 = Class_index_training;
%%    [ neurons_1 topology_1 cpu_time_1 ] = selforgNN(td1,c1,3,5,1,2,2,2,0.3,0.01,8);
%%
%% For training data set 2 (training2.mat), do this
%%		load training2.mat;
%%    td2 = Training_data;
%%    c2 = Class_index_training;
%%    [ neurons_2 topology_2 cpu_time_2 ] = selforgNN(td2,c2,5,20,1,2,2,2,0.3,0.01,8);
%%
%% For training data set 3 (training3.mat), do this
%%		load training3.mat;
%%    td3 = Training_data;
%%    c3 = Class_index_training;
%%    [ neurons_3 topology_3 cpu_time_3 ] = selforgNN(td3,c3,5,20,1,2,2,2,0.3,0.01,8);
%%
%% For training data set 4 (training4.mat), do this
%%		load training4.mat;
%%    td4 = Training_data;
%%    c4 = Class_index_training;
%%    [ neurons_4 topology_4 cpu_time_4 ] = selforgNN(td4,c4,5,20,1,2,2,2,0.3,0.01,8);
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ neurons, twoD_lattice, cpu_time ] = selforgNN( ...
   T, ...
   C, ...
   NNEURON, ...
   NLAYER, ...
   NNEAREST, ...
   NNNEAREST, ...
   NREMOTE, ...
   NSUBSPACECTR, ...
   pss, ...
   idt, ...
   voting_t, ...
   B) 

global ...
   classdist ...
   neurons ...
   operand_size ...
   unary_kernels ...
   binary_kernels ...
   tco_index_set ...
   pick_tco_index_set ...
   tco_marker ...
   id_threshold ...
   voting_threshold ...
   FID;
   
%% Create a log, change string to path where log file is to be created.
mkdirStatus = mkdir('c:\temp','selforgData');
if mkdirStatus == 1
   fprintf('c:\\temp\\selforgData created successfully\n\n');
else
   fprintf('c:\\temp\\selforgData already exists\n\n');
end
fprintf('Creating log c:\\temp\\selforgData\\sim1.log ...\n\n');
FID = fopen('c:\temp\selforgData\sim1.log','w');
%fprintf('Creating log c:\\temp\\selforgData\\results.log ...\n\n');
%FID2 = fopen('c:\temp\selforgData\results.log','w');
neurons = {};
operand_size = B;
classdist = C;
id_threshold = idt;
voting_threshold= voting_t;

%% 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 uniformly approximated by a series expansion
%% based on the kernel set.  Current research in neuroscience suggests that
%% neural computations have origins in specific molecular dynamics mechanisms,
%% e.g. division or gain normalization can be implemented by shunt inhibition
%% [2, 3].  As suggested by Koch, understanding how multiplication is implemented
%% at the molecular level can give us a ground example of how canonical
%% computations could have been implemented by neuronal hardware.
%%
%% Here only a rather limited set of functions were chosen.  At this point,
%% we do not yet attach any physical meanings to the input arguments of these
%% functions as well as the values returned by them.  This will largely be
%% determined by the context of the input data.  For example, training data set 0
%% may be given in units of length, in that case the data could represent diffuse
%% radar images of a flying saucer being intercepted by a different object.  In
%% this case the arguments and return values of these functions may assume the
%% physical dimension of length.  The evolved network then behaves like spatial
%% filters recognizing complex spatial patterns.  In biological or VLSI-based
%% neural networks, the prevailing physical variables are currents and voltages,
%% and so the assumption is that multimodal real world data should have been
%% transduced by specialized sensors to the common language of currents and
%% voltages, before such network can be used.  Often case, these sensors are
%% themselves special kind of neurons that transduce pressure (e.g. sound),
%% photons into appropriate current and voltage signals.  Associating physical
%% meanings to variables introduce physical constraints that could result in
%% better learning and physically feasible solutions [4, 5].
unary_kernels = { 'ident', 'log_2', 'exp_2' };
binary_kernels = { 'f1', 'f2', 'f3', 'f4', 'f5', 'f6', 'f7' };

%% We will restrict our neurons to a 2D square lattice, which stores the IDs
%% of all the generated neurons (except the raw input neurons), column
%% indices represent layers, row indices represent individual neurons per layer.
twoD_lattice = [];

%% nd = dimension of pattern vector; ms = number of samples
[nd,ms] = size(T)

%% Generate dummy neurons for raw inputs.
%% These are considered 'remote' neurons for all layers except the first.
for neuron_id=1:nd,
   neurons{neuron_id} = make_neuron;
    %	added ,2
   neurons{neuron_id}.out = num2cell(T(neuron_id,:),2);
   neurons{neuron_id}.tco = num2cell(ones(1,ms),2);
   neurons{neuron_id}.tcot = num2cell(ones(1,ms),2);
   neurons{neuron_id}.tcoti = num2cell(ones(1,ms),2);
end

%% Generate 2D lattice of dummy neurons based on design parameters NLAYER, NNEURON.
for i_layer=1:NLAYER
   for neuron=1:NNEURON,
      neuron_id = (i_layer-1)*NNEURON + neuron + nd;
      neurons{neuron_id} = make_neuron;
       %	added ,2
      neurons{neuron_id}.out = num2cell(zeros(1,ms),2);
      twoD_lattice(neuron,i_layer) = neuron_id;
   end
end

%% Now choose a feedforward organization for the 2-dimensional lattice of neurons,
%% with preference given to nearest neighbor connectivity.  In 2D lattice, we use
%% Manhattan distance as a distance measure for any given pair of neurons.  Due to
%% the discrete nature of the lattice, it is appropriate to use the following
%% rules to pick nearest neighbors, the rules are based on three design parameters
%% NNEAREST, NNNEAREST, NREMOTE > 0, instead of modeling the choice through the use
%% of pdf.  Essentially, the parameters allow us to explore pdfs of different shapes.
%%
%% 1. We define NNEAREST, NNNEAREST, NREMOTE > 0 as the number of nearest
%%    neighbors, number of next nearest neighbors, and number from the rest of
%%    remaining neurons that each neuron must connect to.
%% 2. In 2D lattice, there is only 1 nearest neighbor from previous layers,
%%    this neighbor will always be picked
%% 3. In finite 2D lattice, there are at most 3 available next nearest neighbors.
%%    If NNNEAREST is < 3, neighbors will be drawn randomly, one at a time,
%%    from the available neighbors; otherwise all available neighbors are picked.
%% 4. All neurons with distance > 2 are considered remote, connections to these
%%    neurons are made the same way as next nearest neighbors based on parameter
%%    NREMOTE.
%% 5. No lateral connections are allowed, i.e. no connection allowed between neurons
%%    in the same layer

for i_layer=1:NLAYER,
   %% first layer is directly connected to raw signal inputs
   if i_layer == 1
      nearest_id = 1:nd;
      avail_scf_neuron_id = 1:nd;
      
      for neuron=1:NNEURON,
         %% generate input organization psuedo-randomly
         
         %% get neuron ID
         neuron_id = twoD_lattice(neuron,i_layer);
         
         pick_nearest_id = randperm(nd);
         pick_nearest_id = randperm(nd);
         pick_nearest_id = randperm(nd);
         pick_nearest_id = randperm(nd);
         for i=1:min(1+NNNEAREST+NREMOTE,nd),
            n = length(neurons{neuron_id}.presynaptic_neurons);
            neurons{neuron_id}.presynaptic_neurons{n+1} = ...
               nearest_id(pick_nearest_id(i));
            m = length(neurons{nearest_id(pick_nearest_id(i))}.postsynaptic_neurons);
            neurons{nearest_id(pick_nearest_id(i))}.postsynaptic_neurons{m+1} = neuron_id;
         end
         
         %% generate subspace control organization pseudo-randomly        
         pick_avail_scf_neuron_id = randperm(nd);
         pick_avail_scf_neuron_id = randperm(nd);
         pick_avail_scf_neuron_id = randperm(nd);
         pick_avail_scf_neuron_id = randperm(nd);
         for i=1:min(NSUBSPACECTR,nd),
            n = length(neurons{neuron_id}.scf_neurons);
            neurons{neuron_id}.scf_neurons{n+1} = ...
               avail_scf_neuron_id(pick_avail_scf_neuron_id(i));
         end
         
      end
   else
      for neuron=1:NNEURON,
         %% Generate input organization psuedo-randomly.
        
         %% Get neuron ID
         neuron_id = twoD_lattice(neuron,i_layer);
         
         %% Get IDs of NEAREST neighbors, there is only 1 in 2D lattice.
         nearest_id = twoD_lattice(neuron,i_layer-1);
         n = length(neurons{neuron_id}.presynaptic_neurons);
         neurons{neuron_id}.presynaptic_neurons{n+1} = nearest_id;
         m = length(neurons{nearest_id}.postsynaptic_neurons);
         neurons{nearest_id}.postsynaptic_neurons{m+1} = neuron_id;
         
         %% Get IDs of NEXT NEAREST neighbors, at most 3 in finite 2D lattice.
         if neuron == 1 | neuron == NNEURON
            if i_layer > 2
               if neuron == 1
                  nnearest_id = [ twoD_lattice(neuron,i_layer-2), ...
                        twoD_lattice(neuron+1,i_layer-1) ];
               else
                  nnearest_id = [ twoD_lattice(neuron,i_layer-2), ...
                        twoD_lattice(neuron-1,i_layer-1) ];
               end
            else
               if neuron == 1
                  nnearest_id = [ twoD_lattice(neuron+1,i_layer-1) ];
               else
                  nnearest_id = [ twoD_lattice(neuron-1,i_layer-1) ];
               end
            end
         else
            if i_layer > 2
               nnearest_id = [ twoD_lattice(neuron,i_layer-2), ...
                     twoD_lattice(neuron-1,i_layer-1), ...
                     twoD_lattice(neuron+1,i_layer-1) ];
            else
               nnearest_id = [ twoD_lattice(neuron-1,i_layer-1), ...
                     twoD_lattice(neuron+1,i_layer-1) ];
            end
         end
         
         %% Draw items from nnearest_id.
         %% Randomly permute the set of next nearest IDs couple times.
         pick_nearest_id = randperm(length(nnearest_id));
         pick_nearest_id = randperm(length(nnearest_id));
         pick_nearest_id = randperm(length(nnearest_id));
         pick_nearest_id = randperm(length(nnearest_id));
         %% pick items according to the permuted IDs
         for i=1:min(NNNEAREST,length(nnearest_id)),
            n = length(neurons{neuron_id}.presynaptic_neurons);
            neurons{neuron_id}.presynaptic_neurons{n+1} = nnearest_id(pick_nearest_id(i));
            m = length(neurons{nnearest_id(pick_nearest_id(i))}.postsynaptic_neurons);
            neurons{nnearest_id(pick_nearest_id(i))}.postsynaptic_neurons{m+1} = neuron_id;
         end
         
         %% Lastly, make connections to the rest of those neighbors with distance > 2.
         %% These simply belong to the set of all IDs generated in previous layers
         %% UNION input neuron IDs 1:nd MINUS the set of IDs in
         %% (nearest_id UNION nnearest_id)
         remote_id = setdiff(1:nd+(i_layer-1)*NNEURON,union(nearest_id,nnearest_id));
         %% Ready to pick neighbors from remote IDs.
         %% Draw items from remote_id.
         %% Randomly permute the set of remote IDs couple times.
         pick_remote_id = randperm(length(remote_id));
         pick_remote_id = randperm(length(remote_id));
         pick_remote_id = randperm(length(remote_id));
         pick_remote_id = randperm(length(remote_id));
         %% Pick items according to permuted IDs
         for i=1:min(NREMOTE+max(0,NNNEAREST-length(nnearest_id)),length(remote_id)),
            n = length(neurons{neuron_id}.presynaptic_neurons);
            neurons{neuron_id}.presynaptic_neurons{n+1} = remote_id(pick_remote_id(i));
            m = length(neurons{remote_id(pick_remote_id(i))}.postsynaptic_neurons);
            neurons{remote_id(pick_remote_id(i))}.postsynaptic_neurons{m+1} = neuron_id;
         end
         
         %% Generate subspace control organization pseudo-randomly.
         avail_scf_neuron_id = 1:nd+(i_layer-1)*NNEURON;
         pick_avail_scf_neuron_id = randperm(length(avail_scf_neuron_id));
         pick_avail_scf_neuron_id = randperm(length(avail_scf_neuron_id));
         pick_avail_scf_neuron_id = randperm(length(avail_scf_neuron_id));
         pick_avail_scf_neuron_id = randperm(length(avail_scf_neuron_id));
         for i=1:min(NSUBSPACECTR,length(avail_scf_neuron_id)),
            n = length(neurons{neuron_id}.scf_neurons);
            neurons{neuron_id}.scf_neurons{n+1} = ...
               avail_scf_neuron_id(pick_avail_scf_neuron_id(i));
         end

      end
   end
end
%% Done prewiring all neurons
%% Write prewired organization to log file.
for i_layer=1:NLAYER,
   for neuron=1:NNEURON,
      neuronID = twoD_lattice(neuron,i_layer);
      %% displayNNInfo(neuronID);
      fprintf(FID,'Neuron ID # %d (layer %d, %d-th neuron)\n', ...
         neuronID, i_layer, neuron);
      fprintf(FID,'   Presynaptic neurons: ');
      for j=1:length(neurons{neuronID}.presynaptic_neurons),
         fprintf(FID,'%d ',neurons{neuronID}.presynaptic_neurons{j});
      end
      fprintf(FID,'\n\n');
   end
end
%%	Plot prewired organization
figure;
plot_prewired_structure(NLAYER,NLAYER,NNEURON,NNEURON,neurons);


%% Subspace learning starts here, this is a supervised learning procedure
%% assuming that the class information is given.  Subspace learning is based
%% on local optimization of mutual entropy between the class and the subspace,
%% the subspace learning concept is equivalent to the concept of clustering, that
%% is why it is also called self organization.  In general, subspace learning
%% idea is to explore the best local transformation (restricted to a given set
%% of kernel functions, their linear combinations and composites) such that
%% classes are optimally separated into subspaces.  Since we are using the
%% optimal local information index, the subspace learned is the best possible
%% one given the available kernel functions (other sets of kernel functions
%% might give even better optimization, but the kernel set we used should
%% be simple enough for VLSI implementation).
figure;
hold on;
%% create master list of 1000 or so TCO indices such that TCOT or TCOTI occurs
%% with probability pss
%% TCO = 1
%% TCOT = 2
%% TCOTI = 3
tco_index_set = ...
   [ ones(1,1000-2*floor(pss*1000)),ones(1,floor(pss*1000))+1,ones(1,floor(pss*1000))+2 ];
pick_tco_index_set = randperm(length(tco_index_set));
pick_tco_index_set = randperm(length(tco_index_set));
tco_marker = 1;

%% start recording the time for training
start = cputime;

   % The learning proceeds as follows:
   % Starting from the first layer of neurons, pick the neuron with the lowest
   % ID value and proceed in order until all neurons in the layer are exhausted.
   % We then proceed to the second layer of neurons and repeat the previous
   % process for the second layer.  Continue through the rest of the layers.
   % Each neuron runs the entire epoch of training samples presented to it with
   % clock selection, before giving up its job to some other neurons.
   
   % Things to do:
   % Explore feedback organization Spring quarter.
   
   for i_layer=1:NLAYER,
      for neuron=1:NNEURON,
         fprintf('Training: layer %d, neuron %d (ID %d)\n', i_layer, neuron, ...
            twoD_lattice(neuron,i_layer));
         %% This neuron got to select which subspace it is going to learn, and
         %% this is done by keeping the following rules in mind:
         %%
         %% For each neuron ID in 'scf_neurons', choose an SCFN sequentially
         %%
         %%    From the chosen SCFN, select TCO, TCOI, or TCOTI subspaces randomly
         %%    such that the combined probability of TCOI and TCOTI would not
         %%    exceed 2*pss, the selected subspace becomes the 'tci' of this neuron,
         %%    at the same time, transfer the corresponding output information
         %%    deficiency of the selected subspace, i.e.,
         %%    (the corresponding delta_o, delta_ot, or delta_oti from the chosen
         %%    SCFN) and this becomes the input information deficiency 'delta_in'
         %%    of this neuron.
         %%
         %%    After 'tci', 'delta_in' of this neuron are all set, we are ready
         %%    to select the 'features' from among the (1+NNNEAREST+NREMOTE) prewired
         %%    'presynaptic_neurons'.  In this project, all kernel functions are
         %%    unary or binary, hence at most 2 coordinates may be selected at a time
         %%    from each sample pattern vector.  A kernel function is then applied to
         %%    the selected feature and the maximal information reduction in the
         %%	   chosen subspace is computed against all possible thresholds.
         %%
         %%    The above process is repeated for all possible combinations of features
         %%    and all kernel functions.  The maximal information reduction in the
         %%    selected subspace can thus be found, and the corresponding function,
         %%    threshold, selected features, 'tco', output information deficiency
         %%    of input space (delta_tco), 'tcot', 'tcoti', output information
         %%    deficiencies of learned subspaces (delta_tcot, delta_tcoti) as well as
         %%    system 'out' can be updated in neuron's memory.
         %%
         %% The entire indented process is repeated for all SCFNs in 'scf_neurons'
         %%
         %% The function EBE will do all the job described by above.
         EBE(twoD_lattice(neuron,i_layer));
      end
      plot_learned_structure(NLAYER,i_layer,NNEURON,NNEURON,neurons);
      pause(1);
   end

%% cpu_time = total training time
cpu_time = cputime - start
fclose(FID);
%dpyAllNNInfo(neurons, 'c:\temp\selforgData\results.log');

%%	end selforgNN


function EBE(neuronID)
global ...
   neurons ...
   unary_kernels ...
   binary_kernels ...
   tco_index_set ...
   pick_tco_index_set ...
   tco_marker ...
   id_threshold ...
   voting_threshold ...
   FID;

scf_neuronIDs = neurons{neuronID}.scf_neurons;
presynaptic_neuronIDs = neurons{neuronID}.presynaptic_neurons;

%% For each SCFN this neuron (neuronID) connects, get the
%% threshold control input ('tci') for analogue processing, the
%% threshold control input ('tci') models digital selection of
%% the subspace this neuron is going to examine, and is obtained
%% by randomly selecting one of three threshold control outputs
%% 'tco', 'tcot', 'tcoti' from the corresponding SCFN.
for i=1:length(scf_neuronIDs),
   
   clk=tco_index_set(pick_tco_index_set(tco_marker));
   neurons{neuronID}.clock{i} = clk 
   
   if clk == 1
      tci = neurons{scf_neuronIDs{i}}.tco;
      delta_in = neurons{scf_neuronIDs{i}}.delta_o;
      fprintf(FID,'Neuron ID# %d selected TCO from Neuron ID# %d as its TCI, with clock 1 and delta_in is %f\n\n',neuronID,scf_neuronIDs{i},delta_in);
   elseif clk == 2
      tci = neurons{scf_neuronIDs{i}}.tcot;
      delta_in = neurons{scf_neuronIDs{i}}.delta_ot; 
      fprintf(FID,'Neuron ID# %d selected TCOT from Neuron ID# %d as its TCI, with clock 1 and delta_in is %f\n\n',neuronID,scf_neuronIDs{i},delta_in);
   elseif clk == 3
      tci = neurons{scf_neuronIDs{i}}.tcoti;
      delta_in = neurons{scf_neuronIDs{i}}.delta_oti;
      fprintf(FID,'Neuron ID# %d selected TCOTI from Neuron ID# %d as its TCI, with clock 1 and delta_in is %f\n\n',neuronID,scf_neuronIDs{i},delta_in);
   end
   
   %% forward marker for next SCFN
   tco_marker = tco_marker + 1;
   if tco_marker > 1000
      tco_marker = tco_marker - 1000;
   end
   

   %% If the selected subspace and its corresponding information deficiency
   %% is less than or equal to the chosen threshold idt, then not much
   %% information can be gained by further dividing the selected subspace,
   %% save the duty of this neuron from learning the selected subspace so
   %% it can move on to learn other selected subspaces right the way.
   if delta_in > id_threshold
      %% neurons{neuronID}.tci = tci;
      %% neurons{neuronID}.delta_in = delta_in;
      %% neurons{neuronID}.tco = neurons{neuronID}.tci;
      %% neurons{neuronID}.delta_o = neurons{neuronID}.delta_in;
      
      %% Now 'tci' is in place, this neuron begins stage 1 learning.  Stage 1
      %% learning entails the followings:
      if  any(cat(2,tci{:}))
         %% First evolve the best local clustering routine by linear separation
         %% based on available 1D input features, store that information in the
         %% memory of this neuron.
         for j=1:length(presynaptic_neuronIDs),
            %% fprintf(FID,'   Neuron ID# %d applies EBE_aux to output from presynaptic neuron ID# %d\n',neuronID,presynaptic_neuronIDs{j});
            %% fprintf(FID,'System output from neuron ID# %d:\n',presynaptic_neuronIDs{j});
            %% fprintf('Neuron ID# %d applies EBE_aux to output from presynaptic neuron ID# %d\n',neuronID,presynaptic_neuronIDs{j});
            %% fprintf('System output from neuron ID# %d:\n',presynaptic_neuronIDs{j});
            %% pause(4)
            %% fprintf(FID,'%f\n',cat(2,neurons{presynaptic_neuronIDs{j}}.out{:}));
            %	modified May 2001  Janusz Starzyk
            %flag = EBE_aux(cat(2,neurons{presynaptic_neuronIDs{j}}.out{:}), ...
            %   neuronID,delta_in,tci);
            flag = EBE_aux(neurons{presynaptic_neuronIDs{j}}.out{:}, neuronID,delta_in,tci);
            if flag
             	disp('flag set simple 1D features');
               delta_in
               neurons{neuronID}.func = 'feval(unary_kernels{1},@1)';
               neurons{neuronID}.features = { presynaptic_neuronIDs{j} };
               neurons{neuronID}.clock_selected{1} = neurons{neuronID}.clock{i};
               neurons{neuronID}.scf_neuron_selected{1} = scf_neuronIDs{i};
               %	added by Janusz Starzyk May 29 2001
               if neurons{neuronID}.clock{i}==1
             	neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tco;
               neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_o;
            	elseif neurons{neuronID}.clock{i}==2
               neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tcot;
               neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_ot;
               elseif neurons{neuronID}.clock{i}==3
               neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tcoti;
               neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_oti;
               end;
            end
      	end
      
      	%% Next see if the neuron can evolve better local clustering routines
      	%% based on the same 1D features, by exploring more complex separating
      	%% functions.  If a better one can be found, the corresponding
         %% information in the memory of this neuron will be updated accordingly.
         
         %%	modified May 18 2001 Janusz Starzyk
         %%	the purpose is to speed up neuron learning so the nonlinear input functions 
         %%	are selected randomly and search is performed on a limited selection of functions
         %%	the range of selection indexes is the same as in the full search
         
      %	for j=1:length(presynaptic_neuronIDs),
      %      for k=1:length(binary_kernels),
      %         for p=1:length(unary_kernels),
      %            for q=1:length(unary_kernels),
      	count_single=5;
         for counter=1:count_single,
            j=find([presynaptic_neuronIDs{:}]==neurons{neuronID}.features{1});
           % j=1+floor(rand*(length(presynaptic_neuronIDs)-eps));
            k=1+floor(rand*(length(binary_kernels)-eps));
            p=1+floor(rand*(length(unary_kernels)-eps));
            q=1+floor(rand*(length(unary_kernels)-eps));

                     %% fprintf(FID,'   Neuron ID# %d applies EBE_aux to outputs from presynaptic neuron ID# %d with binary kernel %s and unary kernels %s %s\n',neuronID,presynaptic_neuronIDs{j},binary_kernels{k},unary_kernels{p},unary_kernels{q});
                     %% fprintf(FID,'System output from neuron ID# %d:\n',presynaptic_neuronIDs{j});
                     %% fprintf('Neuron ID# %d applies EBE_aux to outputs from presynaptic neuron ID# %d\n',neuronID,presynaptic_neuronIDs{j});
                     %% fprintf('with binary kernel %s and unary kernels %s %s\n',binary_kernels{k},unary_kernels{p},unary_kernels{q});
                     %% fprintf('System output from neuron ID# %d:\n',presynaptic_neuronIDs{j});
                     %% pause(4);
                     %% fprintf(FID,'%f\n',cat(2,neurons{presynaptic_neuronIDs{j}}.out{:}));
                     flag = EBE_aux(feval(binary_kernels{k}, ...
                        feval(unary_kernels{p},cat(2,neurons{presynaptic_neuronIDs{j}}.out{:})), ...
                        feval(unary_kernels{q},cat(2,neurons{presynaptic_neuronIDs{j}}.out{:}))), ...
                        neuronID,delta_in,tci);
%                        feval(unary_kernels{p},cat(2,neurons{j}.out{:})), ...
                        %feval(unary_kernels{q},cat(2,neurons{j}.out{:}))), ...
                        if flag
                     	disp('flag set complex 1D features');
                        delta_in
                        neurons{neuronID}.func = ...
                           strcat(strcat(strcat(strcat(strcat(strcat(strcat( ...
                           strcat('feval(binary_kernels{',num2str(k)),'},'), ...
                           'feval(unary_kernels{'),num2str(p)),'},@1),'), ...
                           'feval(unary_kernels{'),num2str(q)),'},@2))');
                        %j
                        neurons{neuronID}.features = ...
                           { presynaptic_neuronIDs{j} presynaptic_neuronIDs{j} };
                        neurons{neuronID}.clock_selected{1} = neurons{neuronID}.clock{i};
                        neurons{neuronID}.scf_neuron_selected{1} = scf_neuronIDs{i};
                        %	added by Janusz Starzyk May 29 2001
                			if neurons{neuronID}.clock{i}==1
             				neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tco;
               			neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_o;
            				elseif neurons{neuronID}.clock{i}==2
               			neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tcot;
               			neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_ot;
               			elseif neurons{neuronID}.clock{i}==3
               			neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tcoti;
               			neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_oti;
            				end;
         				end
 			 end;
                  
        %          end
        %       end
        %    end
        % end
      
      	%% So far the neuron has examined only 1D features, one at a time, out
      	%% of its total number of presynaptic inputs.  Higher dimensional features
      	%% may offer better resolution.  So let's see if the neuron can evolve
      	%% still even better local clustering routines based on 2D features, using
      	%% complex separating functions which are composites of unary and binary
         %% kernel functions.
         
         %%	modified May 18 2001 Janusz Starzyk
         %%	the purpose is to speed up neuron learning so the nonlinear input functions 
         %%	are selected randomly and search is performed on a limited selection of functions
         %%	the range of selection indexes is the same as in the full search

      	count_double=10;
         for counter=1:count_double,

            j1=find([presynaptic_neuronIDs{:}]==neurons{neuronID}.features{1});
            j2=1+floor(rand*(length(presynaptic_neuronIDs)-eps));
            if length(presynaptic_neuronIDs)>1
            while j2==j1
               j2=1+floor(rand*(length(presynaptic_neuronIDs)-eps));
            end;
            end;
            %j1=1+floor(rand*(length(presynaptic_neuronIDs)-eps));
            %j2=1+floor(rand*(length(presynaptic_neuronIDs)-eps));
            k=1+floor(rand*(length(binary_kernels)-eps));
            p=1+floor(rand*(length(unary_kernels)-eps));
            q=1+floor(rand*(length(unary_kernels)-eps));
            
         %for j1=1:length(presynaptic_neuronIDs),
         %   for j2=j1+1:length(presynaptic_neuronIDs),

         %         for k=1:length(binary_kernels),
         %            for p=1:length(unary_kernels),
         %               for q=1:length(unary_kernels),
                           %% fprintf(FID,'   Neuron ID# %d applies EBE_aux to outputs from presynaptic neuron ID#s %d %d with binary kernel %s and unary kernels %s %s\n',neuronID,presynaptic_neuronIDs{j1},presynaptic_neuronIDs{j2},binary_kernels{k},unary_kernels{p},unary_kernels{q});
                           %% fprintf(FID,'System output from neuron ID#s %d %d:\n',presynaptic_neuronIDs{j1},presynaptic_neuronIDs{j2});
                           %% fprintf('Neuron ID# %d applies EBE_aux to outputs from presynaptic neuron ID#s %d %d\n',neuronID,presynaptic_neuronIDs{j1},presynaptic_neuronIDs{j2});
                           %% fprintf('with binary kernel %s and unary kernels %s %s\n',binary_kernels{k},unary_kernels{p},unary_kernels{q});
                           %% fprintf('System output from neuron ID#s %d %d:\n',presynaptic_neuronIDs{j1},presynaptic_neuronIDs{j2});
                           %% pause(4);
                           %% fprintf(FID,'%f\n',cat(2,neurons{presynaptic_neuronIDs{j1}}.out{:}));
                           %% fprintf(FID,'%f\n',cat(2,neurons{presynaptic_neuronIDs{j2}}.out{:}));
                           flag = EBE_aux(feval(binary_kernels{k}, ...
										feval(unary_kernels{p},cat(2,neurons{presynaptic_neuronIDs{j1}}.out{:})), ...
                              feval(unary_kernels{q},cat(2,neurons{presynaptic_neuronIDs{j2}}.out{:}))), ...
                              neuronID,delta_in,tci);
                           %% feval(unary_kernels{p},cat(2,neurons{presynaptic_neuronIDs{j1}}.out{:})), ...
                           % feval(unary_kernels{p},cat(2,neurons{j1}.out{:})), ...
                           if flag
                              disp('flag set 2D features');
                              delta_in
                              neurons{neuronID}.func = ...
                                 strcat(strcat(strcat(strcat(strcat(strcat(strcat( ...
                                 strcat('feval(binary_kernels{',num2str(k)),'},'), ...
                                 'feval(unary_kernels{'),num2str(p)),'},@1),'), ...
                                 'feval(unary_kernels{'),num2str(q)),'},@2))');
                              neurons{neuronID}.features = ...
                                 { presynaptic_neuronIDs{j1} presynaptic_neuronIDs{j2} };
                              neurons{neuronID}.clock_selected{1} = neurons{neuronID}.clock{i};
                              neurons{neuronID}.scf_neuron_selected{1} = scf_neuronIDs{i};
                              %	added by Janusz Starzyk May 29 2001
                					if neurons{neuronID}.clock{i}==1
             						neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tco;
               					neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_o;
            						elseif neurons{neuronID}.clock{i}==2
               					neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tcot;
               					neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_ot;
               					elseif neurons{neuronID}.clock{i}==3
               					neurons{neuronID}.tci = neurons{scf_neuronIDs{i}}.tcoti;
               					neurons{neuronID}.delta_in = neurons{scf_neuronIDs{i}}.delta_oti;
                     			end;
                  			end   
         %              end
         %           end
         %        end

         %  end
         %end
         end;
      end
   end %% Done learning selected TCI
end %% Done all SCFNs

%% Now the EBE has examined all SCFNs.  Has the neuron really learnt anything new?
%% It maybe the case that all TCIs that this neuron picked have already been
%% resolved by previous neurons, and hence no significant gain in information
%% deficiency reduction can be achieved anyway, this is the case when delta_in
%% of all selected subspaces are <= id_threshold.  If nothing has been learnt,
%% then this neuron will simply set all its delta's to 0.
if strcmp(neurons{neuronID}.func,'' )
   neurons{neuronID}.delta_in = 0;
   neurons{neuronID}.delta_o = 0;
   neurons{neuronID}.delta_ot = 0;
   neurons{neuronID}.delta_oti = 0;
end;
%%	end EBE



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Auxiliary EBE function for threshold selection.
%% Somehow we need to understand the structure of the problem of threshold
%% selection more in order to find a more efficient algorithm.  
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function flag = EBE_aux(computed_feature,neuronID,selected_delta_in,selected_tci)
global neurons classdist id_threshold voting_threshold threshold_learning_step;

%% set the step for threshold iteration
%% step = 0.05;
max_cid = max(classdist);
delta_s = zeros(2,1);
delta_t = 1;
delta_ti = 1;
flag = 0;
selected_threshold = 0;
selected_nsc = [];

%% small non-zero number to prevent overflowing log2
eps = 1e-12;
tci = find(cat(2,selected_tci{:}));
%		modified April 22, 2001 by Janusz Satrzyk
   mn=mean(computed_feature(tci));
   mx=max(computed_feature(tci));
%	delt is the interval size for the subspace probability
   delt=std((computed_feature(tci)));
   pickxm=([ mn mn-delt mn+delt]);
   dend=0.05*delt;
   ntot=size(computed_feature(tci),2);
   for c=1:max_cid,
     nclass(c)=sum(classdist(tci) == c);  
   end;
%	calculate class probabilities
probc=nclass/ntot;
while delt>dend
   for c=1:max_cid,
      for i=2:3
      prob(c,i) = sum(computed_feature(tci) <= pickxm(i) & classdist(tci) == c);
   	end;   
   end;

nsubs1 = sum(prob);
nsubs2=ntot-nsubs1;
%	calculate probabilities of subspaces
probs1=nsubs1/ntot+eps;
probs2=nsubs2/ntot+eps;
%	calculate class probabilities in the first subspace
probcs1=prob/ntot+eps;
%	calculate class probabilities in the second subspace
probcs2=probc'*ones(1,size(probcs1,2))-probcs1+eps;
%	calculate relative entropy and information index
delentr=sum(probcs1.*log2(probcs1)+ probcs2.*log2(probcs2)) ...
   - probs1.*log2(probs1) - probs2.*log2(probs2);
inf=1-delentr./sum(probc.*log2(probc));
probct=probc;
%	identify iteratively the pick location
%	start with pick position which corresponds 
%	at mean values of the sample points
   % pickxm contains the locations of 3 points near the pick
   % the maximum pick value is selected and its location passed
   % as pickxm(1)
   [maxinf,ind]=max(inf);
   pick1=pickxm(ind);
   prob(:,1)=prob(:,ind);
   delt=delt/2;
   pickxm=[pick1,pick1-delt,pick1+delt];
end;  % while
%%****end modification 

selected_threshold = pickxm(1);

%% Based on the selected threshold, find the selected joint distribution
for i=1:max_cid,
   selected_nsc(1,i) = ...
      sum(computed_feature(tci) >= selected_threshold & classdist(tci) == i);
   selected_nsc(2,i) = ...
      sum(computed_feature(tci) < selected_threshold & classdist(tci) == i);
end

nc = sum(selected_nsc,1);
ns = sum(selected_nsc,2);
nt = sum(sum(selected_nsc,1),2);

neg_avg_class_entropy = sum((nc+eps)/nt.*log2((nc+eps)/nt));
neg_joint_entropy = sum((selected_nsc+eps)/nt.*log2((selected_nsc+eps)/nt),2);
neg_subspace_entropy = (ns+eps)/nt.*log2((ns+eps)/nt);

delta_s = (neg_joint_entropy - neg_subspace_entropy)/neg_avg_class_entropy;

delta_t = delta_s(1,1);
delta_ti = delta_s(2,1);
%selected_nsc;

%% The above maximum information deficiency reduction is compared to the
%% maximum information deficiency reduction that has so far been obtained
%% by this neuron.  If the maximum information deficiency reduction found
%% during the above iteration process is larger, then update this neuron
%% with the new knowledge gained, also return a flag to the EBE and tell
%% it to record the selected input features and the name of the function
%% that results in better separation of the transformed features.
if (delta_t + delta_ti)*selected_delta_in < ...
      neurons{neuronID}.delta_ot + neurons{neuronID}.delta_oti

   neurons{neuronID}.tci = selected_tci;
   neurons{neuronID}.delta_in = selected_delta_in;
   neurons{neuronID}.tco = selected_tci;
   neurons{neuronID}.delta_o = selected_delta_in;
   neurons{neuronID}.delta_ot = delta_t*selected_delta_in;
   neurons{neuronID}.delta_oti = delta_ti*selected_delta_in;
    %	added ,2
   neurons{neuronID}.tcot = ...
      num2cell((computed_feature >= selected_threshold) & ...
   cat(2,selected_tci{:}),2);
    %	added ,2
   neurons{neuronID}.tcoti = ...
      num2cell((computed_feature < selected_threshold) & ...
      cat(2,selected_tci{:}),2);
   %	added ,2
   neurons{neuronID}.out = num2cell(computed_feature,2);
   neurons{neuronID}.nsc = selected_nsc;
   neurons{neuronID}.threshold = selected_threshold;
   if neurons{neuronID}.delta_ot <= voting_threshold
      neurons{neuronID}.sot = 1;
   else
      neurons{neuronID}.sot = 0;
   end
   if neurons{neuronID}.delta_oti <= voting_threshold
      neurons{neuronID}.soti = 1;
   else
      neurons{neuronID}.soti = 0;
   end
   flag = 1;
end
%%	end EBE_aux



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% make_neuron creates a neuron structure that contains all information
%% necessary for the subspace learning algorithm.  Definitions of the fields
%% within this structure are given below:
%%
%% func 						- function that achieves the best local subspace learning
%% threshold 				- margin that gives the best local clustering, format: analog value
%% features 				- best input features for local subspace learning,
%%            					format: { 4 } or { 4 6 } neuron IDs
%% out 						- feature map generated by the selected function
%% presynaptic_neurons 	- the set of neuron IDs prewired for input
%% postsynaptic_neurons - the set of neuron IDs prewired for output
%% scf_neurons 			- the set of neuron IDs prewired for subspace control
%% scf_neuron_selected 	- neuron IDs selected for subspace control
%% clock 					- prewired clock controlling neurons outputs from 1 (tco) to 3 (tcoti)
%% clock_selected 		- clock controlling neuron output from 1 (tco) to 3 (tcoti)
%% tci 						- input subspace selected from 'scf_neurons' for the best clustering,
%%       						format: digit cell { 0 0 1 0 1 ... }
%% delta_in					- input subspace information deficiency
%%	nsc						- output subspaces number of points from each class ( 2x(number of classes) matrix)
%% delta_o 					- output subspace information deficiency defined by delta_o = delta_in
%% tco 						- output subspace = input subspace
%% delta_ot 				- accumulated threshold subspace information deficiency defined by
%%            					delta_ot = delta_in * delta_t
%% tcot 						- threshold output subspace
%% delta_oti 				- accumulated inverted threshold subspace information deficiency
%%             				defined by delta_oti = delta_in * delta_ti
%% tcoti 					- inverted threshold output subspace
%% sot 						- 0 (default) indicates that the threshold output subspace's accumulated
%%       information deficiency is > idt, in this case further subspace learning
%%       is required in order to gain the best mutual entropy between local
%%       subspace distribution and local class distribution, hence maximizing the
%%       information deficiency reduction in the subsequent learned subspaces
%%     - 1 indicates that the threshold output subspace's accumulated information
%%       deficiency is <= idt, in this case this neuron considers its threshold
%%       output subspace mature enough to participate in voting procedure
%% soti 						- same as sot except that all references to threshold output subspace
%%        						are replaced by the complimentary inverted threshold output subspace
%%
%% Note that the information deficiency reduction index (IDR) is equal to
%% (delta_in - delta_ot - delta_oti) and hence can be computed once subspace
%% learning is done.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function elt = make_neuron()
elt = struct( ...
   'func','', ...
   'threshold',0, ...
   'features',{{}}, ...
   'delta_in',1,'tci',{{}}, 'clock',{{}}, 'clock_selected',{{}}, ...
   'out',{{}}, ...
   'presynaptic_neurons',{{}}, ...
   'postsynaptic_neurons',{{}}, ...
   'scf_neurons',{{}}, 'scf_neuron_selected',{{}},...
   'delta_o',1,'tco',{{}}, ...
   'delta_ot',1,'tcot',{{}}, ...
   'delta_oti',1,'tcoti',{{}}, ...
   'nsc',{{}}, ...
   'sot',0, ...
   'soti',0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% 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
%% linear, 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);