svmfusvmsmallopt.cpp

This is SvmFu, a package for training and testing support vector machines (SVMs). It s written in C

//     This is a part of the SvmFu library, a library for training
//     Support Vector Machines.
//     Copyright (C) 2000  rif and MIT
//
//     Contact: [e-mail]

//     This program is free software; you can redistribute it and/or
//     modify it under the terms of the GNU General Public License as
//     published by the Free Software Foundation; either version 2 of
//     the License, or (at your option) any later version.

//     This program 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.

//     You should have received a copy of the GNU General Public
//     License along with this program; if not, write to the Free
//     Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
//     MA 02111-1307 USA

#include "SvmFuSvmSmallOpt.h"

#ifdef CallTemplate
#undef CallTemplate
#endif

#define CallTemplate(rettype, funcname) \
template<class DataPt, class KernVal> \
rettype SvmSmallOpt<DataPt, KernVal>::funcname

#ifdef InlineTemplate
#undef InlineTemplate
#endif

#define InlineTemplate(rettype, funcname) \
template<class DataPt, class KernVal> \
rettype SvmSmallOpt<DataPt, KernVal>::funcname


template<class DataPt, class KernVal> 
SvmSmallOpt<DataPt, KernVal>::SvmSmallOpt
(int svmSize, const IntVec y, 
 DataPt * trnSetPtr,
 const KernVal (*kernProdFuncPtr)(const DataPt &pt1, const DataPt &pt2),
 double C, double tol, double eps,
 SvmKernCache<DataPt, KernVal> *kernCache)
  : SvmBase<DataPt, KernVal>(svmSize, y, trnSetPtr, 
			     kernProdFuncPtr, C, eps),
    fullCacheOnP_(true), tol_(tol), alphaOffset_(0)
{
  int size = getSize();

  oldAlphas_ = new double[size];
  cachedGradients_ = new double[size];

  workingSet_ = new int[size];
  nonWorkingSet_ = new int[size];
  totalSet_ = new int[size];
  offsetVec_ = new double[size];
      
  for (int i = 0; i < size; i++) {
    oldAlphas_[i] = 0;
    cachedGradients_[i] = -1;

    workingSet_[i] = -1;
    nonWorkingSet_[i] = -1;
    totalSet_[i] = i;

    offsetVec_[i] = 0;
  }

  if (kernCache == 0) {
    kernCache_ = new SvmKernCache<DataPt, KernVal>
      (svmSize, 0, svmSize, trnSetPtr, kernProdFuncPtr);
    kernCacheCreatedP_ = true;
  } else {
    kernCache_ = kernCache;
    kernCacheCreatedP_ = false;
  }
}
  
template<class DataPt, class KernVal> 
SvmSmallOpt<DataPt, KernVal>::~SvmSmallOpt() {
  delete[] oldAlphas_;
  delete[] cachedGradients_;

  delete[] workingSet_;
  delete[] nonWorkingSet_;
  delete[] totalSet_;
    
  delete[] offsetVec_;
  if (kernCacheCreatedP_) {
    delete kernCache_;
  }
}

CallTemplate(void, optimize)() {
  bool done = false;
  bool success;

  stepsTaken_ = 0;
  int numInteriorSteps = 0;

  while(!done) {
    setTopBottomPosNeg();
    chooseExamples();

    success = false;
    if (maxVal_ > getTolerance()) {
      success = takeStep(ex1_, ex2_);

      if (success) {
	stepsTaken_++;

	if (unbndSupVecP(ex1_) && unbndSupVecP(ex2_)) {
	  numInteriorSteps++;
	} else {
	  numInteriorSteps = 0;
	}
	
	if (numInteriorSteps == 100 && fullCacheOnP()) {
	  useUnbndCacheOnly();
	}
      } else {
	cout << maxVal_ << endl;
	cout << ex1_ << " " << ex2_ << endl;
	cout << topPosInd_ << " " << topNegInd_ << " " << bottomPosInd_ << " " << bottomNegInd_ << endl;
	cerr << "ERROR: takeStep failed.  This shouldn't happen." << endl;
	exit(-1);

	done = true;
	break;
      }
    }

    if (success == false) { // We can't take a step under current conditions,
      if (!fullCacheOnP()) {
	useFullCache();
	numInteriorSteps = 0;
      } else {
	done = true;
      }
    }
  }
}

