function [V,Y,OAE,CF]=verhulst2012(sign,fs,fc,spl,varargin)
%VERHULST2012 Process a signal with the cochlear model by Verhulst et. al. 2012
%   Usage: output = verhulst2012(insig,fs,fc,varargin)
%
%   Input parameters:
%      sign             : the input signal to be processed [time channel]. 
%                         Each channel will be independently processed in parallel.
%      fs               : sampling rate (Hz)
%      fc               : list of frequencies specifying the probe positions on
%                         the basilar membrane, or 'all' to probe all 1000
%                         cochlear sections
%      spl              : array of SPLs that correspond to value 1 of the correspondent input channel               
%       
%   Output parameters:
%       V               : velocity of the basilar membrane sections V(time,section,channel) 
%       Y               : displacement of the basilar membrane sections Y(time,section,channel)
%       OAE             : otto-acoustic emissions given by sound pressure at the middle ear
%       CF              : center frequencies of the probed basiliar membrane sections
%
%   This function computes the basilar membrane displacement and the
%   velocity of the movement at different positions employing a faster
%   implementation of the nonlinear time-domain model of cochlea by
%   Verhulsts, Dau, Shera 2012, through the method described in Altoe et
%   al. 2014
%
%   VERHULST2012 accepts the following optional parameters:
%
%     'normalizeRms',n     array to control the normalization of each channel. With value 1 
%                          normalize the energy of the signal, so the SPL corresponds to the RMS of the signal (default 0)
%
%     'subject',s          the subject number controling the cochlear irregulatiries (default 1)
%
%     'irr',i              array that enable (1) or disable (0) irregularities and nonlinearities for each simulation (default 1) 
%
%
%   The processing is implemented as follows:
%
%   1) the input signal is resampled to the 96 kHz sampling rate employed in the cochlea model
%
%   2) the list of frequencies in fc are converted in to probe 
%      positions in a manner that the frequencies are divided evenly into
%      low and high frequency categories. 
%
%   3) the signals are processed in parallel
%
%   4) the values obtained are resampled back to the original sampling
%      rate
%
%
%   Requirements and installation: 
%   ------------------------------
%
%   1) Python >2.6 is required with numpy and scipy packages. On Linux, use sudo apt-get install python-scipy python-numpy
% 
%   2) Compiled files with a C-compiler, e.g. gcc. In amtbase/src/verhulst start make (Linux) or make.bat (Windows)
%
%   
%
%   References:
%     S. Verhulst, T. Dau, and C. A. Shera. Nonlinear time-domain cochlear
%     model for transient stimulation and human otoacoustic emission. The
%     Journal of the Acoustical Society of America, 132(6):3842 -- 3848,
%     2012.
%     
%     A. Altoè, S. Verhulst, and V. Pulkki. Transmission line cochlear
%     models: improved accuracy and efficiency. J. Acoust. Soc. Am.,
%     136:EL302--EL308, 2014.
%     
%            
%   See also: verhulst2015 verhulst2018 demo_verhulst2012
%             verhulst2018_ihctransduction verhulst2015_cn
%             verhulst2015_ic verhulst2018_auditorynerve exp_verhulst2012
%             verhulst2012 verhulst2015
%             verhulst2018 middleearfilter data_takanen2013 takanen2013_periphery
%             exp_osses2022 exp_takanen2013 takanen2013
%
%   Url: http://amtoolbox.org/amt-1.6.0/doc/models/verhulst2012.php


%   #StatusDoc: Perfect
%   #StatusCode: Perfect
%   #Verification: Verified
%   #Requirements: MATLAB M-Signal PYTHON3 C
%   #Author: Alessandro Altoe' (2014)

% This file is licensed unter the GNU General Public License (GPL) either 
% version 3 of the license, or any later version as published by the Free Software 
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and 
% at <https://www.gnu.org/licenses/gpl-3.0.html>. 
% You can redistribute this file and/or modify it under the terms of the GPLv3. 
% This file is distributed without any warranty; without even the implied warranty 
% of merchantability or fitness for a particular purpose. 

