function circuit_demo(s, rulen, norm_method)
% CIRCUIT_DEMO(s, rulen, norm_method)
% 
% Demonstration of basic circuit described in 
% 'Sensory Optimization by Stochastic Tuning'.
%
% Keyboard controls:
%  '1' - location (toggle between methods of input spike combination)
%  '2' - both
%  '3' - frequency
%  'escape' - stop demonstration
%
% Paramters:
%  s - row vector of receptive field sizes
%    (example and default = [1.0, 2.0])
%  rulen - name of rule used for combination of input responses, this
%    parameter determines which aspect of stimulus controls the circuits
%    adaptation
%    (values = 'both' | 'location' | 'frequency')
%  norm_method -
%    - method of weight normalization
%    (values = 'decay' | 'sum')
% 
% Author: Peter Jurica (pjurica@brain.riken.jp, juricap@gmail.com)
% RIKEN Brain Science Institute, Wako-shi, Japan, 2012
%

    PLOT = 1;
    
    if nargout > 0, PLOT = 0; end
    if nargin < 1, s = [1.0, 2.0]; end
    if nargin < 2,
        rulen = 'both';  % 'both' | 'location' | 'frequency'
    end
    if nargin < 3,
        norm_method = 'decay';  % 'decay' | 'sum'
    end
    
    N = length(s);
    tau = 0.2;
    gamma = tau/2;
    
    fs = 1.0./s;
    f0 = -(0.25 + max(fs)); f1 = -f0;
    x0 = -(0.5 + max(s));   x1 = -x0;    % stimulus space boundaries
    
    g  = @(x_, s_) exp(-0.5*(x_./s_).^2);

    col = 'kgr';
    rule = [];
    set_rule(rulen);

    if strcmp(norm_method, 'sum')
        norm_weights = @(w_,w0_) w_./sum(w_);
    elseif strcmp(norm_method, 'decay')
        norm_weights = @(w_,w0_) w_ - tau.*w0_;
    end

    % random initial condition
    ws = rand(1,length(s));
    
    % optional visualization
    if PLOT,
        title('Dynamics of receptive field size.');
        axu = subplot(121);
        ms = mean(s);
        ylim([0, 0.25]);
        plot([ms, ms], ylim(), 'k:'); hold on;
        wm = (ws*s')./sum(ws);
        pm = plot([wm, wm], ylim(), 'r--');
        xlim([1.4, 1.6]);
        xbins = linspace(1.4, 1.6, 51);
        dx = diff(xbins)./2; xx = xbins(1:end-1) + dx;
        hs = zeros(3,50); ns = zeros(3,50);
        b(1) = bar(xx, hs(1,:), 'FaceColor', [69,69,69]./255);
        b(2) = bar(xx, hs(2,:), 'g');
        b(3) = bar(xx, hs(3,:), 'r');
        ylim([0.0, 0.25]);
        set(gca, 'XTick', [1.45, 1.5, 1.55]);
        set(gca, 'YTick', [0.0, 0.1, 0.2]);
        axis('square');
        xlabel('RF size S'); ylabel('p(S_r)');
                
        axw = subplot(122);
        hold on; box on;
        Np = 1000;
        wp = nan(Np, 2, 3);
        for j = 1:3,
            pwp(j) = plot(wp(:,1,j), wp(:,2,j), '.', ...
                          'Color', col(j), 'MarkerSize', 1.0);
        end
        axis([0.2,0.4,0.2,0.4]);
        set(gca, 'XTick', [0.2, 0.3, 0.4]);
        set(gca, 'YTick', [0.2, 0.3, 0.4]);
        axis('square');
        xlabel('w_1'); ylabel('w_2');
        
        %set(gca, 'Units', 'Pixels');
        uicontrol('Style', 'pushbutton', 'String', 'Location', ...
                  'Position', [ 20 20 150 20], ...
                  'Callback', {@set_rule_ev,'location'});
        uicontrol('Style', 'pushbutton', 'String', 'Both', ...
                  'Position', [180 20 150 20], ...
                  'Callback', {@set_rule_ev,'both'});
        uicontrol('Style', 'pushbutton', 'String', 'Frequency', ...
                  'Position', [340 20 150 20], ...
                  'Callback', {@set_rule_ev,'frequency'});
        
        set(gcf, 'KeyPressFcn', @key);
    end
    
    N_stim = 1;   % number of stimuli presented simulatneously
                  % see also first line of the loop
    Nt = 20;      % number of tests
    y = zeros(Nt, length(s));
    c = zeros(size(ws));

    Nbuf = 100;
    wa   = zeros(Nbuf, length(ws));
    
    % simulation loop
    epoch = 1;
    run = true;
    s0 = (ws*s')./sum(ws);
    while run,
        % generate stimuli
        N_stim = randi(5);     % random number of concurrent stimuli
        stim_x = x0 + (x1-x0).*rand(Nt, N_stim);
        stim_f = f0 + (f1-f0).*rand(Nt, N_stim);
        
        o = ones(size(stim_x));
        % detection occured when both space and feature detections are
        % succesful
        for i = 1:length(s),
            loc = g(stim_x, s(i).*o);       % stimulus in space
            ide = g(stim_f, fs(i).*o);      % stimulus in frequency domain
            
            % obtain average firing rate for multiple concurrent stimuli
            loc = mean(loc, 2);
            ide = mean(ide, 2);
            
            % combine results of measurements in the two domain
            y(:,i) = rule(loc, ide);

            % add noise, make sure some spikes are present
            y(:, i) = y(:, i) + 0.1.*randn(size(y(:, i)));
        end
        
        % readout y_r
        u_r = ws*y';
        % adaptive threshold of the output cell
        % here adaptation is simulated by an additive noise
        Theta_r = mean(u_r);
        y_r = u_r > Theta_r + 0.2*randn(size(u_r));

        % compute input-readout coincidence rate
        c = y_r*y./sum(y);
        
        % update rule
        ws_now = ws;
        ws = ws_now + gamma.*c;
        ws = norm_weights(ws, ws_now);
        
        % update histograms
        s1 = (ws*s')./sum(ws);
        if s0 > 1.4 && s0 < 1.6,
            ds = abs(s1 - s0);
            bi = floor(50*(s0-1.4)/0.2)+1;
            hs(wi,bi) = hs(wi,bi) + 100*ds;
            ns(wi,bi) = ns(wi,bi) + 1;
        end
        s0 = s1;
        
        % update visualization
        wp(mod(epoch-1,Np)+1,:,wi) = ws;
        wa(mod(epoch-1,Nbuf)+1,:) = ws;
        if PLOT & mod(epoch, 5) == 1
            wm = (ws*s')./sum(ws);
            set(pm, 'XData', [wm, wm]);
            set(pwp(wi), 'XData', wp(:,1,wi), 'YData', wp(:,2,wi));
            set(b(wi), 'YData', ns(wi,:)./sum(ns(wi,:))); 
            drawnow();
        end

        epoch = epoch + 1;
    end
    
    function set_rule_ev(h,a,rulen)
        set_rule(rulen);
    end
    
    function set_rule(rulen) 
        if strcmp(rulen, 'both')
            rule = @(a, b) a.*b; wi = 2;
        elseif strcmp(rulen, 'location')
            rule = @(a, b) a; wi = 1;
        elseif strcmp(rulen, 'frequency')
            rule = @(a, b) b; wi = 3;
        end
        disp('Current rule:'), disp(rule);
    end

    function key(varargin)
        if strcmp(varargin{2}.Key, 'escape') == 1,
            run = false;
        elseif strcmp(varargin{2}.Key, '1') == 1,
            set_rule('location');
        elseif strcmp(varargin{2}.Key, '2') == 1,
            set_rule('both');
        elseif strcmp(varargin{2}.Key, '3') == 1,
            set_rule('frequency');
        end
    end
end