function varargout = nesim_main(p)
%% Nesim 3.8 Version, Jan. 21, 2004
%% Main routine called by mains/runsim.m
%% Copyright (C) by Charles. H. Anderson and Chris Eliasmith (All Rights Reserved)
%% Dept. Anatomy and Neurobiology
%% Washington Univ. School of Medicine
%% St. Louis, MO
%% cha@wustl.edu
%% eliasmith@uwaterloo.ca)
%% Modified Dec. 5, 2003 CHA
% Original Code 10/31/01 by John Harwell
%
%
%%% Information created by setupsim.m .
global Ensemble Connections Connection_InputMode Connection_extID
%%% Time and noise
global T dt time Noise;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% NEW GLOBALS
global ext_INPUTS InputDimension inputIndex %% Input Signals
%%% Storage for the system variables, (1) Inputs to neurons, (2) Outputs, (3) Function Outputs
global sys_INPUTS In_Start In_End sys_OUTPUTS Out_Start Out_End F_OUTPUTS Fout_Start Fout_End
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Ploting information that tells the routine what to save.
global PlotInfo NumPlotInfo
%%%% Storage for single neuron variables
global SpikeTimes SpikeCnt SpikeEnsOffset Activities ActEnsOffset Currents C_EnsOffset Voltages V_EnsOffset
%%% Created by gui_signal_generator
global NumSignalGenerators SignalGenerators SignalGeneratorDimension
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Global Variables needed for saving the output of the run in runsim.m
global sig Xdirect
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global SimName lastSavedName %% Global Variables for saving and restoring a simulation
SimName = p.simName; %%
fullName = [SimName,'.mat'];
if ~exist(fullName)
DataBase = ['..',filesep,'NeuronData',filesep]; % Directory for storing the database.
fullName = [DataBase,SimName,'.mat'];
if( ~exist(fullName))
errstr= sprintf('Cannot find simulation %s\n',SimName);
error(errstr);
end
end
NumEns = 0;
load(fullName);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tauSyn = zeros(NumEns,NumEns);
for i = 1:length(p.ConnData)
a = p.ConnData{i};
Connections(a(1),a(2)).synapse.tc = a(3); %% Modified 11/8/01 by CHA
tauSyn(a(1),a(2)) = a(3);
end
method = p.method;
for(m=1:NumEns)
Ensemble(m).Method = method(m); %%
end
dt = p.dt;
T = p.T;
Noise = p.noiseLevel;
RandSeed = p.randomSeed;
directOnly = p.directOnly;
%%% We get the input dimensions from the ensemble and connection structures.
MaxDimension=0;
InputDimension=zeros(1,NumEns);
Connection_extID = zeros(1,NumEns); %%% New Globals 8/5/03 -cha
Connection_InputMode = zeros(NumEns,NumEns); %%%
for (i=1:NumEns)
Connection_extID(i) = Connections(i,i).extID;
if(Connection_extID(i) >0)
InputDimension(i)=Ensemble(i).D;
if(InputDimension(i)>MaxDimension)
MaxDimension=InputDimension(i);
end
end
for(j=1:NumEns)
Connection_InputMode(i,j) = Connections(i,j).inputMode;
end
end
if(RandSeed>0) randn('state',RandSeed); rand('state',RandSeed); end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% Set up input signals %%%%%%%%%%%%%%%%%
time = [0:dt:T]; Nt = length(time);
totalD = sum(InputDimension);
ext_INPUTS = zeros(Nt,totalD);
inputIndex = zeros(1,NumSignalGenerators);
m=1;
for (i=1:NumSignalGenerators)
inputIndex(i) = m;
Dim = InputDimension(i);
sigInfo = SignalGenerators(i,:);
sig = sim_genSignal(time, sigInfo,Dim);
ext_INPUTS(:,m:m+Dim-1) = sig(1:Dim,:)';
m=m+Dim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% genWeights;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Direct solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
directMethod = zeros(NumEns,1);
RunSimulation(directMethod,tauSyn);
DirectValues = sys_INPUTS;
D_Start = In_Start;
D_End = In_End;
DirectFunction = F_OUTPUTS;
DF_Start = Fout_Start;
DF_End = Fout_End;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Run the simulation
if (directOnly == 0) %% Run with the neurons
RunSimulation(method,tauSyn);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figureNumber = 1;
numEnsemblePlots = 0;
numExternalPlots = 0;
numOutputPlots = 0;
for i=1:NumPlotInfo
if ( PlotInfo(i).type ) % Neurons
if ((PlotInfo(i).inputValues ~= 0) | ...
(PlotInfo(i).directSolution ~= 0))
if (numEnsemblePlots == 0)
ensemblePlotsFigure = figureNumber;
figureNumber = figureNumber + 1;
end
numEnsemblePlots = numEnsemblePlots + 1;
end
if (PlotInfo(i).outputValues ~= 0)
if (numOutputPlots == 0)
outputPlotsFigure = figureNumber;
figureNumber = figureNumber + 1;
end
if(~isempty(Ensemble(PlotInfo(i).ID).fcoefs{1}));
numOutputPlots = numOutputPlots + 2;
else
numOutputPlots = numOutputPlots + 1;
end
end
elseif (PlotInfo(i).type == 0) % External Inputs
if (PlotInfo(i).directSolution ~= 0)
if (numExternalPlots == 0)
externalPlotsFigure = figureNumber;
figureNumber = figureNumber + 1;
end
numExternalPlots = numExternalPlots + 1;
end
end
end
ensembleSubPlotNumber = 1;
externalSubPlotNumber = 1;
ensembleOutputSubPlotNumber = 1;
for i=1:NumPlotInfo
if (PlotInfo(i).type == 0) % external input
if (PlotInfo(i).directSolution ~= 0)
extNum = PlotInfo(i).ID;
dim = SignalGeneratorDimension(extNum);
plotExternalInput(externalPlotsFigure, numExternalPlots, ...
externalSubPlotNumber, ...
time, sig, extNum, dim);
externalSubPlotNumber = externalSubPlotNumber + 1;
end
elseif (PlotInfo(i).type == 1) % neurons
if ((PlotInfo(i).inputValues ~= 0) | ...
(PlotInfo(i).directSolution ~= 0))
ensNum = PlotInfo(i).ID;
title = Ensemble(i).title;
if(length(PlotInfo(i).sel_components))
comps = str2num(PlotInfo(i).sel_components)-1;
else
comps = [0:PlotInfo(i).D-1];
end
inputValues = sys_INPUTS(:,In_Start(i)+comps)';
directValues = DirectValues(:,D_Start(i)+comps)';
plotEnsembleInput(ensemblePlotsFigure, ...
numEnsemblePlots, ensembleSubPlotNumber, time, ...
inputValues, PlotInfo(i).inputValues, ...
directValues, PlotInfo(i).directSolution, ...
title);
ensembleSubPlotNumber = ensembleSubPlotNumber + 1;
end
if ((~directOnly)&(PlotInfo(i).outputValues ~= 0))
Fil_tc = PlotInfo(i).Filter_tc;
outputValues = FilterSignal(sys_OUTPUTS(:,Out_Start(i)+comps)',time,Fil_tc);
directValues = FilterSignal(DirectValues(:,D_Start(i):D_End(i))',time,Fil_tc);
ensNum = PlotInfo(i).ID;
title = Ensemble(i).title;
switch(method(ensNum))
case 0
title = char(title, 'No Neurons');
case 1
title = char(title, 'Activity');
case 2
title = char(title, 'Spikes');
end
title = strjust(title, 'center');
plotEnsembleOutput(outputPlotsFigure, ...
numOutputPlots, ensembleOutputSubPlotNumber, time, ...
outputValues, PlotInfo(i).outputValues, ...
directValues, PlotInfo(i).directSolution, ...
title);
ensembleOutputSubPlotNumber = ensembleOutputSubPlotNumber + 1;
if(~isempty(Ensemble(ensNum).fcoefs{1}));
Fil_tc = Connections(i,i).synapse.tc;
title = ['Func ' Ensemble(i).title];
switch(method(ensNum))
case 0
title = char(title, 'No Neurons');
case 1
title = char(title, 'Activity');
case 2
title = char(title, 'Spikes');
end
title = strjust(title, 'center');
outputValues = FilterSignal(F_OUTPUTS(:,Fout_Start(i):Fout_End(i))',time,Fil_tc)';
directFunction = FilterSignal(DirectFunction(:,DF_Start(i):DF_End(i))',time,Fil_tc)';
plotEnsembleOutput(outputPlotsFigure, ...
numOutputPlots, ensembleOutputSubPlotNumber, time, ...
outputValues, PlotInfo(i).outputValues, ...
directFunction, PlotInfo(i).directSolution, ...
title);
ensembleOutputSubPlotNumber = ensembleOutputSubPlotNumber + 1;
end
end
end
end
if (directOnly == 0)
for i=1:NumPlotInfo
if (PlotInfo(i).type ~= 0)
ensNum = PlotInfo(i).ID;
numNeuron = PlotInfo(i).numNeurons;
if (PlotInfo(i).rasters ~= 0)
thisFilePlotRasters(figureNumber, ensNum, numNeuron);
figureNumber = figureNumber + 1;
end
if (PlotInfo(i).somaCurrents ~= 0)
plotCurrents(figureNumber, time, ensNum, numNeuron);
figureNumber = figureNumber + 1;
end
if (PlotInfo(i).somaVoltages ~= 0)
plotSomaVoltages(figureNumber, time, ensNum, numNeuron);
figureNumber = figureNumber + 1;
end
if (PlotInfo(i).act ~= 0)
plotActivities(figureNumber, time, ensNum, numNeuron);
figureNumber = figureNumber + 1;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot an external input signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = plotExternalInput(figNum, numSubPlots, subPlotIndex, ...
t, sig, extNum, dim)
global ext_INPUTS inputIndex; %%% CHA 1/28/04 Add this line
if (subPlotIndex == 1)
fig = figure(figNum);
clf;
set(fig, 'Name', 'External Inputs');
end
subplot(numSubPlots, 1, subPlotIndex);
for (i=1:dim) %%% CHA 1/21/04
temp = ext_INPUTS(:,inputIndex(extNum)+i-1); %%% CHA 1/28/04 And this line
plot(t, temp);
hold on;
end
if (subPlotIndex == numSubPlots)
xlabel('time (sec)');
end
ylabel(['Input(' num2str(extNum) ')']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plot input vectors to neurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = plotEnsembleInput(figNum, numSubPlots, subPlotIndex, ...
t, decoded, plotDecoded, ...
direct, plotDirect, title)
if (subPlotIndex == 1)
fig = figure(figNum);
clf;
set(fig, 'Name', 'Input System Values');
end
fig = figure(figNum);
subplot(numSubPlots, 1, subPlotIndex);
if ((plotDecoded ~= 0) & (plotDirect ~= 0))
plot(t, decoded');
hold on;
plot(t, direct', 'r:');
elseif (plotDecoded ~= 0)
plot(t, decoded');
elseif (plotDirect ~= 0)
plot(t, direct','r:');
end
if (subPlotIndex == numSubPlots)
xlabel('time (sec)');
end
ylabel(title);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plot output vectors decoded from neurons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = plotEnsembleOutput(figNum, numSubPlots, subPlotIndex, ...
t, decoded, plotDecoded, ...
direct, plotDirect, title)
if (subPlotIndex == 1)
fig = figure(figNum);
clf;
set(fig, 'Name', 'Output System Values');
end
fig = figure(figNum);
subplot(numSubPlots, 1, subPlotIndex);
if ((plotDecoded ~= 0) & (plotDirect ~= 0))
plot(t, decoded');
hold on;
plot(t, direct', 'r:');
elseif (plotDecoded ~= 0)
plot(t, decoded');
elseif (plotDirect ~= 0)
plot(t, direct','r:');
end
axis tight;
if (subPlotIndex == numSubPlots)
xlabel('time (sec)');
end
ylabel(title);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plot rasters produced by an ensemble
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = thisFilePlotRasters(figNum, ensNumber,dn)
global Ensemble SpikeTimes SpikeCnt SpikeEnsOffset T;
if (Ensemble(ensNumber).Method ~= 2)
return;
end
fig = figure(figNum);clf;
set(fig, 'Name', ['Spike Rasters for ' Ensemble(ensNumber).title]);
markers = ['rk'];
%% SpikeCnts = SpikeData(ensNumber).cnt{1};
%% SpikeTimes = SpikeData(ensNumber).times{1};
N = Ensemble(ensNumber).N;
pSpikeCnt = SpikeCnt(SpikeEnsOffset(ensNumber)+1:SpikeEnsOffset(ensNumber)+N);
pSpikeTimes = SpikeTimes(:,SpikeEnsOffset(ensNumber)+1:SpikeEnsOffset(ensNumber)+N);
%NeuronNum = length(SpikeCnts);
for(n=1:dn:N)
subplot(1,1,1);
numspikes = pSpikeCnt(n)-1;
if(numspikes>0)
s_times = pSpikeTimes(:,n);
index = find(s_times>0);
s_times = s_times(index);
%if(pSpikeTimes(numspikes)>T)
% numspikes = numspikes-1; % Don't plot the last spike if it runs past T
%end
plot(s_times,n,[markers(1+mod(1,2)),'.'],'MarkerSize',4);
hold on;
ytxt = sprintf('Ensemble %s', Ensemble(ensNumber).title);
ylabel(ytxt);
end
end
title('Spike Rasters');
xlabel('Time in seconds');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plot soma currents produced by an ensemble
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = plotCurrents(figNum, t, ensNumber, dn)
global Ensemble Currents C_EnsOffset
if (Ensemble(ensNumber).Method < 1)
return;
end
N = Ensemble(ensNumber).N;
pCurrents = (Currents(:,C_EnsOffset(ensNumber)+1:C_EnsOffset(ensNumber)+N))';
fig = figure(figNum);clf;
set(fig, 'Name', ['SomaCurrents for ' Ensemble(ensNumber).title]);
plot(t, pCurrents(1:dn:N,:));
title('Soma Currents');
ylabel('Current');
xlabel('time (sec)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plot soma voltages produced by an ensemble
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = plotSomaVoltages(figNum, t, ensNumber, dn)
global Ensemble Voltages V_EnsOffset
if (Ensemble(ensNumber).Method ~=2)
return;
end
N = Ensemble(ensNumber).N;
pVoltages = (Voltages(:,V_EnsOffset(ensNumber)+1:V_EnsOffset(ensNumber)+N))';
Nt = length(pVoltages(1,:));
displace = [1:dn:N]'*ones(1,Nt)/dn;
fig = figure(figNum);clf;
titleStr = ['SomaVoltage for ' Ensemble(ensNumber).title];
set(fig, 'Name', titleStr);
plot(t,displace+pVoltages(1:dn:N,:));
title(titleStr);
ylabel('Voltage');
xlabel('time (sec)');
axis tight;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plot Activities produced by an ensemble
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = plotActivities(figNum, t, ensNumber, dn)
global Ensemble Activities ActEnsOffset
if (Ensemble(ensNumber).Method == 1)
titleStr = 'Activities';
else
return;
end
N = Ensemble(ensNumber).N;
pActivities = (Activities(:,ActEnsOffset(ensNumber)+1:ActEnsOffset(ensNumber)+N))';
fig = figure(figNum);clf;
set(fig, 'Name', [titleStr ' for ' Ensemble(ensNumber).title]);
plot(t, pActivities(1:dn:N,:));
title(titleStr);
ylabel('Spikes/sec');
xlabel('time (sec)');
function FilSignal = FilterSignal(Signal,t,TimeConstant)
Nt = length(t);
dt = t(2)-t(1);
if(TimeConstant<dt)
eps = 1;
else
eps = dt/TimeConstant;
end
FilSignal = zeros(size(Signal));
FilSignal(:,1) = Signal(:,1);
for(n=2:Nt);
FilSignal(:,n) = (1-eps)*FilSignal(:,n-1)+eps*Signal(:,n);
end