definput.import={'verhulst2012'};  % load defaults from arg_verhulst2012

[flags,keyvals]  = ltfatarghelper([],definput,varargin);


if isempty(keyvals.normalize),
    [channels,idx]=min(size(sign));
    normalizeRMS=zeros(channels,1);
else
  normalizeRMS=keyvals.normalize; 
end

if isempty(keyvals.irr)
    [channels,idx]=min(size(sign));
    irregularities=ones(1,channels);
else
  irregularities = keyvals.irr;
end

subject=keyvals.subject;
modfs=96000;
sectionsNo=1000;
[channels,idx]=min(size(sign));
if(idx==2) %transpose it (python C-style row major order)
    sign=sign';
end
stim=zeros(channels,length(resample(sign(1,:),modfs,fs)));
for i=1:channels
    stim(i,:)=resample(sign(i,:),modfs,fs);
    if normalizeRMS(i)
        s_rms=rms(stim(i,:));
        stim(i,:)=stim(i,:)./s_rms;
    end
end
% sheraPo=0.061;
if(ischar(fc) && strcmp(fc,'all')) %if probing all sections 1001 output (1000 sections plus the middle ear)
    p=sectionsNo;
else %else pass it as a column vector
    [p,idx]=max(size(fc));
    if(idx==2)
        fc=fc'; 
    end
    fc=round(fc);
end
len=length(stim(1,:));

% probes=fc; Fs=modfs; sheraPo=0.061;
% path=fileparts(which('verhulst2012'));
% path=fileparts(path);
% act_path=pwd;
% cd(strcat(path,'/environments/verhulst2012/')); 
% save('input.mat','stim','Fs','channels','spl','subject','sheraPo','irregularities','probes','-v7');
% system('python run_cochlear_model.py');

in.stim=stim; in.Fs=modfs; in.channels=channels; in.spl=spl; 
in.subject=subject; in.sheraPo=0.061; in.irregularities=irregularities;
in.probes=fc;
out.v=[p len channels];
out.y=[p len channels];
out.E=[len 1 channels];
out.F=[p 1];
output=amt_extern('Python','verhulst2012','run_cochlear_model.py',in,out); 
Vs=output.v;
Ys=output.y;
OAEs=squeeze(output.E);
CF=output.F;

% Vs=zeros(p,len,channels);
% Ys=zeros(p,len,channels);
% OAEs=zeros(len,channels);
% for i=1:channels
%     fname=strcat('out/v',int2str(i),'.np');
%     f=fopen(fname,'r');
%     Vs(:,:,i)=fread(f,[p,len],'double','n');
%     fclose(f);
%     fname=strcat('out/y',int2str(i),'.np');
%     f=fopen(fname,'r');
%     Ys(:,:,i)=fread(f,[p,len],'double','n');
%     fclose(f);
%     fname=strcat('out/E',int2str(i),'.np');
%     f=fopen(fname,'r');
%     OAEs(:,i)=fread(f,[len,1],'double','n');
%     fclose(f);
%     if(i==1)
%       fname=strcat('out/F',int2str(i),'.np');
%       f=fopen(fname,'r');
%       CF=fread(f,[p,1],'double','n');
%       fclose(f);
%     end
% end

rl=length(resample(stim(1,:),fs,modfs));
V=zeros(rl,p,channels);
Y=zeros(rl,p,channels);
OAE=zeros(rl,channels);

for i=1:channels
    V(:,:,i)=resample(squeeze(Vs(:,:,i))',fs,modfs);
    Y(:,:,i)=resample(squeeze(Ys(:,:,i))',fs,modfs);
    OAE(:,i)=resample(squeeze(OAEs(:,i)),fs,modfs);
end
% cd(act_path);