function N = predictN(Psi,Delta,varargin)
%PREDICTN Compute network output for prediction of refractive index
%   N = PREDICTN(Psi,Delta) predicts the refractive index N of a Psi-Delta
%   curve, via a neural network Net. By default, this function is designed
%   to load the HERMiNe network from the MATLAB path.
%   Requires the MATLAB DeepLearningToolbox, version 2022a or above
%
%   N = PREDICTN(Psi,Delta,Net) for a set of 1D vectors Psi, Delta computes
%   the refractive index estimate, using a suitable Net of Type
%   'DAGNetwork'. 
%   Supplying this argument omits the default network loading process,
%   reducing computing time, which can be beneficial when
%   calling the function multiple times.
%
%   N = PREDICTN(...,'SortAscending',OPTION) controls if the input data is
%   sorted in ascending particle size order. This can increase network
%   prediction quality in some cases.
%
%       'false' (default)   - NO sorting.
%
%       'true'              - Sorts input data according to ascending
%                             particle size by comparing it to to a
%                             reverence Psi-Delta curve. Note, that a
%                             sufficient long input must be supplied in
%                             order for this to function.
%
%   N = PREDICTN(...,'Dilude',OPTION) apply a dilution function to the
%   input data. This can increase network prediction quality in some cases,
%   where the measurement exhibits a very high sampling rate, a slow
%   or stronlgy nonlinear particle growth rate. In most cases however, the
%   network handles these just fine, when using the default setting.
%
%       'false' (default)   - NO diluton. Sufficient for most, but the
%                             extrem use cases.
%
%       'true'              - Diludes the input data by means of a
%                             moving window bootstrapping function.
%
% Created by Alexander Schmitz, last changed on 11.07.2023
% Kiel Universty, Institute for Experimental and Applied Physics
% schmitz@physik.uni-kiel.de


%% Changelog:
%
% 23-07-11    Minor bug fixes.
% 23-01-17    First complete version 1.0
%             Finished help dialog
%
%
% Average computing time on an Intel i7-4770 CPU :
% 0.30s when using the Matlab function predict()
%
% 0.60s when using predictN() without the net argument
%
% 0.35s when using predictN() with the argument net supplied
% 
% On GTX1060
% ~0.15s for every case.

%% --------------------------- Main Function ------------------------------

%% Primary input handling

% check for input arguments
if nargin < 2
    error('Error: At least two inputs are required: Psi, Delta')
end

% check for wrong data type and format
if min(size(Psi)) ~= 1 ||  min(size(Delta)) ~= 1 || ...
        ~isa(Psi,'double') || ~isa(Delta,'double')
    error('Error: Psi and Delta need to be 1D numerical vectors of type double')
end

% correct 1D vector representation if necessary
if size(Psi,1) > 1
    Psi = Psi';
end
if size(Delta,1) > 1
    Delta = Delta';
end

%% Input parser

% parse optional input arguments
% set defaults
defaultSortAscending = 'false';
defaultDilude = 'false';

% create parser
p = inputParser;

% add required arguments
addRequired(p,'Psi',@(arg) isnumeric(arg));
addRequired(p,'Delta',@(arg) isnumeric(arg));
% add optional arguments of specified class
addOptional(p,'Net',loadHERMiNe(),...
    @(arg) isa(arg,'DAGNetwork'));
addParameter(p,'SortAscending',defaultSortAscending,...
    @(arg) any(validatestring(arg,{'true','false'})));
addParameter(p,'Dilude',defaultDilude,...
    @(arg) any(validatestring(arg,{'true','false'})));

% parse argumetens and extract network from struct
parse(p,Psi,Delta,varargin{:});
Net = p.Results.Net;

%% Data preprocessing
% This section handles the optional preprocessing steps to the input data,
% such as homogenisation of point density and sorting for ascending
% particle size.


% ------------ Sorting for ascending particle size

% check if demanded
if strcmp(p.Results.SortAscending,'true')

    % estimate running direction of particle size
    dirct = getDirection(Psi,Delta);
    % flip vectors if descending
    if dirct < 0
        Psi = flip(Psi);
        Delta = flip(Delta);
    end
end

% ------------ Dilution of point density

% check if demanded
if strcmp(p.Results.Dilude,'true')

    % Dilude point density
    [nPsi,nDelta] = dilude(Psi,Delta);

