function output=model2018(sign,fs,fc,irregularities,storeflag,subject,sheraPo,IrrPct,non_linear_type,nH,nM,nL,clean,data_folder)

% An fiber model
% input:
%   -sign: stimulus
%   -fs: samplerate
%   -fc: probe frequency or alternatively a string 'all' to probe all cochlear
%        sections or 'half' to probe half sections, 'abr' to store the
%        401 sections used to compute the abr responses in Verhulst et al. 2017
%   -irregularities: decide wether turn on (1) or off (0) irregularities and
%                  nonlinearities of the cochlear model (default 1 )
%   -storeflag: string that sets what variables to store from the
%   computation, each letter correspond to one desired output variable
%   (e.g. 'avhl' to store acceleration, displacement,  high and low spont.
%   fibers.) See "output" for the corresponding characters. Default
%   'avihmlebw'.
%   -subject: number representing the seed to generate the random
%       irregularities in the cochlear sections (default 1)
%   -sheraPo: starting real part of the poles of the cochlear model
%       it can be either an array with one value per BM section,  or a
%       single value far all section (default 0.6)
%   -IrrPct: magnitude of random perturbations on the BM (irregularities, default 0.05=5%)
%   -non_linear_type: select the type of nonlinearity in the BM model.
%       Currently implemented
%           'vel'= instantaneouns nonlinearity based on local BM velocity (see Verhulst et al. 2012)
%           'none'= linear model
%   -nH,nM,nL= number of high, medium and low spont. fibers employed to
%   compute the response of cn and ic nuclei. Default 13,3,3. This
%   parameters can be passed either as a single value for all sections or
%   as an array with each value corresponding to a single CF location (and
%   the length of those array must match the number of the section probed)

%  output:
%   -output.v: BM velocity     (store 'v')
%   -output.y: BM displacement (store 'y')
%   -output.emission: pressure output from the middle ear (store 'e')
%   -output.cf: center frequencies (always stored)
%   -output.ihc: IHC receptor potential (store 'i')
%   -output.anfH: HSR fiber spike probability [0,1] (store 'h')
%   -output.anfM: MSR fiber spike probability [0,1] (store 'm')
%   -output.anfL: LSR fiber spike probability [0,1] (store 'l')
%   -output.an_sum: summation of HSR, MSR and LSR per channel,i.e., the input to the cn model  (storeflag 'b')
%   -output.cn: cochlear nuclei output (storeflag 'b')
%   -output.ic: IC (storeflag 'b')
%   -output.w1,outupt.w3,output.w5: wave 1,3 and 5 (storeflag 'w')
%   -output.fs_bm= sampling frequency of the bm simulations
%   -output.fs_ihc= sample rate of the inner hair cell output
%   -output.fs_an= sample rate of the an output
%   -output.fs_abr= sample rate of the IC,CN and W1/3/5 outputs
%   -clean: delete the temporary files generated by the python call (1) or
%            not (0). Default 1
%   -data_folder: string of folder where to save the temporary files
%               generated by the python call (default current folder)
DECIMATION=5;
[channels,idx]=min(size(sign));
if nargin<13
    data_folder=[pwd(),'/'];
end
if nargin<13
    clean=1;
end
if nargin < 12
    nL=3;
end
if nargin < 11
    nM = 3;
end
if nargin < 10
    nH = 13;
end
if nargin < 9
    non_linear_type = 'vel';
end
if nargin < 8
    IrrPct = 0.05;
end
if nargin < 7
    sheraPo = 0.06;
end
if nargin < 6
    subject = 1;
end
if nargin < 5
    storeflag = 'vihlmeb';
end
if nargin< 4
    irregularities=ones(1,channels);
end


if(numel(irregularities)==1)
    irregularities=irregularities*ones(1,channels);
end

Fs=fs;
sectionsNo=1e3;
[channels,idx]=min(size(sign));
if(idx==2) %transpose it (python C-style row major order)
    sign=sign';
end
stim=sign;
% for i=1:channels
%     stim(i,:)=resample(sign(i,:),Fs,fs);
% end
% sheraPo=0.060;
if(isstr(fc) && strcmp(fc,'all')) %if probing all sections 1001 output (1000 sections plus the middle ear)
    l=sectionsNo;
