function [result pbeDefinition] = solvePbe(pbeDefinition, varargin)
% solve a population balance (growth, nucleation, seeding, bulk states)
%
% Syntax: result = solvePopulationBalance(pbeDefinition)
%
% This routine essentially implements the method of characteristics for a
%   set of input points. In this scope, the routine considers the
%   limitation of the growth rate to the gemetrical boundaries (e.g.
%   validity). It does also recognize dissolving particles.
% Not explicitely contained in the scheme of this function is the
%   integration of the particle population to, e.g. the total particle
%   volume. This functionality is implicitely provided by the definition of
%   bulk properties in the PBE definition structure (pbeDef). However, this
%   library follows the approach of Monte Carlo integral estimates which
%   is, in particular, reflected by the function extractDataFromPbeResult.
% The support for multiple seeding and nucleation was cancelled because the
%   overal scheme with Monte Carlo integral estimates does not yet cover
%   these situations. It would, however, be straightforward to add this
%   functionality. Your turn! ;)
% For input/output see the link below.
%
% Detailed documentation for inpu/output structures and option/value pairs
%   for the function call are available, seperately, <a href="matlab: helpPbeSolver">here</a>.
%   This documentation contains information on:
%     <a href="matlab: help helpPbeSolver>optionValuePairs">List of option/value pairs</a>
%     <a href="matlab: help helpPbeSolver>variables">Nomenclature  for input/output structures</a>
%     <a href="matlab: help helpPbeSolver>inputStruct">Input structure: pbeDefinition</a>
%       > <a href="matlab: help helpPbeSolver>seedingStruct">Structure for a seeding event.</a>
%     <a href="matlab: help helpPbeSolver>outputStruct">Output structure array: result</a>
%
% See also: extractDataFromPbeResult, executePbeCase, loadPbeREsult,
%   limitGrowthRate

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program 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 program 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 program.  If not, see <http://www.gnu.org/licenses/>.

% start timer
tic

%% Input handling
% variable input stuff
optionStruct = mapOptionStruct(varargin, ...
    'odeSolver', @ode45, ...
    'absTolPosition', 1e-12, ...
    'absTolNumber', 1e10, ...
    'absTolBulk', 1e-12, ...
    'absTolDensity', 1e8, ...
    'relTol', 1e-5, ...
    'initialTimeStep', 1, ...
    'distanceLimit', 1e-8, ...
    'binLoad', 1e17, ...
    'gradientRelativeDelta', 1e-3, ...
    'gradientAbsoluteDelta', 1e-10, ...
    'updateUnifiedLimit', 1, ...
    'updateTiming', 1, ...
    'limitGrowthRate', 1, ...
    'useEmbeddingPartitions', 1, ...
    'verbose', 2, ...
    'density', 0);

% allow pbeDef to overwrite some of the options above
if isfield(pbeDefinition, 'optionStruct')
    pbeDefinition.optionStruct = mergeStruct(pbeDefinition.optionStruct, ...
        optionStruct);
else
    pbeDefinition.optionStruct = optionStruct;
end
optionStruct = pbeDefinition.optionStruct;

% prepare additional information
pbeDefinition = checkInputSolver(pbeDefinition);
% pbeDefinition = computeLimitCone(pbeDefinition);
[pbeDefinition.seeding seedingTimeVector] = prepareSeeding(...
    pbeDefinition.seeding, pbeDefinition, optionStruct);

%% Init: event handling
nextSeedIndex = 1;  % index to next seeding population
eventIndex = [];    % index of an event to handle
% set seeding time vector: last event is simulation end
seedingTimeVector = [seedingTimeVector pbeDefinition.timeSpan(end)];
%% Init: solver state
% use a central definition of al properly sorted fields
solverState = obtainCleanSolverState;
% time
solverState.time = pbeDefinition.timeSpan(1);
solverState.evalTime = toc;
% state variables
solverState.bulkStateVector = pbeDefinition.bulkInitial(:);
solverState.hcMatrix = zeros(pbeDefinition.nDecCoordinate + ...
                             pbeDefinition.nExtraCoordinate, 0);
solverState.numberVector = zeros(1, 0);
if optionStruct.density
    solverState.probabilityVector = zeros(1, 0);
