function out = lif( task, opt)
% lif: a sample m-file for leaky integrate-and-fire neuron model
% Usage:
% >>lif init	
%	initialization
% >>lif run
%	run a simulation
% August 1999 by K. Doya

% Declaration of global variables to be used repeatedly
% any variable not declared as "global" will be cleared after each call
global	ts ns vname		% time span, state dim, and variable names
global	v0 v1 v2		% reset, threshold, and spike potential
global	Inp				% current input
global	T Xt It 		% time, state, and input tables
% If you want to check or change these variables from command line, type, e.g.
% >>global ts

% If you just type
% >>lif
% a demo script below will be run.
if nargin == 0
	lif init
	% the above is interpreted as "lif( 'init')" in matlab's hacky syntax.
	lif run
	return
end
% If a task is specified:
switch( task)
case 'init'				% Initialization
	% default time span and initial state
	[ts,x] = lint([],[],'init');
	ns = size(x,1);		% state space dimension
	vname = lint([],[],'name');	% variable names
	lif( 'figs');		% setup figures
	lif( 'inp', 2);	% default stimulus
	T = 0;				% initialize time, state, and input trace
	Xt = x';
	It = 0;
case 'run'				% Run the simulation
	figure(1);			% real-time display of phase protrait
	% options for ode simulation: plot first and last variable
	odeopt = odeset('Events','on','OutputFcn','odeplot');
	% integration of 'non-stiff' differential equation
	% time steps are chosen automatically
	% Inp could be constant or (t,I) table
	if size(Inp,2)==2, ts=Inp([1,end],1); end
	t0 = T(end);
	t = ts(1);
	x = Xt(end,:)';
	while t < ts(end)
		[ T1, X1, te] = ode23( 'lint', [t,ts(end)], x, odeopt, Inp);
		if ~isempty(te)		% fire!
			T1 = [ T1; T1(end); T1(end)];
			X1 = [ X1; v2; v0];
		end
		t = T1(end);
		x = X1(end,:)';
		% concatenate with previous output
		T = [ T; t0+T1];	% time 
		Xt = [ Xt; X1];		% state
		if size(Inp,2)==2	% input
			It = [ It; interp1( Inp(:,1), Inp(:,2), T1)];	% table
		else
			It = [ It; repmat( Inp(1), size(T1))];		% const
		end
	end
	figure(2); lif( 'wave');
	% figure(1); lif( 'phase');
case 'figs'
	ss = get(0,'ScreenSize');
	dx = ss(3)/2; dy = ss(4)/2;		% figure offset
	sx = dx-24; sy = dy-48;			% figure size
	figure(1); set(1,'Position',[10,dy,sx,sy],'Name','Phase'); clf;
	grid('on'); xlabel( vname(1)); ylabel( vname(end));
	figure(2); set(2,'Position',[10+dx,dy,sx,sy],'Name','Wave'); clf;
case 'phase'		% Plot a phase portrait
	plot( Xt(:,1), Xt(:,end));		% plot x1 and xn
	grid('on');  xlabel( vname(1)); ylabel( vname(end));
case 'wave'			% Plot waveforms
	for i = 1:ns
		subplot(ns+1,1,i);	% plot i-th state variable
		plot( T, Xt(:,i));
		ylabel( vname(i,:));
	end
	subplot(ns+1,1,ns+1);	% plot current input
	plot( T, It(:));
	ylabel( 'I'); xlabel( 't (ms)');
case 'inp'			% Current input
	if ischar(opt), opt = eval(opt); end	% strig to value
	Inp = opt				% maybe constant or table
case 'pulse'			% Simple pulse
	if ischar(opt), opt = eval(opt); end	% strig to value
	t = [ 0; 0.2; 0.201; 0.25; 0.251; 1]*(ts(end)-ts(1)) + ts(1);
	i = [ 0; 0; 1; 1; 0; 0]*opt;
	Inp = [ t, i]
case 'save'			% Save all global variables, e.g.: >>lif save 'data1'
	save( opt);
case 'load'			% Relaod all global variables, including T, X
	load( opt);
otherwise
	error( [' invalid task: ', task]);
end

% end of lif.m