function [STMtx,T,V,dV,Ninp]=IDNetSim(SimPar,ConMtx)
% MATLAB wrapper for C program IDNet.c to simulate abritrary biological
% neural networks. The neuron and synapse models are specified in IDNet.c,
% the parameters SimPar in ConfigIDNet.
%
% INPUT:
%   SimPar: A set of parameters (neurons, synapses, control parameters)
%   ConMtx (optionally): A network structure (random connectivity specified 
%                        in SimPar is used if this is empty).
% 
% OUTPUT:
%   STMtx: Cell array of spike times t_sp. Each cell contains t_sp of one
%          neuron
%   T:     Vector of simulation time steps for visualizing voltage traces
%   V:     Cell array of voltage traces for each neuron in the viewlist
%   Ninp:  Number of simulated neurons



% ------------------ Read parameters from SimPar --------------------------

CtrPar = SimPar.CtrPar;                             % control parameters

Nstripes = SimPar.Nstripes;                         % number of stripes
NTypes = SimPar.NTypes;                             % number of neurons of each type
InpNTypes = SimPar.InpNTypes;                       % number of input neurons of each type

CellPar = SimPar.CellPar;                           % (mean) neuron parameters
k_trans = SimPar.k_trans;                           % parameter for inv_transform_distribution
ParCov = SimPar.ParCov;                             % Cell-Array with 14 covariation matrices
N_sig = SimPar.N_sig;                               % SDs of neuron parameters
N_max = SimPar.N_max;                               % maxima of neuron parameters
N_min = SimPar.N_min;                               % minima of neuron parameters

NPList = SimPar.NPList;                             % type of each neuron

STypPar = SimPar.STypPar;                           % synapse parameters

ConPar= SimPar.ConPar;                              % (mean) connection parameters (weights etc.)
pCon = SimPar.pCon;                                 % connectivity (for random networks)
cluster_flag = SimPar.cluster_flag;                 % determines whether to use a common neighbour rule for connectivity
S_sig = SimPar.S_sig;                               % SDs of connection parameters
S_max = SimPar.S_max;                               % maxima of connection parameters
S_min = SimPar.S_min;                               % minima of connection parameters
ConParStripes = SimPar.ConParStripes;               % parameters for inter-stripe connectivity

ViewList = SimPar.ViewList;                         % neurons to be recorded in detail

InpSTtrains = SimPar.InpSTtrains;                   % input spike trains
NoisePar = SimPar.NoisePar;                         % parameters for random background activity (noise)
NoiseDistr = SimPar.NoiseDistr;                     % Distribution of random background activity
V0Par = SimPar.V0Par;                               % parameters for initial conditions

CA = SimPar.CA;                                     % cell assemblies

NeuronGroupsSaveArray = SimPar.NeuronGroupsSaveArray;   % Array specifying neuron groups to save synapse output for

UniqueNum = SimPar.UniqueNum;
% "UniqueNum" is an integer added to IDNet's temporary data files. This
% allows differentiation of files if multiple instances of IDNet are
% desired.


% -------------------- Set overall parameters -----------------------------
% Seed random number generator to a given state
% mtstream = RandStream.getGlobalStream;
% mtstream.State = SimPar.RState;
% RandStream.setGlobalStream(mtstream);

% Set network parameters
N = sum(NTypes);                % # neurons (per stripe)
NTypesN = length(NTypes);       % # neuron types
M = sum(InpNTypes);             % # input neurons (per stripe)
InpNTypesN = length(InpNTypes); % # input neuron types

% Initalize neuron and synapse parameters
NeuPar = zeros(length(CellPar(:,1)),Nstripes*N);
V0 = zeros(size(V0Par,2),Nstripes*N);
SPMtx=zeros(Nstripes*N,Nstripes*N, max(max(ConPar{1})));

NN = 0;
SynPar=[];
Nsyn=0;
dtax_back = [];

ind1 = 1:9;
ind3 = [1 2 3 4 5 6 8 9 10];


% -------------- Set neuron parameters and most synaptic connections -----------------

