matRad_OptimizerIPOPT

0001 classdef matRad_OptimizerIPOPT < matRad_Optimizer
0002 % matRad_OptimizerIPOPT implements the interface for ipopt
0003 %
0004 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % Copyright 2019 the matRad development team.
0007 %
0008 % This file is part of the matRad project. It is subject to the license
0009 % terms in the LICENSE file found in the top-level directory of this
0010 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0011 % of the matRad project, including this file, may be copied, modified,
0012 % propagated, or distributed except according to the terms contained in the
0013 % LICENSE file.
0014 %
0015 % References
0016 %   -
0017 %
0018 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0019     
0020     properties
0021         options
0022         wResult
0023         resultInfo
0024         env
0025     end
0026     
0027     properties (Access = private)
0028         allObjectiveFunctionValues
0029         axesHandle
0030         plotHandle
0031         abortRequested
0032         plotFailed
0033     end
0034     
0035     methods
0036         function obj = matRad_OptimizerIPOPT
0037             %matRad_OptimizerIPOPT
0038             %   Construct an instance of the IPOPT optimizer (mex
0039             %   interface)
0040             
0041             matRad_cfg = MatRad_Config.instance();
0042             
0043             obj.wResult = [];
0044             obj.resultInfo = [];
0045             obj.axesHandle = [];
0046             obj.allObjectiveFunctionValues = [];
0047             obj.abortRequested = false;
0048             
0049             %Set Default Options
0050             if matRad_cfg.logLevel <= 1
0051                 lvl = 0;
0052             elseif matRad_cfg.logLevel <= 2
0053                 lvl = 2;
0054             elseif matRad_cfg.logLevel <= 3
0055                 lvl = 5;
0056             else 
0057                 %There seems to be a problem with higher log levels in
0058                 %IPOPT!
0059                 lvl = 5;
0060             end
0061                 
0062             obj.options.print_level                   = lvl;
0063             obj.options.print_user_options            = 'no';
0064             obj.options.print_options_documentation   = 'no';
0065             
0066             % Termination (C.2)
0067             obj.options.tol                           = 1e-8; % (Opt1)
0068             obj.options.dual_inf_tol                  = 1;    % (Opt2)
0069             obj.options.constr_viol_tol               = 1e-4; % (Opt3)
0070             obj.options.compl_inf_tol                 = 1e-4; % (Opt4), Optimal Solution Found if (Opt1),...,(Opt4) fullfiled
0071             
0072             obj.options.acceptable_iter               = 3;    % (Acc1)
0073             obj.options.acceptable_tol                = 1e10; % (Acc2)
0074             obj.options.acceptable_constr_viol_tol    = 1e10; % (Acc3)
0075             obj.options.acceptable_dual_inf_tol       = 1e10; % (Acc4)
0076             obj.options.acceptable_compl_inf_tol      = 1e10; % (Acc5)
0077             obj.options.acceptable_obj_change_tol     = 1e-3; % (Acc6), Solved To Acceptable Level if (Acc1),...,(Acc6) fullfiled
0078             
0079             obj.options.max_iter                      = matRad_cfg.propOpt.defaultMaxIter;
0080             obj.options.max_cpu_time                  = 3000;
0081             
0082             % Barrier Parameter (C.6)
0083             obj.options.mu_strategy = 'adaptive';
0084             
0085             % Line Sarch (C.8)
0086             %obj.options.accept_every_trial_step = 'yes';
0087             %obj.options.line_search_method = 'cg-penalty';
0088             
0089             % Restoration Phase (C.10)
0090             %obj.options.soft_resto_pderror_reduction_factor = 100;
0091             %obj.options.required_infeasibility_reduction    = 0.9999;
0092             
0093             % Quasi-Newton (C.13)
0094             obj.options.hessian_approximation         = 'limited-memory';
0095             obj.options.limited_memory_max_history    = 6;
0096             obj.options.limited_memory_initialization = 'scalar2';
0097             
0098             % determine once if Matlab or Octave
0099             obj.env = matRad_getEnvironment();
0100             
0101             % for derivate checking
0102             % obj.options.derivative_test              = 'first-order'; % none / first-order / second-order / only-second-order
0103             % obj.options.derivative_test_perturbation = 1e-6; % default 1e-8
0104             % obj.options.derivative_test_tol          = 1e-6;
0105             
0106             if ~matRad_checkMexFileExists('ipopt')
0107                 matRad_cfg.dispError('IPOPT mex interface not available for %s!',obj.env);
0108             end
0109 
0110         end
0111         
0112         function obj = optimize(obj,w0,optiProb,dij,cst)
0113             matRad_cfg = MatRad_Config.instance();
0114             
0115             % set optimization options
0116             
0117             %Set up ipopt structure
0118             ipoptStruct = struct;
0119             
0120             %optimizer options
0121             ipoptStruct.ipopt = obj.options;
0122             
0123             %variable bounds
0124             ipoptStruct.lb = optiProb.lowerBounds(w0);
0125             ipoptStruct.ub = optiProb.upperBounds(w0);
0126             
0127             %constraint bounds;
0128             [ipoptStruct.cl,ipoptStruct.cu] = optiProb.matRad_getConstraintBounds(cst);
0129             
0130             % set callback functions.
0131             
0132             funcs.objective         = @(x) optiProb.matRad_objectiveFunction(x,dij,cst);
0133             funcs.constraints       = @(x) optiProb.matRad_constraintFunctions(x,dij,cst);
0134             funcs.gradient          = @(x) optiProb.matRad_objectiveGradient(x,dij,cst);
0135             funcs.jacobian          = @(x) optiProb.matRad_constraintJacobian(x,dij,cst);
0136             funcs.jacobianstructure = @( ) optiProb.matRad_getJacobianStructure(w0,dij,cst);
0137             funcs.iterfunc          = @(iter,objective,paramter) obj.iterFunc(iter,objective,paramter,ipoptStruct.ipopt.max_iter);
0138             
0139             % Informing user to press q to terminate optimization
0140             matRad_cfg.dispInfo('\nOptimzation initiating...\n');
0141             
0142             % set Callback
0143             qCallbackSet = false;
0144             if ~isdeployed % only if _not_ running as standalone
0145                 switch obj.env
0146                     case 'MATLAB'
0147                         try
0148                             % get handle to Matlab command window
0149                             mde         = com.mathworks.mde.desk.MLDesktop.getInstance;
0150                             cw          = mde.getClient('Command Window');
0151                             xCmdWndView = cw.getComponent(0).getViewport.getComponent(0);
0152                             h_cw        = handle(xCmdWndView,'CallbackProperties');
0153                         
0154                             % set Key Pressed Callback of Matlab command window
0155                             set(h_cw, 'KeyPressedCallback', @(h,event) obj.abortCallbackKey(h,event));
0156                             fprintf('Press q to terminate the optimization...\n');
0157                             qCallbackSet = true;
0158                         catch
0159                             matRad_cfg.dispInfo('Manual termination with q not possible due to failing callback setup.\n');
0160                         end
0161                 end                
0162             end
0163             
0164             %ipoptStruct.options = obj.options;
0165             obj.abortRequested = false;
0166             obj.plotFailed = false;
0167             
0168             % Run IPOPT.
0169             try
0170                 [obj.wResult, obj.resultInfo] = ipopt(w0,funcs,ipoptStruct);
0171             catch ME
0172                 errorString = [ME.message '\nThis error was thrown by the MEX-interface of IPOPT.\nMex interfaces can raise compatability issues which may be resolved by compiling them by hand directly on your particular system.'];
0173                 matRad_cfg.dispError(errorString);
0174             end
0175             
0176             % unset Key Pressed Callback of Matlab command window
0177             if qCallbackSet
0178                 set(h_cw, 'KeyPressedCallback',' ');
0179             end
0180             
0181             obj.abortRequested = false;
0182             % Empty the array of stored function values
0183             obj.allObjectiveFunctionValues = [];
0184         end
0185         
0186         function [statusmsg,statusflag] = GetStatus(obj)
0187             try
0188                 switch obj.resultInfo.status
0189                     case 0
0190                         statusmsg = 'solved';
0191                     case 1
0192                         statusmsg = 'solved to acceptable level';
0193                     case 2
0194                         statusmsg = 'infeasible problem detected';
0195                     case 3
0196                         statusmsg = 'search direction too small';
0197                     case 4
0198                         statusmsg = 'diverging iterates';
0199                     case 5
0200                         statusmsg = 'user requested stop';
0201                     case -1
0202                         statusmsg = 'maximum number of iterations';
0203                     case -2
0204                         statusmsg = 'restoration phase failed';
0205                     case -3
0206                         statusmsg = 'error in step computation';
0207                     case -4
0208                         statusmsg = 'maximum CPU time exceeded';
0209                     case -10
0210                         statusmsg = 'not enough degrees of freedom';
0211                     case -11
0212                         statusmsg = 'invalid problem definition';
0213                     case -12
0214                         statusmsg = 'invalid option';
0215                     case -13
0216                         statusmsg = 'invalid number detected';
0217                     case -100
0218                         statusmsg = 'unrecoverable exception';
0219                     case -101
0220                         statusmsg = 'non-IPOPT exception thrown';
0221                     case -102
0222                         statusmsg = 'insufficient memory';
0223                     case -199
0224                         statusmsg = 'IPOPT internal error';
0225                     otherwise
0226                         statusmsg = 'IPOPT returned unknown status';
0227                 end
0228                 
0229                 if obj.resultInfo.status == 0
0230                     statusflag = 0;
0231                 elseif obj.resultInfo.status > 0
0232                     statusflag = 1;
0233                 else
0234                     statusflag = -1;
0235                 end
0236             catch
0237                 statusmsg = 'No Last IPOPT Status Available!';
0238                 statusflag = -1;
0239             end
0240         end
0241         
0242         function flag = iterFunc(obj,iter,objective,~,~)
0243              
0244             obj.allObjectiveFunctionValues(iter + 1) = objective;
0245             %We don't want the optimization to crash because of drawing
0246             %errors
0247             if ~obj.plotFailed
0248                 try            
0249                     obj.plotFunction();
0250                 catch ME
0251                     matRad_cfg = MatRad_Config.instance();
0252                     %Put a warning at iteration 1 that plotting failed
0253                     matRad_cfg.dispWarning('Objective Function plotting failed and thus disabled. Message:\n%s',ME.message);
0254                     obj.plotFailed = true;
0255                 end                
0256             end
0257             flag = ~obj.abortRequested;
0258         end
0259         
0260         function plotFunction(obj)
0261             % plot objective function output
0262             y = obj.allObjectiveFunctionValues;
0263             x = 1:numel(y);
0264             
0265             if isempty(obj.axesHandle) || ~isgraphics(obj.axesHandle,'axes')
0266                 %Create new Fiure and store axes handle
0267                 hFig = figure('Name','Progress of IPOPT Optimization','NumberTitle','off','Color',[.5 .5 .5]);
0268                 hAx = axes(hFig);
0269                 hold(hAx,'on');
0270                 grid(hAx,'on');
0271                 grid(hAx,'minor');
0272                 set(hAx,'YScale','log');
0273                  
0274                 %Add a Stop button with callback to change abort flag
0275                 c = uicontrol;
0276                 cPos = get(c,'Position');
0277                 cPos(1) = 5;
0278                 cPos(2) = 5;
0279                 set(c,  'String','Stop',...
0280                         'Position',cPos,...
0281                         'Callback',@(~,~) abortCallbackButton(obj));                
0282                 
0283                 %Set up the axes scaling & labels
0284                 defaultFontSize = 14;
0285                 set(hAx,'YScale','log');
0286                 title(hAx,'Progress of Optimization','LineWidth',defaultFontSize);
0287                 xlabel(hAx,'# iterations','Fontsize',defaultFontSize),ylabel(hAx,'objective function value','Fontsize',defaultFontSize);
0288                 
0289                 %Create plot handle and link to data for faster update
0290                 hPlot = plot(hAx,x,y,'xb','LineWidth',1.5,'XDataSource','x','YDataSource','y');
0291                 obj.plotHandle = hPlot;
0292                 obj.axesHandle = hAx;
0293                                 
0294             else %Figure already exists, retreive from axes handle
0295                 hFig = get(obj.axesHandle,'Parent');
0296                 hAx = obj.axesHandle;
0297                 hPlot = obj.plotHandle;
0298             end
0299             
0300             % draw updated axes by refreshing data of the plot handle (which is linked to y and y)
0301             % in the caller workspace. Octave needs and works on figure handles, which
0302             % is substantially (factor 10) slower, thus we check explicitly
0303             switch obj.env
0304                 case 'OCTAVE'
0305                     refreshdata(hFig,'caller');
0306                 otherwise
0307                     refreshdata(hPlot,'caller');
0308             end
0309             drawnow;
0310             
0311             % ensure to bring optimization window to front
0312             figure(hFig);
0313         end
0314         
0315         function abortCallbackKey(obj,~,KeyEvent)
0316             % check if user pressed q
0317             if  get(KeyEvent,'keyCode') == 81
0318                 obj.abortRequested = true;
0319             end
0320         end
0321         
0322         function abortCallbackButton(obj,~,~,~)
0323             obj.abortRequested = true;
0324         end        
0325     end
0326     
0327     methods (Static)
0328         function available = IsAvailable()
0329             available = matRad_checkMexFileExists('ipopt');                   
0330         end
0331     end
0332 end