%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Simulator for Attentional Strategies
%
% By: Kevin Passino
% Version: 4/19/01
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

Tmax=100;  					% Length of simulation time (real time, in seconds)
Tsample=0.01; 				% Sampling period (used to discretize)
Nsteps=(1/Tsample)*Tmax;  	% Number of steps to run the simulation
N=4;         				% Number of predators/prey
T=0*ones(N,Nsteps+1); % T(i,k) will denote the number of steps that predator/prey i has
					  % been ignored as of step k.  This simply allocates memory 
					  % for it.

%-------------------------
% Predator/prey characteristics:
%-------------------------

e=0*ones(N,Nsteps+1); % Allocate memory for predators/prey signals

% Define predators/prey appearance instants (note that f(i,1) refers to predator/prey i, and that predator/prey i
% has priority i, where higher values of i indicate higher priority predators/prey)

f(:,1)=0*ones(N,1); % Allocate memory
f(1,1)=1.2*(1/Tsample);  % Sets that we want an appearance from predator/prey 1 every f(1,1) steps (1.2 sec.)
						 % (this sets the frequency of appearance)
f(2,1)=1.5*(1/Tsample);  
f(3,1)=2*(1/Tsample);  
f(4,1)=1*(1/Tsample);  
%f(5,1)=1.7*(1/Tsample);
fcounter(:,1)=ones(N,1); % Allocate memory for counters 

% Next, define the bounds delta^i on the maximum amount of time between the appearances of each predator/prey
% (note that these must be consistent with the choice for the f(;,1) parameters above that define when the 
% appearances occur).  These are given in real time, not steps.

deltai(1,1)=1.3; % Overbounds a choice of 1.2 above
deltai(2,1)=1.6; 
deltai(3,1)=2.3;
deltai(4,1)=1.1;
%deltai(5,1)=1.75;

%---------------------------
% Organism characteristics:
%---------------------------

% Pick attentional strategy (scheduler)  type

sched=1; % Chooses the policy of choosing the predator/prey that has been ignored for the longest time
%sched=2; % Chooses the policy of choosing the predator/prey that is the highest priority one that has 
      % been ignored longer than the average one
%sched=3; % Chooses the policy of choosing the predators/prey that may be the most difficult to find
          % (of course if the a(:,1) are all the same then this is the same as sched=1)
%sched=4; % Chooses the policy of choosing the predators/prey that may be most difficult to detect
% For sched=4, need weighting factors
w(1,1)=2;
w(2,1)=1;  
w(3,1)=1;  
w(4,1)=1;  
peakest=0*ones(N,1); % Allocate memory for computation of estimate of peak for sched=4

% Set the a_i parameters so that the capacity condition is met:
epsilon=0.1;                    % Used to define a_i parameters next
a(:,1)=epsilon*(1/N)*ones(N,1); % Picks so will meet capacity condition - scale it 
                                      % (other choices are possible)

% Another choice:
%a(1,1)=0.1;
%a(2,1)=0.2;
%a(3,1)=0.3;
%a(4,1)=0.399;

deltas=1*Tsample; % Sets the delay (in steps) that it takes before the organism can 
                 % switch from focusing on one predator/prey to another.
				 % NOTE: This only accounts for part of the "delta" used in the book.
				 % This is the part due to switching from one predator/prey to another.  
deltacounter=0; % Used to implement the deltas delay

decisiontime=1; % Signals that it is a time (if =1) that pick a new predators/prey 
				% to focus on (initially pick one to focus on)
predpreyfocusedon=0*ones(Nsteps+1,1); % Initialize variable that indicates which predators/prey is 
					% currently being focused on (if in a delay period, means
					% the one that it is trying to focus on)
predpreyprocess=0;	% Flag to indicate that switch delay is over
detectedpredprey=0; % Flag that indicates that the predators/prey that is focused on is detected
					% (=0 means not detected, =1 means detected)

%-------------------------------------------------
% Variables for measuring performance:
%-------------------------------------------------

avgT=0*ones(Nsteps+1,1); % avgT(k,1) denotes the average size of T at each step
						 % (average of N values)
maxT=0*ones(Nsteps+1,1); % maxT(k,1) denotes the maximum size of any T at each step
						 % (maximum of N values)
avgP=0; % Time average of priorities

							
%-----------------------------------
% Start simulation loop:
%-----------------------------------


