function [xphase,phaseawaremean,phaseawareL,phaseawareU,werror] = create_phaseaware_forecast(x,varargin)

% Creates a phase-aware ensemble forecast from an original ensemble
% forecast.

%depedencies and where to find download them:
%    rednoise:  https://github.com/grinsted/wavelet-coherence/tree/master/private
%    invcwt:    https://www.mathworks.com/matlabcentral/fileexchange/20821-continuous-wavelet-transform-and-inverse
%    contwt:    https://www.mathworks.com/matlabcentral/fileexchange/20821-continuous-wavelet-transform-and-inverse
%    parseArgs: https://github.com/grinsted/wavelet-coherence/tree/master/private
%    formatts:  https://github.com/grinsted/wavelet-coherence
%    circ_mean: https://www.mathworks.com/matlabcentral/fileexchange/10676-circular-statistics-toolbox--directional-statistics-



%Obligatory inputs:
%      x = An array of ensemble members. Rows are points in time and
%      coloums are the individual ensemble members.


[~,dt] = formatts(x(:,1));

%Lorg is the length of the ensemble members before padding.
[Lorg, Nmember] = size(x);

%----------default arguments for the wavelet transform-----------
Args=struct('Pad',0,...      % pad the time series with zeroes (recommended)
    'Dj',1/12, ...    % this will do 12 sub-octaves per octave
    'S0',2*dt,...    % this says start at a scale of 2 hours, days, months, years, etc.
    'J1',[],...
    'Mother','Morlet', ...
    'MaxScale',[],...   %a more simple way to specify J1
    'Morletwavenumber',6,...
    'Paulorder',4,...
    'Dogderivativeorder',2,...
    'Padlength',200,...
    'sigma_pw',1,... %number of standard deviations used to compute phase-aware spread
    'AR1','auto');

Args = parseArgs(varargin,Args);

if isempty(Args.J1)
    if isempty(Args.MaxScale)
        Args.MaxScale=(Lorg*.17)*2*dt; %automaxscale
    end
    % J1+1 is the number of wavelet scales.
    Args.J1=round(log2(Args.MaxScale/Args.S0)/Args.Dj);
end


% set mother wavelet parameter based on chosen mother wavelet
if strcmpi('Morlet',Args.Mother)
    param = Args.Morletwavenumber;
elseif strcmpi('Paul',Args.Mother)
    param = Args.Paulorder;
elseif strcmpi('Dog',Args.Mother)
    param = Args.Dogderivativeorder;
    warning('The DOG mother wavelet is not recomended for computing phasee-aware statistics')
else
    error('Please select among the Morlet, Paul, and Dog mother wavelets');
end

%--------------------------------------------------------------------------
%pad ends of ensemble members with a realization of a white noise process.
%This padding method will reduce edge effects.

xnoise = rednoise(Args.Padlength,0,1);

%reduces spectral leakage from padded region of wavelet spectrum.
xnoise = xnoise/(30*std(xnoise));