end
% help variables
solverState.indexVector = zeros(1,0);
solverState.upVector = zeros(1, 0);
solverState.dissolvedPivotVector = [];
% property information
solverState.pivotPropertyMatrix = [];
solverState.bulkPropertyVector = [];
% option/structure information
solverState.flagNucleation = 0;
solverState.flagDensity = optionStruct.density;

%% Init: nucleation state
% if nucleation should be swiched on or off corresponds to the third event
%   of the genericEventFunction
[initialEventVector, ~, ~] = genericEventFunction(packState(solverState), ...
    solverState.time, pbeDefinition, solverState);
if initialEventVector(1) > 0
    solverState = addNucleation(solverState, pbeDefinition);
elseif initialEventVector(1) == 0
    error('elk:populationBalance:initialNucleation', ...
        ['The initial nucleation rate is exactly zero. ' ...
         'Please alter the initial nucleation or time slightly to resolve ' ...
         'this issue. If you don''t have nucleation, choose a negative ' ...
         'nucleation rate']);
end

%% Init: result
% initialize field names
result = solverState;
result.hcMatrixDynamic = result.hcMatrix;
result.derivativeStruct = struct([]);
% ensure no data is written
result(:) = [];
% extend to target result vector
result(length(pbeDefinition.timeSpan)).time = [];
% fill time with nan for later filtering
[result(:).time] = deal(nan);

%% solver loop: BEGIN
% 1) Add seeding pivots (if required)
% 2) Setup ODE-state vector:
%    * initial condition
%    * handle for right-hand side
% 3) Integrate to: (1) next seeding event, 
%                  (2) first nucleation bin full // start nucleation
%                  (3) dissolution occures (unknown in which population)
%                  (4) update unified partition data
%                  (5) end of simulation
% 4) Save solver states

