function [output,cf] = verhulst2018(insig,fs,fc_flag,varargin)
%VERHULST2018 Cochlear transmission-line model (improved, incl. brainstem)
%
%   Usage: 
%     output = verhulst2018(insig,fs,fc_flag)
%     [output, cf] = verhulst2018(insig,fs,fc_flag, ...)
%
%   Input parameters:
%
%     insig          : the input signal to be processed. Each column is 
%                      processed in parallel, so it is possible to run several simulations in parallel
%     fs             : sampling rate [Hz]
%     fc_flag        : list of frequencies specifying the probe positions 
%                      along the basilar membrane, or 'all'to probe all 
%                      1000 cochlear sections, or 'abr' to probe 401 locations
%                      between 112 and 12000 Hz.
%
%   Output parameters:
%
%     cf : center frequencies of the probed basiliar membrane sections
%     output : structure with the following fields:
%     v  : velocity of the basilar membrane sections v(time,section,channel)
%     y  : displacement of the basilar membrane sections y(time,section,channel)
%     oae : ottoacoustic emission as sound pressure at the middle ear
%     fs_an : sample rate of the an output
%     fs_abr : sample rate of the brainstem sections (IC,CN and W1/3/5 outputs)
%     ic : output of the IC
%     cn : output of the CN
%     anfH : HSR fiber spike probability [0,1]
%     anfM : MSR fiber spike probability [0,1]
%     anfL : LSR fiber spike probability [0,1]
%     an_summed : sum of HSR, MSR and LSR per channel, i.e., the input to the CN
%     ihc : IHC receptor potential
%     w1: wave 1
%     w3 : wave 3
%     w5 : wave 5
%
%   References:
%     S. Verhulst, A. Altoè, and V. Vasilkov. Functional modeling of the
%     human auditory brainstem response to broadband stimulation.
%     hearingresearch, 360:55--75, 2018.
%     
%
%   License:
%   --------
%
%   This model is licensed under the UGent Academic License. Further usage details are provided 
%   in the UGent Academic License which can be found in the AMT directory "licences" and at 
%   <https://raw.githubusercontent.com/HearingTechnology/Verhulstetal2018Model/master/license.txt>.
%
%   See also: verhulst2015 verhulst2018 demo_verhulst2018 demo_verhulst2012
%             demo_verhulst2015 verhulst2018_ihctransduction verhulst2015_cn
%             verhulst2015_ic verhulst2018_auditorynerve exp_verhulst2012
%             exp_verhulst2018 verhulst2012 verhulst2015
%             verhulst2018 middleearfilter data_takanen2013 takanen2013_periphery
%             exp_osses2022 exp_takanen2013 takanen2013
%
%
%
%
%   Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/models/verhulst2018.php


%   #License: ugent
%   #StatusDoc: Good
%   #StatusCode: Good
%   #Verification: Unknown
%   #Requirements: MATLAB M-Signal PYTHON C
%   #Author: Alejandro Osses (2020): primary implementation based on https://github.com/HearingTechnology/Verhulstetal2018Model
%   #Author: Piotr Majdak (2021): adaptations for the AMT 1.0

amt_info('once');
% ------ Checking of input parameters ------------

if nargin<3
	error('%s: Too few input arguments.',upper(mfilename));
end;

if ~isnumeric(insig)
	error('%s: insig must be numeric.',upper(mfilename));
end;

if ~isnumeric(fs) || ~isscalar(fs) || fs<=0
    error('%s: fs must be a positive scalar.',upper(mfilename));
end;

definput.import={'verhulst2018'};  % load defaults from arg_verhulst2018
[flags,keyvals]  = ltfatarghelper([],definput,varargin);

IrrPct = keyvals.IrrPct;
storeflag = ' ';
if flags.do_oae, storeflag = [storeflag 'e']; end
if flags.do_debug, storeflag = [storeflag 'y']; end
if flags.do_ihc, storeflag = [storeflag 'i']; end

[Nr_signals,idx]=min(size(insig));
irregularities=keyvals.irr_on.*ones(1,Nr_signals);

%%% Fix model parameters so far: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subject = keyvals.subject; % seed
non_linear_type = keyvals.non_linear_type;

% Number of neurones: 13-3-3 => no synaptopathy, any loss of neurones => synaptopathy
numH = keyvals.numH;
numM = keyvals.numM;
numL = keyvals.numL;

if idx == 1
    error('If multiple signals are used as input make sure each column is one signal')
end

