function out = approxNonlinear(hcMatrix, dec, par)
% approximation by nonlinear optimization
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 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/>.

%! DEV: I should allow to control special parameters from here as there
%       are: magnitude of upper/lower bounds; scaling (on/off/factor to
%       use); objective/coordinate tolerances; stage/trial points;
%/

%% Initialize projection data
if strcmpi(par.projectionType, 'confined')
    % this operation is only required for "confined" projections
    sampleProjectionData = generateSampleProjectionData(hcMatrix, dec, par);
else
    sampleProjectionData = [];
end

%% Initialize projection (Grassmann data and first guess)
% we use the PCA-approximation as the initial guess for optimization. 
% pcaProjectionMatrix = approxPca(dataOriginal, par);
% It is used as the origin for the grassmann coordinates. The following 
%   also calculates a coordinate transformation that normalizes the data 
%   points. Therefore, it handles the options 'manualEnclosingCone' and
%   'manualDataScaling' (not set implemented).
grassStruct = obtainGrassStruct(hcMatrix, par);
% we add the handle for calculating projection matrices (private function)
grassStruct.projectionMatrixHandle = @(x)(...
    computeProjectionFromGrass(x, grassStruct));

%% parameters for optimizer fmincon
% algorithm setup - use 'interior-point' or 'sqp' only because they satisfy
%   the constraints in all iterations.
% was 'active-set' before
problem.solver = 'fmincon';
problem.options = optimset('algorithm', 'sqp', ...
                           'display', 'off');
problem.x0 = zeros(grassStruct.grassDim, 1);

% use some upper and lower boundaries as the subspaces are periodic in x.
%   Note, however that simple lower/upper bounds are not the tightest
%   bounds that can be used.
factor = 0.5*pi*sqrt(min(par.targetDim, dec.nC - par.targetDim));
problem.lb = -factor*ones(grassStruct.grassDim, 1);
problem.ub = factor*ones(grassStruct.grassDim, 1);

% prohibit invalid subspaces that have an empty intersection with the
%   validity cone: nonlinear inequality condition f(x) <= 0

%! Dev: removed nonlcon: objective function is bad there, anyway and the
%    affected region should be scaled out successfully.
% problem.nonlcon = @(x)(genericNonlinearCondition(x, grassStruct));
%/
problem.nonlcon = @(x)(deal(max(0, sum(x.^2) - factor^2), []));

% the objective function handles several options: 'projectionType',
%   'errorType', 'sampleWeightVector' and 'measureWeightVector'
problem.objective = @(x)(genericObjectiveFunction(x, hcMatrix, dec, par, ...
      grassStruct, sampleProjectionData));

% tolerance settings - the following should be sufficiently accurate for
%   approximations
problem.options.TolFun = abs(problem.objective(problem.x0)) * 1e-6;
problem.options.TolX = 1e-4;

%% Run optimization
if strcmpi(par.method, 'nonlinearLocal')
    problem.Display = 'iter';
    [x, f, exitFlag, outStruct] = fmincon(problem);
    
    % provide dummy-gs
    gs = [];
    
elseif strcmpi(par.method, 'nonlinear')
    % Configuration and run: GlobalSearch
    gs = GlobalSearch('Display', 'iter', ...
        'StartPointsToRun', 'bounds', ...
        'NumTrialPoints', par.numTrialPoints, ...
        'NumStageOnePoints', par.numStageOnePoints);
    [x dumpA dumpB dumpC outStruct] = run(gs, problem);
else
    error('elk:approximation:wrongInput', ['Optimization method '''
        par.method 'is unknown']);
end

% catch "no result"
if isempty(x)
    error('elk:decomposition:noSolution', ['No optimal solution was ' ...
        'obtained. Likely, this is an internal error and should not happen.']);
end

%% output data
out.projectionMatrix = computeProjectionFromGrass(x, grassStruct);
out.mappingReducedToEmbedding = orth(out.projectionMatrix);
out.nonlinearData.grassSruct = grassStruct;
out.nonlinearData.sampleProjectionData = sampleProjectionData;
out.nonlinearData.fminconStruct = problem;
out.nonlinearData.globalObject = gs;
out.nonlinearData.globalResult = outStruct;
out.nonlinearData.globalX = x;


%! Debug:
% for iLocal = 1:length(allMinima)
%     thisX = allMinima(iLocal).X;
%     thisF = allMinima(iLocal).Fval;
%     plot3(thisX(1), thisX(2), thisF, 'kx', 'markerSize', 12);
%     
%     for iStart = 1:length(allMinima(iLocal).X0)
%     thisX = allMinima(iLocal).X0{iStart};
%     plot3(thisX(1), thisX(2), 0, 'kd', 'markerSize', 8);
%     end
%     
% end
% error('stop!')
%/