xnoise = (ones(Nmember,1)*xnoise')';
x = vertcat(xnoise,x,xnoise);

%-------------------------------------------------------------------------
%remove ensemble members with NaNs.
[~, iy] = find(isnan(x));
x(:,iy) = [];

%Lpad is the length of the ensemble members after padding. Also, the number
%of ensemble members will now be Nmember because  ensemble members with NaNs
%have been removed.
[Lpad, Nmember_noNAN] = size(x);

%-------------------------------------------------------------------------
%compute wavelet and inverse wavelet transformations
if(Nmember_noNAN>0)
    %preallocate memory
    wave = NaN(Args.J1+1,Lpad,Nmember_noNAN);
    modulus = NaN(Args.J1+1,Lpad,Nmember_noNAN);
    phase = NaN(Args.J1+1,Lpad,Nmember_noNAN);
    werror = NaN(Lpad,Nmember_noNAN);
    xnew = NaN(Lpad,Nmember_noNAN);
    
    %compute wavelet and inverse wavelet transformations of each ensemble
    %member
    for jj = 1:Nmember_noNAN
        [wave(:,:,jj),~,SCALE,~,~, PARAMOUT, K] = contwt(x(:,jj),dt,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother,param);
        xnew(:,jj) = invcwt(wave(:,:,jj), Args.Mother, SCALE, PARAMOUT, K);
        modulus(:,:,jj) =  abs(wave(:,:,jj));
        phase(:,:,jj) = angle(wave(:,:,jj));
        %error resulting from the inverse wavelet transformation.
        werror(:,jj) = x(:,jj) - xnew(:,jj);
        
    end
    
    %calculate mean error resulting from the inverse wavelet transformation.
    meanerror = median(werror,2);
    
    %-------------------------------------------------------------------------
    %compute the phase-aware mean
    
    %ensemble mean modulus spectrum
    meanmodulus = mean(modulus,3);
    %ensemble mean phase spectrum
    meanphase = circ_mean(phase,[],3);
    
    %calculate phaseaware ensemble mean
    psi = complex(cos(meanphase),sin(meanphase));
    Z = meanmodulus.*psi;
    phaseawaremean = invcwt(Z, Args.Mother, SCALE, PARAMOUT, K);
    phaseawaremean = phaseawaremean + meanerror';
    %------------------------------------------------------------------------
    %calculate phaseaware spread
    
    modulusdeviation = std(modulus,[],3);
    %upper bound of phaseaware spread envelope
    Umodulus = meanmodulus + Args.sigma_pw*modulusdeviation;
    %lower bound of phaseaware spread envelope.
    Lmodulus = meanmodulus - Args.sigma_pw*modulusdeviation;
    
    Uspectral = Umodulus.*psi;
    Lspectral = Lmodulus.*psi;
    
    %transform the phase-aware spread envelope in spectral space to the
    %phase-aware envelope in physical space.
    phaseawareU = invcwt(Uspectral, Args.Mother, SCALE, PARAMOUT, K);
    phaseawareU = phaseawareU + meanerror';
    
    phaseawareL = invcwt(Lspectral, Args.Mother, SCALE, PARAMOUT, K);
    phaseawareL = phaseawareL + meanerror';
    %----------------------------------------------------------------------------
    % create phaseaware extended forecast by pairing the phase and modulus
    % spectra of different ensemble members.
    
    %num_newmember is the number of new ensemble members that will be
    %created from the phaseaware extenson procedure.
    num_newmembers = Nmember_noNAN^2 - Nmember_noNAN;
    xphase = NaN(Lpad,num_newmembers);
    
    
    % generate combination of vectors to avoid nested for loops
    Z = combvec(1:Nmember_noNAN,1:Nmember_noNAN)';
    index1 = Z(:,1);
    index2 = Z(:,2);
    
    indexdiff = index1 - index2;
    index1(indexdiff==0) = [];
    index2(indexdiff==0) = [];
    
    
    for j = 1:num_newmembers
        psi = complex(cos(phase(:,:,index1(j))),sin(phase(:,:,index1(j))));
        Ztotal = modulus(:,:,index2(j)).*psi;
        xphase_base = invcwt(Ztotal, Args.Mother, SCALE, PARAMOUT, K);
        xphase(:,j) = xphase_base' + meanerror;
    end
    
    %remove padded parts of time series.
    xphase(Lorg+201:Lorg+400,:) = [];
    xphase(1:200,:) = [];
    
    x(Lorg+201:Lorg+400,:) = [];
    x(1:200,:) = [];
    
    phaseawaremean(Lorg+201:Lorg+400) = [];
    phaseawaremean(1:200) = [];
    
    phaseawareL(Lorg+201:Lorg+400) = [];
    phaseawareL(1:200) = [];
    
    phaseawareU(Lorg+201:Lorg+400) = [];
    phaseawareU(1:200) = [];
    
    %Note that no ensemble members will have NaNs in the phase-aware forecast.
    xphase = cat(2,xphase,x);
    
else
    %all ensemble members have NaNs. Can not calculate phaseaware
    %statistics
    xphase = NaN(Lorg,Nmember);
    werror = NaN(1,Lorg);
    phaseawaremean = NaN(1,Lorg);
    phaseawareL = NaN(1,Lorg);
    phaseawareU = NaN(1,Lorg);
    
end

end