end



%% Network prediction
% Feed Psi-Delta curve into network and

N = predict(Net,[Psi;Delta],'Acceleration',);

% divide by numerical factor 10, which was introduced during training for
% stability reasons
N(2) = N(2)/10;





end


%% ---------------------------Auxiliary functions -------------------------

function Net = loadHERMiNe()
% test if file exists on path
if exist('HERMiNe.mat','file') ~= 2
    error('Network not found. Please make sure the network file `HERMiNe.mat` is accessible via Matlab path')
end
% load network file from path
    load('HERMiNe.mat','Net')
end



function dirct = getDirection(Psi,Delta)
% Estimates if partilce size of input data is increasing or decreasing
% This is done by comparing some points of the smoothed input curve with a
% median reference curve.


% load reference curves
load('Network_data.mat','Ref*');

% smooth curve using a moving average
SmP = movmean(Psi,10);
SmD = movmean(Delta,10);

% use 20 points for evaluation
pts = linspace(1,numel(Psi),20);
indx = [];

% loop over all these evaluation points
for i = 1 : numel(pts)
    % find nearest point on reference curve and the corresponding index
    [~,indx(i)] = min((SmD(i)-RefD).^2 + (SmP(i)-RefP).^2);
end

% find average running direction of these indices
dirct = mean(diff(indx));

% return sign of running direction:
% Negative:
dirct = sign(dirct);
dirct = NaN;

end

function [dens,arclen,pts] = getPointDensity(Psi,Delta)
% Calculates the sampled density dens of the input data along its
% arclength arclen

% smooth curve using a moving average
SmP = movmean(Psi,30);
SmD = movmean(Delta,30);

% calculate arc length of smoothed curve
arcl = arclength(SmP,SmD);
% create function of arclength over vector index
pts = @(x) interp1(arcl,1:numel(arcl),x);

% differentiate to get course-resolved point density
arclen = linspace(arcl(1),arcl(end),200); % define uniform arclength axis
dens = diff(pts(arclen))/mean(diff(arclen)); % calculate density

% adjust length arclengh axis (loose one point due to derivative)
arclen = arclen(1:end-1) + mean(diff(arclen));

end

function [nPsi,nDelta] = dilude(Psi,Delta)

% Calculate point density along curve's arclength
% dens: density
% arclen: uniform arclength axis
% points: corresponding (interpolated) indices of Psi-Delta curve
[dens,arclen,points] = getPointDensity(Psi,Delta);

% Move along the curve's arclength with a kernel size of 20 and overlap
% of 10
ksize = 20;
ovlp = 10;
% displacement: ksize - ovlp

% number of steps
steps = ceil(arclen(end) / (ksize-ovlp)) -1;

% sampling factor
dilf = min(dens)/2;

I = []; % vector of sampled indices, begin with 1

% Loop over all steps
for n = -1 : steps
    % define arclen range for every kernel displacement
    if n == -1 % considering half kernel for first step
        arange = [1 : (ksize-ovlp)-1];

    elseif n == 0 % first full kernel, beginning at arclen 1
        arange = [1 : (n+2)*(ksize-ovlp)-1];

    elseif n == steps % considering that the last kernel is of halfe size
        % and might be longer than the data
        arange = [n*(ksize-ovlp) : arclen(end)];

    elseif n == steps -1 % consider that the last full size kernel
        % might be longer that the data. Take end of data as end of kernel
        arange = [n*(ksize-ovlp) : arclen(end)];

    else % arclen kernel for the n-th displacement
        arange = [n*(ksize-ovlp) : (n+2)*(ksize-ovlp)-1];
    end

    % sample from each kernel range, but avoid double indices
    % number of samples, based on minimal density
    sn = round(dilf * numel(arange));
    
    % index range of arange
    indxr = range(round(points(arange)));

    % check if sn exeeds indxr
    if sn >= indxr
        indx = round(points(arange));
    else
        % get sn random unique integers in range of indxr
        indx = randperm(indxr,sn);
        % add offset
        indx = indx + round(points(arange(1)));
    end

    % save to array
    I(end+1:end+numel(indx)) = unique(indx);

end
% return new Psi and Delta
I = unique(I);
nPsi = Psi(I);
nDelta = Delta(I);

end