if optionStruct.verbose,disp('Start solving pbe ..');end
while solverState.time < seedingTimeVector(end)
       
    %% solver loop:event handling
    eventIndex = unique(eventIndex);
    % correct event indexing for dissolved pivot Indexing
    %                original                
    %                (in event function)     
    % nuc on/off     1                       
    % nuc new bin    2                       
    % update up      3   
    % confined error 4
    % dissolution    5..                     
    particleDissolveIndex = eventIndex(eventIndex > 4) - 4;
    while ~isempty(eventIndex)
        
        %% solver loop:switch nucleation on/off
        if eventIndex(1) == 1
            if solverState.flagNucleation
                disp('  ..stopping nucleation');
                solverState = removeNucleation(solverState);
                % new bin event would be superflous now:
                eventIndex(eventIndex == 2) = [];
                % if simultaneaosly pivots grow out then there new
                %   corresponding index is (boundary removed):
                particleDissolveIndex = particleDissolveIndex - 1;
            else
                disp('  ..starting of nucleation');
                solverState = addNucleation(solverState, pbeDefinition, ...
                    optionStruct);
                % if simultaneaosly pivots grow out then there new
                %   corresponding event index is (bin + boundary added):
                particleDissolveIndex = particleDissolveIndex + 2;
            end
            
        %% solver loop:nucleation new bin event
        elseif eventIndex(1) == 2
            disp(['  ..new nucleation bin, now: ' num2str(sum(...
                solverState.indexVector==-1)+2)]);
            solverState = newBinForNucleation(solverState, pbeDefinition, ...
                optionStruct);
            % if simultaneaosly pivots grow out then there new
            %   corresponding event index is (pivot added):
            particleDissolveIndex = particleDissolveIndex + 1;
            
        %% update unified partitions
        elseif eventIndex(1) == 3
            disp('  ..updating unified partition data') 
            [dump solverState.upVector] = ensureProperDecData(...
                solverState.hcMatrix, solverState.upVector, pbeDefinition);
            
        %% confined vector created
        elseif eventIndex(1) == 4
            %! DEV: Yes, this part is not implemented because it was not
            %       necessary for case studies. See the
            %       thesis/documentation for details. The idea is to handle
            %       confined vectors by the closest unconfined partition.
            %       This is not as trivial as it might sound.
            error('elk:pbeSolver:confinedVectorGenerated', ['This ' ...
                'simulation generated an unconfined vector. This can ' ...
                'happen but you might try an implicit solver or different ' ...
                'initial conditions to overcome this issue. An automatic ' ...
                'handling of this case is not yet implemented. Sorry.']);
            %/
            
            
        %% solver loop:dissolution event
        % nucleatio pivot and boundary are treated otherwise, see the end
        %   of genericEventFunction for more information
        elseif eventIndex > 4 && ...
               particleDissolveIndex >= (1+2*solverState.flagNucleation) && ...
               particleDissolveIndex <= length(solverState.numberVector)
            % remove pivot position
            disp(['  ..removing pivot' num2str(particleDissolveIndex(1)) ...
                  ' .. remaining ' ...
                  num2str(length(solverState.numberVector)-1)]);
            % increase count for removed pivots
            thisPivotClass = solverState.indexVector(particleDissolveIndex(1));
            if thisPivotClass > 0
                solverState.dissolvedPivotVector(thisPivotClass) = ...
                      solverState.dissolvedPivotVector(thisPivotClass) + 1;
            end
            solverState.indexVector(particleDissolveIndex(1)) = [];
            solverState.hcMatrix(:, particleDissolveIndex(1)) = [];
            solverState.upVector(particleDissolveIndex(1)) = [];
            % remove pivot numbers
            solverState.numberVector(particleDissolveIndex(1)) = [];
            if solverState.flagDensity
                solverState.probabilityVector(particleDissolveIndex(1)) = [];
            end
            % we start with the first pivot so that ALL remaining pivots
            %   are indexed now as:
            particleDissolveIndex = particleDissolveIndex - 1;
            % remove particle index that was lined up:
            particleDissolveIndex(1) = [];
        %% no valid event
        else
            error('elk:solver:unknownEvent', 'This should not happen.');
        end
        %% remove event
        eventIndex(1) = [];
    end

    %% solver loop: seeding event
    if solverState.time >= seedingTimeVector(nextSeedIndex)
        % shortcuts
        thisHcMatrix = pbeDefinition.seeding(nextSeedIndex).positionMatrix;
        thisNumberVector = pbeDefinition.seeding(nextSeedIndex).numberVector;
        nDecCoordinate = pbeDefinition.nDecCoordinate;
        
        % filter useful points (having some particle number)
        newPivotFilter = thisNumberVector ~= 0;
        nNew = sum(newPivotFilter);
        nVoid = sum(~newPivotFilter);
        
        % assign hcMatrix
        solverState.hcMatrix(:, end + (1:nNew)) = thisHcMatrix(:, ...
            newPivotFilter);
        % assign particle number
        solverState.numberVector(1, end + (1:nNew)) = ...
            pbeDefinition.seeding(nextSeedIndex).numberVector(newPivotFilter);
        
        % assign upVecctor
        [hcMatrixDump solverState.upVector(1, end + (1:nNew)) ~] = ensureProperDecData(...
             thisHcMatrix(:, newPivotFilter), ones(1, nNew), pbeDefinition);
        
        %! Dev: below is obsolete after adding the above
%         if pbeDefinition.dec.isComplete
%             pivotPositionMeasureMatrix = mapToConfinementCone(...
%                 thisHcMatrix(1:nDecCoordinate, newPivotFilter), ...
%                 pbeDefinition.dec, 0);
%             solverState.upVector(1, end + (1:nNew)) = computeUpVector(...
%                 pbeDefinition, pivotPositionMeasureMatrix, ...
%                 elkZerotol);
%         else
%             pivotPositionMeasureMatrix = ...
%                 thisHcMatrix(1:nDecCoordinate, newPivotFilter);
%             solverState.upVector(1, end + (1:nNew)) = 1;
%         end
        %/
        

        solverState.indexVector(1, end + (1:nNew)) = ...
            ones(1, nNew)*nextSeedIndex;
        
        solverState.dissolvedPivotVector(end+1) = nVoid;
        
        if solverState.flagDensity
            solverState.probabilityVector(1, end + (1:nNew)) = ...
                pbeDefinition.seeding(nextSeedIndex).probabilityVector(...
                  newPivotFilter);
        end
        
        nextSeedIndex = nextSeedIndex + 1;
    end
    
