fn_sffs.hpp

Go to the documentation of this file.

// Copyright (C) 2011 the authors listed below
// http://ecocpak.sourceforge.net
//
// Authors:
// - Dimitrios Bouzas (bouzas at ieee dot org)
// - Nikolaos Arvanitopoulos (niarvani at ieee dot org)
// - Anastasios Tefas (tefas at aiia dot csd dot auth dot gr)
//
// This file is part of the ECOC PAK C++ library. It is
// provided without any warranty of fitness for any purpose.
//
// You can redistribute this file and/or modify it under
// the terms of the GNU Lesser General Public License (LGPL)
// as published by the Free Software Foundation, either
// version 3 of the License or (at your option) any later
// version.
// (see http://www.opensource.org/licenses for more info)




double
criterion_fqmi
  (
  const vector<ClassData>& classes_vector,
  const ucolvec& indexPtr,
  const u32 curSubsetSize,
  const double sigma
  )
  {
  // print header
  #ifdef NCURSES_OUTPUT
  mvaddstr(3, 0, "Computing FQMI:      ");
  refresh();
  #endif

  if((classes_vector.size() - curSubsetSize) == 0)
    {
    return 0.0;
    }

  // data samples
  mat samples;

  // sample labels
  ucolvec labels;

  // compute complement subset of current best subset
  ucolvec cSet = complement_sffs(indexPtr, curSubsetSize);

  // create first bipartition
  for(u32 i = 0; i < curSubsetSize; i++)
    {
    samples = join_cols(samples, classes_vector[indexPtr[i]].Data());
    labels =
      join_cols
        (
        labels,
        ones<ucolvec>(classes_vector[indexPtr[i]].Data().n_rows)
        );
    }

  // create second bipartition
  for(u32 i = 0; i < cSet.n_rows; i++)
    {
    samples = join_cols(samples, classes_vector[cSet[i]].Data());
    labels =
      join_cols
        (
        labels,
        2 * ones<ucolvec>(classes_vector[cSet[i]].Data().n_rows)
        );
    }

  return fqmi(samples, labels, sigma);
  }



double
criterion_fldr
  (
  const vector<ClassData>& classes_vector,
  const ucolvec& indexPtr,
  const u32 curSubsetSize
  )
  {
  #ifdef NCURSES_OUTPUT
  mvaddstr(3, 0, "Computing FLDR:     ");
  refresh();
  #endif

  if((classes_vector.size() - curSubsetSize) == 0)
    {
    return 0.0;
    }

  // size of current subset is 1
  if(curSubsetSize == 1)
    {
    // if the number of samples of the current subset is less than 2
    if(classes_vector[indexPtr[0]].Samples() < 2)
      {
      return 0.0;
      }
    }

  // find the complement set of indexPtr
  ucolvec cSet = complement_sffs(indexPtr, curSubsetSize);

  // number of classes and current subset size differ only by 1
  if((classes_vector.size() - curSubsetSize ) == 1)
    {
    // if the number of samples of the complement subset is less than 2
    if(classes_vector[cSet[0]].Samples() < 2)
      {
      return 0.0;
      }
    }

  // ClassData vector for partition 1
  vector<ClassData> newClassVec1;

  // push back the classes that belong to one set
  for(u32 i = 0; i < curSubsetSize; i++)
    {
    newClassVec1.push_back(classes_vector[indexPtr[i]]);
    }

  // ClassData vector for partition 2
  vector<ClassData> newClassVec2;

  // push back classes that belong to the complement set
  for(u32 i = 0; i < classes_vector.size() - curSubsetSize; i++)
    {
    newClassVec2.push_back(classes_vector[cSet[i]]);
    }

  // create data matrices that correspond to the two bi-partitions
  mat first_bipartition = create_samples_matrix(newClassVec1);
  mat second_bipartition = create_samples_matrix(newClassVec2);

  return fldr(first_bipartition, second_bipartition);
  }



double
criterion
  (
  const vector<ClassData>& class_vector,
  const ucolvec& indexPtr,
  const u32 curSubSetSize,
  const u32 criterion_option,
  const double sigma
  )
  {
  // engage according to user enter criterion option
  switch(criterion_option)
    {
    case FQMI:
      {
      // call to fqmi procedure
      return criterion_fqmi
               (
               class_vector,
               indexPtr,
               curSubSetSize,
               sigma
               );

      }

    case FLDR:
      {
      // call to fldr procedure
      return criterion_fldr(class_vector, indexPtr, curSubSetSize);
      }

    case CUSTOM_CRITERION:
      {
      // call to custom criterion procedure
      return criterion_custom(class_vector, indexPtr, curSubSetSize);
      }

    default:
      {
      #ifdef NCURSES_OUTPUT
      endwin();
      #endif
      cerr << "Unknown criterion option" << endl;
      exit(1);
      }

    }

  }