% loop over stripes/columns
reorder = cell(NTypesN,Nstripes);
for ii=1:Nstripes
    
    % loop over neuron types
    t_lat = zeros(NTypesN,1);
    t_lat_LIF = zeros(NTypesN,1);
    for i=1:NTypesN
        
        % Set neuron parameters and initial conditions (one individual set for each neuron)
        iset = (1:NTypes(i))+NN;
        reorder{i,ii} = iset;
        ind_out = iset;
        iii=0;
        while(~isempty(ind_out) && iii<1000)           % comment out this whole block before using update_inv_con_PSP

            % draw random numbers from a multivariate normal distribution
            NeuPar_multi=mvnrnd( zeros(1,length(ParCov{i}(:,1))),ParCov{i},length(ind_out))';

            % invert transformation of marginal distributions
            for j=ind1
                NeuPar_multi(j,:) = inv_transform_distribution2(NeuPar_multi(j,:),k_trans(j,i),CellPar(j+1,i),N_sig(j+1,i),N_min(ind3(j)+1,i));
            end
            
            % complete parameter matrix by "C" and "a"
            for j=ind1
                NeuPar(ind3(j),ind_out) = NeuPar_multi(j,:);
            end
            NeuPar(1,ind_out) = NeuPar(1,ind_out).*NeuPar(2,ind_out);                % C from gL and tau
            NeuPar(7,ind_out) = zeros(size(NeuPar(4,ind_out)));                      % a=0
            
            % detect parameters outside boundaries and draw them from a mulitvariate uniform distribution
            ind_out_para = [];
            for j=ind1
                ind_out_para = [ind_out_para find(NeuPar(ind3(j),iset)<N_min(ind3(j)+1,i) | NeuPar(ind3(j),iset) >N_max(ind3(j)+1,i))];  % check parameter boundaries
            end
            ind_out_V2 = find(NeuPar(9,iset)>=NeuPar(10,iset));    % Vr must be smaller than Vth
            ind_out_tau = find(NeuPar(1,iset)./NeuPar(2,iset)<N_min(12,i) | NeuPar(1,iset)./NeuPar(2,iset)>N_max(12,i)); % check tau boundaries
            ind_out_tcw = find(NeuPar(6,iset)<=NeuPar(1,iset)./NeuPar(2,iset)); % tcw must be larger than tau
            ind_out=iset(unique([ind_out_V2 ind_out_tau ind_out_tcw ind_out_para])); 
            iii=iii+1;
        end
        
        iii=0;
        while (~isempty(ind_out) && iii<1000)
            if iii==0
                disp('Number of trials for the full multivariate distribution has been exceeded. Use multivariate uniform distribution for the rest.');
            end
            % draw random numbers from a multivariate uniform distribution
            NeuPar_multi=multivar_distr(length(ind_out),N_min(ind3+1,i), N_max(ind3+1,i), ParCov{i},1)';
            
            % complete parameter matrix by "C", "a" and "Vup"
            for j=ind1
                NeuPar(ind3(j),ind_out) = NeuPar_multi(j,:);
            end
            NeuPar(1,ind_out) = NeuPar(1,ind_out).*NeuPar(2,ind_out);                % C from gL and tau
            NeuPar(7,ind_out) = zeros(size(NeuPar(4,ind_out)));                      % a=0
            
            % detect parameters outside boundaries and draw them from a mulitvariate uniform distribution
            ind_out_para = [];
            for j=ind1
                ind_out_para = [ind_out_para find(NeuPar(ind3(j),iset)<N_min(ind3(j)+1,i) | NeuPar(ind3(j),iset) >N_max(ind3(j)+1,i))];  % check parameter boundaries
            end
            ind_out_V2 = find(NeuPar(9,iset)>=NeuPar(10,iset));    % Vr must be smaller than Vth
            ind_out_tau = find(NeuPar(1,iset)./NeuPar(2,iset)<N_min(12,i) | NeuPar(1,iset)./NeuPar(2,iset)>N_max(12,i)); % check tau boundaries
            ind_out_tcw = find(NeuPar(6,iset)<=NeuPar(1,iset)./NeuPar(2,iset)); % tcw must be larger than tau
            ind_out=iset(unique([ind_out_V2 ind_out_tau ind_out_tcw ind_out_para])); 
            iii=iii+1;
        end
        
        if (~isempty(ind_out))
            error('Error: Number of trials for univariate uniform distribution has been exceeded.');
        end

