function closed_loop

% 1/28/2012: This function computes the solution to a simple two-unit
% system of sigmoidal units with input to the first, excitation of the
% second by the first, and return excitation of the first by the second (of
% varying strength).
%
% It simulates varying threshold levels applied to evidence accumulated
% under a fixed signal-to-noise ratio.
%
% This simulation is deterministic. A stochastic Euler method should be
% used with the system in the function 'local_closed_loop_ode' if noise is
% not 0.


num_spans = 10;


options = odeset('Events',@events_function);



% Network parameters:
lambda = 4;

% accum_bias = -0.5;
accum_offset = -0.6;
switch_offset = -1.4;

integrator_self_excitation = 1;

latch_self_excitation = 2;

thresh_to_accum_wt = 0;
thresh_to_accum_wt = 0.5;

shutoff_to_accum_wt = -1;
accum_to_thresh_wt = 1.5;
shutoff_to_thresh_wt = -1;
shutoff_to_thresh_wt = -0.5;
thresh_to_shutoff_wt = 1;

% threshold_bias_signal = 0;

threshold_reductions = [-0.555 -0.55 -0.54 -0.5 -0.4 -0.1];
% threshold_reductions = -0.4;
color_list = [ 1 0 0; 0 1 0; 0 0 1 ];
color_list = [ 1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1; 0 1 1; 0 0 0 ];


figure(1);
hold on

figure(2);
hold on


for j = 1:length(threshold_reductions)
    
    T = [];
    Y = [];
    
    t_current = 0;
    
    Y0 = [0 0 0 0];

    threshold_bias_signal = threshold_reductions(j);

    for i = 1:num_spans
        
        [t,y, TE,YE,IE] = ode45(@local_closed_loop_ode,[t_current t_current+300], Y0, options);
        
        
        T = [T; t];
        Y = [Y; y];
        
        Y0 = y(end,:);

        
        % If there's an event, find its time index to help with
        % response-locked plotting:
        if ~isempty(TE)
            idx = find(T>=TE,1);
        end
        
        
        % During one window, set an input signal to some level, else shut it
        % off:
        if i == 2
            Y0(1) = 0.11;
        else
            Y0(1) = 0;
        end
        
        t_current = T(end);
        
    end
        
    figure(1)
    %     plot(T,Y(:,3),'Color',color_list(j,:));
    plot(T,Y(:,3),'Color','k');
    hold on
    
    figure(2)
    %     plot( T(idx-150:idx+50) - T(idx), Y(idx-150:idx+50,3), 'Color', color_list(j,:) );
    plot( T(idx-150:idx+50) - T(idx), Y(idx-150:idx+50,3), 'Color', 'k' );
    hold on
    
end

% keyboard

figure(1)
title('Computed by closed\_loop.m')
% xlim([75 150])

text(80, 0.8, sprintf('accumulator to threshold wt = %f',accum_to_thresh_wt) )
text(80, 0.6, sprintf('threshold to accumulator wt = %f',thresh_to_accum_wt) )





% --------------------------------------------------------------
% The ODE function:
    function dY = local_closed_loop_ode(t,Y)
        
        % This function defines the ode's for a two-unit system forming a closed
        % loop of connections, with one exciting the other, which returns
        % inhibition.      
        
        dY(1) = 0;  % 1/28/12: An input that can be set or unset
        
        dY(2) = -Y(2) + 1./(1+exp(-lambda*(Y(1) + integrator_self_excitation*Y(2) + thresh_to_accum_wt*Y(3) + shutoff_to_accum_wt*Y(4) + accum_offset)));     % 1/28/2012: Creating a simple integrator
        
        dY(3) = -Y(3) + 1./(1+exp(-lambda*(latch_self_excitation*Y(3) + accum_to_thresh_wt*Y(2) + shutoff_to_thresh_wt*Y(4) + switch_offset + threshold_bias_signal)));     % A threshold latch
        
        dY(4) = -Y(4) + 1./(1+exp(-lambda*(latch_self_excitation*Y(4) + thresh_to_shutoff_wt*Y(3) - 1.4)));    % A shutoff switch
        
        dY = dY';
        
%         threshold_bias_signal
        
    end


% --------------------------------------------------------------
    function [value,isterminal,direction] = events_function(t,y)
        % Locate the time when height passes through zero in a
        % decreasing direction and stop integration.
        value = 0.8 - y(4);     % Detect height = 0.8
        isterminal = 1;   % Stop the integration
        % direction = -1;   % Negative direction only
        direction = -1;   % Positive direction only
    end


end