function [MCresult, FORMresult] = profile_fluctuation_probcalc(varargin)
%PROFILE_FLUCTUATION_PROBCALC  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = profile_fluctuation_probcalc(varargin)
%
%   Input: For <keyword,value> pairs call profile_fluctuation_probcalc() without arguments.
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   profile_fluctuation_probcalc
%
%   See also

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2014 Deltares
%       Kees den Heijer
%
%       Kees.denHeijer@deltares.nl
%
%       Deltares
%       P.O. Box 177
%       2600 MH Delft
%       The Netherlands
%
%   This library is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this library.  If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://www.OpenEarth.eu">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute
% your own tools.

%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
% Created: 10 Mar 2014
% Created with Matlab version: 8.1.0.604 (R2013a)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $

%%
OPT = struct(...
    'SaveDir',          '',             ...
    'jarkus_url',       jarkus_url,     ...
    'id',               7003800,        ...
    'station1',         'Den Helder',   ...
    'station2',         'IJmuiden',     ...
    'year',             2013,           ...
    'ProfileVariation', true,           ...
    'VariationVolumeMean',  0,          ...
    'VariationVolumeSTD',   100,        ...
    'MaxRelError',      0.1,            ...
    'MaxNrIterations',  5,              ...   
    'NrSamplesMC',      1000,           ...
    'ReturnFreqs',      [1e-3 1e-4 1e-5]);

% return defaults (aka introspection)
% if nargin==0;
%     varargout = {OPT};
%     return
% end
% overwrite defaults with user arguments
OPT = setproperty(OPT, varargin);

