softDistanceConstraint

 - parent class for soft distance constraints

 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 % - parent class for soft distance constraints
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 
0024 classdef softDistanceConstraint < softConstraint
0025     %% creates and supports soft distance constraints
0026     %
0027     % This class creates softDistanceConstraint objects that measure the
0028     % cost (in terms of minimizing or maximizing a distance of a particular stat)
0029     % associated with the current items in a set, and subsequent to
0030     % swapping a particular item with another items.
0031     %
0032     % Additional functionality and interface requirements are inherited
0033     % from softConstraint.
0034     %
0035     % Currently allows for the following:
0036     %
0037     %   min/max/OrderedMax of Group Means, Paired means, Group Stdev,
0038     %   Paired Stdev, Matching Mean to target, Matching stdev to target
0039     %
0040     %
0041     %PROPERTIES
0042     %     stat % the statistic to calculate
0043     %     paired - logical indicating if paired differences or group differences are bieng used
0044     %     comparison - handle to specific comparison calculator
0045     %     initStats % handle to global stat initialization method
0046     %     swStats % handle to local stat update method
0047     %     acceptsw % handle to local accept swap method
0048     %     rejectsw % handle to local reject swap method
0049     %     s1 - sample1
0050     %     s2 - sample2
0051     %     s1Col - s1 column index
0052     %     s2Col - s2 column index
0053     %     s1ColName - name of column in s1
0054     %     s2ColName - name of column in s2
0055     %     stat1 - current stat of sample 1
0056     %     stat2 - current stat of sample2
0057     %     swstat1 - stat of sample 1 if swap occurs
0058     %     swstat2 - stat of sample2 if swap occurs
0059     %     sumx1sq % sum of the squared values, x, in an array.  Used in stdev calculation
0060     %     meanx1 % square of the mean value of x in an array.  Used in stdev calculation
0061     %     swpsumx1sq % swap sum of the squared values, x, in an array.  Used in stdev calculation
0062     %     swpmeanx1 % swap square of the mean value of x in an array.  Used in stdev calculation
0063     %     sumx2sq % sum of the squared values, x, in an array.  Used in stdev calculation
0064     %     meanx2 % square of the mean value of x in an array.  Used in stdev calculation
0065     %     swpsumx2sq % swap sum of the squared values, x, in an array.  Used in stdev calculation
0066     %     swpmeanx2 % swap square of the mean value of x in an array.  Used in stdev calculation
0067     %     targVal % target value to match for 1 sample-to-value matching
0068     %     zTargVal % normalized target value
0069     %
0070     %METHODS
0071     %   obj =softDistanceConstraint(varargin) - constructor
0072     %   obj = construct2SampleSoftDistanceConstraint(varargin) % creates a constraint object for cases wehre 2 samples are implicated.
0073     %   obj = construct1SampleSoftDistanceConstraint(varargin)  % creates a constraint object for cases where 1 sample is implicated.
0074     %   cost = initCost() - Calculates, saves, and returns the cost value for the currentitems in the sample.
0075     %   swCost = swapCost(targSample,targSampleIndex, feederdf,feederdfIndex) - Calculates the new cost if items from targSample and feederdf were swapped.
0076     %   acceptSwap() - updates variables when the proposed swap is accepted
0077     %   initGroupMeans() - calculates the means of the samples
0078     %   initPairedMeans() % prepares data for paired mean calculation
0079     %   initGroupStdev() %calculates the stdevs of the samples
0080     %   initPairedStdev() % prepares data for paired stdev calculation
0081     %   initSingleMeanandTarg()  % calculates the mean of the sample1 and sets the target as the target stat
0082     %   obj = initSingleStdevandtarg(obj) %calculates the stdevs of the samples
0083     %   swGroupMeans(targSample,targSampleIndex, feederdf,feederdfIndex) - Calculates the new (swap) means if items from targSample and feederd were swapped.
0084     %   swPairedMeans(targSample,targSampleIndex, feederdf,feederdfIndex) % Calculates the new (swap) means if items from targSample and  feederd were swapped.
0085     %   swGroupStdev(targSample,targSampleIndex, feederdf,feederdfIndex) % Calculates the new (swap) stdevs if items from targSample and feederd were swapped.
0086     %   swPairedStdev(targSample,targSampleIndex, feederdf,feederdfIndex) % Calculates the new (swap) stdev if items from targSample and feederdf were swapped.
0087     %   swSingleMean(targSample,targSampleIndex, feederdf,feederdfIndex) % Calculates the new (swap) mean if items from targSample and  feederd were swapped.
0088     %   swSingleStdev(targSample,targSampleIndex, feederdf,feederdfIndex)  % Calculates the new (swap) stdevs if items from targSample and feederd were swapped.
0089     %   cost = minDiff(x1,x2) -  cost function minimized by reduced differences on the statistic
0090     %   cost = maxDiff(x1,x2) - cost function minimized by maximizing differences on the stat.  Many uses may want to use orderedMaxDiff, instead.
0091     %   cost = orderedMaxDiff(x1,x2) - cost function minimized by maximizing a difference between the two groups that adheres to s1.stat < s2.stat cost function minimized by reduced differences on the statistic for paired groups
0092     %   acceptSwapMeans() % updates variables when the proposed swap is accepted for  means
0093     %   acceptSwapStdev()   % updates variables when the proposed swap is accepted for standard deviations
0094     %   cost = acceptSwap() % generic reject function
0095     %   rejectSwapMeans() % updates variables when swap is rejected
0096     %   rejectSwapStdev()  % updates variables when the proposed swap is rejected
0097     %
0098     %METHODS (Static)
0099     %   p = softDistanceConstraintInputParser2Sample() generates an input parser with parameter / value pairs and validators for the constructor args
0100     %
0101     %METHODS (Static, Access = private)
0102     %   p = parse2SampleConstructorArgs(varargin) - parses the constructor args for cases invovling 2 samples
0103     %   p = parse1SampleConstructorArgs(varargin)  % parses the constructor args for cases invovling 1 sample
0104     
0105     %% PROPERTIES
0106     properties
0107         stat % the statistic to calculate
0108         paired % logical indicating if paired differences or group differences are bieng used
0109         comparison % handle to specific comparison calculator
0110         initStats % handle to global stat initialization method
0111         swStats % handle to local stat update method
0112         acceptsw % handle to local accept swap method
0113         rejectsw % handle to local reject swap method
0114         s1 % sample1
0115         s2 % sample2
0116         s1Col %s1 column index
0117         s2Col %s2 column index
0118         s1ColName % name of column in s1
0119         s2ColName % name of column in s2
0120         stat1 % current stat of sample 1
0121         stat2 % current stat of sample2
0122         swstat1 % stat of sample 1 if swap occurs
0123         swstat2 % stat of sample2 if swap occurs
0124         sumx1sq % sum of the squared values, x, in an array.  Used in stdev calculation
0125         meanx1 % square of the mean value of x in an array.  Used in stdev calculation
0126         swpsumx1sq % swap sum of the squared values, x, in an array.  Used in stdev calculation
0127         swpmeanx1 % swap square of the mean value of x in an array.  Used in stdev calculation
0128         sumx2sq % sum of the squared values, x, in an array.  Used in stdev calculation
0129         meanx2 % square of the mean value of x in an array.  Used in stdev calculation
0130         swpsumx2sq % swap sum of the squared values, x, in an array.  Used in stdev calculation
0131         swpmeanx2 % swap square of the mean value of x in an array.  Used in stdev calculation
0132         targVal % target value to match for 1 sample-to-value matching
0133         zTargVal % normalized target value
0134     end
0135     
0136     methods
0137         
0138         %% softDistanceConstraint CONSTRUCTOR
0139         function  obj = softDistanceConstraint(varargin)
0140             % Creates a softDistance Constraint object.  See the
0141             % constructors for the 2Sample and 1Sample case for additional
0142             % parameters
0143             
0144             
0145             %first, decide what type of object to create.  This will be
0146             %based on whether the desired function involves 2 means, or 1.
0147             p = inputParser;
0148             p.addParamValue('fnc','null', ...
0149                  @(fnc)any(strcmp({'min', 'max','orderedMax','match1SampleVal'},fnc)));
0150             p.KeepUnmatched = true;
0151             p.parse(varargin{:});
0152             
0153             %cases involving 2 samples
0154             if any([strcmp(p.Results.fnc,'min'), ...
0155                     strcmp(p.Results.fnc,'max'), ...
0156                     strcmp(p.Results.fnc,'orderedMax')])
0157                 obj = construct2SampleSoftDistanceConstraint(obj,varargin);
0158             % cases invovling 1 sample matched to a value
0159             elseif any(strcmp(p.Results.fnc,'match1SampleVal'))
0160                 obj = construct1SampleSoftDistanceConstraint(obj,varargin);                
0161             else 
0162                 error(['Could not create a soft constraint with <fnc>: ', ...
0163                         p.Results.fnc]);
0164             end
0165              
0166             obj.cost = NaN;
0167             obj.swCost = NaN;
0168                
0169             
0170             verbosePrint('Soft Distance Constraint has been created', ...
0171                     'softDistanceConstraint_Constructor_endObjCreation');
0172         end % constructor
0173         
0174         %% obj = construct2SampleSoftDistanceConstraint(varargin) METHOD
0175         function obj = construct2SampleSoftDistanceConstraint(obj,varargin)
0176             % creates a constraint object for cases where 2 samples are implicated.
0177             % CONSTRUCTOR - Creates a softDistanceConstraint object
0178             %
0179             % CALL:
0180             % softDistanceConstraint(varargin<defined below>)
0181             %
0182             % PARAMETERS:
0183             % REQUIRED:
0184             %   'sosObj'/sos object - the SOS object the constraint will be linked to, and which contains the samples the constraint operates on.
0185             %   'constraintType'/'soft' - the type of contraint - must be 'soft'
0186             %   'fnc'/'min'|'max'|'orderedMax' the distance function to create.  If orderedMax, s1 < s2 on the specified dimension
0187             %   'stat'/'mean'|'stdev' - the statistic to calculate the difference on
0188             %   'sample1'/sample - the first sample
0189             %   'sample2'/sample - the second sample
0190             %   's1ColName'/string - name of column in 1st sample
0191             %   's2ColName'/string - name of column in 2nd sample
0192             %   'paired'/logical - should pairwise or group level distances be calculated
0193             % OPTIONAL:
0194             %   'exponent'/numeric - defaults to 2 (quadratic difference)
0195             %   'weight'/numeric - defaults to 1 (equal weighting of all soft costs)
0196             %
0197             % EXAMPLE:
0198             % mySOS.addConstraint('sosObj',mySOS,'constraintType','soft', ...
0199             %   'fnc','min','stat','mean','sample1',s1,'sample2',s2,'s1ColName','Lg10WF', ...
0200             %   's2ColName','Lg10WF','exponent',2,'paired',false,'weight', 1);
0201             
0202                         
0203             p = softDistanceConstraint.parse2SampleConstructorArgs(varargin{:});
0204             
0205             if(p.Results.sosObj.containsSample(p.Results.sample1) == false)
0206                 error('Cannot create soft distance constraint: sos Object does not contain the sample1');
0207             end
0208             
0209             if(p.Results.sosObj.containsSample(p.Results.sample2) == false)
0210                 error('Cannot create soft distance constraint: sos Object does not contain sample2');
0211             end
0212             
0213             
0214             col1 = p.Results.sample1.colName2colNum(p.Results.s1ColName);           
0215             if(col1 == -1)
0216                 error('Specified column name not found in sample1');
0217             end
0218             
0219             if(strcmp(p.Results.sample1.format{col1},'%f') == 0)
0220                 error('Specified column is not of numeric (%f) format, so cannot use as hard bound');
0221             end           
0222                        
0223             col2 = p.Results.sample2.colName2colNum(p.Results.s2ColName);            
0224             if(col2 == -1)
0225                 error('Specified column name not found in sample');
0226             end
0227             
0228             if(strcmp(p.Results.sample2.format{col2},'%f') == 0)
0229                 error('Specified column is not of numeric (%f) format, so cannot use as hard bound; 2nd sample');
0230             end           
0231             
0232             %if the lengths are not going to be identical, do not run the
0233             %test if it is paired.  The init and swap methods will ensure
0234             %that if practically the data are NaN, the optimization will
0235             %stop.
0236             if(p.Results.paired == true)
0237                 if length(p.Results.sample1.n) ~= length(p.Results.sample2.n)
0238                     error('Sample sizes must be equal if using paired matching');
0239                 end
0240             end
0241             % Assign the comparison function appropriate for the constraint
0242             
0243             if(p.Results.paired ~= true & p.Results.paired ~= false)
0244                 error('<paired> parameter must be either true/false (1/0)');
0245             end
0246             
0247             
0248             %ASSIGN HANDLE TO STAT CALCULATION METHOD
0249             if(strcmp(p.Results.stat,'mean') && p.Results.paired == false)
0250                 obj.initStats = @obj.initGroupMeans;
0251                 obj.swStats = @obj.swGroupMeans;
0252                 obj.acceptsw = @obj.acceptSwapMeans;
0253                 obj.rejectsw = @obj.rejectSwapMeans;
0254             elseif(strcmp(p.Results.stat,'mean') && p.Results.paired == true)
0255                 obj.initStats = @obj.initPairedMeans;
0256                 obj.swStats = @obj.swPairedMeans;    
0257                 obj.acceptsw = @obj.acceptSwapMeans;
0258                 obj.rejectsw = @obj.rejectSwapMeans;
0259             elseif(strcmp(p.Results.stat,'stdev')  && p.Results.paired == false)
0260                 obj.initStats = @obj.initGroupStdev;
0261                 obj.swStats = @obj.swGroupStdev;
0262                 obj.acceptsw = @obj.acceptSwapStdev;  
0263                 obj.rejectsw = @obj.rejectSwapStdev;  
0264             elseif(strcmp(p.Results.stat,'stdev')  && p.Results.paired == true)
0265                 obj.initStats = @obj.initPairedStdev;
0266                 obj.swStats = @obj.swPairedStdev;
0267                 obj.acceptsw = @obj.acceptSwapStdev; 
0268                 obj.rejectsw = @obj.rejectSwapStdev;  
0269             else
0270                error('function not yet supported');
0271             end
0272             
0273             % ASSIGN HANDLE TO DIFFERENCE COMPARISON METHOD
0274             % Group differences:
0275             if(strcmp(p.Results.fnc,'min'))
0276                 obj.comparison = @obj.minDiff;
0277             elseif(strcmp(p.Results.fnc,'max'))
0278                 obj.comparison = @obj.maxDiff;                     
0279             elseif(strcmp(p.Results.fnc,'orderedMax'))
0280                 obj.comparison = @obj.orderedMaxDiff;   
0281             else
0282                 error('function not yet supported');
0283             end
0284             
0285             % parent properties
0286             obj.sosObj = p.Results.sosObj;
0287             obj.constraintType = p.Results.constraintType;
0288             obj.fnc = p.Results.fnc;
0289             
0290             
0291             obj.weight = p.Results.weight;
0292             obj.exp = p.Results.exponent;         
0293             
0294             obj.stat = p.Results.stat;
0295             obj.s1 = p.Results.sample1;
0296             obj.s2 = p.Results.sample2;
0297             obj.s1Col = col1;
0298             obj.s2Col = col2;           
0299             obj.s1ColName = p.Results.s1ColName;
0300             obj.s2ColName = p.Results.s2ColName;           
0301               
0302             obj.paired = p.Results.paired;
0303             
0304             
0305             % add the name and the label
0306             obj.label = [obj.constraintType,'_',obj.fnc,'_',...
0307                     obj.stat,'_',...
0308                     obj.s1.name,'_',obj.s1ColName,'_',...
0309                     obj.s2.name,'_',obj.s2ColName,'_',...
0310                     'p',num2str(obj.paired)','_w',...
0311                     num2str(obj.weight),'_e',num2str(obj.exp)];              
0312             if any(strcmp(p.UsingDefaults,'name'))                 
0313                 obj.name = obj.label;
0314             else
0315                  obj.name = p.Results.name;  
0316             end               
0317         end %construct2SampleSoftDistanceConstraint
0318         
0319         %% obj = construct1SampleSoftDistanceConstraint(varargin) METHOD
0320         function obj = construct1SampleSoftDistanceConstraint(obj,varargin)
0321             % creates a constraint object for cases where 1 sample is implicated.
0322             % CONSTRUCTOR - Creates a softDistanceConstraint object
0323             %
0324             % CALL:
0325             % softDistanceConstraint(varargin<defined below>)
0326             %
0327             % PARAMETERS:
0328             % REQUIRED:
0329             %   'sosObj'/sos object - the SOS object the constraint will be linked to, and which contains the samples the constraint operates on.
0330             %   'constraintType'/'soft' - the type of contraint - must be 'soft'
0331             %   'fnc'/'match1SampleVal' the distance function to create.  If orderedMax, s1 < s2 on the specified dimension
0332             %   'stat'/'mean'|'stdev' - the statistic to calculate the difference on
0333             %   'sample1'/sample - the first sample
0334             %   's1ColName'/string - name of column in 1st sample
0335             % OPTIONAL:
0336             %   'exponent'/numeric - defaults to 2 (quadratic difference)
0337             %   'weight'/numeric - defaults to 1 (equal weighting of all soft costs)
0338             %   'targVal'/numeric - value to match the stat to.
0339             %
0340     
0341             
0342             p = softDistanceConstraint.parse1SampleConstructorArgs(varargin{:});
0343             
0344             if(p.Results.sosObj.containsSample(p.Results.sample1) == false)
0345                 error('Cannot create soft distance constraint: sos Object does not contain the sample1');
0346             end
0347             
0348             
0349             col1 = p.Results.sample1.colName2colNum(p.Results.s1ColName);           
0350             if(col1 == -1)
0351                 error('Specified column name not found in sample1');
0352             end
0353             
0354             if(strcmp(p.Results.sample1.format{col1},'%f') == 0)
0355                 error('Specified column is not of numeric (%f) format, so cannot use as hard bound');
0356             end           
0357                   
0358             
0359             % Assign the comparison function appropriate for the constraint
0360             
0361             %ASSIGN HANDLE TO STAT CALCULATION METHOD
0362             if(strcmp(p.Results.stat,'mean'))
0363                 obj.initStats = @obj.initSingleMeanandTarg;
0364                 obj.swStats = @obj.swSingleMean;
0365                 obj.acceptsw = @obj.acceptSwapMeans;
0366                 obj.rejectsw = @obj.rejectSwapMeans;
0367                 
0368                 obj.comparison = @obj.minDiff;
0369                 
0370             elseif(strcmp(p.Results.stat,'stdev'))
0371                 obj.initStats = @obj.initSingleStdevandTarg;
0372                 obj.swStats = @obj.swSingleStdev;
0373                 obj.acceptsw = @obj.acceptSwapStdev;  
0374                 obj.rejectsw = @obj.rejectSwapStdev;
0375                 
0376                 obj.comparison = @obj.minDiff;
0377             else
0378                error('function not yet supported');
0379             end
0380             
0381             
0382             %Finally, need to convert the raw value supplied by the user
0383             %into the standardized value used for comparison purposes.  We
0384             %do this by conslulting the SOS object's normalization data.
0385             
0386             %find the normalization params for this data column:
0387             
0388             obj.targVal = p.Results.targVal;
0389             obj.zTargVal = NaN;
0390             
0391             % parent properties
0392             obj.sosObj = p.Results.sosObj;
0393             obj.constraintType = p.Results.constraintType;
0394             obj.fnc = p.Results.fnc;
0395             
0396             
0397             obj.weight = p.Results.weight;
0398             
0399             
0400             obj.stat = p.Results.stat;
0401             obj.s1 = p.Results.sample1;
0402             obj.s1Col = col1;   
0403             obj.s2 = NaN;
0404             obj.s1ColName = p.Results.s1ColName;         
0405             obj.exp = p.Results.exponent;           
0406             
0407             % add the name and the label
0408             obj.label = [obj.constraintType,'_',obj.fnc,'_',...
0409                     obj.stat,'_',num2str(obj.targVal),'_',obj.s1.name,'_',...
0410                     obj.s1ColName,'_w',...
0411                     num2str(obj.weight),'_e',num2str(obj.exp)];              
0412             if any(strcmp(p.UsingDefaults,'name'))                 
0413                 obj.name = obj.label;
0414             else
0415                  obj.name = p.Results.name;  
0416             end            
0417             
0418         end %construct1SampleSoftDistanceConstraint
0419         
0420         
0421         
0422         %% cost = initCost() METHOD
0423         function cost = initCost(obj)
0424             % Calculates, saves, and returns the cost value for the current items in the sample.
0425             %
0426             %CALL:
0427             %   <softDistanceConstraint>.initCost();
0428             %
0429             %SYNOPSIS:
0430             % Calculates the cost for the current items in the sample based
0431             % on the specified constraint.  This calculation is done on the
0432             % global dataset, rather than using the differential local
0433             % computation used when calculating swap costs.
0434             %
0435           
0436             %init the stats, then compare them.  Return the calculated cost
0437             obj = obj.initStats(); 
0438             cost = obj.comparison(obj.stat1,obj.stat2);
0439             
0440             obj.swstat1 = NaN;
0441             obj.swstat2 = NaN;
0442             
0443             obj.cost = cost;
0444             obj.swCost = NaN;
0445         end
0446 
0447         %%  swCost(targSample,targSampleIndex, feederdf,feederdfIndex) FUNCTION
0448         function swCost = swapCost(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
0449             % Calculates the new cost if items from targSample and feederdf were swapped.
0450             %
0451             %By definition, if this method is called it means that at least
0452             % one of the two swap objects is implicated in this function
0453             %
0454             %PARAMETERS:
0455             %   targSample - the target sample (i.e., the object that will call it's swapSample() method if a swap later occurs).
0456             %   targSampleIndex - row index of item to swap
0457             %   feederdf - dataframe (sample/pop) containin the other item to swap
0458             %   feederdfIndex - row index of item to swap.
0459                        
0460            if (obj.s1 ~= targSample && obj.s2 ~= targSample ... 
0461                    && obj.s1 ~= feederdf && obj.s2  ~= feederdf)
0462                error('swCost called, but no sample part of this cost function');
0463            end
0464 
0465            % update the stats, then calculate cost of new stats
0466            obj.swStats(targSample,targSampleIndex, feederdf,feederdfIndex);
0467  
0468             swCost = obj.comparison(obj.swstat1,obj.swstat2);
0469             obj.swCost = swCost;
0470               
0471         end
0472         
0473         
0474         %% initGroupMeans() METHOD
0475         function obj = initGroupMeans(obj)
0476             % calculates the means of the samples
0477             obj.stat1= mean(obj.s1.zdata{obj.s1Col});
0478             obj.stat2= mean(obj.s2.zdata{obj.s2Col});
0479             
0480             
0481             if isnan(obj.stat1) || isnan(obj.stat2)
0482                 error('NaN obtained during initstat computation (perhaps there is missing data for an item?)');
0483             end
0484         end
0485 
0486         %% initPairedMeans() METHOD
0487         function obj = initPairedMeans(obj)
0488             % prepares data for paired mean calculation
0489                       
0490             % make sure there are no nan's.
0491             if (any(isnan(obj.s1.zdata{obj.s1Col})))
0492                 error('NaN obtained during initstat computation of sample 1 (perhaps there is missing data for an item?)');
0493             elseif (any(isnan(obj.s2.zdata{obj.s2Col})))
0494                 error('NaN obtained during initstat computation of sample 2 (perhaps there is missing data for an item?)');
0495             end
0496             
0497             % calculate the statistic.  In contrast to the group mean
0498             % statistic, here we'll calculate the difference in each pair
0499             % of observations and sum that in stat2, and leave stat1  ==0
0500             
0501             tmpstat2 = 0;
0502             for i=1:obj.s1.n
0503                 x = obj.s2.zdata{obj.s2Col}(i) - obj.s1.zdata{obj.s1Col}(i);
0504                 
0505                 % when trying to minimize, any difference, in either
0506                 % direction, is bad.  When trying to maximize, then let the
0507                 % sign enter into the equation.
0508                 if strcmp(obj.fnc,'min')
0509                     x = abs(x);
0510                 end
0511 
0512                  tmpstat2 = tmpstat2 + x;               
0513             end
0514             
0515             % calculate the mean of x
0516             obj.stat2 = tmpstat2/obj.s1.n;
0517             obj.stat1 = 0; % always relative to no difference
0518             
0519             
0520         end        
0521         
0522         %% initStdev(obj) METHOD
0523         function obj = initGroupStdev(obj)
0524             %calculates the stdevs of the samples
0525             
0526             obj.sumx1sq = sum((obj.s1.zdata{obj.s1Col}).^2);
0527             obj.meanx1 = mean(obj.s1.zdata{obj.s1Col});
0528             n1 = length(obj.s1.zdata{obj.s1Col});
0529             
0530             %use the computational formula to efficiently calculate the sum
0531             stdevx1 = (1/(n1-1) * obj.sumx1sq ...
0532                         - n1/(n1-1) * obj.meanx1^2)^0.5;
0533                     
0534             obj.sumx2sq = sum((obj.s2.zdata{obj.s2Col}).^2);
0535             obj.meanx2 = mean(obj.s2.zdata{obj.s2Col});
0536             n2 = length(obj.s2.zdata{obj.s2Col});
0537             
0538             %use the computational formula to efficiently calculate the sum
0539             stdevx2 = (1/(n2-1) * obj.sumx2sq ...
0540                         - n2/(n2-1) * obj.meanx2^2)^0.5;
0541                     
0542             obj.stat1 = stdevx1;
0543             obj.stat2 = stdevx2;
0544             
0545             if isnan(obj.stat1) || isnan(obj.stat2)
0546                 error('NaN obtained during initstat computation (perhaps there is missing data for an item?)');
0547             end     
0548             
0549         end
0550             
0551         %% initPairedMeans() METHOD
0552         function obj = initPairedStdev(obj)
0553             % prepares data for paired mean calculation
0554                       
0555             % make sure there are no nan's.
0556             if (any(isnan(obj.s1.zdata{obj.s1Col})))
0557                 error('NaN obtained during initstat computation of sample 1 (perhaps there is missing data for an item?)');
0558             elseif (any(isnan(obj.s2.zdata{obj.s2Col})))
0559                 error('NaN obtained during initstat computation of sample 2 (perhaps there is missing data for an item?)');
0560             end
0561             
0562             % calculate the statistic.  In contrast to the group mean
0563             % statistic, here we'll calculate the difference in each pair
0564             % of observations and sum that in stat2, and leave stat1  ==0
0565             
0566             obj.sumx2sq = 0;
0567             obj.meanx2 = 0;
0568             n = length(obj.s2.zdata{obj.s2Col});
0569             
0570             for i=1:n
0571                 x = obj.s2.zdata{obj.s2Col}(i) - obj.s1.zdata{obj.s1Col}(i);
0572                 obj.sumx2sq = obj.sumx2sq + x^2;
0573                 obj.meanx2 = obj.meanx2 + x/n;
0574                             
0575             end
0576             
0577             %use the computational formula to efficiently calculate the sum
0578             stdevx2 = (1/(n-1) * obj.sumx2sq ...
0579                         - n/(n-1) * obj.meanx2^2)^0.5;            
0580             
0581             
0582             % calculate the mean of x
0583             obj.stat2 = stdevx2;
0584             obj.stat1 = 0; % always relative to no difference
0585             
0586             
0587         end    
0588 
0589         %% initSingleMeanandTarg() METHOD
0590         function obj = initSingleMeanandTarg(obj)
0591             % calculates the mean of the sample1 and sets the target as the target stat
0592             
0593             obj.stat2= mean(obj.s1.zdata{obj.s1Col});
0594             
0595              %normalize the target value to match against
0596             targCol = obj.sosObj.allData.colName2colNum(obj.s1ColName);           
0597             if(targCol == -1)
0598                 error('Normalized data for this value may not yet exist (try normalizing first)');
0599             end    
0600             
0601             obj.zTargVal = (obj.targVal - obj.sosObj.allDataColMean(targCol)) ...
0602                             / obj.sosObj.allDataColStd(targCol);
0603                         
0604            obj.stat1 = obj.zTargVal;
0605                         
0606             if isnan(obj.stat1) || isnan(obj.stat2)
0607                 error('NaN obtained during initstat computation (perhaps there is missing data for an item?)');
0608             end
0609         end
0610 
0611         %% initSingleStdevandtarg(obj) METHOD
0612         function obj = initSingleStdevandTarg(obj)
0613             %calculates the stdevs of the samples
0614                     
0615             obj.sumx2sq = sum((obj.s1.zdata{obj.s1Col}).^2);
0616             obj.meanx2 = mean(obj.s1.zdata{obj.s1Col});
0617             n2 = length(obj.s1.zdata{obj.s1Col});
0618             
0619             %use the computational formula to efficiently calculate the sum
0620             stdevx2 = (1/(n2-1) * obj.sumx2sq ...
0621                         - n2/(n2-1) * obj.meanx2^2)^0.5;
0622             
0623             obj.stat2 = stdevx2;
0624 
0625              %normalize the target value to match against
0626             targCol = obj.sosObj.allData.colName2colNum(obj.s1ColName);           
0627             if(targCol == -1)
0628                 error('Normalized data for this value may not yet exist (try normalizing first)');
0629             end    
0630             
0631             obj.zTargVal = (obj.targVal) ...
0632                             / obj.sosObj.allDataColStd(targCol);
0633                         
0634             obj.stat1 = obj.zTargVal;            
0635             
0636             
0637             if isnan(obj.stat1) || isnan(obj.stat2)
0638                 error('NaN obtained during initstat computation (perhaps there is missing data for an item?)');
0639             end     
0640             
0641         end
0642         
0643         %% swGroupMeans(targSample,targSampleIndex,feederdf,feederdfIndex) FUNCTION
0644         function swGroupMeans(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
0645            % Calculates the new (swap) means if items from targSample and feederd were swapped.
0646            %
0647            % Inputs are the same as for swapCost()
0648          
0649            tempswm1 = obj.stat1;
0650            tempswm2 = obj.stat2;
0651        
0652             %do the adjustments of the means
0653             %basic idea is to take out the old value from that mean and
0654             %then put in the new one.  Do that seperately for each
0655             %dataframe to cover all eventualities (i.e., if either
0656             %targSample or feederdf, or both, are part of the current cost
0657             %function.
0658             
0659             if targSample == obj.s1
0660                 tempswm1 = tempswm1 - (targSample.zdata{obj.s1Col}(targSampleIndex)/length(targSample.zdata{obj.s1Col}));
0661                 tempswm1 = tempswm1 + (feederdf.zdata{obj.s1Col}(feederdfIndex)/length(targSample.zdata{obj.s1Col}));
0662             elseif targSample == obj.s2
0663                 tempswm2 = tempswm2 - (targSample.zdata{obj.s2Col}(targSampleIndex)/length(targSample.zdata{obj.s2Col}));
0664                 tempswm2 = tempswm2 + (feederdf.zdata{obj.s2Col}(feederdfIndex)/length(targSample.zdata{obj.s2Col}));
0665             end
0666                                
0667             if feederdf == obj.s1
0668                 tempswm1 = tempswm1 - (feederdf.zdata{obj.s1Col}(feederdfIndex)/length(feederdf.zdata{obj.s1Col}));
0669                 tempswm1 = tempswm1 + (targSample.zdata{obj.s1Col}(targSampleIndex)/length(feederdf.zdata{obj.s1Col}));
0670             elseif feederdf == obj.s2 
0671                 tempswm2 = tempswm2 - (feederdf.zdata{obj.s2Col}(feederdfIndex)/length(feederdf.zdata{obj.s2Col}));
0672                 tempswm2 = tempswm2 + (targSample.zdata{obj.s2Col}(targSampleIndex)/length(feederdf.zdata{obj.s2Col}));
0673             end           
0674             
0675             obj.swstat1 = tempswm1;
0676             obj.swstat2 = tempswm2;
0677                
0678             if isnan(obj.swstat1) || isnan(obj.swstat2)
0679                 error('NaN obtained during swstat computation (perhaps there is missing data for an item?)');
0680             end
0681         end
0682         
0683 
0684         
0685         
0686         %% swPairedMeans(targSample,targSampleIndex,feederdf,feederdfIndex) FUNCTION
0687         function swPairedMeans(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
0688            % Calculates the new (swap) means if items from targSample and feederd were swapped.
0689            %
0690            % Inputs are the same as for swapCost()
0691          
0692            tempswm2 = obj.stat2;
0693            %tempswm2 should remain 0 at all times
0694        
0695             %do the adjustments of the means
0696             %basic idea is to take out the old value from that mean and
0697             %then put in the new one.  Do that seperately for each
0698             %dataframe to cover all eventualities (i.e., if either
0699             %targSample or feederdf, or both, are part of the current cost
0700             %function.
0701 
0702                        
0703             nobs = length(obj.s1.zdata{obj.s1Col});
0704 
0705             if targSample == obj.s1
0706                 %take the existing contribution to the stat due to the data
0707                 %out.
0708                 curx = (obj.s2.zdata{obj.s2Col}(targSampleIndex) - ...
0709                         targSample.zdata{obj.s1Col}(targSampleIndex))/nobs;
0710                 
0711                 if strcmp(obj.fnc,'min')
0712                     curx = abs(curx);
0713                 end
0714                 
0715                 tempswm2 = tempswm2 - curx;
0716                 
0717                 %add in the new value
0718                 newx = (obj.s2.zdata{obj.s2Col}(targSampleIndex) - ...
0719                           feederdf.zdata{obj.s1Col}(feederdfIndex)) / nobs;
0720                     
0721                 if strcmp(obj.fnc,'min')
0722                     newx = abs(newx);
0723                 end     
0724                 
0725                 tempswm2 = tempswm2 + newx;     
0726                 
0727                 % update is complete
0728             elseif targSample == obj.s2
0729                 
0730                 curx  = (targSample.zdata{obj.s2Col}(targSampleIndex) - ...
0731                             obj.s1.zdata{obj.s1Col}(targSampleIndex))/nobs;
0732                         
0733                 if strcmp(obj.fnc,'min')
0734                     curx = abs(curx);
0735                 end
0736 
0737                 tempswm2 = tempswm2 - curx;         
0738                 
0739                 newx = (feederdf.zdata{obj.s2Col}(feederdfIndex) - ...
0740                         obj.s1.zdata{obj.s1Col}(targSampleIndex))/nobs;
0741 
0742                 if strcmp(obj.fnc,'min')
0743                     newx = abs(newx);
0744                 end     
0745                 
0746                 tempswm2 = tempswm2 + newx;    
0747                 
0748             end
0749             
0750             
0751             if feederdf == obj.s1
0752                 
0753                 curx = (obj.s2.zdata{obj.s2Col}(feederdfIndex) - ...
0754                           feederdf.zdata{obj.s1Col}(feederdfIndex)) / nobs;
0755                 
0756                 if strcmp(obj.fnc,'min')
0757                     curx = abs(curx);
0758                 end                
0759         
0760                 tempswm2 = tempswm2 - curx;       
0761                 
0762                 newx = (obj.s2.zdata{obj.s2Col}(feederdfIndex) - ...
0763                         targSample.zdata{obj.s1Col}(targSampleIndex))/nobs;
0764                 
0765                  if strcmp(obj.fnc,'min')
0766                     newx = abs(newx);
0767                 end     
0768                 
0769                 tempswm2 = tempswm2 + newx;   
0770                 
0771                 
0772             elseif feederdf == obj.s2 
0773                 curx = (feederdf.zdata{obj.s2Col}(feederdfIndex) - ...
0774                         obj.s1.zdata{obj.s1Col}(feederdfIndex))/nobs;
0775                     
0776                  
0777                 if strcmp(obj.fnc,'min')
0778                     curx = abs(curx);
0779                 end                
0780         
0781                 tempswm2 = tempswm2 - curx;       
0782                 
0783                 newx = (targSample.zdata{obj.s2Col}(targSampleIndex) - ...
0784                             obj.s1.zdata{obj.s1Col}(feederdfIndex))/nobs;
0785                 
0786                 if strcmp(obj.fnc,'min')
0787                     newx = abs(newx);
0788                 end     
0789                 
0790                 tempswm2 = tempswm2 + newx;      
0791                 
0792             end
0793             
0794             obj.swstat2 = tempswm2;
0795             obj.swstat1 = obj.stat1; % again, just copy the zero over
0796                
0797             if isnan(obj.swstat1)
0798                 error('NaN obtained during swstat computation (perhaps there is missing data for an item?)');
0799             end
0800         end %swPairedMeans
0801 
0802         %% swGroupStdev(targSample,targSampleIndex, feederdf,feederdfIndex)
0803         function swGroupStdev(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
0804            % Calculates the new (swap) stdevs if items from targSample and feederd were swapped.
0805            %
0806            % Inputs are the same as for swapCost()
0807             
0808             n1 = length(obj.s1.zdata{obj.s1Col});
0809             n2 = length(obj.s2.zdata{obj.s2Col});           
0810                        
0811             obj.swpsumx1sq = obj.sumx1sq;
0812             obj.swpmeanx1 = obj.meanx1;
0813             
0814             obj.swpsumx2sq = obj.sumx2sq;
0815             obj.swpmeanx2 = obj.meanx2;
0816             
0817             if targSample == obj.s1
0818                 %must update 2 params %first, take out the old observation, then add in thew new
0819                 obj.swpsumx1sq = obj.swpsumx1sq - (targSample.zdata{obj.s1Col}(targSampleIndex))^2;
0820                 obj.swpmeanx1 = obj.swpmeanx1 - (targSample.zdata{obj.s1Col}(targSampleIndex))/n1;                
0821                 obj.swpsumx1sq = obj.swpsumx1sq + (feederdf.zdata{obj.s1Col}(feederdfIndex))^2;
0822                 obj.swpmeanx1 = obj.swpmeanx1 + (feederdf.zdata{obj.s1Col}(feederdfIndex))/n1;
0823             elseif targSample == obj.s2
0824                 obj.swpsumx2sq = obj.swpsumx2sq - (targSample.zdata{obj.s2Col}(targSampleIndex))^2;
0825                 obj.swpmeanx2 = obj.swpmeanx2 - (targSample.zdata{obj.s2Col}(targSampleIndex))/n2;                
0826                 obj.swpsumx2sq = obj.swpsumx2sq + (feederdf.zdata{obj.s2Col}(feederdfIndex))^2;
0827                 obj.swpmeanx2 = obj.swpmeanx2 + (feederdf.zdata{obj.s2Col}(feederdfIndex))/n2;     
0828             end
0829             
0830             if feederdf == obj.s1
0831                 obj.swpsumx1sq = obj.swpsumx1sq -  (feederdf.zdata{obj.s1Col}(feederdfIndex))^2;
0832                 obj.swpmeanx1 = obj.swpmeanx1 - (feederdf.zdata{obj.s1Col}(feederdfIndex))/n1;                
0833                 obj.swpsumx1sq = obj.swpsumx1sq + (targSample.zdata{obj.s1Col}(targSampleIndex))^2;
0834                 obj.swpmeanx1 = obj.swpmeanx1 + (targSample.zdata{obj.s1Col}(targSampleIndex))/n1;   
0835             elseif feederdf == obj.s2
0836                 obj.swpsumx2sq = obj.swpsumx2sq -  (feederdf.zdata{obj.s2Col}(feederdfIndex))^2;
0837                 obj.swpmeanx2 = obj.swpmeanx2 - (feederdf.zdata{obj.s2Col}(feederdfIndex))/n2;                
0838                 obj.swpsumx2sq = obj.swpsumx2sq + (targSample.zdata{obj.s2Col}(targSampleIndex))^2;
0839                 obj.swpmeanx2 = obj.swpmeanx2 + (targSample.zdata{obj.s2Col}(targSampleIndex))/n2;  
0840             end
0841      
0842             stdevx1 = (1/(n1-1) * obj.swpsumx1sq ...
0843                         - n1/(n1-1) * obj.swpmeanx1^2)^0.5;
0844                     
0845             stdevx2 = (1/(n2-1) * obj.swpsumx2sq ...
0846                         - n2/(n2-1) * obj.swpmeanx2^2)^0.5;
0847                     
0848             obj.swstat1 = stdevx1;
0849             obj.swstat2 = stdevx2;
0850             
0851             if isnan(obj.swstat1) || isnan(obj.swstat2)
0852                 error('NaN obtained during swstat computation (perhaps there is missing data for an item?)');
0853             end            
0854           
0855         end
0856         
0857         %% swPairedStdev(targSample,targSampleIndex,feederdf,feederdfIndex) FUNCTION
0858         function swPairedStdev(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
0859            % Calculates the new (swap) stdev if items from targSample and feederd were swapped.
0860            %
0861            %NOTE: When repeated tens of thousands of times, this manner of
0862            %updated std via local difference calculations tends to lead to
0863            %some very small imprecisions relative to re-calculating the
0864            %actual values from scratch.  These deviations appear to be
0865            %trivially small in all cases examined to date though.  This
0866            %fact is mentioned here presently only because the
0867            %detail-oriented individual may find values like  0.00000001
0868            %instead of 0, which nevertheless
0869            %likely will have no practical consequence
0870            %for the matching algorithm.
0871            %
0872            % Inputs are the same as for swapCost()
0873 
0874            
0875            
0876             obj.swpsumx2sq = obj.sumx2sq;
0877             obj.swpmeanx2 = obj.meanx2;     
0878             n = length(obj.s1.zdata{obj.s1Col});
0879                        
0880             if targSample == obj.s1
0881                 %take the existing contribution to the stat due to the data
0882                 %out.
0883                 curx = (obj.s2.zdata{obj.s2Col}(targSampleIndex) - ...
0884                         targSample.zdata{obj.s1Col}(targSampleIndex));
0885                                
0886                 obj.swpsumx2sq = obj.swpsumx2sq - curx^2;
0887                 obj.swpmeanx2 = obj.swpmeanx2 - curx/n;
0888                 
0889                 %add in the new value
0890                 newx = (obj.s2.zdata{obj.s2Col}(targSampleIndex) - ...
0891                           feederdf.zdata{obj.s1Col}(feederdfIndex));
0892                     
0893                 obj.swpsumx2sq = obj.swpsumx2sq + newx^2;
0894                 obj.swpmeanx2 = obj.swpmeanx2 + newx/n;                
0895                                 
0896                 % update is complete
0897             elseif targSample == obj.s2
0898                 
0899                 curx  = (targSample.zdata{obj.s2Col}(targSampleIndex) - ...
0900                             obj.s1.zdata{obj.s1Col}(targSampleIndex));
0901                                         
0902                 obj.swpsumx2sq = obj.swpsumx2sq - curx^2;
0903                 obj.swpmeanx2 = obj.swpmeanx2 - curx/n;
0904                 
0905                 newx = (feederdf.zdata{obj.s2Col}(feederdfIndex) - ...
0906                         obj.s1.zdata{obj.s1Col}(targSampleIndex));
0907                     
0908                 obj.swpsumx2sq = obj.swpsumx2sq + newx^2;
0909                 obj.swpmeanx2 = obj.swpmeanx2 + newx/n;       
0910                 
0911             end
0912             
0913             
0914             if feederdf == obj.s1
0915                 
0916                 curx = (obj.s2.zdata{obj.s2Col}(feederdfIndex) - ...
0917                           feederdf.zdata{obj.s1Col}(feederdfIndex));
0918                       
0919                 obj.swpsumx2sq = obj.swpsumx2sq - curx^2;
0920                 obj.swpmeanx2 = obj.swpmeanx2 - curx/n;
0921                 
0922                 newx = (obj.s2.zdata{obj.s2Col}(feederdfIndex) - ...
0923                         targSample.zdata{obj.s1Col}(targSampleIndex));
0924                 
0925                 obj.swpsumx2sq = obj.swpsumx2sq + newx^2;
0926                 obj.swpmeanx2 = obj.swpmeanx2 + newx/n;       
0927                 
0928                 
0929             elseif feederdf == obj.s2 
0930                 curx = (feederdf.zdata{obj.s2Col}(feederdfIndex) - ...
0931                         obj.s1.zdata{obj.s1Col}(feederdfIndex));
0932                     
0933                 obj.swpsumx2sq = obj.swpsumx2sq - curx^2;
0934                 obj.swpmeanx2 = obj.swpmeanx2 - curx/n;
0935                 
0936                 newx = (targSample.zdata{obj.s2Col}(targSampleIndex) - ...
0937                             obj.s1.zdata{obj.s1Col}(feederdfIndex));
0938                         
0939                 obj.swpsumx2sq = obj.swpsumx2sq + newx^2;
0940                 obj.swpmeanx2 = obj.swpmeanx2 + newx/n;       
0941                 
0942             end
0943             
0944             % Special override:
0945             % if it happens that targSample and feederdf both correspond to
0946             % the samples for which stdev is being calculated, and that
0947             % the exact same item is being swapped out in both cases, the
0948             % swpsumx2sq should not change because the two swaps should
0949             % cancel each other out.  However, in the current
0950             % implementation this doesn't happen and instead we get 2x the
0951             % negative subtraction because the second calculation is not
0952             % aware of the shifts completed in the first.  This can,
0953             % however, be fixed in the present case quite simply by not
0954             % changing swpsumx2sq in this specific instance.  This only
0955             % happens in this case because the squaring of the terms leads
0956             % to the sign of the calculation being ignored, and as such the
0957             % calculation of swap means should not suffer in the same
0958             % fashion.
0959             if((targSample == obj.s1 || targSample == obj.s2) && ...
0960                     (feederdf == obj.s1 || feederdf == obj.s2) && ...
0961                 targSampleIndex == feederdfIndex)
0962                 obj.swpsumx2sq = obj.sumx2sq;
0963             end
0964             
0965             % calculate the updated stdev:
0966            %use the computational formula to efficiently calculate the sum
0967             swpstdevx2 = (1/(n-1) * obj.swpsumx2sq ...
0968                         - n/(n-1) * obj.swpmeanx2^2)^0.5;  
0969                     
0970             obj.swstat2 = swpstdevx2;
0971             obj.swstat1 = obj.stat1; % again, just copy the zero over
0972                
0973             if isnan(obj.swstat1)
0974                 error('NaN obtained during swstat computation (perhaps there is missing data for an item?)');
0975             end
0976         end %swPairedStdev
0977        
0978         
0979         %% swSingleMean(targSample,targSampleIndex,feederdf,feederdfIndex) METHOD
0980         function swSingleMean(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
0981            % Calculates the new (swap) mean if items from targSample and feederd were swapped.
0982            %
0983            % Inputs are the same as for swapCost()
0984          
0985            tempswm2 = obj.stat2;
0986        
0987             %do the adjustments of the means
0988             %basic idea is to take out the old value from that mean and
0989             %then put in the new one.  Do that seperately for each
0990             %dataframe to cover all eventualities (i.e., if either
0991             %targSample or feederdf, or both, are part of the current cost
0992             %function.
0993             
0994             if targSample == obj.s1
0995                 tempswm2 = tempswm2 - (targSample.zdata{obj.s1Col}(targSampleIndex)/length(targSample.zdata{obj.s1Col}));
0996                 tempswm2 = tempswm2 + (feederdf.zdata{obj.s1Col}(feederdfIndex)/length(targSample.zdata{obj.s1Col}));
0997             end
0998             
0999             
1000             if feederdf == obj.s1
1001                 tempswm2 = tempswm2 - (feederdf.zdata{obj.s1Col}(feederdfIndex)/length(feederdf.zdata{obj.s1Col}));
1002                 tempswm2 = tempswm2 + (targSample.zdata{obj.s1Col}(targSampleIndex)/length(feederdf.zdata{obj.s1Col}));
1003             end
1004             
1005             obj.swstat2 = tempswm2;
1006             obj.swstat1 = obj.stat1;
1007                
1008             if isnan(obj.swstat1) || isnan(obj.swstat2)
1009                 error('NaN obtained during swstat computation (perhaps there is missing data for an item?)');
1010             end
1011         end
1012         
1013         %% swSingleStdev(targSample,targSampleIndex, feederdf,feederdfIndex)
1014         function swSingleStdev(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
1015            % Calculates the new (swap) stdevs if items from targSample and feederd were swapped.
1016            %
1017            % Inputs are the same as for swapCost()
1018             
1019             n2 = length(obj.s1.zdata{obj.s1Col});           
1020                        
1021             obj.swpsumx2sq = obj.sumx2sq;
1022             obj.swpmeanx2 = obj.meanx2;
1023             
1024             if targSample == obj.s1
1025                 %must update 2 params %first, take out the old observation, then add in thew new
1026                 obj.swpsumx2sq = obj.swpsumx2sq - (targSample.zdata{obj.s1Col}(targSampleIndex))^2;
1027                 obj.swpmeanx2 = obj.swpmeanx2 - (targSample.zdata{obj.s1Col}(targSampleIndex))/n2;                
1028                 obj.swpsumx2sq = obj.swpsumx2sq + (feederdf.zdata{obj.s1Col}(feederdfIndex))^2;
1029                 obj.swpmeanx2 = obj.swpmeanx2 + (feederdf.zdata{obj.s1Col}(feederdfIndex))/n2;
1030             end
1031             
1032             if feederdf == obj.s1
1033                 obj.swpsumx2sq = obj.swpsumx2sq -  (feederdf.zdata{obj.s1Col}(feederdfIndex))^2;
1034                 obj.swpmeanx2 = obj.swpmeanx2 - (feederdf.zdata{obj.s1Col}(feederdfIndex))/n2;                
1035                 obj.swpsumx2sq = obj.swpsumx2sq + (targSample.zdata{obj.s1Col}(targSampleIndex))^2;
1036                 obj.swpmeanx2 = obj.swpmeanx2 + (targSample.zdata{obj.s1Col}(targSampleIndex))/n2;   
1037             end
1038                     
1039             stdevx2 = (1/(n2-1) * obj.swpsumx2sq ...
1040                         - n2/(n2-1) * obj.swpmeanx2^2)^0.5;
1041                     
1042             obj.swstat1 = obj.stat1; %just copy over the target
1043             obj.swstat2 = stdevx2;
1044             
1045             if isnan(obj.swstat1) || isnan(obj.swstat2)
1046                 error('NaN obtained during swstat computation (perhaps there is missing data for an item?)');
1047             end            
1048           
1049         end
1050         
1051         
1052         
1053         %% minDiff() FUNCTION
1054         function cost = minDiff(obj,x1,x2)
1055             % cost function minimized by reduced differences on the statistic
1056             cost = (abs((x2-x1))^obj.exp)*obj.weight;
1057         end % minDiff
1058 
1059         
1060         %% maxDiff() FUNCTION
1061         function cost = maxDiff(obj,x1,x2)
1062             % cost function minimized by maximizing differences on the stat.  Many uses may want to use orderedMaxDiff, instead.
1063             cost = (-(abs((x2-x1))^obj.exp))*obj.weight;
1064         end       
1065         
1066         %% orderedMaxDiff() FUNCTION
1067         function cost = orderedMaxDiff(obj,x1,x2)
1068             % cost function minimized by maximizing a difference between the two groups that adheres to s1.stat < s2.stat
1069             if x1 < x2; cost = (-(abs((x2-x1))^obj.exp))*obj.weight;
1070             else cost = ((abs((x2-x1))^obj.exp))*obj.weight; end;
1071         end
1072         
1073         %% cost = acceptSwap(obj)
1074         function cost = acceptSwap(obj)
1075             %generic acceptSwap function
1076             obj.acceptsw();
1077             cost = acceptSwap@genericConstraint(obj);
1078         end
1079 
1080         %% cost = rejectSwap(obj)
1081         function cost = rejectSwap(obj)
1082             %generic reject function
1083             obj.rejectsw();
1084             cost = rejectSwap@genericConstraint(obj);
1085         end
1086         
1087         %% acceptSwapMeans() METHOD
1088         function acceptSwapMeans(obj)
1089             % updates variables when the proposed swap is accepted for
1090             % means
1091             
1092             if(isnan(obj.swstat1) && isnan(obj.swstat2))
1093                 %do nothing, this method does not need to update
1094             elseif(isnan(obj.swstat1) == false && isnan(obj.swstat2) == false)
1095               obj.stat1=obj.swstat1;
1096               obj.stat2=obj.swstat2;
1097          
1098               obj.swstat1 = NaN;
1099               obj.swstat2 = NaN;
1100             else
1101                 error('This should not have happened.  Maybe NaN in data?');
1102             end            
1103         end
1104 
1105         %% rejectSwapMeans() METHOD
1106         function rejectSwapMeans(obj)
1107             % updates variables when swap is rejected
1108               obj.swstat1 = NaN;
1109               obj.swstat2 = NaN;       
1110         end
1111         
1112          %% acceptSwapStdev() METHOD
1113         function acceptSwapStdev(obj)
1114             % updates variables when the proposed swap is accepted for
1115             % standard deviations
1116             
1117             if(isnan(obj.swstat1) && isnan(obj.swstat2))
1118                 %do nothing, this method does not need to update
1119             elseif(isnan(obj.swstat1) == false && isnan(obj.swstat2) == false)
1120                 
1121 %               disp('clearing swap memory');
1122               obj.stat1=obj.swstat1;
1123               obj.stat2=obj.swstat2;
1124          
1125           
1126               obj.sumx1sq = obj.swpsumx1sq;
1127               obj.meanx1 = obj.swpmeanx1;
1128               
1129               obj.sumx2sq = obj.swpsumx2sq;
1130               obj.meanx2 = obj.swpmeanx2;      
1131               
1132               obj.swpsumx1sq = NaN;
1133               obj.swpmeanx1 = NaN;
1134 
1135               obj.swpsumx2sq = NaN;
1136               obj.swpmeanx2 = NaN;   
1137               
1138               obj.swstat1 = NaN;
1139               obj.swstat2 = NaN;
1140             else
1141                 error('This should not have happened');
1142             end            
1143         end       
1144         
1145         %% rejectSwapStdev() METHOD
1146         function rejectSwapStdev(obj)
1147             % updates variables when the proposed swap is rejected
1148               
1149               obj.swpsumx1sq = NaN;
1150               obj.swpmeanx1 = NaN;
1151 
1152               obj.swpsumx2sq = NaN;
1153               obj.swpmeanx2 = NaN;   
1154               
1155               obj.swstat1 = NaN;
1156               obj.swstat2 = NaN;       
1157         end                 
1158     end
1159         
1160 
1161     
1162     
1163     methods (Static)
1164         %% p = softDistanceConstraintInputParser2Sample() STATIC METHOD
1165         function p = softDistanceConstraintInputParser2Sample()
1166             % generates an input parser with parameter / value pairs and validators for the constructor args
1167             %
1168             % See constructor help/doc for more details
1169             
1170             p = inputParser;
1171 
1172             %NOTE: though these technically are 'optional' according to
1173             %MATLAB's definition of what an input parser does, the fact
1174             %that most of their default values will fail validation in all cases
1175             %makes them de facto required parameters.
1176             
1177             p.addParamValue('sosObj','null',@(sosObj)strcmp(class(sosObj),'sos'));
1178             p.addParamValue('constraintType', 'null', ...
1179                 @(constraintType)any(strcmp({'soft'},constraintType)));
1180             p.addParamValue('fnc','null', ...
1181                  @(fnc)any(strcmp({'min' 'max','orderedMax'},fnc)));
1182             p.addParamValue('stat','null', ...
1183                  @(stat)any(strcmp({'mean' 'stdev'},stat)));
1184             p.addParamValue('sample1','null',@(sample1)strcmp(class(sample1),'sample'));
1185             p.addParamValue('sample2','null',@(sample2)strcmp(class(sample2),'sample'));
1186             p.addParamValue('s1ColName','',@(s1ColName)ischar(s1ColName));
1187             p.addParamValue('s2ColName','',@(s2ColName)ischar(s2ColName));
1188             p.addParamValue('paired','null',@(paired)islogical(paired));
1189             p.addParamValue('exponent',2,@(exponent)isnumeric(exponent));
1190             p.addParamValue('weight',1,@(weight)isnumeric(weight));
1191             p.addParamValue('name','noname',@(name)ischar(name));
1192             
1193         end % softDistanceConstraintInputParser()
1194         
1195         
1196         
1197         function p = softDistanceConstraintInputParser1Sample()
1198             % generates an input parser with parameter / value pairs and validators for the constructor args
1199             %
1200             % See constructor help/doc for more details
1201             
1202             p = inputParser;
1203 
1204             %NOTE: though these technically are 'optional' according to
1205             %MATLAB's definition of what an input parser does, the fact
1206             %that most of their default values will fail validation in all cases
1207             %makes them de facto required parameters.
1208             
1209             p.addParamValue('sosObj','null',@(sosObj)strcmp(class(sosObj),'sos'));
1210             p.addParamValue('constraintType', 'null', ...
1211                 @(constraintType)any(strcmp({'soft'},constraintType)));
1212             p.addParamValue('fnc','null', ...
1213                  @(fnc)any(strcmp({'match1SampleVal'},fnc)));
1214             p.addParamValue('stat','null', ...
1215                  @(stat)any(strcmp({'mean' 'stdev'},stat)));
1216             p.addParamValue('sample1','null',@(sample1)strcmp(class(sample1),'sample'));
1217             p.addParamValue('s1ColName','',@(s1ColName)ischar(s1ColName));
1218             p.addParamValue('exponent',2,@(exponent)isnumeric(exponent));
1219             p.addParamValue('weight',1,@(weight)isnumeric(weight));
1220             p.addParamValue('targVal','null',@(targVal)isnumeric(targVal));
1221             p.addParamValue('name','noname',@(name)ischar(name));
1222             
1223         end % softDistanceConstraintInputParser()
1224             
1225     end
1226         
1227     
1228     methods (Static, Access = private)
1229         
1230         %% p = parse2SampleContructorArgs(varargin) STATIC PRIVATE FUNCTION
1231         function p = parse2SampleConstructorArgs(varargin)
1232             % parses the constructor args
1233             %
1234             % See constructor help/doc for more info
1235                         
1236             varargin = varargin{1};
1237             p = softDistanceConstraint.softDistanceConstraintInputParser2Sample();
1238             p.parse(varargin{:});
1239         end
1240         
1241         %% p = parse1SampleContructorArgs(varargin) STATIC PRIVATE FUNCTION
1242         function p = parse1SampleConstructorArgs(varargin)
1243             % parses the constructor args
1244             %
1245             % See constructor help/doc for more info
1246                         
1247             varargin = varargin{1};
1248             p = softDistanceConstraint.softDistanceConstraintInputParser1Sample();
1249             p.parse(varargin{:});
1250         end
1251     end
1252     
1253 end % end class
1254