%run a complete simulation of a double-WTA configuration of a WTA
%
%params are:
%N -> size of WTA
%W, P1, P2 -> weight matrices (see below)
%Tall -> thresholds for C2
%TtnAll -> thresholds for C1
%inputAmp -> ext input amp for init
%symbInpAmp -> ext input amp for symbols
%dt -> integration timeconstant
%till -> # integration steps
%symbsToProcess -> symbols (string) to process
%symbToTNs -> mapping of symbols to transitions
%
%the weight matrices are:
%W -> describes uniform connectivity on each map
%P1 -> describes connectivity between C1->M2
%P2 -> describes connectivity between C2->C1
%
%the maps/WTAs are:
%Mx -> all units on one WTA, including the inhibitory neuron
%C2 -> WTA that holds the states
%C1 -> WTA that contains the TNs
%
%
%urut/sept08 initial
%urut/oct10 modified to allow noise input to C1 map.
%sm/dec14 reduced to two WTA

function [C2,C1,inpHistory, inpHistorySymb] = runDoubleWTA(N, W, W_TN, P1, P2, Tall, TtnAll, inputAmp,symbInpAmp, dt, till, symbsToProcess, symbToTNs, C1_noise, noiseUpdateFrequency,tau)

if nargin<14
    C1_noise=[];
end

nextSymbOn=7000; %when to start processing
nextSymbOff=0;
nextSymbInd=1;
symbOnTime=1930;
symbOffTime=3000;

%pre-allocate
symbsDone=0;
inp=zeros(N,1);
symbInp=zeros(N,1);
C2=zeros(N,till);
C1=zeros(N,till);
inpHistory=zeros(N,till);
inpHistorySymb=zeros(N,till);
temp = 0;
noiseUpdate = [];
for i=2:till
    inpHistory(:,i)=inp;
    inpHistorySymb(:,i)=symbInp;


    %this could be written more succint, but is intentionally split-up into
    %steps to make it easier to understand. The aim is not to make it
    %optimized for runtime speed!
    dC2 = -C2(:,i-1) + W*C2(:,i-1) + P1*C1(:,i-1) + inp - Tall;       
    dC1 = -C1(:,i-1) + W_TN*C1(:,i-1) + P2*C2(:,i-1) + symbInp - TtnAll;
        
    dC2= (dC2)*(dt/tau) ;
    dC1= (dC1)*(dt/tau) ;

    %apply rectification funct
    C2(:,i) = max( (C2(:,i-1)+dC2 ), 0);
    C1(:,i) = max( (C1(:,i-1)+dC1 ), 0);
    
    %external input on/off (initial reset pulse)
    isSet=0;
    if i>1 & i<4000
        inp=inputAmp;
        isSet=1;
    end


    if ~isSet
        inp=zeros(N,1);
    end

    %external inpt (symbols), to C1
    if i == nextSymbOn & ~symbsDone
        symbInp=zeros(N,1);
        symb = symbsToProcess( nextSymbInd );
        
        %activate them
        TNsToActivate = symbToTNs{symb};
        symbInp(TNsToActivate) = symbInpAmp;
        
        nextSymbOff = i+symbOnTime;
        nextSymbInd = nextSymbInd+1;
        nextSymbOn = nextSymbOff + symbOffTime;
        
        if nextSymbInd>length(symbsToProcess)
            symbsDone=1;
        end
    elseif i == nextSymbOff
        symbInp=zeros(N,1);
    end
    
    if i>15000 & ~isempty(C1_noise)
        %add noise to the C1 map            
        if mod(i,noiseUpdateFrequency)==0  %noise update frequency
            symbInp =  C1_noise(:,i);
        end 
            
        
    end
        
    if mod(i,10000)==0
        disp(['simulation step ' num2str(i) ' of ' num2str(till) ]);
    end
end
