clear all

addpath 'Skuld_2_2' 'Urd_v1_2' 'VerdandiLite_v1_0'

global Network

disp('Preparing simulations')
tic 

% --- Time integration ---
Nsample = 20; %reduces nr of output values, don't change, is hardcoded in C++ code
Network.dt = 0.01; % 0.01 Integration step (in ms)
simdur     = .5e3; % ms of simulation time
Network.nStep = round(simdur/Network.dt); % Number of time steps
tsim = (1:Network.nStep)*Network.dt;

Network.OutputVm       = 0; % 1=Output membrane potentials to Vm.txt, 0=no such output
Network.OutputSp       = 1; % 1=Output spike times to Spikes.txt, 0=no such output
Network.OutputSynState = 1; % 1=Output synapse state to SynState.txt, 0=no such output

fInit = fopen('Inits.txt','w');
fprintf(fInit,'%u\t%g\n%u\t%u\t%u',Network.nStep,Network.dt,Network.OutputVm,Network.OutputSp,Network.OutputSynState);
fclose(fInit);


std_cellparam(1) = 0.005;
std_cellparam(2) = 0.005;
p.offsetg = 0;          % add additional constant synaptic input to the neurons


%% --- Populations ---
   Ncellex = 100;
   Ncellin = 99;
   Lambda = 600;
   Iinput = 0;
   gsynex = randn(1,Ncellex)*std_cellparam(1);
   gsynex = gsynex - mean(gsynex) + p.offsetg;
   gsynin = randn(1,Ncellin)*std_cellparam(1);
   gsynin = gsynin - mean(gsynin) + p.offsetg;
   Ek = -95;

%%
AddPopulation(...
    'Name','Ex',...
    'nCell',Ncellex,...
    'Type','HH',... 
    'Iinput',Iinput ,'Cm',10.0,'g_na',100.0,'g_naL',0.0175,'g_k',40.0,'g_kL',0.05,'g_clL',0.05,'phi',3.0,'E_k',Ek,'E_na',53,'E_cl',-82,...
    'gsyntest', gsynex, 'Esyntest', 50,...
    'V',linspace(-65,-65,Ncellex),'n',0.07,'h',0.97,...
    'Position',repmat(linspace(0,10,Ncellex)',1,3)...
);

AddPopulation(...
    'Name','In',...
    'nCell',Ncellin,...
    'Type','HH',... 
    'Iinput',Iinput ,'Cm',10.0,'g_na',100.0,'g_naL',0.0175,'g_k',40.0,'g_kL',0.05,'g_clL',0.05,'phi',3.0,'E_k',Ek,'E_na',53,'E_cl',-82,...
    'gsyntest', gsynin, 'Esyntest', 50,...
    'V',linspace(-65,-65,Ncellin),'n',0.07,'h',0.97,...
    'Position',repmat(linspace(0,10,Ncellin)',1,3)...
);

AddPopulation(...
    'Name','Input_ex',...
    'nCell', Ncellex, ...
    'Type','PoissonStep',... 
    'Lambda',Lambda,'Lambda2',Lambda/1.5,'StepTime',0.50,...
    'Position',rand(Ncellex,3));

%--- Connections --- 
gammaEx = 1/45;
gammaIn = 1/34;


% give each cell a unique Poisson spike train input
AddConnectivity(...
    'From','Input_ex',...
    'To','Ex',...
    'Type','AlphaLiley',...
    'gamma',gammaEx,'E',50,'g',4*0.5*0.0027,... 
    'Connections','Matrix', eye(Ncellex)...
    ,'Delay','Function',@(r) 0*r);

P = 0.5; % fraction of cells a cell has a connection with
ConMat2 = double(rand(Ncellex,Ncellex)<P);

AddConnectivity(...
    'From','Ex',...
    'To','Ex',...
    'Type','AlphaLiley',...
    'gamma',gammaEx,'E',50,'g',0.5*0.0057/Ncellex/P,...
    'Connections','Matrix', ConMat2...
    ,'Delay','Function',@(r) 0*r);  

P = 0.5; % fraction of cells a cell has a connection with
ConMat3 = double(rand(Ncellin,Ncellex)<P);

AddConnectivity(...
    'From','Ex',...
    'To','In',...
    'Type','AlphaLiley',...
    'gamma',gammaIn,'E',50,'g',0.5*1.32*6*0.0057/Ncellex/P,... %0.015
    'Connections','Matrix', ConMat3...
    ,'Delay','Function',@(r) 0*r);

P = 0.5; % fraction of cells a cell has a connection with
ConMat4 = double(rand(Ncellex,Ncellin)<P);

AddConnectivity(...
    'From','In',...
    'To','Ex',...
    'Type','AlphaLiley',...
    'gamma',gammaEx,'E',-82,'g',15*5*0.5/1.32*0.015/Ncellin/P,... %0.045
    'Connections','Matrix', ConMat4...
    ,'Delay','Function',@(r) 0*r);  

P = 0.5; % fraction of cells a cell has a connection with
ConMat5 = double(rand(Ncellin,Ncellin)<P);

AddConnectivity(...
    'From','In',...
    'To','In',...
    'Type','AlphaLiley',...
    'gamma',gammaIn,'E',-82,'g',0.5*0.015/Ncellin/P,...
    'Connections','Matrix', ConMat5...
    ,'Delay','Function',@(r) 0*r);

%% --- Output to files 
WriteCells;
WriteConnectivity;
save('Network.mat')
toc
%% --- run simulation in C++ environment
disp('Running simulations')
tic 

if(ismac)
    system('VerdandiLite_v1_0/VerdandiLite_OSX');
elseif(ispc)
    system('VerdandiLite_v1_0\VerdandiLite_CPP\Verdandi_Zandt2014.exe');
else
    error('No binaries available for VerdandiLite. Please compile C++ source files and add modify these Matlab lines.')
end

toc
%% generate data for use in Skuld
SkuldOut

%% convert spike timings to seconds for use in Skuld
Spikes = load('Spikes.txt');
Spikes(:,1) = Spikes(:,1)/1e3;
fid = fopen('Spikes.txt','w');
fprintf(fid,'%6.6f  %i\n',Spikes');
fclose(fid);
%%
run_meanfieldapprox_makefigures   % 
break
% to further analyze the spiking dynamics:
Skuld('SkuldSet.mat')