%         for j = 2:length(SimPar.CellPar(:,1))      % use this block instead of the previous one before using update_inv_con_PSP
%             NeuPar(j-1,ind_out) = rand_par(length(ind_out), CellPar(j,i), N_sig(j,i), N_min(j,i), N_max(j,i), 0);
%         end;
%         NeuPar(3,[2:3 5:11 13:17])  = -60;                     % uniform EL above GABA reversal potential
%         NeuPar(5,[2:3 5:11 13:17])  = 10001; % no spiking in simulated neurons (except for the first one)
%         NeuPar(10,[2:3 5:11 13:17]) = 10000;

                   
        % Compute measures for electrophysiological classification of cells
        % and parameters for depolarization block 
        t_lat_act = zeros(size(iset));
        t_lat_LIF_act = zeros(size(iset));
        adapt_act = zeros(size(iset));
        for j=1:length(iset)

            % meaningful model parameter names
            [C,gL,EL,sf,~,~,~,~,Vr,Vth]=names(NeuPar(:,iset(j)));
            
            % compute parameters for depolarization block:
            % I_ref: Input current beyond which the neuron does not react
            %        any more
            % v_dep: Membrane potential that is used to compute the (input-
            %        independent) differential equation for the membrane
            %        potential during the depolarization block
            IRheo=gL*(Vth-EL)+gL*sf;
            OPTIONS=optimset('TolFun',1e-3,'TolX',1e-3,'Display','off');
            I_ref = fminsearch(@Define_I_ref,IRheo+100,OPTIONS,NeuPar(:,iset(j)));
            NeuPar(11,iset(j)) = I_ref;     % I_ref
            NeuPar(12,iset(j)) = Vr;        % v_dep
            
            % compute latency of first spike
            I   = 500;
            t_lat_act(j) = integral(@(V) C./(I - gL*(V-EL) + gL*sf*exp((V-Vth)/sf) ), EL, Vth/2);
            t_lat_LIF_act(j) = C*log(I/(I+gL*(EL-Vth/2)))/gL;
