% Galvis et al. 2018 (written by Daniel Galvis)
% Class Burster
% This class is used for all Bursting neural ensembles

% These bursters activate to max activation in response to an excitatory
% input, stay active for burst_length regardless of if the input is
% removed, then inactivate and stay inactivated until the refractory period
% has passed. They can then be reactivated by a new stimulus

%PROTOCOL FOR USE!!!
% Bursters must be initialized using X = Burster([intra_delay,extra_delay]...
% ,burst_length,refrac_period,dom_frac); Where burst_length/refrac_period
% are in number of timesteps. Dom_frac describes the weight of extra
% nucleus connections for HVC bursters (usually 0.9 for dom/0.1 for non
% dom)

% At each time step, run X = X.Step([...]) for all Spikers and Bursters, 
% then run X = X.Syn_step([...]) for all Spikers and Bursters.
% See the main code or below for the inputs into each of these functions


classdef Burster
    properties
     %Each of these variables is described in the constructor
        V;
        V_now;
        
        
        S_close;
        S_far;

        
        dom_frac;
        
        
        burst_count;
        refrac_count;
    end
    properties(GetAccess=private)
        %Each of these variables is described in the constructor
        delays;
        burst_length;
        refrac_period;
        on; %used to construct arrays of objects
    end
    
    
    
    methods
        
        %Constructor/Initialization
        function  obj = Burster(init_delays,init_burst_length,...
                                init_refrac_period,init_dom_frac)
            if nargin>0
            %The initializer has 4 inputs
            %-"delays" should have two values: (1) intra-nucleus delay
            %(timesteps) and (2) extra-nucleus delay (timesteps)
            
            %"Burst length" and "Refractory_Period" in timesteps
            
            %Dom_frac >0.5 for dominant hemisphere 
            %dom_frac for dom and non-dom hemispheres should sum to 1
            
            obj.delays = init_delays;
            obj.burst_length = init_burst_length;
            obj.refrac_period = init_refrac_period;
            obj.dom_frac = init_dom_frac;
           
            % V is the value of the Burster, but also tracks previous
            % values. The value of synaptic connections is
            % based on previous values of V because of delay
            
            % At each step entry k of this array is updated with the new
            % value of V. k-1...1,n,n-1,...k+1 is the ordering of past
            % values of V (assuming V is an n sized array).
            obj.V = zeros(max(obj.delays),1);
            
            
            %V_now is the index of the immediate value of V
            
            %obj.V_now = 0; (initialized in Step Function)

            %These are values of the synaptic connections at this timestep
            % S_close is intra/S_far is extra 
            obj.S_close = 0;
            obj.S_far = 0;
            
            
            %These are counters that make sure the Burster stays active
            %for burst_length when activated, then stays silent for at 
            %least refrac_period afterwards
            obj.burst_count = 0;
            obj.refrac_count = obj.refrac_period;
            end     
        end
        
        
        %Step function. Finds the new value of the burster given it's
        %current value and the values of its excitatory and inhibitory
        %inputs
        function    obj = Step(obj,step_times,step_exc,step_inh)
            
            
            %update index of V at this timestep
            obj.V_now = (1+mod(step_times,size(obj.V,1)));
            
            
            %value of V at the previous timestep
            aux_V = obj.V(1+mod(step_times-1,size(obj.V,1)));
            
            
            %The following if statements update the value of V at this
            %timestep
            

            % The only way V = 1 is if 
            %           1) step_exc>0 and burst_count = 0
            %           2) burst_count<burst_period, V = 1 at the
            %           previous timestep
            %           AND!!!!
            %           refractory_count>=refractory_period
            %           step_inh = 0
            
            if(step_inh==0)
                if(obj.refrac_count>=obj.refrac_period && obj.burst_count<...
                        obj.burst_length)
                    if((obj.burst_count ==0&&step_exc>0)||obj.burst_count~=0)  
                        obj.V(obj.V_now) = 1;
                    else
                        obj.V(obj.V_now) = 0;
                    end
                else
                    obj.V(obj.V_now) = 0;
                end
                
            else
                obj.V(obj.V_now) = 0;
            end
            %End updating the value of V at this timestep
            
            
            
            
            % The following updates refractory_count and burst_count
            
            % If V just turned off refrac starts counting up from 1
            % burst count resets to 0 (for the next time it bursts)
            if(aux_V>0 && obj.V(obj.V_now)==0)
                obj.refrac_count = 1;
                obj.burst_count = 0; 
            end
            
            % If V has been off keep adding 1 to the refrac_count
            % burst count stays at 0
            if(aux_V==0&&obj.V(obj.V_now)==0)
               obj.refrac_count = obj.refrac_count + 1; 
               obj.burst_count = 0;
            end
            
            % If V is on keep adding 1 to the burst count
            if(obj.V(obj.V_now)>0)
               obj.burst_count = obj.burst_count + 1;
            end
          %End updating refrac and burst  
     
        end
        
        %This function updates the synaptic weights for both intra and
        %extra nucleus connections. This function should be run AFTER step
        %this function updates these values given the values this timestep
        %for the next timestep
        function    obj = Syn_Step(obj,syn_times)
            
            %moves back the current index by delays(x)-1 units. On the next
            %timestep, this will contain the values of V from the time
            %point delays(x) in the past.
            del_1 = 1+mod(syn_times-obj.delays(1)+1,size(obj.V,1));
            del_2 = 1+mod(syn_times-obj.delays(2)+1,size(obj.V,1));
            
            %Update synaptic weights for intra-nucleus connections
            %Dominance is instantiated here!!!!!
            
               obj.S_close = obj.V(del_1);
               obj.S_far   = obj.dom_frac*obj.V(del_2);
        end
        
        
        
    end
    
    
    
end