CallTemplate(void, setTopBottomPosNeg)() {
  int *setPtr;
  int setSize;

  topPosInd_ = bottomPosInd_ = topNegInd_ = bottomNegInd_ = -1;
    
  if (fullCacheOnP()) {
    setPtr = totalSet_;
    setSize = getSize();
  } else {
    setPtr = workingSet_;
    setSize = workingSetSize_;
  }

  for (int i = 0; i < setSize; i++) {
    int ex = setPtr[i];
    double cG = cachedGradients_[ex];

    double C = getC(ex);
    double alpha = getAlpha(ex);

    double eps = getEpsilon();

    if (getY(ex) == 1) { // a Positive Example
      if ((topPosInd_ == -1 || cG > topPos_) &&
	  alpha > tol_){
	topPos_ = cG;
	topPosInd_ = ex;
      }
	
      if ((bottomPosInd_ == -1 || cG < bottomPos_) &&
	  alpha <= C - tol_) {
	bottomPos_ = cG;
	bottomPosInd_ = ex;
      }
    } else { // a Negative Example
      if ((topNegInd_ == -1 || cG > topNeg_) && 
	  alpha > tol_) {
	topNeg_ = cG;
	topNegInd_ = ex;
      }
	
      if ((bottomNegInd_ == -1 || cG < bottomNeg_) &&
	  alpha <= C - tol_) {
	bottomNeg_ = cG;
	bottomNegInd_ = ex;
      }
    }
  }      
}
  
CallTemplate(void, chooseExamples)() {
  maxVal_ = 0;
    
  if ((topPosInd_ != -1 && bottomPosInd_ != -1) &&
      (topPos_ - bottomPos_) >= maxVal_) {
    ex1_ = topPosInd_; ex2_ = bottomPosInd_;
    maxVal_ = topPos_ - bottomPos_;
  }
    
  if ((topNegInd_ != -1 && bottomNegInd_ != -1) &&
      (topNeg_ - bottomNeg_) >= maxVal_) {
    ex1_ = topNegInd_; ex2_ = bottomNegInd_;
    maxVal_ = topNeg_ - bottomNeg_;
  }

  if ((topPosInd_ != -1 && topNegInd_ != -1) &&
      (topPos_ + topNeg_) >= maxVal_) {
    ex1_ = topPosInd_; ex2_ = topNegInd_;
    maxVal_ = topPos_ + topNeg_;
  }
    
  if ((bottomPosInd_ != -1 && bottomNegInd_ != -1) &&
      -(bottomPos_ + bottomNeg_) >= maxVal_) {
    ex1_ = bottomPosInd_; ex2_ = bottomNegInd_;
    maxVal_ = -(bottomPos_ + bottomNeg_);
  }
}