%             t_lat_act(j) = 1;       % replace the two above lines by these ones before running update_inv_con_PSP
%             t_lat_LIF_act(j) = 1;
            
            
            % compute adaptation ratio
            I = 25:25:300;
            f1 = zeros(length(I),1);
            f  = zeros(length(I),1);
            for k=1:length(I)
                f1(k) = FRsimpAdEx(NeuPar(:,iset(j)),I(k),0,NeuPar(3,iset(j)),[]);
                f(k)  = FRsimpAdEx(NeuPar(:,iset(j)),I(k),[],[],[]); 
            end;
            if isnan(median(f1(f>0) ./ f(f>0)))
                adapt_act(j) = 0;
            else
                adapt_act(j) = median(f1(f>0) ./ f(f>0));
            end;
            
        end;
        t_lat(iset) = t_lat_act;
        t_lat_LIF(iset) = t_lat_LIF_act;
        adapt(iset) = adapt_act;
        
        % finally, set initial conditions
        V0(2,iset)=zeros(1,length(iset));
        for j=1:length(iset)
            V0(1,iset(j))=NeuPar(3,iset(j));   
        end
        
        NN = NN+NTypes(i);
    end;
    
    
    % Redistribute neuron types: I-L -> I-L-d             % comment out the redistribution before using update_inv_con_PSP
    ind_L23 = [reorder{2,ii} reorder{3,ii}];
    ind_L23_Ld = ind_L23((t_lat(ind_L23) - t_lat_LIF(ind_L23))>0);
    ind_L23_L = ind_L23(~ismember(ind_L23, ind_L23_Ld));
    NTypes(2) = length(ind_L23_L);
    NTypes(3) = length(ind_L23_Ld);
    reorder{2,ii} = ind_L23_L;
    reorder{3,ii} = ind_L23_Ld;
    
    ind_L5  = [reorder{9,ii} reorder{10,ii}];
    ind_L5_Ld = ind_L5((t_lat(ind_L5) - t_lat_LIF(ind_L5))>0);
    ind_L5_L = ind_L5(~ismember(ind_L5, ind_L5_Ld));
    NTypes(9) = length(ind_L5_L);
    NTypes(10) = length(ind_L5_Ld);
    reorder{9,ii} = ind_L5_L;
    reorder{10,ii} = ind_L5_Ld;
    
    
    % Redistribute neuron types: I-CL -> I-CL-AC
    ind_L23 = [reorder{4,ii} reorder{5,ii}];
    ind_L23_Ld = ind_L23(adapt(ind_L23)>1.5834);
    ind_L23_L = ind_L23(~ismember(ind_L23, ind_L23_Ld));
    NTypes(4) = length(ind_L23_L);
    NTypes(5) = length(ind_L23_Ld);
    reorder{4,ii} = ind_L23_L;
    reorder{5,ii} = ind_L23_Ld;
    
    ind_L5  = [reorder{11,ii} reorder{12,ii}];
    ind_L5_Ld = ind_L5(adapt(ind_L5)>1.5834);
    ind_L5_L = ind_L5(~ismember(ind_L5, ind_L5_Ld));
    NTypes(11) = length(ind_L5_L);
    NTypes(12) = length(ind_L5_Ld);
    reorder{11,ii} = ind_L5_L;
    reorder{12,ii} = ind_L5_Ld;

    

    % Set network connectivity
    if isempty(ConMtx)              % connect randomly if ConMtx is empty

        % intra-celltype connections
        for i=1:NTypesN
            if cluster_flag(i,i) == 1                                                                       % set random connectivity submatrix for the current input->output pair
                [X,idc] = SetCon_CommonNeighbour_Recur(Nsyn, length(reorder{i,ii}), length(reorder{i,ii}), pCon(i,i), 0.47, 1);  % with...
            else
                [X,idc] = SetCon(Nsyn, length(reorder{i,ii}), length(reorder{i,ii}), pCon(i,i));            % ... or without common neighbour rule
            end;
            
            for k = 1:length(ConPar{2}{i,i}{1})                                                             % loop over all synapse types
                l = ConPar{2}{i,i}{1}(k);                                                                   % determine which parameter set to use for each type
                ST = ConPar{2}{i,i}{2}(k);                                                                  % get synapse type ID (for SynPar)
                SPMtx(reorder{i,ii}, reorder{i,ii},k)=(X)+(k-1)*length(idc);                                % fit submatrix in
                
                idc_back = idc;
                [SynPar,Nsyn,idc]=SetSyn(SynPar,Nsyn,idc,ST,ConPar{2}{i,i}{l+2}, ...                        % set synapse parameters
                    ConPar{3}{i,i}, S_sig(i,i,:), S_max(i,i,:), S_min(i,i,:));
                if k==1
                    dtax_back = SynPar(6,idc_back);
                else
                    SynPar(6,idc_back) = dtax_back;
                    dtax_back = [];
                end;
            end;
        end;
        
        % inter-celltype connections
        for i=1:NTypesN                              % loop over output neurons
            for j=setdiff(1:NTypesN,i)               % loop over input neurons (except output neuron)
                if cluster_flag(i,j) == 2                                                                       % set random connectivity submatrix for the current input->output pair
                    [X,idc] = SetCon_CommonNeighbour_Cross(Nsyn, SPMtx(length(reorder{i,ii}),length(reorder{i,ii})), SPMtx(length(reorder{j,ii}),length(reorder{j,ii})), pCon(i,j), [], [], 1);  % with...
                    X=X';
                else
                    [X,idc] = SetCon(Nsyn, length(reorder{i,ii}), length(reorder{j,ii}), pCon(i,j));                                                                             % ... or without common neighbour rule
                end;
                
                for k = 1:length(ConPar{2}{i,j}{1})                                                             % loop over all synapse types
                    l = ConPar{2}{i,j}{1}(k);                                                                   % determine which parameter set to use for each type
                    ST = ConPar{2}{i,j}{2}(k);                                                                  % get synapse type ID (for SynPar)
                    SPMtx(reorder{i,ii}, reorder{j,ii},k)=(X)+(k-1)*length(idc);                                % fit submatrix in

                    idc_back = idc;
                        [SynPar,Nsyn,idc]=SetSyn(SynPar,Nsyn,idc,ST,ConPar{2}{i,j}{l+2}, ...                    % set synapse parameters
                            ConPar{3}{i,j}, S_sig(i,j,:), S_max(i,j,:), S_min(i,j,:));
                    if k==1
                        dtax_back = SynPar(6,idc_back);
                    else
                        SynPar(6,idc_back) = dtax_back;
                        dtax_back = [];
                    end;
                end;
            end;
        end;

    else                            % use ConMtx if available
        
        for i=1:NTypesN                 % loop over input neurons
            for j=1:NTypesN             % loop over output neurons
                for k = 1:length(ConPar{2}{i,j}{1})                                     % loop over all (different) synapse types
                    X=ConMtx(reorder{i,ii}, reorder{j,ii},1);                           % Import connectivity submatrix for the current input->output pair
                    kk=find(X);                                                         % Extract non-zero entries...
                    idc=Nsyn+(1:length(kk));                                            % ... and create the idc vector from them
                    
                    X(kk)=idc;                                                          % Fill submatrix with idc entries...
                    X(X==0)=(k-1)*length(idc);                                          % ... and fill up the rest (for compatibitliy with former code only)
                    SPMtx(reorder{i,ii}, reorder{j,ii},k)=(X);                          % fit submatrix into full connectivity matrix
                    
                    l = ConPar{2}{i,j}{1}(k);                                           % determine which parameter set to use for each type
                    ST = ConPar{2}{i,j}{2}(k);                                          % get synapse type ID (for SynPar)

                    idc_back = idc;
                    [SynPar,Nsyn,idc]=SetSyn(SynPar,Nsyn,idc,ST,ConPar{2}{i,j}{l+2}, ...
                        ConPar{3}{i,j}, S_sig(i,j,:), S_max(i,j,:), S_min(i,j,:));      % set synapse parameters
                    if k==1
                        dtax_back = SynPar(6,idc_back);
                    else
                        SynPar(6,idc_back) = dtax_back;
                        dtax_back = [];
                    end;
                end;
            end;
        end;
    end;

    
    % Add noise synapses (one for each synapse type)
    if max(max(NoisePar))>0
        for i=1:length(NoisePar(:,1))
            idc_back = idc;
            [SynPar,Nsyn,idc]=SetSyn(SynPar,Nsyn,length(SynPar(1,:))+1,i,NoisePar(i,:), {1.0,{[1 0 0], [0 0 0]}}, zeros(1,length(NoisePar(1,:))), zeros(1,length(NoisePar(1,:))),zeros(1,length(NoisePar(1,:))));
            if k==1
                dtax_back = SynPar(6,idc_back);
            else
                SynPar(6,idc_back) = dtax_back;
                dtax_back = [];
            end;
        end;
    end