%     disp('DBG: ...')
%     solverState
%     derivativeState = genericDerivative(solverState, 0, pbeDefinition)
%     disp('..DBG/')
    
    %% solve a step
    thisStateVector = packState(solverState);
    % break if state vector is empty
    if isempty(thisStateVector)
        warning('elk:solver:emptyState', ...
            'Solver cancelled calculations because of empty state vector');
        return
    end
    
    thisDerivativeFunction = @(time, state)(genericDerivativeStateVector(...
        state, time, pbeDefinition, solverState));

    thisEventFunction = @(time, state)(genericEventFunction(state, time, ...
        pbeDefinition, solverState));
    
    thisTimeSpan = [solverState.time seedingTimeVector(nextSeedIndex)];
    % limit maximum step to updateTiming from total interval
    thisTimeSpan(2) = min(thisTimeSpan(2), thisTimeSpan(1) + ...
        optionStruct.updateTiming*...
        (seedingTimeVector(end) - seedingTimeVector(1)));
    thisTimeSpan = [thisTimeSpan(1), ...
                    pbeDefinition.timeSpan((pbeDefinition.timeSpan > thisTimeSpan(1)) & ...
                                           (pbeDefinition.timeSpan < thisTimeSpan(2))), ...
                    thisTimeSpan(2)];
    
    thisAbsTol = obtainAbsTolVector(solverState, ...
        optionStruct.absTolPosition, ...
        optionStruct.absTolNumber, ...
        optionStruct.absTolBulk, ...
        optionStruct.absTolDensity);
    
    
    [thisTimeVector, thisStateMatrix, eventTime, eventState, eventIndex] = ...
        feval(optionStruct.odeSolver,thisDerivativeFunction, thisTimeSpan, thisStateVector, ...
          odeset(...
            'absTol', thisAbsTol, ...
            'relTol', optionStruct.relTol, ...
            'events', thisEventFunction, ...
            'InitialStep', optionStruct.initialTimeStep ...
          ));
    if ~any([1 eventIndex] == 3)
        eventIndex = [eventIndex 3];
    end

    solverState = unpackState(thisStateMatrix(end, :)', thisTimeVector(end), ...
        solverState);
    if isempty(eventIndex)
        disp(['..stopped ode45 at time ' num2str(solverState.time)])
    else
        disp(['..stopped ode45 for event at time ' num2str(solverState.time) ...
            ' with event index ' num2str(eventIndex(:)')])
    end
    
    %% save unpacked states
    for iTime = 1:length(thisTimeVector)
        % hide event times where the solver also stops
        thisWriteIndex = find(thisTimeVector(iTime) == pbeDefinition.timeSpan, 1);
        if isempty(thisWriteIndex)
            % this is no time to be written
            continue
        end
        
        % initialize thisResult
        thisResult = unpackState(thisStateMatrix(iTime, :)', thisTimeVector(iTime), ...
            solverState);
        % save state values of hc-matrix
        thisResult.hcMatrixDynamic = thisResult.hcMatrix;
        % try to use previous unified index vector for speedup
        if thisWriteIndex > 1 && ~isnan(result(thisWriteIndex - 1).time)
            thisResult.upVector = result(thisWriteIndex - 1).upVector;
        end
        % compute derivative before it is too late (hcMatrix changed below)
        thisResult.derivativeStruct = genericDerivative(thisResult, ...
            thisTimeVector(iTime), pbeDefinition);

        % ensure proper data
        [thisResult.hcMatrix thisResult.upVector] = ...
            ensureProperDecData(thisResult.hcMatrix, thisResult.upVector, ...
            pbeDefinition);
        
        % pivot properties
        thisResult.pivotPropertyMatrix = pbeDefinition.pivotProperty(...
            thisResult.hcMatrix, thisResult.upVector);
        
        % bulk properties
        thisResult.bulkPropertyVector = pbeDefinition.bulkProperty(...
            thisResult.pivotPropertyMatrix, ...
            thisResult.numberVector, thisResult.bulkStateVector, ...
            thisResult.time);
        
        % assign
        thisResult.evalTime = toc;
        result(thisWriteIndex) = thisResult;
    end
    
end
disp('integration complete')

%% cleanup results
removeFilter = isnan([result(:).time]);
if any(removeFilter)
    warning('elk:pbeSolver:missingOutput', ['Not all requested output '...
        'times were generated. Check [result(:).time] for ' ...
        num2str(sum(removeFilter)) ' missing entries.'])
end
result(removeFilter) = [];
