function [outX, outY] = extractDataFromPbeResult(taskString, taskParameter, ...
    resultStruct, pbeDef)
% extract information from solver result structure
%
% Syntax: [outX outY] = extractDataFromPbeResult(taskString, propertyHandle, ...
%                         resultStruct, pbeDefinition)
%
% This function generates vectors or matrices that can be used to plot data
%   obtained by solvePbe. Therefore, the taskString specifies the type of
%   output where the task parameter selects a state variable, a property or
%   gives a custom function handle.
%
% The following tasks return a time vector and a value vector of some kind:
%   bulkProperty - for taskParameter being an integer, the bulk property
%       defined in pbeDefinition is used. The taskParameter can also be a
%       function handle used for bulk properties in solvePbe @(P, N, x, t).
%   bulkState - the taskParameter is interpretet as an index to the state
%       vector.
%   integralProperty - taskParameter is the string 'number' or a string of
%     a measure contained in pbeDef.dec.measureData. The integral of that
%     measure (first row of Y) ingluding its standard deviation (second row
%     of Y) is calculated.
%   genericIntegral - the taskParameter is interpreted as a function handle
%     @(H, up), output is given as above.
%   propertyDistribution - task parameter is a cell array where {1} gives
%     the desired measure used on the y-axis of the distribution and {2}
%     gives the measure on the x-axis. {3} contains a matrix where the
%     first row represents the lower bound of each bin and the second row
%     the upper bound. Whether these bins overlap, are continuous between
%     two limits or where the pivot point is is not considered.
%
% See also: solvePbe, executePbeCase, loadPbeCase

% 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/>.

%% catch empty result
if isempty(resultStruct)
    error('elk:pbeSolver:wrongInput', 'PBE result structure is empty')
end

%% for measure computation
% mappingToMeasComputation = eye(pbeDef.dec.nC);
% if pbeDef.dec.isProper || ~pbeDef.optionStruct.useEmbeddingPartitions
%     dec = pbeDef.dec;
% elseif pbeDef.dec.isComplete
%     dec = pbeDef.dec.properData;
%     mappingToMeasComputation = pbeDef.dec.properData.mappingNewToProper;
% else
%     % incomplete dec
%     dec = pbeDef.dec;
% end
dec = pbeDef.dec;