end;


% ------------------------- Define inter-stripe connections ---------------------------------

for i=1:length(ConParStripes{1})
    i_act = ConParStripes{1}(i,1);
    j_act = ConParStripes{1}(i,2);
    
    for j=1:Nstripes
        for kk=1:length(ConParStripes{3}{i})
            
            % Set target stripe, dist_act and p_act
            target_act = j + ConParStripes{3}{i}(kk);
            while target_act<=0
                target_act = target_act + Nstripes;
            end;
            while target_act>Nstripes
                target_act = target_act - Nstripes;
            end;
            dist_act = abs(ConParStripes{3}{i}(kk));
            p_act = pCon(i_act,j_act)*exp(-dist_act/ConParStripes{2}{i}(1));

            
            % insert matrix and parameter definitions (as for intra-stripe connections)
            if isempty(ConMtx)              % connect randomly if ConMtx is empty
                if cluster_flag(i_act,j_act) == 2                                                                       % set random connectivity submatrix for the current input->output pair
                    [X,idc] = SetCon_CommonNeighbour_Cross(Nsyn, SPMtx(length(reorder{i_act,target_act}),length(reorder{i_act,target_act})), SPMtx(length(reorder{j_act,j}),length(reorder{j_act,j})), p_act, 1);   % with...
                    X=X';
                else
                    [X,idc] = SetCon(Nsyn, length(reorder{i_act,target_act}), length(reorder{j_act,j}), p_act);                                                                                                     % ... or without common neighbour rule
                end;
                
                for k = 1:length(ConPar{2}{i_act,j_act}{1})                                                             % loop over all synapse types
                    l = ConPar{2}{i_act,j_act}{1}(k);                                                                   % determine which parameter set to use for each type
                    ST = ConPar{2}{i_act,j_act}{2}(k);                                                                  % get synapse type ID (for SynPar)
                    SPMtx(reorder{i_act,target_act}, reorder{j_act,j},k)=(X)+(k-1)*length(idc);                         % fit submatrix in
                    
                    % Define synapse parameters from horizontal distance
                    wgt_act = ConPar{2}{i_act,j_act}{l+2}(1)*exp(-dist_act/ConParStripes{2}{i}(1)); 
                    dtax_act = ConPar{2}{i_act,j_act}{l+2}(2) + ConParStripes{2}{i}(2)*dist_act;
                    S_sig_act(1) = S_sig(i_act,j_act,1)*exp(-dist_act/ConParStripes{2}{i}(1));
                    S_sig_act(2) = S_sig(i_act,j_act,2) + ConParStripes{2}{i}(2)*dist_act;
                    
                    idc_back = idc;
                    [SynPar,Nsyn,idc]=SetSyn(SynPar,Nsyn,idc,ST,[wgt_act dtax_act], ...                                 % set synapse parameters
                        ConPar{3}{i_act,j_act}, S_sig_act, S_max(i_act,j_act,:), S_min(i_act,j_act,:));
                    if k==1
                        dtax_back = SynPar(6,idc_back);
                    else
                        SynPar(6,idc_back) = dtax_back;
                        dtax_back = [];
                    end;
                end;
            else                            % use ConMtx if available

                for k = 1:length(ConPar{2}{i_act,j_act}{1})                                                 % loop over all (different) synapse types
                    X=ConMtx(reorder{i_act,target_act}, reorder{j_act,j},1);                                % Import connectivity submatrix for the current input->output pair
                    kkk=find(X);                                                                            % Extract non-zero entries...
                    idc=Nsyn+(1:length(kkk));                                                               % ... and create the idc vector from them
                    
                    X(kkk)=idc;                                                                             % Fill submatrix with idc entries...
                    X(X==0)=(k-1)*length(idc);                                                              % ... and fill up the rest (for compatibitliy with former code only)
                    SPMtx(reorder{i_act,target_act}, reorder{i_act,j},k)=(X);                               % fit submatrix into full connectivity matrix
                    
                    l = ConPar{2}{i_act,j_act}{1}(k);                                                       % determine which parameter set to use for each type
                    ST = ConPar{2}{i_act,j_act}{2}(k);                                                      % get synapse type ID (for SynPar)

                    % Define synapse parameters from horizontal distance
                    wgt_act = ConPar{2}{i_act,j_act}{l+2}(1)*exp(-dist_act/ConParStripes{2}{i}(1)); 
                    dtax_act = ConPar{2}{i_act,j_act}{l+2}(2) + ConParStripes{2}{i}(2)*dist_act;
                    S_sig_act(1) = S_sig(i_act,j_act,1)*exp(-dist_act/ConParStripes{2}{i}(1));
                    S_sig_act(2) = S_sig(i_act,j_act,2) + ConParStripes{2}{i}(2)*dist_act;
                    
                    idc_back = idc;
                    [SynPar,Nsyn,~]=SetSyn(SynPar,Nsyn,idc,ST,[wgt_act dtax_act], ...
                        ConPar{3}{i_act,j_act}, S_sig_act, S_max(i_act,j_act,:), S_min(i_act,j_act,:));      % set synapse parameters
                    if k==1
                        dtax_back = SynPar(6,idc_back);
                    else
                        SynPar(6,idc_back) = dtax_back;
                        dtax_back = [];
                    end;
                end;
            end;
        end;
    end;