ucolvec
sffs
  (
  const vector<ClassData>& classes_vector,
  const u32 crit
  )
  {
  // sigma for fqmi()
  double sigma = 1.0;

  // if criterion is FQMI compute sigma
  if(crit == 1)
    {
    // compute sigma
    sigma = compute_sigma(classes_vector);
    }

  // --- declarations of constants --- //

  // full feature set size
  const u32 n = classes_vector.size();

  // desired feature subset size
  const u32 d = n - 1;

  // set floating search stopping constraint
  int delta = 1;

  // set generalization level
  int genLevel = 1;

  // a multiplier of estimated delta
  const double DELTAMUL = 1.0;

  // a constant to be added to estimated delta
  const int DELTAADD = 1;

  // double overflow limit
  const double SAFEUP = 5e300;

  // ---declarations of variables --- //

  // of size n, stores 0/1 information on currently selected features
  ucolvec binPtr;

  // of size d, it is computed from bin, stores indexPtrs of selected
  // features
  ucolvec indexPtr;

  // temporary best subset array
  ucolvec tempBestSetPtr;

  // best value
  double tempBestValue;

  // so far best subsets
  ucolvec globalBestSetPtr;

  // so far best values
  colvec globalBestValuePtr;

  // indicates, if a better subset has been found during last step
  bool betterFound;

  // criterion function return value
  double value;

  // current subset size
  int curSubSetSize;

  // current subset size in integers
  int curSubSetSizeInt;

  // how to change curSubSetSize (when changing direction of search):
  // -2,-1,+1,+2 steps
  int curSubsetChange;

  // initial generalization level
  int initGenLevel;

  // auxiliary variable
  int aux;

  // beginning of block of identifiers
  int begin;

  // "pivot" - last identifier in a block
  int pivot;

  // identifies end of internal exhaustive search
  bool stopExhaustive;

  // identifies end of it all
  bool stopAll;

  // identifies the dimension, for which the algorithm should end
  int stopDimension = 0;

  // identifier 0 exhanged in case of backward search
  int id0 = 0;

  // identifier 2 exhanged in case of backward search
  int id2 = 2;

  // 0=sfs, 1=sbs current search direction
  bool sbs = 0;

  // stores current backtracking depth (+forward,-backward)
  int backtrackDepth = 0;

  // predicted delta estimate
  float predDeltaEst = 0.0;

  // number of direction changes (needed for delta averaging)
  int dirCount = 0;

  // number of steps
  long int stepsNum;

  // set stopDimension according to delta value
  if(delta == 0)
    {
    stopDimension = n;
    }
  else
    {
    if(delta > 0)
      {
      stopDimension = d + delta;

      if(stopDimension > n)
        {
        stopDimension = n;
        }

      }

    }

  // initialize back tracking depth
  backtrackDepth = 0;

  // initialize predicted delta estimate
  predDeltaEst = 0.0;

  // initialize direction counter
  dirCount = 0;

  // allocate memory for feature code word
  binPtr = zeros<ucolvec>(n + 1);

  // allocate memory for indices of selected features
  indexPtr = zeros<ucolvec>(n);

  // allocate memory for current best set
  tempBestSetPtr = zeros<ucolvec>(n);

  // allocate memory for so far best subsets
  globalBestSetPtr = zeros<ucolvec>((n * ( n + 1 ) ) / 2);

  // allocate memory for so far best subset values
  globalBestValuePtr.set_size(n);
  globalBestValuePtr.fill(-SAFEUP);

  // set initial generalization level to current generalization level
  initGenLevel = genLevel;

  // --- Beginning of printing info block according to type of    --- //
  // --- floating search                                          --- //

  // set current subset size equal to initial generalization level
  curSubSetSize = initGenLevel;

  // initialize stop centinel
  stopAll = false;

  // do while floating search ends
  do
    {
    // initialize temporary best set
    tempBestSetPtr.zeros();

    // initialize best value with underflow limit
    tempBestValue = -SAFEUP;

    // initialize binPtr for exhaustive step
    aux = initGenLevel;
    for(u32 i = 0; ((aux > 0) && ( i < n )); i++)
      {
      if(binPtr[i] == id0)
        {
        binPtr[i] = id2;
        aux--;
        }

      }

    // size of current subset size in number integer (usually 4 byte per
    // int)
    curSubSetSizeInt = curSubSetSize * sizeof( int );

    // --- estimate the number of steps --- //
    stepsNum = 1;

    // backward step sbs = 1
    if(sbs)
      {
      aux = curSubSetSize + initGenLevel;
      }
     else // else forward step
      {
      aux = n - ( curSubSetSize - initGenLevel );
      }

    for(u32 i = 0; i < initGenLevel; i++)
      {
      stepsNum *= aux - i;
      }

    for(u32 i = 2; i <= initGenLevel; i++)
      {
      stepsNum /= i;
      }
    // --- end of steps estimation --- //

    stopExhaustive = false;

    // do while internal exhaustive search ends
    do
      {
      // --- convert "binPtr" to "indexPtr" --- //
      aux = 0;
      for(u32 i = 0; i < curSubSetSize; i++)
        {
        while(binPtr[aux] <= 0)
          {
          aux++;
          }

        indexPtr[i] = aux;
        aux++;
        }
      // --- end of conversion --- //

      // evaluate criterion function
      value = criterion
                (
                classes_vector,
                indexPtr,
                curSubSetSize,
                crit,
                sigma
                );

      // if the new value is better than the best value already
      if(value > tempBestValue)
        {
        tempBestSetPtr = indexPtr;
        tempBestValue = value;
        }

      // iterate through pBegin to find the beginning of the block of
      // selected features
      for(begin = 0; binPtr[begin] != id2; begin++)
        {
        // just iterate
        }

      // iterate through pBegin to find the pivot
      for(pivot = begin; (pivot < n) && (binPtr[pivot] != id0); pivot++)
        {
        // just iterate
        }

      if(pivot == n)
        {
        stopExhaustive = true;
        }
      else
        {
        // remember the position of first 0 on the right
        aux = pivot;

        // find a real pivot
        do
          {
          pivot--;
          }
        while(binPtr[pivot] != id2);

        // shift pivot to the right
        binPtr[pivot] = id0;
        binPtr[aux] = id2;

        aux = 0;

        // --- run "aux" from left, "piv" from right. the 0,2     --- //
        // --- pairs found are changed to 2,0                     --- //
        do
          {
          pivot--;
          }
        while((pivot > 0) && (binPtr[pivot] != id2));

        while((aux < pivot) && (binPtr[aux] != id0))
          {
          aux++;
          }

        while((pivot - aux) > 0)
          {
          binPtr[pivot] = id0;
          binPtr[aux] = id2;

          do
            {
            pivot--;
            }
          while((pivot > 0) && (binPtr[pivot] != id2));

          while((aux < pivot) && (binPtr[aux] != id0))
            {
            aux++;
            }

          }
        // --- end of run --- //

        }

      }
    while(!stopExhaustive);

    // update if best value is found
    if(tempBestValue > globalBestValuePtr[curSubSetSize - 1])
      {
      memcpy
        (&globalBestSetPtr[(curSubSetSize * (curSubSetSize - 1)) / 2],
        tempBestSetPtr.memptr(),
        curSubSetSizeInt
        );

      globalBestValuePtr[curSubSetSize - 1] = tempBestValue;
      betterFound = true;
      }
    else
      {
      betterFound = false;
      }

    // if last step was sbs
    if(sbs)
      {
      // a better subset was found
      if(betterFound)
        {
        // --- convert according to format --- //

        // convert to sbs format
        if(curSubSetSize > genLevel)
          {
          curSubsetChange = -1;

          for(u32 i = 0; i < n; i++)
            {
            binPtr[i] = -1;
            }

          // freeze changes
          for(u32 i = 0; i < curSubSetSize; i++)
            {
            binPtr[tempBestSetPtr[i]] = 2;
            }

          }
        else // nothing to remove, convert to sfs format
          {
          for(u32 i = 0; i < n; i++)
            {
            binPtr[i] = 0;
            }

          for(u32 i = 0; i < curSubSetSize; i++)
            {
            binPtr[tempBestSetPtr[i]] = 1;
            }

          // change to sfs
          sbs = false;

          // change "binPtr" identifiers 0 and 2
          id0 = 0;
          id2 = 2;

          curSubsetChange = 1;
          }

        // --- end of conversion --- //

        }
      else // no improvement during last step (sbs)
        {
        // change to sfs
        sbs = false;

        // change to "binPtr" identifiers 0 and 2
        id0 = 0;
        id2 = 2;

        // forget last step and perform one new sfs step
        curSubsetChange = 2;
        }

      }
    else // last step was sfs
      {
      // if removing is possible, prepare sbs
      if(curSubSetSize > genLevel)
        {
        for(u32 i = 0; i < n; i++)
          {
          binPtr[i] = -1;
          }

        for(u32 i = 0; i < curSubSetSize; i++)
          {
          binPtr[tempBestSetPtr[i]] = 2;
          }

        // change to sbs
        sbs = true;

        // change to "binPtr" identifiers 0 and 2
        id0 = 2;
        id2 = 0;

        curSubsetChange = -1;

        }
      else // otherwise stay by sfs
        {
        curSubsetChange = 1;

        for(u32 i = 0; i < n; i++)
          {
          binPtr[i] = 0;
          }

        // freeze changes
        for(u32 i = 0; i < curSubSetSize; i++)
          {
          binPtr[tempBestSetPtr[i]] = 1;
          }

        }

      }

    // renew curSubSetSize and initGenLevel
    if(curSubsetChange == -1)
      {
      if((curSubSetSize == d) && (d % genLevel != 0))
        {
        curSubSetSize = d - ( d % genLevel );
        initGenLevel = d % genLevel;
        }
      else
        {
        curSubSetSize -= genLevel;
        initGenLevel = genLevel;
        }

      if(backtrackDepth < 0)
        {
        // continue
        backtrackDepth -= initGenLevel;
        }
      else
        {
        // begin
        backtrackDepth = -initGenLevel;
        }

      }
    else // else curSubsetChange == 1 or 2
      {
      if((curSubSetSize < d) && ((curSubSetSize + genLevel) > d))
        {
        curSubSetSize = d;
        initGenLevel = d % genLevel;
        }
      else
        {
        curSubSetSize += genLevel;
        initGenLevel = genLevel;
        }

      // if end of going down
      if(backtrackDepth < 0)
        {
        if(delta == -1)
          {
          // averaging
          predDeltaEst =
            ( ( dirCount * predDeltaEst ) - backtrackDepth ) /
            ( dirCount + 1 );

          dirCount++;
          }
        else // else maximize
          {
          if((-backtrackDepth) > predDeltaEst)
            {
            predDeltaEst = -backtrackDepth;
            }

          }

        }

      backtrackDepth = 0;
      }

    //  once more and renew "binPtr"
    if(curSubsetChange == 2)
      {
      // change to sfs format now with actualized "curSubSetSize"
      for(u32 i = 0; i < n; i++)
        {
        binPtr[i] = 0;
        }

      aux = (curSubSetSize * (curSubSetSize - 1)) / 2;

      for(u32 i = 0; i < curSubSetSize; i++)
        {
        binPtr[globalBestSetPtr[aux + i]] = 1;
        }

      if((curSubSetSize < d) && ((curSubSetSize + genLevel) > d))
        {
        curSubSetSize = d;
        initGenLevel = d % genLevel;
        }
      else
        {
        curSubSetSize += genLevel;
        initGenLevel = genLevel;
        }

      // no direction change
      }

    if(delta < 0)
      {
      stopDimension = int(d + (DELTAMUL * predDeltaEst) + DELTAADD);

      if(stopDimension > n)
        {
        stopDimension = n;
        }

      }

    // end if delta reached
    if(curSubSetSize > stopDimension)
      {
      stopAll = true;
      }

    }
  while(!stopAll);

  // ---  find the best value and best subset from all possible   --- //
  // --- subsets                                                  --- //
  double bestValue = globalBestValuePtr[0];
  int bestValueIndex = 0;

  // find the index of the best value
  for(u32 i = 1; i < n - 1; i++)
    {
    if(globalBestValuePtr[i] > bestValue)
      {
      bestValue = globalBestValuePtr[i];
      bestValueIndex = i;
      }

    }
  // --- end of best value and subset --- //

  // auxiliary variable
  aux = ((bestValueIndex + 1) * bestValueIndex) / 2;

  return globalBestSetPtr.rows(aux, aux + bestValueIndex);
  }