CallTemplate(bool, takeStep)(int ex1, int ex2) {

  // cout << "SvmSmallOpt::takeStep(): " << ex1 << " " << ex2 << endl;

  if (ex1 == ex2) {
    cerr << "WARNING: ex1 == ex2 in takeStep.  This shouldn't happen."
	 << endl;
    return false;
  }

  double a1 = getAlpha(ex1);
  double a2 = getAlpha(ex2);
  int y1 = getY(ex1);
  int y2 = getY(ex2);
  double C1 = getC(ex1);
  double C2 = getC(ex2);

  double err1 = y1*cachedGradients_[ex1];
  double err2 = y2*cachedGradients_[ex2];
  
  int s = y1*y2;
    
  double La2, Ha2, gamma;
    
  if (s == -1) {
    gamma = a2-a1;
    La2 = gamma < 0 ? 0 : gamma;
    Ha2 = C2 < gamma+C1 ? C2 : gamma+C1;
  } else {
    gamma = a1+a2;
    La2 = gamma < C1 ? 0 : gamma - C1;
    Ha2 = gamma < C2 ? gamma : C2;
  }

  double k11 = kernCache_->cachedKernProd(ex1, ex1);
  double k12 = kernCache_->cachedKernProd(ex1, ex2);
  double k22 = kernCache_->cachedKernProd(ex2, ex2);
    
  double a1new, a2new, a2want=0;
  double eta = 2*k12-k11-k22;
  double eps = getEpsilon();
    
  if (eta < -eps) {
    a2want = a2-y2*(err1-err2)/eta;
    a2new = a2want < La2 ? La2 : a2want;
    a2new = a2new > Ha2 ? Ha2 : a2new;
  } else {
    // This can possibly be changed, see old version of code.
    cerr << "WARNING: eta too small in takeStep.  This shouldn't happen."
	 << endl;
    cerr << "Are examples " << kernCache_->getColToTrnSet(ex1) << " and " << kernCache_->getColToTrnSet(ex2)
	 << " (counting from 0) identical points from opposite classes?" << endl;

    return false;
  }
    
  a1new = a1 + s*(a2-a2new);
    
  if (fabs(a1new-a1) < eps*(a1new+a1+1+eps)) {
    cerr << "WARNING: step too small in takeStep.  This shouldn't happen."
	 << endl;
    cerr << a1new << " " << a1 << " " << a2 << " " << a2new << " " << a2want << " " << k11 << " " << k12 << " " << k22 << " " << eta << " " << err1 << " " << err2 << " " << y1 << " " << y2 << endl;
    return false;
  }
    
  setAlpha(ex1, a1new);
  setAlpha(ex2, a2new);
  
  return true;
}
  
CallTemplate(void, setAlpha)(int ex, double newAlpha) {
  double oldAlpha = getAlpha(ex);
  double alphaDiff = newAlpha - oldAlpha;

  if (fabs(alphaDiff) > getEpsilon()) {
    const KernVal * kernRow = kernCache_->cachedKernProdRowPtr(ex);

    int y = getY(ex);
      
    int *setPtr;
    int setSize;

    if (fullCacheOnP()) {
      setPtr = totalSet_;
      setSize = getSize();
    } else {
      setPtr = workingSet_;
      setSize = workingSetSize_;
    }
      
    for (int i = 0; i < setSize; i++) {
      int ex2 = setPtr[i];
      int y2 = getY(ex2);
	
      cachedGradients_[ex2] += y*y2*alphaDiff*kernRow[ex2];
    }

    SvmBase<DataPt, KernVal>::setAlpha(ex, newAlpha);
  }
}
  
CallTemplate(void, useUnbndCacheOnly)() {
  if (!fullCacheOnP_) {
    cerr << "WARNING: Redundant attempt to switch to unbndCache only; no effect." << endl;
    return;
  }
   
  // double dOF = dualObjFunc();
    
  DoubleVec alphas = getAllAlphas();
  DoubleVec C_Vec = getCVec();

  int size = getSize();
  workingSetSize_ = nonWorkingSetSize_ = 0;

  double eps = getEpsilon();

  for (int i = 0; i < size; i++) {
    //      if ((alphas[i] > eps && alphas[i] < C_Vec[i] - eps) ||
    //(alphas[i] <= eps && cachedGradients_[i] <= tol) ||
    // (alphas[i] >= C_Vec[i]-eps && cachedGradients_[i] >= -tol)) {
    if ((alphas[i] > eps && alphas[i] < C_Vec[i] - eps)) {
      workingSet_[workingSetSize_++] = i;
    } else {
      nonWorkingSet_[nonWorkingSetSize_++] = i;
    }

    oldAlphas_[i] = alphas[i];
  }

  fullCacheOnP_ = false;

  delete[] alphas;
  delete[] C_Vec;
}