end;


% ------------------------- Define connections to input neurons ---------------------------------

for ii=1:Nstripes
    reorder{NTypesN+1,ii} = NN+1:NN+InpNTypes(1);
    reorder{NTypesN+2,ii} = NN+InpNTypes(1)+1:NN+M;
    
    if isempty(ConMtx)              % connect randomly if ConMtx is empty
        
        % inter-celltype connections
        for i=1:NTypesN                              % loop over output neurons
            for j=NTypesN+1:NTypesN+InpNTypesN       % loop over input neurons (literally)
                if cluster_flag(i,j) == 2                                                                       % set random connectivity submatrix for the current input->output pair
                    [X,idc] = SetCon_CommonNeighbour_Cross(Nsyn, SPMtx(length(reorder{i,ii}),length(reorder{i,ii})), SPMtx(length(reorder{j,ii}),length(reorder{j,ii})), pCon(i,j), [], [], 1);  % with...
                    X=X';
                else
                    [X,idc] = SetCon(Nsyn, length(reorder{i,ii}), length(reorder{j,ii}), pCon(i,j));                                                                             % ... or without common neighbour rule
                end;
                
                for k = 1:length(ConPar{2}{i,j}{1})                                                             % loop over all synapse types
                    l = ConPar{2}{i,j}{1}(k);                                                                   % determine which parameter set to use for each type
                    ST = ConPar{2}{i,j}{2}(k);                                                                  % get synapse type ID (for SynPar)
                    SPMtx(reorder{i,ii}, reorder{j,ii},k)=(X)+(k-1)*length(idc);                                % fit submatrix in
                    
                    idc_back = idc;
                    [SynPar,Nsyn,idc]=SetSyn(SynPar,Nsyn,idc,ST,ConPar{2}{i,j}{l+2}, ...                        % set synapse parameters
                        ConPar{3}{i,j}, S_sig(i,j,:), S_max(i,j,:), S_min(i,j,:));
                    if k==1
                        dtax_back = SynPar(6,idc_back);
                    else
                        SynPar(6,idc_back) = dtax_back;
                        dtax_back = [];
                    end;
                    
                end;
            end;
        end;
        
    else                            % use ConMtx if available
        
        for i=1:NTypesN                              % loop over input neurons
            for j=NTypesN+1:NTypesN+InpNTypesN       % loop over input neurons (literally)
                for k = 1:length(ConPar{2}{i,j}{1})                                     % loop over all (different) synapse types
                    X=ConMtx(reorder{i,ii}, reorder{j,ii},1);                           % Import connectivity submatrix for the current input->output pair
                    kk=find(X);                                                         % Extract non-zero entries...
                    idc=Nsyn+(1:length(kk));                                            % ... and create the idc vector from them
                    
                    X(kk)=idc;                                                          % Fill submatrix with idc entries...
                    X(X==0)=(k-1)*length(idc);                                          % ... and fill up the rest (for compatibitliy with former code only)
                    SPMtx(reorder{i,ii}, reorder{j,ii},k)=(X);                          % fit submatrix into full connectivity matrix
                    
                    l = ConPar{2}{i,j}{1}(k);                                           % determine which parameter set to use for each type
                    ST = ConPar{2}{i,j}{2}(k);                                          % get synapse type ID (for SynPar)
                    
                    idc_back = idc;
                    [SynPar,Nsyn,~]=SetSyn(SynPar,Nsyn,idc,ST,ConPar{2}{i,j}{l+2}, ...
                        ConPar{3}{i,j}, S_sig(i,j,:), S_max(i,j,:), S_min(i,j,:));      % set synapse parameters
                    if k==1
                        dtax_back = SynPar(6,idc_back);
                    else
                        SynPar(6,idc_back) = dtax_back;
                        dtax_back = [];
                    end;
                end;
            end;
        end;
    end;
