% BlackboxModelOptimizer: Class to interact with globaloptim.bmo.solver and
% perform optimization iterations.

%   Copyright 2018-2023 The MathWorks, Inc.

classdef BlackboxModelOptimizer
    %
    properties(Hidden) % SetAccess=private
        % Reference to the solver is cleaned up by onCleanup object
        trialRequest
        statusRequest
        resultsRequest
        saveRequest
        state = struct('exitflag', [], 'msg', '', 'startTime', [], ...
            'displayState',[] ,'elapsedtime',0, 'funccount',0, 'iter',0, ... 
            'rngstate', [], 'firstEvalOK',false, 'forceHeaderDisplay',false, ...
            'stallCount', 0, 'elapsedTimeSinceCheckpoint',0);
        
        info = struct('doOutput',false, 'nvars', [], ...
            'outputfcn',[], 'doPlot',false, 'plotfcn',[], 'rngstate', [], ...
            'timeoffset',0);
        
        Verbosity
        
        % Struct containing solver status
        solverStatus
        % Results structure
        results
        % optimoptions: a separate object to hold properties
        options
        % Handle to the expensive model
        expensive
        % Reference to the modelMgr
        modelMgr
        % Reference to the solver used for optimization
        solverRef
        
    end
    
    properties(Hidden, Constant)
        varName = 'X';
        changingVars = {'state', 'results'};
        
        constantVars = {'Verbosity', 'solverRef', 'options', 'trialRequest', ... 
            'statusRequest', 'resultsRequest', 'solverStatus', 'saveRequest', ...
            'expensive', 'options','results','info'};
        EFLAG_FIXED = 10;
        EFLAG_LINEAR_INFEAS = -2;
        STALL_NO_POINTS = 0;
        
    end
    
    properties(Hidden) % , SetAccess=private, GetAccess=public
        myCleanUp
        interrupt = false
        savedOnce = false
    end
    
    methods(Access = public)
        
        function self = BlackboxModelOptimizer(expensive,lb,ub,intcon, ... 
                Aineq,bineq,Aeq,beq,options)
            
            if ~isstruct(expensive)
                error(message('globaloptim:surrogateopt:modelInvalidStruct'))
            end
            % Save the current rng state now
            self.info.rngstate = rng;
            
            if nargin == 2
                options = lb;
            elseif nargin < 9
                % options is supposed to be a struct.
                options = prepareOptionsForSolver(optimoptions(@Surrogateopt), ...
                    'surrogateopt');
            end
            
            if isfield(expensive, 'CheckpointFile')
                % Restore problem and the solver
                self = globaloptim.bmo.restoreSolver(self, ... 
                    expensive.CheckpointFile,options);
                
                % Setup a cleanup object for the solver and cancel any pending
                % job on workers.
                self.myCleanUp = onCleanup( @() deleteSolver(self.solverRef.id, ... 
                    self.modelMgr));
                
                % Save some data to the file now so if there are any errors
                % we stop now.
                if ~isempty(self.saveRequest)
                    self.saveState(true,false);
                end                
                % Mock solver run from restored state
                self.interrupt = self.pastIterations();

                % Start the time counter now.
                self.state.startTime = tic;

            else
                % Setup new problem and the solver
                self = globaloptim.bmo.createSolver(self,expensive,lb,ub, ... 
                    intcon,Aineq,bineq,Aeq,beq,options);
                if ~isempty(self.solverRef)
                    self.myCleanUp = onCleanup( @() deleteSolver(self.solverRef.id, ... 
                        self.modelMgr));
                end
            end
            
            if ~isempty(self.state.exitflag)
                % Check for bound infeasibility.
                return
            end
            
        end
        
        %==================================================================
        function self = optimize(self)
            
            if ~isempty(self.state.exitflag)% == EFLAG_LINEAR_INFEAS, EFLAG_FIXED
                % if exitflag is non-empty, return now.
                return
            end
            
            % Check for solver setup if that was ok. 
            self.solverStatus = globaloptim.bmo.solver(self.statusRequest);
            if self.solverStatus.retcode == -100
                % This is unlikely to happen because self.createSolver
                % should catch any potential error.
                error(self.solverStatus.messageID,self.solverStatus.message);
            end
            
            if self.info.doPlot || self.info.doOutput
                if isempty(self.state.displayState)
                    self.state.displayState = 'init';
                end
                self.interrupt = self.outputAndPlots();
            end
            userDisplay = self.Verbosity > 0 ||  self.info.doPlot || self.info.doOutput;
            stop = self.stopTests();
            
            % Checkpoint to file
            if ~isempty(self.saveRequest)
                self.saveState(true,false);
                self.savedOnce = true;
            end

            while ~stop

                if self.modelMgr.availableWorkers() > 0
                    nextTrial = globaloptim.bmo.solver(self.trialRequest);
                    
                    if nextTrial.retcode == -100
                        % Ensure no error
                        error(nextTrial.messageID,nextTrial.message);
                    elseif nextTrial.retcode == self.STALL_NO_POINTS
                        % Did solver failed to create a trial?
                        self.state.stallCount = self.state.stallCount + 1;
                        % If stalled, we must clear the evaluation queue
                        % when using vectorized.
                        self.modelMgr.evaluate_vectorized_now();
                    else
                        self.state.stallCount = 0;
                        % Non-blocking (if using parfeval)
                        self.modelMgr.evaluate(nextTrial);
                    end
                    
                end
                
                % Check if ModelManager queue can be flushed
                if self.modelMgr.isTrialInQueue() && self.modelMgr.responseAvailable()
                    % Update the solver with response.
                    trialResponse = self.modelMgr.getNext();                   
                    response = globaloptim.bmo.solver(trialResponse);
                    
                    % Ensure no errors
                    if response.retcode == -100
                        error(response.messageID,response.message);
                    end
                    % Update counters for display/stopping
                    self.state.elapsedtime = toc(self.state.startTime) + self.info.timeoffset;
                    self.state.funccount = self.modelMgr.EvaluationCount;                    
            
                    self.solverStatus = globaloptim.bmo.solver(self.statusRequest);
                    if self.solverStatus.exitflag == -99
                        % Check if initial point and response are consistent
                       error(message('globaloptim:surrogateopt:InconsistentFunctionResponse'))
                    end
                    
                    if self.solverStatus.surrogateReset
                        self.modelMgr.cancelPending();
                    end
                    if userDisplay
                        self.state.displayState = 'iter';
                        self = self.iterativeUpdate(); % Iterative updates
                    end
                    % If objective/constraints requested a stop then stop
                    self.interrupt = self.interrupt || trialResponse.Stop;
                    
                    % Iteration count updated after the output/plot
                    % function is called so that iteration is one less than
                    % funccount (after funccount > 0). This is to match
                    % fmincon convention.
                    self.state.iter = self.state.iter + 1;
                    
                end
                % Check stopping criteria
                stop = self.stopTests();
                
                % Checkpoint to file
                if ~isempty(self.saveRequest)
                    self.saveState(false,stop);
                end

            end
            
            if userDisplay
                self.state.displayState = 'done';
                self.solverStatus = globaloptim.bmo.solver(self.statusRequest);
                self = self.iterativeUpdate(); % Iterative updates
            end
        end
        
        %==================================================================
        function [x,fval,exitflag,output,trials] = solution(self)
            
            exitflag = [];
            tolCon = self.options.ConstraintTolerance;
            
            if ~isempty(self.state.exitflag)
                % Stopped in pre-processing
                
                if self.state.exitflag == self.EFLAG_FIXED
                    % All variables are fixed or unique solution.
                    
                    fixedX = globaloptim.bmo.solver(self.trialRequest);
                    % Perform model evaluation
                    self.modelMgr.evaluate(fixedX);
                    fixedResponse = self.modelMgr.getNext();
                    % Update output 
                    fval = fixedResponse.Fval;
                    ineq = fixedResponse.Ineq;
                    x = fixedResponse.(self.varName);
                    trials.X = x;
                    trials.Fval = fval;
                    trials.Ineq = ineq(:)';
                    output.funccount = self.modelMgr.EvaluationCount;                    
                    output.ineq = ineq;
                    output.constrviolation = max([0, ineq]);
                    
                    if output.constrviolation > tolCon
                        % Override exitflag/message for infeasible ineq
                       self.state.exitflag = -2;
                       self.state.msg = ... 
                           getString(message('globaloptim:surrogateopt:FixedBoundInfeasibleIneq'));
                    end
                elseif self.state.exitflag == self.EFLAG_LINEAR_INFEAS
                    % Infeasible bounds/linear
                    x = [];
                    fval = [];
                    trials = [];
                    output.funccount = 0;
                    output.ineq = [];
                    output.constrviolation = [];
                
                end
                output.elapsedtime = toc(self.state.startTime) + ...
                    self.info.timeoffset;
                
                exitflag = self.state.exitflag;
                
            else % Normal optimization run
                self.results = globaloptim.bmo.solver(self.resultsRequest);
                trials = self.results.trialData;                
                output.elapsedtime = toc(self.state.startTime) + ...
                    self.info.timeoffset;                
                output.funccount = self.state.funccount;
                output.constrviolation = self.results.constrviolation;
                output.ineq = self.results.ineq;
                % Get fval and x as non-empty even if x is infeasible
                fval = self.results.fval;
                x = self.results.(self.varName);
                stalled = self.state.stallCount >= self.options.MaxStallIterations;

                
                if isempty(x) && stalled
                    % This may happen if the solver cannot find even one
                    % integer feasible point.
                    self.state.msg = sprintf('%s\n%s', ... 
                        getString(message('globaloptim:surrogateopt:InfeasibleIntegerStall')), ... 
                        getString(message('globaloptim:surrogateopt:StalledMsg')));
                    exitflag = -2;
                    
                elseif ~isempty(x) && output.constrviolation > tolCon && ...
                        ~self.interrupt
                    % Infeasible
                    exitflag = -2;
                    % Set message
                    if stalled
                        % infeasible and stopped due to stalling
                        self.state.msg = sprintf('%s\n%s', ... 
                            getString(message('globaloptim:surrogateopt:InfeasibleIneqStalledMsg')), ...
                            getString(message('globaloptim:surrogateopt:StalledMsg')));
                    else
                        self.state.msg = ... 
                            getString(message('globaloptim:surrogateopt:InfeasibleIneq'));
                    end
                    
                elseif ~isempty(x) && output.constrviolation <= tolCon && ...
                       stalled && ~self.interrupt

                    exitflag = 3;
                        self.state.msg = sprintf('%s\n%s', ... 
                            getString(message('globaloptim:surrogateopt:FeasibleStalledMsg')), ...
                            getString(message('globaloptim:surrogateopt:StalledMsg')));
                end
            end
            
            % Get exitflag and messages
            if ~isempty(fval) && fval < self.options.ObjectiveLimit 
                % Do this check first so we can use this exitflag even when
                % solver exited with eflag = 10 (all fixed variables).
                
                if (~isempty(output.constrviolation) && output.constrviolation > tolCon)
                    % override.
                    exitflag = -2;
                    output_msg = getString(message('globaloptim:surrogateopt:InfeasibleIneq'));
                else
                    exitflag = 1;
                    output_msg = getString(message('globaloptim:surrogateopt:ObjectiveLimitReached'));
                end
                
            elseif ~isempty(exitflag)
                % Cases when exitflag is already assigned; populate output
                output_msg = self.state.msg;
                
            elseif self.state.elapsedtime > self.options.MaxTime
                exitflag = 0;
                output_msg = getString(message('globaloptim:surrogateopt:MaxTimeExceeded'));
                
            elseif self.modelMgr.noMoreEvaluations() 
                exitflag = 0;
                output_msg = getString(message('globaloptim:surrogateopt:MaxFvalExceeded'));
                
            elseif self.interrupt
                exitflag = -1;
                output_msg = getString(message('globaloptim:surrogateopt:FcnInterrupt'));
                
            else % Should not reach here
                exitflag = -100; % What happened here?
                output_msg = 'Unknown Error';
            end
            
            output.rngstate = self.info.rngstate;
            output.message = output_msg;            
            if self.Verbosity > 0
                disp(output.message);
            end
        end
        
        %==================================================================
    end
    
    methods(Hidden)
        
        function self = copyOptions(self, opts)
            % Keep some options as member variables (easy access)
            self.Verbosity = opts.Verbosity;
            
            % Handle output functions
            fcn = opts.OutputFcn;
            if isempty(fcn)
                self.info.doOutput = false;
                self.info.outputfcn = {};
            else
                self.info.doOutput = true;
                % Need to parse OutputFcn to support cell-array syntax for OutputFcn.
                self.info.outputfcn = createCellArrayOfFunctions(fcn,'OutputFcn');
            end
            
            % Handle plot functions
            fcn = opts.PlotFcn;
            if isempty(fcn)
                self.info.doPlot = false;
                self.info.plotfcn = {};
            else
                self.info.doPlot = true;
                % Need to parse PlotFcn to support cell-array syntax for PlotFcns.
                self.info.plotfcn = createCellArrayOfFunctions(fcn,'PlotFcn');
            end
            
        end

    end
    
    methods(Access=private)
        
        %==================================================================
        function TF = stopTests(self)
            fval = self.solverStatus.fval;
            constrviolation = max([0, self.solverStatus.ineq]);
            tolCon = self.options.ConstraintTolerance;
           
            % Test for exitflag == 0 or -1/-2
            TF = self.modelMgr.noMoreEvaluations() || ...            
                self.state.elapsedtime > self.options.MaxTime || ...
                self.interrupt || ~isempty(self.state.exitflag) || ...
                self.state.stallCount >= self.options.MaxStallIterations;
           if TF
               return
           end
           
            % Test for exitflag == 1
             TF = constrviolation <= tolCon && ...
                 (~isempty(fval) && fval < self.options.ObjectiveLimit);

        end
        
        %==================================================================
        function self = iterativeUpdate(self)
                        
            % Handle iterative display differently when the first few
            % points fails (returns NaN for fval/ineq).
            if (self.Verbosity > 0 || ... 
                self.info.doPlot   || ... 
                self.info.doOutput)&& ~self.state.firstEvalOK
                % Is fval NaN? current iterate is always non-empty but best
                % will be empty if function returned NaN (for fval or
                % ineq).
                fvalFailed = isempty(self.solverStatus.fval) && ...
                    ~isempty(self.solverStatus.current.fval) && ...
                    isnan((self.solverStatus.current.fval));
                % If ineq NaN?
                ineqFailed = isempty(self.solverStatus.ineq) && ...
                    ~isempty(self.solverStatus.current.ineq) && ...
                    any(isnan((self.solverStatus.current.ineq)));
                
                self.state.firstEvalOK =  ~fvalFailed && ~ineqFailed;
                
                if ~self.state.firstEvalOK               
                    fprintf('%s\n',getString(message('globaloptim:surrogateopt:InitEvalFailed')));
                    self.state.forceHeaderDisplay = true;
                end
            end
            
            if self.Verbosity > 0 && self.state.firstEvalOK
                % Print iterative display (table style) when we have at
                % least one successful function evaluation.
                self = self.printIterativeData();
            end
            
            if self.info.doPlot || self.info.doOutput
                % Call plot functions irrespective of the failure. We
                % expect that plots can handle NaN and we do not set [] for
                % fval/ineq
                self.interrupt = self.outputAndPlots();
            end
            
        end
        
        %==================================================================
        function self = printIterativeData(self)
            
            if strcmp(self.state.displayState,'done') || self.state.iter < 0
                % Skip iterative display if done iterating.
                return
            end
            
            if self.state.iter == 0 || self.state.forceHeaderDisplay
                self.showProbelmInfo();
            end
            
            self.showHeaders();
            
            if self.state.forceHeaderDisplay
                % If header was displayed then turn this flag off.
                self.state.forceHeaderDisplay = false;
            end

            if self.Verbosity >= 2
               % Prepare data to print depending on the problem
                if isempty(self.solverStatus.fval)
                    % feasibility problem
                    best = sprintf('%- 3g',self.solverStatus.nFeasPoints);
                    current = sprintf('%- 8.4e',self.solverStatus.current.constrviolation);
                elseif isempty(self.solverStatus.ineq)
                    % no ineq
                    best = sprintf('%- 8.4e',self.options.ObjectiveSenseManager.updateFval(self.solverStatus.fval));
                    current = sprintf('%- 8.4e',self.options.ObjectiveSenseManager.updateFval(self.solverStatus.current.fval));
                else
                    % fval and ineq
                    if self.solverStatus.constrviolation > self.options.ConstraintTolerance
                        % infeasible
                        best = sprintf('(%-s, %- 8.4e)','-',self.solverStatus.constrviolation);
                    else
                        % feasible
                        best = sprintf('(%- 8.4e, %-s)',self.options.ObjectiveSenseManager.updateFval(self.solverStatus.fval),'-');
                    end
                    
                    if self.solverStatus.current.constrviolation > self.options.ConstraintTolerance
                        % infeasible
                        current = sprintf('(%-s, %- 8.4e)','-',self.solverStatus.current.constrviolation);
                    else
                        % feasible
                        current = sprintf('(%- 8.4e, %-s)',self.options.ObjectiveSenseManager.updateFval(self.solverStatus.current.fval),'-');
                    end                    
                end
                
            end
            
            if self.Verbosity > 2
                if self.solverStatus.surrogateReset
                    surrogateReset = 'Yes';
                else
                    surrogateReset = '';
                end                
                
                fprintf('%4d        %- 6.2f     %- 10s     %- 10s   %- 7s    %3s\n',self.state.funccount, ...
                    self.solverStatus.elapsedtime, best, current, ...
                    self.solverStatus.current.flag, surrogateReset);
                
                if self.solverStatus.checkpointResume                    
                    fprintf('%s\n',getString(message('globaloptim:surrogateopt:ResumeFromCheckPointMsg')));
                end
                
            elseif self.Verbosity == 2
                fprintf('%4d        %- 6.2f     %- 10s     %- 10s   %- 7s\n',self.state.funccount, ...
                    self.solverStatus.elapsedtime, best, current, ...
                    self.solverStatus.current.flag);
                
                if self.solverStatus.checkpointResume
                    fprintf('%s\n',getString(message('globaloptim:surrogateopt:ResumeFromCheckPointMsg')));
                end
            end
            
        end
        
        %==================================================================
        function stop = outputAndPlots(self)
            
            optimValues.currentX = self.solverStatus.current.X;
            optimValues.currentFval = self.solverStatus.current.fval;
            optimValues.currentIneq = self.solverStatus.current.ineq;
            optimValues.currentConstrviolation = self.solverStatus.current.constrviolation;            
            optimValues.currentFlag = self.solverStatus.current.flag;
            
            optimValues.incumbentX = self.solverStatus.incumbent.X;
            optimValues.incumbentFval = self.solverStatus.incumbent.fval;
            optimValues.incumbentIneq = self.solverStatus.incumbent.ineq;
            optimValues.incumbentConstrviolation = self.solverStatus.incumbent.constrviolation;
            optimValues.incumbentFlag = self.solverStatus.incumbent.flag;

            xCurrent = self.solverStatus.X;
            optimValues.fval = self.solverStatus.fval;
            optimValues.ineq = self.solverStatus.ineq; % Inequalities
            optimValues.constrviolation = self.solverStatus.constrviolation;
            optimValues.flag = self.solverStatus.flag;
            optimValues.iteration = self.state.iter;
            optimValues.elapsedtime = self.solverStatus.elapsedtime;
            optimValues.funccount = self.state.funccount;
            
            optimValues.surrogateReset = self.solverStatus.surrogateReset;
            optimValues.surrogateResetCount = self.solverStatus.surrogateResetCount;
            
            optimValues.checkpointResume = self.solverStatus.checkpointResume;
            optimValues.checkpointResumeCount = self.solverStatus.checkpointResumeCount;
            
            stop = false;            
            fcn = self.info.outputfcn;
            if ~isempty(fcn)
                switch self.state.displayState
                    case 'init'
                        stop = callAllOptimOutputFcns(fcn,xCurrent,optimValues,self.state.displayState) || stop;
                        checkForPlotFigure();
                    case {'iter','interrupt'}
                        stop = callAllOptimOutputFcns(fcn,xCurrent,optimValues,self.state.displayState) || stop;
                    case 'done'
                        callAllOptimOutputFcns(fcn,xCurrent,optimValues,self.state.displayState);
                end
            end
            
            if ~self.state.firstEvalOK
                % Override optimValues (so fval/ineq are non-empty) to
                % prevent any plot failures. These values are NaN or
                % useless anyway.
                xCurrent = self.solverStatus.current.X;
                optimValues.fval = self.solverStatus.current.fval;
                optimValues.ineq = self.solverStatus.current.ineq; % Inequalities
                optimValues.constrviolation = self.solverStatus.current.constrviolation;
                optimValues.flag = self.solverStatus.current.flag;
                
                optimValues.incumbentX = xCurrent;
                optimValues.incumbentFval = optimValues.fval;
                optimValues.incumbentIneq = optimValues.ineq;
                optimValues.incumbentConstrviolation = optimValues.constrviolation;
                optimValues.incumbentFlag = optimValues.flag;
                
            end
            
            % Call plot functions
            fcn = self.info.plotfcn;
            if ~isempty(fcn)
                switch self.state.displayState
                    case 'init'
                        stop = callAllOptimPlotFcns(fcn,xCurrent,optimValues,self.state.displayState) || stop;
                        checkForPlotFigure();
                    case 'iter'
                        stop = callAllOptimPlotFcns(fcn,xCurrent,optimValues,self.state.displayState) || stop;
                    case 'done'
                        callAllOptimPlotFcns(fcn,xCurrent,optimValues,self.state.displayState);
                end
            end
            
            stop = stop || self.interrupt;

            function checkForPlotFigure()
                % Add options as userdata in the figure for plot
                % functions to use if needed.
                fname = getString(message('MATLAB:optimfun:funfun:optimplots:WindowTitle'));
                fig = findobj(0,'Type','figure','name',fname);
                if ~isempty(fig)
                    set(fig,'UserData',self.options)
                end
            end
        end
        
        %==================================================================
        function showHeaders(self)
            
            if rem(self.state.iter,30) == 0 || self.state.forceHeaderDisplay
                if isempty(self.solverStatus.fval)
                    % feasibility problem
                    if self.Verbosity > 2
                        fprintf('\nF-count    Time       Number of           Current   Trial    Surrogate\n')
                        fprintf('        (seconds)     Feas Pts          Max Infeas  type      reset\n')
                    elseif self.Verbosity == 2
                        fprintf('\nF-count    Time       Number of           Current   Trial\n')
                        fprintf('        (seconds)     Feas Pts          Max Infeas  type\n')
                    end
                elseif isempty(self.solverStatus.ineq)
                    % no ineq
                    if self.Verbosity > 2
                        fprintf('\nF-count       Time        Best          Current        Trial    Surrogate\n')
                        fprintf('             (seconds)    Fval          Fval           type      reset\n')
                    elseif self.Verbosity == 2
                        fprintf('\nF-count       Time        Best          Current         Trial\n')
                        fprintf('             (seconds)    Fval          Fval            type\n')
                    end
                else
                    if self.Verbosity > 2
                        fprintf('\nF-count    Time            Best               Current          Trial    Surrogate\n')
                        fprintf('          (seconds)    (Fval, Infeas)       (Fval, Infeas)     type      reset\n')
                    elseif self.Verbosity == 2
                        fprintf('\nF-count    Time            Best               Current          Trial\n')
                        fprintf('          (seconds)    (Fval, Infeas)       (Fval, Infeas)     type\n')
                    end
                end
                
                
            end
            
        end
        
        function showProbelmInfo(self)
            if self.Verbosity >= 2                
                if numel(self.solverStatus.fval) > 0                    
                    fprintf('%s\n',getString(message('globaloptim:surrogateopt:ScalarObjectiveFunction')));
                else
                    fprintf('%s\n',getString(message('globaloptim:surrogateopt:NoObjective')));
                end
                
                fprintf('%s\n',getString(message('globaloptim:surrogateopt:NumberOfVariables',self.info.nvars)));
                fprintf('%s\n',getString(message('globaloptim:surrogateopt:NumberOfIntVariables', ... 
                    self.info.nInt)));
                fprintf('%s\n',getString(message('globaloptim:surrogateopt:NumberOfLinearInequalities', ... 
                    self.info.mLinIneq)));
                fprintf('%s\n',getString(message('globaloptim:surrogateopt:NumberOfLinearEqualities', ... 
                    self.info.mLinEq)));
                
                fprintf('%s\n',getString(message('globaloptim:surrogateopt:NumberOfInequalities', ... 
                    numel(self.solverStatus.ineq))));
            end
        end
        

        %==================================================================
        function saveState(self, saveConstVarsOnly, lastSave)
            
            if isempty(self.saveRequest)
                return;
            end

            if ~self.savedOnce
                BlackboxModelConstantData = struct();
                for ii = 1:length(self.constantVars)
                    var = self.constantVars{ii};
                    BlackboxModelConstantData.(var) = self.(var);
                end
                if exist(self.options.CheckpointFile,'file')
                    CheckPointData = struct(); loadCheckPointFile()                                        
                    CheckPointData.BlackboxModelConstantData = BlackboxModelConstantData;
                else
                    CheckPointData.BlackboxModelConstantData = BlackboxModelConstantData;                    
                end                
                globaloptim.bmo.checkpointToFile(self.options.CheckpointFile, CheckPointData);
            end

            if saveConstVarsOnly
                return;
            end

            % Save only after checkpointInterval interval has passed. test
            % after funccount > 1 so we have saved at least one data point.
            if ~lastSave
                if self.state.funccount > 1 && ...
                        (self.state.elapsedtime - self.state.elapsedTimeSinceCheckpoint) < ...
                        self.options.CheckpointInterval
                    return;
                end
            end

            self.state.rngstate = rng;
            % Save state
            BlackboxModelOptimizerData = struct();
            for ii = 1:length(self.changingVars)
                var = self.changingVars{ii};
                BlackboxModelOptimizerData.(var) = self.(var);
            end

            loadCheckPointFile()
            CheckPointData.BlackboxModelOptimizerData = BlackboxModelOptimizerData;
            globaloptim.bmo.checkpointToFile(self.options.CheckpointFile, CheckPointData);

            % Request to save solver's state
            response = globaloptim.bmo.solver(self.saveRequest);
            if response.retcode == -100
                 error(response.messageID,response.message);
            end

            % Record the checkpoint time
            self.state.elapsedTimeSinceCheckpoint = ... 
                        toc(self.state.startTime) + self.info.timeoffset;

            function loadCheckPointFile()
                % Load 'BlackboxModelOptimizerData from the file (assuming there is no data)
                warnState = warning('query','MATLAB:load:variableNotFound');
                warning('off','MATLAB:load:variableNotFound')
                try
                    load(self.options.CheckpointFile,'CheckPointData')
                    warning(warnState)
                catch ME
                    warning(warnState)
                    rethrow(ME)
                end
            end
        end
        
        %==================================================================
        function stop = pastIterations(self)            

            userDisplay = self.Verbosity > 0 ||  self.info.doPlot || self.info.doOutput;
            stop = false;
            if ~userDisplay
                return;
            end
            
            % Reset the display state 
            self.state.displayState = [];
            nCount = self.state.funccount;
            % Call the solver to get status for past iterations
            for ii = 0:nCount

                self.statusRequest.iter = ii;                
                self.solverStatus = globaloptim.bmo.solver(self.statusRequest);
                self.state.iter = ii-1;
                self.state.funccount = ii;
                self.state.elapsedtime =  self.solverStatus.elapsedtime;

                
                if self.state.iter < 0
                    self.state.displayState = 'init';
                else
                    self.state.displayState = 'iter';
                end
                
                if self.state.iter == 0
                    % Reset this flag to false until fval/ineq is checked.
                    self.state.firstEvalOK = false;
                end
                
                if ii == nCount
                    self.solverStatus.checkpointResume = true;
                    self.solverStatus.checkpointResumeCount = ... 
                        self.solverStatus.checkpointResumeCount + 1; 
                end
                
                self = self.iterativeUpdate(); % Iterative updates
                stop = self.interrupt;
                if stop
                    break;
                end

            end
            
        end
    end
    
end

%==========================================================================
function deleteSolver(id, modelMgr)

if ~isempty(id)
    globaloptim.bmo.solver(struct('id',id,'clear',true));
end

% Cancel any remaining task
modelMgr.cancelAll();

end