switch lower(taskString)    
    %% bulk property
    case 'bulkproperty'
        if isa(taskParameter, 'function_handle')
            % given function handle evaluates bulk property
            thisFunHandle = @(info)(taskParameter(info.P, info.N, ...
                info.x, info.time));
            
        elseif (taskParameter > 0) && (mod(taskParameter, 1) == 0)
            % taskString indicates existing property
            thisFunHandle = @(info)(subsref(...
                pbeDef.bulkProperty(info.P, info.N, info.x, info.time), ...
                struct('type', '()', 'subs', {{taskParameter}})));
          
        else
            error('elk:pbeSolver:wrongInput', ['task parameter ' ...
                'must be a function handle or a positive integer']);
        end
        [outX, outY] = obtainTimeEvolution(...
            resultStruct, thisFunHandle);
        
    %% bulk state
    case 'bulkstate'
        if (taskParameter > 0) && (mod(taskParameter, 1) == 0)
            % taskString indicates existing state
            thisFunHandle = @(info)(info.x(taskParameter));
        else
            error('elk:pbeSolver:wrongInput', ['task parameter ' ...
                'must be a positive integer']);
        end
        [outX, outY] = obtainTimeEvolution(...
            resultStruct, thisFunHandle);
        
    %% master integral evaluation
    case 'masterintegral'
        % integral: note that the core integral is already evaluated as
        %   extra information so that it is easily reused for the error
        %   estimates. The real integral is, therefore, given in the last
        %   line of this case.
        thisFunHandle = @(info)([info.extra(end).data; ...
            sqrt(info.nSample/(info.nSample-1) * (...
                sum((...
                  info.extra(end-1).data .* ...
                  (ones(size(info.extra(end-1).data, 1), 1)*info.N) - ...
                  (info.extra(end).data*...
                   ones(1, size(info.extra(end-1).data, 2)))/info.nSample ...
                  ).^2, 2) + ...
                  info.nZero/info.nSample.^2*info.extra(end).data.^2 ...
                ))]);
            
        % allow for multiple parameter input
        if ~iscell(taskParameter)
            taskParameter = {taskParameter};
        end
        % call the subroutine that crawls through the time steps in
        %   resultStruct. This line also includes the evaluation of the
        %   integral
        [outX, outY] = obtainTimeEvolution(...
            resultStruct, thisFunHandle, taskParameter{:}, ...
            @(info)(info.extra(end).data*info.N'));
    
    %% integral property (generic)
    case 'genericintegralproperty'
        % integrand
        integrandHandle = @(info)(taskParameter(info.H, info.up));
        % evaluate
        [outX, outY] = extractDataFromPbeResult('masterIntegral', ...
            integrandHandle, resultStruct, pbeDef);
    
    %% integal property (string, index)
    case 'integralproperty'
        % integrand
        if strcmpi(taskParameter, 'number')
            integrandHandle = @(info)(ones(size(info.N)));
        elseif ischar(taskParameter)
            integrandHandle = @(info)(computeMeasureDec(dec, 'hC', ...
                {info.up, info.H}, taskParameter, [], 0, 1));
        elseif iscell(taskParameter)
            integrandHandle = @(info)(computeMeasureDec(dec, 'hC', ...
                {info.up, info.H}, {taskParameter{1}, taskParameter{2}}, [], 0, 1));
        elseif isnumeric(taskParameter) && taskParameter == 0
            integrandHandle = @(info)(info.H);
        elseif isnumeric(taskParameter) && taskParameter < 0
            disp('yup')
            integrandHandle = @(info)(info.H(-1*taskParameter, :));
        else
            integrandHandle = @(info)(info.P(taskParameter, :));
        end
        % evaluate
        [outX, outY] = extractDataFromPbeResult('masterIntegral', ...
            integrandHandle, resultStruct, pbeDef);
        
    %% property distribution
    case 'propertydistribution'
        % yAxis property
        if ~ischar(taskParameter{1})
            yAxisHandle = @(info)(info.P(taskParameter{1}, :));
        elseif strcmpi(taskParameter{1}, 'number')
            yAxisHandle = @(info)(ones(size(info.N)));
        else
            yAxisHandle = @(info)(computeMeasureDec(dec, 'hC', ...
                {info.up, info.H(1:pbeDef.nDecCoordinate,:)}, taskParameter{1}, [], 0, 1));
        end
        % xAxis property
        if ~ischar(taskParameter{2})
            xAxisHandle = @(info)(info.P(taskParameter{2}, :));
        elseif strcmpi(taskParameter{2}, 'number')
            xAxisHandle = @(info)(ones(size(info.N)));
        else
            xAxisHandle = @(info)(computeMeasureDec(dec, 'hC', ...
                {info.up, info.H(1:pbeDef.nDecCoordinate,:)}, taskParameter{2}, [], 0, 1));
        end
        % integrand: time evolution will evaluate yAxisProperty and
        %   xAxisProperty in info.extra(1).data and info.extra(2) so that
        %   the integrand can be constructed on these information. The
        %   integrand looks like:
        %     yProperty*(xProperty > lowerBound)*(xProperty <= upperBound)
        %   Now consider that yProperty is the vector as given by
        %   yAxisHandle. While yProperty is a row vector, the boundary
        %   filters can be written as matrices that are multiplied from the
        %   right:
        %     yProperty*boundaryFilterMatrix
        %   so that the result of the integral will become a column vector.
        % Hence, we construct the boundary filter, first.
        boundaryFilter = @(info)(...
            bsxfun(@gt, info.extra(2).data, taskParameter{3}(1,:)') .* ...
            bsxfun(@le, info.extra(2).data, taskParameter{3}(2,:)'));
        % integrand      
        integrandHandle = @(info)(info.extra(3).data .* ...
            (ones(size(info.extra(3).data, 1), 1)*info.extra(1).data));
        % call master integral
        [outX, outY] = extractDataFromPbeResult('masterIntegral', ...
            {yAxisHandle, xAxisHandle, boundaryFilter, integrandHandle}, ...
            resultStruct, pbeDef);
        % make density values
        size(outY)
        size(taskParameter{3}(2,:))
        outY = bsxfun(@rdivide, outY, ...
            [(taskParameter{3}(2,:) - taskParameter{3}(1,:))';
             (taskParameter{3}(2,:) - taskParameter{3}(1,:))']);
        
              
    otherwise
        error('elk:populationBalance:wrongInput', 'Specified task is unknown.')
end

end