end



% ----------------------------- Define cell assemblies --------------------------------------
if ~isempty(CA)
    SPMtxCA = SPMtx(CA(:,1),CA(:,1),:);

    for i=1:length(CA(:,1))
        CA_ind = find(SPMtxCA(:,i,1)>0);
        if(~isempty(CA_ind))
            for k=1:2       % exc. synapses only
                SynPar(5,SPMtxCA(CA_ind,i,k)) = CA(i,2)*SynPar(5,SPMtxCA(CA_ind,i,k));
            end;
        end;
    end;
end;


% ------------------------- Save full parameter set  ---------------------------------
if ~isempty(SimPar.fnOut)
    EvtMtx=SimPar.EvtMtx;
    EvtTimes=SimPar.EvtTimes;
    save(SimPar.fnOut,'SimPar', 'CtrPar','NeuPar','NPList', ...
        'STypPar','SynPar','SPMtx','EvtMtx','EvtTimes','ViewList','V0', '-v7.3');
end;


% -------------- Run the actual simulation (described in IDNet.c) -------------------------

[Ninp,~]=IDNet(CtrPar,NeuPar,NPList,STypPar,SynPar,SPMtx,SimPar.EvtMtx,SimPar.EvtTimes,ViewList,InpSTtrains,NoiseDistr,V0,UniqueNum,NeuronGroupsSaveArray);


% -------------- Convert data into binary MATLAB format -------------------------
% check if there are neurons where the input oscillates around I_ref (not necessary any more)
% idx_osc=find(N_osc>0);
% if ~isempty(idx_osc)
%     disp('The following neurons might have a problem with their input (oscillates around I_ref):');
%     disp(num2str(idx_osc));
% end

% Read in temporary files and delete them
STMtx=[];
if CtrPar(5)>0
    for i=1:N*Nstripes+M
        fnST=['ISIu' num2str(i-1) '_' num2str(UniqueNum) '.dat'];
        if exist(fnST,'file')
            ST=load(fnST);
            STMtx{i}=ST';
            delete(fnST);
        else
            STMtx{i}=[];
        end;
    end;
end;

T=[]; V=[];dV=[];
if ~isempty(ViewList)
    X=load(['IDN_' num2str(UniqueNum) '.dat']);
    nv=length(ViewList);
    T=X(1:nv:end,1);
    for i=1:nv
        V{i}=X(i:nv:end,3:end);
    end;
    delete(['IDN_' num2str(UniqueNum) '.dat'])
end;

if exist(['IDN2_' num2str(UniqueNum) '.dat'],'file')
    dX=load(['IDN2_' num2str(UniqueNum) '.dat']);
    for Counter1 = 1:3*length(NeuronGroupsSaveArray)
        dV{Counter1} = dX(:,Counter1);
    end
    delete(['IDN2_' num2str(UniqueNum) '.dat'])
end

if exist(['IDN3_' num2str(UniqueNum) '.dat'],'file')
    XX=importdata(['IDN3_' num2str(UniqueNum) '.dat']);
    delete(['IDN3_' num2str(UniqueNum) '.dat'])
end


% ------------------------- Save results  ---------------------------------
if ~isempty(SimPar.fnOut)
    EvtMtx=SimPar.EvtMtx;
    EvtTimes=SimPar.EvtTimes;
    save(SimPar.fnOut,'Ninp','T','V','dV','XX','STMtx', '-append', '-v7.3');
end;



% -------------------------------------------------------------------------
% -------------- Functions to generate connectivity -----------------------
% -------------------------------------------------------------------------

function [X,idc]=SetCon(Nsyn, Nin,Nout,pCon)
% Generates a random connectivity matrix X[ConPar(1), ConPar(2)] with
% connection prob. ConPar(3). Nsyn (number of synapses) is continuated
% from input matrices. idc gives the running indices of the generated
% connections (reshaped in one dimension)
% 
% INPUT:
%   Nsyn:   Number of synapses that are already created
%   Nin:    Number of input neurons
%   Nout:   Number of output neurons
%   pcon:   Connection probability
% 
% OUTPUT:
%   X:      Connection matrix (Nin x Nout matrix) containing idc on 
%           non-zero entries
%   idc:    Synapse indices (starting with Nsyn)