%%% Stage 1A. Outer ear: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if flags.do_outerear
    % Same outer ear filter as in other peripheral models
    hp_fir = headphonefilter(fs);% Getting the filter coefficients at fs
    N = ceil(length(hp_fir)/2);  % group delay for a FIR filter of order length(hp_fir)
    M = size(insig,2);
    insig = [insig; zeros(N,M)]; % group delay compensation: step 1 of 2.
    insig = filter(hp_fir,1,insig); % filtering
    insig = insig(N+1:end,1:M); % group delay compensation: step 2 of 2
end

%%% Stage 1B. Middle ear: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[b,a] = middleearfilter(fs,'verhulst2018');
if flags.do_middleear
    insig = filter(b,a,insig);
end
if flags.do_no_middleear
    % If no middle-ear filter is applied, a gain/attenuation equivalent to
    % the level in the filter band-pass is applied to the input signal. This
    % is important to get an appropriate compression in the subsequent filter
    % bank.
    K = 8192; % arbitrary
    h = ones(K,1);
    h = h.*freqz(b,a,K);

    me_gain_TF = max( 20*log10(abs(h)) ); % max of the filter response
    insig = gaindb(insig,me_gain_TF);
end

%%% Stage 2. Cochlear filter bank - Transmission line model: %%%%%%%%%%%%%%
%         2.1 Signal preprocessing:
fs_in = fs;
if fs ~= keyvals.fs_up
    insig = resample(insig,keyvals.fs_up,fs);
    fs = keyvals.fs_up;
end
stim=transpose(insig); % transpose it (python C-style row major order)

if(ischar(fc_flag)) %if probing all sections 1001 output (1000 sections plus the middle ear)
    switch  fc_flag
        case 'all'
            Nr_sections = 1000;
        case 'half'
            Nr_sections = 500;
        case 'abr'
            Nr_sections = 401;
            % Nothing to do, it is correctly formatted
        otherwise
            error('fc flag not recognised it should be either ''all'',''abr'',''half'', or a numeric array');
    end
    fc_str = fc_flag;
else
    % then fc is a numeric array that will be passed as a column vector
    fc_flag=round(fc_flag(:));
    fc_str = 'custom CF(s)';
    Nr_sections = length(fc_flag);
end
probes=fc_flag;

  % Load poles for the selected hearing profile
poles=amt_load('verhulst2018','Poles.mat','Poles');
if ~isfield(poles.Poles,keyvals.hearing_profile)
  error(['Hearing profile ' keyvals.hearing_profile ' not available.']);
end
sheraPo=poles.Poles.(keyvals.hearing_profile);

dir_model = fullfile(amt_basepath,'environments','verhulst2018',filesep);
dir_data  = fullfile(dir_model,'out',filesep);
version_year = 2018; 

% if ~exist(fullfile(dir_model,'tridiag.so'),'file')
% 	error('/environments/verhulst2018/tridiag.so library is missing. Run amt_mex');
% end

channels = Nr_signals;
L_samples = length(stim(1,:));

amt_disp('VERHULST 2018: The following parameters will be passed to Python:',flags.disp);
amt_disp(['  Number of signals to be processed: ',num2str(channels),' channels'],flags.disp);
amt_disp(['  Cochlear hearing profile: ',keyvals.hearing_profile,' (poles between ',num2str(sheraPo(1)),' and ',num2str(sheraPo(end)),')'],flags.disp);
amt_disp(['  Number of auditory nerve fibres: (',num2str([numH,numM,numL]),') (HSR,MSR,LSR)'],flags.disp);
amt_disp(['  Number of cochlear sections to be stored: flag ',num2str(fc_str),' (all=1000, half=500, abr=401)'],flags.disp);
amt_disp(['  Seed number: ',num2str(subject),' (subject ''variable'')'],flags.disp);
amt_disp(['  Irregularities=',num2str(irregularities(1)),' (1=on,0=off), IrrPrct=',num2str(IrrPct)],flags.disp);
amt_disp(['  Non linear type=',num2str(non_linear_type),'(''vel'' or ''disp'')'],flags.disp);
amt_disp(['  Extra variables to be saved: ',num2str(storeflag),' (storeflag)'],flags.disp);

%%% Storing input, change dir, run the model, and come back
amt_disp('Cochlear processing...','volatile');
% act_path=pwd;
% cd(dir_model);
% fname = 'input.mat';
% save(fname,'stim','fs','channels','subject','sheraPo','irregularities','probes', ... % this line: same parameters as verhulst2012.m
%     'dir_data','storeflag','IrrPct','non_linear_type','version_year','-v7');
% [status,res] = system(sprintf('python %srun_cochlear_model.py',dir_model),'-echo');
% cd(act_path);
% if status ~= 0
%     %%% Then the simulation did not succeed
% 	amt_disp();
%     amt_disp(res); % it re-prints the Python warning...
%     error('Something went wrong (see the Python warning above)')
% end

