%
% This is a short demo which shows how to use 'csim_lifnet.c'
%
% We implement a network with one spike input neuron,
% one analog input neuron, one exzitatory and
% one inhibitory neuron with simple connections between them:
%
%                          syn4
%  analog input(4) ------------------------+
%                                          |
%                   syn1           syn2    v
%  spike input(1) ------> exz(2) ------> inh(3) -+
%                           o                    |
%                           |        syn3        |
%                           +--------------------+
%               
% Author: Thomas Natschlaeger, 15/11/2001
%

clear all

%
% define constants for synapse and neuron types 
%
csim_defs

%
% define the neurons
%

% the input neuron (index 1)
Neurons(1,1) = SPIKE_INPUT;

% a leaky integrate and fire neuron (index 2)
Neurons(1,2) = LIF_NEURON; % type
Neurons(2,2) = 0.015;      % threshold  
Neurons(3,2) = 0.003;      % refractory period
Neurons(4,2) = 0.030;      % tau_m
Neurons(5,2) = 0.005;      % V_reset
Neurons(6,2) = 0.002;      % I_back
Neurons(7,2) = 0.001;      % V_init

% another leaky integrate and fire neuron (index 3)
Neurons(1,3) = LIF_NEURON; % type
Neurons(2,3) = 0.010;      % threshold  
Neurons(3,3) = 0.002;      % refractory period
Neurons(4,3) = 0.050;      % tau_m
Neurons(5,3) = 0.000;      % V_reset
Neurons(6,3) = 0.000;      % I_back
Neurons(7,3) = 0.000;      % V_init

% an analog input wave form unit
Neurons(1,4) = ANALOG_INPUT;

%
% set up the connections
%

% exzitatory synapse 1 (a static one) connects neuron 1 (the input) to neuron 2;
Synapses(1,1) = STATIC_SYN; % type
Synapses(2,1) = 1;          % pre
Synapses(3,1) = 2;          % post
Synapses(4,1) = 0.5;        % weight
Synapses(5,1) = 0.001;      % delay
Synapses(6,1) = 0.003;      % tau_s
  
% exzitatory synapse 2 (a dynamic one) connects neuron 2 to neuron 3;
Synapses(1,2)  = DYNAMIC_SYN; % type
Synapses(2,2)  = 2;           % pre
Synapses(3,2)  = 3;           % post
Synapses(4,2)  = 0.5;         % weight
Synapses(5,2)  = 0.001;       % delay
Synapses(6,2)  = 0.003;       % tau_s
Synapses(7,2)  = 0.00;        % -- not used --
Synapses(8,2)  = 0.2;         % U
Synapses(9,2)  = 1.0;         % D
Synapses(10,2) = 0.3;         % F
  
% inhibitory synapse 3 (a dynamic one) connects neuron 3 back to neuron 2;
Synapses(1,3)  = DYNAMIC_SYN; % type
Synapses(2,3)  = 3;           % pre
Synapses(3,3)  = 2;           % post
Synapses(4,3)  = -0.1;        % weight
Synapses(5,3)  = 0.001;       % delay
Synapses(6,3)  = 0.03;        % tau_s
Synapses(7,3)  = 0.02;        % -- not used --
Synapses(8,3)  = 0.5;         % U
Synapses(9,3)  = 0.5;         % D
Synapses(10,3) = 1.0;         % F
  
% exzitatory synapse 4 (an analog onw) connects neuron 3 back to neuron 2;
Synapses(1,4)  = ANALOG_SYN; % type
Synapses(2,4)  = 4;          % pre
Synapses(3,4)  = 3;          % post
Synapses(4,4)  = 0.018;      % weight
Synapses(5,4)  = 0.0;        % -- not used --
Synapses(6,4)  = 0.0;        % -- not used --
Synapses(7,4)  = 0.01;      % noise
  

%
% We want to plot the Vm's and PSC's of nrns 2 and 3 and of all syns
%
Synapses(1,:)  = bitset(Synapses(1,:),REC_BIT_NMR);
Neurons(1,2:3) = bitset(Neurons(1,2:3),REC_BIT_NMR);

%
% Set up the spiking input: first we have to initialize
% the 'InSpikes' cell array and then place the spikes at
% the right place (neuron 1 emmits a Poisson spike train
% (20 spikes) with a mean rate of 25 Hz:
%
InSpikes{size(Neurons,2)}  = [];
InSpikes{1} = cumsum(exponentialrnd(1/25,1,20));;

%
% Set up the analog input:
%
dt_analog = 1e-3; 
t_analog=0:dt_analog:1;                    % 1.  define the time base
AnalogIn{size(Neurons,2)}  = [];           % 2. init 'AnalogIn'
AnalogIn{4}  = 1+sin(2*pi*6*t_analog);     % 3. place a wave form at right place (neuron 4)

%
% simulation parameters
%
dt      = 1e-4;  % integration time step [sec]
Tsim    = 0.5;   % simulation time [sec]
dt_out  = 0.002; % intervals at which VMs and PSCs are recorded [sec]
dt_disp = 0.1;   % intervals at which to print status to console [sec]
noise   = 0.03;  % the amount of noise [e.g. nA^2]

%
% run the simulation
%
[ST,VM,PSC] = csim_lifnet(Neurons,Synapses,InSpikes,dt,Tsim,...
                          dt_out,dt_disp,noise,AnalogIn,dt_analog);


%
% now we plot the output
%
figure(1); clf reset;

%
% first the spikes
%
subplot(3,1,1); cla reset;
if ~isempty(ST)
  % if there are some spikes we want to plot spikes of neuron 2 in green
  % and spikes from neuron 3 in red so we look for these individual
  % spikes ...
  
  s2=find(ST(1,:)==2);
  s3=find(ST(1,:)==3);

  % and plot them 
  plot(ST(2,s2),ST(1,s2),'b.'); hold on
  plot(ST(2,s3),ST(1,s3),'g.'); hold on
  title('output spikes','fontweight','bold');
else
  title('NO output spikes','fontweight','bold');
end
% also plot the input in black 
% (input spikes are not contained in ST)
plot(InSpikes{1},ones(size(InSpikes{1})),'k.');
legend('neuron 2','neuron 3','input');
xlabel('time [sec]');
ylabel('neuron #');
set(gca,'Xlim',[0 Tsim],'Ylim',[0 4]);

%
% plot membrane voltages
%
subplot(3,1,2); cla reset;

t=0:dt_out:Tsim;                 % these are the times at which VMs are recorded
plot(t,VM(1,:),t,VM(2,:));       % since only neuron 2 and three have the recording
legend('neuron 2','neuron 3');   % flag set, 'VM' only contains this two traces.

title('membrane voltages','fontweight','bold');
xlabel('time [sec]');
ylabel('Vm(t)');
set(gca,'Xlim',[0 Tsim]);

%
% plot PSCs
%
subplot(3,1,3); cla reset;

% plot the PSCs of all four synapses in the net
plot(t,PSC(1,:),'k',t,PSC(2,:),'b',t,PSC(3,:),'g',t,PSC(4,:),'r');
legend('syn 1','syn 2','syn 3','syn 4');
				 
title('PSCs','fontweight','bold');
xlabel('time [sec]');
ylabel('psc(t)');
axis tight
set(gca,'Xlim',[0 Tsim]);