function [output,cf] = verhulst2018(insig,fs,fc_flag,varargin)
%verhulst2018 Cochlear transmission-line model including a model of the brainstem
%
%   Usage: 
%     output = verhulst2018(insig,fs,fc_flag)
%     [output, cf] = verhulst2018(insig,fs,fc_flag, ...)
%
%   Input parameters:
%
%     insig          : Acoustic signal (in Pa). Size: (*time x ear*).
%     fs             : Sampling rate (in Hz).
%     fc_flag        : Either a vector with frequencies specifying the probe positions
%                      along basilar-membrane sections, or 'all' to probe all
%                      1000 cochlear sections, or 'abr' to probe 401 locations
%                      between 112 and 12000 Hz.
%
%
%   Output parameters:
%
%     cf        : Center frequencies (in Hz) of the probed basiliar membrane sections.
%     output    : Structure with the following fields:
%
%                 - fs_an*: Sample rate (in Hz) of the output.
%
%                 - fs_abr*: Sample rate (in Hz) of the brainstem sections (IC, CN, W1, W3, and W5).
%
%                 - w1*: Wave 1, output of the AN model.
%
%                 - w3*: Wave 3, output of the CN model.
%
%                 - w5*: Wave 5, output of the IC model.
%
%                 - an_summed*: Sum of HSR, MSR and LSR responses (per ear*) and the input 
%                   to the CN (modelled by VERHULST2015_CN). Provided by default. 
%                   Can be disabled by the flag 'no_an'.
%
%                 - ihc*: IHC receptor potential. Provided by default. 
%                   Can be disabled by the flag 'no_ihc'.
%
%                 - cn*: Detailed output of the CN. Provided by default. 
%                   Can be disabled by the flag 'no_cn'.
%
%                 - ic*: Detailed output of the IC. Provided by default. 
%                   Can be disabled by the flag 'no_ic'.
%
%
%   VERHULST2018(..) computes the basilar membrane displacement. 
%   The following additinal flags can be specified:
%
%     'an'        Calculate the AN responses. Default.
%
%     'no_an'     Disable the caclculation of the AN and following stages.
%
%     'cn'        Calculate the CN responses. Default.
%
%     'no_cn'     Disable the caclculation of the CN and following stages.
%
%     'ic'        Calculate the IC responses. Default.
%
%     'no_ic'     Disable the caclulation of the IC and following stages.
%
%     'no_mfb'    Determines the output of W3 and W5: The output of only a single 
%                 CN filter is returned, output.w3 and output.w5 are numeric. Default.
%
%     'mfb'       Determines the output of W3 and W5: The output of the whole modulation filterbank 
%                 is returned, output.w3 and output.w5 are cell arrays. 
%
%     'progress'  Display the calculation progress. Default. 
%
%     'silent'    No progress display. Opposite of 'progress'. 
%
%     'no_debug'  No display of debugging information. Default. 
%
%     'debug'     Display additional debugging information. 
%                 Opposite of 'no_debug'.
%
%
%   VERHULST2018(..) takes also many additional parameters, see arg_VERHULST2018 for more details. 
%
%
%   The output can optionally provide the following information:
%
%     anfH*: Responses of the high-SR fibers. Optional, only when called with flag 'anfH'.
%
%     anfM*: Responses of the medium-SR fibers. Optional, only when called with flag 'anfM'.
%
%     anfL*: Responses of the low-SR fibers. Optional, only when called with flag 'anfL'.
%
%     v*: Optional velocity of the basilar membrane sections. Size: (*time x section x ear*). This is 
%     also the input to the AN modelled by VERHULST2018_IHCTRANSDUCTION. Provided only
%     when called with flag 'v'.
%
%     y*: Optional displacement of the basilar membrane sections. Size: (*time x section x ear*). 
%     Provided only when called with flag 'y'.
%
%     oae*: Optional ottoacoustic emission as sound pressure at the middle ear. Provided only
%     when called with flag 'oae'.
%
%   License
%   --------
%
%   This model is licensed under the UGent Academic License. For non-commercial academic research, 
%   you can use this file and/or modify it under the terms of that license. Further usage details 
%   are provided in the in the AMT directory "licences".
%
%   See also: demo_verhulst2018 
%             verhulst2018_ihctransduction verhulst2015_cn
%             verhulst2015_ic verhulst2018_auditorynerve 
%             middleearfilter 
%
%   References:
%     S. Verhulst, A. Altoè, and V. Vasilkov. Functional modeling of the
%     human auditory brainstem response to broadband stimulation. Hearing
%     Research, 360:55--75, 2018.
%     
%
%
%   Url: http://amtoolbox.org/amt-1.6.0/doc/models/verhulst2018.php


%   #License: ugent
%   #StatusDoc: Perfect
%   #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.
%   #Author: Piotr Majdak (2024): Documentation improvements.

% This file is licenced under the terms of the UGent Academic License, which details can be found in the AMT directory "licences" and at <https://raw.githubusercontent.com/HearingTechnology/Verhulstetal2018Model/master/license.txt>.
% For non-commercial academic research, you can use this file and/or modify it under the terms of that license. This file is distributed without any warranty; without even the implied warranty of merchantability or fitness for a particular purpose. 

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_y  , 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
if ~flags.do_silent
  amt_disp(''); % start of volatile display
  amt_disp('Cochlear processing...','volatile');
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.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_y,   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); 

v = out.v;

if flags.do_no_v, out = rmfield(out,'v'); end % reduce memory usage if v not used.

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
    if ~flags.do_silent
      amt_disp(['Neural processing ' num2str(ii) ' of ' num2str(Nr_signals) '...'],'volatile');
    end

    v_here = squeeze(v(:,:,ii))';

    output(ii).cf=squeeze(out.cf(:,1,ii));
    if flags.do_v,   output(ii).v=resample(v_here,fs_in,fs); end
    if flags.do_y  , 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: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if flags.do_ihc || flags.do_an || flags.do_cn || flags.do_ic
        outsig = verhulst2018_ihctransduction(v_here,fs,version_year,keyvals); % unscaled input
        output(ii).ihc = resample(outsig,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: %%%%%%%%%%%%%%%%%%%%%%%%%
        % 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);
      
      if numH ~= 0 % HSR neurones
          anfH = verhulst2018_auditorynerve(Vm_res,fs_abr,keyvals.kSR_H,keyvals.kmax_H,version_year);
      else
          error('There should be at least one HSR neurone, set numH to a non-null value...')
      end
      
      if numM ~= 0 % MSR neurones
          anfM = verhulst2018_auditorynerve(Vm_res,fs_abr,keyvals.kSR_M,keyvals.kmax_M,version_year);
      else
          anfM = zeros(size(Vm_res));
      end
        
      if numL ~= 0 % LSR neurones
          anfL = verhulst2018_auditorynerve(Vm_res,fs_abr,keyvals.kSR_L,keyvals.kmax_L,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:
      M1 = keyvals.M1;
      M3 = keyvals.M3;
      M5 = keyvals.M5;

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

      output(ii).w1 = cal_factor*M1*sum(outsig,2);
    end
    
    if flags.do_cn || flags.do_ic
      %%% 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 = local_sumcell(outsig);
        end
      end

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

    if flags.do_ic
      %%% 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 = local_sumcell(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(local_sumcell(outsig),2);
      end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    output(ii).fs_an  = fs_an;
    output(ii).fs_abr = fs_abr;
end

output(1).keyvals = keyvals;
amt_disp(); %end of volatile display

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function inoutsig = local_sumcell(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