k1 = randperm(Nin*Nout);
k = k1(1:round(pCon*Nin*Nout));
X=zeros(Nin,Nout);
idc=Nsyn+(1:length(k));
X(k)=idc;


function [SynPar,Nsyn,idc]=SetSyn(SynPar,Nsyn,idc,ST,ConPar,ConParSTSP,S_sig,S_max,S_min)
% Sets (i.e. continuates) synaptic parameter matrix SynPar(6,idc), where
% idc is the running index of connections (see SetCon). 
% 
% INPUT:
%   SynPar:         Synapse parameters (excluding those computed here)
%   Nsyn:           Number of synapses that are already set
%   idc:            Indices of synapses that are already set
%   ST:             Synapse type
%   ConPar:         Mean connection parameters (synaptic weight, delay and 
%                   failure probability)
%   ConParSTSP:     Short-term synaptic plasticity parameters
%   S_sig:          Standard deviation of connection parameters
%   S_max:          Maximum of connection parameters
%   S_min:          Minimum of connection parameters
% 
% OUTPUT:
%   SynPar:         Synapse parameters (concatinating old and new columns)
%   Nsyn:           Updated number of synapses
%   idc:            Updated synapse indices

wgt=ConPar(1)';
dtax=ConPar(2);
SynPar(1,idc) = ST;                                                        % Synapse type (AMPA, GABA and NMDA here)


mean_wgt = log((wgt^2)/sqrt(S_sig(1)^2+wgt^2));
std_wgt = sqrt(log(S_sig(1)^2/(wgt^2)+1));
if wgt==0 || S_sig(1)==0
    SynPar(5,idc) = rand_par(length(idc), wgt, S_sig(1),   S_min(1)*wgt,     S_max(1)*wgt,    0);            % wgt
else
    SynPar(5,idc) = rand_par(length(idc), mean_wgt,    std_wgt,    S_min(1)*wgt,     S_max(1)*wgt,    2);    % wgt
end

SynPar(6,idc) = rand_par(length(idc), dtax,        S_sig(2),   S_min(2)*dtax,    S_max(2)*dtax,   0);        % dtax
SynPar(7,idc) = ConPar(3);                                                                                   % p_fail  

if ~isempty(idc)
    % Deal with rounding errors
    idc_frac = round(ConParSTSP{1}*length(idc));
    if sum(idc_frac) > length(idc)
        [~,ind] = min(ConParSTSP{1}*length(idc) - floor(ConParSTSP{1}*length(idc)));
        idc_frac(ind) = idc_frac(ind)-1;
    else
        if sum(idc_frac) < length(idc)
            [~,ind] = max(idc_frac - floor(idc_frac));
            idc_frac(ind) = idc_frac(ind)+1;
        end;
    end;
    
    idc_act = Nsyn+(1:idc_frac(1));
    Nsyn_act = Nsyn+idc_frac(1);
    for i=1:length(ConParSTSP{1})
        use    = ConParSTSP{i+1}{1}(1);
        tc_rec = ConParSTSP{i+1}{1}(2);
        tc_fac = ConParSTSP{i+1}{1}(3);
        std_use    = ConParSTSP{i+1}{2}(1);
        std_tc_rec = ConParSTSP{i+1}{2}(2);
        std_tc_fac = ConParSTSP{i+1}{2}(3);
        
        % set maxima and standard deviations to zero before using update_inv_con_PSP
        SynPar(2,idc_act) = rand_par(length(idc_act), use,    std_use,    0,    1, 0);    % use
        SynPar(3,idc_act) = rand_par(length(idc_act), tc_rec, std_tc_rec, 0, 1500, 0);    % tc_rec
        SynPar(4,idc_act) = rand_par(length(idc_act), tc_fac, std_tc_fac, 0, 1500, 0);    % tc_fac
        
        if i<length(ConParSTSP{1}) && length(ConParSTSP{1})>1
            idc_act=Nsyn_act+(1:idc_frac(i+1));
            Nsyn_act=Nsyn_act+idc_frac(i+1);
        end;
    end;
end;

Nsyn=Nsyn+length(idc);
idc=Nsyn+(1:length(idc));


% (c) 2016 J. Hass, L. Hertaeg and D. Durstewitz,
% Central Institute of Mental Health, Mannheim University of Heidelberg 
% and BCCN Heidelberg-Mannheim