Home > src > softEntropyConstraint.m

# softEntropyConstraint

## PURPOSE

- soft entropy constraint object

## SYNOPSIS

This is a script file.

## DESCRIPTION

     - soft entropy constraint 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.

## CROSS-REFERENCE INFORMATION

This function calls:

- genericConstraint - parent class for constraints
- softConstraint - parent class for soft constraints
- softEntropyConstraint - soft entropy constraint object
- verbosePrint - master print function that regulates cmd line output based on verbosity flags

This function is called by:

- genericConstraint - parent class for constraints
- softEntropyConstraint - soft entropy constraint object

## SUBFUNCTIONS

- function obj = softEntropyConstraint(varargin)
- function cost = initCost(obj)
- function initEnt(obj)
- function cost = minEnt(obj,curEnt)
- function cost = maxEnt(obj,curEnt)
- function swCost = swapCost(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
- function swEntropy(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
- function cost = acceptSwap(obj)
- function cost = rejectSwap(obj)
- function plotDistribution(obj)

## SOURCE CODE

    0001 % - soft entropy constraint 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
    0024 classdef softEntropyConstraint < softConstraint
    0025     %% creates and supports soft entropy constraints
    0026     %
    0027     % Objects of this class measure the cost in terms of minimizing or
    0028     % maximizing a measure of entropy, or degree to which items have been
    0029     % randomly distributed on a particular dimension.
    0030     %
    0031     % Additional functionality is inherited from parent softConstraint
    0032     %
    0033     % The particular measure of entropy used in the class takes its
    0034     % inspiration from the standard version of Gibbs entropy:
    0035     %
    0036     %   -k*[Sigma(p)*ln(p)];
    0037     %
    0038     %  In this formula, k is a positive constant, and p is the probability
    0039     %  of choosing an item with a particular value on the
    0040     %  dimension of interest.  The sum of p*ln(p) is summed across all
    0041     %  values of the dimension of interest.
    0042     %
    0043     %  In this standard version of entropy, the values of Entropy increase
    0044     %  as the 'randomness' of the distribution of items increases.  That
    0045     %  is, the highest values of entropy are obtained when a randomly
    0046     %  selected item has an equal chance of having any one of the values on
    0047     %  the dimension of interest.  This also corresponds to having a
    0048     %  uniform distribution of items across values on the dimension of
    0049     %  interest.  In contrast, when items are not randomly distributed in
    0050     %  this uniform fashion, lower values of entropy are obtained.  The
    0051     %  lowest possible entropy would occur when all items shared the same
    0052     %  value.
    0053     %
    0054     %  Maximizing entropy thus provides one means of randomly distributing
    0055     %  items across a range of values, whereas minimizing entropy can serve
    0056     %  to force units to all take on the same value.
    0057     %
    0058     %  Though the standard Gibbs entropy formula has the desirable
    0059     %  characteristics outlined above, for present purposes in one main
    0060     %  respect: the domain of the results of the entropy formula are not
    0061     %  well constrainted to a particular range.  This problem is due to two
    0062     %  main aspects of the formula.  First, if a particular value of a
    0063     %  dimension literally has no item representing it, its probability
    0064     %  would be zero, and ln(0) = -Inf.  Second, the amount of entropy in
    0065     %  the standard formula varies as a function of the number of items.
    0066     %  Neither of these characteristics are desirable in the present
    0067     %  context where it is useful to have the standard output of a cost
    0068     %  function consistently have standardized values in the range -1 - 1
    0069     %   (or at least within an order of magnitude thereof).
    0070     %
    0071     %  This additional desirable characteristic resulted in the development
    0072     %  of a new 'entropy' formula, as follows:
    0073     %
    0074     %     ent = (-1*Sigma(1:N)[p*log((p+1)')/Sum(N)] -1/(N)*ln((1/N)+1)) /
    0075     %           (ln(2)/N - 1/N * ln((1/N) +1));
    0076     %
    0077     %  Don't let the complexity of this formula scare you though, as in
    0078     %  principle it's quite straightforward.  First, to avoid the case
    0079     %  where values with p=0 generated -Inf, it was necessary to add some
    0080     %  constant value to the ln() term.  A value of 1 was selected as this
    0081     %  results in only positive numbers being generated from the equation,
    0082     %  since log(1) == 0 and log(1+positiveN) > 0.  Next, this equation is
    0083     %  divided by N, the number of different possible items, which renders
    0084     %  the output of the equation relatively insensitive to the number of
    0085     %  items for which entropy is being calculated.  The theoretical
    0086     %  minimum entropy for this particular dataset is then subtracted away
    0087     %  from the result of hte previous step (1/N*ln((1/N)+1) corresponding
    0088     %  to the case where all items are distributed uniformally across
    0089     %  values) giving the equation a bound a bound of zero when maximum
    0090     %  entropy is acheived.  The result of this entire calculation
    0091     %  subsequently serves as the dividend in a division which
    0092     %  standardizes the equation so that the maximum possible entropy that
    0093     %  could result will always have a value of 1.0, thus bounding the
    0094     %  formula to the 0-1 domain.  Finally, to conserve the property that
    0095     %  entropy is maximized as randomness increases, the result of the
    0096     %  previous calculation is multiplied by -1, thus yielding an entropy
    0097     %  formula that is bounded at [-1, 0], where -1 corresponds to the
    0098     %  least randomly distributed state, and 0 corresponds to the most
    0099     %  randomly distributed state.
    0100     %
    0101     %
    0102     %     s1 % sample entropy is calculated on
    0103     %     s1Col   % data column entropy is calculated on
    0104     %     nbin    % number of bins to subdivide data into
    0105     %     pdSpread % type of spread of scores desired.  Options currently are 'sample' and 'allItems'.
    0106     %     ent % current entropy value
    0107     %     swEnt % swap entropy value
    0108     %     comparison %handle to comparison calculator (to minimize or max entropy)
    0109     %     initStats % handle to method that initializes stats
    0110     %     swStats % handle to method that calculates swap stats
    0111     %     scores % copy of scores in the column entropy is calucated for.
    0112     %     swScores % copy of scores if a swap occured
    0113     %     bins % mid-point of bins data is divided into
    0114     %     swBins % mid-point of bins data is divided into after a swap
    0115     %     pd % probability distribution of scores
    0116     %     swpd % probability distribution of swap socres
    0117     %     minScore % min score in set
    0118     %     swMinScore % min score in swap set
    0119     %     maxScore % max score in set
    0120     %     swMaxScore % max score in swap set.
    0121     %
    0122     %PROPERTIES
    0123     %     s1 % sample entropy is calculated on
    0124     %     s1Col   % data column entropy is calculated on
    0125     %     nbin    % number of bins to subdivide data into
    0126     %     pdSpread % type of spread of scores desired.  Options currently are 'sample' and 'allItems'.
    0127     %     ent % current entropy value
    0128     %     swEnt % swap entropy value
    0129     %     comparison %handle to comparison calculator (to minimize or max entropy)
    0130     %     initStats % handle to method that initializes stats
    0131     %     swStats % handle to method that calculates swap stats
    0132     %     scores % copy of scores in the column entropy is calucated for.
    0133     %     swScores % copy of scores if a swap occured
    0134     %     bins % mid-point of bins data is divided into
    0135     %     swBins % mid-point of bins data is divided into after a swap
    0136     %     pd % probability distribution of scores
    0137     %     swpd % probability distribution of swap socres
    0138     %     minScore % min score in set
    0139     %     swMinScore % min score in swap set
    0140     %     maxScore % max score in set
    0141     %     swMaxScore % max score in swap set.
    0142     %
    0143     %PROPERTIES (Constant)
    0144     %    s2 = NaN; % s2 = NaN; for consistency with other methods, there is a swap2 property, but it is just set to a null-like value.
    0145     %   smallVal = 0.00000001 % small value added to avoid case where a value is exactly equal the min or max range of the distribution
    0146     %
    0147     %METHODS
    0148     %   obj = softEntropyConstraint(varargin) % CONSTRUCTOR
    0149     %   cost = initCost() METHOD
    0150     %   initEnt() % calculates the initial entropy value for the scores in sample1<sample1Col>
    0151     %   cost = minEnt(curEnt) % calculates the cost associated with curEnt, when minimal Entropy is desired
    0152     %   cost = maxEnt(curEnt)  % calculates the cost associated with curEnt, when maximal entropy is desired
    0153     %   swCost = swapCost(targSample,targSampleIndex, feederdf,feederdfIndex)  % Calculates the new cost if items from targSample and feederdf were swapped.
    0154     %   swEntropy(obj,targSample,targSampleIndex, feederdf,feederdfIndex)  % Calculates the new (swap) entropy if items from targSample and feederdf were swapped.
    0155     %   cost = acceptSwap()  % alter internal variables to reflect accepted swap
    0156     %   cost = rejectSwap() % set internal variables to reflect rejection of proposed swap
    0157     %   plotDistribution() % plots the probability distribution used to calculate entropy
    0158
    0159
    0160     %% Properties
    0161     properties
    0162         s1 % sample entropy is calculated on
    0163         s1Col   % data column entropy is calculated on
    0164         s1ColName %name of column entropy is being calculated on
    0165         nbin    % number of bins to subdivide data into
    0166         pdSpread % type of spread of scores desired.  Options currently are 'sample' and 'allItems'.
    0167         ent % current entropy value
    0168         swEnt % swap entropy value
    0169         comparison %handle to comparison calculator (to minimize or max entropy)
    0170         initStats % handle to method that initializes stats
    0171         swStats % handle to method that calculates swap stats
    0172         scores % copy of scores in the column entropy is calucated for.
    0173         swScores % copy of scores if a swap occured
    0174         bins % mid-point of bins data is divided into
    0175         swBins % mid-point of bins data is divided into after a swap
    0176         pd % probability distribution of scores
    0177         swpd % probability distribution of swap socres
    0178         minScore % min score in set
    0179         swMinScore % min score in swap set
    0180         maxScore % max score in set
    0181         swMaxScore % max score in swap set.
    0182         figHandle % handle to figure for plotting entropy
    0183     end
    0184
    0185     %% Properties (Constant)
    0186     properties (Constant)
    0187         s2 = NaN; % s2 = NaN; for consistency with other methods, there is a swap2 property, but it is just set to a null-like value.
    0188         smallVal = 0.00000001 % small value added to avoid case where a value is exactly equal the min or max range of the distribution
    0189     end
    0190
    0191
    0192     methods
    0193         %% obj = softEntropyConstraint(varargin) CONSTRUCTOR
    0194         function obj = softEntropyConstraint(varargin)
    0195             % Constructs a softEntropy constraint object
    0196             %
    0197             % PARAMETERS:
    0198             % REQUIRED:
    0199             %   'sosObj'/sos object - the SOS object the constraint will be linked to, and which contains the sample the constraint operates on.
    0200             %   'constraintType'/'soft' - the type of contraint - must be 'soft'
    0201             %   'fnc'/'minEnt'|'maxEnt' the entropy cost to calculate.
    0202             %   'sample1'/sample - the first sample
    0203             %   's1ColName'/string - name of column in 1st sample
    0204             %   'pSpread'/'sample'|'allItems'  Should entropy be maximized relative to sample items only, or to the theoretical min and max values in the population?
    0205             %
    0206             %
    0207             % OPTIONAL:
    0208             %   'nbin'/integer - number of bins to divide data into for
    0209             %       purposes of calculating probability distribution.  Defaults
    0210             %       to number of items in the sample.  Must be greater than 2
    0211             %       and <= number of items in the sample
    0212             %   'exponent'/numeric - defaults to 2 (quadratic difference)
    0213             %   'weight'/numeric - defaults to 1 (equal weighting of all soft costs)
    0214
    0215             p = inputParser;
    0216
    0217             p.addParamValue('sosObj','null',@(sosObj)strcmp(class(sosObj),'sos'));
    0218             p.addParamValue('constraintType', 'null', ...
    0219                 @(constraintType)any(strcmp({'soft'},constraintType)));
    0220             p.addParamValue('fnc','null', ...
    0221                  @(fnc)any(strcmp({'minEnt' 'maxEnt'},fnc)));
    0222             p.addParamValue('nbin',2,@(nbin)validateattributes(nbin, {'numeric'}, ...
    0223                 {'scalar', 'integer', 'positive', '>', 0}));
    0224             p.addParamValue('pdSpread','null', ...
    0225                  @(pdSpread)any(strcmp({'sample' 'allItems'},pdSpread)));
    0226             p.addParamValue('sample1','null',@(sample1)strcmp(class(sample1),'sample'));
    0227             p.addParamValue('s1ColName','',@(s1ColName)ischar(s1ColName));
    0228             p.addParamValue('exponent',2,@(exponent)isnumeric(exponent));
    0229             p.addParamValue('weight',1,@(weight)isnumeric(weight));
    0230             p.addParamValue('name','noname',@(name)ischar(name));
    0231
    0232             p.parse(varargin{:});
    0233
    0234             % check additional constraints on values submitted to the
    0235             % constructor
    0236             if(p.Results.sosObj.containsSample(p.Results.sample1) == false)
    0237                 error('Cannot create soft distance constraint: sos Object does not contain the sample1');
    0238             end
    0239
    0240             col1 = p.Results.sample1.colName2colNum(p.Results.s1ColName);
    0241             if(col1 == -1)
    0242                 error('Specified column name not found in sample1');
    0243             end
    0244
    0245             if(strcmp(p.Results.sample1.format{col1},'%f') == 0)
    0246                 error('Specified column is not of numeric (%f) format, so cannot use as hard bound');
    0247             end
    0248
    0249             if any(strcmp(p.UsingDefaults,'nbin'))
    0250                 obj.nbin = p.Results.sample1.n;
    0251             else
    0252                 if p.Results.nbin <= p.Results.sample1.n
    0253                     obj.nbin = p.Results.nbin;
    0254                 else
    0255                     error('Number of bins must be <= number of observations in the sample');
    0256                 end
    0257             end
    0258
    0259             if obj.nbin < 2
    0260                 error ('There must be at least 2 bins in the entropy calculation');
    0261             end
    0262
    0263
    0264             %assign the appropriate handles for the calculation.
    0265             obj.sosObj = p.Results.sosObj;
    0266             obj.constraintType = p.Results.constraintType;
    0267             obj.fnc = p.Results.fnc;
    0268             obj.weight = p.Results.weight;
    0269             obj.s1 = p.Results.sample1;
    0270             obj.s1ColName = p.Results.s1ColName;
    0271             obj.s1Col = col1;
    0272             obj.exp = p.Results.exponent;
    0273
    0274             obj.initStats = @obj.initEnt;
    0275             obj.swStats = @obj.swEntropy;
    0276
    0277             if(strcmp(p.Results.fnc,'minEnt'))
    0278                 obj.comparison = @obj.minEnt;
    0279             elseif(strcmp(p.Results.fnc,'maxEnt'))
    0280                obj.comparison = @obj.maxEnt;
    0281             else
    0282                 error('Entropy <fnc> must be either "minEnt" or "maxEnt"');
    0283             end
    0284
    0285             if(any(strcmp({'sample' 'allItems'},p.Results.pdSpread)))
    0286                 obj.pdSpread = p.Results.pdSpread;
    0287             else
    0288                 error('pdSpread must be either "sample" or "allItems"');
    0289             end
    0290
    0291             obj.cost = NaN;
    0292             obj.swCost = NaN;
    0293
    0294             % add the name and the label
    0295             obj.label = [obj.constraintType,'_',obj.fnc,...
    0296                     '_pd_',obj.pdSpread,'_nbin',num2str(obj.nbin),'_',...
    0297                     obj.s1.name,'_',...
    0298                     obj.s1ColName,'_w',...
    0299                     num2str(obj.weight),'_e',num2str(obj.exp)];
    0300             if any(strcmp(p.UsingDefaults,'name'))
    0301                 obj.name = obj.label;
    0302             else
    0303                  obj.name = p.Results.name;
    0304             end
    0305
    0306
    0307             verbosePrint('Soft Entropy Constraint has been created', ...
    0308                     'softEntropyConstraint_Constructor_endObjCreation');
    0309         end % constructor
    0310
    0311
    0312         %% cost = initCost() METHOD
    0313         function cost = initCost(obj)
    0314             % Calculates, saves, and returns the cost value for the current items in the sample.
    0315
    0316             %init the stats, then compare them.  Return the calculated cost
    0317             obj.initStats();
    0318             cost = obj.comparison(obj.ent);
    0319
    0320             obj.swEnt = NaN;
    0321             obj.swScores = NaN;
    0322             obj.swpd = NaN;
    0323             obj.swMinScore = NaN;
    0324             obj.swMaxScore = NaN;
    0325
    0326             obj.cost = cost;
    0327
    0328         end %initCost
    0329
    0330         %% initEnt() METHOD
    0331         function initEnt(obj)
    0332             % calculates the initial entropy value for the scores in
    0333             % sample1<sample1Col>
    0334
    0335             obj.scores = (obj.s1.zdata{obj.s1Col})';
    0336
    0337
    0338             % get min and max values for the probability distribution:
    0339             %this procedure is different depending on whether we are using
    0340             %the sample as the range or allItems
    0341
    0342             if strcmp(obj.pdSpread,'sample')
    0343                 minVal = min(obj.scores);
    0344                 maxVal = max(obj.scores);
    0345             elseif strcmp(obj.pdSpread,'allItems')
    0346                 %check in the population
    0347                 minVal = min([obj.s1.population.zdata{obj.s1Col}]);
    0348                 maxVal = max([obj.s1.population.zdata{obj.s1Col}]);
    0349
    0350                 %check in all samples associated with the population (this
    0351                 %includes the original one
    0352                 for i=1:length(obj.s1.population.samples)
    0353                     minVal = min([minVal; ...
    0354                         obj.s1.population.samples(i).zdata{obj.s1Col}]);
    0355                     maxVal = max([maxVal; ...
    0356                         obj.s1.population.samples(i).zdata{obj.s1Col}]);
    0357                 end
    0358             else
    0359                 error('Unsupported pdSpread in initEnt.  pdSpread must be "sample" or "allItems"');
    0360             end
    0361
    0362             % add a smallVal to avoid boundary cases
    0363             minVal = minVal - obj.smallVal;
    0364             maxVal = maxVal + obj.smallVal;
    0365
    0366             %now have min and max values, determine size of each bin:
    0367             obj.minScore = minVal;
    0368             obj.maxScore = maxVal;
    0369
    0370             spread  = maxVal - minVal;
    0371             binSize = spread/obj.nbin;
    0372
    0373             %smallVal added so that the last score is included (simulating
    0374             %<= rather than <
    0375             obj.bins = (minVal+0.5*binSize):binSize:(maxVal-0.5*binSize)+obj.smallVal;
    0376
    0377             % scale values to 0-1
    0378             obj.pd = hist(obj.scores,obj.bins)/length(obj.scores);
    0379
    0380
    0381             % calculate entropy
    0382             obj.ent = obj.pd*log((obj.pd+1)')/length(obj.pd);
    0383             obj.ent = -1* (obj.ent - 1/length(obj.pd)*log(1/length(obj.pd)+1)) / ...
    0384                 (1*log(2)/length(obj.pd) - 1/length(obj.pd)*log(1/length(obj.pd)+1));
    0385
    0386         end % initEnt
    0387
    0388         %% cost = minEnt(curEnt) METHOD
    0389         function cost = minEnt(obj,curEnt)
    0390             % calculates the cost associated with curEnt, when minimal entropy is desired
    0391             cost = -(abs((curEnt))^obj.exp)*obj.weight;
    0392         end
    0393
    0394         %% cost = maxEnt(curEnt)
    0395         function cost = maxEnt(obj,curEnt)
    0396             % calculates the cost associated with curEnt, when maximal entropy is desired
    0397             cost = (abs((curEnt))^obj.exp)*obj.weight;
    0398         end
    0399
    0400         %%  swCost(targSample,targSampleIndex, feederdf,feederdfIndex) METHOD
    0401         function swCost = swapCost(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
    0402             % Calculates the new cost if items from targSample and feederdf were swapped.
    0403             %
    0404             %By definition, if this method is called it means that at least
    0405             % one of the two swap objects is implicated in this function
    0406             %
    0407             %PARAMETERS:
    0408             %   targSample - the target sample (i.e., the object that will call it's swapSample() method if a swap later occurs).
    0409             %   targSampleIndex - row index of item to swap
    0410             %   feederdf - dataframe (sample/pop) containin the other item to swap
    0411             %   feederdfIndex - row index of item to swap.
    0412
    0413
    0414            if (obj.s1 ~= targSample && obj.s1 ~= feederdf)
    0415                error('swCost called, but no sample part of this cost function');
    0416            end
    0417
    0418            % update the stats, then calculate cost of new stats
    0419            obj.swStats(targSample,targSampleIndex, feederdf,feederdfIndex);
    0420
    0421            swCost = obj.comparison(obj.swEnt);
    0422
    0423            obj.swCost = swCost;
    0424
    0425         end
    0426
    0427         %% swEnt(targSample,targSampleIndex,feederdf,feederdfIndex) METHOD
    0428         function swEntropy(obj,targSample,targSampleIndex, feederdf,feederdfIndex)
    0429            % Calculates the new (swap) entropy if items from targSample and feederdf were swapped.
    0430            %
    0431            % Inputs are the same as for swapCost()
    0432
    0433
    0434
    0435            obj.swScores = obj.scores;
    0436            obj.swpd = obj.pd;
    0437            obj.swMinScore = obj.minScore;
    0438            obj.swMaxScore = obj.maxScore;
    0439            obj.swBins = obj.bins;
    0440
    0441            % re-extract the raw, unnormalized entropy value from the
    0442            % current entropy
    0443            rawEnt = obj.ent;
    0444            rawEnt = rawEnt *(1*log(2)/length(obj.pd) - 1/length(obj.pd)*log(1/length(obj.pd)+1));
    0445            rawEnt = rawEnt*-1;
    0446            rawEnt = rawEnt + 1/length(obj.pd)*log(1/length(obj.pd)+1);
    0447
    0448            spread  = obj.swMaxScore - obj.swMinScore;
    0449            binSize = spread/obj.nbin;
    0450
    0451            recalculate = false;
    0452
    0453             if targSample == obj.s1
    0454
    0455                 obj.swScores(targSampleIndex) = ...
    0456                     feederdf.zdata{obj.s1Col}(feederdfIndex);
    0457
    0458                 % if the min and max bounds may change and we're only
    0459                 % calculating the entropy of the sample and letting the
    0460                 % range of the sample vary dynamically, recalculate entropy
    0461                 % from scratch
    0462                 if(strcmp(obj.pdSpread,'sample') && ...
    0463                     (lt(obj.swScores(targSampleIndex)-obj.smallVal, obj.minScore) || ...
    0464                      gt(obj.swScores(targSampleIndex)+obj.smallVal, obj.maxScore) || ...
    0465                      obj.scores(targSampleIndex)-obj.smallVal == obj.minScore || ...
    0466                      obj.scores(targSampleIndex)+obj.smallVal == obj.maxScore))
    0467
    0468                     recalculate = true;
    0469                 else
    0470                    %we can do a local update
    0471
    0472                    found = false;
    0473
    0474                    % find the old bin:
    0475                    for i=1:(length(obj.swBins))
    0476
    0477
    0478                        if (obj.scores(targSampleIndex) >= obj.swBins(i) - 0.5*binSize ...
    0479                                && obj.scores(targSampleIndex) < obj.swBins(i)+0.5*binSize)
    0480
    0481                            %remove the old score
    0482                            rawEnt = rawEnt - obj.swpd(i)*log(obj.swpd(i)+1)/length(obj.pd);
    0483                            obj.swpd(i) = obj.swpd(i) - 1/length(obj.scores);
    0484                            rawEnt = rawEnt + obj.swpd(i)*log(obj.swpd(i)+1)/length(obj.pd);
    0485                            found = true;
    0486                            break;
    0487                        end
    0488                    end
    0489
    0490                    if found == false
    0491                        error('unable to find the swap bin');
    0492                    end
    0493
    0494                    found = false;
    0495
    0496                     for i=1:(length(obj.swBins))
    0497                        if obj.swScores(targSampleIndex) >= obj.swBins(i) - 0.5*binSize ...
    0498                                && obj.swScores(targSampleIndex) < obj.swBins(i) + 0.5*binSize
    0499                            %add the new score
    0500
    0501                            rawEnt = rawEnt - obj.swpd(i)*log(obj.swpd(i)+1)/length(obj.pd);
    0502                            obj.swpd(i) = obj.swpd(i) + 1/length(obj.scores);
    0503                            rawEnt = rawEnt + obj.swpd(i)*log(obj.swpd(i)+1)/length(obj.pd);
    0504
    0505                             found = true;
    0506                            break;
    0507                        end
    0508                     end
    0509
    0510                    if found == false
    0511                        error('unable to find the swap bin');
    0512                    end
    0513
    0514                 end
    0515              end
    0516
    0517              if feederdf == obj.s1
    0518
    0519                 obj.swScores(feederdfIndex) =  ...
    0520                     targSample.zdata{obj.s1Col}(targSampleIndex);
    0521
    0522                 % if using the sample's spread and the new score could
    0523                 % exceed the existing bounds, recalculate from sractch
    0524                 if(strcmp(obj.pdSpread,'sample') && ...
    0525                         (lt(obj.swScores(feederdfIndex)-obj.smallVal, obj.minScore) || ...
    0526                          gt(obj.swScores(feederdfIndex)+obj.smallVal, obj.maxScore) || ...
    0527                          obj.scores(feederdfIndex)-obj.smallVal == obj.minScore || ...
    0528                          obj.scores(feederdfIndex)+obj.smallVal == obj.maxScore))
    0529
    0530                     recalculate = true;
    0531                 else
    0532                    %we can do a local update
    0533
    0534                    found = false;
    0535
    0536
    0537
    0538                    % find the old bin:
    0539                    for i=1:(length(obj.swBins))
    0540                        if obj.scores(feederdfIndex) >= obj.swBins(i) - 0.5*binSize ...
    0541                                && obj.scores(feederdfIndex) < obj.swBins(i) + 0.5*binSize
    0542                            %remove the old score
    0543                            rawEnt = rawEnt - obj.swpd(i)*log(obj.swpd(i)+1)/length(obj.pd);
    0544                            obj.swpd(i) = obj.swpd(i) - 1/length(obj.scores);
    0545                            rawEnt = rawEnt + obj.swpd(i)*log(obj.swpd(i)+1)/length(obj.pd);
    0546
    0547                            found = true;
    0548                            break;
    0549                        end
    0550                    end
    0551
    0552                    if found == false
    0553                        error('Unable to find the old bin');
    0554                    end
    0555
    0556                    found = false;
    0557                     for i=1:(length(obj.swBins))
    0558                        if obj.swScores(feederdfIndex) >= obj.swBins(i) - 0.5*binSize ...
    0559                                && obj.swScores(feederdfIndex) < obj.swBins(i) + 0.5*binSize
    0560                            %add the new score
    0561                            rawEnt = rawEnt - obj.swpd(i)*log(obj.swpd(i)+1)/length(obj.pd);
    0562                            obj.swpd(i) = obj.swpd(i) + 1/length(obj.scores);
    0563                            rawEnt = rawEnt + obj.swpd(i)*log(obj.swpd(i)+1)/length(obj.pd);
    0564
    0565                            found = true;
    0566                            break;
    0567                        end
    0568                     end
    0569
    0570                     if found == false
    0571                        error('Unable to find an appropriate swbin for the data');
    0572                     end
    0573
    0574                 end
    0575             end
    0576
    0577             % a local update is possible if neither of the swap scores
    0578             % exceeds the current spread.
    0579
    0580             %recalculate = true;
    0581
    0582             if recalculate == true
    0583                 %re-derive all measures
    0584                 minVal = min(obj.swScores);
    0585                 maxVal = max(obj.swScores);
    0586
    0587                 minVal = minVal - obj.smallVal;
    0588                 maxVal = maxVal + obj.smallVal;
    0589
    0590                 %now have min and max values, determine size of each bin:
    0591                 obj.swMinScore = minVal;
    0592                 obj.swMaxScore = maxVal;
    0593
    0594                 spread  = maxVal - minVal;
    0595                 binSize = spread/obj.nbin;
    0596                 % remove the small -1.0e-14 value to avoid rounding errors
    0597                 % that result in non-overlapping bins.  So long as this
    0598                 % number is < e-16 then this appears to prevent these
    0599                 % rounding errors.  And so long as this number is
    0600                 % substantially smaller than smallVal it will not interfere
    0601                 % with the actual bin sizes (currently smallVal is <
    0602                 % 1.0-e10)
    0603                 obj.swBins = (minVal + 0.5*binSize):binSize:(maxVal - 0.5*binSize)+obj.smallVal;
    0604
    0605                 for i=1:length(obj.swBins)
    0606                     obj.swBins(i) = obj.swBins(i) - 1.0e-14;
    0607                 end
    0608
    0609                 obj.swpd = hist(obj.swScores,obj.swBins)/length(obj.scores);
    0610                 obj.swEnt = obj.swpd*log((obj.swpd+1)')/length(obj.pd);
    0611                 obj.swEnt = -1* (obj.swEnt - 1/length(obj.pd)*log(1/length(obj.pd)+1)) / ...
    0612                     (1*log(2)/length(obj.pd) - 1/length(obj.pd)*log(1/length(obj.pd)+1));
    0613
    0614             else
    0615                 % cannot exceed bounds because min and max were taken from
    0616                 % the population.  Proceed with local update
    0617                 obj.swEnt = -1* (rawEnt - 1/length(obj.pd)*log(1/length(obj.pd)+1)) / ...
    0618                     (1*log(2)/length(obj.pd) - 1/length(obj.pd)*log(1/length(obj.pd)+1));
    0619
    0620             end
    0621
    0622             if isnan(obj.swEnt)
    0623                 error('NaN obtained during swEnt computation (perhaps there is missing data for an item?)');
    0624             end
    0625
    0626         end %swEntropy()
    0627
    0628
    0629        %% cost = acceptSwap() METHOD
    0630         function cost = acceptSwap(obj)
    0631             % alter internal variables to reflect accepted swap
    0632
    0633
    0634             if isnan(obj.swEnt)
    0635                 %do nothing, no need to swap
    0636             else
    0637                 obj.ent = obj.swEnt;
    0638                 obj.swEnt = NaN;
    0639                 obj.scores = obj.swScores;
    0640                 obj.swScores = NaN;
    0641
    0642                obj.pd = obj.swpd;
    0643                obj.swpd = NaN;
    0644
    0645                obj.minScore = obj.swMinScore;
    0646                obj.swMinScore = NaN;
    0647                obj.maxScore = obj.swMaxScore;
    0648                obj.swMaxScore = NaN;
    0649
    0650                obj.bins = obj.swBins;
    0651                obj.swBins = NaN;
    0652
    0653             end
    0654
    0655             cost = acceptSwap@genericConstraint(obj);
    0656         end
    0657
    0658         %% cost = rejectSwap() METHOD
    0659         function cost = rejectSwap(obj)
    0660             % set internal variables to reflect rejection of proposed swap
    0661             obj.swEnt = NaN;
    0662             obj.swScores = NaN;
    0663             obj.swpd = NaN;
    0664             obj.swMinScore = NaN;
    0665             obj.swMaxScore = NaN;
    0666
    0667             obj.swBins = NaN;
    0668
    0669             cost = rejectSwap@genericConstraint(obj);
    0670         end
    0671
    0672         %% plotDistribution() METHOD
    0673         function plotDistribution(obj)
    0674             % plots the probability distribution used to calculate entropy
    0675
    0676             % make the entropy figure the active figure
    0677             if isempty(obj.figHandle)
    0678                 obj.figHandle = figure('Name',['Entropy - ',obj.s1.name, ...
    0679                     '|',obj.s1ColName],'NumberTitle','off');
    0680
    0681             end
    0682
    0683             figure(obj.figHandle);
    0684
    0685             if isempty(obj.cost) == false && isnan(obj.cost) == false %the constraints have not been initialized
    0686                 clf();
    0687
    0688
    0689                 if strcmp(obj.pdSpread,'sample')
    0690                     minVal = min(obj.s1.data{obj.s1Col});
    0691                     maxVal = max(obj.s1.data{obj.s1Col});
    0692                 elseif strcmp(obj.pdSpread,'allItems')
    0693                     %check in the population
    0694                     minVal = min([obj.s1.population.data{obj.s1Col}]);
    0695                     maxVal = max([obj.s1.population.data{obj.s1Col}]);
    0696
    0697                     %check in all samples associated with the population (this
    0698                     %includes the original one
    0699                     for i=1:length(obj.s1.population.samples)
    0700                         minVal = min([minVal; ...
    0701                             obj.s1.population.samples(i).data{obj.s1Col}]);
    0702                         maxVal = max([maxVal; ...
    0703                             obj.s1.population.samples(i).data{obj.s1Col}]);
    0704                     end
    0705                 else
    0706                     error('Unsupported pdSpread in initEnt.  pdSpread must be "sample" or "allItems"');
    0707                 end
    0708
    0709                 % add a smallVal to avoid boundary cases
    0710                 minVal = minVal - obj.smallVal;
    0711                 maxVal = maxVal + obj.smallVal;
    0712
    0713                 spread  = maxVal - minVal;
    0714                 binSize = spread/obj.nbin;
    0715
    0716                 %smallVal added so that the last score is included (simulating
    0717                 %<= rather than <
    0718                 tmpbins = (minVal+0.5*binSize):binSize:(maxVal-0.5*binSize)+obj.smallVal;
    0719
    0720                 % scale values to 0-1
    0721                 tmppd = hist(obj.s1.data{obj.s1Col},obj.nbin)/length(obj.s1.data{obj.s1Col});
    0722
    0723                 bar(tmpbins,tmppd,'hist');
    0724                 ylim([0 1]);
    0725                 title(['Entropy - ',obj.s1.name, ...
    0726                     '|',obj.s1ColName]);
    0727                  xlabel([obj.s1ColName, ]);%' - each bar is centered on bin mid-point']);
    0728                 ylabel('p(item from a particular bin)');
    0729             else
    0730                 clf();
    0731                 ylim([0 1]);
    0732                 title(['Entropy - ',obj.s1.name, ...
    0733                     '|',obj.s1ColName]);
    0734                 xlabel([obj.s1ColName, ]);%' - each bar is centered on bin mid-point']);
    0735                 ylabel('p(item from a particular bin)');
    0736                 text(0.5,0.5,'Cost must be initialized before data is plotted', ...
    0737                     'HorizontalAlignment','center');
    0738             end
    0739
    0740
    0741
    0742
    0743         end
    0744
    0745     end
    0746
    0747 end
    0748