elseif(isstr(fc) && strcmp(fc,'half')) %if probing half sections sections 1001 output (1000 sections plus the middle ear)
    l=sectionsNo/2;
elseif(isstr(fc) && strcmp(fc,'abr')) %if probing half sections sections 1001 output (1000 sections plus the middle ear)
    l=401;
else %else pass it as a column vector
    [l,idx]=max(size(fc));
    if(idx==2)
        fc=fc'; 
    end
    fc=round(fc);
end
probes=fc;

act_path=pwd;

save('input.mat','stim','Fs','channels','subject','sheraPo','irregularities','probes',...
    'sectionsNo','data_folder','storeflag','IrrPct','non_linear_type','nH','nM','nL','-v7');
system('python run_model2018.py');

cd(data_folder);

for i=1:channels

    p=length(stim(i,:));
    p2=ceil(length(stim(i,:))/DECIMATION);
    
    fname = strcat('cf',int2str(i),'.mat');
    tmp = load(fname,'cf');
    output(i).cf = tmp.cf;

    if strfind(storeflag,'v')
        fname = strcat('v',int2str(i),'.mat');
        tmp = load(fname,'Vsolution');
        output(i).v = tmp.Vsolution;
    end

    if strfind(storeflag,'y')
        fname = strcat('y',int2str(i),'.mat');
        tmp = load(fname,'Ysolution');
        output(i).y = tmp.Ysolution;
    end
    
    if strfind(storeflag,'e')
        fname = strcat('emission',int2str(i),'.mat');
        tmp = load(fname,'oto_emission');
        output(i).e = tmp.oto_emission;
    end

    if strfind(storeflag,'i')
        fname = strcat('ihc',int2str(i),'.mat');
        tmp = load(fname,'Vm');
        output(i).ihc = tmp.Vm;
    end

    if strfind(storeflag,'h')
        fname = strcat('anfH',int2str(i),'.mat');
        tmp = load(fname,'anfH');
        output(i).anfH = tmp.anfH;
    end
    
    if strfind(storeflag,'m')
        fname = strcat('anfM',int2str(i),'.mat');
        tmp = load(fname,'anfM');
        output(i).anfM = tmp.anfM;
    end

    if strfind(storeflag,'l')
        fname = strcat('anfL',int2str(i),'.mat');
        tmp = load(fname,'anfL');
        output(i).anfL = tmp.anfL;
    end

    if strfind(storeflag,'b')
        fname = strcat('cn',int2str(i),'.mat');
        tmp = load(fname,'cn');
        output(i).cn = tmp.cn;
        
        fname = strcat('ic',int2str(i),'.mat');
        tmp = load(fname,'ic');
        output(i).ic = tmp.ic;
        
        fname = strcat('AN',int2str(i),'.mat');
        tmp = load(fname,'anSummed');
        output(i).an_summed = tmp.anSummed;
    end

    if strfind(storeflag,'w')
        fname = strcat('1w',int2str(i),'.mat');
        tmp = load(fname,'w1');
        output(i).w1 = tmp.w1;

        fname = strcat('3w',int2str(i),'.mat');
        tmp = load(fname,'w3');
        output(i).w3 = tmp.w3;
        
        fname = strcat('5w',int2str(i),'.mat');
        tmp = load(fname,'w5');
        output(i).w5 = tmp.w5;
    end

    output(i).fs_abr=fs/DECIMATION;
    output(i).fs_an=fs/DECIMATION;
    output(i).fs_bm=fs;
    output(i).fs_ihc=fs;
end

if clean==1
    delete *.mat
end
save('input.mat','stim','Fs','channels','subject','sheraPo','irregularities','probes',...
    'sectionsNo','data_folder','storeflag','IrrPct','non_linear_type','nH','nM','nL','-v7');

cd(act_path);

%The model code and interface was written by Alessandro Altoč and Sarah Verhulst (copyright 2012,2014,2015,2016,2018) 
%and is licensed under the UGent acadamic license (see details in license file that is part of this repository). 
%The Verhulstetal2018Model consists of the following files: 
%tridiag.so, cochlea_utils.c, run_model2018.py, model2018.m, cochlear_model2017.py, inner_hair_cell2018.py, auditory_nerve2017.py, ic_cn2017.py, ExampleSimulation.m, ExampleAnalysis.m, the HI profiles in the Poles folder.