function [model_output] = perceptual_rivalry_simulation(model_output)
% Function to implement the perceptual rivalry simulations of the prosciptive
% integration model.

% define functions
f =  @(x) x^2 / (x^2 + 1^2);
gauss = @(x, mu, sigma) exp(-0.5 * ((x - mu)./sigma).^2);

% define parameters 
h = .5;           % step size
time_steps = 12^3;% simulation duration
tau_x = 1;        % timescale of inhibition
tau_a = 125;      % timescale of adaptation
g_i = 7;          % gain of inhibition
g_a = 7;          % gain of adaptation
noise_std = .005; % stadard deviation of stochastic vaiability

% define initial bistable response
if nargin == 0 || isempty('model_output')
    % generate bimodal distribution
    mu = [12 24];
    sigma = 2;
    nPool = 37;
    I = gauss(1:nPool,mu(1),sigma)+gauss(1:nPool,mu(2),sigma);
else
    % use response from proscriptive_integration_model.m output
    I = model_output.readout_layer_response(end,:);
    I = I*1/max(I);
    nPool = numel(I);
end
[~,mu] = findpeaks(I);
[lm, rm] = meshgrid(linspace(-pi,pi,nPool));
l_c = max(0,cos(abs(lm+rm)));
l_c = l_c/sum(l_c(1,:));

% generate empty matricies
X   = zeros(time_steps,nPool);
S_X = zeros(time_steps,nPool);
A   = zeros(time_steps,nPool);

% simulate dynamic rivalry
for t = 1:time_steps
    % compute S(X1) and S(X2)
    for n = 1:nPool
        S_X(t,n) = f(X(t,n));
    end
    for n = 1:nPool
        % adaptation
        A(t+1,n) = A(t,n) + h * (-A(t,n) + g_a * S_X(t,n)) / tau_a;
    end
    for n = 1:nPool
        % neural populations
        X(t+1,n) = max(0, X(t,n) + h * (I(n) - (1 + A(t,n)) * S_X(t,n) - g_i * sum(S_X(t,:).*l_c(n,:)) + randn*noise_std) / tau_x);
    end
end
model_output.rivalry_timecourse = X(:,mu);