sos

 - SOS optimization object

 copyright 2009-2012 Blair Armstrong, Christine Watson, David Plaut

    This file is part of SOS

    SOS is free software: you can redistribute it and/or modify
    it for academic and non-commercial purposes
    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.  For commercial or for-profit
    uses, please contact the authors (sos@cnbc.cmu.edu).

    SOS 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.

0001 % - SOS optimization object
0002 %
0003 % copyright 2009-2012 Blair Armstrong, Christine Watson, David Plaut
0004 %
0005 %    This file is part of SOS
0006 %
0007 %    SOS is free software: you can redistribute it and/or modify
0008 %    it for academic and non-commercial purposes
0009 %    under the terms of the GNU General Public License as published by
0010 %    the Free Software Foundation, either version 3 of the License, or
0011 %    (at your option) any later version.  For commercial or for-profit
0012 %    uses, please contact the authors (sos@cnbc.cmu.edu).
0013 %
0014 %    SOS is distributed in the hope that it will be useful,
0015 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
0016 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0017 %    GNU General Public License for more details.
0018 
0019 %    You should have received a copy of the GNU General Public License
0020 %    along with SOS (see COPYING.txt).
0021 %    If not, see <http://www.gnu.org/licenses/>.
0022 
0023 classdef sos < genericStatTest
0024     % Creates and supports optimization objects
0025     %
0026     % Inherits from handle so that by-reference passing can be used with
0027     % the object
0028     %
0029     % PROPERTIES
0030     %     samples % array of samples linked to SOS object
0031     %     hardConstraints % array of hard constraint objects
0032     %     softConstraints % array of soft constraint objects
0033     %     metaConstraints % array of meta constraint objects
0034     %     cost % current cost
0035     %     sCandObj % object used to identify the target item in a sample which may be swapped
0036     %     feederCandObj % object used to identify the item in a population/other sample to be swapped with the target sample item
0037     %     pSwapObj  % object which determines the probability of making a swap based on deltaCost
0038     %     annealObj % object used to anneal temperature in pSwap
0039     %     maxIt % maximum number of iterations that an optimization will run for
0040     %     allData % merge of all data from all samples/ pops linked to sosobj
0041     %     allDataColName % name of columns in allData
0042     %     allDataColMean % mean for each column in allData
0043     %     allDataColStd  % stdev for each column in all data
0044     %     reportInterval % number of iterations between progress reports
0045     %     curIt % current optimization iteration
0046     %     curFreezeIt % number of iterations cost has been frozen at a particular value
0047     %     stopFreezeIt % number of frozen iterations before algorithm stops
0048     %     sosstattests % array of handles to sos  stat tests requested by user
0049     %     statInterval % number of iterations between stat reports
0050     %     statTestReportStyle % style of report ('short' or 'full')
0051     %     deltaCost % delta cost from soft constraints from last attempted flip that met all hard constraints
0052     %     blockSize % number of iterations in a 'block' of trials for which average flip history /minDeltaCost stats are calculated
0053     %     nFlip % number of flips so far in the block
0054     %     deltaCostLog % log of deltaCost values for delta cost distribution analysis
0055     %     oldNumIt % number of iterations requested on the last optimization.  This is the default number for future optimizations.
0056     %     plotObj % sos object responsible for plotting
0057     %     histObj % object that contains detailed optimization history data
0058     %
0059     % METHODS:
0060     %   sos(varargin) -  Constructor
0061     %   addSample(sample) adds a sample to the SOS object so that its items can be optimized.
0062     %   present = containsSample(sample) evaluates whether the sos object already contains the specified sample
0063     %   constraint = addConstraint(varargin) % Adds a constraint to be optimized and returns it
0064     %   initFillSamples() % fills all samples to their specified capacity
0065     %   validItem = checkHardConstraintsFilling(sample,sItemIndex,newItemDataFrame,newItemIndex) returns a flag indicating if newIem can be added to sample during the initial filling of samples.
0066     %   swCost = checkHardConstraintsSwapping(sample,sItemIndex,newItemDataFrame,newItemIndex) returns the hard constraint cost of swapping item from newItemDataframe with item in sample.
0067     %   swCost = checkSoftConstraintsSwapping(targSample,targSampleIndex,feederdf, feederdfIndex)  % calculates the cost of a swap on the soft constraints
0068     %   swCost = checkMetaConstraintsSwapping(~,~,~,~)  % calcualtes the cost of a swap on the meta constraints
0069     %   [colName,colMean,colStd,allData] = normalizeData() % Normalizes the data in the samples and populations linked to the SOS object to provide some coarse balancing of initial cost values.
0070     %   cost = initCost()  %initializes and reports hard/soft/meta cost contraints
0071     %   setpSwapFunction(obj,swFunctionName) %sets the pSwap function for the sos obj.
0072     %   setSampleCandidateSelectionMethod(methodName) % sets the sample Candidate selection method
0073     %   setFeederdfCandidateSelectionMethod(methodName)  % sets the feeder candidate selection method
0074     %   setAnnealSchedule(varargin) % sets the anneal schedule for the SOS object
0075     %   goodEnding = optimize()  % Optimizes samples based on specified constraints
0076     %   doReport(tStart)  % displays a progress report from the last iteration in optimization
0077     %   writeSamples() % writes the data from the samples to text files specified in their 'outFile' property
0078     %   writePopulations() % writes the data from the populations to text files specified in their 'outFile' property
0079     %   writeAll() write all populations and samples associated with SOS object to disk
0080     %   addttest(obj, varargin) %adds a t-test to the list of analyses to run
0081     %   addztest(varargin) %adds a z-test to the list of analyses to run
0082     %   pass = doStatTests(varargin) % runs the stat tests and indicates if all passed the user hypotheses.
0083     %   present = containsConstraint(constraint) % determines whether <constraint> is a constraint associated with the current SOSobj
0084     %   cost = dispCost()  %reports hard/soft/meta cost contraints
0085     %   deltaCostPercentiles() % displays the breakdown of deltaCost values per decile
0086     %   createPlots(dispIt) % creates plots for several optimization parameters
0087     %   updatePlots(curIt,cost,deltaCost,pFlip,temp) %updates the contents of the plot
0088     %   createHistory() % creates plots for several optimization parameters
0089     %   updateHistory(curIt,cost,deltaCost,pFlip,temp) %updates the contents of the plot
0090     %   writeHistory(fileName) %writes all history saved to date to file 'fileName'
0091     %   setBufferedHistoryOutfile(outFile) % writes the history on-line, one update at a time, to outfile
0092     %   enableBufferedHistoryWrite() % enables writing of buffered history.
0093     %   disbleBufferedHistoryWrite() % disables writing of buffered history.
0094     %
0095     % METHODS (STATIC,Acess = private)
0096     %   p = parseConstructorArgs(varargin) - parses the sos constructor arguments
0097     
0098     %% PROPERTIES
0099     properties 
0100         samples % array of samples linked to SOS object
0101         hardConstraints % array of hard constraint objects
0102         softConstraints % array of soft constraint objects
0103         metaConstraints % array of meta constraint objects
0104         cost % current cost
0105         sCandObj % object used to identify the target item in a sample which may be swapped
0106         feederCandObj % object used to identify the item in a population/other sample to be swapped with the target sample item
0107         pSwapObj  % object which determines the probability of making a swap based on deltaCost
0108         annealObj % object used to anneal temperature in pSwap
0109         maxIt % maximum number of iterations that an optimization will run for
0110         allData % merge of all data from all samples/ pops linked to sosobj
0111         allDataColName % name of columns in allData
0112         allDataColMean % mean for each column in allData
0113         allDataColStd  % stdev for each column in all data
0114         reportInterval % number of iterations between progress reports
0115         curIt % current optimization iteration
0116         curFreezeIt % number of iterations cost has been frozen at a particular value
0117         stopFreezeIt % number of frozen iterations before algorithm stops
0118         sosstattests % array of handles to sos  stat tests requested by user
0119         statInterval % number of iterations between stat reports
0120         statTestReportStyle % style of report ('short' or 'full')
0121         deltaCost % delta cost from soft constraints from last attempted flip that met all hard constraints
0122         blockSize % number of iterations in a 'block' of trials for which average flip history /minDeltaCost stats are calculated
0123         nFlip % number of flips so far in the block
0124         deltaCostLog % log of deltaCost values for delta cost distribution analysis
0125         oldNumIt % number of iterations requested on the last optimization.  This is the default number for future optimizations.
0126         plotObj % object responsible for plotting
0127         histObj % object that contains detailed optimization history data
0128         targSampleCandSelectMethod % name of sample replacement method
0129         feederdfCandSelectMethod % name of neighbor method
0130         queryStopOptimize % how often to probe the gui to see if the stop button was pressed
0131     end
0132     
0133     %% PROPERTIES (Constant)
0134     properties (Constant)
0135         maxIt_def = 10000; % default
0136         pSwapFunction_def = 'logistic'; % default
0137         targSampleCandSelectMethod_def = 'random'; % default
0138         feederdfCandSelectMethod_def = 'randomPopulationAndSample'; % default
0139         reportInterval_def = 1000; % default
0140         stopFreezeIt_def = 10000; % default
0141         statInterval_def = 10000; % default
0142         statTestReportStyle_def = 'short'; % default
0143         blockSize_def = 1000; % default
0144         queryStopOptimize_def = 100; % default
0145     end
0146          
0147     methods
0148         
0149         %% sos CONSTRUCTOR
0150         function obj = sos(varargin)
0151             % Constructor
0152             %
0153             % SYNOPSIS:
0154             % Creates an SOS object
0155             %
0156             %PARAMETERS:
0157             %Optional:
0158             %   'maxIt'/integer - maximum number of iterations to run the optimizer
0159             %   'pSwapFunction'/string - name of pSwapFunction to use.
0160             %   'targSampleCandSelectMethod'/string - name of target Sample Candidate selection method
0161             %   'feederdfCandSelectMethod'/string - name of feeder dataframe candidate selection method
0162             %   'reportInterval'/int - number of iterations between general cost reports
0163             %   'stopFreezeIt'/int - operationalized number of sequential iterations cost value must remain the same for state to be considered 'frozen'
0164             %   'statInterval'/int - number of iterations between stat reports
0165             %   'statTestReportStyle'/string - style of stat reports ('short' or 'full')
0166             %   'blockSize'/int - number of iterations in a block; used to determine length of deltaCostLog
0167             
0168             verbosePrint([char(10) 'Creating sos Object'], ...
0169                 'sos_constructor_startObjCreation');
0170 
0171             p = sos.parseConstructorArgs(varargin);
0172             
0173             %potentially user-specified variables
0174             obj.maxIt = p.Results.maxIt;
0175             obj.oldNumIt = obj.maxIt;
0176             obj.reportInterval = p.Results.reportInterval;
0177             obj.stopFreezeIt = p.Results.stopFreezeIt; 
0178             obj.statInterval = p.Results.statInterval;
0179             obj.statTestReportStyle = p.Results.statTestReportStyle;
0180             obj.blockSize = p.Results.blockSize;
0181             
0182             obj.queryStopOptimize = obj.queryStopOptimize_def;
0183             
0184             setpSwapFunction(obj,p.Results.pSwapFunction);
0185             setSampleCandidateSelectionMethod(obj,...
0186                 p.Results.targSampleCandSelectMethod);
0187             setFeederdfCandidateSelectionMethod(obj,...
0188                 p.Results.feederdfCandSelectMethod);
0189             
0190             obj.targSampleCandSelectMethod = p.Results.targSampleCandSelectMethod;
0191             obj.feederdfCandSelectMethod = p.Results.feederdfCandSelectMethod;
0192             
0193             
0194             obj.curIt = 0;            
0195             obj.samples = [];       
0196             obj.hardConstraints = [];
0197             obj.softConstraints = [];
0198             obj.metaConstraints = [];
0199             obj.sosstattests = [];
0200             obj.curFreezeIt = 0;
0201             obj.deltaCost = NaN;
0202             obj.nFlip = NaN;
0203             obj.deltaCostLog = nan(obj.blockSize,1);
0204             
0205             %set default anneal schedule to greedy
0206             obj.setAnnealSchedule();
0207 
0208             verbosePrint(['Creation of sos object complete' char(10)], ...
0209                 'sos_constructor_endObjCreation');
0210         end
0211 
0212         
0213         %% addSample METHOD
0214         function obj = addSample(obj, sample)
0215             % adds a sample to the SOS object so that its items can be optimized.
0216             %
0217             % CALL:
0218             % <sosObj>.addSample(<sampleObj>)
0219             %
0220             % PARAMETERS:
0221             % sample - a sample object
0222             
0223             p = inputParser;
0224             p.addRequired('obj');
0225             p.addRequired('sample', ...
0226                  @(sample)strcmp(class(sample),'sample'));
0227             p.parse(obj,sample);
0228                                             
0229             if(isempty(sample.population)==1)
0230                 verbosePrint('Warning: Sample has not been linked with a population', ...
0231                     'sos_addSample_NoPopForSampleWarn');
0232             end
0233           
0234             present = obj.containsSample(sample);
0235             
0236             if(present == 1)
0237                 verbosePrint('Warning: Sample has already been added.  It cannot be added again', ...
0238                     'sos_addSample_SampleAlreadyAddedWarn');              
0239             else
0240                 obj.samples = [obj.samples sample];  
0241                 verbosePrint('Adding new Sample...', ...
0242                     'sos_addSample_sampleAdded');
0243             end
0244             
0245         end %addSample
0246         
0247         
0248         %% present = containsSample(sample) METHOD
0249         function present = containsSample(obj, sample)
0250             % evaluates whether the sos object already contains the specified sample
0251             %
0252             % PARAMETERS:
0253             %   sample - a sample object
0254             
0255             p = inputParser;
0256             p.addRequired('obj');
0257             p.addRequired('sample', ...
0258                  @(sample)strcmp(class(sample),'sample'));
0259             p.parse(obj,sample);
0260                                             
0261             present = max(ismember(obj.samples,sample));
0262  
0263         end %addSample
0264         
0265         
0266         %% addConstraint(varargin) METHOD
0267         function constraint = addConstraint(obj, varargin)
0268            % Adds a constraint to be optimized
0269            %
0270            % PARAMETERS:
0271            % varies dependig on object to be created.  See
0272            % genericConstraint and the specific constraint you wish to
0273            % create for options.
0274            
0275            constraint = genericConstraint.createConstraint('sosObj',obj,varargin{:});
0276            
0277            if(strcmp(constraint.constraintType,'hard'))
0278                obj.hardConstraints = [obj.hardConstraints {constraint}];
0279            elseif(strcmp(constraint.constraintType,'soft'))
0280                obj.softConstraints = [obj.softConstraints {constraint}];          
0281            elseif(strcmp(constraint.constraintType,'meta'))
0282                obj.metaConstraints = [obj.metaConstraints {constraint}];
0283            else
0284                error('The type of the new constraint is not supported by the sos object.  Supported types are hard/soft/meta');
0285            end
0286          
0287         end
0288 
0289         %% present = containsConstraint(constraint) METHOD
0290         function present = containsConstraint(obj, constraint)
0291             % determines whether <constraint> is a constraint associated with the current SOSobj
0292             %
0293             % Returns 1 if <constraint> is present, 0 otherwise
0294             
0295             present = 0;
0296 
0297             for i=1:length(obj.hardConstraints)
0298                 if(obj.hardConstraints{i} == constraint)
0299                     present = 1;
0300                     return;
0301                 end
0302             end
0303             
0304             for i=1:length(obj.softConstraints)
0305                 if(obj.softConstraints{i} == constraint)
0306                     present = 1;
0307                     return;
0308                 end
0309             end
0310             
0311             for i=1:length(obj.metaConstraints)
0312                 if(obj.metaConstraints{i} == constraint)
0313                     present = 1;
0314                     return;
0315                 end
0316             end
0317             
0318         end
0319         
0320         %% initFillSamples()  METHOD
0321         function obj = initFillSamples(obj)
0322             % fills all samples to their specified capacity
0323             %
0324             
0325             verbosePrint('Filling All Samples...','sos_initFillSamples_start');
0326             
0327             if(isempty(obj.samples))
0328                  error('Error: No Samples in SOS Object!');
0329             end
0330             
0331             %total add stores references to the number of items left to
0332             %fill in to all of the different samples.
0333             totalAdd = [];
0334             for i=1:length(obj.samples)
0335                %find out how many observations are needed,
0336                %how many have already been filled,
0337                %and then fill the remainder
0338                
0339                n = obj.samples(i).n;
0340                
0341                if(isempty(obj.samples(i).data))
0342                    curN= 0;
0343                else
0344                 curN = length(obj.samples(i).data{1});
0345                end
0346                
0347                toAdd = n-curN;               
0348                totalAdd = [totalAdd ones(1,toAdd)*i]; %#ok<AGROW>
0349                    
0350             end
0351             
0352            %now add all of those observations from the appropriate
0353            %population, enforcing hard bounds
0354 
0355            while(isempty(totalAdd) ==0)
0356                %select a sample at random
0357                toAddIndex=floor((length(totalAdd)*rand)+1);
0358                sIndex= totalAdd(toAddIndex);
0359 
0360                %randomly select an item from it's population:
0361                 %make sure that the target population has at least one
0362                 %item still in it
0363                 if(isempty(obj.samples(sIndex).population))
0364                     error(['Sample ''',obj.samples(sIndex).name,''' is not filled to capacity, but has not been linked to a population to fill it']);
0365                 end
0366                         
0367                     
0368                 if(isempty(obj.samples(sIndex).population.data)==1)
0369                     error(['sos.initFillSamples Failed because ', ...
0370                         'a population was depeleted of items but ',...
0371                         'more items were still needed']);                       
0372                 else
0373                     pItemIndex = ...
0374                      floor(length(obj.samples(sIndex).population.data{1})*rand+1);
0375                 end
0376    
0377                 if(isempty(obj.samples(sIndex).data))
0378                     sItemIndex = 1;
0379                 else
0380                     sItemIndex = length(obj.samples(sIndex).data{1});
0381                 end
0382                 
0383                 %pass it the current sample, where the item will go in it,
0384                 %then the population, and the item from the pop's index.
0385                 validItem = obj.checkHardConstraintsFilling( ...
0386                     obj.samples(sIndex),sItemIndex, ...
0387                     obj.samples(sIndex).population,pItemIndex);
0388                               
0389                 if(validItem == 1)
0390                     %add that population item to the sample list; remove it
0391                     %from the population.
0392                     totalAdd(toAddIndex)=[]; %#ok<AGROW>
0393                     pItem = obj.samples(sIndex).population.popItem(pItemIndex);
0394                     obj.samples(sIndex).appendItem(pItem);
0395                 end
0396            end            
0397         end % initFillSamples
0398             
0399         
0400         %% validItem = checkHardConstraintsFilling(obj,sample,sItemIndex,newItemDataFrame,newItemIndex)
0401         function validItem = checkHardConstraintsFilling(obj,sample,sItemIndex,newItemDataFrame,newItemIndex)
0402             % returns a flag indicating if newItem can be added to sample during the initial filling of samples.
0403             %
0404             % Should only be invoked when initially filling the samples
0405             %
0406             %
0407             %   sample - the target sample
0408             %   sItemIndex - row index of item to swap
0409             %   newItemDataframe - dataframe (sample/pop) containin the item to fill with.
0410             %   newItemIndex - row index of item to swap.
0411             validItem = 1;
0412             
0413             for i=1:length(obj.hardConstraints)
0414                 if(obj.hardConstraints{i}.s1 == sample)
0415                     if (obj.hardConstraints{i}.itemCostFilling(...
0416                             sample,sItemIndex, ...
0417                             newItemDataFrame,newItemIndex) == 1)   
0418                         
0419                         validItem = 0;
0420                     end
0421                 end
0422             end        
0423         end % checkHardConstraintsFilling
0424         
0425         
0426         %% swCost = checkHardConstraintsSwapping(obj,sample,sItemIndex,newItemDataFrame,newItemIndex)
0427         function swCost = checkHardConstraintsSwapping(obj,sample,sItemIndex,newItemDataFrame,newItemIndex)
0428             % returns the hard constraint cost of swapping item from newItemDataframe with item in sample.
0429             %
0430             % Generally, if swCost > 0 the swap should not be executed
0431             %   sample - the target sample
0432             %   sItemIndex - row index of item to swap
0433             %   newItemDataframe - dataframe (sample/pop) containin the item to fill with.
0434             %   newItemIndex - row index of item to swap.
0435             curCost = 0;
0436             newCost = 0;
0437             for i=1:length(obj.hardConstraints)
0438                 if(obj.hardConstraints{i}.s1 == sample)
0439                         
0440                     curCost = curCost + obj.hardConstraints{i}.cost;
0441                     newCost = newCost + obj.hardConstraints{i}.swapCost(...
0442                         sample,sItemIndex,newItemDataFrame,newItemIndex);
0443                 elseif (obj.hardConstraints{i}.s1 == newItemDataFrame)
0444                     curCost = curCost + obj.hardConstraints{i}.cost;
0445                     newCost = newCost + obj.hardConstraints{i}.swapCost(...
0446                         newItemDataFrame,newItemIndex,sample,sItemIndex);
0447                     
0448                 end        
0449             end
0450                        
0451             swCost = newCost - curCost ;
0452                         
0453         end
0454         
0455         %% swCost =  checkSoftConstraintsSwapping(targSample,targSampleIndex,feederdf, feederdfIndex)
0456         function swCost = checkSoftConstraintsSwapping(obj,targSample,targSampleIndex,feederdf, feederdfIndex)
0457             % calculates the cost of a swap on the soft constraints
0458             %
0459             % PARAMETERS:
0460             %   targSample - the target sample
0461             %   targSampleIndex - row index of item to swap
0462             %   feederdf - dataframe (sample/pop) containin the item to fill with.
0463             %   feederdfIndex - row index of item to swap.
0464             %
0465             % RETURNS:
0466             % swCost - cost of making the swap
0467             
0468             swCost = 0;
0469             
0470             for j=1:length(obj.softConstraints)
0471                 if ((obj.softConstraints{j}.s1 == targSample ||  ...
0472                         obj.softConstraints{j}.s2 == targSample || ...
0473                         obj.softConstraints{j}.s1 == feederdf || ...
0474                         obj.softConstraints{j}.s2 == feederdf ))
0475                     swCost = swCost + ...
0476                       obj.softConstraints{j}.swapCost(...
0477                        targSample,targSampleIndex,feederdf, feederdfIndex);
0478                 else
0479                     swCost = swCost +obj.softConstraints{j}.cost;
0480                 end
0481             end       
0482         end
0483          
0484         %% swCost = checkMetaConstraintsSwapping(~,~,~,~) METHOD
0485         function swCost = checkMetaConstraintsSwapping(obj,~,~,~,~)
0486             % calcualtes the cost of a swap on the meta constraints
0487             %
0488             %Parameters:
0489             % parameters are the same as for hard and soft constraints, but
0490             % are ignored.  They are accepted as arguments merely to
0491             % increase the consistency of the argument.
0492             
0493             %checking meta constraints is extremely computationally cheap,
0494             %so they are calculated every time.
0495             
0496            swCost = 0; 
0497            
0498            for j=1:length(obj.metaConstraints)
0499                swCost = swCost + obj.metaConstraints{j}.swapCost();
0500            end
0501         end
0502              
0503              
0504         %% [colName,colMean,colStd,allData] = normalizeData() METHOD
0505         function [colName,colMean,colStd,allData] = normalizeData(obj)
0506             % Normalizes the data in the samples and populations linked to the SOS object to provide some coarse balancing of initial cost values.
0507             
0508             verbosePrint('Normalizing Population and Sample Data', ...
0509                         'sos_normalizeData_start');
0510          
0511             %empty data frame
0512             allData = dataFrame();
0513             
0514             alreadyMerged = {};
0515             
0516             %go through each of the samples, and their populations.
0517             for i=1:length(obj.samples)
0518                 %first, check the population:
0519                 
0520                 popPresent = 0;
0521                 if(isempty(alreadyMerged))
0522                     popPresent = 0;
0523                 else
0524                     for j=1:length(alreadyMerged)
0525                         if(alreadyMerged{j} == obj.samples(i).population)
0526                             popPresent = 1;
0527                         end
0528                     end
0529                 end
0530                 
0531                 %population has not already been added, so add it
0532                 if(popPresent == 0)
0533                     %only add the population if the sample does in fact
0534                     %have a population
0535                     if isempty(obj.samples(i).population) ~= 1
0536                         allData = dataFrame.aContainsb(allData,obj.samples(i).population);
0537                         allData = dataFrame.aContainsbData(allData,obj.samples(i).population);
0538                         alreadyMerged = [alreadyMerged {obj.samples(i).population}]; %#ok<AGROW>
0539                     end
0540                 end
0541                 
0542                 %now add the sample
0543                 allData = dataFrame.aContainsb(allData,obj.samples(i));
0544                 allData = dataFrame.aContainsbData(allData,obj.samples(i));
0545                 %we don't need to check that the sample hasn't been added
0546                 %stricly speaking, because samples can only appear once.
0547                 %However, this may be necessary in a future version of
0548                 %the code...
0549                 alreadyMerged = [alreadyMerged {obj.samples(i)}]; %#ok<AGROW>
0550                
0551                 
0552             end
0553             
0554             %we now have all of the data in a single dataframe; each column
0555             %can now be normalized and descriptive statistics stored.
0556             
0557             colName = cell(1,length(allData.data));
0558             colMean = NaN(1,length(allData.data));
0559             colStd =  NaN(1,length(allData.data));
0560             
0561             
0562             for i=1:length(allData.data)
0563                 %if it's a numeric vector, calculate mean/stdev
0564                 colName{i} = allData.header{i}; 
0565                 
0566                 if(strcmp(allData.format{i},'%f'))
0567                    colMean(i) = nanmean(allData.data{i});
0568                    colStd(i) = nanstd(allData.data{i});
0569                 end    
0570             end
0571                            
0572             %configure the z-data field in each sample and population to
0573             %reflect the new changes.
0574             
0575             for i=1:length(alreadyMerged)
0576                %link the current normalization to this object
0577                 alreadyMerged{i}.sosObj = obj; %#ok<AGROW>
0578                 alreadyMerged{i}.zdata = {}; %#ok<AGROW>
0579                
0580                for j=1:length(alreadyMerged{i}.data)
0581                    col = alreadyMerged{i}.data{j};
0582                    
0583                   if(strcmp(alreadyMerged{i}.format{j},'%f'))
0584                       m= NaN;
0585                       std = NaN;
0586                       %get the relevant coefficients:
0587                       for k=1:length(colName)
0588                          if (strcmp(colName{k},alreadyMerged{i}.header{j}))                           
0589                             m = colMean(k);
0590                             std = colStd(k);
0591                              break;
0592                          end
0593                       end
0594                       
0595                       %check that m and std exist
0596                       if(isnan(m) || isnan(std))
0597                           error('Normalization failed because either the mean or stdeviation was NaN');
0598                       end
0599                       
0600                       %they exist, so normalization can proceed
0601                       
0602                       
0603                       if std == 0
0604                           verbosePrint(['Warning: No variance in column: ', alreadyMerged{i}.header{j}],...
0605                             'sos_normalizeData_noVarWarn');
0606                           std = 1;
0607                       end
0608                       
0609                       col = (col-m)/std;
0610                                           
0611                   end
0612 
0613                   alreadyMerged{i}.zdata{j} = col; %#ok<AGROW>
0614                           
0615                end                
0616             end    
0617             
0618             obj.allData = allData;
0619             obj.allDataColName = colName;
0620             obj.allDataColMean = colMean;
0621             obj.allDataColStd = colStd;
0622             
0623         end %normalizeData
0624         
0625         
0626         %% initCost() METHOD
0627         function cost = initCost(obj)
0628             %initializes and reports hard/soft/meta cost contraints
0629 
0630             verbosePrint([char(10),'Initializing Constraints',char(10)],...
0631                 'sos_initCost_startInit');
0632             
0633             
0634             
0635             % hard constraints
0636             hardConstraintCost = 0;
0637             for j = 1:length(obj.hardConstraints)
0638                 try
0639                     tempCost = obj.hardConstraints{j}.initCost();
0640                 catch exception
0641                     if strcmp(exception.identifier, ...
0642                             'MATLAB:cellRefFromNonCell')
0643                         error(['ERROR: Data needed to calculate hard constraint cost #' num2str(j),' is missing.',...
0644                             char(10),'       Have samples been filled and normalized?']);
0645                   
0646                     else
0647                         throw exception
0648                     end
0649                     
0650                 end
0651                 
0652                 verbosePrint(['   Hard Constraint #: ',num2str(j),...
0653                     ' Cost: ' num2str(tempCost),'   ',obj.hardConstraints{j}.name], ...
0654                         'sos_initCost_indivHardCost');
0655                 hardConstraintCost = hardConstraintCost + tempCost;
0656             end
0657  
0658             verbosePrint(['Hard Constraint Total: ', num2str(hardConstraintCost)], ...
0659                         'sos_initCost_totalHardCost');
0660                     
0661             % soft constraints
0662             softConstraintCost = 0;
0663             for j=1:length(obj.softConstraints)
0664                 try
0665                     tempCost = obj.softConstraints{j}.initCost();    
0666                 catch exception
0667                     if strcmp(exception.identifier, ...
0668                             'MATLAB:cellRefFromNonCell')
0669                         error(['ERROR: Data needed to calculate soft constraint cost #' num2str(j),' is missing.',...
0670                             char(10),'       Have samples been filled and normalized?']);
0671                   
0672                     else
0673                         throw exception
0674                     end
0675                     
0676                 end
0677                 
0678                 verbosePrint(['   Soft Constraint #: ',num2str(j),...
0679                     ' Cost: ' num2str(tempCost),'   ',obj.softConstraints{j}.name], ...
0680                         'sos_initCost_indivSoftCost');
0681                softConstraintCost = softConstraintCost + tempCost; 
0682             end
0683             
0684             verbosePrint(['Soft Constraint Total: ', num2str(softConstraintCost)], ...
0685                         'sos_initCost_totalSoftCost');
0686             
0687             % meta constraints
0688             metaConstraintCost = 0;
0689             for j=1:length(obj.metaConstraints)
0690                 try
0691                     tempCost = obj.metaConstraints{j}.initCost();    
0692                 catch exception
0693                     if strcmp(exception.identifier, ...
0694                             'MATLAB:cellRefFromNonCell')
0695                         error(['ERROR: Data needed to calculate meta constraint cost #' num2str(j),' is missing.',...
0696                             char(10),'       Have samples been filled and normalized?']);
0697                   
0698                     else
0699                         throw exception
0700                     end
0701                     
0702                 end
0703                 
0704                 verbosePrint(['   Meta Constraint #: ',num2str(j),...
0705                     ' Cost: ' num2str(tempCost),'   ',obj.metaConstraints{j}.name], ...
0706                         'sos_initCost_indivMetaCost');
0707                metaConstraintCost = metaConstraintCost + tempCost; 
0708             end
0709             
0710             verbosePrint(['Meta Constraint Total: ', num2str(metaConstraintCost)], ...
0711                         'sos_initCost_totalMetaCost');            
0712             
0713             
0714             cost = softConstraintCost + metaConstraintCost;
0715             obj.cost=cost;
0716             
0717             verbosePrint([char(10),'TOTAL COST (soft + meta): ', num2str(cost), char(10)], ...
0718                 'sos_initCost_totalSoftMetaCost');             
0719         end
0720         
0721         %% dispCost() METHOD
0722         function cost = dispCost(obj)
0723             %reports hard/soft/meta cost contraints
0724 
0725             verbosePrint([char(10),'Displaying Constraint Costs',char(10)],...
0726                 'sos_dispCost_start');
0727             
0728             % hard constraints
0729             hardConstraintCost = 0;
0730 
0731             for j = 1:length(obj.hardConstraints)
0732                 tempCost = obj.hardConstraints{j}.cost;    
0733                 verbosePrint(['   Hard Constraint #: ',num2str(j),...
0734                     ' Cost: ' num2str(tempCost),'   ',obj.hardConstraints{j}.name], ...
0735                         'sos_dispCost_indivHardCost');
0736                 hardConstraintCost = hardConstraintCost + tempCost;
0737             end
0738             
0739             verbosePrint(['Hard Constraint Total: ', num2str(hardConstraintCost)], ...
0740                         'sos_dispCost_totalHardCost');
0741                     
0742             % soft constraints
0743             softConstraintCost = 0;
0744             
0745            
0746             for j=1:length(obj.softConstraints)
0747                 tempCost = obj.softConstraints{j}.cost;    
0748                 verbosePrint(['   Soft Constraint #: ',num2str(j),...
0749                     ' Cost: ' num2str(tempCost),'   ',obj.softConstraints{j}.name], ...
0750                         'sos_dispCost_indivSoftCost');
0751                softConstraintCost = softConstraintCost + tempCost; 
0752             end
0753 
0754             
0755             
0756             verbosePrint(['Soft Constraint Total: ', num2str(softConstraintCost)], ...
0757                         'sos_dispCost_totalSoftCost');
0758             
0759             % meta constraints
0760             metaConstraintCost = 0;
0761             
0762            
0763             for j=1:length(obj.metaConstraints)
0764                 tempCost = obj.metaConstraints{j}.cost;    
0765                 verbosePrint(['   Meta Constraint #: ',num2str(j),...
0766                     ' Cost: ' num2str(tempCost),'   ',obj.metaConstraints{j}.name], ...
0767                         'sos_dispCost_indivMetaCost');
0768                metaConstraintCost = metaConstraintCost + tempCost; 
0769             end
0770 
0771             
0772             verbosePrint(['Meta Constraint Total: ', num2str(metaConstraintCost)], ...
0773                         'sos_dispCost_totalMetaCost');            
0774             
0775             
0776             cost = softConstraintCost + metaConstraintCost;
0777             %obj.cost=cost;
0778             
0779             verbosePrint([char(10),'TOTAL COST (soft + meta): ', num2str(cost), char(10)], ...
0780                 'sos_dispCost_totalSoftMetaCost'); 
0781             
0782             if isnan(cost)
0783                 verbosePrint([char(10), ...
0784                     'Warning: Cost was NaN.  This should not happen if cost terms ',...
0785                     char(10),'         are present and cost has been initialized',...
0786                     char(10)], 'sos_dispCost_warnCostNaN');
0787             end
0788             
0789         end %dispCost
0790 
0791         %% setAnnealSchedule(varargin) METHOD
0792         function setAnnealSchedule(obj,varargin)
0793             % sets the anneal schedule for the SOS object
0794             %
0795             % PARAMETERS:
0796             % Required:
0797             % 'schedule'/scheduleName - param/value pair indicating the name of the anneal schedule to use.  Defaults to greedy'
0798             %
0799             % Optional:
0800             %  - As required by specific schedule to create.  See its
0801             %  constructor for details.
0802             
0803             p = inputParser;
0804             p.addRequired('obj', @(obj)strcmp(class(obj),'sos'));
0805             p.addParamValue('schedule','greedy', ...
0806                 @(schedule)ischar(schedule));            
0807             p.KeepUnmatched = true;
0808             
0809             % in case no parameters are passed to varargin, must make the
0810             % following decision
0811             if(isempty(varargin) == false)
0812                  p.parse(obj,varargin{:});
0813             else
0814                 p.parse(obj);
0815             end
0816             
0817             scheduleName = p.Results.schedule;
0818              
0819             if any(strcmp(p.UsingDefaults,'schedule'))
0820                 verbosePrint(['    Defaulting to ', scheduleName, ' annealing...'], ...
0821                     'sos_setAnnealSchedule_defaultAnneal'); 
0822             end
0823             
0824            if(strcmp(scheduleName, 'greedy') == 1)
0825                obj.annealObj = greedyAnneal(varargin{:});
0826            elseif (strcmp(scheduleName, 'exp') == 1)
0827                 obj.annealObj = expAnneal(varargin{:});
0828            else
0829                error(['There is no schedule named: ' schedule])
0830            end                 
0831         end
0832             
0833         
0834         %% setSampleCandidateSelectionMethod(methodName) METHOD
0835         function setSampleCandidateSelectionMethod(obj,methodName)
0836             % sets the sample Candidate selection method
0837             %
0838             %PARAMETERS:
0839             % Name of a sampleCandidateSelectionMethod
0840             
0841             if (ischar(methodName) == false) 
0842                 error('{methodName} must be a string');
0843             end   
0844             
0845            if (strcmp(methodName, 'random') == 1)
0846                 obj.sCandObj = randSampleCandidateSelection(obj);
0847            else
0848                error(['Specified methodName ',methodName,' does not exist']); 
0849            end
0850         end
0851         
0852         
0853         %% setFeederdfCandidateSelectionMethod(obj,methodName)
0854         function setFeederdfCandidateSelectionMethod(obj,methodName)
0855             % sets the feeder candidate selection method
0856             %
0857             %PARAMETERS:
0858             % name of a feederCandidateSelectionMethod
0859             
0860             if (ischar(methodName) == false) 
0861                 error('{methodName} must be a string');
0862             end   
0863             
0864            if (strcmp(methodName, 'randomPopulation') == 1)
0865                 obj.feederCandObj= randPopulationCandidateSelection();
0866                 obj.feederdfCandSelectMethod = methodName;
0867            elseif (strcmp(methodName, 'randomPopulationAndSample') == 1)
0868                 obj.feederCandObj= randPopulationAndSampleCandidateSelection(obj);
0869                 obj.feederdfCandSelectMethod = methodName;
0870            else
0871                 error('Specified {methodName} does not exist'); 
0872            end
0873             
0874         end
0875         
0876          %% setpSwapFunction(methodName) METHOD
0877          function setpSwapFunction(obj,swFunctionName)
0878              %sets the pSwap function for the sos obj.
0879              %
0880              %PARAMETERS:
0881              %  swFunctionName - name of swap function.  Currently, 'logicistic' is the only supported function.
0882 
0883             if (ischar(swFunctionName) == false) 
0884                 error('SOS:setSwapFunction {swFunctionName} is not a valid string');
0885             end   
0886             
0887            if (strcmp(swFunctionName, 'logistic'))
0888                 obj.pSwapObj = logisticpSwapFunction();
0889            else
0890                 error('SOS:setpSwapFunction: specified function does not exist');
0891            end            
0892         end %setpSwapFunction
0893         
0894         
0895         %% optimize() METHOD
0896         function goodEnding = optimize(obj,varargin)
0897             % Optimizes samples based on specified constraints
0898             %
0899             % PARAMETERS
0900             %   numIt - number of iterations to run (otherwises uses
0901             %           default from obj init)
0902             %   'isGui' / 1 or 0 - flag denoting whether there is an active gui or not
0903             % RETURNS:
0904             %   1 if ended because statistical constraints reached, zero
0905             %   otherwise
0906             
0907             % It would be possible to seperate the normalization and
0908             % costInitialization procedures so that they need not be run
0909             % every time you begin optimizing.  In fact, the current code
0910             % probably enforces such a constraint and should generate an
0911             % error if normalization is attempted anyways.  Nevertheless,
0912             % it costs little to enforce it at this point in time given
0913             % that most time should be spent trying to make swaps, so it's
0914             % left that way in this version of the code.  In the future, it
0915             % may be useful to avoid any possible redundancy here, and also
0916             % consider adding the fillSamples() command to the
0917             % pre-optimization procedure.
0918                
0919             verbosePrint(['Setting up optimization procedure...', ...
0920                         char(10)], 'sos_optimize_begin');    
0921             
0922             
0923             %get all parameters passed to the alogirthm
0924             
0925             p = inputParser;
0926             p.addOptional('numIt',-1);
0927             p.addParamValue('isGui',NaN);
0928             p.parse(varargin{:});
0929             
0930             if any(strcmp(p.UsingDefaults,'numIt'))
0931                 % do nothing, allow defaults to be used
0932             else
0933                 numIt = p.Results.numIt;
0934                 validateattributes(numIt, {'numeric'}, ...
0935                     {'scalar', 'integer', 'positive', '>', 0})
0936                 obj.maxIt = obj.curIt + numIt; 
0937                 obj.oldNumIt = numIt;
0938             end
0939             
0940             if any(strcmp(p.UsingDefaults,'isGui'))
0941                 % do nothing, no gui
0942                 isGui = 0;
0943             else
0944                 isGui = p.Results.isGui;
0945                 
0946                 if isGui == 1
0947                     mainWindowHandle = sos_gui; %this connects/ launches the gui
0948                     mainWindowData = guidata(mainWindowHandle);
0949                     set(mainWindowData.pushbutton_stopOptimize,'UserData',0)
0950                 elseif isGui == 0
0951                     % do nothing
0952                 else
0953                     error('isGui should be a flag with value 1 or 0).');
0954                 end
0955             end
0956             
0957              
0958             % flag indicating whether the optimization should be
0959             % stopped or not, as inidicated in the GUI
0960             stopOptimize = 0;
0961                     
0962             goodEnding = 0;
0963             
0964             % make sure that there are samples linked with the object
0965             % before running.
0966             
0967             if isempty(obj.samples) 
0968                 error('Cannot start optimization - there are no samples to optimize');
0969             end
0970             
0971            if isempty(obj.pSwapObj) 
0972                 error('Cannot start optimization - the annealing schedule must be set (e.g., greedy, exp)');
0973             end
0974             
0975             % next, make sure that all of the samples have been filled.
0976             
0977             for i=1:length(obj.samples)
0978                 if isempty(obj.samples(i).data)
0979                     error(['Cannot start optimization - sample ''',...
0980                         obj.samples(i).name,''' is empty']);
0981                 elseif length(obj.samples(i).data{1}) < obj.samples(i).n
0982                     error(['Cannot start optimization - sample ''',...
0983                         obj.samples(i).name,''' is missing items.  Did you initFillSamples()?']);
0984                 end
0985             end
0986               
0987             
0988             if(obj.curIt == obj.maxIt)
0989                 obj.maxIt = obj.maxIt+obj.oldNumIt;
0990             elseif obj.curIt > obj.maxIt
0991                 error('Cannot start Optimization - curIt should never be greater than maxIt.  Increase by <SOSobj>.maxIt = <larger val>?');
0992             end
0993                     
0994             
0995             if(obj.curFreezeIt == obj.stopFreezeIt)
0996                 obj.curFreezeIt = 0;
0997                 verbosePrint('curFreezeIt was at maximum value at start of optimization; resetting to 0', ...
0998                 'sos_optimize_resetFreezeIt');
0999             end
1000             
1001             obj.normalizeData();
1002             obj.initCost();
1003  
1004             % set up the sample and neighbor replacement methods
1005             obj.setSampleCandidateSelectionMethod(obj.targSampleCandSelectMethod);
1006             obj.setFeederdfCandidateSelectionMethod(obj.feederdfCandSelectMethod);
1007             
1008             
1009             
1010             verbosePrint(['Optimizing for ', num2str(obj.maxIt-obj.curIt), ' iterations'], ...
1011                 'sos_optimize_startOptimization');
1012             
1013             tStart = tic;
1014             startIt = obj.curIt;
1015             allStatsPass = NaN;
1016             
1017             reportHeader = [char(10), ' Iteration              Cost  %Complete       Elapsed   Remaining'];
1018             verbosePrint(reportHeader,'sos_optimize_reportHeader');
1019             
1020             while obj.curIt<obj.maxIt && ...
1021                     obj.curFreezeIt < obj.stopFreezeIt && ...
1022                     allStatsPass ~= 1 && ...
1023                     stopOptimize == 0
1024                 obj.curIt=obj.curIt+1;
1025                
1026                 if mod(obj.curIt,obj.reportInterval) == 0  || obj.curIt == 1
1027                     obj.doReport(tStart,startIt);
1028                      obj.doStatTests('reportStyle','none');
1029                     
1030                     if isempty(obj.histObj) == 0
1031                         %updatePlots
1032                         pFlip = obj.nFlip/obj.reportInterval;
1033                         obj.updateHistory(obj.curIt,obj.cost,...
1034                               obj.deltaCost,pFlip,obj.annealObj.getTemp());
1035                         obj.nFlip = 0;
1036                         
1037                         if isempty(obj.plotObj) == 0
1038                             obj.updatePlots();
1039                         end
1040                     end
1041                 end
1042                           
1043                 
1044                 %select items to swap
1045                 [targSample,targSampleIndex] = obj.sCandObj.getCandidateIndex(); 
1046                 [feederdf, feederdfIndex] = obj.feederCandObj.getCandidateIndex(targSample);
1047                 
1048 
1049                 % check to see how swap fares with regard to hard
1050                 % constraints.  Abort swap if it would lead to an
1051                 % (increased) violation of hard constraints.
1052                 sumHardswCost = obj.checkHardConstraintsSwapping(...
1053                        targSample,targSampleIndex,feederdf, feederdfIndex);
1054                 
1055                 %if hard constraint was not violated or improved continue.
1056                 %A value < 0 indicates that hard constraint was improved so
1057                 %there will be a swap no matter what the soft cost.
1058                 if sumHardswCost <= 0 
1059                    curCost = obj.cost;
1060                        
1061                     sumSoftswCost = obj.checkSoftConstraintsSwapping(...
1062                        targSample,targSampleIndex,feederdf, feederdfIndex);
1063                    
1064                     sumMetaswCost = obj.checkMetaConstraintsSwapping(...
1065                        targSample,targSampleIndex,feederdf, feederdfIndex);
1066                     
1067                     swCost = sumSoftswCost + sumMetaswCost;
1068                     %sign flipped so that if swCost < curCost, we get a
1069                     %negative value, and negative values descend in cost
1070                     %space
1071                     
1072                     deltaCost = -1*(curCost - swCost); %#ok<PROP>
1073                     
1074                     % if a swap is possible based on hard constraints, but
1075                     % not necessary, base decision to swap on soft
1076                     % constraints
1077                     if sumHardswCost == 0
1078                         shouldSwap = obj.pSwapObj.shouldSwap(...
1079                             deltaCost,obj.annealObj.getTemp()); %#ok<PROP>
1080                         
1081                         obj.deltaCost = deltaCost; %#ok<PROP>
1082                         
1083                         
1084                         
1085                     else %swap is necessary based on hard constraints.
1086                         shouldSwap = 1;
1087                         
1088                     end
1089                     
1090                     if shouldSwap == 1
1091                         
1092                         obj.nFlip = obj.nFlip + 1;
1093                         
1094                         
1095                         %we need to swap items
1096                         targSample.swapItems(...
1097                         targSampleIndex,feederdf, feederdfIndex,obj);
1098 
1099                         %update constraints
1100                         for j = 1:length(obj.hardConstraints)
1101                             obj.hardConstraints{j}.acceptSwap();
1102                         end
1103 
1104                         for j=1:length(obj.softConstraints)
1105                             obj.softConstraints{j}.acceptSwap();
1106                         end
1107 
1108                         for j=1:length(obj.metaConstraints)
1109                             obj.metaConstraints{j}.acceptSwap();
1110                         end
1111 
1112                         % see if the set is frozen
1113                         if sumHardswCost < 0
1114                             obj.curFreezeIt = 0;
1115                         elseif obj.cost ~= swCost
1116                             obj.curFreezeIt = 0;
1117                         else % we're swapping items with the same cost
1118                             obj.curFreezeIt = obj.curFreezeIt +1;
1119                         end
1120                         
1121                         obj.cost = swCost;
1122                     else
1123                         % reject the swap:
1124                         for j = 1:length(obj.hardConstraints)
1125                             obj.hardConstraints{j}.rejectSwap();
1126                         end
1127 
1128                         for j=1:length(obj.softConstraints)
1129                             obj.softConstraints{j}.rejectSwap();
1130                         end
1131 
1132                         for j=1:length(obj.metaConstraints)
1133                             obj.metaConstraints{j}.rejectSwap();
1134                         end
1135                         
1136                         obj.curFreezeIt = obj.curFreezeIt + 1;
1137                         
1138                     end
1139                 else
1140                     for j = 1:length(obj.hardConstraints)
1141                             obj.hardConstraints{j}.rejectSwap();
1142                     end
1143                     
1144                     obj.curFreezeIt = obj.curFreezeIt + 1;
1145                 end
1146                 %on some subset of the intervals, display relevant
1147                 %statistics.
1148                 
1149                 %mod indexing needed to have deltaCostLog loop after
1150                 %blockSize updates.
1151                 obj.deltaCostLog(mod(obj.curIt-1,obj.blockSize)+1) = obj.deltaCost;
1152                 
1153                         
1154                 obj.annealObj.anneal(obj.curIt,obj.cost,obj.deltaCost);
1155                 
1156                 % see if stats constraint has been met
1157                 if mod(obj.curIt,obj.statInterval) == 0
1158                     allStatsPass = ...
1159                         obj.doStatTests('reportStyle',obj.statTestReportStyle);
1160                     verbosePrint(reportHeader,'sos_optimize_reportHeader');
1161                 end
1162                 
1163                 %see if GUI interrupt has been requested
1164                 if isGui
1165                     % check to see if the stop button has been pressed
1166                     % every so often.  This should not be done too
1167                     % frequently as the drawnow event could take some time.
1168                     if mod(obj.curIt,obj.queryStopOptimize) == 0
1169                         drawnow;
1170                         if  get(mainWindowData.pushbutton_stopOptimize,'UserData') ~= 0
1171                             stopOptimize = 1;
1172                         end
1173                     end
1174                 end
1175                 
1176                 
1177             end %end iterations
1178             
1179             %avoid a redundant report
1180             if mod(obj.curIt,obj.reportInterval) ~= 0 && obj.curIt ~= 1
1181                 obj.doReport(tStart,startIt);
1182                 obj.doStatTests('reportStyle','none');
1183                 
1184                 if isempty(obj.histObj) == 0
1185                     %updatePlots
1186                     pFlip = obj.nFlip/obj.reportInterval;
1187                     obj.updateHistory(obj.curIt,obj.cost,...
1188                           obj.deltaCost,pFlip,obj.annealObj.getTemp());
1189                     obj.nFlip = 0;
1190 
1191                     if isempty(obj.plotObj) == 0
1192                         obj.updatePlots();
1193                     end
1194                 end
1195             end
1196             
1197             if mod(obj.curIt,obj.statInterval) ~= 0
1198                     allStatsPass = ...
1199                         obj.doStatTests('reportStyle',obj.statTestReportStyle);
1200             end
1201              
1202             % state the reason why optimization was stopped
1203            if (allStatsPass == 1)
1204                 verbosePrint([[char(10), 'Optimization ended after all statistical tests passed defined criteria', char(10)], ...
1205                         char(10)], 'sos_optimize_endallStatsPass');
1206                 goodEnding = 1;
1207            elseif obj.curIt == obj.maxIt
1208                 verbosePrint([char(10), 'Optimization ended after reaching the maximum number of iterations', ...
1209                         char(10)], 'sos_optimize_endMaxIt');
1210            elseif obj.curFreezeIt >= obj.stopFreezeIt
1211                 verbosePrint([char(10), 'Optimization ended after cost value was frozen for specified number of iterations', ...
1212                         char(10)], 'sos_optimize_endstopFreezeIt');
1213            elseif isGui
1214                if get(mainWindowData.pushbutton_stopOptimize,'UserData') == 1
1215                verbosePrint([char(10),'Optimization ended after user interrupt (from graphical interface)',...
1216                         char(10)],'sos_optimize_endUserGuiInterrupt');
1217                end
1218            end
1219 
1220             if isGui == 1
1221                 mainWindowHandle = sos_gui; %this connects/ launches the gui
1222                 mainWindowData = guidata(mainWindowHandle);
1223                 set(mainWindowData.pushbutton_optimize,'Enable','on')
1224                 set(mainWindowData.pushbutton_stopOptimize,'Enable','off')
1225             end
1226                     
1227         end
1228         
1229         
1230         %% doReport (obj,tstart) METHOD
1231         function doReport(obj,tStart,startIt)
1232             % displays a progress report from the last iteration in optimization
1233             %
1234             %PARAMETERS:
1235             %   tStart - tic() object indicating start of optimization
1236 
1237             percentComplete = (obj.curIt/obj.maxIt)*100;
1238             elapsedTime = toc(tStart);
1239             
1240             startPercent = (startIt/obj.maxIt)*100;
1241             
1242             remainingTime = (elapsedTime/ ...
1243                (percentComplete-startPercent))*(100-startPercent)-elapsedTime;     
1244             
1245             hElapsedTime = seconds2human(elapsedTime,'full');
1246             tempLen = length(hElapsedTime);
1247             hElapsedTime = strjust([blanks(8-tempLen) hElapsedTime],'right');
1248 
1249             hRemainingTime = seconds2human(remainingTime,'full');
1250             tempLen = length(hRemainingTime);
1251             hRemainingTime = strjust([blanks(8-tempLen) hRemainingTime],'right');
1252 
1253 
1254             report = sprintf('%9.0f) \t %15.5f \t %5.2f%% \t %s \t %s', ...
1255                     obj.curIt,obj.cost,percentComplete, ...
1256                     hElapsedTime,hRemainingTime);
1257             verbosePrint(report, 'sos_optimize_report');
1258 
1259         end
1260         
1261         %% createPlots() METHOD
1262         function createPlots(obj,dispIt)
1263             % creates plots for several optimization parameters, as
1264             % follows:
1265             %   - cost
1266             %   - deltaCost
1267             %   - sosStatTestpvals
1268             %   - temperature
1269             %   - pFlipHistory
1270             %
1271             % PARAMETERS:
1272             %   dispIt - number of iterations to show on the plot screen.
1273             %       All other datapoints plotted since the creation of the
1274             %       object will still be accessible by scrolling the x-axis
1275             
1276             
1277             if isempty(obj.histObj)
1278                 obj.createHistory();
1279             end
1280             
1281             if exist('dispIt','var') == 0
1282                 dispIt = 100000;
1283             else
1284                 validateattributes(dispIt, {'numeric'}, ...
1285                 {'scalar', 'integer', 'positive', '>', 0});
1286             end
1287                 
1288                            
1289             if isempty(obj.plotObj)            
1290                 obj.plotObj = sosPlots(obj,obj.histObj,dispIt, ...
1291                                         obj.curIt);
1292             else
1293                 verbosePrint('Warning: Plots already exist so cannot be created',...
1294                     'sos_createPlots_alreadyCreated');
1295             end            
1296         end
1297         
1298         %% updatePlots() METHOD
1299         function updatePlots(obj)
1300             %updates the contents of the plot
1301             
1302             obj.plotObj.updatePlots();
1303             
1304         end
1305 
1306         
1307         %% createHistory() METHOD
1308         function createHistory(obj)
1309             % creates plots for several optimization parameters, as
1310             % follows:
1311             %   - cost
1312             %   - deltaCost
1313             %   - sosStatTestpvals
1314             %   - temperature
1315             %   - pFlipHistory
1316             
1317             if isempty(obj.histObj)
1318             
1319                 obj.histObj = sosHistory(obj);
1320             else
1321                 verbosePrint('Warning: History already exists so cannot be created',...
1322                     'sos_createHistory_alreadyCreated');
1323             end
1324             
1325         end
1326 
1327         
1328         %% updateHistory() METHOD
1329         function updateHistory(obj,curIt,cost,deltaCost,pFlip,temp)
1330             %updates the contents of the plot
1331             
1332             % get the information about the p-vals.
1333             testNames = {};
1334             testps = [];
1335             
1336             for i=1:length(obj.sosstattests)
1337                 testNames = horzcat(testNames,obj.sosstattests{i}.label); %#ok<AGROW>
1338                 testps = horzcat(testps,obj.sosstattests{i}.lastp); %#ok<AGROW>
1339             end
1340             
1341             obj.histObj.updateHistory(curIt,cost,deltaCost,pFlip, ...
1342                             testNames,testps,temp);
1343             
1344         end
1345         
1346         %% function writeHistory(outFile) METHOD
1347         function writeHistory(obj,outFile)
1348             %writes all history saved to date to file 'outFile'
1349                 
1350              if exist('outFile','var') == 0
1351                  error('"Outfile" argument was not supplied to writeHistory()');
1352              end
1353              
1354              if(ischar(outFile) == false)
1355                 error('Outfile has not been set to a string != "null".');
1356              end
1357     
1358             if isempty(obj.histObj) == 0
1359                 obj.histObj.writeHistory(outFile);
1360             else
1361                 error('No saved history to write - did you run createHistory() before optimizing?');
1362             end
1363         end
1364         
1365         %% setbufferedHistoryOutfile(outFile) METHOD
1366         function setBufferedHistoryOutfile(obj,outFile)
1367             % writes the history on-line, one update at a time, to outfile
1368             % If outfile exists, it will be overridden.
1369             
1370             if exist('outFile','var') == 0
1371                 error('"Outfile" argument was not supplied to writeHistory()');
1372             end
1373 
1374             if(ischar(outFile) == false)
1375                 error('Outfile has not been set to a string != "null".');
1376             end
1377 
1378             if ischar(outFile) == false || strcmp(outFile,'null')
1379                 error('Outfile has not been set to a string != "null".');
1380             end
1381             
1382             if isempty(obj.histObj) == 0
1383                 obj.histObj.setBufferedHistoryOutfile(outFile);
1384             else
1385                 error('No saved history to write - did you start recording history with createHistory()?');
1386             end                                    
1387         end
1388 
1389         %% enableBufferedHistoryWrite() METHOD
1390         function enableBufferedHistoryWrite(obj)
1391             % enables writing of buffered history.  Enabled automatically
1392             % when bufferedHistoryWrite is first called; cannot be called
1393             % until that method has first been called to specify an
1394             % outfile.  Subsequent enablings continue to write to that
1395             % file.
1396             
1397             if isempty(obj.histObj) == 0
1398                 obj.histObj.enableBufferedHistoryWrite();
1399             else
1400                 error('History is not being recorded so buffered writing cannot be enabled - run createHistory() first');
1401             end     
1402         end
1403         
1404         %% disableBufferedHistoryWrite() METHOD
1405         function disableBufferedHistoryWrite(obj)
1406             % disables writing of buffered history.
1407             
1408             if isempty(obj.histObj) == 0
1409                 obj.histObj.disableBufferedHistoryWrite();
1410             else
1411                 error('History is not being recorded so buffered writing cannot be disabled - run createHistory() first');
1412             end     
1413         end
1414         
1415         
1416         %% writeSamples()  METHOD
1417         function writeSamples(obj)
1418            % writes the data from the samples to text files specified in their 'outFile' property
1419            
1420            verbosePrint('Writing all samples to disk...',...
1421                             'sos_writeSamples_begin');
1422            
1423            for i=1:length(obj.samples)
1424                obj.samples(i).writeData();
1425            end
1426         end %writeSamples()
1427         
1428         %% writePopulations()  METHOD
1429         function writePopulations(obj)
1430            % writes the data from the populations to text files specified in their 'outFile' property
1431            
1432            verbosePrint('Writing all populations to disk...',...
1433                             'sos_writePopulations_begin');
1434            
1435            alreadyWritten = {};
1436            popWritten = 0;
1437            
1438            for i=1:length(obj.samples)
1439                if isempty(alreadyWritten)
1440                    popWritten = 0;
1441                else
1442                 for j = 1:length(alreadyWritten)
1443                     if (alreadyWritten{j} == obj.samples(i).population)
1444                         popWritten = 1;
1445                     end
1446                 end
1447                end
1448                
1449                if(popWritten == 0)
1450                 obj.samples(i).population.writeData();
1451                 alreadyWritten = [alreadyWritten {obj.samples(i).population}]; %#ok<AGROW>
1452                end
1453            end
1454         end %writeSamples()
1455             
1456         %% writeAll() METHOD
1457         function writeAll(obj)
1458             %write all populations and samples associated with SOS object to disk
1459             obj.writeSamples();
1460             obj.writePopulations();
1461             
1462         end % writeAll()
1463         
1464         %% addttest(sample1, sample2, s1ColName, s2ColName, type) METHOD
1465         function addttest(obj, varargin)
1466                 %adds a t-test to the list of analyses to run
1467                 newTest = genericStatTest.createStatTest('ttest','sosObj',obj,varargin{:});
1468                 obj.sosstattests = [obj.sosstattests {newTest}];   
1469                 
1470                 % if there is a history object at work, let it know that a
1471                 % new stat test has been added.
1472                 if isempty(obj.histObj) == 0
1473                     obj.histObj.addStatTestName(newTest.name);
1474                 end
1475                 
1476              verbosePrint('t-test added...',  ...
1477                     'sos_addttest_end');    
1478         end
1479         
1480         %% addtzest(sample1, sample2, s1ColName, s2ColName, type, ...) METHOD
1481         function addztest(obj, varargin)
1482                 %adds a z-test to the list of analyses to run
1483                 
1484                 % NOTE: If additional, more general types of z-tests are
1485                 % added, this will need to be conditionalized so that
1486                 % sosCorrelTest is not chosen by default...
1487                 newTest = genericStatTest.createStatTest('matchCorrel','sosObj',obj,varargin{:});
1488                 obj.sosstattests = [obj.sosstattests {newTest}];   
1489                 
1490                 % if there is a history object at work, let it know that a
1491                 % new stat test has been added.
1492                 if isempty(obj.histObj) == 0
1493                     obj.histObj.addStatTestName(newTest.name);
1494                 end
1495                 
1496              verbosePrint('z-test added...',  ...
1497                     'sos_addztest_end');    
1498         end        
1499 
1500         
1501         %% addtzest(sample1, sample2, s1ColName, s2ColName, type, ...) METHOD
1502         function addkstest(obj, varargin)
1503                 %adds a kstest to the list of analyses to run
1504                 
1505                 % NOTE: If additional, more general types of kstests are
1506                 % added, this will need to be conditionalized so that
1507                 % sosCorrelTest is not chosen by default...
1508                 newTest = genericStatTest.createStatTest('matchUniform','sosObj',obj,varargin{:});
1509                 obj.sosstattests = [obj.sosstattests {newTest}];   
1510                 
1511                 % if there is a history object at work, let it know that a
1512                 % new stat test has been added.
1513                 if isempty(obj.histObj) == 0
1514                     obj.histObj.addStatTestName(newTest.name);
1515                 end
1516                 
1517              verbosePrint('ks-test added...',  ...
1518                     'sos_addkstest_end');    
1519         end    
1520         
1521         
1522         
1523         function pass = doStatTests(obj,varargin)
1524             % runs the stat tests and indicates if all passed the user
1525             % hypotheses.
1526             %
1527             %PARAMETERS:
1528             %   'reportStyle'/'short'|'full' - style of report to be printed.
1529             %                               Either short or long
1530             %
1531             % RETURNS:
1532             %   NaN if user had no hypotheses about any of the tests that
1533             %   were run.  1 if there were hypotheses and all hypotheses
1534             %   passed their tests.  0  otherwise.
1535             
1536            p = inputParser;            
1537            p.addParamValue('reportStyle','short', ...
1538                     @(reportStyle)any([strcmp(reportStyle,'short') ...
1539                                     strcmp(reportStyle,'full') ...
1540                                     strcmp(reportStyle,'none')]));
1541            p.parse(varargin{:});
1542            
1543            if strcmp(p.Results.reportStyle,'none') == 0
1544                 verbosePrint([char(10),'Running all stat tests:'],...
1545                             'sos_doStatTests_begin');            
1546            end              
1547             pass = NaN;
1548             
1549             %
1550             
1551             
1552             for i=1:length(obj.sosstattests)
1553                 userHypothesis = obj.sosstattests{i}.runTest(varargin{:});
1554                 
1555                 if userHypothesis == 0
1556                     pass = 0;
1557                 % if a hypothesis passes (as opposed to being NaN), set the
1558                 % condition to pass;  If a hypothesis has already failed,
1559                 % it cannot pass
1560                 elseif userHypothesis == 1 && isnan(pass)
1561                     pass = 1;
1562                 end
1563             end
1564         end % doStatTests
1565         
1566         %% function deltaCostPercentiles() METHOD
1567         function deltaCostPercentiles(obj)
1568             % displays the breakdown of deltaCost values per decile and
1569             % other important percentiles
1570             
1571             %one over 100 so that we get the 100th decile
1572             deciles = 0:10:101;
1573             scores = prctile(obj.deltaCostLog, deciles);
1574             
1575             %let the user know how many observations entered into the delta
1576             %cost report
1577             verbosePrint(['Deciles for deltaCost in last block:',char(10), ...
1578                 '(',num2str(min(length(obj.deltaCostLog),obj.curIt)),...
1579                 ' iterations total; numIt < blockSize in first block)',char(10)], ...
1580                         'sos_deltaCostPercentiles_Header');
1581                     
1582                     
1583            values = '';
1584            for i=1:length(deciles)
1585                line=sprintf('\t%1.0f:  \t %+10.8f \n',deciles(i),scores(i)); 
1586                values = [values line]; %#ok<AGROW>
1587            end 
1588                
1589            verbosePrint(values, 'sos_deltaCostPercentiles_Scores');
1590              
1591            
1592            %print out information that could be helpful for configuring the
1593            %expAnneal function
1594            
1595            percentiles = [2.5 97.5];
1596            scores = prctile(obj.deltaCostLog, percentiles);
1597                 
1598            deltaCost95 = scores(2) - scores(1);
1599            
1600            msg = sprintf('\n\n 97.5th - 2.5th percentile deltaCost: %10.8f\n\n', ...
1601                         deltaCost95);
1602                     
1603            verbosePrint(msg, 'sos_deltaCostPercentiles_deltaCost95');
1604                 
1605              
1606         end
1607         
1608     end %end methods
1609     
1610     methods (Static)
1611         
1612         %% p = sosInputParser() STATIC METHOD
1613         function p = sosInputParser()
1614             % returns an input parser for the sos constructor args
1615             %
1616             %CALL:
1617             % p = sos.sosInputParser()
1618             %
1619             %SYNOPSIS:
1620             % returns an input parser for the sos constructor args
1621             %
1622             %PARAMETERS:
1623             %   Same as object constructor
1624             %
1625             %EXAMPLE:
1626             % p = sos.sosInputParser()
1627             %
1628             
1629              p = inputParser;
1630 
1631              %use NaN as null, since matlab doesn't support standard
1632              %NULL
1633              
1634              p.addParamValue('maxIt',  sos.maxIt_def, ...
1635                 @(maxIt)validateattributes(maxIt, {'numeric'}, ...
1636                 {'scalar', 'integer', 'positive', '>', 0}));
1637             p.addParamValue('pSwapFunction',sos.pSwapFunction_def, ...
1638                 @(pSwapFunction)ischar(pSwapFunction));
1639             p.addParamValue('targSampleCandSelectMethod',sos.targSampleCandSelectMethod_def, ...
1640                 @(targSampleCandSelectMethod)ischar(targSampleCandSelectMethod));
1641             p.addParamValue('feederdfCandSelectMethod',sos.feederdfCandSelectMethod_def, ...
1642                 @(feederdfCandSelectMethod)ischar(feederdfCandSelectMethod));
1643             p.addParamValue('reportInterval',sos.reportInterval_def, ...
1644                 @(reportInterval)validateattributes(reportInterval, {'numeric'}, ...
1645                 {'scalar', 'integer', 'positive', '>', 0}));
1646             p.addParamValue('stopFreezeIt',sos.stopFreezeIt_def, ...
1647                 @(stopFreezeIt)validateattributes(stopFreezeIt, {'numeric'}, ...
1648                 {'scalar', 'integer', 'positive', '>', 0}));
1649             p.addParamValue('statInterval',sos.statInterval_def, ...
1650                 @(statInterval)validateattributes(statInterval, {'numeric'}, ...
1651                 {'scalar', 'integer', 'positive', '>', 0})); 
1652             p.addParamValue('statTestReportStyle',sos.statTestReportStyle_def, ...
1653                 @(statTestReportStyle)any( ...
1654                     [strcmp(statTestReportStyle,'short'), ...
1655                      strcmp(statTestReportStyle,'full')]));
1656             p.addParamValue('blockSize',sos.blockSize_def, ...
1657                 @(statInterval)validateattributes(statInterval, {'numeric'}, ...
1658                 {'scalar', 'integer', 'positive', '>', 0})); 
1659             
1660       
1661         end
1662     end
1663     
1664     methods (Static, Access = private)
1665         %% parseConstructorArgs PRIVATE STATIC METHOD
1666         function p = parseConstructorArgs(varargin)
1667             %parses the constructor arguments
1668             %
1669             % CALL:
1670             % p = sos.parseConstructorArgs(varargin);
1671             %
1672             %parses the arguments from the  constructor.  Default
1673             %values are substituted where appropriate.  Returns a struct
1674             %with the parsed args
1675             %
1676             %PARAMETERS:
1677             % SAME as sos CONSTRUCTOR
1678             %
1679             %RETURNS:
1680             %    p - parsed constructor input
1681             
1682              varargin = varargin{1};
1683 
1684              p = sos.sosInputParser();
1685              p.parse(varargin{:});
1686              
1687              if any(strcmp(p.UsingDefaults,'maxIt'))
1688                  verbosePrint(['    Defaulting to ',  ...
1689                      num2str(p.Results.maxIt), ' maximum iterations'], ...
1690                     'sos_parseConstructor_defaultMaxIt');
1691              end
1692              
1693              if any(strcmp(p.UsingDefaults,'pSwapFunction'))
1694                  verbosePrint(['    Defaulting to ', ...
1695                      (p.Results.pSwapFunction), ' pSwapFunction'], ...
1696                     'sos_parseConstructor_pSwapFunction');                         
1697              end
1698 
1699              if any(strcmp(p.UsingDefaults,'targSampleCandSelectMethod'))
1700                  verbosePrint(['    Defaulting to ', ...
1701                      (p.Results.targSampleCandSelectMethod), ...
1702                      ' targSample Candidate Selection Method'], ...
1703                     'sos_parseConstructor_targSampleCandSelectMethod');      
1704              end
1705              
1706              if any(strcmp(p.UsingDefaults,'feederdfCandSelectMethod'))
1707                  verbosePrint(['    Defaulting to ', ...
1708                      (p.Results.feederdfCandSelectMethod), ...
1709                      ' feederdf Candidate Selection Method'], ...
1710                     'sos_parseConstructor_feederdfCandSelectMethod');      
1711              end
1712              
1713              if any(strcmp(p.UsingDefaults,'reportInterval'))
1714                  verbosePrint(['    Defaulting to ', ...
1715                      num2str(p.Results.reportInterval), ...
1716                      ' iterations between reports'], ...
1717                     'sos_parseConstructor_reportInterval');      
1718              end
1719              
1720              if any(strcmp(p.UsingDefaults,'stopFreezeIt'))
1721                  verbosePrint(['    Defaulting to ', ...
1722                      num2str(p.Results.stopFreezeIt), ...
1723                      ' of same frozen cost before stopping'], ...
1724                     'sos_parseConstructor_stopFreezeIt');      
1725              end
1726              
1727              
1728              if any(strcmp(p.UsingDefaults,'statInterval'))
1729                  verbosePrint(['    Defaulting to ', ...
1730                      num2str(p.Results.statInterval), ...
1731                      ' intervals between stat reports'], ...
1732                     'sos_parseConstructor_statReports');      
1733              end
1734              
1735              if any(strcmp(p.UsingDefaults,'statTestReportStyle'))
1736                  verbosePrint(['    Defaulting to ', ...
1737                      p.Results.statTestReportStyle, ...
1738                      ' style stat reports'], ...
1739                     'sos_parseConstructor_statReports');      
1740              end             
1741                          
1742         end % parseConstructorArgs
1743     end        
1744     
1745 end