in.stim=stim; in.fs=fs; in.channels=channels; in.subject=subject;
in.sheraPo=sheraPo; in.irregularities=irregularities;
in.probes=probes; in.storeflag=storeflag; in.IrrPct=IrrPct;
% in.dir_data=dir_data; 
in.non_linear_type=non_linear_type; in.version_year=version_year;
out.cf=[Nr_sections,1,Nr_signals];
out.v=[Nr_sections,L_samples,Nr_signals];
if flags.do_debug, out.y=[Nr_sections,L_samples,Nr_signals]; end
if flags.do_oae, out.e=[L_samples,1,Nr_signals]; end
out=amt_extern('Python','verhulst2018','run_cochlear_model.py',in,out); 

% if status ~= 0
%     %%% Then the simulation did not succeed
% 	amt_disp();
%     amt_disp(res); % it re-prints the Python warning...
%     error('Something went wrong (see the Python warning above)')
% end


% output = [];
% fname2clean = []; % temporal files will be removed after simulating
% bClean = 1;
% 
% for i=1:Nr_signals
% 
%     fname=[dir_data 'cf',int2str(i),'.np'];
%     if i == 1
%         % Characteristic frequencies are read only once (the same for all Nr_signals)
%         f=fopen(fname,'r');
%         cf=fread(f,[Nr_sections,1],'double','n');
%         fclose(f);
%     end
%     if bClean; fname2clean{end+1} = fname; end
% 
%     fname=[dir_data 'v' int2str(i) '.np'];
%     f=fopen(fname,'r');
%     output(i).v=fread(f,[Nr_sections,L_samples],'double','n')';
%     fclose(f);
%     if bClean; fname2clean{end+1} = fname; end
% 
%     if flags.do_debug
%         fname=[dir_data 'y' int2str(i) '.np'];
%         f=fopen(fname,'r');
%         output(i).ys=fread(f,[Nr_sections,L_samples],'double','n');
%         fclose(f);
%         if bClean; fname2clean{end+1} = fname; end
%     end
% 
%     if flags.do_oae
%         fname=[dir_data 'e',int2str(i),'.np'];
%         f=fopen(fname,'r');
%         output(i).e=fread(f,[L_samples,1],'double','n');
%         fclose(f);
%         if bClean; fname2clean{end+1} = fname; end
%     end
%     output(i).fs_bm=fs_in;
% end

fs_abr = keyvals.subfs; % default subfs = 20000 Hz
fs_an  = keyvals.subfs;
DECIMATION = fs/keyvals.subfs;

output(1:Nr_signals) = struct('fs_bm',fs_in);

for ii=1:Nr_signals
  
    amt_disp(['Neural processing ' num2str(ii) ' of ' num2str(Nr_signals)],'volatile');

