function plotHandles = visualizeSplittingFunction(obj,numberOfLocalModels,inputDimensions,axesHandle)
% This function visualizes the splitting functions without normalization.
%
% visualizeSplittingFunction(obj,numberOfLocalModels,inputDimensions,axesHandle)
%
% INPUTS
%
% obj:                      local model network object
%
% numberOfLocalModels:      Integer value equals the number of splitting
%                           functions, that will be visualized.
%
% inputDimensions:          Integer value or [1 x 2] Vector, containing the
%                           number of the input dimensions, that should be
%                           visualized.
%
% axesHandle:               Axes handle to determine in which axes the
%                           visualization should be plotted in.


% Test if all needed information is given to visualize the splitting
% functions
if ~exist('obj','var') || isempty(obj)
    fprintf('\n No local model network object passed.\n \n');
    return
end
if ~exist('numberOfLocalModels','var') || isempty(numberOfLocalModels)
    numberOfLocalModels = obj.suggestedNet;
end
if ~exist('inputDimensions','var') || isempty(inputDimensions)
    if size(obj.input,2) == 2
        inputDimensions = [1 2];
    else
        inputDimensions = 1;
    end
end
if size(inputDimensions,2)>2
    fprintf(['\n Number of chosen inputs is greater than two. Because only\n two ',...
        'dimensions can be visualized only the first two entries\n (%1.0f %1.0f) of ',...
        'the corresponding vector are utilized.\n \n'],inputDimensions(1),inputDimensions(2));
    inputDimensions = inputDimensions(1:2);
end
if ~exist('axesHandle','var') || isempty(axesHandle)
    axesHandle = gca;
end

% Options for the visualization
pointsPerAxis2d = 40;
pointsPerAxis1d = 1001;

% 1d or 2d are the only possibilities
inputDimensionality = size(inputDimensions,2);

% Minima and Maxima
minima = min(obj.unscaledInput);
maxima = max(obj.unscaledInput);

% Determine the leaf models for the chosen number of local models
leafModels = find(obj.history.leafModelIter{1,numberOfLocalModels});

% Distinguish between sigmoidal and gaussian splitting functions
if isa(obj,'gaussianOrthoGlobalModel')
    
    splittingFcn = @(object,zRegressor) object.calculateMSF(zRegressor);
    
elseif isa(obj,'sigmoidGlobalModel')
    
    splittingFcn = @(object,zRegressor) object.calculatePsi(zRegressor);
    
end

% Hold currently chosen axis
hold(axesHandle,'on');
grid(axesHandle,'on');

% Initialize cell array, that contains all plot handles
plotHandles = cell(1,numberOfLocalModels);

if inputDimensionality == 1
    
    % Create Data
    if size(obj.input,2) == 1
        zReg = linspace(minima(inputDimensions(1)),maxima(inputDimensions(1)),pointsPerAxis1d)';
    else
        zReg = ones(pointsPerAxis1d,1) * mean([minima;maxima]);
        zReg(:,inputDimensions) = linspace(minima(inputDimensions(1)),maxima(inputDimensions(1)),pointsPerAxis1d)';
    end
    
    for ii=1:numberOfLocalModels
        
        % Evaluate values of the current splitting function
        yTmp = splittingFcn(obj.localModels(1,leafModels(ii)),zReg);
        
        % Plot the splitting function for every local model
        plotHandles{1,ii} = plot(axesHandle,zReg(:,inputDimensions(1)),yTmp);
        
    end
    
    set([plotHandles{:}],'Color',[0 0 0],'LineWidth',2,'LineStyle',':');
    
elseif inputDimensionality == 2
    
    % Create Grid Data
    [z1g,z2g] = meshgrid(linspace(minima(inputDimensions(1)),maxima(inputDimensions(1)),pointsPerAxis2d),linspace(minima(inputDimensions(2)),maxima(inputDimensions(2)),pointsPerAxis2d));
    
    % create an input
    if size(obj.input,2) == 2
        zReg = [z1g(:) z2g(:)];
    else
        zReg = ones(length(z1g(:)),1) * mean([minima;maxima]);
        zReg(:,inputDimensions) = [z1g(:) z2g(:)];
    end
    
    for ii=1:numberOfLocalModels
        
        if isprop(obj.localModels(1,leafModels(ii)),'splittingParameter') && isempty(obj.localModels(1,leafModels(ii)).splittingParameter)
            % No parameters are here in the first iteration of an hilomot
            % model
            yTmp = ones(pointsPerAxis2d*pointsPerAxis2d,1);
        else
            % Evaluate values of the current splitting function
            yTmp = splittingFcn(obj.localModels(1,leafModels(ii)),zReg);
        end
        
        yG = reshape(yTmp,pointsPerAxis2d,pointsPerAxis2d);
        
        % Plot the splitting function for every local model
        plotHandles{1,ii} = surf(axesHandle,z1g,z2g,yG);
        
    end
    
end

end