approximation by nonlinear optimization THIS IS NO USER FUNCTION
0001 function out = approxNonlinear(hcMatrix, dec, par) 0002 % approximation by nonlinear optimization 0003 % 0004 % THIS IS NO USER FUNCTION 0005 0006 % The elk-library: convex geometry applied to crystallization modeling. 0007 % Copyright (C) 2013 Alexander Reinhold 0008 % 0009 % This program is free software: you can redistribute it and/or modify it 0010 % under the terms of the GNU General Public License as published by the 0011 % Free Software Foundation, either version 3 of the License, or (at your 0012 % option) any later version. 0013 % 0014 % This program is distributed in the hope that it will be useful, but 0015 % WITHOUT ANY WARRANTY; without even the implied warranty of 0016 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0017 % General Public License for more details. 0018 % 0019 % You should have received a copy of the GNU General Public License along 0020 % with this program. If not, see <http://www.gnu.org/licenses/>. 0021 0022 %! DEV: I should allow to control special parameters from here as there 0023 % are: magnitude of upper/lower bounds; scaling (on/off/factor to 0024 % use); objective/coordinate tolerances; stage/trial points; 0025 %/ 0026 0027 %% Initialize projection data 0028 if strcmpi(par.projectionType, 'confined') 0029 % this operation is only required for "confined" projections 0030 sampleProjectionData = generateSampleProjectionData(hcMatrix, dec, par); 0031 else 0032 sampleProjectionData = []; 0033 end 0034 0035 %% Initialize projection (Grassmann data and first guess) 0036 % we use the PCA-approximation as the initial guess for optimization. 0037 % pcaProjectionMatrix = approxPca(dataOriginal, par); 0038 % It is used as the origin for the grassmann coordinates. The following 0039 % also calculates a coordinate transformation that normalizes the data 0040 % points. Therefore, it handles the options 'manualEnclosingCone' and 0041 % 'manualDataScaling' (not set implemented). 0042 grassStruct = obtainGrassStruct(hcMatrix, par); 0043 % we add the handle for calculating projection matrices (private function) 0044 grassStruct.projectionMatrixHandle = @(x)(... 0045 computeProjectionFromGrass(x, grassStruct)); 0046 0047 %% parameters for optimizer fmincon 0048 % algorithm setup - use 'interior-point' or 'sqp' only because they satisfy 0049 % the constraints in all iterations. 0050 % was 'active-set' before 0051 problem.solver = 'fmincon'; 0052 problem.options = optimset('algorithm', 'sqp', ... 0053 'display', 'off'); 0054 problem.x0 = zeros(grassStruct.grassDim, 1); 0055 0056 % use some upper and lower boundaries as the subspaces are periodic in x. 0057 % Note, however that simple lower/upper bounds are not the tightest 0058 % bounds that can be used. 0059 factor = 0.5*pi*sqrt(min(par.targetDim, dec.nC - par.targetDim)); 0060 problem.lb = -factor*ones(grassStruct.grassDim, 1); 0061 problem.ub = factor*ones(grassStruct.grassDim, 1); 0062 0063 % prohibit invalid subspaces that have an empty intersection with the 0064 % validity cone: nonlinear inequality condition f(x) <= 0 0065 0066 %! Dev: removed nonlcon: objective function is bad there, anyway and the 0067 % affected region should be scaled out successfully. 0068 % problem.nonlcon = @(x)(genericNonlinearCondition(x, grassStruct)); 0069 %/ 0070 problem.nonlcon = @(x)(deal(max(0, sum(x.^2) - factor^2), [])); 0071 0072 % the objective function handles several options: 'projectionType', 0073 % 'errorType', 'sampleWeightVector' and 'measureWeightVector' 0074 problem.objective = @(x)(genericObjectiveFunction(x, hcMatrix, dec, par, ... 0075 grassStruct, sampleProjectionData)); 0076 0077 % tolerance settings - the following should be sufficiently accurate for 0078 % approximations 0079 problem.options.TolFun = abs(problem.objective(problem.x0)) * 1e-6; 0080 problem.options.TolX = 1e-4; 0081 0082 %% Run optimization 0083 if strcmpi(par.method, 'nonlinearLocal') 0084 problem.Display = 'iter'; 0085 [x, f, exitFlag, outStruct] = fmincon(problem); 0086 0087 % provide dummy-gs 0088 gs = []; 0089 0090 elseif strcmpi(par.method, 'nonlinear') 0091 % Configuration and run: GlobalSearch 0092 gs = GlobalSearch('Display', 'iter', ... 0093 'StartPointsToRun', 'bounds', ... 0094 'NumTrialPoints', par.numTrialPoints, ... 0095 'NumStageOnePoints', par.numStageOnePoints); 0096 [x dumpA dumpB dumpC outStruct] = run(gs, problem); 0097 else 0098 error('elk:approximation:wrongInput', ['Optimization method ''' 0099 par.method 'is unknown']); 0100 end 0101 0102 % catch "no result" 0103 if isempty(x) 0104 error('elk:decomposition:noSolution', ['No optimal solution was ' ... 0105 'obtained. Likely, this is an internal error and should not happen.']); 0106 end 0107 0108 %% output data 0109 out.projectionMatrix = computeProjectionFromGrass(x, grassStruct); 0110 out.mappingReducedToEmbedding = orth(out.projectionMatrix); 0111 out.nonlinearData.grassSruct = grassStruct; 0112 out.nonlinearData.sampleProjectionData = sampleProjectionData; 0113 out.nonlinearData.fminconStruct = problem; 0114 out.nonlinearData.globalObject = gs; 0115 out.nonlinearData.globalResult = outStruct; 0116 out.nonlinearData.globalX = x; 0117 0118 0119 %! Debug: 0120 % for iLocal = 1:length(allMinima) 0121 % thisX = allMinima(iLocal).X; 0122 % thisF = allMinima(iLocal).Fval; 0123 % plot3(thisX(1), thisX(2), thisF, 'kx', 'markerSize', 12); 0124 % 0125 % for iStart = 1:length(allMinima(iLocal).X0) 0126 % thisX = allMinima(iLocal).X0{iStart}; 0127 % plot3(thisX(1), thisX(2), 0, 'kd', 'markerSize', 8); 0128 % end 0129 % 0130 % end 0131 % error('stop!') 0132 %/