%     output(ii).fs_bm=fs_in;
    output(ii).cf=squeeze(out.cf(:,1,ii));
    output(ii).v=squeeze(out.v(:,:,ii))';
    if flags.do_debug, output(ii).y=resample(squeeze(out.y(:,:,ii))',fs_in,fs); end
    if flags.do_oae, output(ii).oae=resample(squeeze(out.e(:,:,ii)),fs_in,fs); end
    
    
    %%% Stage 3. Inner hair cell model: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    outsig = verhulst2018_ihctransduction(output(ii).v,fs,version_year,keyvals); % unscaled input

    if flags.do_ihc, output(ii).ihc = outsig; end

    % Reducing the sampling frequency: Decimation.
    n = 30;
    for jj = 1:size(outsig,2)
        Vm_res(:,jj) = decimate(outsig(:,jj),DECIMATION,n,'fir'); % n-th order FIR filter before decimation
    end
    Vm_res(1:DECIMATION,:) = repmat(outsig(1,:),DECIMATION,1);

    % Now that Vm_res has been decimated we can reduce the sampling frequency
    % of v and ihc (if do_ihc):
    output(ii).v = resample(output(ii).v,fs_in,fs);
    if flags.do_ihc
        output(ii).ihc = resample(output(ii).ihc,fs_in,fs);
        output(ii).fs_ihc=fs_in;
    end

    if flags.do_an || flags.do_cn || flags.do_ic
        %%% Stage 4: Auditory nerve model and Wave I: %%%%%%%%%%%%%%%%%%%%%%%%%
        % HSR neurones
        if numH ~= 0
            kSR  = keyvals.kSR_H;  %   68.5 spikes/s
            kmax = keyvals.kmax_H; % 3000 spikes/s, peak exocytosis rate
            anfH = verhulst2018_auditorynerve(Vm_res,fs_abr,kSR,kmax,version_year);
        else
            error('There should be at least one HSR neurone, set numH to a non-null value...')
        end

        % MSR neurones
        if numM ~= 0
            kSR  = keyvals.kSR_M;  %   10 spikes/s
            kmax = keyvals.kmax_M; % 1000 spikes/s, peak exocytosis rate
            anfM = verhulst2018_auditorynerve(Vm_res,fs_abr,kSR,kmax,version_year);
        else
            anfM = zeros(size(Vm_res));
        end
        % LSR neurones
        if numL ~= 0
            kSR  = keyvals.kSR_L;  %    1 spike/s
            kmax = keyvals.kmax_L; %  800 spikes/s, peak exocytosis rate
            anfL = verhulst2018_auditorynerve(Vm_res,fs_abr,kSR,kmax,version_year);
        else
            anfL = zeros(size(Vm_res)); % empty array if no anfL, saves some computation power
        end

        if flags.do_anfH
            output(ii).anfH = anfH;
        end
        if flags.do_anfM
            output(ii).anfM = anfM;
        end
        if flags.do_anfL
            output(ii).anfL = anfL;
        end

        outsig = numL*anfL+numM*anfM+numH*anfH;
        if flags.do_an
            output(ii).an_summed = outsig;
        end

        % Loading scaling constants for Wave I, III, and V:
        model_version = keyvals.model_version;

        M1 = keyvals.M1;
        M3 = keyvals.M3;
        M5 = keyvals.M5;

        switch Nr_sections
            case {401,500}
                % same cochlear tuning in both configurations:
                cal_factor = 1;
            case 1000
                % more ''dense'' cochlear resolution (twice as many channels):
                cal_factor = 0.5;
            otherwise
                cal_factor = 1;
        end

        output(ii).w1 = cal_factor*M1*sum(outsig,2);

        %%% Stage 5: Cochlear nucleus: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Tex = keyvals.Tex_cn;
        Tin = keyvals.Tin_cn;
        dly = keyvals.dly_cn;
        A   = keyvals.Acn;
        S   = keyvals.Scn;
        outsig = verhulst2015_cn(outsig,fs_abr,Tex,Tin,dly,A,S);
        if flags.do_cn
            if ~iscell(outsig)
                output(ii).cn = outsig;
            else
                output(ii).cn_mfb = outsig;
                output(ii).cn = il_sum_cell(outsig);
            end
        end

        if flags.do_no_mfb
            % Only one CN filter, 'outsig' is numeric:
            output(ii).w3 = cal_factor*M3*sum(outsig,2);
        end
        if flags.do_mfb
            % CN from the modulation filter bank: outsig is a cell variable
            output(ii).w3 = cal_factor*M3*sum(il_sum_cell(outsig),2);
        end

        %%% Stage 6: Inferior Colliculus: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Tex = keyvals.Tex_ic;
        Tin = keyvals.Tin_ic;
        dly = keyvals.dly_ic;
        A   = keyvals.Aic;
        S   = keyvals.Sic;
        outsig = verhulst2015_ic(outsig,fs_abr,Tex,Tin,dly,A,S);
        if flags.do_ic
            if ~iscell(outsig)
                output(ii).ic = outsig;
            else
                output(ii).ic_mfb = outsig;
                output(ii).ic = il_sum_cell(outsig);
            end

        end
        if flags.do_no_mfb
            % Only one IC filter, 'outsig' is numeric:
            output(ii).w5 = cal_factor*M5*sum(outsig,2);
        end
        if flags.do_mfb
            % IC from the modulation filter bank: outsig is a cell variable
            output(ii).w5 = cal_factor*M5*sum(il_sum_cell(outsig),2);
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    output(ii).fs_an  = fs_an;
    output(ii).fs_abr = fs_abr;
%     output(ii).cf     = cf;
end
output(1).keyvals = keyvals;

% if bClean==1
%     for i = 1:length(fname2clean)
%         try
%             delete(fname2clean{i});
%         end
%     end
% end
% cd(act_path);

amt_disp();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function inoutsig = il_sum_cell(inoutsig)

for i = 2:length(inoutsig)

    inoutsig{1} = inoutsig{1}+inoutsig{i};

end
inoutsig = inoutsig{1}; % keeps only the first cell, and converts it into a double array