%% code
for iYear = 1:numel(OPT.year)
    [stochast, Lambda1] = profile_fluctuation_stochast(...
        'jarkus_url', OPT.jarkus_url, ...
        'id',       OPT.id,         ...
        'station1', OPT.station1,   ...
        'station2', OPT.station2,   ...
        'ProfileVariation', OPT.ProfileVariation, ...
        'VariationVolumeMean', OPT.VariationVolumeMean, ...
        'VariationVolumeSTD', OPT.VariationVolumeSTD, ...
        'year',     OPT.year(iYear));
    
    tr = jarkus_transects(...
        'year', OPT.year(iYear),       ...
        'url',  OPT.jarkus_url, ...
        'id',   OPT.id);
    
    nnid    = ~isnan(tr.altitude);
    crid    = squeeze(nnid(1,1,:));
    z       = squeeze(tr.altitude(nnid));
    x       = tr.cross_shore(crid)';
    
    % importance sampling structure
    sampling = struct(                              ...
        'Name', {                                   ...
        'Ph'                                        ...
        },                                          ...
        'Method', {                                 ...
        @prob_is_uniform                            ...
        },                                          ...
        'Params', {                                 ...
        num2cell([0 6])                             ...
        }                                           ...
        );
    
    xref = max(findCrossings(x,z, x, 5*ones(size(x))));
    
    RetreatDistanceRange = 0:300;
    
    MCresult(iYear,:) = {OPT.year(iYear), MC(stochast,...
        'x2zFunction', @profile_fluctuation_LSF,...
        'NrSamples', OPT.NrSamplesMC,...
        'x2zVariables', {'x' x 'z' z 'lambda' Lambda1 'station1' OPT.station1 'station2' OPT.station2 'critical_rsp' xref - RetreatDistanceRange},...
        'IS', sampling)};
    
    save(fullfile(OPT.SaveDir,[num2str(OPT.id) '_MCresult_' num2str(OPT.year(iYear)) '_PV.mat']),'MCresult');
    
    %% First FORM iteration
    idx = NaN(numel(OPT.ReturnFreqs),1);
    for iReturnFreq = 1:numel(OPT.ReturnFreqs)
        [MinVal, ~] = min(abs(MCresult{iYear,2}.Output.P_f - OPT.ReturnFreqs(iReturnFreq)));
        TempIdx = find(abs(MCresult{iYear,2}.Output.P_f - OPT.ReturnFreqs(iReturnFreq)) == MinVal);
        
        if numel(TempIdx) > 1
            % In case there are multiples; select the one that has a
            % neighbour closest to desired Pf
            if sign(MCresult{iYear,2}.Output.P_f(min(TempIdx)-1) - OPT.ReturnFreqs(iReturnFreq)) == sign(MCresult{iYear,2}.Output.P_f(min(TempIdx)) - OPT.ReturnFreqs(iReturnFreq))
                idx(iReturnFreq) = max(TempIdx);
            else
                idx(iReturnFreq) = min(TempIdx);
            end
        else
            idx(iReturnFreq) = TempIdx;
        end
        
        FORMresult{iYear,iReturnFreq,1}     = OPT.year(iYear);
        FORMresult{iYear,iReturnFreq,2}     = OPT.ReturnFreqs(iReturnFreq);
        FORMresult{iYear,iReturnFreq,3}     = FORM(stochast,...
            'x2zFunction', @profile_fluctuation_LSF,...
            'x2zVariables', {'x' x 'z' z 'lambda' Lambda1 'station1' OPT.station1 'station2' OPT.station2 'critical_rsp' MCresult{iYear,2}.settings.x2zVariables{12}(idx(iReturnFreq))});
        RetreatDistances(iReturnFreq) = MCresult{iYear,2}.settings.x2zVariables{12}(idx(iReturnFreq));
        Pfs(iReturnFreq)              = FORMresult{iYear,iReturnFreq,3}.Output.P_f;
    end
    nFORMiteration  = 1;
    
    %% Check whether results are accurate enough, else another iteration
    while any(CheckConvergence(FORMresult, nFORMiteration, OPT, iYear)) && nFORMiteration < OPT.MaxNrIterations
        fitRDs          = FitFORM(Pfs, RetreatDistances, OPT);
        if any(fitRDs > max(xref - RetreatDistanceRange))
            idMax           = find(fitRDs > max(xref - RetreatDistanceRange));
            fitRDs(idMax)   = max(xref - RetreatDistanceRange);
        elseif any(fitRDs < min(xref - RetreatDistanceRange))
            idMin           = find(fitRDs < min(xref - RetreatDistanceRange));
            fitRDs(idMin)   = min(xref - RetreatDistanceRange);
        end
        NotConverged    = CheckConvergence(FORMresult, nFORMiteration, OPT, iYear);
        idx2            = find(NotConverged);
        nFORMiteration  = nFORMiteration + 1;
        for iRD = 1:sum(NotConverged)
            FORMresult{iYear,idx2(iRD),2+nFORMiteration} = FORM(stochast,...
                'startU',   FORMresult{iYear,idx2(iRD),1+nFORMiteration}.Output.designpoint.finalU, ...
                'x2zFunction', @profile_fluctuation_LSF,...
                'x2zVariables', {'x' x 'z' z 'lambda' Lambda1 'station1' OPT.station1 'station2' OPT.station2 'critical_rsp' fitRDs(idx2(iRD))});
            RetreatDistances(numel(RetreatDistances)+1) = fitRDs(idx2(iRD));
            Pfs(numel(Pfs)+1)                           = FORMresult{iYear,idx2(iRD),2+nFORMiteration}.Output.P_f;
            display(['** ' num2str(OPT.year(iYear)) ' RD=' num2str(fitRDs(idx2(iRD))) ' FORM iteration ' num2str(nFORMiteration) ' for return freq ' num2str(OPT.ReturnFreqs(idx2(iRD))) ' completed **'])
        end
    end
    
    save(fullfile(OPT.SaveDir,[num2str(OPT.id) '_FORMresult_' num2str(OPT.year(iYear)) '_PV.mat']),'FORMresult');
    
    clear Pfs RetreatDistances
end
end

function fitRDs = FitFORM(Pfs, RDs, OPT)
fitRDs  = NaN(1,numel(OPT.ReturnFreqs));
[~, iUa, iUc]   = unique(RDs);
RDs             = RDs(iUa);
Pfs             = Pfs(iUa);
% PolyFit         = polyfitn(log10(Pfs),RDs,1);
% fitRDs          = polyvaln(PolyFit,log10(OPT.ReturnFreqs));
for iFreq = 1:numel(OPT.ReturnFreqs)
    idx             = isort(abs(log10(Pfs)-log10(OPT.ReturnFreqs(iFreq))));
    PolyFit         = polyfitn(log10(Pfs(idx(1:min(numel(idx),3)))),RDs(idx(1:min(numel(idx),3))),1);
    fitRDs(1,iFreq) = polyvaln(PolyFit,log10(OPT.ReturnFreqs(iFreq)));
end
end

function anotherIteration = CheckConvergence(FORMresult, nFORMiteration, OPT, iYear)
for iFreq = 1:numel(OPT.ReturnFreqs)
    if ~isempty(FORMresult{iYear,iFreq,2+nFORMiteration})
        anotherIteration(iFreq) = abs((FORMresult{iYear,iFreq,2+nFORMiteration}.Output.P_f - OPT.ReturnFreqs(iFreq))/OPT.ReturnFreqs(iFreq)) > OPT.MaxRelError;
    else
        anotherIteration(iFreq) = false;
    end
end
end