for k=1:Nsteps
	
	for i=1:N	 % First generate the appearance signals at step k
		if fcounter(i,1)>=f(i,1); % Check if have waited long enough to appear
								  %(notice here have >= so since due to how the frequency of
								  % occurences of appearances are defined so that f can be  fraction)
			e(i,k+1)=1; % Appear
			fcounter(i,1)=1; % Reset counter
		else
			e(i,k+1)=0; % Do not appear
			fcounter(i,1)=fcounter(i,1)+1; % Increment counter
		end
	end

	% Compute some performance measures that can be used in decision making:
	% (and for plotting later)
	avgT(k,1)=(1/N)*sum(T(:,k));
	maxT(k,1)=max(T(:,k));

	if decisiontime~=1; 
		predpreyfocusedon(k,1)=predpreyfocusedon(k-1,1); % If not a decision time then same focus as last time
	end
	
	% Decide which predators/prey to focus on:
	if decisiontime==1; % Test if it is a decision time

		% Policy of choosing the predator/prey that has been ignored for the longest time
		if sched==1
		[val,predpreyfocusedon(k,1)]=max(T(:,k)); % Pick predator/prey that has largest T at time k
		end

		% Policy of choosing the predator/prey that is the highest priority one that has 
      	% been ignored longer than the average one		
		if sched==2
		for j=1:N  % Consider predators/prey numbered 1,..,N with N being the highest priority predator/prey
			if T(j,k)>=avgT(k,1)  % Pick the predator/prey encountered that is above the average (so
									% this picks the highest priority one that is above the average)
									% Of course there are some wasted computations here...
				predpreyfocusedon(k,1)=j;
			end
		end
		end
	
		% Policy of choosing the predator/prey that may be the most difficult to find
		if sched==3
		[val,predpreyfocusedon(k,1)]=max(a(:,1).*T(:,k)); % Scales the T(i,k) 
														% elements by a(i,1) in making choice
		end
		
		% Policy of choosing the predators/prey that may be most difficult to detect
		if sched==4 
			for j=1:N
			peakest(j,1)=w(j,1)*(a(j,1)/(1-a(j,1)))*(T(j,k)+deltas+deltai(j,1));
			end
		[val,predpreyfocusedon(k,1)]=max(peakest); % Picks the predator/prey expected to have the highest peak
		end
		
		decisiontime=0; % Do not decide again until detect that predators/prey (consider it
						% "detected" when T(i,.) goes to zero)
						
	end

	if deltacounter>=deltas  % First, implement delay due to switching from one predator/prey to another
		predpreyprocess=1;  % Flag to indicate that we are now going to focus on the predator/prey

	else
		deltacounter=deltacounter+Tsample; % Increment delay since computing T(i,k+1) (note that
										% at every time we are focusing on something so
										% at every time we need to compute this delay)
										% Also, note that below the test for the delay 
										% accounts for incrementing the counter past delta

	end

	if predpreyprocess==1  % If switch delay is over, then can start trying to detect predator/prey
		if e(predpreyfocusedon(k,1),k)==1 % Consider this test for "detection" to only indicate one pulse
			detectedpredprey=1;  % If detect predators/prey set flag
%			D=k; % The time that an predators/prey pulse was found
		end
	end

	if predpreyprocess==1 & detectedpredprey==1
		deltaT=((1-a(predpreyfocusedon(k,1),1))/a(predpreyfocusedon(k,1),1))*Tsample;
		if T(predpreyfocusedon(k,1),k)-deltaT>0;  % Must test if will hit zero (to keep nonneg)
			T(predpreyfocusedon(k,1),k+1)=T(predpreyfocusedon(k,1),k)-deltaT; 
		else 
			T(predpreyfocusedon(k,1),k+1)=0; % Set it to zero since it finishes (an approximation
							% that is accurate in the sense that it is conservative)
			deltacounter=0; % Reset the delay counter to prepare for next predator/prey
			decisiontime=1; % Signal a decision time since just detected an predator/prey
			detectedpredprey=0; % Reset for use in detection of next predator/prey
			predpreyprocess=0;  % Reset for next loop
		end
		
	else
		T(predpreyfocusedon(k,1),k+1)=T(predpreyfocusedon(k,1),k)+Tsample; % Increment amount of time 
																		% since last detect
	end

	for i=1:N 
		if i~=predpreyfocusedon(k,1); % If predator/prey that is focused on is not i
			T(i,k+1)=T(i,k)+Tsample; % Update the amounts of time since last focused on each predator/prey
		end
	end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

time=0:Tsample:Tmax;

figure(1)
clf
subplot(411)
stem(time,e(1,:))
title('Predator/prey appearances, N=4')
subplot(412)
stem(time,e(2,:))
subplot(413)
stem(time,e(3,:))
subplot(414)
stem(time,e(4,:))
xlabel('Time, sec.')

figure(2)
clf
subplot(311)
stairs(time,predpreyfocusedon)
title('Predator/prey focused on and length of time since last detect (i=1,2)')
subplot(312)
plot(time,T(1,:))
subplot(313)
plot(time,T(2,:))
xlabel('Time, sec.')

figure(3)
clf
subplot(311)
stairs(time,predpreyfocusedon)
title('Predator/prey focused on and length of time since last detect (i=3,4)')
subplot(312)
plot(time,T(3,:))
subplot(313)
plot(time,T(4,:))
xlabel('Time, sec.')

% Compute the time average of the priorities (like a measure of control energy).
% This gives a measure of how well the scheduling policy did on focusing on higher
% priority pedators/prey 

disp('----------------------------------------------------------------')
sched
disp('----------------------------------------------------------------')
disp('The time average of the priorities is=')
avgP=(1/(Nsteps+1))*sum(predpreyfocusedon)

figure(4)
clf
subplot(211)
stairs(time,predpreyfocusedon)
title('Predator/prey focused on')
xlabel('Time, sec.')
hold on
plot(time,avgP,'r.')
hold off
subplot(212)
plot(time,T)
xlabel('Time, sec.')
title('Lengths of times since last detected each predator/prey')

% Compute some performance measures:

disp('The time average of the average values of the lengths of times waited=')
(1/(Nsteps+1))*sum(avgT)  % Computes the time average of the average values

disp('The maximum of the average values of the lengths of times waited=')
max(avgT)  % Computes the maximum of the average values

disp('The time average of the maximum values of the lengths of times waited=')
(1/(Nsteps+1))*sum(maxT)  % Computes the time average of the max values

disp('The maximum of the maximum values of the lengths of times waited=')
max(maxT)  % Computes the maximum of the max values
disp('----------------------------------------------------------------')


figure(5)
clf
subplot(211)
plot(time,avgT,time,(1/(Nsteps+1))*sum(avgT),'r.')
title('Average length of time since last detect')
xlabel('Time, sec.')
subplot(212)
plot(time,maxT,time,(1/(Nsteps+1))*sum(maxT),'r.')
xlabel('Time, sec.')
title('Maximum length of time since last detect')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