CallTemplate(void, useFullCache)() {
  if (fullCacheOnP_) {
    cerr << "WARNING: Redundant attempt to switch to fullCache; no effect." << endl;
    return;
  }

  DoubleVec alphas = getAllAlphas();

  for (int i = 0; i < workingSetSize_; i++) {
    int wEx = workingSet_[i];
    if (alphas[wEx] != oldAlphas_[wEx]) {
      double alphadiff = alphas[wEx]-oldAlphas_[wEx];
      const KernVal *kernRow = kernCache_->cachedKernProdRowPtr(wEx);
      int yW_Ex = getY(wEx);
      for (int j = 0; j < nonWorkingSetSize_; j++) {
	int nEx = nonWorkingSet_[j];
	cachedGradients_[nEx] += getY(nEx)*yW_Ex*alphadiff*kernRow[nEx];
      }
    }
  }
    
  fullCacheOnP_ = true;
  delete[] alphas;
}

CallTemplate(double, dualObjFunc)() const {
  if (!fullCacheOnP()) {
    cerr << "WARNING: Full cache is not on; computed dual objective may be inaccurate." << endl;
  }
  
  double result = 0.0;
  int size = getSize();
    
  for (int i = 0;  i < size; i++) {
    // double alpha = getAlpha(i);
    // result += alpha;
    // result -= .5*alpha*(cachedGradients_[i]+1);
    result += getAlpha(i)*(1-cachedGradients_[i]+getY(i)*getOffset(i));
  }
    
  result /= 2;
  result += alphaOffset_;

  return result;
}

InlineTemplate(bool, fullCacheOnP)() const {
  return fullCacheOnP_;
}

CallTemplate(double, outputAtTrainingExample)(int ex) const {
  if (!fullCacheOnP()) {
    cerr << "WARNING: Full cache is not on in outputAtTrainingExample; computed output may be inaccurate." << endl;
  }
    
  double output = getY(ex)*(cachedGradients_[ex]+1) + getB() + getOffset(ex);
  
  return output;
}

CallTemplate(void, fixB)(bool printInfo) {
  const IntVec USVs = getUnbndSupVecIDsPtr();
  int numUSVs = getNumUnbndSupVecs();

  // bM will be mean, bV will be variance.
  // We'll compute bV in one pass using Var(X) = E[X^2] - E[X]^2
  double bM = 0, bV = 0;
  for (int i = 0; i < numUSVs; i++) {
    double bEst = -getY(USVs[i])*cachedGradients_[USVs[i]]; 
      
    bM += bEst;
    bV += bEst*bEst;
  }

  bM /= numUSVs;
  bV = (bV/numUSVs)-bM*bM;
    
  setB(bM);
    
  if (printInfo) {
    cout << "Computed b using " << numUSVs << " USVs: b = " << bM << 
      ", var(b) = " << bV << endl;
  }
}

InlineTemplate(void, addToCachedGradient)(int ex, double amt) {
  cachedGradients_[ex] += amt;
}

CallTemplate(void, setTolerance)(double tol) {
  // This cannot change whether a point is an SV or not.
  // It CAN affect whether a point violates the KKT conditions, but
  // we don't save that info, so it's no problem.
  tol_ = tol;
}

CallTemplate(void, setOffsetVec)(const DoubleVec offsetVec) {
  int size = getSize();

  for (int i = 0; i < size; i++) { 
    // offsetVec_[i] = offsetVec[i]; 
    addToCachedGradient(i, getY(i)*offsetVec[i]);
  }
}

InlineTemplate(double, getOffset)(int ex) const { return offsetVec_[ex]; }

InlineTemplate(double, getCachedGradient)(int ex) const {
  return cachedGradients_[ex];
}

InlineTemplate(void, setAlphaOffset)(double alphaOffset) {
  alphaOffset_ = alphaOffset;
}
  
InlineTemplate(double, getAlphaOffset)() const {
  return alphaOffset_;
}
  
CallTemplate(int, getStepsTaken)() const {
  return stepsTaken_;
}

InlineTemplate(double, getTolerance)() const { return tol_; }


#include "SvmFuSvmDataPoint.h"
#define IterateTypes(datatype, kerntype) \
    template class SvmSmallOpt<DataPoint<datatype>, kerntype>;
#include "SvmFuSvmTypes.h"